基于hessian矩阵的线条中心亚像素提取（matlab）

clear ;close ;
file='';%%自己加地址
im=imread(file);

im=im2double(im);
imshow(im);hold on;

%采用hessian法进行检测
Sigma =4;
[a,b,c,d,e]=Hessian2D(im,Sigma);
n=1;
for i=1:size(im,1)
    for j=1:size(im,2)

        if im(i,j)>200/255
            [V,D]=eig([c(i,j),d(i,j);d(i,j),e(i,j)]);
  if abs(D(1, 1)) >= abs(D(2, 2))
       norm = V(:, 1) ;
       lamda1 = D(2, 2);        
       lamda2 = D(1, 1);     
    else
       norm = V(:, 2) ; 
       lamda1 = D(1, 1);     
       lamda2 = D(2, 2);     
    end
    lamda = [lamda1; lamda2];  
   
    t = -(norm(1)*a(i,j) + norm(2)*b(i,j))/ ...
         (norm(1)^2*c(i,j) + 2*norm(1)*norm(2)*d(i,j) + norm(2)^2*e(i,j)); 

    if abs(t*norm(1))<=1/2&&abs(t*norm(2))<=1/2        

        sss(n,1)=i+t*norm(1);sss(n,2)=j+t*norm(2);
          n=n+1;

    end 

    end
end
end

plot(sss(:,2),sss(:,1),'*');

%hessian子函数
function [Dx,Dy,Dxx,Dxy,Dyy] = Hessian2D(I,Sigma)
[X,Y]   = ndgrid(-round(3*Sigma):round(3*Sigma));
Dx=( - X / (Sigma*Sigma)).*exp(-1 * (X.^2 + Y.^2) / (2 * Sigma*Sigma))*(1 / (2 * pi*Sigma^2));
Dy=( - Y / (Sigma*Sigma)).*exp(-1 * (X.^2 + Y.^2) / (2 * Sigma*Sigma))*(1 / (2 * pi*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 = 1/(2*pi*Sigma^4) * (Y.^2/Sigma^2 - 1) .* exp(-(X.^2 + Y.^2)/(2*Sigma^2));
Dx = imfilter(I,Dx,'conv');
Dy = imfilter(I,Dy,'conv');
Dxx = imfilter(I,DGaussxx,'conv');
Dxy = imfilter(I,DGaussxy,'conv');
Dyy = imfilter(I,DGaussyy,'conv');
end