algorithm - How to normalize elliptic fourier coefficients? -

i'm writing program find elliptic fourier coefficients (efc) of closed planar curve, , have trouble coefficients normalization.

the closed planar polyline p described set of points m_points: m_points[i][0] keeps xi-coordinates, m_pointsi keeps yi-coordinates. start form 0 m_num_points-1.

the efc of polyline calculated algorithm (coefficients in efd array).

// calc accumulated curve length, delta x, , delta y dt[0] := 0; dx[0] := 0; dy[0] := 0; tp[0] := 0; t := 0; i:=1 p.m_num_points-1 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  := 0; cn  := 0; i:=0 p.m_num_points-1 begin     := + (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] := / t; efd[0][1] := 0; efd[0][2] := cn / t; efd[0][3] := 0; // next n:=1 m_num_efd begin     tn  := tpi / (n*n);     pitn  := pit * n;      := 0;     bn  := 0;     cn  := 0;     dn  := 0;     i:=1 p.m_num_points-1 begin         := + 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] := * tn;     efd[n][1] := bn * tn;     efd[n][2] := cn * tn;     efd[n][3] := dn * tn; end; 

to restore polyline set of harmonics use algorithm.

// restore outline resp := typolyline.create(); i:=0 p.m_num_points-1 begin     xn := efd[0][0];     yn := efd[0][2];     n:=1 m_num_efd 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; 

now need normalize efd , place them in new array nefd. so:

// normalize efd := efd[0][0]; bn := efd[0][1]; cn := efd[0][2]; dn := efd[0][3];  n:=0 m_num_efd 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 ); n:=1 m_num_efd+1 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 semi-major axis parallel x-axis angl_w := arctan( nefd[1][2] / nefd[1][0] ); cs   := cos(angl_w); ss   := sin(angl_w); n:=1 m_num_efd 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] ); n:=0 m_num_efd 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; 

when try make semi-major axes parallel x , rotate coefficeints, ugly result. looks restored shape rotated around z-axis. (source shape on left, restored on right.)

what wrong code?

update

after successful answer of @mbo following modifications of code necessary:

// modified: change n-1 n // rotate starting point angl_o := 0.5 * arctan2( 2*(an*bn + cn*dn), an*an + cn*cn - bn*bn - dn*dn ); n:=1 m_num_efd+1 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 efd // make semi-major axis parallel x-axis angl_w := arctan( nefd[1][2] / nefd[1][0] ); cs   := cos(angl_w); ss   := sin(angl_w); n:=1 m_num_efd 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 nefd n:=0 m_num_efd 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; 

possible reason:
in cycle use modified value of nefd[n][0] calculate new value of nefd[n][2] (the same nefd[n][1]). seems have store , use unmodified values.

for n:=1 m_num_efd 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;