光条中心线提取-Steger算法（基于Hessian矩阵）类似于ArcScan提取道路中心线

ゝ一纸荒年。 2023-03-13 02:45 6阅读 0赞

采用结构光进行扫描检测时，需要提取激光条纹的中心线，本文采用经典的Steger算法提取光条中心。

Steger算法原理  
Steger算法基于Hessian矩阵，能够实现光条中心亚像素精度定位：首先通过Hessian矩阵能够得到光条的法线方向，然后在法线方向利用泰勒展开得到亚像素位置。  
对于图像中激光条纹上的任意一点(x,y)(x,y)，Hessian矩阵可以表示为：  
H(x,y)=\[rxxrxyrxyryy\]  
H(x,y)=\[rxxrxyrxyryy\]

其中rxxrxx表示图像沿x的二阶偏导数，其他参数类似。需要注意的是在求Hessian矩阵之前需要对图像进行高斯滤波，高斯滤波时，根据文献\[1\]中，设置高斯方差σ<w3√σ<w3，其中ww为光条宽度。Hessian矩阵最大特征值对应的特征向量对应于光条的法线方向，用(nx,ny)(nx,ny) 表示，以点(x0,y0)(x0,y0) 为基准点，则光条中心的亚像素坐标为：  
(px,py)=(x0+tnx，y0+tny)(px,py)=(x0+tnx，y0+tny)  
公式中  
t=−nxrx+nyryn2xrxx+2nxnyrxy+n2yryyt=−nxrx+nyrynx2rxx+2nxnyrxy+ny2ryy  
如果(tnx,tny)∈\[−0.5,0.5\]×\[−0.5,0.5\](tnx,tny)∈\[−0.5,0.5\]×\[−0.5,0.5\] ，即一阶导数为零的点位于当前像素内，且(nx,ny)(nx,ny) 方向的二阶导数大于指定的阈值，则该点(x0,y0)(x0,y0) 为光条的中心点，(px,py)(px,py) 则为亚像素坐标。  
示例代码  
根据Steger算法的原理，借助opencv的Mat数据结构，实现Steger算法如下：

void StegerLine()
    { 
        Mat img0 = imread("image_0.png", 1);
        Mat img;
        cvtColor(img0, img0, CV_BGR2GRAY);
        img = img0.clone();
    
        //高斯滤波
        img.convertTo(img, CV_32FC1);
        GaussianBlur(img, img, 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(img, dx, CV_32FC1, m1);
        filter2D(img, 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(img, dxx, CV_32FC1, m3);
        filter2D(img, dyy, CV_32FC1, m4);
        filter2D(img, 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)
                    {
                        Pt.push_back(i);
                        Pt.push_back(j);
                    }
                }
            }
        }
    
        for (int k = 0;k<Pt.size()/2;k++)
        {
            Point rpt;
            rpt.x = Pt[2 * k + 0];
            rpt.y = Pt[2 * k + 1];
            circle(img0, rpt, 1, Scalar(0, 0, 255));
        }
    
        imshow("result", img0);
        waitKey(0);
    }

### 示例结果 ###

![这里写图片描述][SouthEast]

参考博客：

其中测绘可能会用到的工具ArcToolBox进行道路中心线提取：[ArcMap：矢量线数据提取（道路）中心线][ArcMap]

[SouthEast]: /images/20230312/d1914dacbff44505ac00b0b034c12681.png
[ArcMap]: https://blog.csdn.net/xinjiang666/article/details/78481059