NSBezierPathとNSPointの距離 その２（それなりにズルしてません）

まず

結果

で

動機
前回のを解析的(?)には難しそうなので数値計算で（？）やってみることを考えた．

方針
まず方針は前回の文中の追記に書いたものとする．

[x,y]=[X(t), Y(t)], D(t, px, py)=sqrt((X(t)-px)^2+(Y(t)-py)^2), d( D(t, px, py) ) / dt = 0 となる t

具体的には

１．ベジエ曲線の媒介変数表示(X(t), Y(t))と点(px, py)との距離の関数(D)を作成
２．これの微分したものが0となるtを求めることで距離が極値となる点を求める(tの解のなかで虚部が0となるもののみ．多くの場合５次だから一つはあるはず(?)）．
３．該当する点が複数あるケースを一応考慮すると距離が最小となるものを利用する．
４．虚部が0でもtが0より小さいや1より大きい場合はそれぞれt=0, t=1として処理して曲線端部を割り当てる．
５．ベジエ曲線の構成点を-px,-pyずらして式を簡略化しD(t)=sqrt((X(t))^2+(Y(t))^2)とする．（平行移動でありtを元の構成点での式に入れればも元のベジエ曲線での位置になる．）

求めたい式
手作業でやるのは苦労なことは目に見えているのでmaxima(portで入れた)に登場いただく．

f(t):= (t^3*cx4+3*t^2*(1-t)*cx3+3*t*(1-t)^2*cx2+(1-t)^3*cx1);

g(t):= (t^3*cy4+3*t^2*(1-t)*cy3+3*t*(1-t)^2*cy2+(1-t)^3*cy1);

v(t):=sqrt((f(t))^2+(g(t))^2);

rhs(solve(expand(diff(v(t),t)),t)[1]);

で出てきたのがこれ．


このままだとプログラムに使えないし，書き写すと間違うことには自信がある．ので，探したらありました(“付録”のところ)．そういうことをしてくれるlispファイルが．ありがとうございます．というかmaxima自体にあっても良いんじゃないかしらね．

ということで前述のlispファイルのあるフォルダでmaximaを起動して以下のようになります．

load("cform.lisp");

cform(rhs(solve(expand(diff(v(t),t)),t)[1]));

ーー出力ーー
(pow(cy4,2)+cy2*(6*cy4-18*cy3)-6*cy3*cy4+cy1*(-2*cy4+6*cy3-6*cy2)+9*
pow(cy3,2)+9*pow(cy2,2)+pow(cy1,2)+pow(cx4,2)+cx2*(6*cx4-18*cx3)-6*cx3*
cx4+cx1*(-2*cx4+6*cx3-6*cx2)+9*pow(cx3,2)+9*pow(cx2,2)+pow(cx1,2))*
pow(t,5)+(cy1*(5*cy4-20*cy3+25*cy2)+5*cy3*cy4+cy2*(45*cy3-10*cy4)-15*
pow(cy3,2)-30*pow(cy2,2)-5*pow(cy1,2)+cx1*(5*cx4-20*cx3+25*cx2)+5*cx3*
cx4+cx2*(45*cx3-10*cx4)-15*pow(cx3,2)-30*pow(cx2,2)-5*pow(cx1,2))*
pow(t,4)+(cy2*(4*cy4-36*cy3)+cy1*(-4*cy4+24*cy3-40*cy2)+6*pow(cy3,2)+36*
pow(cy2,2)+10*pow(cy1,2)+cx2*(4*cx4-36*cx3)+cx1*(-4*cx4+24*cx3-40*cx2)+6*
pow(cx3,2)+36*pow(cx2,2)+10*pow(cx1,2))*pow(t,3)+(cy1*(cy4-12*cy3+30*cy2)+9*
cy2*cy3-18*pow(cy2,2)-10*pow(cy1,2)+cx1*(cx4-12*cx3+30*cx2)+9*cx2*cx3-18*
pow(cx2,2)-10*pow(cx1,2))*pow(t,2)+(cy1*(2*cy3-10*cy2)+3*pow(cy2,2)+5*
pow(cy1,2)+cx1*(2*cx3-10*cx2)+3*pow(cx2,2)+5*pow(cx1,2))*t+cy1*
cy2-pow(cy1,2)+cx1*cx2-pow(cx1,2);

でこのtに関する方程式を解けば良いわけだ．

５次方程式/数値計算
で，５次方程式になっちゃったわけですが５次方程式は一般的な代数的解法は必ずしも存在しない/代数的な根の公式は存在しない そうです．特別な場合なら公式を得ることができたり，楕円関数を使えば…とかいろいろ書いてあったのですが，特殊かどうか判断できるほど頭よくないのと，楕円うんぬんは最終的に手順がみつけられず，さらには楕円関数のライブラリとかってどれなんだろう？とかいろいろ考えて結局，ライブラリを使うんだったら，数値計算しましょうということにしました．ライブラリは簡単に入れられそうだった(?)GSLを使うことにします(portで入れた)．
でGSLでの代数方程式のサンプルさがしてたまたま巡り会ったのがこちら．ということで準備が整いました．

コード
NSPoint(p0,c1,c2,p1)を取得するのは前回を参考にしてください．でNSCurveToBezierPathElementのケース部だけ書きます．

#include <gsl/gsl_complex.h>
#include <gsl/gsl_poly.h>

   //NSPoint p0= ...;
   //NSPoint c1= ...;
   //NSPoint d2= ...;
   //NSPoint p1= ...;
   
   double cx4= (double)p0.x-p.x;
   double cx3= (double)c1.x-p.x;
   double cx2= (double)c2.x-p.x;
   double cx1= (double)p1.x-p.x;
   double cy4= (double)p0.y-p.y;
   double cy3= (double)c1.y-p.y;
   double cy2= (double)c2.y-p.y;
   double cy1= (double)p1.y-p.y;
   
   double a[6]; 
   gsl_complex z[5];
   a[5]= (pow(cy4,2)+cy2*(6*cy4-18*cy3)-6*cy3*cy4+cy1*(-2*cy4+6*cy3-6*cy2)+9*pow(cy3,2)+9*pow(cy2,2)+pow(cy1,2)+pow(cx4,2)+cx2*(6*cx4-18*cx3)-6*cx3*cx4+cx1*(-2*cx4+6*cx3-6*cx2)+9*pow(cx3,2)+9*pow(cx2,2)+pow(cx1,2));
   a[4]= (cy1*(5*cy4-20*cy3+25*cy2)+5*cy3*cy4+cy2*(45*cy3-10*cy4)-15*pow(cy3,2)-30*pow(cy2,2)-5*pow(cy1,2)+cx1*(5*cx4-20*cx3+25*cx2)+5*cx3*cx4+cx2*(45*cx3-10*cx4)-15*pow(cx3,2)-30*pow(cx2,2)-5*pow(cx1,2));
   a[3]= (cy2*(4*cy4-36*cy3)+cy1*(-4*cy4+24*cy3-40*cy2)+6*pow(cy3,2)+36*pow(cy2,2)+10*pow(cy1,2)+cx2*(4*cx4-36*cx3)+cx1*(-4*cx4+24*cx3-40*cx2)+6*pow(cx3,2)+36*pow(cx2,2)+10*pow(cx1,2));
   a[2]= (cy1*(cy4-12*cy3+30*cy2)+9*cy2*cy3-18*pow(cy2,2)-10*pow(cy1,2)+cx1*(cx4-12*cx3+30*cx2)+9*cx2*cx3-18*pow(cx2,2)-10*pow(cx1,2));
   a[1]= (cy1*(2*cy3-10*cy2)+3*pow(cy2,2)+5*pow(cy1,2)+cx1*(2*cx3-10*cx2)+3*pow(cx2,2)+5*pow(cx1,2));
   a[0]= cy1*cy2-pow(cy1,2)+cx1*cx2-pow(cx1,2);  

   
   gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc(6);
   gsl_poly_complex_solve(a, 6, w, z);
   gsl_poly_complex_workspace_free(w);
   
   int k;
   float zd= FLT_MAX;
   NSPoint zp= NSZeroPoint;
   for(k= 0; k<5; k++) {
    if (GSL_IMAG(z[k])==0) {
     float t= (float)GSL_REAL(z[k]);
     t= (t<0.0)?0.0:t;
     t= (t>1.0)?1.0:t;
     float x= f(p0.x, c1.x, c2.x, p1.x, t);
     float y= f(p0.y, c1.y, c2.y, p1.y, t);
     NSPoint tmp_zp= NSMakePoint(x, y);
     float tmp_zd= length(tmp_zp, p);
     if (tmp_zd<zd) {
      zd= tmp_zd;
      zp= tmp_zp;
     }
    }
   }
   tmp_d= zd;
   tmp_ret= zp;

ということでtmp_d(zd)が最短距離でtmp_ret(zp)がベジエ曲線上の最近傍点(?)となります．

コード補足
gslのライブラリの使い方は基本的にzlibの使い方と同じです．/opt/local/lib/libgsl.dylibをプロジェクトに追加します．あとヘッダ検索パス/opt/local/includeを追加しました．

色々やったけど自信はないですw．