Steger算法提取曲线中心线

Posted 2023-05-15 孟祥度与小岱的博客

找中心线在机器视觉上的应用可以说几乎没有，那我为什么要写这篇文章呢？主要原因是这种方法和这一篇文章《9点拟合梯度边缘亚像素方法 - 兜尼完 - 博客园 (cnblogs.com)》功能是相近的。都是根据一个整数像素点拟合一个亚像素值，都是求一个二元函数在某个方向上的极值点，但是方法却不一样。本文侧重理论上的推导。多年前我曾验证过这个算法，但是那些代码早就没了。原始论文链接：An Unbiased Detector of Curvilinear Structures (tum.de)。相关知识点参见：

• 知乎的介绍：基于Hessian矩阵的光条中心提取算法 - 知乎 (zhihu.com)
• 二元泰勒展开的介绍：Steger算法原理详解_Jieckiee的博客-CSDN博客

已知条件是整数像素点$ \\left( x_0,y_0 \\right) $。将图像看成一个二元函数$ f \\left( x,y \\right) $，那么在这点附近可进行泰勒展开。为了突出海森矩阵，我们把它写成矩阵形式：

$$ f \\left( x_0+\\Delta x,y_0+\\Delta y \\right) \\approx f \\left( x_0,y_0 \\right) + \\beginpmatrix \\Delta x & \\Delta y \\endpmatrix \\beginpmatrix f^\'_x \\left( x_0,y_0 \\right) \\\\ f^\'_y \\left( x_0,y_0 \\right) \\endpmatrix+\\frac12 \\beginpmatrix \\Delta x & \\Delta y \\endpmatrix \\beginpmatrix f^\'\'_xx \\left( x_0,y_0 \\right) & f^\'\'_xy \\left( x_0,y_0 \\right) \\\\ f^\'\'_xy \\left( x_0,y_0 \\right) & f^\'\'_yy \\left( x_0,y_0 \\right) \\endpmatrix \\beginpmatrix \\Delta x \\\\ \\Delta y \\endpmatrix $$

我们并不需要知道函数$ f \\left( x,y \\right) $的具体形式。从上式可知只要知道它的一阶、二阶偏导数就可以对它进行处理。而图像的一阶、二阶偏导数确实可以通过卷积运算得到。接下来就是确定拟合亚像素的方向。论文是利用海森矩阵的特征值和特征向量取得这个方向的：海森矩阵较大的那个特征值对应的特征向量是函数数值变化率最大的方向；较小的那个特征值对应的特征向量是函数数值变化率最小的方向。因此这里使用的是较大特征值对应的特征向量$ \\left( e_x,e_y \\right) $（可单位化此向量）。计算特征值、特征向量可以用OpenCV中的cv::eigen(...)函数。我们对上式加上方向限制，令$ \\Delta x= te_x,\\Delta y=te_y $，则有：

$$ f \\left( x_0+te_x,y_0+te_y \\right) \\approx f \\left( x_0,y_0 \\right) + \\beginpmatrix te_x & te_y \\endpmatrix \\beginpmatrix f^\'_x \\left( x_0,y_0 \\right) \\\\ f^\'_y \\left( x_0,y_0 \\right) \\endpmatrix+\\frac12 \\beginpmatrix te_x & te_y \\endpmatrix \\beginpmatrix f^\'\'_xx \\left( x_0,y_0 \\right) & f^\'\'_xy \\left( x_0,y_0 \\right) \\\\ f^\'\'_xy \\left( x_0,y_0 \\right) & f^\'\'_yy \\left( x_0,y_0 \\right) \\endpmatrix \\beginpmatrix te_x \\\\ te_y \\endpmatrix $$

可以看出来上式已经变成了一个关于t的一元函数。要计算它的极值只需要对t求导令其为0即可：

