matlab - 繰り返される x 位置の y 点によるカーブフィッティング (Galaxy Spiral arm)

Question

私は現在、銀河から追跡された渦巻き腕のRGB画像を取得し、最大の腕成分を選択してそれのみをプロットするMATLABプログラムを持っています。

曲線フィッティングの準備が整ったスパイラルアーム

matlab のビルトインカーブフィッティングツールとスムージングスプラインを使用してフィットさせてみましたが、次の結果が得られました。

MATLAB を使用したスパイラルアーム

パラメトリックフィッティングで interp1 を使用して、悪い結果しか得られないようにしました。

このタイプの曲線に合わせる方法はありますか?

score 3 · Accepted Answer

あなたの誤った結果は、2D 曲線を関数として処理しないという事実によるものであり (y同じに対して複数の値がありxます)、右側でフィッティングが失敗するのはそのためです (非関数領域にヒットしたとき)。）。

これを改善するには、カーブフィットを各次元に分離する必要があります。したがって、各軸を個別の関数として適合させることができます。そのためには、別の関数パラメーター (x ではない) を使用する必要があります。何らかの方法でポイントを並べ替える場合 (たとえば、始点からの曲線の距離、極角など) は、そのような関数パラメーターとしてポイントインデックスを使用できます。

だからあなたはこのようなことをしました：

y(x) = fit((x0,y0),(x1,y1),(x2,y2)...)

の多項式を返しますy(x)。代わりに、次のようにする必要があります。

x(t) = fit(( 0,x0),( 1,x1),( 2,x2)...)
y(t) = fit(( 0,y0),( 1,y1),( 2,y2)...)

t順序付けられたリストのポイントの順序に締められた新しいパラメーターはどこにありますか。ほとんどの曲線は、範囲内のパラメータを使用t=<0.0,1.0>して、計算と使用を容易にします。したがって、Nポイントを取得した場合は、次のようにポイントインデックスi=<0,N-1>をカーブパラメータに変換できますt。

t=i/(N-1);

プロットするときは、変更する必要があります

plot(x,y(x))

に

plot(x(t),y(t))

