

%Reading the image


im=imread('or13.tif');;

%size(im)        1136        1298


r=5; 
red1=[81:r:1055];
red2=[81:r:1217];


rad=40;

num=1;
mfd=[];

for ii=1:max(size(red1))
    ii
    for jj=1:max(size(red2))
  
      a=[red1(ii) red2(jj)];
     val1=a(1)-2*rad;
     val2=a(1)+2*rad;
     val3=a(2)-2*rad;
     val4=a(2)+2*rad;
     
     peq=im(val1:val2,val3:val4);
   vg=[0];
   %vg=[];
   at=sum(sum(peq));
   
   for k=2:50 
    
   
   vg=[vg 1 - sum(sum(imopen(peq,strel('square',k))))/at];
   

   end

   mfd(num,:)=vg;
   num=num+1;
    end
end
   
   plot(1:50,mfd') 

    
	%dlmwrite('mfdall_orto.txt',mfd,' ');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%After segmentation in R
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


eall=dlmread('eallorto4.txt'); %AA k=4

imll=im;

r=5; 
red1=[81:r:1055];
red2=[81:r:1217];


rad2=2;

num=1;
mfd=[];

for ii=1:max(size(red1))
    ii
    for jj=1:max(size(red2))
  
	if eall(num) == 1 
      a=[red1(ii) red2(jj)];
     val1=a(1)-2*rad2;
     val2=a(1)+2*rad2;
     val3=a(2)-2*rad2;
     val4=a(2)+2*rad2;
     
     imll(val1:val2,val3:val4)=0;
  end
   
	
	if eall(num) == 3
      a=[red1(ii) red2(jj)];
     val1=a(1)-2*rad2;
     val2=a(1)+2*rad2;
     val3=a(2)-2*rad2;
     val4=a(2)+2*rad2;
     
     imll(val1:val2,val3:val4)=255;
  end
	
	if eall(num) == 4 
      a=[red1(ii) red2(jj)];
     val1=a(1)-2*rad2;
     val2=a(1)+2*rad2;
     val3=a(2)-2*rad2;
     val4=a(2)+2*rad2;
     
     imll(val1:val2,val3:val4)=128;
  end
	
	num=num+1;
    end
end
   
	
	imshow(imll(81:1055,81:1217));
	


	Lc=imll(81:1055,81:1217);
	
	imshow(o(81:1055,81:1217));
	imcontrast
	
	
bw=Lc==128;
B2 = bwboundaries(bw);


hold on

for k = 1:length(B2)    
boundary = B2{k};   

 plot(boundary(:,2), boundary(:,1), 'w--', 'LineWidth', 2) 


end
	


bw=Lc==0;
B3 = bwboundaries(bw);


hold on

for k = 1:length(B3)    
boundary = B3{k};   

 plot(boundary(:,2), boundary(:,1), 'w', 'LineWidth', 2) 


end

	bw=Lc==255;
B = bwboundaries(bw);


hold on

for k = 1:length(B)    
boundary = B{k};   

 plot(boundary(:,2), boundary(:,1), 'k', 'LineWidth', 2) 


end
	
	

%%%%%%%%%%%%%%%%%%%%%%%%
%After segmentation in R
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


eall=dlmread('eallortok4.txt'); %k-means k=4

imll=im;

r=5; 
red1=[81:r:1055];
red2=[81:r:1217];


rad2=2;

num=1;
mfd=[];

for ii=1:max(size(red1))
    ii
    for jj=1:max(size(red2))
  
	if eall(num) == 2
      a=[red1(ii) red2(jj)];
     val1=a(1)-2*rad2;
     val2=a(1)+2*rad2;
     val3=a(2)-2*rad2;
     val4=a(2)+2*rad2;
     
     imll(val1:val2,val3:val4)=128;
  end
   
	
	if eall(num) == 3
      a=[red1(ii) red2(jj)];
     val1=a(1)-2*rad2;
     val2=a(1)+2*rad2;
     val3=a(2)-2*rad2;
     val4=a(2)+2*rad2;
     
     imll(val1:val2,val3:val4)=0;
  end
	
	if eall(num) == 4 
      a=[red1(ii) red2(jj)];
     val1=a(1)-2*rad2;
     val2=a(1)+2*rad2;
     val3=a(2)-2*rad2;
     val4=a(2)+2*rad2;
     
     imll(val1:val2,val3:val4)=255;
  end
	
	num=num+1;
    end
end
   
	imshow(imll(81:1055,81:1217));
	


	Lc=imll(81:1055,81:1217);
	
	imshow(o(81:1055,81:1217));
	imcontrast
	
	
bw=Lc==128;
B2 = bwboundaries(bw);


hold on

for k = 1:length(B2)    
boundary = B2{k};   

 plot(boundary(:,2), boundary(:,1), 'w--', 'LineWidth', 2) 


end
	


bw=Lc==0;
B3 = bwboundaries(bw);


hold on

for k = 1:length(B3)    
boundary = B3{k};   

 plot(boundary(:,2), boundary(:,1), 'w', 'LineWidth', 2) 


end

	bw=Lc==255;
B = bwboundaries(bw);


hold on

for k = 1:length(B)    
boundary = B{k};   

 plot(boundary(:,2), boundary(:,1), 'k', 'LineWidth', 2) 


end
	
	
