Steger 算法技术细节推导

面试问到了 steger 算法中非常细的几个问题：

• 为什么要用海森矩阵计算光条中心点？
• 往哪个方向展开海森矩阵？
• 光条的方向怎么计算得到的？

下面从算法原理讲起，回答这些问题。

光条提取究竟在解决什么问题？

对于一张包含激光条的图像，在激光条覆盖的像素邻域内，灰度分布有两个核心特征：

1. 垂直于光条的方向（法向 \(\boldsymbol{n}\)）：激光能量集中在光条中心，向法向 \(\boldsymbol{n}\) 两侧快速衰减，灰度 \(g\) 分布严格符合高斯型亮脊结构（中间亮、两侧暗），因此灰度变化最剧烈、二阶曲率最大，光条的真实中心就位于这个方向的灰度极值点。
2. 沿光条延伸的方向（切向 \(\boldsymbol{n}^{\perp}\)）：激光能量沿光条方向 \(\boldsymbol{n}^{\perp}\) 均匀分布，灰度变化平缓、二阶曲率极小，仅存在缓慢的亮度波动，不存在明显的极值点。

简单来说：光条是一条「沿切向延伸、法向呈高斯峰」的线状结构，我们的任务就是精准定位这条线的中心。

基于上述特性，光条提取需要解决两个核心问题：

方向检测：对每个像素，精准判断哪个方向是光条的法向（灰度变化最剧烈的方向），哪个是切向（灰度变化平缓的方向）；

亚像素定位：在法向方向上，精准计算出光条中心的亚像素坐标（精度可达 0.01~0.1 像素，远高于普通像素级定位），同时过滤噪声、反光、杂散光的干扰。

steger 算法是如何做的？

针对前面所要解决的两个核心问题，steger 算法是这么做的：

方向检测：

• 二阶导数描述局部灰度曲率：将图像考虑成一个离散二元函数，求某个像素沿着不同方向的灰度值变化剧烈程度，本质是求该方向的局部曲率，而二阶导数正是描述曲率的核心指标，这也是所谓的海森矩阵的核心作用。
• 二阶方向导数描述灰度曲率变化方向：在所有的灰度曲率变化方向中，只有垂直于光条方向的灰度值变化剧烈程度最大，且沿着光条方向灰度值变化剧烈程度最小，二者方向是相互垂直的。因此，在已知局部曲率的条件下，需要找出一个方向导数，使得该方向的曲率变化最大。而海森矩阵的特征值本质上就是对应特征向量方向上的二阶导数，绝对值越大，代表该方向的曲率变化越剧烈。于是在海森矩阵计算出的两个特征值中，绝对值较大的那个对应的特征向量就是光条法向，而绝对值较小（接近于0）的那个就是光条切向。

亚像素定位：

• 法向灰度分布获取：得到光条法向后，可获取沿该方向的灰度值分布，其分布近似高斯型亮脊，光条中心对应其灰度极值点所在位置。
• 泰勒展开求亚像素：图像坐标是整数离散的，要得到精确的极值点位置，需利用二阶泰勒展开对局部灰度做多项式近似，再通过令一阶导数为零，求解得到极值点的亚像素坐标。同时，会通过设置二阶方向导数阈值、约束亚像素偏移量在当前像素邻域内等方式，过滤噪声、反光、杂散光的干扰，保证定位的鲁棒性。

详细推导过程

海森矩阵计算

令包含激光条的灰度图像为 \(I(x,y)\) ，考虑任意位置的某个像素 \((x_0,y_0)\) ，其一阶偏导为 \(r_x(x_0,y_0)\)、\(r_y(x_0,y_0)\)，二阶偏导为 \(r_{xx}(x_0,y_0)\) 、\(r_{xy}(x_0,y_0)\)、\(r_{yx}(x_0,y_0)\)、\(r_{yy}(x_0,y_0)\)。因此该处的海森矩阵为：

\[H(x_0,y_0)=\left( \begin{array} {ll}{r_{xx}(x_0,y_0)}&{r_{xy}(x_0,y_0)}\\ {r_{yx}(x_0,y_0)}&{r_{yy}(x_0,y_0)}\end{array} \right) \]

