在重复的x个位置上用y点进行曲线拟合(Galaxy螺旋臂) [英] 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曲线作为函数来处理，而不是这种情况(对于相同的x，您会获得更多的y值)这就是为什么拟合在右侧失败的原因(当您击中非功能区域时).
Your fail is due to that you handle your 2D curve as function which is not the case (you got more y values for the same x) and that is why the fit 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 tighted 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 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,yssize of image in pixelsp[y][x].ddis pixel at (x,y) position as 32 bit integer typep[y][x].db[4]is pixel access by color bands (r,g,b,a)p.load(filename),p.save(filename)guess what ... loads/saves imagep.bmp->Canvasis GDI bitmap access so I can use GDI stuff too
拟合是由我的近似搜索类完成的:
The fitting is done by my approximation search class from:
因此，只需从此处复制class approx.
So just copy the class approx from there.
List<T>模板只是动态数组(列表):
The List<T> template is just dynamic array (list):
-
List<int> q;与int q[];相同
-
q.num保存内部的元素数量 -
q.add()向列表末尾添加新的空元素 -
q.add(10)将10作为新元素添加到列表末尾
List<int> q;is the same asint q[];q.numholds the number of elements insideq.add()adds new empty element to end of listq.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 ...
如果您需要 BEZIER 而不是插值多项式，则可以直接转换控制点，请参见:
If you need BEZIER instead 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:
-
对于
- 曲线始于
p1(t=0.0)，结束于p2(t=1.0).点p0,p3可通过(t=-1.0和t=2.0)到达，并确保面片之间的连续性条件.因此，对于所有相邻面片，在p0和p1处的推导都相同.就像将BEZIER补丁合并在一起一样.
4输入点p0,p1,p2,p3的- for
4input pointsp0,p1,p2,p3the curve starts atp1(t=0.0) and ends atp2(t=1.0). the pointsp0,p3are reachble through (t=-1.0andt=2.0) and ensure the continuity condition between patches. So the derivation atp0andp1are the same for all neighboring patches. It is the same as merging BEZIER patches together.
多项式拟合很容易:
-
我将
p1,p2设置为螺旋端点
I set the
p1,p2to the spiral endpoints
因此曲线的起点和终点应
so the curve starts and ends where it should
我在p1,p2附近搜索p0,p3达一定距离m
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螺旋臂)的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