$$ f^\' \\left( x_0+te_x,y_0+te_y \\right) \\approx e_x f^\'_x \\left( x_0,y_0 \\right) + e_y f^\'_y \\left( x_0,y_0 \\right) + t \\left( e^2_x f^\'\'_xx \\left( x_0,y_0 \\right) + 2e_xe_yf^\'\'_xy \\left( x_0,y_0 \\right) + e^2_y f^\'\'_yy \\left( x_0,y_0 \\right) \\right)=0 $$

即可解得：

$$ t=\\frace_x f^\'_x \\left( x_0,y_0 \\right) + e_y f^\'_y \\left( x_0,y_0 \\right)e^2_x f^\'\'_xx \\left( x_0,y_0 \\right) + 2e_xe_yf^\'\'_xy \\left( x_0,y_0 \\right) + e^2_y f^\'\'_yy \\left( x_0,y_0 \\right) $$

因此，亚像素值为：

$$ \\left\\ \\beginmatrix x_sub=x_0+te_x \\\\ y_sub=y_0+te_y \\endmatrix \\right. , 需要te_x\\in \\left[ 0.5,0.5 \\right], te_y\\in \\left[ 0.5,0.5 \\right] $$

理论如上。可以看出来，利用海森矩阵求二元函数在某个方向上的极值点比用两曲面方程列方程组然后消元求解要方便一点。因它不需要区分特殊情况。

MATLAB | 垂距法提取离散坐标数据特征点(矢量曲线压缩)

垂距法是指根据中间顶点到其前、后两相邻顶点连线的距离的大小，来确定是否保留该顶点的一种 线要素顶点压缩 算法。当求得的距离大于给定的限差（阈值）时，保留该顶点，否则删除该顶点(如下图所示)。

一般使用所有点到直线距离的 中位数 作为阈值。同时，本人所编写的工具函数不仅可以处理二维数据点，三维甚至更高维度的数据点依旧可以处理。

PART.0 工具函数

怕大家找不到工具函数，这里放在最前面啦：

function [newPntSet,vertPnt]=getFeaturePnt(pntSet)
% @author:slandarer
% newPntSet : 特征点
% vertPnt   : 原始曲线垂足
[rows,cols]=size(pntSet);
if rows<4||cols<2,error('数据点过少或维度异常');end

