最近の季の節で風邪をひいてへたってますが、先日書いた「ベジェ曲線と点の距離」の続きを書きます。計算のやり方としては、2次ベジェ関数から導いた距離の数式を微分して解を求める、という感じ。
さて、2次ベジェのx,y座標を、tの関数とすると、それぞれ
// 始点 point0 // 終点 point1 // コントロール点 control //x座標の関数 fx(t) = point0.x*tp*tp + 2*control.x*t*tp + point1.x*t*t; //y座標の関数 fy(t) = point0.y*tp*tp + 2*control.y*t*tp + point1.y*t*t;
というような2次式になります。ここで、任意の点(x,y)からの距離を考えると、距離はtの関数で
//点(x,y)と曲線の距離平方の関数 fd2(t) = (fx(t)-x)*(fx(t)-x) + (fy(t)-y)*(fy(t)-y); //点(x,y)と曲線の距離の関数 fd(t) = Math.sqrt( fd2(t) );
になります。
ただ、この式だと展開した後面倒なので、point0、point1、controlをそれぞれ(-x,-y)だけoffsetして、単純に原点と曲線の距離の問題になるように考えます。
で、任意点からの近傍点を求めるわけですが、これはfd2(t)が最小になるtを求めることになります。距離は近づいたり離れたり、連続的な変化をしてるので、fd2(t)を微分した値の傾きが0の時が、距離の山か谷、ということになります。
fd2(t)は4次式なので、微分した関数fd2(t)’は3次式です。これの解を求めます。
ただ、3次式の解には虚数が含まれるので計算はやっかいで、この問題に限っては次の条件で限定的な計算をします。条件というのは、
・tが0から1の値
・fd2(t)’の3次式は、t^3の係数が正なので、山・谷の順の関数になる。
※fd2(t)が2次式の二乗を足しているからです。
・実数解が得られる場合は、それが最近傍
・実数解が得られない場合は、数値計算で求める
・解が0から1でない場合は、端点の近いほうを解とする
なにやら面倒ですね…
実数解が得られない場合の数値計算をする場合、fd2(t)’をもう一度微分して、その山と谷にあたるt値を求めて、
・t=0~山の値
・t=谷~1
の間でfd2(t)’が0になるt値を再帰的に求めます。それぞれの関数はこんな曲線になります。
下部の青地のところが関数のグラフです。青が距離の関数fd(t)、黄が微分した関数fd2(t)’、赤がさらに微分したもの。なんとなくイメージつきますかね…。スクリプトは次のような感じ。ちょっとぬるいので、どこかで発散しそうな気もしますが…。
/** * 2次ベジェ曲線と点との最短距離 t値 * @param ベジェ始点 * @param ベジェ終点 * @param ベジェコントロール点 * @param Point * @return t値 */ function nearestFromPoint( point0:Point, point1:Point, control:Point, p:Point ):Number{ //ベジェ曲線を(-p)でオフセットした座標 var x0:Number = point0.x - p.x; var y0:Number = point0.y - p.y; var x1:Number = point1.x - p.x; var y1:Number = point1.y - p.y; var cx:Number = control.x - p.x; var cy:Number = control.y - p.y; //計算パラメータ var xx:Number = x0+x1; var yy:Number = y0+y1; var cc:Number = cx*cx + cy*cy; //曲線と点の距離平方の微分計算定数 // fd(t) = A*t^3 + B*t^2 + C*t + D var A:Number = 4*( xx*xx + yy*yy + 4*cc - 4*cx*xx - 4*cy*yy ); var B:Number = -12*( x0*xx + y0*yy + 2*cc + (-3*x0-x1)*cx + (-3*y0-y1)*cy ); var C:Number = 4*( (3*x0+x1-6*cx)*x0 + (3*y0+y1-6*cy)*y0 + 2*cc ); var D:Number = -4*(y0*y0-cy*y0+x0*x0-cx*x0); var TA:Number = 27*A*A*D*D + (4*B*B*B-18*A*B*C)*D + 4*A*C*C*C - B*B*C*C; if( TA>0 ){ //実数解 TA = Math.sqrt(TA)/( 6*Math.sqrt(3)*A*A ) - (27*A*A*D - 9*A*B*C + 2*B*B*B )/(54*A*A*A); if( TA<0 ){ TA = -Math.pow(-TA,1/3); }else{ TA = Math.pow(TA,1/3); } return Math.max( 0.0, Math.min( 1.0, (TA - (3*A*C-B*B)/(9*A*A*TA) - B/(3*A)) ) ); }else if( TA<0 ){ //虚数解あり、数値計算で //距離平方の微分 ※2次関数の平方の和なので、Aは正→山・谷の3次式 var diff:Function = function(t:Number):Number{ return (A*t*t*t + B*t*t + C*t + D); }; //2回微分「3*A*t*t + 2*B*t + C」の解 var tt0:Number = (-B - Math.sqrt(B*B-3*A*C))/(3*A); var tt1:Number = (-B + Math.sqrt(B*B-3*A*C))/(3*A) //最適解 var t0:Number = 0; var t1:Number = 1.0; var ot0:Number; var ot1:Number; var dt:Number; if( tt0>0 && diff(0)<0 ){ //数値計算 ot0 = 0.0; ot1 = tt0; t0 = (ot0 + ot1)/2; while( Math.abs( dt=diff(t0) )>ERR_T ){ if( dt>0 ){ ot1 = t0; t0 = (ot0 + ot1)/2; }else{ ot0 = t0; t0 = (ot0 + ot1)/2; } } } if( tt1<1.0 && diff(1.0)>0 ){ //数値計算 ot0 = tt1; ot1 = 1.0; t1 = (ot0 + ot1)/2; while( Math.abs( dt=diff(t1) )>ERR_T ){ if( dt<0 ){ ot0 = t1; t1 = (ot0 + ot1)/2; }else{ ot1 = t1; t1 = (ot0 + ot1)/2; } } } var dv0:Point = new Point( t0*t0*(x1+x0-2*cx)+t0*(2*cx-2*x0)+x0, t0*t0*(y1+y0-2*cy)+t0*(2*cy-2*y0)+y0 ); var dv1:Point = new Point( t1*t1*(x1+x0-2*cx)+t1*(2*cx-2*x0)+x0, t1*t1*(y1+y0-2*cy)+t1*(2*cy-2*y0)+y0 ); //とりあえず近いほう。同じ距離は無視しとく return ( dv0.length<dv1.length ) ? t0 : t1; }else{ //コントロール点が、始点と終点の中点 return -D/C; } }
実際は、解が2つある場合、それぞれ同じ距離って場合もあるのですが、そこはまぁええか、という感じです。
距離を使って衝突判定ライクなことをすると、こんな感じです。あんまりちゃんとしてないので、抜けたりしますけど、そこはご愛嬌です。あるいは距離を測らなくても、軌道直線と曲線の交点を求めればいいような気もしますが…
やはりというか、使う場面は限られているようにも思えますが…。僕の場合は、曲線と曲線の交点を出す場合に、t値ではなく座標がでてくる感じなので、座標からt値への変換で使うつもりでつくった、というとこと。