Optimise scaling factor for 3D points

I have two matrices of point positions in 3D space, and I want to find the optimal scaling in x,y and z directions of one set of points to minimise the difference between the sets of points (so one scaling value for each direction for all points).

I currently have written a script to do this, and a separate script to optimise the rotation of the points. The scripts iterate through several different scaling factors/rotations then plot the difference and differentiate to find the minimum.

Now I want to try and use the optimisation toolbox to optimise both simultaneously but I am struggling... I have come across the Procrustes Analysis in Matlab, but really need separate scaling factors in x, y and z.

Any help would be fab, thanks.

[Walter Roberson]

Are the matrices the same size? Is there to be a one-to-one correspondence between points, or just an overall shape?

[awinka]

Yes the matrices are the same size, and have a one-to-one correspondence.

Answer 1

[Alan Weiss]

If I understand you correctly, the transformations you are willing to perform are rotation in 3-D space, and linear scaling in each coordinate direction. Are you also willing to recenter the points, meaning add an arbitrary translation? Are you willing to reflect the points, meaning use a rotation with negative determinant?

A rotation matrix has three parameters, the three rotational angles. Translation takes another three parameters, one for each coordinate direction, if you are willing to consider this. And scaling in each coordinate direction takes another three parameters. So you might have up to 9 parameters.

Your code would look something like this for the objective function:

function f = objfn(x,M1,M2) % The two matrices are M1 and M2

theta1 = x(1);
theta2 = x(2);
theta3 = x(3);
% Do the three rotations here, multiplying M1 by the
% three rotation matrices
% Then do the three translations:
M1(:,1) = M1(:,1) + x(4);
M1(:,2) = M1(:,2) + x(5);
M1(:,3) = M1(:,3) + x(6);
% Now the three dilations
M1(:,1) = M1(:,1)*x(7);
%etc
% Then compare the two point matrices
y = norm(M1 - M2)

Optimize the function as follows:

fun = @(x)objfn(x,M1,M2);
% Make an initial point x0. Then:
x = fminunc(fun,x0);

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

[awinka]

Thanks Alan, this looks great and I can mostly follow what it's doing, but after the very first line I get this error;

>> function f = objfn(x,M1,M2) function f = objfn(x,M1,M2) | Error: Function definitions are not permitted in this context.

Any ideas? Many thanks!

[Torsten]

Did you put both parts of Alan's code in separate .m-files ?

Best wishes

Torsten.

[awinka]

No I didn't, sorry!

I've tried putting them in two separate .m files, broken down as Alan did in his post, then running the second .m file, but it states that M1 is undefined, then it runs into errors within fminunc...

Thanks, I am learning as I go!

[awinka]

I am still struggling to get this working so if anyone has any tips that would be great.

At the moment I have a file saved as "objfn.m"

function = objfn(x,M1,M2) % The two matrices are M1 and M2
x = zeros(6,1);
theta1=x(1);
theta2=x(2);
theta3=x(3);
scale1=x(4);
scale2=x(5);
scale3=x(6);
rot1=[1,0,0;0,cos(x(1)),-sin(x(1));0,sin(x(1)),cos(x(1))];
rot2=[cos(x(2)),0,sin(x(2));0,1,0;-sin(x(2)),0,cos(x(2))];
rot3=[cos(x(3)),-sin(x(3)),0;sin(x(3)),cos(x(3)),0;0,0,1];
% M1 gets modified to M3 - rotation
M1 = [0,0,0;1,1,1;2,2,2;4,4,4];
M2 = [0,0,0; 0.5,0.5,0.5;2.2,2.2,2.2;3.8,3.8,3.8];
M2 = M1 * rot1 * rot2 * rot3;
% now do scaling
M2(:,1)=M2(:,1)*x(4);
M2(:,2)=M2(:,2)*x(5);
M2(:,3)=M2(:,3)*x(6);
% compare the matrices
y=norm(M2-M1);

end

And a file "optimise.m";

    fun = @(x)objfn(x,M1,M2);
x0=[0,0,0,1,1,1];
% Make an initial point x0. Then:
% minimum x of function described in fun. 
x = fminunc(fun,x0);
  end

But when I try and run the second file it doesn't work... thanks for any suggestions.

[Alan Weiss]

You seem to have misunderstood how you pass data. Either define M1 and M2 in your workspace before you begin, or change the function definition not to require them, but as it is you both pass them in and then overwrite them in your function, which is nonsensical.

You should probably have x0 as a column vector, not a row vector.

Why do you have all those definitions of thetai and scalei when you don't use those variables? If they are for your information, comment them out.

And now for the heart of the matter: your line

M2 = M1 * rot1 * rot2 * rot3;

is incorrect. You don't want to rotate and scale M1, and then compare it to M1. Obviously, the best thing to do in that case is not to rotate or scale at all. Instead, you should use

M2 = M2 * rot1 * rot2 * rot3;

And please double-check that you are multiplying the correct way; I am not sure if you should instead be writing

M2 = rot3 * rot2 * rot1 * M2;

Alan Weiss

MATLAB mathematical toolbox documentation

[awinka]

Hi Alan, I had since noticed this error with the M2 = M1 * rot1... Yes I was hoping to define M1 and M2 in my workspace first, but this didn't seem to work so I thought maybe I needed it in the code instead.

I thought I was required to define thetai and scalei so it knew it was to be optimised? sorry if I have misunderstood? I have commented them out now. Sorry I am very new to Matlab!

I still get these errors

Error: File: objfn.m Line: 1 Column: 10 The expression to the left of the equals sign is not a valid target for an assignment.

Error in tuesday (line 8) x = fminunc(objfn,x0);

Many thanks for your help, it is much appreciated.

[Alan Weiss]

The first line should say

function y = objfn(x,M1,M2) % The two matrices are M1 and M2

Alan Weiss

MATLAB mathematical toolbox documentation

[awinka]

Many thanks, I think I finally have this working.

My final question is related to options; I currently get this message

Solver stopped prematurely.
fminunc stopped because it exceeded the function evaluation limit,
options.MaxFunEvals = 600 (the default value).

So I know I have to change the number for MaxFunEvals. I have tried doing this by adding the line above the optimisation function

options = optimoptions('fminunc', 'MaxFunEvals', 1000);

and changing the final line to

x = fminunc(opt,x0,options);

But I get the error message

Undefined function or variable 'options'.

Error in fun (line 7) options = optimset(options, 'MaxFunEvals', '1000');

Any further help would be great, thanks very much.