Matlab calculate 3D similarity transformation. fitgeotrans for 3D

Question

How can I calculate in MatLab similarity transformation between 4 points in 3D? I can calculate transform matrix from

T*X = Xp,

but it will give me affine matrix due to small errors in points coordinates. How can I fit that matrix to similarity one? I need something like fitgeotrans, but in 3D

Thanks

Answer 1

The answer by @rayryeng is correct, given that you have a set of up to 3 points in a 3-dimensional space. If you need to transform m points in n-dimensional space (m>n), then you first need to add m-n coordinates to these m points such that they exist in m-dimensional space (i.e. the a matrix in @rayryeng becomes a square matrix)... Then the procedure described by @rayryeng will give you the exact transformation of points, you then just need to select only the coordinates of the transformed points in the original n-dimensional space.

As an example, say you want to transform the points:

(2 -2 2) -> (-3 5 -4)
(2 3 0) -> (3 4 4)
(-4 -2 5) -> (-4 -1 -2)
(-3 4 1) -> (4 0 5)
(5 -4 0) -> (-3 -2 -3)

Notice that you have m=5 points which are n=3-dimensional. So you need to add coordinates to these points such that they are n=m=5-dimensional, and then apply the procedure described by @rayryeng.

I have implemented a function that does that (find it below). You just need to organize the points such that each of the source-points is a column in a matrix u, and each of the target points is a column in a matrix v. The matrices u and v are going to be, thus, 3 by 5 each.

WARNING:

• the matrix A in the function may require A LOT of memory for moderately many points nP, because it has nP^4 elements.

• To overcome this, for square matrices u and v, you can simply use T=v*inv(u) or T=v/u in MATLAB notation.

• The code may run very slowly...

In MATLAB:

u = [2 2 -4 -3 5;-2 3 -2 4 -4;2 0 5 1 0]; % setting the set of source points
v = [-3 3 -4 4 -3;5 4 -1 0 -2;-4 4 -2 5 -3]; % setting the set of target points
T = findLinearTransformation(u,v); % calculating the transformation

You can verify that T is correct by:

I = eye(5);
uu = [u;I((3+1):5,1:5)]; % filling-up the matrix of source points so that you have 5-d points
w = T*uu; % calculating target points
w = w(1:3,1:5); % recovering the 3-d points
w - v % w should match v ... notice that the error between w and v is really small

The function that calculates the transformation matrix:

function [T,A] = findLinearTransformation(u,v)
% finds a matrix T (nP X nP) such that T * u(:,i) = v(:,i)
% u(:,i) and v(:,i) are n-dim col vectors; the amount of col vectors in u and v must match (and are equal to nP)
%
    if any(size(u) ~= size(v))
        error('findLinearTransform:u','u and v must be the same shape and size n-dim vectors');
    end
    [n,nP] = size(u); % n -> dimensionality; nP -> number of points to be transformed
    if nP > n % if the number of points to be transform exceeds the dimensionality of points
        I = eye(nP);
        u = [u;I((n+1):nP,1:nP)]; % then fill up the points to be transformed with the identity matrix
        v = [v;I((n+1):nP,1:nP)]; % as well as the transformed points
        [n,nP] = size(u);
    end
    A = zeros(nP*n,n*n);
    for k = 1:nP
        for i = ((k-1)*n+1):(k*n)
            A(i,mod((((i-1)*n+1):(i*n))-1,n*n) + 1) = u(:,k)';
        end
    end
    v = v(:);
    T = reshape(A\v, n, n).';
end