Субпиксельная корреляция

Материал из Техническое зрение
Перейти к: навигация, поиск

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

4-2-16.jpg

Обозначение систем координат в методе субпиксельной корреляции

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

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

Обозначим через $f(x,y)$ распределение яркостей на участке левого изображения, который будет эталоном. Начало прямоугольной системы координат $\langle x,y\rangle$ поместим в середину центрального пиксела эталона. Обозначим через $g(x_1 ,y_1 )$ распределение яркостей на участке правого изображения, который соответствует эталону. Форма этого участка отличается от формы эталона по причине перспективных искажений наблюдаемой сцены. Начало системы координат $\langle x_1 ,y_1 \rangle$ поместим в центре участка правого изображения (рис. 16).


Координаты $\langle x,y\rangle$ и $\langle x_1 ,y_1 \rangle$ связаны неизвестным преобразованием (13), где ${\rm {\textbf p}}$ - вектор параметров преобразования (например аффинных или др.): $$ \begin{gather}\tag{1} {\begin{array}{*{20}c} {x_1 =x_1 (x,y,{\rm {\textbf p}}),} \\ {y_1 =y_1 (x,y,{\rm {\textbf p}})} \\ \end{array} }_{.} \end{gather} $$

Вектор параметров ${\rm {\textbf p}}$ найдем путем максимизации меры сходства изображений, за которую принимается функция нормализованной корреляции (14), $$ \begin{gather}\tag{2} K({\rm {\textbf p}})= \frac {\sum\limits_{\langle x,y\rangle} {f(x,y)g(x_1 ,y_1 )-N{\bar f}{\bar g}}} {{\big(\sum\limits_{\langle x,y\rangle} f^2(x,y)-N{\bar f}^2\big)}^{1/2} {\big(\sum\limits_{\langle x,y\rangle} g^2(x_1 , y_1 )-N {\bar g}^2\big)}^{1/2}}. \end{gather} $$

В этой формуле знак $\sum\limits_{\langle x,y\rangle} $ обозначает суммирование по всем пикселам эталона; $N$ - количество пикселов. Средняя яркость находится по формулам $$ \begin{gather}\tag{3} {\bar f}= \frac{1}{N} \sum\limits_{\langle x,y\rangle} f(x,y); \quad {\bar g}= \frac{1}{N} \sum\limits_{\langle x,y\rangle} {g(x_1 ,y_1 )}. \end{gather} $$

Теперь задачу стереоотождествления можно сформулировать следующим образом: найти вектор параметров преобразования из условия $$ \begin{gather}\tag{4} {\rm {\textbf p}}^{\ast} =\mathop {\arg \max }\limits_{\rm {\textbf p}} K({\rm {\textbf p}}). \end{gather} $$

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

В дальнейшем для упрощения будем предполагать, что среднее значение яркости эталона равно нулю - $\bar {f}=0$. При этом яркости всех пикселов эталона преобразуются следующим образом: $$ f(x,y)\to f(x,y)-\bar {f}. $$

Преобразование (13) учитывает изменение формы участка текущего изображения, соответствующего эталону. В качестве этого преобразования могут быть выбраны различные модели с большим количеством параметров. Однако следует учесть, что субпиксельная точность оценивания может быть достигнута и в случае простых моделей искажений, тогда как использование сложных и гибких моделей, например полиномиальных, приводит в последствии к большим вычислительным затратам и проблемам сходимости процедур оценки модельных параметров. Поэтому будем, как в использовать в качестве (13) аффинное преобразование вида $$ x_1 = a_1 + a_2 x + a_3 y, $$

$$ \begin{gather}\tag{5} \end{gather} $$

$$ y_1 = b_1 + b_2 x + b_3 y, $$ которое обладает достаточной гибкостью и, как правило, приводит к хорошим результатам.

Таким образом, для решения задачи (16) необходимо найти вектор параметров $$ \textbf{p}=[a_1 a_2 a_3 b_1 b_2 b_3 ]^{T}. $$ Примем начальное приближение вектора параметров в виде $[a^* 1 0 b^* 0 1]^{T}$, которое можно найти с помощью классического корреляционного метода. Обозначим $g^*(x,y)$ распределение яркости на участке правого изображения, положение которого задается начальным вектором параметров. Линеаризуем неизвестную функцию $g(x_1 ,y_1 )$ относительно $g^\ast (x,y)$ по параметрам преобразования (17). В форме приращений (17) имеет вид $$ { \Delta x_1 =\Delta a_1 +x\Delta a_2 +y \Delta a_3 ,} $$

