It's working!!!

alignShapes.m

@@ -1,4 +1,4 @@
-function [alignedShapes, avgDiff] = alignShapes(allShapes, scaling)
+function alignedShapes = alignShapes(allShapes, scaling)
 % ALIGNSHAPES uses Procrustes analysis to align a set of shapes (with or without
 % scaling).
 %
@@ -10,7 +10,6 @@
 %
 %	OUTPUT
 %       alignedShapes: The realigned shapes. Same shape as totalShapes
-%       avgDiff: The average difference (error?) between each shape and the mean shape
 %
 %   Shape analysis techniques based on this paper:
 %   Cootes, T. F., Taylor, C. J., Cooper, D. H., & Graham, J. (1995).
@@ -25,7 +24,6 @@
 % Pre-allocate
 n_shapes = size(allShapes,2);
 alignedShapes = zeros(size(allShapes));
-totalDiff = 0;
 
 % Mean shape across all subjects (assuming each shape in totalShapes is a different subj)
 meanShape = mean(allShapes,2); % x1, y1, x2, y2, ..., x20, y20
@@ -37,15 +35,11 @@
     iShape = [allShapes(1:2:end,n_shape) allShapes(2:2:end,n_shape)];
 
     % Do the Procrustes alignment
-    [d, iShapeAligned] = procrustes(meanShape,iShape,'scaling',scaling,'reflection','best');
-    totalDiff = totalDiff + d; % Total difference between iShape and the mean
+    [~, iShapeAligned] = procrustes(meanShape,iShape,'scaling',scaling,'reflection','best');
 
     % Store the aligned shape in a similar manner as how totalShapes is passed in            
     alignedShapes(1:2:end,n_shape) = iShapeAligned(:,1)';
     alignedShapes(2:2:end,n_shape) = iShapeAligned(:,2)';            
 end
 
-% On average, how far off was each shape from the mean shape?
-avgDiff = totalDiff/n_shapes;
-
 end % End of main

buildGrayLevelModel.m

@@ -19,13 +19,13 @@
 tic
 
 % Multi-resolution
-downsample_factors = 4:-1:2;
+downsample_factors = 6:-1:2;
 n_resolutions = numel(downsample_factors);
 
 % Change patterns based on resolution
 interp_step_sizes = 1*ones(n_resolutions,1); % Probably just leave this as 1
 filter_sigs = linspace(0.6,0.3,n_resolutions);
-region_size = linspace(8,3,n_resolutions);
+region_size = linspace(10,3,n_resolutions);
 
 % Build a 2D gray-level profile for each landmark at each resolution
 for n_resolution = 1:n_resolutions

buildShapeModel.m

@@ -11,10 +11,14 @@
 %       V: Eigenvectors (decreasing energy)
 %       D: Eigenvalues  (decreasing energy)
 %
+% TODO: Do I need to account for resolution somehow with the eigenvectors and values?
+%
 % John W. Miller
 % 25-Apr-2017
 
-x = alignShapes(unalignedShapes,0);
+% Align shapes from training images using Procrustes
+scaling = 0; % Yes or no, remove scale differences between training images
+x = alignShapes(unalignedShapes,scaling);
 
 % Use PCA to create model
 xBar = mean(x,2);  % Mean shape
@@ -23,10 +27,7 @@
 D = sort(diag(D),'descend');
 V = fliplr(V);
 
-% Store as struct
-
-%TODO: Do I need to account for resolution somehow with the eigenvectors and values?
-
+% Store model as a struct
 shapeModel = struct();
 shapeModel.meanShape = xBar;
 shapeModel.eVectors = V;

findFace.m

@@ -7,16 +7,16 @@
 %
 %	OUTPUT
 %
-%   TODO: Add varargin for selecting parameters (e.g. dist_metric)
+%   TODO: Add varargin for selecting parameters (e.g. dist_metric, n_evolutions, etc.)
 %
 % John W. Miller
 % 25-Apr-2017
 
+% The initial placement has to be very close to the center of the face...
+
 % Use Mahalanobis distance or gray-level PCA space for search?
 dist_metric = 'pca'; % Options: 'maha' or 'pca'
 
-% The initial placement has to be very close to the center of the face...
-
 % Inverse of covariance matrix sometimes badly scaled
 warning('off','MATLAB:nearlySingularMatrix');
 
@@ -36,10 +36,8 @@
 
     % Downsample everything that needs to be downsampled
     downsampleFactor = GM.Info.downsampleFactor;
-    im = imresize(im,1/downsampleFactor);           % Resize image
-    x_mean = shapeModel.meanShape./downsampleFactor;                % Resize mean shape
-    %         V = V./downsampleFactor;    % Is this a legit method for downsampling V and D?
-    %         D = D./downsampleFactor;
+    im = imresize(im,1/downsampleFactor);            % Resize image
+    x_mean = shapeModel.meanShape./downsampleFactor; % Resize mean shape    
 
     % Gray-level profiles model
     [g_mean, g_S, g_eVecs, g_eVals] = deal(GM.meanProfile, GM.covMatrix, GM.eVectors, GM.eValues);
@@ -62,7 +60,7 @@
 
     % Place mean shape over face (this should be done automatically)
     if n_resolution == 1
-        [x_original_estimate,h_im] = placeShape(im,x_mean);
+        [x_original_estimate,T_model_to_image, h_im] = placeShape(im,x_mean);
         %         [x_original_estimate] = asm_multiResolution(im,x_mean,GM.Info.SmoothingFilter);
         h_im = gcf; % I should be able to use h_im from above, there is a weird error
         x_aligned = x_original_estimate;
