眼底图像血管增强与分割--（3）基于PCA的增强算法原理及实现

偶然间在网上发现的一种眼底血管分割算法，没有论文出版，但是效果不错，我这里整理一下作者的思路，并测试。

算法步骤主要分为：

（1）输入一幅彩色图像，首先将其转换到Lab空间；

（2）将a,b通道分量置0，对L分量向量化然后PCA分析；

（3）对得到的SCORE矩阵L分量（第一通道）归一化，然后进行CLAHE直方图均衡；

（4）图像背景处理。具体为：对（3）的均衡结果中值平滑，然后与平滑前图像做差，得到细节图像Z；然后对Z阈值处理得到二值化图像BW，利用bwareaopen函数（类似于形态学操作）将BW中面积小的区域去除，得到BW2;

（5）对BW2取反色，并将整体背景设为黑色，就能看到血管边缘了。

具体代码如下：

% Retinal Blood Vessel Segmentation Test
clc; clear all; close all;
I = imread('.\Test Images\17_right.jpeg');
figure(1), imshow(I),title('original image');
% Resize image for easier computation
[len,wid,channel]=size(I);
B = imresize(I, [len/4 wid/4]); %简化计算量
% B = imresize(I, [600 600]);

% Convert RGB to Gray via PCA
lab = rgb2lab(im2double(B));
%get wlab's method1
% f = 0;
% A_F=cat(3,1-f,f/2,f/2);
% B_F=bsxfun(@times,A_F,lab);
% wlab = reshape(B_F,[],3);

%get wlab's method2,same to method1
lab(:,:,2)=0;lab(:,:,3)=0;
wlab = reshape(lab,[],3); %向量化

[C,S] = pca(wlab); %主成分分析
S = reshape(S,size(lab));%S为PCA后新坐标下的矩阵
S = S(:,:,1);
gray = (S-min(S(:)))./(max(S(:))-min(S(:))); %归一化

%% Contrast Enhancment of gray image using CLAHE
J = adapthisteq(gray,'numTiles',[8 8],'nBins',256); %CLAHE直方图均衡

%% Background Exclusion
% Apply Average Filter
h = fspecial('average', [11 11]);
JF = imfilter(J, h);
% Take the difference between the gray image and Average Filter
Z = imsubtract(JF, J); %相减得到高频细节信息
figure(2),
subplot(211),imshow(JF),title('CLAHE后均值滤波');
subplot(212), imshow(Z),title('CLAHE后均值滤波差值');
%% Threshold using the IsoData Method
%level=isodata(Z) ; % threshold level
level = graythresh(Z);
%% Convert to Binary
BW = im2bw(Z, level-0.008);
%% Remove small pixels
BW2 = bwareaopen(BW, 50); %删除二值图像BW中面积小于50的对象，默认情况下使用8邻域，参数可调
figure(3),
subplot(211),imshow(BW),title('BW');
subplot(212),imshow(BW2),title('BW2_1');
%% Overlay
BW2 = imcomplement(BW2); %对图像数据进行取反运算（可以实现底片效果）
out = imoverlay(B, BW2, [0 0 0]); %底片整体颜色设置为黑色
figure(4),
subplot(211),imshow(out),title('最终结果');
subplot(212),imshow(BW2),title('BW2_2');

实验结果：









整个测试项目下载地址：http://download.csdn.net/download/piaoxuezhong/10046407

附录：

matlab命令解释：bsxfun,repmat,reshape,
http://blog.sina.com.cn/s/blog_6ca002a50100wvu1.html
matlab滚动条设置：
http://blog.csdn.net/fdybit/article/details/7976542
matlab pca讲解：
http://blog.csdn.net/watkinsong/article/details/38536463 http://blog.csdn.net/watkinsong/article/details/8234766
coyefilter：http://cn.mathworks.com/matlabcentral/fileexchange/50839-a-novel-retinal-blood-vessel-segmentation-algorithm-for-fundus-images