Iterative algorithm for offsets, scale and rotate estimation for television image superposition with additive and multiplicative noise

Full Text

Введение

Совмещение телевизионных изображений является важной задачей в научных исследованиях и технических приложениях. Вопросы совмещения решаются в медицине (совмещение гистологических изображений, трехмерная реконструкция и т. д.) [1], в авиационной технике (совмещение подстилающей поверхности) [2; 3], в железнодорожной сфере (поиск эталонных изображений в видеопотоке для оценки смещения рельсов относительно шпальной решетки), а также в различных системах индексации данных (сопоставление портретных фотографий, поиск изображения по фрагменту и т. д.) [4–7].

Из-за большого объема данных, содержащихся в изображениях, использование алгоритмов оценки параметров, основанные на методе полного перебора, приводит к существенному снижению производительности алгоритма совмещения.

Оценка параметров происходит по метрике (как правило, по коэффициенту корреляции). Количество метрик равно количеству сочетаний оцениваемых параметров. При этом снижение производительности связано с преобразованием изображений, которых будет столько же, сколько и сочетаний параметров.

По этой причине актуальной проблемой является разработка принципов оценки параметров, отличных от принципов полного перебора. Они должны обеспечивать более высокую производительность. Одним из таких принципов является итерационный метод оценки параметров. В данной статье представлен алгоритм совмещения, в основе которого лежит итерационная процедура.

Обзор существующих работ

Задача совмещения изображений известна очень давно и по этой причине разработаны различные подходы, позволяющие повысить скорость обработки.

Также следует отметить, что ряд подходов не являются универсальными и позволяют решить задачу, если выполняются определенные ограничения.

Если совмещаемые изображения связаны между собой только плоскопараллельным смещением, то для повышения скорости вычислений может быть использовано:

– расчёт свёрток на базе преобразования Фурье [4; 7–9];

– метод «пирамид изображений» для уменьшения размерности данных в 2ⁿ × 2ⁿ раз [10; 11].

Если совмещаемые изображения связаны между собой плоскопараллельным смещением и смещение сравнимо по величине с шагом дискретизации, то для повышения скорости вычислений может быть использован метод линеаризации [12]. Производиться разложение сигналов в ряд и параметры оцениваются в результате решения системы линейных уравнений.

Также существует модификация этого метода для учета поворота, если угол не превышает 15°. Авторы статьи [13] предлагают кроме разложения сигналов в ряд вводить замену тригонометрических функций синуса и косинуса на приближенные значения. Повышение скорости вычислений также обеспечивается за счет решения в аналитическом виде.

Если совмещаемые изображения связаны между собой смещением, масштабом и поворотом, то для повышения скорости вычислений может быть использован метод сопоставления реперных точек. Реперные точки – это, как правило, локальные экстремумы. Сопоставление реперных точек проводиться за счет совмещения фрагментов в окрестности этих точек. Фрагменты содержат меньшее количество элементов, чем исходные изображения. За счет этого получается выигрыш в скорости обработки. Этот метод часто ассоциируется с дескрипторами [14–16], который представляет собой вектор, описывающий фрагмент изображения, малочувствительный к масштабным искажениям и повороту. Серьезным недостатком этого метода является неверное сопоставление отдельных фрагментов и, соответственно, последующее исправление данной ситуации.

Если смещение не превышает +/–10 % исходного изображения, масштаб +/–20 %, а поворот +/–30°, то можно использовать метод раздельной оценки [4–7]. В методе отдельно друг от друга оцениваются две группы параметров:

– смещения вдоль координатных осей;

– масштаб и поворот.

При этом используется декартовая и логарифмически-полярная система координат. Метод позволяет уменьшить число преобразований изображений (в сравнении с методом проверки всех возможных комбинаций), тем самым обеспечивая повышение скорости обработки.

