可展面作成アルゴリズム Ver 0.2

先日書いたアルゴリズムからいろいろ変更。まず、一般的なベクター描画ソフト、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レンダリングする