Определение аффинного преобразования методом наименьших квадратов

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

Аффинное преобразование задается одной из форм линейного преобразования: $$H' = AH+B$$ Здесь $H$ - вектор заданных точек, $H'$ - вектор уточненных, $A$ - матрица 3x3, $B$ - вектор-столбец.

Требуется найти такие коэффициенты матрицы $A$ и вектора $B$, при которых для всех пар точек выполняется наилучшее преобразование.

Задача характеризуется тем, что как исходные значения $H$, так и целевые значения $H'$ заданы с некоторыми погрешностями и, вообще говоря, не существует одного такого преобразования с матрицей $A$ и вектором $B$, которое бы в точности перевело все исходные точки в целевые уточненные. Какое бы преобразование ни использовалось, в итоге какие-то точки разойдутся или сильно, или очень сильно.

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

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

Далее нас интересуют значения коэффициентов матрицы $A$ и вектора $B$, при которых такая сумма минимальна. Если сумма квадратов отклонений минимальна, то набор соответствующих частных производных по коэффициентам $A_{ij}$ и $B_i$ должны быть равны нулям.

Положим, что пары точек заданы точками $h(x,y,x)$ и $h'(x',y',z')$ и эти значения $x$, $y$, $z$, $x'$, $y'$ и $z'$ заданы для всего набора, для всех пар.

Положим, что аффинное преобразование отображает: $$\left\{ \begin{array}{c} x' = a_{11}x + a_{12}y + a_{13}z + b_1 \\ y' = a_{21}x + a_{22}y + a_{23}z + b_2 \\ z' = a_{31}x + a_{32}y + a_{33}z + b_3 \end{array} \right.$$ Мы должны минимизировать значение $$\begin{array}{c} \sum\limits_n\left(a_{11}x + a_{12}y + a_{13}z + b_1 - x'\right)^2+\\ \sum\limits_n\left(a_{21}x + a_{22}y + a_{23}z + b_2 - y'\right)^2+\\ \sum\limits_n\left(a_{31}x + a_{32}y + a_{33}z + b_3 - z'\right)^2 = min\\ \end{array}$$ Для упрощения представим эту сумму в виде: $$\sum\limits_nS_x^2 + \sum\limits_nS_y^2 + \sum\limits_nS_z^2 = min$$ Берем частные производные по $a_{ij}$ и $b_i$ и приравниваем их нулям и используем промежуточное равенство: $$\frac{\partial}{\partial a_{ij}}\sum\limits_n S_k^2 = 2\sum\limits_n S_k \frac{\partial}{\partial a_{ij}}S_k$$ Таким образом, нужно опустить 2 и найти чему равна соответствующая частная производная по $a_{ij}$ и $b_i$. Для трехмерного пространства матрица $A$ содержит в общем случае $3\cdot 3=9$ коэффициентов и вектор $b$ содержит $3$ коэффициента. Итого мы должны получить систему уравнений с $9+3=12$ неизвестных, и решить её.

При взятии частных производных по коэффициентам $a_{ij}$ от величин $S_k$ задача существенно обьлегчается тем, что производная $S_k$ по $a_{ij}$ или $b_i$ либо равна 0 если $S_k$ не зависит от соответствующего коэффициента, либо равна одной из координат $x$, $y$, $z$, либо равна единице.