$$ \begin{gather}\tag{6} \end{gather} $$

$$ { \Delta y_1 =\Delta b_1 +x\Delta b_2 +y \Delta b_3 . } $$ Линеаризация $g(x_1 ,y_1 )$ дает $$ \begin{gather}\tag{7} g(x_1 ,y_1 )\approx g^{\ast} (x,y)+\frac{\partial g^{\ast} (x,y)}{\partial x}\Delta x_1 +\frac{\partial g^{\ast} (x,y)}{\partial y}\Delta y_1 . \end{gather} $$ Для краткости в дальнейшем в выражениях вида $g^\ast (x,y)$, $f(x,y)$ и т.д. аргументы $(x,y)$ опускаются. Обозначим также $$ g_x^\ast =\frac{\partial g\ast } {\partial x},\quad {g_y^\ast } =\frac{\partial g\ast }{\partial y}. $$

Из (18) и (19) имеем $$ \begin{gather}\tag{8} g(x_1 ,y_1 )\approx g^\ast + g_x^\ast \Delta a_1 + g_x^\ast x \Delta a_2 + g_x^\ast y\Delta a_3 + g_y^\ast \Delta b_1 + g_y^\ast x \Delta b_2 + g_y^\ast y \Delta b =\textbf{\textrm{g}}^\textrm{T}\Delta \textbf{\textrm{p}}, \end{gather} $$ где $$ {\rm{\bf g}}^\textrm{T}= [{ {g^\ast } {g_x^\ast } {xg_x^\ast } {yg_x^\ast } {g_y^\ast } {xg_y^\ast } {yg_y^\ast }} ], $$

$$ \begin{gather}\tag{9} \Delta {\rm {\textbf p}}^\textrm{T}=[ 1 {\Delta a_1 } {\Delta a_2 } {\Delta a_3 } {\Delta b_1 } {\Delta b_2 } {\Delta b_3 }]. \end{gather} $$

Производные функции $g(x,y)$ вычисляются с помощью их дискретных аппроксимаций

$g_x (x,y)=0{,}5(g(x+1,y)-g(x-1,y))$ - для внутренних точек изображения;

$g_x (x,y)=g(x+1,y)-g(x,y)$ - на левом краю изображения;

$g_x (x,y)=g(x,y)-g(x-1,y)$ - на правом краю изображения.

Аналогичные формулы применяются для координаты $y$.

Выражение (15) для среднего значения яркости участка правого изображения имеет вид $$ \begin{gather}\tag{10} \bar {g}\approx \frac{1}{N}\sum\limits_{\langle x,y\rangle} {{\rm {\bf g}}^\textrm{T}\Delta {\rm {\textbf p}}={\rm {\bf \bar {g}}}^\textrm{T}\Delta {\rm {\textbf p}}}. \end{gather} $$

Подставим (21), (22) в выражение для коэффициента корреляции (14): $$ \begin{gather}\tag{11} K(\Delta {\rm\bf p})= \frac { \sum\limits_{\langle x,y\rangle} f{\rm\bf g}^{\textrm{T}} \Delta {\rm\bf p} } { \big({\sum\limits_{\langle x,y\rangle} f^2}\big)^{1/2} \big({\sum\limits_{\langle x,y\rangle} \Delta {\rm\bf p}^{\textrm{T}} {\rm\bf g}{\rm\bf g}^{\textrm{T}} \Delta {\rm\bf p} - N \Delta {\rm\bf p}^{\textrm{T}} {\rm\bf \bar g} {\rm\bf \bar g} ^{\textrm{T}} \Delta {\rm\bf p}}\big)^{1/2} } . \end{gather} $$

Тогда исходная задача (16) сведется к задаче

$$ \begin{gather}\tag{12} \Delta {\rm {\textbf p}}^\ast =\arg \max K(\Delta {\rm {\textbf p}}). \end{gather} $$

Дисперсия яркости эталона $$ D_f =\frac{1}{N}\sum\limits_{\langle x,y\rangle} {f^2} $$ является константой. Поэтому задача (23) после равносильных преобразований $$ {K}'(\Delta {\rm {\textbf p}})=(ND_f K(\Delta {\rm {\textbf p}}))^2 $$ имеет вид $$ \begin{gather}\tag{13} K' (\Delta \textbf{p})= \frac {\Delta \textbf{p}^{\textrm{T}} \big({\sum\limits_{\langle x,y\rangle} f\textbf{g} }\big) \big({\sum\limits_{\langle x,y\rangle} f\textbf{g}^{\textrm{T}} }\big)\Delta \textbf{p}} {\Delta \textbf{p}^{\textrm{T}} \big({\sum\limits_{\langle x,y\rangle} \textbf{gg}^{\textrm{T}} - N {\bar {\textbf{g}}}{\bar {\textbf{g}}}^{\textrm{T}} }\big)\Delta \textbf{p}} =\frac {\Delta \textbf{p}^{\textrm{T}} \textbf{A} \Delta \textbf{p}} {\Delta \textbf{p}^{\textrm{T}} \textbf{B} \Delta \textbf{p}}. \end{gather} $$

