Как нормализовать эллиптические коэффициенты Фурье?

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

Замкнутая плоская ломаная p описывается набором точек m_points: m_points[i][0] сохраняет координаты xi, m_points i сохраняет координаты yi. Я начинаю с формы 0 до m_num_points-1.

EFC полилинии рассчитывается по этому алгоритму (коэффициенты находятся в массиве EFD).

// calc accumulated curve length, delta x, and delta y
dt[0] := 0;
dx[0] := 0;
dy[0] := 0;
tp[0] := 0;
T := 0;
for i:=1 to p.m_num_points-1 do begin
    va := VectorAffineSubtract(p.m_points[i], p.m_points[i-1]);
    dt[i] := VectorLength( va );
    dx[i] := p.m_points[i][0] - p.m_points[i-1][0];
    dy[i] := p.m_points[i][1] - p.m_points[i-1][1];
    tp[i] := tp[i-1] + dt[i];
    T := tp[i];
end;

Tpi := T / (2*PI*PI);
piT := PI*2 / T;

// find elliptic fourier coefficients
// first
An  := 0;
Cn  := 0;
for i:=0 to p.m_num_points-1 do begin
    An := An + (p.m_points[i][0] + p.m_points[i-1][0]) * dt[i] / 2;
    Cn  := Cn + (p.m_points[i][1] + p.m_points[i-1][1]) * dt[i] / 2;
end;
EFD[0][0] := An / T;
EFD[0][1] := 0;
EFD[0][2] := Cn / T;
EFD[0][3] := 0;
// next
for n:=1 to m_num_efd do begin
    Tn  := Tpi / (n*n);
    piTn  := piT * n;
    An  := 0;
    Bn  := 0;
    Cn  := 0;
    Dn  := 0;
    for i:=1 to p.m_num_points-1 do begin
        An := An + dx[i]/dt[i]*( cos(piTn*tp[i]) - cos(piTn*tp[i-1]) );
        Bn := Bn + dx[i]/dt[i]*( sin(piTn*tp[i]) - sin(piTn*tp[i-1]) );
        Cn := Cn + dy[i]/dt[i]*( cos(piTn*tp[i]) - cos(piTn*tp[i-1]) );
        Dn := Dn + dy[i]/dt[i]*( sin(piTn*tp[i]) - sin(piTn*tp[i-1]) );
    end;
    EFD[n][0] := An * Tn;
    EFD[n][1] := Bn * Tn;
    EFD[n][2] := Cn * Tn;
    EFD[n][3] := Dn * Tn;
end;

Для восстановления полилинии из набора гармоник я использую этот алгоритм.

// restore outline
resP := TYPolyline.create();
for i:=0 to p.m_num_points-1 do begin
    xn := EFD[0][0];
    yn := EFD[0][2];
    for n:=1 to m_num_efd do begin
        xn := xn + EFD[n][0] * cos(piT*n*tp[i]) + EFD[n][1] * sin(piT*n*tp[i]);
        yn := yn + EFD[n][2] * cos(piT*n*tp[i]) + EFD[n][3] * sin(piT*n*tp[i]);
    end;
    resP.add_point( xn, yn );
end;

Теперь мне нужно нормализовать EFD и поместить их в новый массив NEFD. Я делаю так:

// Normalize EFD
An := EFD[0][0];
Bn := EFD[0][1];
Cn := EFD[0][2];
Dn := EFD[0][3];

for n:=0 to m_num_efd do begin
    NEFD[n][0]  := EFD[n][0];
    NEFD[n][1]  := EFD[n][1];
    NEFD[n][2]  := EFD[n][2];
    NEFD[n][3]  := EFD[n][3];
end;

// rotate starting point
angl_o := 0.5 * ArcTan2( 2*(An*Bn + Cn*Dn), An*An + Cn*Cn - Bn*Bn - Dn*Dn );
for n:=1 to m_num_efd+1 do begin
    NEFD[n-1][0]:= EFD[n-1][0]*cos((n)*angl_o) + EFD[n-1][1]*sin((n)*angl_o);
    NEFD[n-1][1]:=-EFD[n-1][0]*sin((n)*angl_o) + EFD[n-1][1]*cos((n)*angl_o);

    NEFD[n-1][2]:= EFD[n-1][2]*cos((n)*angl_o) + EFD[n-1][3]*sin((n)*angl_o);
    NEFD[n-1][3]:=-EFD[n-1][3]*sin((n)*angl_o) + EFD[n-1][3]*cos((n)*angl_o);
