MATLAB Answers

Automated detection of diabetic retinopathy

조회 수: 7(최근 30일)
Anas Bilal
Anas Bilal 2020년 4월 23일
댓글: Mrutyunjaya Hiremath 2020년 4월 28일
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
Anas Bilal
Anas Bilal 2020년 4월 25일
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

댓글을 달려면 로그인하십시오.

채택된 답변

Mrutyunjaya Hiremath
Mrutyunjaya Hiremath 2020년 4월 26일
Hello Anas Bilal,
Code tested working on 2012.
  댓글 수: 3
Mrutyunjaya Hiremath
Mrutyunjaya Hiremath 2020년 4월 28일
@Anas,
Which version of Matlab are you using?

댓글을 달려면 로그인하십시오.

Community Treasure Hunt

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

Start Hunting!

Translated by