Окончательно, задача (16) принимает вид $$ \begin{gather}\tag{14} \Delta {\rm {\textbf p}}=\arg \max K'(\Delta {\rm {\textbf p}}), \end{gather} $$ где ${\rm {\textbf B}}$ - матрица размером 7$\times $7, ${\rm {\bf r}}=\sum\limits_{\langle x,y\rangle} {f{\rm {\bf g}}^\textrm{T}} $ - вектор размерности 7, ${\rm {\textbf B}}=\sum\limits_{\langle x,y\rangle} {{\rm {\bf gg}}^\textrm{T}-N{\rm {\bf \bar {g}}}{\rm {\bf \bar {g}}}^\textrm{T}}$ - матрица размером 7$\times $7.

Матрица \textbf{B} - симметрическая и положительно определенная. Последнее следует из того, что знаменатель формулы (24) есть величина, пропорциональная дисперсии яркости участка правого изображения. Для реальных изображений матрицу \textbf{B} можно считать невырожденной.

Из матричной алгебры известно что задача максимизации отношения квадратичных форм $$ \begin{gather}\tag{15} \lambda =\frac{\rm{\mathbf{\textbf{x}^\textrm{T}A\textbf{x}}}}{\rm{\mathbf{\textbf{x}^\textrm{T}B\textbf{x}}}}\to \mathop {\max }\limits_\textbf{x} , \end{gather} $$ где \textbf{x} - вектор неизвестных параметров, сводится к эквивалентной задаче на обобщенные собственные значения $$ \begin{gather}\tag{16} {\rm {\bf A\textbf{x}}}=\lambda {\rm {\bf B\textbf{x}}}. \end{gather} $$

Это легко показать, векторно продифференцировав по ${\rm {\bf x}}$\textbf{ }обе части выражения (26): $$ \frac{\partial}{\partial{\rm {\bf x}}}\left(\frac{{\rm {\bf x}}^\textrm{T}{\rm {\bf Ax}}}{{\rm {\bf x}}^\textrm{T}{\rm {\textbf B}}{\rm {\bf x}}}\right)=\frac{({\rm {\bf x}}^\textrm{T}{\rm {\textbf B}}{\rm {\bf x}}){\rm {\bf Ax}}-({\rm {\bf x}}^\textrm{T}{\rm {\bf Ax}}){\rm {\textbf B}}{\rm {\bf x}}}{({\rm {\bf x}}^\textrm{T}{\rm {\bf Bx}})^2}=\frac{{\rm {\bf Ax}}-\lambda {\rm {\textbf B}}{\rm {\bf x}}}{{\rm {\bf x}}^\textrm{T}{\rm {\textbf B}}{\rm {\bf x}}}=0. $$ Таким образом, решением задачи (26) будет собственный вектор, отвечающий максимальному собственному значению в задаче (27).

Рассмотрим теперь решение эквивалентной задачи методом Холецкого [166].

В основе такого решения лежит вспомогательная лемма, подробное доказательство которой приводится в [263].

Из доказательства леммы следует эффективный метод нахождения собственного вектора \textbf{x}, отвечающего единственному ненулевому собственному значению задачи (27), а именно, вектор \textbf{х} можно найти как решение системы линейных уравнений $$ {\rm {\textbf B}}^{\rm {\bf x}}={\rm {\bf r}}, $$ которая может быть решена методом Холецкого, требующим примерно $n^3/3+2n^2$ операций.

Последовательность действий при использовании этого метода состоит в следующем.


  1. Получаем разложение матрицы ${\rm {\textbf B}}={\rm {\bf LL}}^\textrm{T}$ .

Обозначим ${\rm\bf{B}}= \Vert b_{ij} \Vert$, ${\rm\bf L}=\Vert l_{ij}\Vert$.