Разработанный итерационный алгоритм наиболее близок к алгоритмам раздельной оценки параметров совмещения. Но он характеризуется менее строгими органическими по масштабу (+/–50 %) и повороту (+/–180°).

Исходные данные

На рис. 1 представлены изображения, которые будут использованы для совмещения.

Рис. 1. Исходные изображения: f(x, y) (а), g(x, y) (б)

Fig. 1. Initial images: f(x, y) (a), g(x, y) (b)

На рис. 1, а показано изображение f(x, y), которое является фрагментом изображения g(x, y). На рис. 1, б на изображении g(x, y) показана область, соответствующая f(x, y).

Параметры совмещения:

– смещение вдоль оси абсцисс: h = 200 пикселей;

– смещение вдоль оси ординат: p = 20 пикселей;

– масштаб: a = 0,8;

– поворот: j = 30°;

– мультипликативная помеха: l = 1,2;

– аддитивная помеха: g = 20 у.е.

Примечание: у.е. описывает яркости пикселя изображения, кодируемая числом от 0 до 255 у.е.

Эти сигналы будут использованы для проведения численного моделирования для тестирования разработанного алгоритма.

Исходные данные

Изображение представляет собой двумерный сигнал с равномерным шагом дискретизации.

Параметрическая модель, связывающая между собой изображения, имеет следующий вид:

