function [X, template] = pcalib_GPA(varargin)
% function [X] = pcalib_GPA(varargin)
%
% Generalized procrustes analysis
%
% figure(1); clf;
% subplot(1,2,1);
% template = [0 1 1 0; 0 0 1 1];
% template = template(:);
% X = repmat(template, 1, 100);
% X = X + rand(size(X))/10;
% D = pcalib_randDeformMatrix2D(X);
% pcalib_plotshapematrix(D, 'g', gca);
% opts.scaling = 0;
% opts.rotation = 1;
% opts.translation = 1;
% [X, template] = pcalib_GPA('data',D, 'template', template, 'opts', opts);
% subplot(1,2,2);
% pcalib_plotshapematrix(X, 'r', gca);
%
% Dr. A. I. Hanna (2006)
error(nargchk(0,inf,nargin));
template = [];
template_flag = 0;
X = [];
opts.scaling = 1;
opts.translation = 1;
opts.rotation = 1;
maxiter = 100;
dim = 2;
if mod(length(varargin),2) ~= 0
    % input args have not com in pairs, woe is me
    error('Arguments to LISTDLG must come param/value in pairs.')
end
for i=1:2:length(varargin)
    switch lower(varargin{i})
        case 'data'
            X = varargin{i+1};
        case 'template'
            template = varargin{i+1};
            template_flag = 1;
            maxiter = 0;
        case 'opts'
            opts = varargin{i+1};
        case 'dimension'
            dim = varargin{i+1};
        otherwise
            error(['Unknown parameter name passed to pcalib_GPA.  Name was ' varargin{i}])
    end
end
if isempty(X)
    error('Data parameter is required.')
end
if size(template,2)==1
    template = reshape(template, dim, length(template)/dim);
end
% if the template is empty, then use the first shape in the data set
if isempty(template)
	template = reshape(X(:,1), dim, length(X(:,1))/dim);
end
[template, temp_centroid] = pcalib_centershape(template(:), dim);
scaling = opts.scaling;
rotation = opts.rotation;
translation = opts.translation;
temp_scale  = pcalib_shapesizemetric(template, dim);
dt = realmax;
thr = 10e-4;
iter = 0;
while (dt>thr) && (iter <= maxiter)
    for i=1:size(X,2)
        x = X(:,i);
        [x, mu] = pcalib_centershape(x, dim);
        % opts(1) is the scaling flag
        if scaling==1
           [x] = pcalib_scaleshape(x, temp_scale);
           %[x] = pcalib_scaleshape(x, 1);
        end
        % opts(2) is the rotation flag. Here R refers to the rotation
        % matrix [cos(theta) -sin(theta); sin(theta) cos(theta)] to align
        % the points.
        if rotation==1
            [x, R] = pcalib_rotateshape(x, template);
        end
        % opts(3) is the translation flag
        if translation==1
            %[x] = pcalib_transshape(x, zeros(dim,1));
            [x] = pcalib_transshape(x, temp_centroid);
        else
            [x] = pcalib_transshape(x, mu);
        end
        X(:,i) = x(:);
    end
    if template_flag==0
        [new_template, temp_centroid] = pcalib_estimate_template(X, dim);
        dt = sum(sqrt(sum((template - new_template).^2,1)));
        template = new_template;
    end
    iter = iter + 1;
end
[new_template, temp_centroid] = pcalib_estimate_template(X, dim);
dt = sum(sqrt(sum((template - new_template).^2,1)));
fprintf('pcalib_GPA: Converged to %2.4f in %d iterations\n', dt, iter);
template = mean(X,2);
template = reshape(template, dim, length(template)/dim);
return;

