MATLAB: Feature Extraction from Breast ROI

Hello,

I am trying to extract features from breast images using the code below:

clc;clear;close all
%% Getting Image
i=imread('B16.jpg');
figure(1)
imshow(i);title('Original Photo')
% if image is rgb
try
    i=rgb2gray(i);
end
%% Crop The Breast
z=im2bw(i,0.1);
figure(2)
imshow(z);title('Original B&W')
info=regionprops(z);
a=cat(1,info.Area);
[m,l]=max(a);
X=info(l).Centroid;
bw2=bwselect(z,X(1),X(2),8);
i=immultiply(i,bw2);
figure(3)
imshow(i);
title('Getting the Breast and Muscle')
%% Deleting Black Ground
% We will delete the black corners
% So that we can select the muscle
% using bwselect
% convert to B&W first time
[x,y]=size(z);
tst1=zeros(x,y);
% detect empty rows
r1=[];
m=1;
for j=1:x
    if z(j,:)==tst1(j,:)
        r1(m)=j;
        m=m+1;
    end
end
% detect empty columns
r2=[];
m=1;
for j=1:y
    if z(:,j)==tst1(:,j)
        r2(m)=j;
        m=m+1;
    end
end
% Deleting
i(:,r2)=[];
i(r1,:)=[];
figure(4)
imshow(i);title('after deleting background');
%% Deleting the Muscle
if i(1,1)~=0
    c=3;
    r=3;
else
    r=3;
    c=size(i,2)-3;
end
z2=im2bw(i,0.5);
bw3=bwselect(z2,c,r,8);
bw3=~bw3;
ratio=min(sum(bw3)/sum(z2));
if ratio>=1
    i=immultiply(i,bw3);
else
    z2=im2bw(i,0.75);
    bw3=bwselect(z2,c,r,8);
    ratio2=min(sum(bw3)/sum(z2));
    if round(ratio2)==0
        lvl=graythresh(i);
        z2=im2bw(i,1.75*lvl);
        bw3=bwselect(z2,c,r,8);
        bw3=~bw3;
        i=immultiply(i,bw3);
    else
        bw3=~bw3;
        i=immultiply(i,bw3);
    end
end
figure(5)
imshow(i)
title('Getting only the Breast')
%% Adaptive Median Filter
% clc;
% clear;
% close all;
% a=imread('M1.jpg');
% b=rgb2gray(i);
J = imnoise(i,'salt & pepper', 0.02);
NoisyImage=J;
[R C P]=size(NoisyImage);
OutImage=zeros(R,C);
figure;
% imshow(J);
Zmin=[];
Zmax=[];
Zmed=[];
for i=1:R
    for j=1:C
               if (i==1 & j==1)
                   
                   
                        
                        
          % for right top corner[8,7,6]
                elseif (i==1 & j==C)
                        
                        
			% for bottom left corner[2,3,4]
                elseif (i==R & j==1)
                        
                        
        		     % for bottom right corner[8,1,2]
                elseif (i==R & j==C)
                        
                        
        		    %for top edge[8,7,6,5,4]
                elseif (i==1)
		      			 	
		      			 
                     		      			% for right edge[2,1,8,7,6]
                elseif (i==R)
		      			 
		      			% // for bottom edge[8,1,2,3,4]
		      			 
                elseif (j==C)
		      			 
						 
                         %// for left edge[2,3,4,5,6]
						 
                elseif (j==1)
                    
						 
                            
                         
						 
               else
                    
                   
						 	
                             SR1 = NoisyImage((i-1),(j-1));
 							 SR2 = NoisyImage((i-1),(j));
    						 SR3 = NoisyImage((i-1),(j+1));
 							 SR4 = NoisyImage((i),(j-1));
                             SR5 = NoisyImage(i,j);
 							 SR6 = NoisyImage((i),(j+1));
 							 SR7 = NoisyImage((i+1),(j-1));
 							 SR8 = NoisyImage((i+1),(j));
 							 SR9 = NoisyImage((i+1)),((j+1));
                             TempPixel=[SR1,SR2,SR3,SR4,SR5,SR6,SR7,SR8,SR9];
                             Zxy=NoisyImage(i,j);
                             Zmin=min(TempPixel);
                             Zmax=max(TempPixel);
                             Zmed=median(TempPixel);
                             A1 = Zmed - Zmin;
                             A2 = Zmed - Zmax;
                         
                             if A1 > 0 && A2 < 0
                             
                                  %   go to level B
                                  B1 = Zxy - Zmin;
                                  B2 = Zxy - Zmax;
                                  if B1 > 0 && B2 < 0
                                      OutImage(i,j)= Zxy;
                                  else
                                      OutImage(i,j)= Zmed;
                                  end
                             else
                                 
                                 if ((R > 4 && R < R-5) && (C > 4 && C < C-5))
                                 
                                 S1 = NoisyImage((i-1),(j-1));
                                 S2 = NoisyImage((i-2),(j-2));
                                 S3 = NoisyImage((i-1),(j));
                                 S4 = NoisyImage((i-2),(j));
                                 S5 = NoisyImage((i-1),(j+1));
                                 S6 = NoisyImage((i-2),(j+2));
                                 S7 = NoisyImage((i),(j-1));
                                 S8 = NoisyImage((i),(j-2));
                                 S9 = NoisyImage(i,j);
                                 S10 = NoisyImage((i),(j+1));
                                 S11 = NoisyImage((i),(j+2));
                                 S12 = NoisyImage((i+1),(j-1));
                                 S13 = NoisyImage((i+2),(j-2));
                                 S14 = NoisyImage((i+1),(j));
                                 S15 = NoisyImage((i+2),(j));
                                 S16 = NoisyImage((i+1)),((j+1));
                                 S17 = NoisyImage((i+2)),((j+2));
                                 TempPixel2=[S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,S12,S13,S14,S15,S16,S17];
                                 Zmed2=median(TempPixel2);
                                 OutImage(i,j)= Zmed2;
                                 else
                                     
                                 OutImage(i,j)= Zmed;
                                 end
                             end         
                     
          end    
      
    end