注意：二阶混合偏导 \(r_{xy} \equiv r_{yx}\) 。

偏导计算

这里的一阶导数和二阶导数计算可以用使用高斯偏导核做卷积，可以获得更好的噪声抑制效果。并且通过调整高斯核的 \(\sigma\) 的大小，可以控制平滑尺度，适配不同粗细的光条（粗光条用大 \(\sigma\)，细光条用小 \(\sigma\)）。

二维高斯核可以表示为：

\[G(x,y,\sigma) = \frac{1}{2\pi\sigma^2} e^{-\frac{x^2+y^2}{2\sigma^2}} \]

图像的高斯卷积（后续实际没使用到）表示为：

\[r(x,y) = I(x,y)*G(x,y,\sigma) \]

于是图像的一阶偏导可以表示为：

\[\begin{cases} r_x(x,y) = \frac{\partial I}{\partial x} = I(x,y) * G_x(x,y,\sigma) \\ r_y(x,y) = \frac{\partial I}{\partial y} = I(x,y) * G_y(x,y,\sigma) \end{cases} \]

其中 \(G_x = \frac{\partial G}{\partial x}\)、\(G_y = \frac{\partial G}{\partial y}\) 是一阶高斯偏导核（也叫高斯导数滤波器）。

具体地：

\[\begin{align*} G_x(x,y,\sigma) &= \frac{\partial G}{\partial x} = -\frac{x}{2\pi\sigma^4} e^{-\frac{x^2+y^2}{2\sigma^2}} \\ G_y(x,y,\sigma) &= \frac{\partial G}{\partial y} = -\frac{y}{2\pi\sigma^4} e^{-\frac{x^2+y^2}{2\sigma^2}} \end{align*} \]

图像的二阶偏导可表示为：

\[\begin{cases} r_{xx}(x,y) = \frac{\partial^2 I }{\partial x^2} = I(x,y) * G_{xx}(x,y,\sigma) \\ r_{xy}(x,y) = \frac{\partial^2 I }{\partial x \partial y} = I(x,y) * G_{xy}(x,y,\sigma) \\ r_{yy}(x,y) = \frac{\partial^2 I }{\partial y^2} = I(x,y) * G_{yy}(x,y,\sigma) \end{cases} \]

其中 \(G_{xx}, G_{xy}, G_{yy}\) 是二阶高斯偏导核 。

具体地：

\[\begin{align*} G_{xx}(x,y,\sigma) &= \frac{\partial^2 G}{\partial x^2} = \frac{x^2 - \sigma^2}{2\pi\sigma^6} e^{-\frac{x^2+y^2}{2\sigma^2}} \\ G_{xy}(x,y,\sigma) &= \frac{\partial^2 G}{\partial x \partial y} = \frac{xy}{2\pi\sigma^6} e^{-\frac{x^2+y^2}{2\sigma^2}} \\ G_{yy}(x,y,\sigma) &= \frac{\partial^2 G}{\partial y^2} = \frac{y^2 - \sigma^2}{2\pi\sigma^6} e^{-\frac{x^2+y^2}{2\sigma^2}} \end{align*} \]

关键补充

归一化问题

高斯偏导核的积分不为1，因此卷积后图像的灰度会发生缩放，工程中需：

• 对核做归一化（保证核的积分和为1，避免图像整体变亮/变暗）；
• 或在后续计算中忽略绝对灰度，仅用相对值（如Steger算法中仅用偏导的比值、符号）。

核尺寸选择

高斯核是无限支撑的，工程中需截断为有限尺寸，通常取：

\[\text{核尺寸} = 2 \times \lceil 3\sigma \rceil + 1 \]

（\(3\sigma\) 原则：截断后核的能量保留99.7%，误差可忽略）

steger 论文指出，若光条中心到边缘的宽度为 \(w\) ，那么当 满足 \(\sigma \ge \frac{w}{\sqrt{3}}\) 时，线点中心处的一阶导数和二阶导数的取值特点才稳定存在。

可分离性加速

二维高斯偏导核是可分离的，可拆分为「一维行核 + 一维列核」，将卷积复杂度从 \(O(N^2)\) 降为 \(O(N)\)，大幅提升计算效率（Steger算法工程实现的核心优化点）。

