sr_ridge.m

读入一副图片

%% Read image from file 
inImg = imread('C:\Documents and Settings\cc\桌面\毕业设计\retinal images\RI1五点中值滤波.bmp');
% inImg = imread('E:\defect images\w1.jpg');


%figure,imshow(inImg);

% inImg = double(inImg);
% inImg = im2double(rgb2gray(inImg));
inImg = double(inImg);

[m n] = size(inImg);

inImg1 = imresize(inImg, [64, 64], 'bilinear');
%figure,imshow(inImg);

%% Spectral Residual
myFFT = fft2(inImg1); 
% myAmplitude = abs(myFFT);
%figure,imshow(myAmplitude,[ ]);

myLogAmplitude = log(abs(myFFT));
myPhase = angle(myFFT);
% figure,imshow(myLogAmplitude,[ ]);
% figure,imshow(myPhase);


mySpectralResidual = myLogAmplitude - imfilter(myLogAmplitude, fspecial('average', 3), 'replicate'); 
% figure,imshow(mySpectralResidual);
saliencyMap = abs(ifft2(exp(mySpectralResidual + i*myPhase))).^2; %复数操作
% saliency  n.突起

%% After Effect
saliencyMap = imfilter(saliencyMap, fspecial('disk', 3));

MeanSaliency = mean(saliencyMap(:));
Threshold = MeanSaliency*3;
Threshold = 0.2;

saliencyMap = mat2gray(saliencyMap);
saliencyMap = imresize(saliencyMap, [m, n], 'bilinear');
% figure,imshow(saliencyMap);


objectMap = zeros(m,n);
inImg1 = zeros(m,n);
for i=1:m
    for j=1:n
        if saliencyMap(i,j) >Threshold
            objectMap(i,j) = 255;
            inImg1(i,j) = inImg(i,j);
        end
    end
end
% objectMap= uint8(objectMap);
% figure,imshow(objectMap);

inImg1= uint8(inImg1);
%figure,imshow(inImg1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% gradient image
inImg= double(inImg);
[Ix Iy] = gradient(inImg);
imImg = Ix.^2 + Iy.^2;

% imImg = mat2gray(imImg);
% figure,imshow(imImg);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ridge detector
[m n] = size(inImg);

[Ix Iy] = gradient(inImg);
[Ixx Ixy] = gradient(Ix);
[Iyx Iyy] = gradient(Iy);

vec = zeros(m,n,2);
val = zeros(m,n);
Q = zeros(m,n);

for i=1:m
    for j=1:n
        
        if objectMap(i,j)>0
            H(1,1) = Ixx(i,j);
            H(1,2) = Ixy(i,j);
            H(2,1) = Ixy(i,j);
            H(2,2) = Iyy(i,j);
            [v d] = eig(H);
        
            if abs(d(1,1))>abs(d(2,2))
                vec(i,j,1) = v(1,1);
                vec(i,j,2) = v(2,1);
                val(i,j) = d(1,1);
            else
                vec(i,j,1) = v(1,2);
                vec(i,j,2) = v(2,2);
                val(i,j) = d(2,2);
            end
        
             Q(i,j) = vec(i,j,1)*Ix(i,j) + vec(i,j,2)*Iy(i,j);
        end
        
    end
end

RidgeMap = ones(m,n)*255;

for i=2:m-1
    for j=2:n-1
        
        if val(i,j)<0
            for k=-1:1
                for l=-1:1
                    %%%%%%%%%%%%%%%%
                    if k~=0 && l~=0
                        Qy = Q(i+k,j+l);
                        P1 = vec(i,j,1)*vec(i+k,j+l,1) + vec(i,j,2)*vec(i+k,j+l,2);
                        P1 = Q(i,j)*Qy*P1;
                        if P1<0                 
                            P2 = Q(i,j)* (vec(i,j,1)*Ix(i+k,j+l) + vec(i,j,2)*Iy(i+k,j+l));
                            if P2<0
                                RidgeMap(i,j) = 0;
                            end
                        end
                    end
                    %%%%%%%%%%%%%%%%
                end
            end
        end
        
        
    end
end

RidgeMap = uint8(RidgeMap);
figure,imshow(RidgeMap);
%imwrite(RidgeMap, 'C:\Documents and Settings\cc\桌面\毕业设计\image_processed\ridge\matlab_ridge\RI1_guass80_SRridge.bmp');