end
imshow(OutImage,[]);
title('Adaptive Median Filter')
disp('exit');
%%GLCM Feature Extraction
% Y=rgb2gray(OutImage);
Y=double(OutImage);
% Y=OutImage;
% ShapeTexture=wlt4(Y);
% statsa = GLCM_Features4(Y,ShapeTexture);
% ExtractedFeatures1=statsa;
% glcm2=graycomatrix(Y,'Offset',[2 0;0 2]);
% glcm2=graycomatrix(Y);
% stats = GLCM_Features1(glcm2,0);
% ExtractedFeatures1=stats;
% statsTable = struct2table(stats);
% statsArray = table2array(statsTable);
% statsArray'
        k=2; % k: number of regions
        g=2; % g: number of GMM components
        beta=1; % beta: unitary vs. pairwise
        EM_iter=10; % max num of iterations

        MAP_iter=10; % max num of iterations
            [X,GMM,ShapeTexture]=image_kmeans(Y,k,g);
        [X,Y,GMM]=HMRF_EM(X,Y,GMM,k,g,EM_iter,MAP_iter,beta);
        Y=Y*80;
        Y=uint8(Y);
        Y =double(Y);
        statsa = glcm(Y,0,ShapeTexture);
        ExtractedFeatures1=statsa;
         imshow(Y,[]);

However, I am getting the following error:

Error using reshape

To RESHAPE the number of elements must not change.

Error in MRF_MAP (line 10)

y=reshape(Y,[m*n 3]);

Error in HMRF_EM (line 8)

[X sum_U(it)]=MRF_MAP(X,Y,GMM,k,g,MAP_iter,beta,0);

Error in weiner (line 275)

[X,Y,GMM]=HMRF_EM(X,Y,GMM,k,g,EM_iter,MAP_iter,beta);

I have attached the relevant files above. Any help would be very much appreciated. Thank you.

Best Answer

Evidently Y does not have the same number of elements. Put this right before the error and what does it show in the command window.

whos Y
m
n
numberOfPixels = numel(Y)
mn3 = m * n * 3

Don't use semicolons. Tell us what you see. You'll see that numberOfPixels does not equal mn3 and since they're not the same, some pixels won't get used, or you're missing some pixels. Either way, a reshape cannot happen.

Best Answer

Related Question