最近の季の節で風邪をひいてへたってますが、先日書いた「ベジェ曲線と点の距離」の続きを書きます。計算のやり方としては、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値への変換で使うつもりでつくった、というとこと。