在matlab中计算闭合曲线（或多边形）的曲率

Question

考虑以下一组要点

x = [1.34, 0.92, 0.68, 0.25, -0.06, -0.34, -0.49, -0.72, -0.79, -0.94, -1.35, -0.35, 0.54, 0.68, 0.84, 1.20, 1.23, 1.32, 1.34];
y = [0.30, 0.43, 0.90, 1.40, 1.13, 1.08, 1.14, 1.23, 0.52, 0.21, -0.20, -0.73, -0.73, -0.82, -0.71, -0.76, -0.46, -0.13, 0.30];

给出闭合曲线（或多边形）：

figure(1)
hold on
plot(x,y,'k');
scatter(x,y,'r');
xlim([-2 2]);
ylim([-2 2]);
axis equal

我想计算曲线上各点的曲率（尽可能准确）。

到目前为止我所拥有的是切向量（一阶导数）和曲率（二阶导数）的简单计算：

dsx = diff(x);
dsy = diff(y);
ds = sqrt(dsx.^2+dsy.^2);
Tx = dsx./ds;
Ty = dsy./ds;
ds2 = 0.5*(ds(1:end-1)+ds(2:end));
Hx = diff(Tx)./ds2;
Hy = diff(Ty)./ds2;

但我得到的曲率非常不准确：

figure(1)
quiver(x(1:end-2),y(1:end-2),Hx,Hy,'b','autoscalefactor',1.2); 
xlim([-2 2]); ylim([-2 2]);
axis equal

这应该是一个简单的计算，但是它不起作用，请指教：如何找到最简单近似的曲率并且在方向和大小上具有合理的精度？

Answer 1

曲率计算正确，是绘图不对。请注意，diff 计算后续元素之间的差异，生成一个少了一个元素的向量。它估计样本对之间的导数。如果你重复这个，你会得到样本的二阶导数，但不是第一个或最后一个样本（你现在少了 2 个元素）。

您确实注意到了这一点，因为除了一个顶点之外，您正在绘制曲率。

所以您需要做的就是在一阶导数之后复制一个点（我将最后一个点添加到开头，因此元素的顺序与输入数组中的顺序相同）。索引语句 Tx([end,1:end]) 就是这样做的。

在下面的代码中，我还绘制了黑色的法线 (Ty,-Tx)。

x = [1.34, 0.92, 0.68, 0.25, -0.06, -0.34, -0.49, -0.72, -0.79, -0.94, -1.35, -0.35, 0.54, 0.68, 0.84, 1.20, 1.23, 1.32, 1.34];
y = [0.30, 0.43, 0.90, 1.40, 1.13, 1.08, 1.14, 1.23, 0.52, 0.21, -0.20, -0.73, -0.73, -0.82, -0.71, -0.76,-0.46, -0.13, 0.30];

% First derivative
dsx = diff(x);
dsy = diff(y);
ds = sqrt(dsx.^2+dsy.^2);
Tx = dsx./ds;
Ty = dsy./ds;

% Second derivative & curvature
ds2 = 0.5*(ds([end,1:end-1])+ds);
Hx = diff(Tx([end,1:end]))./ds2;
Hy = diff(Ty([end,1:end]))./ds2;

% Plot
clf
hold on
plot(x,y,'ro-');
x = x(1:end-1);
y = y(1:end-1); % remove repeated point
quiver(x+dsx/2,y+dsy/2,Ty,-Tx,'k','autoscalefactor',0.3); 
quiver(x,y,Hx,Hy,'b','autoscalefactor',1.2); 
set(gca,'xlim',[-2 2],'ylim',[-1.5 2]);
axis equal

Answer 2

严格来说，这个问题没有意义。

多边形在顶点处没有曲率，沿边的曲率为零。

您需要更具体地说明要评估的数量，例如解释用途。

在matlab中计算闭合曲线（或多边形）的曲率

Calculating the curvature of a closed curve ( or polygon) in matlab

matlab

computational-geometry