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.

Answer 1

[Image Analyst]

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.

[Warid Islam]

Hi@ImageAnalyst,

After putting the above lines of code, I get the following result:

Name Size Bytes Class Attributes

Y 605x686 3320240 double

m =

605

n =

686

numberOfPixels =

415030

mn3 =

1245090

[Image Analyst]

So why are you trying to reshape this thing into 3 columns? If you put all m*n pixels into column 1, then there are none left to be put into columns 2 and 3. But why do you want to do that anyway? It's not a color image. What does 3 have to do with anything??? Maybe you want repmat() to replicate the pixels into 3 columns so that you have 3 identical columns? I have no idea what you want to do with the reshaping.

[Warid Islam]

Hi @ImageAnalyst,

You are right. I wanted to replicate the pixels into 3 columns. I used the wrong function as reshape(). However, I have some other issues. After applying repmat(), I get the following error:

Error using bsxfun

Non-singleton dimensions of the two input arrays must match each other.

Error in MRF_MAP (line 31)

yi=bsxfun(@minus,y,mu);

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);

[Image Analyst]

So, what are y and mu? Are they the same size? They need to be if you're going to subtract them, unless mu is a scalar.

whos y
whos mu

Why not just do

yi = y - mu;

[Warid Islam]

Hi Image Analyst,

Name Size Bytes Class Attributes

Y 605x686 3320240 double

m =

605

n =

686

numberOfPixels =

415030

mn3 =

1245090

Name Size Bytes Class Attributes

y 1815x2058 29882160 double

Name Size Bytes Class Attributes

mu 1x3 24 double

After doing yi = y - mu; I get the following error:

Matrix dimensions must agree.

Error in MRF_MAP (line 33)

yi = y - mu;

[Image Analyst]

Yes, so you see y is a 1815x2058 matrix, and mu is not. It's not even a scalar. mu is a 1x3 row vector. It's not the same shape as y or a constant (scalar). So what are you thinking here? Why do you think this should work? Please explain your thought process so we can guide you.

[Warid Islam]

Hi @Image Analyst,

I tried to implement the GMM-Based Hidden Markov Random Field from the following link:

https://www.mathworks.com/matlabcentral/fileexchange/39553-gmm-hmrf

But, I guess I followed the wrong process as it is for color segmentation and my image is grayscale.

[ibtihaj abusnyma]

Can you tell me what method you used for segmention here pleas?

%% 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')

[AHT.fatima]

hello everyone, can someone explain the code for removal of the pectoral muscle to me please .... it gives me good results but I did not understand its principle

[Image Analyst]

If you can't figure out that code, look at my well commented code here:

https://www.mathworks.com/matlabcentral/answers/418369-how-to-remove-the-pectoral-muscle-from-the-mammogram#comment_2155590