Subversion Repositories seema-scanner

function [alignment] = groupwiseOrthogonalProcrustes(P, initialAlign)
%groupwiseOrthogonalProcrustes Computes rototranslations to bring N
% matched point sets into alignment using non-linear optimization.
%
% Input:
%   P: (Np x Nc) cell array containing point coordinates in matched order.
%       empty cells indicate that the corresponding point set does not
%       include a point matched to that cluster.
%   initialAlignment: a N length struct with fields Rotation (3x3 matrix)  
%       and Translation (3-vector) for each of the point sets.
%
% Output:
%   alignment: a struct with fields Rotation (3x3 matrix) and Translation
%       (3-vector) for each of the point sets.
%

Np = size(P, 1);
assert(length(initialAlign) == Np);

xInit = zeros(6*Np, 1);

for i=1:Np
    [wi, thetai] = rmat2axis(initialAlign(i).Rotation);
    xInit((i-1)*6+1:(i-1)*6+3) = wi*thetai;
    
    xInit((i-1)*6+4:(i-1)*6+6) = initialAlign(i).Translation;
end

options = optimset('Algorithm', 'levenberg-marquardt', 'Display', 'iter-detailed', 'OutputFcn', @outfun, 'MaxIter', [], 'TolFun', 0, 'TolX', 0);
[x, ~, ~] = lsqnonlin(@orthProcFun, xInit, [], [], options);

alignment = struct('Rotation', {}, 'Translation', {});
for i=1:Np
    wi = x((i-1)*6+1:(i-1)*6+3);
    Ri = axis2rmat(wi, norm(wi));
    Ti = x((i-1)*6+4:(i-1)*6+6);
    alignment(i).Rotation = Ri;
    alignment(i).Translation = Ti;
end


    % objective function
    function e = orthProcFun(x)

        Np = size(P, 1);
        Nc = size(P, 2);

        % transform all points according to current x
        Pbar = cell(size(P));
        for i=1:Np
            wi = x((i-1)*6+1:(i-1)*6+3);
            Ri = axis2rmat(wi, norm(wi));
            Ti = x((i-1)*6+4:(i-1)*6+6);
            for j=1:Nc
                if(not(isempty(P{i,j})))
                    Pbar{i,j} = P{i,j}*Ri' + Ti';
                end
            end
        end

        % include all pairwise distances, normalizing for each point sets
        e = [];
        for i=1:Np   
            % number of points in current point set
            Nci = size(cat(2, Pbar{i,:}), 2);
            for j=1:Nc
                if(not(isempty(Pbar{i,j})))
                    for k=setxor(1:Np, i)
                        if(not(isempty(Pbar{k,j})))
                            e = [e; 1/Nci * norm(Pbar{k,j} - Pbar{i,j})];
                            %e = [e; norm(Pbar{k,j} - Pbar{i,j})];
                        end      
                    end   
                end
            end
        end

    end

    % output function called at every iteration
    function stop = outfun(x, optimValues, ~)
        stop = false;
    end

end

function [w, theta] = rmat2axis(R)

w = zeros(3, size(R,3));
theta = zeros(1, size(R,3));

[V,D] = eig(R);
[~,ix] = min(abs(diag(D)-1)); 

w(:) = V(:,ix); 
t = [R(3,2)-R(2,3),R(1,3)-R(3,1),R(2,1)-R(1,2)];

theta = atan2(t*w(:),trace(R(:,:))-1);

if theta<0
    theta = -theta(1); 
    w(:) = -w(:); 
end
end

function R = axis2rmat(w, theta)

if(theta > 0.0001)
    w = w./norm(w); % w should be a unit vector
end

P = w*transpose(w);
Q = [0 -w(3) w(2);
     w(3) 0 -w(1);
     -w(2) w(1) 0];

% using Rodigues' rotation formula
R = P + (eye(3) - P)*cos(theta) + Q*sin(theta);
end