先日書いたアルゴリズムからいろいろ変更。まず、一般的なベクター描画ソフト、Illustrator とか Inkscape とかはベジエ曲線を使っているので、曲面生成にもベジエ曲線を使う。まず描画ソフトで一本だけベジエ曲線をイメージどおりに書いて、これを svg フォーマットで保存。このフォーマットは html みたいに可読性が高い。これを面生成コードで入力として読んでコード内部でベジエ曲線を再現できるようにする。
ベジエ曲線はy=f(x) という形ではなく x(t),y(t) という風に位置が与えられるので、Case 2 と Case 2 ' を区別しなくていい。Case 2 のみ考える。
//ベジエ曲線のクラス class B3 { public: //三次式の係数 double a0[2], a1[2], a2[2], a3[2]; #define LENDIV 1000 double len; double lbuf[LENDIV]; //制御点の座標で初期化 B3(double *r0, double *r1, double *r2, double *r3); ~B3(); // value void pos(double t, double *d); void dpos(double t, double *d); void ddpos(double t, double *d); }; B3::B3(double *r0, double *r1, double *r2, double *r3) { int i; //三次式の係数を計算 for(i=0;i<2;i++) { // it^3 r0 + 3t it^2 r1 + 3t^2 it r2 + t^3 r3 // = r0 | 1 -3t +3tt -ttt // +r1 3t (1-2t+tt) | 3t -6tt +3ttt // +r2 3tt (1-t) | 3tt -3ttt // +r3 ttt | ttt a3[i]= r3[i]-3*r2[i]+3*r1[i]-r0[i]; a2[i]= 3*r2[i] -6*r1[i]+3*r0[i]; a1[i]= 3*r1[i]-3*r0[i]; a0[i]= r0[i]; } double l=0; double dt=1.0/LENDIV; for(i=0;i<LENDIV;i++) { double t=dt*i; double v[2]; dpos(t, v); l += sqrt(v[0]*v[0]+v[1]*v[1])*dt; lbuf[i]=l; } len = l; } B3::~B3() { } void B3::pos(double t, double *d) { int i; for(i=0;i<2;i++) d[i]= a0[i]+a1[i]*t + a2[i]*t*t + a3[i]*t*t*t; } // 速度 d(x,y)/dt void B3::dpos(double t, double *d) { int i; for(i=0;i<2;i++) d[i]= a1[i] + 2*a2[i]*t + 3*a3[i]*t*t; }; // 加速度 void B3::ddpos(double t, double *d) { int i; for(i=0;i<2;i++) d[i]= 2*a2[i] + 6*a3[i]*t; };
Case1: 二つの曲線が平行な平面にある場合
それぞれのレールをなるべく均等に分割し点を置く。t 等分割でいいかも。曲線は３次ベジエを複数つないだものになっているので、各部分で tを等分割するが、各部分の長さを近似的に計算しておき、それにほぼ比例した数で分割する。分割した各点で、接線の向き (vx, vy)=(dx/dt, dy/dt) と曲がる方向を計算。接線の向きは atan2(vy, vx) で評価。曲がる向きは速度(vx,vy)と加速度 (ax,ay)=d(vx,vy)/dt の外積 vx*ay - vy*ax。これは符号だけ分かればいい。各点にセグメント番号を振る。始点で０でスタート。順に点を見ていって曲がる向きの符号が変わったら番号を＋１する。変化しなければ前と同じ番号。終点ではさらに＋１する。
レール１、２上の点v1, v2 を始点から順にたどっていく。
- 0: 初期値としてv1, v2両方共起点を選ぶ。v1, v2 を結ぶ線を引く。
- 1: ループ開始。v1,v2共にレールの終点に達していたら終了
- 2: v1, v2 のセグメント番号を比較
- セグメント番号が異なる場合はv1,v2でセグメント番号が小さいほうを一歩進める。そしてv1,v2の間に線を引いて１へ戻る
- 3: セグメント番号が同じ場合。v1, v2 上での接線の向き s1, s2 を計算。曲がる向きの符号が負ならそれぞれ-1倍する。
- s1＞s2 であれば、v2を一歩進める。そしてv1,v2の間に線を引いて１へ戻る
- そうでなければ、v1を一歩進める。そしてv1,v2の間に線を引いて１へ戻る
結んだ線は順番に記録しておき、前回の線と、今回一歩進んだ点からなる三角形を面として登録。
Case2: 二つの曲線が平行でない平面にある場合
二つの平面が交わる直線をy軸とする。二つのレール上の点を結べる条件は、各点の接線を伸ばしてy軸に交わった切片が双方で一致すること。y軸に垂直な、曲線1の横軸座標をx1とする。 x1=0での切片は B1= y1-x1*(vy1/vx1)と計算できる。2も同様。B1=B2 となる点同士を結べばいい。 dB/dt= (vy*ax - vx*ay)/(vx*vx+vy*vy) であるから、接線の角度 が単調増加ならB1 も単調増加。したがって平行な場合と同じセグメント番号を使う。
アルゴリズムは Case 1と同じで、ステップ３において傾きの代わりにB1,B2 を比較すればいい。ただし分母が０になる場合があるので、代わりに atan2(vx, vx*y-vy*x)で評価。
展開図作成
すでに平面での座標を求めた２つの頂点v1,v2 があり、三番目の頂点v3を追加して三角形v1,v2,v3を平面に書くことを考える。v1,v2の三次元での距離をL12などと書けば、平面でv1とv3 の距離がL13, v2とv3の距離がL23となる点をを探せばいい。余弦定理から角v2-v1-v3の角度をθとおくと cosθ=(L12*L12+L13*L13-L23*L23)/(2*L12*L13) となる。sin θは sqrt(1-cosθ^2) として、辺v1-v2 のどっちにv3を置くか決めておいて符号を決める。あとは座標変換でv3の位置は分かる。
これをすべての三角形について繰り返せば展開図ができる。
展開図出力
SVG フォーマットで二次元の点をつないだ線を出力。サンプル：
<?xml version="1.0" standalone="no"?>
<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 20001102//EN" "svg-20001102.dtd">
<svg width="744" height="1052" viewBox="0 0 744 1052">
<g style="fill:none; color:black; stroke:currentColor; stroke-width:1.00">
<path d='M100,100 L200,100 L200,200 L100,200 Z100,100 '></path>
</g>
</svg>
TODO
VTK ライブラリを使って３Dレンダリングする
