私は閉曲線の楕円フーリエ係数(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に平行とcoefficeintsを回転させるようにしよう
は、私は醜い結果を得ます。復元された形状がz軸の周りを回転しているように見えます。 (右側に復元された左の光源形状、。)私のコードが間違っていますか?
UPDATE
@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 –