在重复的 x 位置上用 y 点进行曲线拟合（Galaxy Spiral arms）
Curve fitting with y points on repeated x positions (Galaxy Spiral arms)
我目前有一个 MATLAB 程序，它从星系中获取追踪旋臂的 RGB 图像并选择最大的旋臂组件并仅绘制它。
我尝试使用 matlab 的内置曲线拟合工具和平滑样条来拟合它，我得到以下结果：
我试过将 interp1 与参数拟合一起使用，结果却很糟糕。
有没有办法完全适应这种类型的曲线？
您的错误结果是由于您将 2D 曲线作为函数处理，但事实并非如此（对于相同的 x，您有多个 y 值），这就是为什么右侧拟合失败（当您点击非功能区域时）。
要解决这个问题，您需要将曲线拟合与每个维度分开。因此，您可以将每个轴作为单独的函数进行拟合。为此，您需要使用不同的函数参数（不是 x）。如果您以某种方式对您的点进行排序（例如，通过与起点的曲线距离，或通过极角或任何其他方式），那么您可以使用点索引作为此类函数参数。
所以你做了这样的事情：
y(x) = fit((x0,y0),(x1,y1),(x2,y2)...)
其中 returns 是 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 中为您的任务制作了一个简单的三次插值的简单示例，以便您更好地理解我的意思：
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 中提供的输入图像。这是阶段输出
螺旋点选择：
点按极角排序：
最终拟合结果：
如您所见，拟合度不是很好，因为：
- 我对整个手臂使用单个立方体（它可能需要更大次数的多项式或细分到补丁）
- 我使用从图像中提取点并且曲线很粗所以我每个极角有多个点（你有原始点所以这应该不是问题）我应该使用细化算法但懒得添加它...
在 C++ 示例中，我使用了自己的图像 class 所以这里有一些成员：
xs,ys 图像大小（以像素为单位）p[y][x].dd 是 (x,y) 位置的像素为 32 位整数类型p[y][x].db[4] 是色带 (r,g,b,a) 的像素访问p.load(filename),p.save(filename)猜猜看...loads/saves图片p.bmp->Canvas 是 GDI 位图访问所以我也可以使用 GDI 东西
拟合是通过我的近似搜索 class 完成的：
所以只需从那里复制 approx class。
List<T>模板只是一个动态数组（列表）类型：
List<int> q; 等同于 int q[];q.num保存里面的元素个数q.add() 添加新的空元素到列表末尾q.add(10) 添加 10 作为新元素到列表末尾
[注释]
因为你已经有了点列表，所以你不需要扫描输入图像来寻找点...所以你可以忽略那部分代码...
如果你需要Bézier而不是插值多项式那么你可以直接转换控制点，见：
- Interpolation cubic vs. Bézier cubic
如果目标曲线形式不固定，那么你也可以尝试用一些参数化的圆直接拟合螺旋方程，比如带有移动中心和变半径的方程。这应该更精确，并且大多数参数无需拟合即可计算。
[Edit1] 更好地描述了我的多项式拟合
我正在使用上面 link 的插值立方体，这些属性：
- 对于
4 输入点 p0,p1,p2,p3 曲线开始于 p1 (t=0.0) 并结束于 p2 (t=1.0)。点 p0,p3 可通过（t=-1.0 和 t=2.0）到达，并确保补丁之间的连续性条件。所以 p0 和 p1 的推导对于所有相邻的补丁都是相同的。这与将贝塞尔补丁合并在一起是一样的。
多项式拟合很简单：
- 我将
p1,p2 设置为螺旋端点
所以曲线在它应该开始和结束的地方
- 我在
p1,p2 附近搜索 p0,p3 直到一段距离 m
同时记住多项式曲线与原始点的最接近匹配。您可以为此使用平均或最大距离。 approx class 完成您需要的所有工作，只需在每次迭代中计算距离 ee。
对于 m 我使用图像大小的倍数。如果太大你会失去精度（或者需要更多的递归和减慢速度），如果太小你可以限制控制点应该在的区域并且拟合会变形。
迭代起始步 dm 是 m 的一部分，如果太小计算会很慢。如果太高，您可能会错过局部 min/max 解决方案导致不合适的地方。
为了加快计算速度，我只使用了从点中均匀选择的 25 个点（不需要全部使用）步骤在 di
维数分离x,y是一样的，你只要把x全部换成y，否则代码是一样的
我目前有一个 MATLAB 程序，它从星系中获取追踪旋臂的 RGB 图像并选择最大的旋臂组件并仅绘制它。
我尝试使用 matlab 的内置曲线拟合工具和平滑样条来拟合它，我得到以下结果：
我试过将 interp1 与参数拟合一起使用，结果却很糟糕。
有没有办法完全适应这种类型的曲线？
您的错误结果是由于您将 2D 曲线作为函数处理，但事实并非如此（对于相同的 x，您有多个 y 值），这就是为什么右侧拟合失败（当您点击非功能区域时）。
要解决这个问题，您需要将曲线拟合与每个维度分开。因此，您可以将每个轴作为单独的函数进行拟合。为此，您需要使用不同的函数参数（不是 x）。如果您以某种方式对您的点进行排序（例如，通过与起点的曲线距离，或通过极角或任何其他方式），那么您可以使用点索引作为此类函数参数。
所以你做了这样的事情：
y(x) = fit((x0,y0),(x1,y1),(x2,y2)...)
其中 returns 是 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 中为您的任务制作了一个简单的三次插值的简单示例，以便您更好地理解我的意思：
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 中提供的输入图像。这是阶段输出
螺旋点选择：
点按极角排序：
最终拟合结果：
如您所见，拟合度不是很好，因为：
- 我对整个手臂使用单个立方体（它可能需要更大次数的多项式或细分到补丁）
- 我使用从图像中提取点并且曲线很粗所以我每个极角有多个点（你有原始点所以这应该不是问题）我应该使用细化算法但懒得添加它...
在 C++ 示例中，我使用了自己的图像 class 所以这里有一些成员：
xs,ys图像大小（以像素为单位）p[y][x].dd是 (x,y) 位置的像素为 32 位整数类型p[y][x].db[4]是色带 (r,g,b,a) 的像素访问p.load(filename),p.save(filename)猜猜看...loads/saves图片p.bmp->Canvas是 GDI 位图访问所以我也可以使用 GDI 东西
拟合是通过我的近似搜索 class 完成的：
所以只需从那里复制 approx class。
List<T>模板只是一个动态数组（列表）类型：
List<int> q;等同于int q[];q.num保存里面的元素个数q.add()添加新的空元素到列表末尾q.add(10)添加 10 作为新元素到列表末尾
[注释]
因为你已经有了点列表，所以你不需要扫描输入图像来寻找点...所以你可以忽略那部分代码...
如果你需要Bézier而不是插值多项式那么你可以直接转换控制点，见：
- Interpolation cubic vs. Bézier cubic
如果目标曲线形式不固定，那么你也可以尝试用一些参数化的圆直接拟合螺旋方程，比如带有移动中心和变半径的方程。这应该更精确，并且大多数参数无需拟合即可计算。
[Edit1] 更好地描述了我的多项式拟合
我正在使用上面 link 的插值立方体，这些属性：
- 对于
4输入点p0,p1,p2,p3曲线开始于p1(t=0.0) 并结束于p2(t=1.0)。点p0,p3可通过（t=-1.0和t=2.0）到达，并确保补丁之间的连续性条件。所以p0和p1的推导对于所有相邻的补丁都是相同的。这与将贝塞尔补丁合并在一起是一样的。
多项式拟合很简单：
- 我将
p1,p2设置为螺旋端点
所以曲线在它应该开始和结束的地方
- 我在
p1,p2附近搜索p0,p3直到一段距离m
同时记住多项式曲线与原始点的最接近匹配。您可以为此使用平均或最大距离。 approx class 完成您需要的所有工作，只需在每次迭代中计算距离 ee。
对于 m 我使用图像大小的倍数。如果太大你会失去精度（或者需要更多的递归和减慢速度），如果太小你可以限制控制点应该在的区域并且拟合会变形。
迭代起始步 dm 是 m 的一部分，如果太小计算会很慢。如果太高，您可能会错过局部 min/max 解决方案导致不合适的地方。
为了加快计算速度，我只使用了从点中均匀选择的 25 个点（不需要全部使用）步骤在 di
维数分离x,y是一样的，你只要把x全部换成y，否则代码是一样的