$\left\{\begin{array}{l}f\left(x_i,y_i\right)=s\left(x_i,y_i\right)+k\left(x_i,y_i\right);\\ g\left(x_i,y_i\right)=\lambda s\left(u_i,w_i\right)+\gamma +m\left(x_i,y_i\right),\end{array}\right.$(1)

$\left\{\begin{array}{l}{u}_i=\alpha x_i\cos \left(\phi \right)-\alpha y_i\sin \left(\phi \right)+h;\\ w_i=\alpha x_i\sin \left(\phi \right)+\alpha y_i\cos \left(\phi \right)+p,\end{array}\right.$(2)

где f(x, y), g(x, y) – совмещаемые изображения, h, p – смещения вдоль оси абсцисс и ординат, $\alpha$ – масштаб, $\phi$ – поворот, $\lambda , \gamma$  – мультипликативная и аддитивная помеха.

Совмещение определяется параметрами {h, p, $\alpha , \phi , \lambda , \gamma \}.$ Для оценки параметров будет использован критерий максимума коэффициента корреляции:

$\begin{array}{l}R\left(\theta \right)=\left\{\left(\sum _{i=1}^{N}g\left(x_i,y_i\right)\cdot f\left(x_i,y_i,\theta \right)\right)/N-\right.\\ \left.-\left(\left(\sum _{i=1}^{N}g\left(x_i,y_i\right)\right)/N\right)\left(\left(\sum _{i=1}^{N}f\left(x_i,y_i,\theta \right)\right)/N\right)\right\}/\\ /\left\{{\left(\sum _{i=1}^N{g}^2\left(x_i,y_i\right)-{\left(\sum _{i=1}^{N}g\left(x_i,y_i\right)\right)}^2/N\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right.\times \\ \left.\times {\left(\sum _{i=1}^N{f}^2\left(x_i,y_i,\theta \right)-{\left(\sum _{i=1}^{N}f\left(x_i,y_i,\theta \right)\right)}^2/N\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right\},\end{array}$ (3)

$\widehat{\theta }=\underset{\theta }{\mathrm{argmax}}\left(R\left(\theta \right)\right),$ (4)

где $\theta =\{h, p, \alpha , \phi , \lambda , \gamma \},$ N – это количество пикселей изображения.

Критерий, основанный на вычислении коэффициента корреляции, не зависит от параметров $\lambda$ и $\gamma .$ В этом можно убедиться, рассчитав коэффициент корреляции для двух произвольных векторов {I_i, Y_i} и для $\{Z_i=\lambda I_i+\gamma , Y_i\}.$ Коэффициент корреляции для {I_i, Y_i} будет равен коэффициенту для {Z_i, Y_i}.

По этой причине, если оценены параметры $\{\widehat{h}, \widehat{p}, \widehat{\alpha }, \widehat{\phi }\},$ то параметры $\{\lambda , \gamma \}$ могут быть оценены по методу наименьших квадратов:

 $\begin{array}{l}〈\widehat{\lambda },\widehat{\gamma }〉=\\ =\underset{\lambda ,\gamma }{\mathrm{argmin}}\left(\sum _{i=1}^N{\left(\lambda f\left(x_i,y_i,\theta \right)+\gamma -g\left(x_i,y_i\right)\right)}^2\right).\end{array}$(5)

$\widehat{\lambda }=\frac{Sg\cdot Sf-N\cdot Sfg}{S{f}^2-N\cdot Ef},$ (6)

$\widehat{\gamma }=\frac{Sg\cdot Sfg-Sg\cdot Ef}{S{f}^2-N\cdot Ef},$ (7)

где 

$Sf=\sum _{i=1}^{N}f\left(x_i,y_i,\theta \right), Sg=\sum _{i=1}^{N}g\left(x_i\right),$

$Sfg=\sum _{i=1}^{N}f\left(x_i,y_i,\theta \right)g\left(x_i\right),$

$Ef=\sum _{i=1}^N{f}^2\left(x_i,y_i,\theta \right), \theta =\left\{\widehat{h},\widehat{p},\widehat{\alpha },\widehat{\phi }\right\},$

$f\left(x_i,y_i,\theta \right)$ – сигнал после преобразования с учетом найденный смещений, масштаба и угла поворота.

Таким образом, количество параметров, которые достаточно оценить уменьшается с шести до четырех.

Оценка масштаба и поворота

Если известно соответствие между двумя точками на изображениях, то можно оценить масштаб и поворот, используя логарифмически-полярную систему координат.

Пусть точка (x₀, y₀) изображения f соответствует точке (u₀, w₀) изображения g и пусть эти точки связаны между собой выражением (2), т. е. индекс i = 0. Пусть некоторая точка (x, y) изображения f соответствует точке (u, w) изображения g.

Обозначим координату точки (x, y) относительно (x₀, y₀) в логарифмически-полярной системе координат как (r, ang). Их значения будут определяться по формулам:

$\rho =\sqrt{{\left(x-x_0\right)}^2+{\left(y-y_0\right)}^2},$ (8)

$r={\text{log}}_2\left(\rho \right),$ (9)

$ang=\text{ atan}2\left(x-x_0,y-y_0\right).$ (10)

Примечание: atan2 – это функция с двумя параметрами, возвращает значение арктангенса выражения (y – y₀)/(x – x₀) в радианах, в отличие от арктангенса имеет область значений (–π; π) (в арктангенса – (–π/2; π/2)).

Координаты точки (u, w) относительно (u₀, w₀) в логарифмически-полярной системе координат обозначим как (r', ang'). Ниже представлен вывод выражений для них.

$u-u_0=\alpha \left(\cos \right(\phi )x-\sin (\phi \left)y\right).$

$w-w_0=\alpha \left(\sin \right(\phi )x+\cos (\phi \left)y\right).$

$\begin{array}{l}{\rho }^{\text{'}}=\sqrt{{\left(u-u_0\right)}^2+{\left(w-w_0\right)}^2}=\\ =\left|\alpha \right|\sqrt{{\left(x-x_0\right)}^2+{\left(y-y_0\right)}^2}.\end{array}$

С учетом того, что $\alpha >\mathrm{0,}$ то: 

${\rho }^{\text{'}}=\alpha \sqrt{{\left(x-x_0\right)}^2+{\left(y-y_0\right)}^2}=\alpha \rho .$

$r\text{'}={\log }_2\left(\rho \text{'}\right)={\log }_2\left(\alpha \right)+{\log }_2\left(\rho \right)={\log }_2\left(\alpha \right)+r.$ (11)

Рассмотрим выражение:

$\begin{array}{l}\text{tg}\left(an{g}^{\text{'}}\right)=\frac{w-w_0}{u-u_0}=\\ =\frac{\sin \left(\phi \right)\left(x-x_0\right)+\cos \left(\phi \right)\left(y-y_0\right)}{\cos \left(\phi \right)\left(x-x_0\right)-\sin \left(\phi \right)\left(y-y_0\right)}.\end{array}$

Если умножить числитель и знаменатель дроби на выражение $\frac{1}{\cos \left(\phi \right)(x-x_0)},$ то выражение tg(ang') можно записать в виде:

$\begin{array}{l}\text{tg}\left(an{g}^{\text{'}}\right)=\frac{\text{tg}\left(\phi \right)+\frac{\left(y-y_0\right)}{\left(x-x_0\right)}}{1-g\left(\phi \right)\frac{\left(y-y_0\right)}{\left(x-x_0\right)}}=\\ =\frac{\text{tg}\left(\phi \right)+\text{tg}\left(ang\right)}{1-g\left(\phi \right)\text{tg}\left(ang\right)}=\text{tg}\left(\phi +ang\right).\end{array}$

$an{g}^{\text{'}}=ang+\phi .$ (12)

Таким образом, можно записать, что:

$(r^{\text{'}},an{g}^{\text{'}})=(r+{\log }_2(\alpha ),ang+\phi ).$ (13)

На рис. 2 показан пример оценки масштаба и поворота по изображениям в логарифмически-полярной системе координат. На рис. 2, а показано изображение g. На рис. 2, б показан фрагмент, соответствующий изображению f. Он выделен прямоугольником. Смещения вдоль осей определяют масштаб и поворот.

Рис. 2. Оценка масштаба и поворота

Fig. 2. Evaluation of scale and rotation

При известном значении масштаба и поворота изображение f можно преобразовать таким образом, чтобы оно являлось фрагментом изображения g, которое можно совместить плоскопараллельным смещением. Оценить смещения можно известными способами [8–11].

Сложность реализации описанного подхода заключается в том, что не известны реперные точки (x₀, y₀) и (u₀, w₀) на изображениях f и g, которые позволили бы оценить сначала масштаб и поворот, а затем смещения.

Однако, если координаты этих точек известны с некоторой погрешностью, то разработанный алгоритм позволяет уточнить первоначально выбранные координаты реперных точек, тем самым обеспечивая совмещение телевизионных сигналов.

Описание алгоритма

Исходными данными являются: f(x_i, y_i), i = 1, .. N – первое изображение, g(u_j, w_j), j = 1, .. K – второе изображение.

Примечание: координаты пикселов сигнала g переобозначены как (u, w) вместо (x, y), чтобы было очевидно, какое изображение обрабатывается на том или ином шаге алгоритма.

Выходными данными являются: $\widehat{h}, \widehat{p}$ – смещения вдоль координатных осей; $\widehat{\phi }, \widehat{\alpha }$ – угол поворота и масштаб; $\widehat{\gamma }, \widehat{\lambda }$  – значение аддитивной и мультипликативной помехи.

Алгоритм состоит из следующих шагов.

1. Загрузка изображений f(x_i, y_i) и g(u_j, w_j).
2. Предварительный выбор точек (x₀, y₀) и (u₀, w₀).

Примечание: так как изображение f является фрагментом изображения g, то координаты (x₀, y₀) фиксируются и не меняются; для удобства представления в логарифмически-полярной системе они равняются координатам центрально пикселя (т. е. при размере 128 × 128 пикселей, (x₀, y₀) = (64, 64)); координаты (u₀, w₀) неизвестны, по этой причине используется метод перебора, однако шаг может выбираться из условия 0,1 размера изображения, при размерах 128 × 128 пикселей шаг был выбран равным 10 пикселей (т. е. немного меньше, чем 12,8); другими словами проверялись координаты u₀ = 1, 11 .. 381, w₀ = 1, 11, .. 281 (размер изображения g равняется 288 × 384).

3. Выбор количества итераций L (примеч.: в работе L = 4).
4. Инициализация матрицы преобразования в соответствии с количеством итераций:

$M_1=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right), M_2=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right) M{}_L=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right).$

5. Цикл по количеству итераций, k = 1.

5.1. Определение промежуточной матрицы преобразования Mres:

если k = 1, то Mres = M₁,

если k = 2, то Mres = M₂M₁,

…

если k = L, то Mres = M_L···M₂M₁.

5.2. Вычисление положение точек первого изображения (x_j', y_j') в соответствии с матрицей преобразования Mres:

x'_i = x_iMres(1,1) + y_iMres(1,2) + Mres(1,3),

y'_i = x_iMres(2,1) + y_iMres(2,2) + Mres(2,3),

5.3. Вычисление координат реперных точек, относительно которых происходит преобразование из декартовой в логарифмически-полярную систему координат:

X₀ = x₀Mres(1,1) + y₀Mres(1,2) + Mres(1,3),

Y₀ = x₀Mres(2,1) + y₀Mres(2,2) + Mres(2,3),

U₀ = X₀,

W₀ = Y₀.

5.4. Вычисление матрицы M_k по данным f(x'_i, y'_i), (X₀, Y₀), g(u_j, w_j), (U₀, W₀) согласно процедуре раздельной оценки параметров (процедура описана ниже).

5.5. Проверка на выход из цикла: k = k +1, если k ≥ L, то выход из цикла, в противном случае переход к п. 5.1.

6. Оценка параметров $\widehat{h}, \widehat{p}, \widehat{\phi }, \widehat{\alpha }$  по матрице Mres:

$\widehat{h}=Mres\left(\mathrm{1,3}\right), \widehat{p}=Mres\left(\mathrm{2,3}\right),$

$\widehat{\alpha }=\sqrt{Mre{s}^2\left(\mathrm{1,1}\right)+Mre{s}^2\left(\mathrm{1,2}\right)},$

$\widehat{\phi }=\arcsin \left(\frac{Mres\left(\mathrm{2,1}\right)}{\sqrt{Mre{s}^2\left(\mathrm{1,1}\right)+Mre{s}^2\left(\mathrm{1,2}\right)}}\right)180/\pi .$

7. Оценка параметров $\widehat{\gamma }, \widehat{\lambda }$ по формулам (6-7)
8. Выход из программы.

Ниже описана процедура раздельной оценки параметров (п. 5.4).

Входные данные: f(x'_i, y'_i), i = 1, .. N – первое изображение, g(u_j, w_j), j = 1, .. K – второе изображение. (X₀, Y₀), (U₀, W₀) – реперные точки.

Выходные данные: M – матрица преобразования размером [3 × 3].

1. 1. Преобразование координат (x'_i, y'_i) в логарифмически-полярную систему координат относительно точки (X₀, Y₀).
   2. Преобразование координат (u_j, w_j) в логарифмически-полярную систему координат относительно точки (U₀, W₀).
   3. Совмещение сигналов f и g в логарифмически-полярной системе. В результате оцениваются смещения ${\log }_2\left(\alpha r\right)$ и $\phi r$ (см. рис. 2).
   4. Вычисление промежуточной матрицы преобразования:

$\begin{array}{l}mA=\left[\begin{array}{cc}\alpha r\cos \left(\phi r\right)& -\alpha r\sin \left(\phi r\right)\\ \alpha r\sin \left(\phi r\right)& \alpha r\cos \left(\phi r\right)\\ 0& 0\end{array}\right.\\ \left.\begin{array}{c}{U}_0-\left(x_0\cos \left(\phi r\right)-y_0\sin \left(\phi r\right)\right)\\ W_0-\left(x_0\sin \left(\phi r\right)+y_0\cos \left(\phi r\right)\right)\\ 1\end{array}\right].\end{array}$

5. Преобразование координат (x'_i, y'_i) в соответствии с матрицей mA. Формирование (x''_i, y''_i).
6. Совмещение сигналов f(x''_i, y''_i) и g(u, w) в декартовой системе координат. В результате оцениваются смещения вдоль координатах осей hr и pr.
7. Вычисление промежуточной матрицы преобразования:

$mB=\left[\begin{array}{ccc}1& 0& hr\\ 0& 1& pr\\ 0& 0& 1\end{array}\right].$

8. Вычисление выходной матрицы М:

$M=mB\cdot mA.$

Численное моделирование

Для определения погрешности оценки параметров совмещения было проведено численное моделирование, которое заключалось в добавлении помехи с гауссовым распределением к ординатам одномерных сигналов.

При выбранном значении мощности (дисперсии помехи, ${\sigma }_n^2)$ проводилась серия опытов, по которым оценивалась погрешность оценки параметра.

В качестве погрешности выступало значение среднеквадратичной величины (RMS):

$RMS=\sqrt{\frac{{\displaystyle \sum _{i=1}^n{\left(a_i-a\right)}^2}}{n}},$ (14)

где a – действительное значение параметра; a_i – измеренное значение параметра в i-м опыте; n – количество опытов.

По результатам численного моделирования был построен график зависимости RMS от среднеквадратичного отклонения (СКО) помехи $\left({\sigma }_n\right).$

Результаты численного моделирования приведены на рис. 3.

Рис. 3. Погрешность оценки параметров

Fig. 3. Error in parameter estimation

Разработанный алгоритм сравнивался с алгоритмом полного перебора: проверялись все возможные варианты точки (u₀, w₀), т. е. u₀ = 1,2 .. 384, w₀ = 1,2 .. 288. И для каждой точки происходило преобразование в логарифмически-полярную систему координат. Погрешности у обоих алгоритмов почти совпадают. Различия наблюдаются при мощности шума ${\sigma }_n>20$ у.е. Различия возникают вследствие того, что итерационный процесс не приводит к верной оценке параметров за четыре итерации. При увеличении количества итераций до шести графики погрешностей у сравниваемых алгоритмов совпадают. Однако, увеличение итераций, ожидаемо, приводит к снижению скорости обработки и оправдано только при высоком уровне шума, который, как правило, свидетельствует о нештатной работе телевизионного оборудования.

Скорость вычисления у разработанного алгоритма при четырех итерациях в 25 раз выше, чем у алгоритма полного перебора.

Таким образом, разработанный алгоритм характеризуется погрешностью оценки параметров, как в алгоритме полного перебора, но обладает более высоким быстродействием.

Заключение

Разработанный алгоритм был апробирован для поиска изображений в потоковом видео в видеосистеме вагона-путеизмерителя в рамках задачи индексации («быстрого поиска»).

Разработанный алгоритм может быть использован в медицинских целях, в аэрокосмической отрасли, охранных системах и для расширения функциональных возможностей многокамерных систем технического зрения.

Разработанный алгоритм обеспечивает высокую скорость обработки видеоданных и может быть использован в задачах совмещения изображений, ориентированных на работу в режиме реального времени.

About the authors

Rinat R. Diyazitdinov

Author for correspondence.
Email: rinat.diyazitdinov@gmail.com
ORCID iD: 0000-0001-6360-0351

Candidate of Technical Sciences, associate professor of the Department of Networks and Communication Systems

References

Supplementary files

Supplementary Files

Action

1. JATS XML

Download

2. Fig. 1. Initial images: f(x, y) (a), g(x, y) (b)

Download (184KB)

Indexing metadata

3. Fig. 2. Evaluation of scale and rotation

Download (199KB)

Indexing metadata

4. Fig. 3. Error in parameter estimation

Download (1MB)

Indexing metadata