Тогда $$ l_{kk} =(b_{kk} -\sum\limits_{p=1}^{k-1} {b_{kp}^2 )^{1/2}} ,\quad k=1, \ldots , n, $$ $$ l_{ik} =\frac{b_{ik} -\sum\limits_{p=1}^{k-1} {b_{ip} b_{kp} } }{l_{kk} }, \quad i=k+1,\ldots ,n . $$ Этот шаг требует примерно $n^3/3$ операций.

  1. Решаем систему ${\rm {\bf Bx}}={\rm {\bf r}}$, т.е.

последовательно решаем две системы уравнений с треугольными матрицами: $$ {\rm {\bf Ly}}={\rm {\bf r}}, \quad {\rm {\bf L}}^{\textrm{T}}{\rm {\bf x}}={\rm {\bf y}}. $$ Этот шаг требует примерно $2n^2$ операций.


Решением задачи (25) является любой вектор вида $\alpha \Delta {\rm {\textbf p}}$, где $\Delta {\rm {\textbf p}}$ - найденный вектор поправок, $\alpha $ - константа, которая управляет сходимостью метода. Поэтому является естественным выбирать такое $\alpha $, которое при данном векторе $\Delta {\rm {\textbf p}}$ максимизирует коэффициент корреляции (14), т.е. задача принимает вид

$$ \begin{gather}\tag{17} K(\alpha \Delta {\rm {\textbf p}})\to \mathop {\max }\limits_\alpha . \end{gather} $$

Умножение вектора параметров на константу соответствует изменению масштаба исходного преобразования (13). В качестве преобразования (13) используется следующее, зависящее от одного параметра $\alpha $: $$ x_2 =\alpha x_1 (x,y,{\rm {\textbf p}}), \quad y_2 =\alpha y_1 (x,y,{\rm {\textbf p}}), $$ которое в форме приращений имеет вид $$ \Delta x_2 =x_1 \Delta \alpha , \quad y_2 =y_1 \Delta \alpha . $$

Обратимся вновь к линеаризованному выражению для ${\rm {\bf g}}(x_2 ,y_2 )$ (20): $$ {\rm\bf g}(x_2 ,y_2 )\approx {\rm\bf g}(x_1 ,y_1 )+({\rm \bf g}_x (x_1 ,y_1 )x_1 + {\rm \bf g}_y (x_1 ,y_1 ) y_1 )\Delta \alpha = \textrm{\textbf{g}}_\alpha^{\textrm{T}} \alpha , $$ где $$ \boldsymbol\alpha = \left[{ \begin{matrix} 1 \cr \Delta \alpha \end{matrix} }\right] , \quad {\rm\bf g}_a = \left[{ \begin{matrix} {\rm\bf g}(x_1 ,y_1 ) \cr {\rm\bf g}_x x_1 +{\rm\bf g}_y y_1 \end{matrix} }\right] . $$

Преобразованная задача максимизации коэффициента корреляции имеет вид, аналогичный (25): $$ \begin{gather}\tag{18} K''\left( {\Delta \alpha } \right)= \frac {\boldsymbol\alpha ^{\textrm {т}} \left( \sum\limits_{(x,y)}f{\rm {\bf g}}_\alpha \right) \left(\sum\limits_{(x,y)} f{\rm {\bf g}}_\alpha ^{\textrm {т}}\right)\boldsymbol\alpha } {\boldsymbol\alpha^{\textrm {т}} \left(\sum\limits_{(x,y)} {\rm{\bf g}}_\alpha{\rm{\bf g}}_\alpha^{\textrm {т}} - N{\rm{\bf \bar{g}}}_\alpha {\rm{\bf \bar{g}}}_\alpha^{\textrm {т}} \right)\boldsymbol\alpha } =\frac {\boldsymbol\alpha ^{\textrm {т}}A_\alpha \boldsymbol\alpha } {\boldsymbol\alpha ^{\textrm {т}}B_\alpha \boldsymbol\alpha }, \end{gather} $$ $\Delta \alpha =\arg \max K^{\prime\prime}(\Delta \alpha )$, где $$ {\rm\bf A}_{\rm\bf\alpha} ={\rm\bf r}_{\rm\bf\alpha} {\rm\bf r}_{\rm\bf\alpha}^{\rm T}, \quad \textrm{\textbf{r}}_\alpha =\sum\limits_{(x,y)} {f{\rm\bf g}}_\alpha ^{\rm T}= \left[{ \begin{matrix} {r_1 } \cr {r_2 } \end{matrix} }\right] , \quad \rm{\bf{B}}_{\alpha} = \sum\limits_{(x,y)} {\rm\bf g}_{\alpha} {\rm\bf g}_{\alpha}^{\rm T} - N {\bar {\rm\bf g}}_{\alpha} {\bar {\rm\bf g} }_{\alpha}^{\rm T} = \left[{ \begin{matrix} {b_{11} } & {b_{12} } \cr {b_{12} } & {b_{22} } \end{matrix} }\right] . $$