あなたのタスクのためにC++/VCLで3 次の単純な単一補間の簡単な例を作成したので、私の意味がよくわかります。

    picture pic0,pic1;
        // pic0 - source
        // pic1 - output
    int x,y,i,j,e,n;
    double x0,x1,x2,x3,xx;
    double y0,y1,y2,y3,yy;
    double d1,d2,t,tt,ttt;
    double ax[4],ay[4];
    approx a0,a3; double ee,m,dm; int di;
    List<_point> pnt;
    _point p;

    // [extract points from image]
    pic0.load("spiral_in.png");
    pic1=pic0;
    // scan image
    xx=0.0; x0=pic1.xs;
    yy=0.0; y0=pic1.ys;
    for (y=0;y<pic1.ys;y++)
     for (x=0;x<pic1.xs;x++)
      // select red pixels
      if (DWORD(pic1.p[y][x].dd&0x00008080)==0)     // low blue,green
       if (DWORD(pic1.p[y][x].dd&0x00800000)!=0)    // high red
        {
        // recolor to green (just for visual check)
        pic1.p[y][x].dd=0x0000FF00;
        // add found point to a list
        p.x=x;
        p.y=y;
        p.a=0.0;
        pnt.add(p);
        // update bounding box
        if (x0>p.x) x0=p.x;
        if (xx<p.x) xx=p.x;
        if (y0>p.y) y0=p.y;
        if (yy<p.y) yy=p.y;
        }
    // center of bounding box for polar sort origin
    x0=0.5*(x0+xx);
    y0=0.5*(y0+yy);
    // draw cross (for visual check)
    x=x0; y=y0; i=16;
    pic1.bmp->Canvas->Pen->Color=clBlack;
    pic1.bmp->Canvas->MoveTo(x-i,y);
    pic1.bmp->Canvas->LineTo(x+i,y);
    pic1.bmp->Canvas->MoveTo(x,y-i);
    pic1.bmp->Canvas->LineTo(x,y+i);
    pic1.save("spiral_fit_0.png");
    // cpmpute polar angle for sorting
    for (i=0;i<pnt.num;i++)
        {
        xx=atan2(pnt[i].y-y0,pnt[i].x-x0);
        if (xx>0.75*M_PI) xx-=2.0*M_PI; // start is > -90 deg
        pnt[i].a=xx;
        }
    // bubble sort by angle (so index in point list can be used as curve parameter)
    for (e=1;e;)
     for (e=0,i=1;i<pnt.num;i++)
      if (pnt[i].a>pnt[i-1].a)
        {
        p=pnt[i];
        pnt[i]=pnt[i-1];
        pnt[i-1]=p;
        e=1;
        }
    // recolor to grayscale gradient (for visual check)
    for (i=0;i<pnt.num;i++)
        {
        x=pnt[i].x;
        y=pnt[i].y;
        pic1.p[y][x].dd=0x00010101*((250*i)/pnt.num);
        }
    pic1.save("spiral_fit_1.png");

    // [fit spiral points with cubic polynomials]
    n =6;                               // recursions for accuracy boost
    m =fabs(pic1.xs+pic1.ys)*1000.0;    // radius for control points fiting
    dm=m/50.0;                          // starting step for approx search
    di=pnt.num/25; if (di<1) di=1;      // skip most points for speed up
    // fit x axis polynomial
    x1=pnt[0          ].x;  // start point of curve
    x2=pnt[  pnt.num-1].x;  // endpoint of curve
    for (a0.init(x1-m,x1+m,dm,n,&ee);!a0.done;a0.step())
    for (a3.init(x2-m,x2+m,dm,n,&ee);!a3.done;a3.step())
        {
        // compute actual polynomial
        x0=a0.a;
        x3=a3.a;
        d1=0.5*(x2-x0);
        d2=0.5*(x3-x1);
        ax[0]=x1;
        ax[1]=d1;
        ax[2]=(3.0*(x2-x1))-(2.0*d1)-d2;
        ax[3]=d1+d2+(2.0*(-x2+x1));
        // compute its distance to points as the fit error e
        for (ee=0.0,i=0;i<pnt.num;i+=di)
            {
            t=double(i)/double(pnt.num-1);
            tt=t*t;
            ttt=tt*t;
            x=ax[0]+(ax[1]*t)+(ax[2]*tt)+(ax[3]*ttt);
            ee+=fabs(pnt[i].x-x);                   // avg error
//          x=fabs(pnt[i].x-x); if (ee<x) ee=x;     // max error
            }
        }
    // compute final x axis polynomial
    x0=a0.aa;
    x3=a3.aa;
    d1=0.5*(x2-x0);
    d2=0.5*(x3-x1);
    ax[0]=x1;
    ax[1]=d1;
    ax[2]=(3.0*(x2-x1))-(2.0*d1)-d2;
    ax[3]=d1+d2+(2.0*(-x2+x1));
    // fit y axis polynomial
    y1=pnt[0          ].y;  // start point of curve
    y2=pnt[  pnt.num-1].y;  // endpoint of curve
    m =fabs(y2-y1)*1000.0;
    di=pnt.num/50; if (di<1) di=1;
    for (a0.init(y1-m,y1+m,dm,n,&ee);!a0.done;a0.step())
    for (a3.init(y2-m,y2+m,dm,n,&ee);!a3.done;a3.step())
        {
        // compute actual polynomial
        y0=a0.a;
        y3=a3.a;
        d1=0.5*(y2-y0);
        d2=0.5*(y3-y1);
        ay[0]=y1;
        ay[1]=d1;
        ay[2]=(3.0*(y2-y1))-(2.0*d1)-d2;
        ay[3]=d1+d2+(2.0*(-y2+y1));
        // compute its distance to points as the fit error e
        for (ee=0.0,i=0;i<pnt.num;i+=di)
            {
            t=double(i)/double(pnt.num-1);
            tt=t*t;
            ttt=tt*t;
            y=ay[0]+(ay[1]*t)+(ay[2]*tt)+(ay[3]*ttt);
            ee+=fabs(pnt[i].y-y);                   // avg error
//          y=fabs(pnt[i].y-y); if (ee<y) ee=y;     // max error
            }
        }
    // compute final y axis polynomial
    y0=a0.aa;
    y3=a3.aa;
    d1=0.5*(y2-y0);
    d2=0.5*(y3-y1);
    ay[0]=y1;
    ay[1]=d1;
    ay[2]=(3.0*(y2-y1))-(2.0*d1)-d2;
    ay[3]=d1+d2+(2.0*(-y2+y1));
    // draw fited curve in Red
    pic1.bmp->Canvas->Pen->Color=clRed;
    pic1.bmp->Canvas->MoveTo(ax[0],ay[0]);
    for (t=0.0;t<=1.0;t+=0.01)
        {
        tt=t*t;
        ttt=tt*t;
        x=ax[0]+(ax[1]*t)+(ax[2]*tt)+(ax[3]*ttt);
        y=ay[0]+(ay[1]*t)+(ay[2]*tt)+(ay[3]*ttt);
        pic1.bmp->Canvas->LineTo(x,y);
        }
    pic1.save("spiral_fit_2.png");