Беря частные производные по коэффициентам $a_{ij}$ и $b_i$, получаем уравнения: $$\frac{\partial}{\partial a_{11}}:\ \sum(a_{11}xx+a_{12}yx+a_{13}zx+b_1x)=\sum x'x$$ $$\frac{\partial}{\partial a_{12}}:\ \sum(a_{11}xy+a_{12}yy+a_{13}zy+b_1y)=\sum x'y$$ $$\frac{\partial}{\partial a_{13}}:\ \sum(a_{11}xz+a_{12}yz+a_{13}zz+b_1z)=\sum x'z$$ $$\frac{\partial}{\partial b_1}:\ \sum(a_{11}x+a_{12}y+a_{13}z+b_1)=\sum x'$$ $$\frac{\partial}{\partial a_{21}}:\ \sum(a_{21}xx+a_{22}yx+a_{23}zx+b_2x)=\sum y'x$$ $$\frac{\partial}{\partial a_{22}}:\ \sum(a_{21}xy+a_{22}yy+a_{23}zy+b_2y)=\sum y'y$$ $$\frac{\partial}{\partial a_{23}}:\ \sum(a_{21}xz+a_{22}yz+a_{23}zz+b_2z)=\sum y'z$$ $$\frac{\partial}{\partial b_2}:\ \sum(a_{21}x+a_{22}y+a_{23}z+b_2)=\sum y'$$ $$\frac{\partial}{\partial a_{31}}:\ \sum(a_{31}xx+a_{32}yx+a_{33}zx+b_3x)=\sum z'x$$ $$\frac{\partial}{\partial a_{32}}:\ \sum(a_{31}xy+a_{32}yy+a_{33}zy+b_3y)=\sum z'y$$ $$\frac{\partial}{\partial a_{33}}:\ \sum(a_{31}xz+a_{32}yz+a_{33}zz+b_3z)=\sum z'z$$ $$\frac{\partial}{\partial b_3}:\ \sum(a_{31}x+a_{32}y+a_{33}z+b_3)=\sum z'$$ Далее надо раскрыть все знаки сумм и получить систему уравнений относительно коэффициентов $a_{ij}$ и $b_i$. $$a_{11}\sum\limits_nxx+a_{12}\sum\limits_nyx+a_{13}\sum\limits_nzx+ b_1\sum\limits_nx=\sum\limits_nx'x$$ $$a_{11}\sum\limits_nxy+a_{12}\sum\limits_nyy+a_{13}\sum\limits_nzy+ b_1\sum\limits_ny=\sum\limits_nx'y$$ $$a_{11}\sum\limits_nxz+a_{12}\sum\limits_nyz+a_{13}\sum\limits_nzz+ b_1\sum\limits_nz=\sum\limits_nx'z$$ $$a_{11}\sum\limits_nx+a_{12}\sum\limits_ny+a_{13}\sum\limits_nz+ b_1n=\sum\limits_nx'$$ $$a_{21}\sum\limits_nxy+a_{22}\sum\limits_nyy+a_{23}\sum\limits_nzy+ b_2\sum\limits_ny=\sum\limits_ny'x$$ $$a_{21}\sum\limits_nxz+a_{22}\sum\limits_nyz+a_{23}\sum\limits_nzz+ b_2\sum\limits_nz=\sum\limits_ny'y$$ $$a_{21}\sum\limits_nxz+a_{22}\sum\limits_nyz+a_{23}\sum\limits_nzz+ b_2\sum\limits_nz=\sum\limits_ny'z$$ $$a_{21}\sum\limits_nx+a_{22}\sum\limits_ny+a_{23}\sum\limits_nz+ b_2n=\sum\limits_ny'$$ $$a_{31}\sum\limits_nxx+a_{32}\sum\limits_nyx+a_{33}\sum\limits_nzx+ b_3\sum\limits_nx=\sum\limits_nz'x$$ $$a_{31}\sum\limits_nxy+a_{32}\sum\limits_nyy+a_{33}\sum\limits_nzy+ b_3\sum\limits_ny=\sum\limits_nz'y$$ $$a_{31}\sum\limits_nxz+a_{32}\sum\limits_nyz+a_{33}\sum\limits_nzz+ b_3\sum\limits_nz=\sum\limits_nz'z$$ $$a_{31}\sum\limits_nx+a_{32}\sum\limits_ny+a_{33}\sum\limits_nz+ b_3n=\sum\limits_nz'$$ Здесь использовался тот факт, что $$\sum\limits_nb_3=b_3\sum\limits_n1=b_3n$$ Несмотря на то, что в каждом уравнении слева от знака $=$ стоит лишь 4 числа, совокупность этих уравнений задаёт систему 12x12 линейных уравнений для 12-ти коэффициентов, обозначенных через $a_{ij}$ и $b_i$.

Для приведения системы уравнений к СЛАУ нужно зафиксировать их порядок и в них найти коэффициенты матрицы, ну скажем $C$ и $D$, такие что $$CX=D$$ где $C$ - квадратная матрица, $D$ - вектор, $X$ - искомый вектор.

