Subversion Repositories seema-scanner

function [alignment] = groupwiseOrthogonalProcrustes(P, initialAlign)

%groupwiseOrthogonalProcrustes Computes rototranslations to bring Np

% 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 Np 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);

% parameter vector contains relative axis angle and translations for all Np

% point sets

xInit = zeros(6*Np, 1);

% for i=1:Np

%     [w, theta] = rmat2axis(initialAlign(i).Rotation);

%     xInit((i-1)*6+1:(i-1)*6+3) = w*theta;

%     xInit((i-1)*6+4:(i-1)*6+6) = initialAlign(i).Translation;

% end

figure;

%options = optimset('Algorithm', 'levenberg-marquardt', 'Display', 'iter-detailed', 'OutputFcn', @outfun, 'MaxIter', 100, 'TolFun', 0, 'TolX', 0);

options = optimoptions('lsqnonlin', 'Display', 'iter-detailed');

[x, ~, ~] = lsqnonlin(@orthProcFun, xInit, [], [], options);

alignment = struct('Rotation', {}, 'Translation', {});

for i=1:Np

    wi = x((i-1)*6+1:(i-1)*6+3);

    Ri = initialAlign(i).Rotation * axis2rmat(wi, norm(wi));

    Ti =  initialAlign(i).Translation + x((i-1)*6+4:(i-1)*6+6);

    alignment(i).Rotation = Ri;

    alignment(i).Translation =Ti;

end

    % objective function

    function e = orthProcFun(x)

        Nc = size(P, 2);

        % transform all points according to current x

        Pbar = cell(size(P));

        for j=1:Np

            wi = x((j-1)*6+1:(j-1)*6+3);

            Ri = initialAlign(j).Rotation * axis2rmat(wi, norm(wi));

            Ti = initialAlign(j).Translation + x((j-1)*6+4:(j-1)*6+6);

            for k=1:Nc

                if(not(isempty(P{j,k})))

                    Pbar{j,k} = P{j,k}*Ri' + Ti';

                end

            end

        end

        % include all pairwise distances

        e = [];

        for j=1:Nc   

            % points in current cluster

            c = cat(1, Pbar{:,j});

            Ncj = size(c, 1);

            if(Ncj < 2)

                continue;

            end

            for k=1:Ncj-1

                for l=k+1:Ncj

                    dVec = c(k,:) - c(l,:);

                    e(end+1:end+3) = dVec;

                end

            end

        end

    end

    % output function called at every iteration

    function stop = outfun(x, optimValues, ~)

        %display(x);

        hold('on');

        plot(optimValues.residual);

        drawnow;

        stop = false;

    end

end

function [w, theta] = rmat2axis(R)

w = zeros(3, 1);

%theta = zeros(1, 1);

[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; 

    w(:) = -w(:); 

end

end

function R = axis2rmat(w, theta)

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);

% ensure orthonormal matrix

[U,~,V] = svd(R);

R = U*V';

end