@@ -77,7 +75,7 @@
     end
 
     % Evolve estimate of face location, adjusting landmarks w/in model space
-    n_evolutions = 5;
+    n_evolutions = 4;
     for n_evolution = 1:n_evolutions;
         title(sprintf('Downsample: 1/%d. %d remaining.\nEvolution %d/%d',...
             downsampleFactor,n_resolutions-n_resolution,n_evolution,n_evolutions),'fontsize',FS), drawnow('expose')
@@ -96,7 +94,7 @@
             % Shift the 2D profile around the current landmark
             for c = -(search_size/2):(search_size/2)
                 for r = -(search_size/2):(search_size/2)
-                    n_shift = n_shift+1; g_new = [];
+                    n_shift = n_shift+1;
                     iPixel = iPixel_startingPosition+[r c]'; % Coordinates of current pixel
 %                     plot(iPixel(1),iPixel(2),'bs'), hold on
 
@@ -119,7 +117,7 @@
                     % Sigmoid equalization
                     im_region_filt = im_region_filt./(abs(im_region_filt)+C);
 
-                    % Store 2D profile as a 1D vector
+                    % Store the current 2D profile as a 1D vector
                     g_new = reshape(im_region_filt,size(im_region_filt,1).^2,1);
 
                     % Compute the 'distance' from the current profile to the mean profile                    
@@ -163,47 +161,30 @@
 
         end % Looping through landmarks
 
-        %% Remove translation, rotation, and scaling (pose parameters)
-        % from suggested shape. In other words, we're using Procrustes to project the
-        % suggested shape from image space into the shape model space, allowing us to
-        % subtract between the two. We will project back to the image space at the
-        % end of this loop.
-
-        % Need to change array dimensions for transformation
-        xm = [x_mean(1:2:end) x_mean(2:2:end)]; xs = [x_suggested(1:2:end) x_suggested(2:2:end)];
-
-        % Remove pose parameters from suggested shape
-        [~, xs, T_image_to_shape] = procrustes(xm,xs);
-        x_suggested_pose_removed = zeros(size(x_mean));
-        x_suggested_pose_removed(1:2:end) = xs(:,1); x_suggested_pose_removed(2:2:end) = xs(:,2);
-
-        % Determine weights to adjust shape w/in model space
-        % TODO: Do I somehow need to account for resolution scaling w/ the
-        % eigenvectors and values?
-        b = P'*(x_suggested_pose_removed-x_mean); % What should the subtraction be between?
-        %         [~,xc] = procrustes(x_mean,x_current);
-        %         dx = x_current-x_suggested_pose_removed;
-        dx = x_mean-x_suggested_pose_removed;
-        db = P'*dx;
-        b = b+db;
+        % Update pose parameters towards suggested shape
+        xs = [x_suggested(1:2:end) x_suggested(2:2:end)]; xc = [x_current(1:2:end) x_current(2:2:end)];        
+        [~,xp] = procrustes(xs,xc);        
+        x_posed = zeros(size(x_mean));
+        x_posed(1:2:end) = xp(:,1);
+        x_posed(2:2:end) = xp(:,2);
 
-        % Limit the model parameters based on what is considered a nomal
-        % contour, using the eigenvalues of the PCA-model
-        b=max(min(b,maxb),-maxb);
+        % Deform shape towards suggested points        
+        b_suggested = P'*(x_posed-x_mean);
+        b=max(min(b_suggested,maxb),-maxb); % Keep adjustments within model limits
 
-        % Generate new shape within constraints of model
-        %Wb = diag(D); % Eq. 26 (can also just be identity)
-        x_new = x_mean + P*(b);
+        % Generate new shape (within model space)
+        x_new = x_mean + P*b;
 
-        % Restore the pose parameters
-        xn = [x_new(1:2:end) x_new(2:2:end)];
-        xn = xn./T_image_to_shape.b/T_image_to_shape.T - T_image_to_shape.c;
+        % Transfer x_new to image space (for some reason we need to change the array size)        
+        xn = [x_new(1:2:end) x_new(2:2:end)];                        
+        [~,xn] = procrustes(xs,xn);
         x_new = zeros(size(x_mean));
-        x_new(1:2:end) = xn(:,1); x_new(2:2:end) = xn(:,2);
-        x_current = x_new; % Update the current shape position
+        x_new(1:2:end) = xn(:,1); x_new(2:2:end) = xn(:,2);                
+        x_current = x_new;
 
-        imshow(im,[]), plotLandmarks(x_current,'hold',1), hold on
-        plot(x_suggested(1:2:end),x_suggested(2:2:end),'ro','linewidth',2)
+        % View the current evolution
+        imshow(im,[]), plotLandmarks(x_current,'hold',1), hold on        
+        plot(x_suggested(1:2:end),x_suggested(2:2:end),'bo','linewidth',2)
     end % End looping through evolutions
 end % End looping through resolution levels
 

shapeModel_CootesTaylor.m

@@ -43,14 +43,14 @@
 
 %% Edge detection using ASMs
 imDir = fullfile(projectDir,'Faces','faces_B');
-n_im = 40;
+n_im = 35;
 imFile = sprintf('B_%02d_0.jpg',n_im);
 im = imread(fullfile(imDir,imFile));
 figure(34), hold on, imshow(im,[]), hold on
 
 %% Multi-resolution
 % Roughly align mean shape to face in image using multi-resolution
-% x_aligned = asm_multiResolution(im,alignedShapes);
+% x_aligned = asm_multiResolution(im,shapeModel.meanShape);
 x_aligned = placeShape(im,xBar);
 
 %% Find a face!