im_class_mle.m

一个利用Matlab实现PCA（主成分分析）

function class=im_class_MLE(im,plot);
% function class=im_class_MLE(im,plot);
%
% routine for performing classification of multispectral images using
% Maximum Likelihood Estimation algorithm
%
% the input is a single [n m] or multispectral image set [m n p]
%
% user is prompt to select ROIs from image that mostly describes different
% classes
% 
% the output is a classified map CLASS
%
% created by K.Artyushkova
% 022004
%
% Kateryna Artyushkova
% Postdoctoral Scientist
% Department of Chemical and Nuclear Engineering
% The University of New Mexico
% [phone]
% [e-mail] 


% the user is prompt to select the number of classes
h = msgbox('You need to train the MLE classifier by specifying Regions of Interest on the images within areas mostly representing particular class');
uiwait(h)
N = inputdlg('How many training ROIs you want to spesify?');
n=str2double(N);
[a,b,c]=size(im);
if c==1;
% selection of roi from single image
    for i=1:n;
          [J,BW] = roifill(uint8(im));
          roi(:,:,i)=BW;
      end
else
% selection of roi from multispectral images - have to choose the image to
% draw ROIs on
      R = inputdlg('Which image from the array to use for ROIs drawing?');
      R=str2double(R);
      h = msgbox('Using mouse draw boxes on the image. These boxes will be used as traning regions');
      uiwait(h)
   
      for i=1:n;
          [J,BW] = roifill(uint8(im(:,:,R)));
          roi(:,:,i)=BW; %BW is a binary image with 1's for pixels corresponding to the ROI and 0's elsewhere.
      end
  end
close(figure(1))
h = waitbar(0,'Please wait while MLE classification is performed...');
  
for i=1:n;
    [I,J] = find(roi(:,:,i)==1); %find indices of all pixels within ROIs
    [m(i),p(i)]=size(I);
    clear v
    for j=1:m(i)
    v(j,:)=im(I(j),J(j),:); % creating ROI matrices containing pixels within ROIs
    end
    M(i,:)=mean(v); % calculating mean vectors for each ROI matrix 
    siz = [m(i) 1]; 
    clear M1
    M1=M(i,:);
    [p,q] = size(M1); 
    mind1 = (1:p)'; nind1 = (1:q)'; 
    mind1 = mind1(:,ones(1,siz(1)));
    nind1 = nind1(:,ones(1,siz(2)));
    clear Me
    Me = M1(mind1,nind1); % creating ROI mean matrix - each column has the same vector of means
    clear Vc
    Vc = v - Me; % Remove the mean matrix from the ROI matrix data
    clear cov
    cov = Vc' * Vc / (m(i)-1); % calculating covariance matrices for each ROI
    Det(i)=det(cov); % calculating the determinant of the covariance matrix for each ROI
    Inv(:,:,i)=(inv(cov)); %calculating invrse of covariance matrix for each ROI
end

% calculating discriminant function for each pixel
for i1=1:a
    for j1=1:b
        x=im(i1,j1,:);
        x=reshape(x,[1 q]); % extracting pixel value or vector
        for i2=1:n;
            g(i2) = - log(Det(i2))-(x'-M(i2,:)')' *Inv(:,:,i2)*(x'-M(i2,:)'); % ln(det(Cov)-(x-Mean)T*inv(Cov)*(x-Mean)
        end
        G=sum(g); % sum of all discriminant functions
        for i=1:n
            P(i)=g(i)/G; % calculating probability vector

        end
        [r,s]=max(P); % determining the max probbality - index s =class number

        % Class assignment
        for i=1:n
            A(i)=250/(n-i+1); % claculating vector of colors to assing (equally spaced numbers between 0 and 250)
            if s==i
             class(i1,j1)=A(i); %class assignment
         end
        end         
     end
 end
close(h)
% display classification results
if plot==1
H.Position=[421 330 369 337];
figure(H);
imagesc(class,[0 260])
else
end