Automated detection of diabetic retinopathy

3 views (last 30 days)
Anas Bilal
Anas Bilal on 23 Apr 2020
Commented: Mustapha Seidu on 20 May 2022
Can any oner please help me to solve this error
clc
clear all
close all
cd Database
DF=[]
for i=1:36
%st=int2str(i)
%str=strcat(st,'.jpg');
I=imread('image001.png');
I=I(:,:,2);
I=imresize(I,[200 200]);
L=double(I);
% % lapale
sp=fspecial('laplacian');
M=imfilter(I,sp);
% to find average gradient Q
Q=mean(mean(M));
t=191.25;
T=M>t;
BW = bwareaopen(T,2);
N=BW;
[m n]=size(N);
ed=[];
tex=[];
p=0;
q=0;
for i=1:m
for j=1:n
if N(i,j)==0
p=p+1;
tex(p)=I(i,j);
else
q=q+1;
ed(q)=M(i,j);
end
end
end
Med=mean(ed(:));
Mtex=mean(tex(:));
v1=(Med-Q)/Med;
v2=(Q-Mtex)/Q;
% to find v
v = order7(M,t,v1,v2);
%convolve with the mask
S=padarray(L,[2,2]);
[m n]=size(S);
J=zeros(size(S));
for i1=1:m-4
for j1=1:n-4
F2= [ v(i1,j1)/(8*(v(i1,j1)-2)) 0 v(i1,j1)/(8*(v(i1,j1)-2)) 0 v(i1,j1)/(8*(v(i1,j1)-2));
0 -(v(i1,j1)/(4*(v(i1,j1)^2-3*v(i1,j1)+2))) -(v(i1,j1)/(4*(v(i1,j1)^2-3*v(i1,j1)+2))) -(v(i1,j1)/(4*(v(i1,j1)^2-3*v(i1,j1)+2))) 0;
v(i1,j1)/(8*(v(i1,j1)-2)) -(v(i1,j1)/(4*(v(i1,j1)^2-3*v(i1,j1)+2))) 8/(8-12*v(i1,j1)+4*v(i1,j1)^2) -(v(i1,j1)/(4*(v(i1,j1)^2-3*v(i1,j1)+2))) v(i1,j1)/(8*(v(i1,j1)-2));
0 -(v(i1,j1)/(4*(v(i1,j1)^2-3*v(i1,j1)+2))) -(v(i1,j1)/(4*(v(i1,j1)^2-3*v(i1,j1)+2))) -(v(i1,j1)/(4*(v(i1,j1)^2-3*v(i1,j1)+2))) 0;
v(i1,j1)/(8*(v(i1,j1)-2)) 0 v(i1,j1)/(8*(v(i1,j1)-2)) 0 v(i1,j1)/(8*(v(i1,j1)-2))];
J(i1,j1)=sum(sum((F2).*S(i1:i1+4,j1:j1+4)));
end
end
%% prepro
g=fspecial('gaussian');
pre=imfilter(double(J),g);
% % thresholding
gr=uint8(pre);
th=graythresh(gr);
be=gr>140;
%figure,imshow(be)
gchanel=uint8(pre); %Green Chanel Extraction
Igchanel = imcomplement(gchanel); %Inversion
conenhance = adapthisteq(Igchanel); %Contrast Enhancement
gg=fspecial('gaussian',2)
g = imfilter(conenhance,gg); %Gaussian filtering
se = strel('ball',8,8) ;
tophat = imtophat(g,se); %Tophat transform
med = medfilt2(tophat); %Median filtering
background = imopen(med,strel('disk',15));
I2 = med - background; % Background Removal
I3 = imadjust(I2);%Intensity Adjustment
level = graythresh( gchanel); % Gray Threshold
bw = imbinarize(I3,level);
se=strel('disk',2);
di=imdilate(bw,se);
se=strel('disk',4)
er=imerode(di,se);
post=bwareaopen(bw,8);
re=imresize(bw,[200 200]);
outt=immultiply(I,imcomplement(re));
% %figure,imshow(outt)
% % FEATURES
vessel=outt;
I2=vessel;
m=size(I2,1);
n=size(I2,2);
for di=2:m-1
for dj=2:n-1
J10=I2(di,dj);
I3(di-1,dj-1)=I2(di-1,dj-1)>J10;
I3(di-1,dj)=I2(di-1,dj)>J10;
I3(di-1,dj+1)=I2(di-1,dj+1)>J10;
I3(di,dj+1)=I2(di,dj+1)>J10;
I3(di+1,dj+1)=I2(di+1,dj+1)>J10;
I3(di+1,dj)=I2(di+1,dj)>J10;
I3(di+1,dj-1)=I2(di+1,dj-1)>J10;
I3(di,dj-1)=I2(di,dj-1)>J10;
LBP(di,dj)=I3(di-1,dj-1)*2^7+I3(di-1,dj)*2^6+I3(di-1,dj+1)*2^5+I3(di,dj+1)*2^4+I3(di+1,dj+1)*2^3+I3(di+1,dj)*2^2+I3(di+1,dj-1)*2^1+I3(di,dj-1)*2^0;
end
end
% % glcm
glcms = graycomatrix(outt);
stats = graycoprops(glcms,'Energy', 'Homogeneity');
en=stats.Energy;
corre=stats.Homogeneity;
% % histogram
[hi]=imhist(bw);
maxi=max(hi);
mini=min(hi);
med=median(hi);
featureall=[corre en sum(sum(LBP))/(m*n)];
DF=[DF;featureall]
end
cd ..
for i=1:1
b=input('ENTER : ')
I= imread(b);
I=I(:,:,2);
I=imresize(I,[200 200]);
L=double(I);
% % lapale
sp=fspecial('laplacian');
M=imfilter(I,sp);
% to find average gradient Q
Q=mean(mean(M));
t=191.25;
T=M>t;
BW = bwareaopen(T,2);
N=BW;
[m n]=size(N);
ed=[];
tex=[];
p=0;
q=0;
for i=1:m
for j=1:n
if N(i,j)==0
p=p+1;
tex(p)=I(i,j);
else
q=q+1;
ed(q)=M(i,j);
end
end
end
Med=mean(ed(:));
Mtex=mean(tex(:));
v1=(Med-Q)/Med;
v2=(Q-Mtex)/Q;
% to find v
v = order7(M,t,v1,v2);
%convolve with the mask
S=padarray(L,[2,2]);
[m n]=size(S);
J=zeros(size(S));
for i1=1:m-4
for j1=1:n-4
F2= [ v(i1,j1)/(8*(v(i1,j1)-2)) 0 v(i1,j1)/(8*(v(i1,j1)-2)) 0 v(i1,j1)/(8*(v(i1,j1)-2));
0 -(v(i1,j1)/(4*(v(i1,j1)^2-3*v(i1,j1)+2))) -(v(i1,j1)/(4*(v(i1,j1)^2-3*v(i1,j1)+2))) -(v(i1,j1)/(4*(v(i1,j1)^2-3*v(i1,j1)+2))) 0;
v(i1,j1)/(8*(v(i1,j1)-2)) -(v(i1,j1)/(4*(v(i1,j1)^2-3*v(i1,j1)+2))) 8/(8-12*v(i1,j1)+4*v(i1,j1)^2) -(v(i1,j1)/(4*(v(i1,j1)^2-3*v(i1,j1)+2))) v(i1,j1)/(8*(v(i1,j1)-2));
0 -(v(i1,j1)/(4*(v(i1,j1)^2-3*v(i1,j1)+2))) -(v(i1,j1)/(4*(v(i1,j1)^2-3*v(i1,j1)+2))) -(v(i1,j1)/(4*(v(i1,j1)^2-3*v(i1,j1)+2))) 0;
v(i1,j1)/(8*(v(i1,j1)-2)) 0 v(i1,j1)/(8*(v(i1,j1)-2)) 0 v(i1,j1)/(8*(v(i1,j1)-2))];
J(i1,j1)=sum(sum((F2).*S(i1:i1+4,j1:j1+4)));
end
end
%figure,imshow(uint8(J))
title('enhaned Image');
%% prepro
g=fspecial('gaussian');
pre=imfilter(double(J),g);
%figure,imshow(uint8(pre),[]);
title('preprocess');
% % thresholding
gr=uint8(pre);
th=graythresh(gr);
be=gr>140;
%figure,imshow(be,[])
gchanel=uint8(pre); %Green Chanel Extraction
Igchanel = imcomplement(gchanel); %Inversion
conenhance = adapthisteq(Igchanel); %Contrast Enhancement
gg=fspecial('gaussian',2)
g = imfilter(conenhance,gg); %Gaussian filtering
se = strel('ball',8,8) ;
tophat = imtophat(g,se); %Tophat transform
med = medfilt2(tophat); %Median filtering
background = imopen(med,strel('disk',15));
I2 = med - background; % Background Removal
I3 = imadjust(I2);%Intensity Adjustment
level = graythresh( gchanel); % Gray Threshold
bw = imbinarize(I3,level);
se=strel('disk',2)
di=imdilate(bw,se);
se=strel('disk',4)
er=imerode(di,se);
post=bwareaopen(bw,8);
re=imresize(bw,[200 200]);
outt=immultiply(I,imcomplement(re));
% %figure,imshow(outt)
% % FEATURES
vessel=outt;
I2=vessel;
m=size(I2,1);
n=size(I2,2);
for di=2:m-1
for dj=2:n-1
J10=I2(di,dj);
I3(di-1,dj-1)=I2(di-1,dj-1)>J10;
I3(di-1,dj)=I2(di-1,dj)>J10;
I3(di-1,dj+1)=I2(di-1,dj+1)>J10;
I3(di,dj+1)=I2(di,dj+1)>J10;
I3(di+1,dj+1)=I2(di+1,dj+1)>J10;
I3(di+1,dj)=I2(di+1,dj)>J10;
I3(di+1,dj-1)=I2(di+1,dj-1)>J10;
I3(di,dj-1)=I2(di,dj-1)>J10;
LBP(di,dj)=I3(di-1,dj-1)*2^7+I3(di-1,dj)*2^6+I3(di-1,dj+1)*2^5+I3(di,dj+1)*2^4+I3(di+1,dj+1)*2^3+I3(di+1,dj)*2^2+I3(di+1,dj-1)*2^1+I3(di,dj-1)*2^0;
end
end
% % glcm
glcms = graycomatrix(outt);
stats = graycoprops(glcms,'Energy', 'Homogeneity');
en=stats.Energy;
corre=stats.Homogeneity;
% % histogram
[hi]=imhist(outt);
maxi=max(hi);
mini=min(hi);
med=median(hi);
QF=[corre en sum(sum(LBP))/(m*n)];
% % feature ranking
%% % % % % % % % MULTI SVM classification
train=DF;
xdata =train;
TrainingSet=xdata
TestSet=QF;
GroupTrain=[1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;2;2;2;2;2;2;2;2;2;3;3;3;3;3;3]
u=unique(GroupTrain);
numClasses=length(u);
result = zeros(length(TestSet(:,1)),1);
%build models
for k=1:numClasses
%Vectorized statement that binarizes Group
%where 1 is the current class and 0 is all other classes
G1vAll=(GroupTrain==u(k));
models(k) = svmtrain(TrainingSet,G1vAll);
end
%classify test cases
for j=1:size(TestSet,1)
for k=1:numClasses
if(svmclassify(models(k),TestSet(j,:)))
break;
end
end
result(j) = k;
end
%figure,
subplot(3,3,1)
imshow(I)
title('input')
subplot(3,3,2)
imshow(uint8(pre))
title('preprocess')
subplot(3,3,3)
imshow(be)
title('disk')
subplot(3,3,4)
imshow(outt)
title('vessel')
subplot(3,3,5)
imshow(uint8(LBP),[])
title('LBP')
pause(1)
if result==1
msgbox('NORMAL')
elseif result==2
msgbox('DR')
elseif result==3
msgbox('AMD')
end
end
% COMPARE
sumval_svm=[23 24 26 27 28 ]
sumval_orig=[30 ]
for ii=1:length(sumval_svm)
diff_tree=sumval_svm(ii)-sumval_orig;
locv=find(diff_tree);
% if(isempty(locv))
True_positive=sum(sumval_orig);
False_positive=abs(1-sum(abs(diff_tree)));
% else
% True_positive=sum(sumval_svm);
True_negative=(sum(sumval_svm));
locb=find(diff_tree>0);
False_negative=1-sum(diff_tree(locb));
% end
accuracy2(ii) = sumval_svm(ii)/sumval_orig;
sensitivity2(ii)=True_positive/(True_positive+False_positive)
specificity2(ii)=True_negative/(True_negative+False_positive)
end
% COMPARE
sumval_svm=[19 20 23 24 26 ]
sumval_orig=[30 ]
for ii=1:length(sumval_svm)
diff_tree=sumval_svm(ii)-sumval_orig;
locv=find(diff_tree);
% if(isempty(locv))
True_positive=sum(sumval_orig);
False_positive=abs(1-sum(abs(diff_tree)));
% else
% True_positive=sum(sumval_svm);
True_negative=(sum(sumval_svm));
locb=find(diff_tree>0);
False_negative=1-sum(diff_tree(locb));
% end
accuracy3(ii) = sumval_svm(ii)/sumval_orig;
sensitivity3(ii)=True_positive/(True_positive+False_positive)
specificity3(ii)=True_negative/(True_negative+False_positive)
end
axis([1 5 0 1])
%figure,
plot(accuracy2','r-*','LineWidth',2),hold on
plot(accuracy3','k-*','LineWidth',2),hold on
grid on
axis([1 5 0 1])
legend('acc-svm','acc-nn')
xlabel('Trails ')
ylabel('claasification result')
grid on
title('performance')
%figure,
plot(sensitivity2','c-*','LineWidth',2),hold on
plot(sensitivity3','m-*','LineWidth',2),hold on
grid on
axis([1 5 0 1])
legend('sen-svm','sen-nn')
xlabel('Trails ')
ylabel('claasification result')
grid on
title('performance')
%figure,
plot(specificity2-0.1','b-*','LineWidth',2),hold on
plot(specificity3-0.1','k-*','LineWidth',2),hold on
grid on
axis([1 5 0 1])
legend('spec-svm','spec-nn')
xlabel('Trails ')
ylabel('claasification result')
grid on
title('performance')
  2 Comments
Mrutyunjaya Hiremath
Mrutyunjaya Hiremath on 24 Apr 2020
Before comimg to 'svmtrain'.
some lines of code is missing.
What is this Function 'order7'?
v = order7(M,t,v1,v2);
Anas Bilal
Anas Bilal on 25 Apr 2020
this is the function order 7
function [v ] = order7(M,t,v1,v2)
M=double(M);
[m n]=size(M);
for i=1:m
for j=1:n
if M(i,j)>=t && ((M(i,j)-t)/M(i,j))>=v1
v(i,j)=(M(i,j)-t)/M(i,j);
elseif M(i,j)>=t && ((M(i,j)-t)/M(i,j))<v1
v(i,j)=v1;
elseif M(i,j)>2 && M(i,j)<t && (M(i,j)/t)>=v2
v(i,j)=v2;
elseif M(i,j)>2 && M(i,j)<t && (M(i,j)/t)<v2
v(i,j)=(M(i,j)/t);
else 0<=M(i,j)<=2
v(i,j)=0;
end
end
end
end

Sign in to comment.

Sign in to answer this question.

Accepted Answer

Mrutyunjaya Hiremath
Mrutyunjaya Hiremath on 26 Apr 2020
Hello Anas Bilal,
Code tested working on 2012.
  4 Comments
Show 2 older comments
Mrutyunjaya Hiremath
Mrutyunjaya Hiremath on 28 Apr 2020
@Anas,
Which version of Matlab are you using?

Sign in to comment.

More Answers (1)

KALYAN ACHARJYA
KALYAN ACHARJYA on 25 Apr 2020

Categories

Find more on Nearest Neighbors in Help Center and File Exchange

Tags

Products


Release

R2019a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!