MAIN PROGRAM
close all;%close all the windows
clc;%clear command
I = imread('kgp.jpg');%to read an image
figure,imshow(I);title('original image');%to display an image
% Resize image for easier computation
B = imresize(I, [584 565]);
% Read image
im = im2double(B);
% Convert RGB to Gray via PCA
lab = rgb2lab(im);
f = 0;
wlab = reshape(bsxfun(@times,cat(3,1-f,f/2,f/2),lab),[],3);
[C,S] = pca(wlab);
S = reshape(S,size(lab));
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',128);
%% Background Exclusion
% Apply Average Filter
h = fspecial('average', [9 9]);
JF = imfilter(J, h);
figure(1), imshow(JF);title('Filtering after enhancement of retina image');
% Take the difference between the gray image and Average Filter
Z = imsubtract(JF, J);
figure, imshow(Z);
title('difference between the gray image and Average Filter');
%% Threshold using the IsoData Method
level=isodata(Z); % this is our threshold level
%% Convert to Binary
%BW = im2bw(Z, level-.008);
%RM Modified line 42 to line 44
BW = im2bw(Z, level-.0118/8);
%% Remove small pixels
%BW2 = bwareaopen(BW, 100);
BW2 = bwareaopen(BW, 200);
%% Overlay
BW2 = imcomplement(BW2);
out = imoverlay(B, BW2, [0 0 0]);
figure, imshow(out)
BB2= imcomplement(out);
title('overlay image');
figure, imshow(BB2)
title('out');
% % RM Modification
% %Extract Blood Vessels
% Threshold = 2.5;
% inImg = double(rgb2gray(out));
% bloodVessels = VesselExtract(inImg, Threshold);
% figure;
% subplot(121);imshow(inImg);title('Input Image');
% subplot(122);imshow(bloodVessels);title('Extracted Blood Vessels');
Threshold = 2.5;
inImg = double(rgb2gray(out));
bloodVessels = VesselExtract(inImg, Threshold);
bw=im2bw(bloodVessels);
q=bwboundaries(bw);
boundry=q{1};
figure,imshow(bw);
title('Input Image');
hold on;
plot(boundry(:,2),boundry(:,1),'g');
CROPPING PROGRAM
figure,imshow(bw);
hr=imfreehand(gca);
position = getPosition(hr);
bb=createMask(hr);
bw(bb == 0)=0;
figure,imshow(bw);
title('cropping');
hold off
ISOLATE PROGRAM
function level=isodata(I)
% ISODATA Compute global image threshold using iterative isodata method.
% LEVEL = ISODATA(I) computes a global threshold (LEVEL) that can be
% used to convert an intensity image to a binary image with IM2BW. LEVEL
% is a normalized intensity value that lies in the range [0, 1].
% This iterative technique for choosing a threshold was developed by Ridler
and Calvard .
% The histogram is initially segmented into two parts using a starting
threshold value such as 0 = 2B-1,
% half the maximum dynamic range.
% The sample mean (mf,0) of the gray values associated with the foreground
pixels and the sample mean (mb,0)
% of the gray values associated with the background pixels are computed. A
new threshold value 1 is now computed
% as the average of these two sample means. The process is repeated, based
upon the new threshold,
% until the threshold value does not change any more.
%
% Class Support
% -------------
% The input image I can be of class uint8, uint16, or double and it
% must be nonsparse. LEVEL is a double scalar.
%
% Example
% -------
% I = imread('blood1.tif');
% level = graythresh(I);
% BW = im2bw(I,level);
% imshow(BW)
%
% See also IM2BW.
%
% Reference :T.W. Ridler, S. Calvard, Picture thresholding using an iterative
selection method,
% IEEE Trans. System, Man and Cybernetics, SMC-8 (1978) 630-632.
% Convert all N-D arrays into a single column. Convert to uint8 for
% fastest histogram computation.
I = im2uint8(I(:));
% STEP 1: Compute mean intensity of image from histogram, set T=mean(I)
[counts,N]=imhist(I);
i=1;
mu=cumsum(counts);
T(i)=(sum(N.*counts))/mu(end);
T(i)=round(T(i));
% STEP 2: compute Mean above T (MAT) and Mean below T (MBT) using T from
% step 1
mu2=cumsum(counts(1:T(i)));
MBT=sum(N(1:T(i)).*counts(1:T(i)))/mu2(end);
mu3=cumsum(counts(T(i):end));
MAT=sum(N(T(i):end).*counts(T(i):end))/mu3(end);
i=i+1;
% new T = (MAT+MBT)/2
T(i)=round((MAT+MBT)/2);
% STEP 3 to n: repeat step 2 if T(i)~=T(i-1)
while abs(T(i)-T(i-1))>=1
mu2=cumsum(counts(1:T(i)));
MBT=sum(N(1:T(i)).*counts(1:T(i)))/mu2(end);
mu3=cumsum(counts(T(i):end));
MAT=sum(N(T(i):end).*counts(T(i):end))/mu3(end);
i=i+1;
T(i)=round((MAT+MBT)/2);
Threshold=T(i);
end
% Normalize the threshold to the range [i, 1].
level = (Threshold - 1) / (N(end) - 1);
BLOOD VESSELS EXTRACTION PROGRAM
%Program for Retinal Blood Vessel Extraction
%AUTHOR
%MATLAB,
%P.E.S 1College of Engineering
%Erode, KARNNATAKA, India.
%Program Description
%This program extracts blood vessels from a retina image using Kirsch's
Templates.
%Spatial Filtering of the input retina image is done with the Kirsch's
%Templates in different orientations. Followed by thresholding, results in
%the extracted blood vessels. The threshold can be varied to fine tune the
%output.
function bloodVessels = VesselExtract(inImg, threshold)
%Kirsch's Templates
h1=[5 -3 -3;
5 0 -3;
5 -3 -3]/15;
h2=[-3 -3 5;
-3 0 5;
-3 -3 5]/15;
h3=[-3 -3 -3;
5 0 -3;
5 5 -3]/15;
h4=[-3 5 5;
-3 0 5;
-3 -3 -3]/15;
h5=[-3 -3 -3;
-3 0 -3;
5 5 5]/15;
h6=[ 5 5 5;
-3 0 -3;
-3 -3 -3]/15;
h7=[-3 -3 -3;
-3 0 5;
-3 5 5]/15;
h8=[ 5 5 -3;
5 0 -3;
-3 -3 -3]/15;
%Spatial Filtering by Kirsch's Templates
t1=filter2(h1,inImg);
t2=filter2(h2,inImg);
t3=filter2(h3,inImg);
t4=filter2(h4,inImg);
t5=filter2(h5,inImg);
t6=filter2(h6,inImg);
t7=filter2(h7,inImg);
t8=filter2(h8,inImg);
s=size(inImg);
bloodVessels=zeros(s(1),s(2));
temp=zeros(1,8);
for i=1:s(1)
for j=1:s(2)
temp(1)=t1(i,j);temp(2)=t2(i,j);temp(3)=t3(i,j);temp(4)=t4(i,j);
temp(5)=t5(i,j);temp(6)=t6(i,j);temp(7)=t7(i,j);temp(8)=t8(i,j);
if(max(temp)>threshold)
bloodVessels(i,j)=max(temp);
end
end
end
VESSEL WIDTH CALCULATION
w=csvread('E:\project\cr\p10.dat');
wt=imrotate(w,0)
[x,y]=find(w);
cnt=1;temp_cnt=1;vxx=0;
for i=1:1:size(x,1)-1
if(y(i+1,1)==y(i,1))
vxx(cnt)=x(i,1);
cnt=cnt+1;
else
wx(temp_cnt)=x(i,1);
wxx(temp_cnt)=x(i,1);
if(cnt~=1)
wx(temp_cnt)=min(vxx);
wxx(temp_cnt)=max(vxx);
end
wy(temp_cnt)=y(i);
temp_cnt=temp_cnt+1;
cnt=1;
clear vxx;
end
end
%figure(1);imshow(wt);hold on;plot(wy,wx,'.r');hold on;plot(wy,wxx,'.m');
xs=size(wx);
seg_len=xs(2);
sub_seg_size=round(seg_len/4);
wpx1=wx(1,1);
wpxx1=wxx(1,1);
wpx2=wx(1,1+sub_seg_size);
wpxx2=wxx(1,1+sub_seg_size);
wpx3=wx(1,1+2*sub_seg_size);
wpxx3=wxx(1,1+2*sub_seg_size);
wpx4=wx(1,1+3*sub_seg_size);
wpxx4=wxx(1,1+3*sub_seg_size);
w1=abs(wpx1-wpxx1);
w2=abs(wpx2-wpxx2);
w3=abs(wpx3-wpxx3);
w4=abs(wpx4-wpxx4);
wii=1:1:seg_len;
for wi=1:1:seg_len
width(wi)=abs(wx(1,wi)-wxx(1,wi));
end
avg_width=round(sum(width)/seg_len);
figure(2);plot(wii,width);
DIFFERENT ENHANCEMENT TECHNIQUE
clc;
clear all;
im=imread('kgp.jpg');%to read an image
subplot(4,3,1);
imshow(im);%to display an image
title('original image');
t=rgb2gray(im);%to convert from color image to gray image or extraction of
gray image
subplot(4,3,2);
imshow(t);
title('gray image');
a=histeq(t);%enhance constrast using histogram equalization enhancement
technique
subplot(4,3,3);
imshow(a);
title('histogram equalization');
b=adapthisteq(t);%constrast-limited adaptation equilization enhancement
technique avoid the amplifiying noise
subplot(4,3,4);
imshow(b);
title('adaptive equalization');
h=fspecial('gaussian',15,7);%create predefined 2-D filter
c=imfilter(t,h);
subplot(4,3,5);
imshow(c);
title('imfilter');
d=imadjust(t);%adjust the image intensity values or colourmap
subplot(4,3,6);
imshow(d);
title('imadjust technique');
b=adapthisteq(c);%constrast-limited adaptive histogram equalization
subplot(4,3,7);
imshow(b);
title('adptation eqlzation apply to fltrd image');
b=adapthisteq(d);%constrast-limited adaptive histogram equalization
subplot(4,3,8);
imshow(b);
title('adption eqlzation apply to imadjust image');
g=medfilt2(t);%reduce the salt and pepper noisee and 2-D median filtering
subplot(4,3,9);
imshow(g);
title('median filter');
gms=medfilt2(t,'symmetric');%to reduce the black point sorrounding the image
subplot(4,3,10);
imshow(gms);
title('median filter with symmetric');
