结合血管信息的视盘定位

结合血管信息的视盘定位

%% matched filter2，测试程序
clc,close all,clear all;
%% 读取图像
img=imread('E:\糖网\预处理\shipantiqu\shipan4\ROI4.jpg');
% img=imread('E:\HEI_MED\images\HEI_MED40.jpg');

%% 转换成灰度图像
if length(size(img))==3
    img0=rgb2gray(img);
end

%% 自适应直方图均衡化
img1 = adapthisteq(img0, 'NumTiles', [15 15], 'ClipLimit', 0.05);
figure,imshow(img1);

%% 高斯匹配滤波
[g,bg]=matchedFilter2(img1);
figure,imshow(bg);

%% 腐蚀处理
he = strel('disk',2);
img2=imerode(bg,he);
figure;imshow(img2);

%% 提取连通域
img3=xueguanLianTongYu(img2);
figure,imshow(img3);

%% 形态学膨胀得到较为完整的包含渗出物的候选区域
he = strel('disk',1);
img3=imdilate(img3,he);
figure,imshow(img3);

img0=double(img0);
img4=img0.*img3;

% img0=double(img0);
% img4=img0.*(1-img3);
% figure,imshow(uint8(img4));
% 
% img5=imreconstruct(img4,img0);
% figure,imshow(uint8(img5));
% 
% img6=img0-img5;
% figure,imshow(img6);
%% 视盘定位
N1=[135  120  120  105  90  75  60  60  45;
    150  135  120  105  90  75  60  45  30;
    150  135  135  120  90  60  45  30  30;
    165  165  150  135  90  45  30  15  15;
    0  0  0  0  0  0  0  0  0;
    15  15  30  45  90  135  150  165  165;
    30  30  45  60  90  120  135  150  150;
    30  45  60  75  90  105  120  135  150;
    45  60  60  75  90  105  120  120  135];
N2=[2 2 2 2 2 2 2 2 2;
    2 2 2 2 2 2 2 2 2;
    2 2 3 3 3 3 3 2 2;
    2 2 3 3 3 3 3 2 2;
    2 2 3 3 3 3 3 2 2;
    2 2 3 3 3 3 3 2 2;
    2 2 3 3 3 3 3 2 2;
    2 2 2 2 2 2 2 2 2;
    2 2 2 2 2 2 2 2 2];
N3=N1.*N2;
y=imfilter(img4,N3,'replicate');
figure,
subplot(121),imshow(y);
[n]=max(max(y));
[p q]=find(n==y);
m=y>=n;
% imwrite(m,'E:\糖网\预处理\shipantiqu\shipan6\xiangying6.jpg')
subplot(1,2,2),imshow(m);
hold on;
plot(q, p, 'g*'); 

需要用到的一些function

function [g,bg]=matchedFilter2(f)
 
f=double(f);
% subplot(1,2,1);
% imshow(uint8(f));
% title('origin image');
% mean filter
f=medfilt2(f,[5,5]);
% f=medfilt2(f,[21 1]);
% f=medfilt2(f,[1,7]);
% 参数
os=12;  % 角度的个数
sigma=2;
tim=3;
L=9;
t=120; % 全局阈值,需要多次尝试
 
thetas=0:(os-1);
thetas=thetas.*(180/os);
N1=-tim*sigma:tim*sigma;
N1=-exp(-(N1.^2)/(2*sigma*sigma));
N=repmat(N1,[2*floor(L/2)+1,1]);
r2=floor(L/2);
c2=floor(tim*sigma);
[m,n]=size(f);
RNs=cell(1,os);  % rotated kernals 旋转角尺
MFRs=cell(1,os); % filtered images 
g1=f;
 
% matched filter
for i=1:os
    theta=thetas(i);
    RN=imrotate(N,theta);
    %去掉多余的0行和零列
    RN=RN(:,any(RN));
    RN=RN(any(RN'),:);
    meanN=mean2(RN);
    RN=RN-meanN;
    RNs{1,i}=RN;
    MFRs{1,i}=imfilter(f,RN,'conv','symmetric');
end
 
% get the max response
g=MFRs{1,1};
for j=2:os
    g=max(g,MFRs{1,j});
end
bg=g>t;
 
% subplot(1,2,2);
% imshow(bg);
% title('filterd image');
end

function [img]=xueguanLianTongYu(I)
if length(size(I))>2
    I = rgb2gray(I);
end
if ~islogical(I)
    imBw = im2bw(I);                        %转换为二值化图像
else
    imBw = I;
end
imBw = im2bw(I);                        %转换为二值化图像
imLabel = bwlabel(imBw);                %对各连通域进行标记
stats = regionprops(imLabel,'Area');    %求各连通域的大小
area = cat(1,stats.Area);
index = find((area >=1000)&(area<=10000000));        %求最大连通域的索引
img = ismember(imLabel,index);          %获取最大连通域图像

原图实际效果