И решить это уравнение, найдя $X$.

Оставив уравнения в том же порядке, получим что соответственно $$D_1 = \sum x'x$$ $$D_2=\sum x'y$$ и так далее.

Упорядочив коэффициенты $a_{ij}$ и $b_i$ в определенном порядке и рассматривая их как элементы вектора $X$, выберем сопоставление: $$\begin{array}{cccc} X_1 = a_{11}, & X_2 = a_{12}, & X_3 = a_{13}, & X_4 = b_1, \\ X_5 = a_{21}, & X_6 = a_{22}, & X_7 = a_{23}, & X_8 = b_2, \\ X_9 = a_{31}, & X_{10}=a_{32}, & X_{11} = a_{33}, & X_{12} = b_3 \end{array}$$ После этого сопоставления можем определить по вышеприведенным уравнениям значения коэффициентов $C_{ij}$, например: $$C_{1,1} = \sum xx$$ $$C_{1,2} = \sum yx$$ $$C_{8,10} = 0$$ $$C_{4,1} = \sum x$$ $$C_{7,7} = \sum zz$$ $$C_{12,12} = n$$

Далее следует банальное программирование. Заводим матрицу 12x12 и вектор на 12. Всем коэффициентам присваиваем нулевые значения. В цикле по всем парам точек расписываем по вышеприведенным уравнениям значения $x$, $y$, $z$, $x'$, $y'$ и $z'$ или их произведения.

Решив систему уравнений (например, методом Гаусса), получаем в качестве результата вектор из искомых значений $a_{ij}$ и $b_i$.

Также можно отметить, что вышеприведенная система уравнений может решаться как система для 12-ти переменных, так и может быть разделена на 3 системы по 4 уравнения для вычисления векторов из 4-х коэффициентов: $$\left( \begin{array}{cccc} a_{11} & a_{12} & a_{13} & b_1 \end{array} \right)$$ $$\left( \begin{array}{cccc} a_{21} & a_{22} & a_{23} & b_2 \end{array} \right)$$ $$\left( \begin{array}{cccc} a_{31} & a_{32} & a_{33} & b_3 \end{array} \right)$$ Если расписать в виде матричного произведения эти отдельные 3 системы для 4-х неизвестных каждая, то получим: $$\left( \begin{array}{cccc} \sum xx & \sum yx & \sum zx & \sum x \\ \sum xy & \sum yy & \sum zy & \sum y \\ \sum xz & \sum yz & \sum zz & \sum z \\ \sum x & \sum y & \sum z & n \end{array} \right) \left( \begin{array}{c} a_{11} \\ a_{12} \\ a_{13} \\ b_1 \end{array} \right)= \left( \begin{array}{c} \sum x'x \\ \sum x'y \\ \sum x'z \\ \sum x' \end{array} \right)$$ $$\left( \begin{array}{cccc} \sum xx & \sum yx & \sum zx & \sum x \\ \sum xy & \sum yy & \sum zy & \sum y \\ \sum xz & \sum yz & \sum zz & \sum z \\ \sum x & \sum y & \sum z & n \end{array} \right) \left( \begin{array}{c} a_{21} \\ a_{22} \\ a_{23} \\ b_2 \end{array} \right)= \left( \begin{array}{c} \sum y'x \\ \sum y'y \\ \sum y'z \\ \sum y' \end{array} \right)$$ $$\left( \begin{array}{cccc} \sum xx & \sum yx & \sum zx & \sum x \\ \sum xy & \sum yy & \sum zy & \sum y \\ \sum xz & \sum yz & \sum zz & \sum z \\ \sum x & \sum y & \sum z & n \end{array} \right) \left( \begin{array}{c} a_{31} \\ a_{32} \\ a_{33} \\ b_3 \end{array} \right)= \left( \begin{array}{c} \sum z'x \\ \sum z'y \\ \sum z'z \\ \sum z' \end{array} \right)$$ Соответственно, если решать эти системы методом Гаусса, то надо три решения, по одной на каждую систему. А если решать методом нахождения обратной матрицы, то такую матрицу нужно искать лишь один раз в силу того что все три системы имеют одну и ту же матрицу.