OPで提供した入力画像を使用しました。ここにステージ出力があります

らせん点の選択:

極角順のポイント:

最終適合結果:

ご覧のとおり、適合はあまり良くありません。理由は次のとおりです。

腕全体に 1 つの 3 次を使用します (多項式の次数が大きくなるか、パッチに分割する必要がある場合があります)。
画像からのポイント抽出を使用し、曲線が太いため、単一の極角ごとに複数のポイントがあります（元のポイントを取得したので問題ありません）間引きアルゴリズムを使用する必要がありますが、追加するのが面倒です...

C++の例では、独自の画像クラスを使用しているため、ここにいくつかのメンバーがあります。

xs,ysピクセル単位の画像サイズ
p[y][x].dd32 ビット整数型としての (x,y) 位置のピクセル
p[y][x].db[4]カラーバンド (r、g、b、a) によるピクセルアクセスです。
p.load(filename),p.save(filename)何を推測する...画像の読み込み/保存
p.bmp->CanvasGDIビットマップアクセスなので、GDIのものも使用できます

フィッティングは、次の近似検索クラスによって行われます。

近似検索の仕組み

approxそこからクラスをコピーするだけです。

List<T>テンプレートは単なる動的配列 (リスト) 型です。

List<int> q; と同じですint q[];
q.num内部の要素数を保持します
q.add()リストの最後に新しい空の要素を追加します
q.add(10)リストの最後に新しい要素として 10 を追加します

[ノート]

すでにポイントリストがあるので、ポイントの入力画像をスキャンする必要はありません...そのため、コードのその部分を無視できます...

補間多項式の代わりにベジエが必要な場合は、制御点を直接変換できます。次を参照してください。

補間 3 次対ベジエ 3 次

ターゲットの曲線の形が固定されていない場合は、中心が移動し、半径が可変の方程式のようなパラメトリックな円によって、らせん方程式を直接当てはめることもできます。これははるかに正確なはずであり、ほとんどのパラメーターはフィッティングなしで計算できます。

[Edit1]鉱山多項式フィッティングのより良い説明

上記のリンクから補間キュービックを使用しています。これらのプロパティは次のとおりです。

4入力ポイントの場合p0,p1,p2,p3、曲線は ( ) で始まり( p1)t=0.0で終わります。ポイントは (および) を介して到達可能であり、パッチ間の連続性条件を保証します。したがって、およびでの導出は、隣接するすべてのパッチで同じです。ベジエパッチをマージするのと同じです。p2t=1.0p0,p3t=-1.0t=2.0p0p1

多項式の当てはめは簡単です：

p1,p2をらせんの端点に設定します

そのため、曲線は本来あるべき場所で開始および終了します

p0,p3近くp1,p2から遠くまで探すm

多項式曲線の元の点への最も近い一致を記憶しながら。これには平均距離または最大距離を使用できます。クラスは、各反復approxで距離を計算するために必要なすべての作業を行います。ee

m画像サイズの倍数を使用するためです。大きすぎると、精度が失われます (または、より多くの再帰が必要になり、速度が低下します)。小さすぎると、コントロールポイントが必要な領域が制限され、フィットが変形します。

反復開始ステップdmは一部でmあり、小さすぎると計算が遅くなります。高すぎると、ソリューションが不適切な適合をもたらすローカルの最小値/最大値を見逃す可能性があります。

計算を高速化するために、ポイントから均等に選択された25ポイントのみを使用しています（すべてを使用する必要はありません）ステップはdi

次元の分離x,yは同じです。すべてのxfor を変更するだけyです。それ以外のコードは同じです。

matlab - 繰り返される x 位置の y 点によるカーブ フィッティング (Galaxy Spiral arm)

1 に答える 1

Related

Reference

matlab - 繰り返される x 位置の y 点によるカーブフィッティング (Galaxy Spiral arm)