Ненормализация коэффициентов эллипса после прямого выравнивания эллипса

Question

Я пытаюсь понять шаги нормализации и «ненормализации» в алгоритме прямого эллипсового метода наименьших квадратов, разработанном Фицгиббоном, Пилу и Фишером (улучшенный Галиром и Флюссером).

EDITED: Более подробная информация о теории добавлена. Является ли проблема с собственными значениями, из-за которой возникает путаница?

Краткая теория:

Эллипс представлен неявным полином второго порядка (общая коническая уравнения):

где:

Чтобы ограничить это общую конику к эллипсу, коэффициенты должны удовлетворять квадратичное ограничение:

, которое эквивалентно:

где С представляет собой матрицу нулей, за исключением :

Конструкция матрица D состоит из всех точек данных х юга я.

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

Обозначая:

Теперь у нас есть система:

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

Код:

фрагментов кода здесь непосредственно из кода MATLAB, представленной авторами: http://research.microsoft.com/en-us/um/people/awf/ellipse/fitellipse.html

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

% normalize data 
% X and Y are the vectors of data points, not normalized 
mx = mean(X); 
my = mean(Y); 
sx = (max(X)-min(X))/2; 
sy = (max(Y)-min(Y))/2; 

x = (X-mx)/sx; 
y = (Y-my)/sy; 

% Build design matrix 
D = [ x.*x x.*y y.*y x y ones(size(x)) ]; 

% Build scatter matrix 
S = D'*D; %' 

% Build 6x6 constraint matrix 
C(6,6) = 0; C(1,3) = -2; C(2,2) = 1; C(3,1) = -2; 

[gevec, geval] = eig(S,C); 

% Find the negative eigenvalue 
I = find(real(diag(geval)) < 1e-8 & ~isinf(diag(geval))); 

% Extract eigenvector corresponding to negative eigenvalue 
A = real(gevec(:,I)); 

После этого, нормализация восстанавливается на коэффициенты:

par = [ 
    A(1)*sy*sy, ... 
    A(2)*sx*sy, ... 
    A(3)*sx*sx, ... 
    -2*A(1)*sy*sy*mx - A(2)*sx*sy*my + A(4)*sx*sy*sy, ... 
    -A(2)*sx*sy*mx - 2*A(3)*sx*sx*my + A(5)*sx*sx*sy, ... 
    A(1)*sy*sy*mx*mx + A(2)*sx*sy*mx*my + A(3)*sx*sx*my*my ... 
    - A(4)*sx*sy*sy*mx - A(5)*sx*sx*sy*my ... 
    + A(6)*sx*sx*sy*sy ... 
    ]'; 

На данный момент, я не знаю, что случилось. Почему ненормализация последних трех коэффициентов A (d, e, f) зависит от первых трех коэффициентов? Как вы математически показываете, откуда берутся эти уравнения ненормализации?

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

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

Любая помощь приветствуется. Благодаря!

Answer 1

Во-первых, позвольте мне формализовать проблему в однородном пространстве (как и в Ричард Хартли и книги Эндрю Зиссерман в Multiple View Geometry):
Предположим, что

P=[X,Y,1]' 

наша точка в ненормированного пространстве и

p=lambda*[x,y,1]' 

наша точка в нормализованном пространстве, где лямбда является неважной свободной шкалой (в однородном пространстве [х, у, 1] = [10 * х, 10 * у, 10] и так далее).

Теперь ясно, что мы можем написать

x = (X-mx)/sx; 
y = (Y-my)/sy; 

в качестве простого матричного уравнения вида:

p=H*P; %(equation (1)) 

где

H=[1/sx, 0,  -mx/sx; 
    0,  1/sy, -my/sy; 
    0,  0,   1]; 

Кроме того, мы знаем, что эллипс с уравнением

A(1)*x^2 + A(2)*xy + A(3)*y^2 + A(4)*x + A(5)*y + A(6) = 0 %(first representation) 

