最近一直在做线结构光扫描三维成像方面的内容，采用结构光进行扫描检测时，需要提取激光条纹的中心线，我利用Steger算法提取激光中心线的亚像素坐标，在Matlab 2018b 软件上运行。
联系方式：QQ：936874728
邮箱：[email protected]
[email protected]
欢迎读者指正错误，希望和大家一起交流线结构光扫描三维成像相关的内容

下面先介绍一下steger算法的运行思路：(借鉴的别人的总结，直接截图白嫖过来了)

根据算法思路，我们
①输入图像
②图像处理
③利用高斯模板与图像卷积构造hessian矩阵
④求解hessian的特征值得到光条的法线方向
⑤法线方向泰勒展开获得亚像素位置
⑥选择合适的像素点，此时获得的是像素坐标（x0，y0），和对应的（tnx,tny）相加获得亚像素坐标。
根据上述思路编写matlab代码：（代码主题是嫖的网上的开源码，然后结合自己需要进行的修改）。

第一部分先是主函数啦
代码解释：
1.bmp是我采集的激光线的图片，会附在下面。
然后选取感兴趣的区域
最后plot函数绘图
保存激光中心线的亚像素坐标

img = imread('1.bmp');
 laserPixel = findLaserCenter(img,1,960,301,957,0);
 plot(laserPixel(:,1),laserPixel(:,2),'o');
 axis equal
 save C:\Users\17806\Desktop\Steger算法\1.txt -ascii laserPixel

第二部分是调用函数部分
①对感兴趣区域的图像处理
这一部分我就不解释了，都注释在代码里面了

function realPixel = findLaserCenter(frame, startRow, endRow, startColumn, endColumn, plotStatue)
if nargin < 2, plotStatue = 0;    %plotStatue,默认为0
end
frame_gray = rgb2gray(frame);%将文件的第i帧转换为灰度图像
frame_gray_roi = frame_gray(startRow:endRow, startColumn:endColumn);%感兴趣范围
[row, column]=size(frame_gray);     % 图像行列数
keral = fspecial('gaussian',[3,3],0.5);%fspecial函数用于创建预定义的滤波算子，格式：h=fspecial(type,parameters,sigma)
%gaussian,高斯低通滤波器，参数有两个，n表示模板尺寸，默认值为[3，3]，sigma为滤波器的标准差，单位为像素，默认值为0.5
frame_gray_gauss = imfilter(frame_gray_roi,keral,'replicate');%imfilter对任意类型数组或多维图像进行滤波，imfilter(区域，类型，边界)
%replicate图像大小通过复制外边界的值来扩展。
Pixel=steger(frame_gray_gauss);%对滤波之后的图像运用steger算法求得像素坐标。
realPixel = [ Pixel(:, 2) + startColumn - 1,  Pixel(:, 1) + startRow - 1];%真实像素坐标
%Generat a binary image of equal size to show the middle line of laser
%生成大小相同的二值图像，以显示激光器的中心线
if plotStatue ~= 0
    image_1=zeros(row,column);      % 等大的全黑背景
    [max, t] = size(Pixel);
    for index = 1 : max
        image_1(Pixel(index, 1) + startRow - 1, Pixel(index, 2) + startColumn - 1) = 255;
    end
    figure;
    imshow(image_1);
end
end

②调用steger函数
完全根据steger算法思路进行编写的，简单也不用解释了，看代码注释

function linePixel = steger(EinRDR)
%采用steger算法提取光条中心
 im = double(EinRDR);
 im = im/max(im(:)); %%归一化
 bw=im2bw(im, 0.4);  %%二值化
 bw = bwmorph(bw, 'clean', 1);    % 去孤立点
[x0,y0]=find(bw==1);  %%找到白点坐标
sigma=25/sqrt(3);%求Hessian矩阵之前需要对图像进行高斯滤波，高斯滤波时，设置高斯方差σ<（w/√3），w是光条宽度
[Dx,Dy,Dxx,Dxy,Dyy] = Hessian2D(im,sigma);
[eigenvalue1, eigenvalue2, eigenvectorx, eigenvectory]=eig2image(Dxx, Dxy, Dyy); 
%eigenvectorx,eigenvectory:Hessian矩阵最大特征值对应的特征向量对应于光条的法线方向 
t = -(Dx.*eigenvectorx + Dy .* eigenvectory) ./...  
    (Dxx .* eigenvectorx.^2 + 2*Dxy.*eigenvectorx.*eigenvectory + Dyy.*eigenvectory.^2 );   
