线条中心提取(Steger)

Steger（线条中心提取算法）

听风者868

4734人浏览 · 2022-05-10 09:13:17

听风者868 · 2022-05-10 09:13:17 发布

Steger（线条中心提取算法）

1、算法原理

Steger算法基于Hessian矩阵，能够实现对线条中心的亚像素提取。对于图像上任一点(x,y)，Hessian矩阵可以表示为：
$H(x,y)=\left[ \begin{matrix} h_{xx} & h_{xy}\\ h_{xy} & h_{yy}\\ \end{matrix} \right]$
该矩阵由图像某点沿不同方向的四个偏导数组成，该矩阵可以确定z(x,y)的二阶方向导数最大绝对值的方向。
Steger提取线条中心的步骤如下：
①高斯平滑；
②基于Hessian矩阵得到该像素点的法向；
③法向上泰勒展开拟合亚像素坐标。
Hessian 矩阵最大特征值对应的特征向量对应于线条的法线方向，用 $n_x,n_y)$ 表示，以点 $x_0,y_0)$ 为基准点，
图像上每个点都可以计算 Hessian矩阵，由于Hessian矩阵的特征值就是形容其在该点附近特征向量方向的凹凸性,特征值越大,凸性越强。线条中心点的Hessian矩阵特征值相比其他位置更大，所以可以设置阈值进行筛选，线条中心的亚像素坐标为：
$(p_x,p_y) = (x_0 + tn_x,y_0+tN_y) \\其中，\\ t=-\frac{r_xn_x+r_yn_y}{r_{xx}n_x^2+2r_{xy}n_xn_y+r_{yy}n_y^2} \\若存在点满足\\ (tn_x,tn_y)\epsilon[-0.5,0.5]*[-0.5,0.5]$
即一阶导数为零的点位于当前像素内（找拐点附近的点，因为实际中心截面不是规则的曲线），且方向的二阶导数大于指定的阈值（确定是极值点），则该点 $x_0,y_0)$ 为线条的中心点， $p_x,p_y)$ 则为亚像素坐标。 $tn_x,tn_y)$ 是中心点与亚像素的偏移值，所以其绝对值不能超过0.5，超过0.5就超出该像素的范围。

2、参考代码

void StegerCV()
{ 
    Mat srcImg = imread("srcImg.png", 1);
    Mat grayImg;
    cvtColor(srcImg, grayImg, CV_BGR2GRAY);

    // 高斯滤波
    grayImg.convertTo(grayImg, CV_32FC1);
    GaussianBlur(grayImg, grayImg, Size(0, 0), 6, 6);

    // 一阶偏导数
    Mat m1, m2;
    m1 = (Mat_<float>(1, 2) << 1, -1);  // x偏导
    m2 = (Mat_<float>(2, 1) << 1, -1);  // y偏导

    Mat dx, dy;
    filter2D(grayImg, dx, CV_32FC1, m1);
    filter2D(grayImg, dy, CV_32FC1, m2);

    //二阶偏导数
    Mat m3, m4, m5;
    m3 = (Mat_<float>(1, 3) << 1, -2, 1);        // 二阶x偏导
    m4 = (Mat_<float>(3, 1) << 1, -2, 1);        // 二阶y偏导
    m5 = (Mat_<float>(2, 2) << 1, -1, -1, 1);    // 二阶xy偏导

    Mat dxx, dyy, dxy;
    filter2D(grayImg, dxx, CV_32FC1, m3);
    filter2D(grayImg, dyy, CV_32FC1, m4);
    filter2D(grayImg, dxy, CV_32FC1, m5);

    // hessian矩阵
    double maxD = -1;
    int imgcol = img.cols;
    int imgrow = img.rows;
    vector<double> Pt;
    for (int i=0;i<imgcol;i++)
    {
        for (int j=0;j<imgrow;j++)
        {
            if (img0.at<uchar>(j,i)>200)
            {
                Mat hessian(2, 2, CV_32FC1);
                hessian.at<float>(0, 0) = dxx.at<float>(j, i);
                hessian.at<float>(0, 1) = dxy.at<float>(j, i);
                hessian.at<float>(1, 0) = dxy.at<float>(j, i);
                hessian.at<float>(1, 1) = dyy.at<float>(j, i);

                Mat eValue;
                Mat eVectors;
                eigen(hessian, eValue, eVectors);

                double nx, ny;
                double fmaxD = 0;
                if (fabs(eValue.at<float>(0,0))>= fabs(eValue.at<float>(1,0)))  //求特征值最大时对应的特征向量
                {
                    nx = eVectors.at<float>(0, 0);
                    ny = eVectors.at<float>(0, 1);
                    fmaxD = eValue.at<float>(0, 0);
                }
                else
                {
                    nx = eVectors.at<float>(1, 0);
                    ny = eVectors.at<float>(1, 1);
                    fmaxD = eValue.at<float>(1, 0);
                }

                double t = -(nx*dx.at<float>(j, i) + ny*dy.at<float>(j, i)) / (nx*nx*dxx.at<float>(j,i)+2*nx*ny*dxy.at<float>(j,i)+ny*ny*dyy.at<float>(j,i));

                if (fabs(t*nx) <= 0.5 && fabs(t*ny) <= 0.5)
                {
                    double x_sub = i + t*nx;
                    double y_sub = j + t*ny;
                    Pt.push_back(x_sub);
                    Pt.push_back(x_sub);
                }
            }
        }
    }

    for (int k = 0;k<Pt.size()/2;k++)
    {
        Point rpt;
        rpt.x = Pt[2 * k + 0];
        rpt.y = Pt[2 * k + 1];
        circle(srcImg, rpt, 1, Scalar(0, 0, 255));
    }

    imshow("result", srcImg);
    waitKey(0);
}