在重复的 x 位置上使用 y 点进行曲线拟合(Galaxy Spiral 臂) [英] Curve fitting with y points on repeated x positions (Galaxy Spiral arms)

问题描述

我目前有一个 MATLAB 程序，它从星系中获取跟踪的旋臂的 RGB 图像，并选择最大的旋臂组件并仅绘制它.

I currently have a MATLAB program which takes RGB images of traced spiral arms from galaxies and selects the biggest arm component and plots only that.

我尝试使用 matlab 的内置曲线拟合工具和平滑样条来拟合它，我得到以下结果:

I have tried using matlab's built in curve fitting tool with smoothing spline to fit it and I get the following result:

我曾尝试将 interp1 与参数拟合一起使用，但结果却很糟糕.

I have tried using interp1 with parametric fitting to only get bad results.

有没有办法完全拟合这种类型的曲线?

Is there a way to fit this type of curve at all?

推荐答案

您的错误结果是由于您将 2D 曲线作为函数处理，而事实并非如此(您有多个 y相同 x 的值)，这就是为什么拟合在右侧失败(当您击中非功能区域时).

Your erroneous result is due to the fact that you handle your 2D curve as function which is not the case (you have multiple y values for the same x) and that is why the fitting fails on the right side (when you hit the non-function area).

为了解决您需要将曲线拟合到每个维度的问题.因此，您可以将每个轴作为单独的函数进行拟合.为此，您需要使用不同的函数参数(不是 x).如果您以某种方式对点进行排序(例如，通过距起点的曲线距离、极角或其他任何方式)，那么您可以使用点索引作为此类函数参数.

To remedy that you need to separate the curve fit to each dimension. So you can fit each axis as separate function. For that you need to use different function parameter (not x). If you order your points somehow (for example by curve distance from start point, or by polar angle or by what ever) then you can use point index as such function parameter.

所以你做了这样的事情:

So you have done something like this:

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

返回 y(x) 的多项式.相反，您应该这样做:

which returns a polynomial for y(x). Instead you should do something like this:

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 ，如下所示:

where t is your new parameter tightened to the order of point in the ordered list. Most curves use parameter in range t=<0.0,1.0> to ease up the computation and usage. So if you got N points then you can convert point index i=<0,N-1> to curve parameter t like this:

t=i/(N-1);

绘图时您需要更改您的

plot(x,y(x))

到

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

我在 C++/VCL 中为您的任务制作了一个简单的三次方插值简单示例，以便您更好地理解我的意思:

I made a simple example of simple single interpolation of a cubic in C++/VCL for your task so you better see what I mean:

    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 中提供的输入图像.这是阶段输出

I used your input image you provided in OP. Here are the stages outputs

螺旋点选择:

按极角排序的点:

最终拟合结果:

如您所见，贴合度不是很好，因为:

• 我对整个手臂使用单三次(它可能需要更大的多项式或细分到补丁)
• 我使用图像中的点提取，曲线很粗，所以每个极角都有多个点(你得到了原始点，所以这应该不是问题)我应该使用细化算法但懒得添加它...

在 C++ 示例中，我使用自己的图像类，因此这里有一些成员:

In the C++ example I use my own image class so here some members:

• xs,ys 图像的像素大小
• p[y][x].dd 是 (x,y) 位置的像素，为 32 位整数类型
• p[y][x].db[4] 是按色带 (r,g,b,a) 的像素访问
• p.load(filename),p.save(filename) 猜猜...加载/保存图像
• p.bmp->Canvas 是 GDI 位图访问，所以我也可以使用 GDI 的东西

• xs,ys size of image in pixels
• p[y][x].dd is pixel at (x,y) position as 32 bit integer type
• p[y][x].db[4] is pixel access by color bands (r,g,b,a)
• p.load(filename),p.save(filename) guess what ... loads/saves image
• p.bmp->Canvas is GDI bitmap access so I can use GDI stuff too

拟合是由我的近似搜索类完成的:

The fitting is done by my approximation search class from:

