光条中心提取-steger代码

        Steger光条中心提取原理可查看之前的文章《光条中心提取-steger》，以下Steger光条中心提取为代码，该代码需修改图像路径，根据实际需求替换高斯核，替换完后即可正常运行。（欢迎进Q群交流：874653199）

clc;
clear;
%% 读取光条图片
srcImage=imread('C:\Users\admin\Desktop\test\6.jpg');

if(numel(size(srcImage))>2)
    srcImage=rgb2gray(srcImage);
end

%% 设置标准差构造高斯核及二值化
sigma=3;
thresh= graythresh(srcImage);%otsu  
srcImage=double(srcImage);

[m,n]=size(srcImage);
ky=[-1,1];
kx=[-1;1];
kyy=[1,-2,1];
kxx=[1;-2;1];
kxy=[1,-1;-1,1];


gausFilter = fspecial('gaussian',17,sigma);   %高斯滤波
dstImage=imfilter(srcImage,gausFilter,'replicate');  

dx=imfilter(dstImage,kx);
dy=imfilter(dstImage,ky);
dxx=imfilter(dstImage,kxx);
dyy=imfilter(dstImage,kyy);
dxy=imfilter(dstImage,kxy);

%% 计算Hessian矩阵
hessian=zeros(2,2);
points=zeros(m*n,2);

for i=1:m
    for j=1:n
        if(srcImage(i,j)/255>thresh)
            hessian(1,1)=dxx(i,j);
            hessian(1,2)=dxy(i,j);
            hessian(2,1)=dxy(i,j);
            hessian(2,2)=dyy(i,j);
            [eigenvectors,eigenvalues]=eig(hessian);
            
            if(abs(eigenvalues(1,1))>= abs(eigenvalues(2,2)))
            nx=eigenvectors(1,1);
            ny=eigenvectors(2,1);
            fmax_dist=eigenvalues(1,1);
            else
            nx=eigenvectors(1,2);
            ny=eigenvectors(2,2);
            fmax_dist=eigenvalues(2,2);
                
            end
            
            t=-(nx*dx(i,j)+ny*dy(i,j))/(nx*nx*dxx(i,j)+2 * nx*ny*dxy(i,j)+ny*ny*dyy(i,j));
            
            if(abs(t*nx) <= 0.5 && abs(t*ny) <= 0.5)
                points((i-1)*m+j,:)=[ j+t*ny,i+t*nx];
            end
            
            
        end
        
    end
end

index = find(points(:,1)==0);
points(index,:) = [];
index = find(points(:,2)==0);
points(index,:) = [];

figure(1);
imshow(uint8(srcImage));
title('Steger光条中心提取','color','b')
hold on
plot(points(:,1),points(:,2),'r.','MarkerSize',1)
hold off

Steger光条中心提取结果：