Wie normalisieren elliptische Fourier-Koeffizienten?

Question

Ich schreibe ein Programm, um einen Elliptischen Fourier Koeffizienten (EFC) einer geschlossenen planaren Kurve zu finden, und habe Probleme mit der Koeffizientennormierung.

Die geschlossene planare Polylinie p wird durch die Menge der Punkte m_points beschrieben: m_points [i] [0] behält xi-Koordinaten, m_points i behält yi-Koordinaten. Ich beginne Form 0 bis m_num_points-1.

Der EFC der Polylinie wird mit diesem Algorithmus berechnet (Koeffizienten sind im EFD-Array).

// 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; 

Um Polylinie aus dem Satz von Oberschwingungen wiederherzustellen, verwende ich diesen Algorithmus.

// 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; 

Jetzt muss ich EFD normalisieren und sie in einem neuen Array NEFD platzieren. Ich mache das so:

// 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; 

Wenn ich versuche, große Halbachsen X parallel zu machen und drehen coefficeints, erhalte ich ein hässliches Ergebnis. Es sieht so aus, als ob die wiederhergestellte Form um die z-Achse gedreht wird. (Quelle Form auf der linken Seite, auf der rechten Seite wieder hergestellt.)

Was mit meinem Code falsch?

UPDATE

Nach der erfolgreichen Antwort von @MBo folgende Änderungen des Codes notwendig:

// 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; 

Answer 1

Möglicher Grund:
In diesem Zyklus Sie modifizierte Wert von NEFD[n][0] verwenden, um Berechnen Sie den neuen Wert NEFD[n][2] (das gleiche für NEFD[n][1]). Es scheint, dass Sie unmodifizierte Werte speichern und verwenden müssen.

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;