По лемме решение (28) имеет вид

$$ \begin{gather}\tag{19} \alpha =\left[{ \begin{matrix} {\alpha _1 } \cr {\alpha _2 } \end{matrix} }\right] ={\rm{\textbf B}}_\alpha ^{-1} {\rm{\bf r}}_\alpha =\frac{1}{b_{11} b_{22} -b_{12}^2 } \left[{ \begin{matrix} {r_1 b_{22} -r_2 b_{12} } \cr {r_2 b_{11} -r_1 b_{12} } \end{matrix} }\right]. \end{gather} $$

Таким образом, решение задачи нахождения максимума коэффициента корреляции сводится к двум шагам:


  1. нахождение вектора параметров $\Delta {\rm {\textbf p}}$ аффинного преобразования (17), определяющего направление в 6-мерном пространстве параметров;
  2. нахождение оптимального шага $\alpha \Delta {\rm {\textbf p}}$ (30) вдоль этого вектора.


Коэффициент корреляции, соответствующий найденному вектору поправок параметров $\Delta {\rm {\textbf p}}$, находится по формуле (29): $$ K\left( {\Delta \textrm{\textbf{p}}} \right)=\sqrt {\frac{\textrm{\textbf{r}}^\textrm{T}\Delta \textrm{\textbf{p}}}{ND_f }} . $$ Аналогично, коэффициент корреляции, соответствующий найденному вектору $\alpha $, находится по формуле $$ K(\Delta \textbf{\textrm{p}})=\sqrt {\frac{r_1 \alpha _{{\kern 1pt}1} +r_2 \alpha _2 }{ND_f }}. $$

Поскольку исходная задача (14) была нелинейной, ее финальное решение находится последовательными итерациями. На каждой итерации решается линеаризованная задача и определяется очередное приближение положения максимума коэффициента корреляции. Начальным приближением является функция $g^\ast $, определяемая начальным вектором параметров, который находится одним из методов грубого стереоотождествления. На $m$-й итерации сначала находится вектор поправок аффинного преобразования $\Delta \alpha \Delta {\rm {\textbf p}}_{(m)} $ и уточняется само преобразование (из (18)): $$ \begin{array}{c} x_1 =x+\Delta x_1 ,\\ y_1 =y+\Delta y_1 . \end{array} $$

Это преобразование дает смещение центра участка правого изображения. На следующей итерации это смещение является новым началом прямоугольной системы координат $\langle x, y\rangle$. Так как новая система координат может быть сдвинута на нецелое число пикселов, то для вычисления производных при линеаризации яркость в точках с нецелыми координатами получается при помощи билинейной интерполяции по значениям яркости в соседних точках с целыми координатами.

Процесс прекращается при достижении заданного количества итераций или когда абсолютные значения поправок к параметрам сдвига становятся меньше заданных констант: $$ \begin{array}{c} \left| {\Delta a_{1(m)} -\Delta a_{1(m-1)} } \right|<C_1 ,\\ %\quad \left| {\Delta b_{1(m)} -\Delta b_{1(m-1)} } \right|<C_2 . \end{array} $$ Результаты использования метода субпиксельной корреляции иллюстрируются многочисленными вычислительными экспериментами (рис. 17 - 19).

На рис. 17 - 18 показаны траектории сходимости восьми начальных приближений. Показаны два случая: близкое $(d=2)$ и далекое $(d=4)$ начальное приближение. На рис. 19 показан пример ложной сходимости некоторых траекторий. Толстым штрихом показан пиксел, найденный классическим корреляционным алгоритмом, от которого отсчитывались начальные приближения. Вид этих траекторий наглядно демонстрирует более быструю сходимость метода с выбором шага оптимизации (рис. 18, 19$\textit{б}$).

Метод субпиксельной корреляции может также применяться с учетом предварительной сегментации изображений, в частности, в модели стереоотождествления.

4-2-17.jpg

Траектории сходимости различных начальных приближений $d=2$

4-2-18.jpg

Траектории сходимости различных начальных приближений $d=4$

4-2-19.jpg

Случай ложной сходимости некоторых траекторий $d=5$

Полезные ссылки

  1. ☝ К началу
  2. ☜ Сравнение и привязка изображений. Стереоотождествление
Личные инструменты
Пространства имён

Варианты
Действия
Навигация
Инструменты