можно записать в матричной форме как:

p'*C*p=0 %you can easily verify this by matrix multiplication 

где

C=[A(1), A(2)/2, A(4)/2; 
    A(2)/2, A(3), A(5)/2; 
    A(4)/2, A(5)/2, A(6)]; %second representation 

и

p=[x,y,1] 

, и ясно, что эти два представления эллипса точно так же, и эквивалент.

Также известно, что вектор A = [A (1), A (2), A (3), A (4), A (5), A (6)] является представлением типа 1 эллипс в нормализованном пространстве.
Таким образом, мы можем написать:

p'*C*p=0 

где р нормированный точки В и С, как определено ранее. Теперь мы можем использовать "уравнение (1): р = HP", чтобы получить хороший результат:

(H*P)'*C*(H*P)=0 

=====>

P'*H'*C*H*P=0 

=====>

P'*(H'*C*H)*P=0 

=====>

P'*(C1)*P=0 %(equation (2)) 

Мы видим, что Уравнение (2) является уравнением эллипса в ненормированного пространстве, где C1 вл етс представление 2-го типа эллипса, и мы знаем, что:

C1=H'*C*H 

Ans Кроме того, так как уравнение (2) является уравнением нулю, получаем может умножить его на любое ненулевое число.Таким образом, мы умножаем его на SX^2 * сы^2 и мы можем написать:

C1=sx^2*sy^2*H'*C*H 

И, наконец, мы получим результат

C1=[               A(1)*sy^2,               (A(2)*sx*sy)/2,                                        (A(4)*sx*sy^2)/2 - A(1)*mx*sy^2 - (A(2)*my*sx*sy)/2; 
                 (A(2)*sx*sy)/2,                A(3)*sx^2,                                        (A(5)*sx^2*sy)/2 - A(3)*my*sx^2 - (A(2)*mx*sx*sy)/2; 
    -(- (A(4)*sx^2*sy^2)/2 + (A(2)*my*sx^2*sy)/2 + A(1)*mx*sx*sy^2)/sx,  -(- (A(5)*sx^2*sy^2)/2 + A(3)*my*sx^2*sy + (A(2)*mx*sx*sy^2)/2)/sy,  (mx*(- (A(4)*sx^2*sy^2)/2 + (A(2)*my*sx^2*sy)/2 + A(1)*mx*sx*sy^2))/sx + (my*(- (A(5)*sx^2*sy^2)/2 + A(3)*my*sx^2*sy + (A(2)*mx*sx*sy^2)/2))/sy + A(6)*sx^2*sy^2 - (A(4)*mx*sx*sy^2)/2 - (A(5)*my*sx^2*sy)/2] 

, который может быть преобразован в тип-2 эллипса и получить точный результат, который мы искали:

[ A(1)*sy^2, A(2)*sx*sy, A(3)*sx^2, A(4)*sx*sy^2 - 2*A(1)*mx*sy^2 - A(2)*my*sx*sy, A(5)*sx^2*sy - 2*A(3)*my*sx^2 - A(2)*mx*sx*sy, A(2)*mx*my*sx*sy + A(1)*mx*my*sy^2 + A(3)*my^2*sx^2 + A(6)*sx^2*sy^2 - A(4)*mx*sx*sy^2 - A(5)*my*sx^2*sy] 

Если вам интересно, как мне удалось caculate этих трудоемких уравнений я могу дать вам MatLab код, чтобы сделать это для вас следующим образом:

syms sx sy mx my 
syms a b c d e f 
C=[a, b/2, d/2; 
    b/2, c, e/2; 
    d/2, e/2, f];  
H=[1/sx, 0, -mx/sx; 
    0,  1/sy, -my/sy; 
    0,  0,  1]; 
C1=sx^2*sy^2*H.'*C*H 
a=[Cp(1,1), 2*Cp(1,2), Cp(2,2), 2*Cp(1,3), 2*Cp(2,3), Cp(3,3)]