px = t.*eigenvectorx;  
py = t.*eigenvectory;
[candidateX1, candidateY1] = find(px >= -0.5 & px <= 0.5 & py >= -0.5 & py <= 0.5 & bw==1);  
%判断：如果(px,py)∈[?0.5,0.5]×[?0.5,0.5]
%即一阶导数为零的点位于当前像素内，
%且(nx,ny)方向的二阶导数大于指定的阈值，
%则该点(candidataX1,candidataY1）为光条的中心点，candidataX1+px,candidataY1+py 则为亚像素坐标。
linePixel_t = [candidateX1, candidateY1];
for i=1:size(candidateX1,1)
    m1=candidateX1(i,1);
    n1=candidateY1(i,1);
    px1(i,1)=px(m1,n1);
    py1(i,1)=py(m1,n1);
    x1(i,1)=m1+px(m1,n1);
    x2(i,1)=py(m1,n1)+n1;
end
linePixel = [x1,x2];
%获得激光线像素数点坐标
end

③构造hessian矩阵
利用高斯模板和图像卷积获得图像沿x，y方向上的二阶导数

function [Dx,Dy,Dxx,Dxy,Dyy] = Hessian2D(I,Sigma)
%构造高斯模板
if nargin < 2, Sigma = 1; end
[X,Y]   = ndgrid(-round(3*Sigma):round(3*Sigma));
DGaussx  = 1/(2*pi*Sigma^4)*(-X).* exp(-(X.^2 + Y.^2)/(2*Sigma^2));
DGaussy  = 1/(2*pi*Sigma^4)*(-Y).* exp(-(X.^2 + Y.^2)/(2*Sigma^2));
DGaussxx = 1/(2*pi*Sigma^4) * (X.^2/Sigma^2 - 1) .* exp(-(X.^2 + Y.^2)/(2*Sigma^2));
DGaussxy = 1/(2*pi*Sigma^6) * (X .* Y)           .* exp(-(X.^2 + Y.^2)/(2*Sigma^2));
DGaussyy = DGaussxx';
%卷积
Dx  = imfilter(I,DGaussx,'conv');
Dy  = imfilter(I,DGaussy,'conv');
Dxx = imfilter(I,DGaussxx,'conv');
Dxy = imfilter(I,DGaussxy,'conv');
Dyy = imfilter(I,DGaussyy,'conv');
end

④计算黑森矩阵的特征值

function [Lambda1,Lambda2,Ix,Iy]=eig2image(Dxx,Dxy,Dyy)
tmp = sqrt((Dxx - Dyy).^2 + 4*Dxy.^2);
v2x = 2*Dxy; v2y = Dyy - Dxx + tmp;
% Normalize
%标准化
mag = sqrt(v2x.^2 + v2y.^2); i = (mag ~= 0);
v2x(i) = v2x(i)./mag(i);
v2y(i) = v2y(i)./mag(i);
% The eigenvectors are orthogonal
%实对称矩阵性质：不同特征值对应的特征向量是正交的
v1x = -v2y; 
v1y = v2x;
% Compute the eigenvalues
%计算特征值
mu1 = 0.5*(Dxx + Dyy + tmp);
mu2 = 0.5*(Dxx + Dyy - tmp);
% Sort eigen values by absolute value abs(Lambda1)<abs(Lambda2)
%按绝对值abs(Lambda1)<abs(Lambda2)对特征值排序
check=abs(mu1)>abs(mu2);
Lambda1=mu1; Lambda1(check)=mu2(check);
Lambda2=mu2; Lambda2(check)=mu1(check);
Ix=v1x; Ix(check)=v2x(check);
Iy=v1y; Iy(check)=v2y(check);
end

附上程序所需的图片：
我也不知道什么原因，可能网页原因或者我的网络原因吧，原图没传上去。。。。下面截图，大小发生了变化，读者还原过程时需要重新选择感兴趣区域，，，

试验结果：


部分坐标值：

可能叙述和代码还有些许问题，请大家加以指正，目前在做线结构光扫描三维成像相关的内容，有相关课题的可以加好友讨论。
QQ：936874728
邮箱：[email protected]
[email protected]