end;

// make the semi-major axis parallel to the x-axis
angl_w := ArcTan( NEFD[1][2] / NEFD[1][0] );
cs   := cos(angl_w);
ss   := sin(angl_w);
for n:=1 to m_num_efd do begin
    NEFD[n][0]  := cs*NEFD[n][0] + ss*NEFD[n][2];
    NEFD[n][1]  := cs*NEFD[n][1] + ss*NEFD[n][3];

    NEFD[n][2]  :=-ss*NEFD[n][0] + cs*NEFD[n][2];
    NEFD[n][3]  :=-ss*NEFD[n][1] + cs*NEFD[n][3];
end;

// size invariant
R   := sqrt( NEFD[1][0]*NEFD[1][0] );
for n:=0 to m_num_efd do begin
    NEFD[n][0]  := NEFD[n][0] / R;
    NEFD[n][1]  := NEFD[n][1] / R;
    NEFD[n][2]  := NEFD[n][2] / R;
    NEFD[n][3]  := NEFD[n][3] / R;
end;

Когда я пытаюсь сделать большие полуоси параллельными X и поворачивать их, я получаю ужасный результат. Похоже, что восстановленная форма вращается вокруг оси Z. (Исходная форма слева, восстановлена ​​справа.)

Что не так с моим кодом?

ОБНОВИТЬ

После успешного ответа @MBo необходимы следующие модификации кода:

// MODIFIED: change n-1 to n
// rotate starting point
angl_o := 0.5 * ArcTan2( 2*(An*Bn + Cn*Dn), An*An + Cn*Cn - Bn*Bn - Dn*Dn );
for n:=1 to m_num_efd+1 do begin
    NEFD[n][0]:= EFD[n][0]*cos((n)*angl_o) + EFD[n][1]*sin((n)*angl_o);
    NEFD[n][1]:=-EFD[n][0]*sin((n)*angl_o) + EFD[n][1]*cos((n)*angl_o);

    NEFD[n][2]:= EFD[n][2]*cos((n)*angl_o) + EFD[n][3]*sin((n)*angl_o);
    NEFD[n][3]:=-EFD[n][3]*sin((n)*angl_o) + EFD[n][3]*cos((n)*angl_o);
end;

// MODIFIED: change left NEFD to EFD
// make the semi-major axis parallel to the x-axis
angl_w := ArcTan( NEFD[1][2] / NEFD[1][0] );
cs   := cos(angl_w);
ss   := sin(angl_w);
for n:=1 to m_num_efd do begin
    EFD[n][0]  := cs*NEFD[n][0] + ss*NEFD[n][2];
    EFD[n][1]  := cs*NEFD[n][1] + ss*NEFD[n][3];

    EFD[n][2]  :=-ss*NEFD[n][0] + cs*NEFD[n][2];
    EFD[n][3]  :=-ss*NEFD[n][1] + cs*NEFD[n][3];
end;

// ADDED: place normalized EFD into NEFD
for n:=0 to m_num_efd do begin
    NEFD[n][0]  := EFD[n][0];
    NEFD[n][1]  := EFD[n][1];
    NEFD[n][2]  := EFD[n][2];
    NEFD[n][3]  := EFD[n][3];
end;
Источник

1 ответ

Решение

Возможная причина:
В этом цикле вы используете измененное значение NEFD[n][0] рассчитать новое значение NEFD[n][2] (то же самое для NEFD[n][1]). Кажется, вы должны хранить и использовать неизмененные значения.

for n:=1 to m_num_efd do begin
    NEFD[n][0]  := cs*NEFD[n][0] + ss*NEFD[n][2];
    NEFD[n][1]  := cs*NEFD[n][1] + ss*NEFD[n][3];

                      vvvvvvvvv
    NEFD[n][2]  :=-ss*NEFD[n][0] + cs*NEFD[n][2];

                      vvvvvvvvv   
    NEFD[n][3]  :=-ss*NEFD[n][1] + cs*NEFD[n][3];
end;
Источник
Другие вопросы по тегам