| INDEX |

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

コメントは受け付けていません。

| INDEX |