innerPntSet=pntSet(2:end-1,:);             % 内部点
adjPnt1V=innerPntSet-pntSet(1:end-2,:);    % 当前点与前一点向量
adjPnt2V=pntSet(3:end,:)-pntSet(1:end-2,:);% 当前点的两个相邻点向量
adjPnt2V=adjPnt2V./vecnorm(adjPnt2V')';    % 归一化

lVert=sum(adjPnt1V.*adjPnt2V,2);   
vertPnt=adjPnt2V.*lVert+pntSet(1:end-2,:); % 获取垂足
vertL=vecnorm((innerPntSet-vertPnt)');     % 计算垂线距离
vertL(isinf(vertL)|isnan(vertL))=0;        % 修正/0的情况
innerPntSet(vertL<median(vertL),:)=[];     % 删掉特征性不强的点
newPntSet=[pntSet(1,:);innerPntSet;pntSet(end,:)];
end

PART.1 基础使用(二维)

这里随机生成一组二维数据(两列)，取出特征点后并绘图：

% 随机构造数据
X=linspace(0,25,10)';
Y=randi([0,10],[10,1]);
pntSet=[X,Y];

% 获取特征点
[nPntSet,vertPnt]=getFeaturePnt(pntSet);

% 坐标区域修饰
hold on
ax=gca;
ax.DataAspectRatio=[1,1,1];

% 绘制原始数据曲线
plot(pntSet(:,1),pntSet(:,2),'Color',[0 0.4470 0.7410],'LineWidth',2,'Marker','*');
% 绘制新数据曲线
plot(nPntSet(:,1),nPntSet(:,2),'Color',[0.6350 0.0780 0.1840 .7],'LineWidth',2,'Marker','s');

legend('original-curve','feature-curve')

PART.2 使用并修饰绘图(二维)

就把辅助线也画上并各种加属性呗，需要注意的是，垂足并不一定在相邻两点的连线上，有时候是在其延长线上，但因为怕麻烦这里延长线就没有画：

% 随机构造数据
X=linspace(0,25,10)';
Y=randi([0,10],[10,1]);
pntSet=[X,Y];

% 获取特征点
[nPntSet,vertPnt]=getFeaturePnt(pntSet);

% 坐标区域修饰
hold on
ax=gca;
ax.YLim=[0,10];
ax.DataAspectRatio=[1,1,1];
ax.Color=[1,1,1];
ax.XColor=[1,1,1].*.3;
ax.YColor=[1,1,1].*.3;
ax.LineWidth=1.5;
ax.FontName='cambria';

% 绘制原始数据曲线
l1=plot(pntSet(:,1),pntSet(:,2),'Color',[0 0.4470 0.7410],'LineWidth',2,'Marker','*');

% 绘制辅助线及垂线
innerPntSet=pntSet(2:end-1,:);
l2=plot([innerPntSet(:,1),vertPnt(:,1)]',[innerPntSet(:,2),vertPnt(:,2)]','Color',[.3,.3,.3],'LineWidth',1.2,'LineStyle','-.');
plot([pntSet(3:end,1),pntSet(1:end-2,1)]',[pntSet(3:end,2),pntSet(1:end-2,2)]','Color',[.3,.3,.3],'LineWidth',1.2,'LineStyle','-.')

% 绘制新数据曲线
l3=plot(nPntSet(:,1),nPntSet(:,2),'Color',[0.6350 0.0780 0.1840 .7],'LineWidth',2,'Marker','s');


% 增添图例
legend([l1,l2(1),l3],'original-curve','auxiliary-line','feature-curve');

PART.3 三维数据

也是几乎完全相同的使用方式，不过绘图用的plot3：

% 随机构造数据
X=linspace(0,25,10)';
Y=randi([0,10],[10,1]);
Z=randi([0,10],[10,1]);
pntSet=[X,Y,Z];

% 获取特征点
[nPntSet,vertPnt]=getFeaturePnt(pntSet);

% 坐标区域修饰
hold on
grid on
ax=gca;
ax.YLim=[0,10];
ax.ZLim=[0,10];
ax.DataAspectRatio=[1,1,1];
ax.Color=[1,1,1];
ax.XColor=[1,1,1].*.3;
ax.YColor=[1,1,1].*.3;
ax.ZColor=[1,1,1].*.3;
ax.LineWidth=1.5;
ax.FontName='cambria';

% 绘制原始数据曲线
l1=plot3(pntSet(:,1),pntSet(:,2),pntSet(:,3),'Color',[0 0.4470 0.7410],'LineWidth',2,'Marker','*');

% 绘制辅助线
innerPntSet=pntSet(2:end-1,:);
l2=plot3([innerPntSet(:,1),vertPnt(:,1)]',[innerPntSet(:,2),vertPnt(:,2)]',...
    [innerPntSet(:,3),vertPnt(:,3)]','Color',[.3,.3,.3],'LineWidth',1.2,'LineStyle','-.');
plot3([pntSet(3:end,1),pntSet(1:end-2,1)]',[pntSet(3:end,2),pntSet(1:end-2,2)]',...
    [pntSet(3:end,3),pntSet(1:end-2,3)]','Color',[.3,.3,.3],'LineWidth',1.2,'LineStyle','-.')

% 绘制新数据曲线
l3=plot3(nPntSet(:,1),nPntSet(:,2),nPntSet(:,3),'Color',[0.6350 0.0780 0.1840 .7],'LineWidth',2,'Marker','s');

% 增添图例
legend([l1,l2(1),l3],'original-curve','auxiliary-line','feature-curve');
view(3)