findshapeinimage.m

cootes提出的ASM算法的matlab实现。在cootes用于课堂使用的源代码的基础上改写。改正了原来代码中的一些错误

function TerminateNotContinue=FindShapeInImage(StartShape,P,tEigenValues,W,ContoursEndingPoints,...
   MnNrmDrvProfiles,ProfilesCov,TrnPntsBelow,TrnPntsAbove,MaxNumPyramidLevels);
%function TerminateNotContinue=FindShapeInImage(StartShape,P,tEigenValues,W,ContoursEndingPoints,...
%   MnNrmDrvProfiles,ProfilesCov,TrnPntsBelow,TrnPntsAbove,MaxNumPyramidLevels);

%MnNrmDrvProfiles(level,landmark,mn_nrm_grd_profile)
%ProfilesCov{level,landmark}
%  global DEBUG_00;

ButtonName='Yes';
%def={'1','0','0','0',num2str(3*TrnPntsAbove),num2str(3*TrnPntsBelow),MaxNumPyramidLevels,'40'};
def={'1','0','0','0',num2str(3*TrnPntsAbove),num2str(3*TrnPntsBelow),num2str(MaxNumPyramidLevels),'40'};
while ButtonName=='Yes',
   
   %%%%%%load the image file. 'the shape is hiding some where in the image file'
   NumEigen=size(P,2);
   FileName=0;PathName=0;
   [FileName,PathName]=uigetfile('*.bmp','ASM: choose image file to find shape in');
   if FileName==0
      TerminateNotContinue = 1;
      return;
   end
   Img=double(imread([PathName,FileName]));      
   
   %%%%%%%%%%% data for search %%%%%%%%%%%%%
   SatisfactoryInit='No';
   while (strcmp(SatisfactoryInit,'No'))
      prompt={'Enter the initial scaling factor (>0):',...
            'Enter the initial rotation in degrees (0-360):',...
            'Enter the initial x(right) translation :',...
            'Enter the initial y(down)  translation :',...
            ['Enter number of above points for search (above+below > ',num2str(TrnPntsBelow+TrnPntsAbove),'):'],...
            'Enter number of below points for search :',...
            ['Enter number MR levels (<= ',num2str(MaxNumPyramidLevels),'):'],...
            'Enter the max number of loops for image search (>1):'};
      TheTitle='ASM';
      lineNo=[1,1,1,1,1,1,1,1];
      answer=inputdlg(prompt,TheTitle,lineNo,def);
      def=answer;
      if isempty(answer)
         TerminateNotContinue = 1;
         return;
      end   
      in_s=str2num(answer{1});
      in_Theta=str2num(answer{2})*pi/180;
      in_xc=str2num(answer{3});
      in_yc=str2num(answer{4});
      SrchPntsAbove  = str2num(answer{5});
      SrchPntsBelow  = str2num(answer{6});
      UsedPyramidLevels= str2num(answer{7});
      MAX_SEARCH_LOOPS=str2num(answer{8});
      %%%%%%%%%%% end of data for search %%%%%%%     
      
      
      %%% displaying the chosen image with the initiliazation      
      % ..1.. start with x (ex.MeanShape)
      x = StartShape;
      % ..2.. choose s,Theta,Xc  &  b=zeros
      s = in_s;
      Theta = in_Theta;
      Xc = [in_xc * ones(length(x)/2,1) ; in_yc * ones(length(x)/2,1)];
      b=zeros(NumEigen,1);
      
      X=ScaleRotateTranslate(x,s,Theta,Xc(1),Xc(end));
      %rm
      figure
      imagesc(Img);
      colormap('gray');
      PlotShapes(X,'as if we started in HI res',ContoursEndingPoints);
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      
      
      %%% Asking the user if the initialization looks OK
      SatisfactoryInit=questdlg('Is the initialization satisfactory?','ASM: Searching Images');
      if(strcmp(SatisfactoryInit,'Cancel'))
         TerminateNotContinue = 1;
         return; 
      end      
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      
   end %SatisfactoryInit
   
   
   ImagePyramid=GetImagePyramid(Img,UsedPyramidLevels);
   X=X/(2^UsedPyramidLevels);
   
   for ind2=UsedPyramidLevels:-1:1,%for each level
      disp(['Entering MR Level: ',num2str(ind2)]);
      s = 1;
      Theta = 0;
      Xc = [zeros(length(x)/2,1) ; zeros(length(x)/2,1)];
      mrImg=ImagePyramid{ind2};
      X=X*2;
      x=X;
      %rm
% %       if DEBUG_00==1
          figure
          subplot(1,2,1)
          imagesc(mrImg);
          PlotShapes(X,'Progress in MR',ContoursEndingPoints);
          subplot(1,2,2)
          imagesc(Img);
          PlotShapes(X*2^(ind2-1),'Progress in HiR',ContoursEndingPoints);
          colormap('gray');
          drawnow 