海森矩阵特征值和特征向量计算

根据海森矩阵，求其特征值，具体为：

\[\begin{align*} \lambda_1 &= \frac{1}{2}(r_{xx}(x_0,y_0) + r_{yy}(x_0,y_0) + \sqrt{(r_{xx}(x_0,y_0)-r_{yy}(x_0,y_0))^2 + 4r_{xy}(x_0,y_0)^2}) \\ \lambda_2 &= \frac{1}{2}(r_{xx}(x_0,y_0) + r_{yy}(x_0,y_0) - \sqrt{(r_{xx}(x_0,y_0)-r_{yy}(x_0,y_0))^2 + 4r_{xy}(x_0,y_0)^2}) \end{align*} \]

对应的特征向量为：

\[\begin{align*} \boldsymbol{v_1} &= (v_{1x},v_{1y})^T \\ &=(r_{xx}(x_0,y_0) + r_{yy}(x_0,y_0) + \sqrt{(r_{xx}(x_0,y_0)-r_{yy}(x_0,y_0))^2 + 4r_{xy}(x_0,y_0)^2} , 2r_{xy}(x_0,y_0))^T \\ \boldsymbol{v_2} &= (v_{2x},v_{2y})^T \\ &=(r_{xx}(x_0,y_0) - r_{yy}(x_0,y_0) -\sqrt{(r_{xx}(x_0,y_0)-r_{yy}(x_0,y_0))^2 + 4r_{xy}(x_0,y_0)^2} , 2r_{xy}(x_0,y_0))^T \end{align*} \]

法向和切向确定

这里特征向量需要做单位化:

\[\begin{align*} {\hat{\boldsymbol{v}}_i} &= (v_{ix},v_{iy})^T \\ &= \frac{\boldsymbol{v_i}}{\|\boldsymbol{v_i}\|}\\ &= \frac{\boldsymbol{v_i}}{\sqrt{v_{ix}^2+v_{iy}^2}}\\ \end{align*} \]

根据特征值大小确定法向和切向，即

• 当 \(|\lambda_{1}| \gt |\lambda_{2}|\) 时，\(\boldsymbol{n} = \boldsymbol{\hat{v}_1}\)，\(\boldsymbol{n}^{\perp} = \boldsymbol{\hat{v}_2}\)；
• 当 \(|\lambda_{1}| \lt|\lambda_{2}|\) 时，\(\boldsymbol{n} = \boldsymbol{\hat{v}_2}\)，\(\boldsymbol{n}^{\perp} = \boldsymbol{\hat{v}_1}\)。

注意：每一个像素都可以求出这样的切向和法向，但是只有光条中心位置的沿着法向的曲率变化是最大的，而非光条中心位置曲率变化会明显更小一些，因此可以通过限制最大特征值（对应特征向量为最大曲率变化方向）必须高于某个阈值 \(m\) 来滤除这些非光条中心点。即 \(max(|\lambda_{1}|,|\lambda_{2}|) > m\) 。

沿法向做二阶泰勒展开

对于图像\(I(x,y)\)，在任意位置的某个像素 \((x_0,y_0)\) 处的二阶泰勒展开为：

\[{ I \left(x,y\right) \approx \\ r(x_0,y_0) + \begin{pmatrix} {x-x_0} & {y-y_0} \end{pmatrix} \begin{pmatrix} r_x(x_0,y_0) \\ r_y(x_0,y_0)\end{pmatrix} +\frac{1}{2} \begin{pmatrix} {x-x_0} & {y-y_0}\end{pmatrix} \begin{pmatrix} r_{xx}(x_0,y_0) & r_{xy}(x_0,y_0) \\ r_{xy}(x_0,y_0) & r_{yy}(x_0,y_0)\end{pmatrix} \begin{pmatrix} {x-x_0} \\ {y-y_0} \end{pmatrix} } \]

光条法向 \(\boldsymbol{n} = (n_x,n_y)^T\) ，那么沿着这方向的任意一点 \((x,y)\) 应满足：

\[\begin{align*} x-x_0 &= tn_x \\ y-y_0 &= tn_y \end{align*} \]

其中 \(t\) 是一个可变参数。