• 近似搜索的工作原理

所以只需从那里复制 approx 类.

So just copy the approx class from there.

List 模板只是一个动态数组(列表)类型:

The List<T> template is just a dynamic array (list) type:

• Listq; 与 int q[];
相同
• q.num 保存里面的元素个数
• q.add() 向列表末尾添加新的空元素
• q.add(10) 将 10 作为新元素添加到列表末尾

• List<int> q; is the same as int q[];
• q.num holds the number of elements inside
• q.add() adds new empty element to end of list
• q.add(10) adds 10 as new element to end of list

[注释]

因为你已经有了点列表，所以你不需要扫描输入图像的点......所以你可以忽略那部分代码......

As you already have the point list then you do not need to scan input image for points... so you can ignore that part of code ...

如果您需要 Bézier 而不是插值多项式，那么您可以直接转换控制点，请参阅:

If you need Bézier instead of an interpolation polynomial then you can convert the control points directly, see:

• 插值三次与贝塞尔三次方

如果目标曲线形式不是固定的，那么你也可以尝试通过一些参数圆，如移心变半径方程，直接拟合螺旋方程.这应该更精确，并且大多数参数无需拟合即可计算.

If the target curve form is not fixed then you can also try to fit directly the spiral equation by some parametric circle like equation with shifting center and variable radius. That should be far more precise and most of the parameters can be computed without fitting.

[Edit1] 对矿井多项式拟合的更好描述

我正在使用上面链接中的插值三次，并具有这些属性:

I am using interpolation cubic from the above link whit these properties:

• 对于 4 个输入点 p0,p1,p2,p3 曲线从 p1 (t=0.0>) 并在 p2 (t=1.0) 处结束.点p0,p3是通过(t=-1.0 and t=2.0)可达的，保证了patch之间的连续性条件.所以 p0 和 p1 的推导对于所有相邻的补丁都是相同的.这与合并 Bézier 补丁相同.

• for 4 input points p0,p1,p2,p3 the curve starts at p1 (t=0.0) and ends at p2 (t=1.0). the points p0,p3 are reachable through (t=-1.0 and t=2.0) and ensure the continuity condition between patches. So the derivation at p0 and p1 are the same for all neighboring patches. It is the same as merging Bézier patches together.

多项式拟合很简单:

1. 我将 p1,p2 设置为螺旋端点

所以曲线在它应该的位置开始和结束

so the curve starts and ends where it should

1. 我在 p1,p2 附近搜索 p0,p3 到一定距离 m

1. I search p0,p3 near p1,p2 up to some distance m

同时记住多项式曲线与原始点的最接近匹配.您可以为此使用平均或最大距离.approx 类完成您需要的所有工作，只需计算每次迭代中的距离 ee.

while remembering the closest match of polynomial curve to original points. You can use average or max distance for this. The approx class do all the work you need just to compute the distance ee in each iteration.

对于 m 我使用图像大小的倍数.如果太大你会失去精度(或者需要更多的递归和减慢速度)，如果太低你可以限制控制点应该在的区域并且拟合会变形.

for m I use multiple of image size. If too big you will lose precision (or need more recursions and slow things down), if too low you can bound out the area where the control points should be and the fit will be deformed.

迭代开始步骤dm 是m 的一部分，如果太小计算会很慢.如果太高，您可能会错过解决方案导致错误拟合的局部最小值/最大值.

Iteration starting step dm is part of m and if too small computation will be slow. If to high you can miss local min/max where the solution is resulting in wrong fit.

为了加速计算，我只使用了从点中均匀选择的 25 个点(不需要使用所有点)，步骤在 di

To speed up the computation I am using only 25 points evenly selected from the points (no need to use all of them) the step is in di

维度分离x,y是一样的，你只要把x全部改成y，否则代码是一样的.

The dimension separation x,y is the same, you just change all the x for y, otherwise the code is the same.

这篇关于在重复的 x 位置上使用 y 点进行曲线拟合(Galaxy Spiral 臂)的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！