% %       end
      
      hwtbar = waitbar(0,['Searching level ',num2str(ind2),'. Please wait...']);
      ind1=0;
      Converged=0;
      while (ind1<MAX_SEARCH_LOOPS & Converged==0),%search loops in each level
          disp(['  search loops ',num2str(ind1)])
         ind1=ind1+1;
         % ..2.. find X=M(s,Theta)[x]+Xc (the initial shape guess)      
         X=ScaleRotateTranslate(x,s,Theta,Xc(1),Xc(end))+P*b;
%          X=ScaleRotateTranslate(x+P*b,s,Theta,Xc(1),Xc(end));
         
         %rm
%          if DEBUG_00==1
             clf
             subplot(1,2,1)
             imagesc(mrImg);
             PlotShapes(X,'Progress in MR',ContoursEndingPoints);
             subplot(1,2,2)
             imagesc(Img);
             PlotShapes(X*2^(ind2-1),'Progress in HiR',ContoursEndingPoints);
             colormap('gray');
             drawnow
%          end
         
         % ..3.. find needed changes dX to be applied to X to better fit the image. 2 ways to find dX
         [dX,Converged]=GetdX(X,mrImg,MnNrmDrvProfiles,ProfilesCov,...
            SrchPntsAbove,SrchPntsBelow,...
            TrnPntsBelow,ContoursEndingPoints,ind2);
         
         %rm
         if Converged==1
             disp('    Converged.');
         end
         
         
%          dX=LimitTheJump(dX);
         XPdX=X+dX;
         
         %rm
%          if DEBUG_00==1
%              figure
%              imagesc(mrImg);
%              colormap('gray');
%              PlotShapes(XPdX,'X+dX',ContoursEndingPoints);
%              title('pur:X+dX');
%             
%          end
         
         
         % ..4.. find 1+ds,dTheta,dXc that fit X to X+dX
         [x2New,y2New,dsP1,dTheta,tx,ty]=AlignShapeToShape(...
            XPdX(1:end/2),...
            XPdX(end/2+1:end),...
            X(1:end/2),...
            X(end/2+1:end),...
            W);
         dXc=[tx*ones(length(X)/2,1);ty*ones(length(X)/2,1)];
         
%          if DEBUG_00==1
%              hold on;
%              XTEMP=ScaleRotateTranslate(X,dsP1,dTheta,dXc(1),dXc(end));
%              xxtemp=XTEMP(1:39);
%              yytemp=XTEMP(40:78);
%              plot(xxtemp,yytemp,'y.');
%              
%              plot(x2New,y2New,'y-');
%              title('pur:X+dX ; yel:M_dsP1_dTheta[X]+dXc');
%          end
         
         
         
         % ..5.. find dx to make X fit exactly to XPdX
         % find dx to make M(s(1+ds),(Theta+dTheta))[x+dx]+(Xc+dXc) == X+dX , here X = M(s,Theta)[x]+dX
         dx=find_dx(s,dsP1,Theta,dTheta,x,dX,Xc,dXc);
%          dx=find_dx(s,dsP1,Theta,dTheta,X,dX,dXc);

%          if DEBUG_00==1
%              hold on;
%              XTEMP=ScaleRotateTranslate(x+dx,s*dsP1,Theta+dTheta,Xc(1)+dXc(1),Xc(end)+dXc(end));
%              xxtemp=XTEMP(1:39);
%              yytemp=XTEMP(40:78);
%              plot(xxtemp,yytemp,'r-');
%          end
             
             
         % ..6..find allowable dx i.e. db: changes allowed by varying the shape parameters b
         
             
         db=find_db(dx,P);         
         % ..7.. updating pose and shape parameters
         wt=1;wth=1;ws=1;Wb=eye(length(db));
         Xc = Xc + wt * dXc;
         Theta = Theta + wth * dTheta;
         %s = s * (1 + ws * ds);  %ds=dsP1-1
         s  = s * (1 + ws * (dsP1-1));
         
         b = b + Wb * db;
         
         b=LimitTheB(b,tEigenValues);
         
         waitbar((UsedPyramidLevels-ind2+1)/UsedPyramidLevels);
      end %for ind1=1:NumSearchLoops,
      X=ScaleRotateTranslate(x,s,Theta,Xc(1),Xc(end))+P*b;
%       X=ScaleRotateTranslate(x+P*b,s,Theta,Xc(1),Xc(end));
      close(hwtbar);
   end%for ind2=UsedPyramidLevels:-1:1
   
   %prepare data for return
   ShapeX=X(1:end/2);
   ShapeY=X(end/2+1:end);
   
   figure
   imagesc(Img);
   PlotShapes([ShapeX;ShapeY],'ASM: Shape resulting from image search and original image',ContoursEndingPoints);
   colormap('gray');
   
   ButtonName=questdlg('Do you want to try again?','ASM: Searching Images');
   if(strcmp(ButtonName,'Cancel'))
      TerminateNotContinue = 1;
      return; 
   elseif (strcmp(ButtonName,'No'))
      TerminateNotContinue = 0;
      return;
   end
   keyboard   
end%while button yes