将上式代入二阶泰勒展开，即可得到沿着法向的二阶泰勒展开：

\[{ r_n \left(x,y\right) \approx \\ r(x_0,y_0) + \begin{pmatrix} {tn_x} & {tn_y} \end{pmatrix} \begin{pmatrix} r_x(x_0,y_0) \\ r_y(x_0,y_0)\end{pmatrix} +\frac{1}{2} \begin{pmatrix} {tn_x} & {tn_y}\end{pmatrix} \begin{pmatrix} r_{xx}(x_0,y_0) & r_{xy}(x_0,y_0) \\ r_{xy}(x_0,y_0) & r_{yy}(x_0,y_0)\end{pmatrix} \begin{pmatrix} {tn_x} \\ {tn_y} \end{pmatrix} } \]

极值点亚像素坐标计算

令 $ I'_n \left(x,y\right)=0$ ，则有：

\[{ I'_n \left(x,y\right) \approx n_{x} r_x(x_0,y_0) + n_{y} r_y(x_0,y_0) + t \left( n^{2}_{x} r_{xx}(x_0,y_0) + 2n_{x}n_{y} r_{xy}(x_0,y_0)+ n^{2}_{y} r_{yy}(x_0,y_0) \right)=0 } \]

于是可变参数 \(t\) 的表达式为:

\[t=-\frac{r_{x}(x_0,y_0) n_{x}+r_{y}(x_0,y_0) n_{y}}{n_{x}^{2}r_{xx}(x_0,y_0)+2n_{x} n_{y}r_{xy}(x_0,y_0) +n_{y}^{2}r_{yy}(x_0,y_0)} \]

于是精确的极值点坐标为 \((x_0+tn_x, y_0+tn_y)\) 。因为像素是离散型变量，即使求精确解也不应该超过当前像素所在的区域，因此若 \((tn_x,tn_y)\in [-\frac{1}{2}, \frac{1}{2}] \times [-\frac{1}{2}, \frac{1}{2}]\) ，则该极值点为线条中心点。

后续步骤

上述内容只涉及到光条中心点的提取（或者称作线点检测），后续还有线点连接、线宽检测、偏差消除等一系列操作。在此不做详细推导。详见 https://www.cnblogs.com/gshang/p/18864955/steger_note 。

意外发现

python 有 开源的 ridge-detector 库，安装：

pip install ridge-detector

使用案例

from ridge_detector import RidgeDetector
import numpy as np
det = RidgeDetector(
    line_widths=np.array([3, 5, 9]),  # Line widths to detect
    low_contrast=50,  # Lower bound of intensity contrast, decrease this value if ridges are missed out
    high_contrast=128,  # Higher bound of intensity contrast, decrease this value if ridges are missed out
    min_len=10, # Ignore ridges shorter than this length
    max_len=0, # Ignore ridges longer than this length, set to 0 for no limit
    dark_line=False, # Set to True if detecting black ridges in white background, False otherwise
    estimate_width=True, # Estimate width for each detected ridge point
    extend_line=True, # Tend to preserve ridges near junctions if set to True
    correct_pos=False,  # Correct ridge positions with asymmetric widths if set to True
)
det.detect_lines("laser.png")
det.show_results()
det.save_results("./")  # Comment out if you want to save the detection results

源码可查看，这位作者代码写的很不错，基本符合 steger 算法原论文，并且还做了光条宽度多尺度自适应

def detect_lines(self, image):
    # 兼容图像与图像路径读取图像数据
    self.image = iio.imread(image) if isinstance(image, str) else image

    # 图像数据转换到 uint8 格式, 如果有需要
    if self.image.dtype != np.uint8:
        self.image = ((image - image.min()) / (image.max() - image.min()) * 255).astype(np.uint8)

    # 转成灰度图
    self.gray = cv2.cvtColor(self.image, cv2.COLOR_RGB2GRAY) if self.image.ndim == 3 else self.image

    # 图像求导
    self.apply_filtering()

    # 线点计算
    self.compute_line_points()

    # 线点连接
    self.compute_contours()

    # 线宽计算
    if self.estimate_width:
        self.compute_line_width()

    # 偏差消除
    self.prune_contours()