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

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

m=0;n=0;m1=0;n1=0;
for i=1:size(im,1)
    for j=1:size(im,2)
        m=m+i*im(i,j);m1=m1+im(i,j);
        n=n+j*im(i,j);n1=n1+im(i,j);
    end
end
m0=m/m1;
n0=n/n1;
m00=0;n00=0;m10=0;n10=0;
for i=1:size(im,1)
    for j=1:size(im,2)
        m00=m00+i*im(i,j).^2;m10=m10+im(i,j).^2;
        n00=n00+j*im(i,j).^2;n10=n10+im(i,j).^2;
    end
end
m000=m00/m10;
n000=n00/n10;
figure
imshow(im);hold on;
plot(m0,n0,'x');hold on;%质心法
plot(m000,n000,'+');hold on;%加权质心法

% 下面用hessian矩阵实现
Sigma =3/sqrt(3);
[a,b,c,d,e]=Hessian2D(im,Sigma);
for i=1:size(im,1)
    for j=1:size(im,2)
            [V,D]=eig([c(i,j),d(i,j);d(i,j),e(i,j)]);
            im2(i,j)=D(1,1)*D(2,2);
    end
end
im3=im2;
for ii=1:60
    [X Y] = find (im3 == max(max(im3)));
for nnn=1:size(X,1)
    i=X(nnn,1);j=Y(nnn,1);
s=(b(i,j)*d(i,j)-a(i,j)*e(i,j))/(c(i,j)*e(i,j)-d(i,j)^2);
t=(a(i,j)*d(i,j)-b(i,j)*c(i,j))/(c(i,j)*e(i,j)-d(i,j)^2);
if abs(s)<=1/2&&abs(t)<=1/2
    sss(1,1)=i+s;sss(1,2)=j+t;
plot(sss(:,1),sss(:,2),'*');hold on;
break
else
im3(X,Y)=0;
end
end
end
f1=0;f2=0;f3=0;g=0;%%判断标准
for i=1:60
    for j=1:60
if im(i,j)>200/255
    f3=f3+sqrt((i-sss(1))^2+(j-sss(2))^2);
        f1=f1+sqrt((i-m0)^2+(j-n0)^2);    
        f2=f2+sqrt((i-m000)^2+(j-n000)^2);
    g=g+1;
end
    end
end
lv1=f1/g;
lv2=f2/g;
lv3=f3/g;
% 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 = DGaussxx';
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