Rotation, Translation, Scaling of 2 clouds of 3D points

Hi

I have two sets of corresponding 3D points and want to calculate the rotation, translation and scaling between those two sets. The scaling is a result of errors made while placing the points (so the points are not scaled "in general").

My test data:

The sets each contain of 6 corresponding points which are a generating set of the R^3. They are transformed by a rotation matrix which consists of 3 rotations around the x, y and z axis:

http://www.wolframalpha.com/input/?i=%7B%7B0.93969%2C+-0.34202%2C+0%7D%2C+%7B0.34202%2C+0.93969%2C+0%7D%2C+%7B0%2C+0%2C+1%7D%7D+.+%7B%7B0.99619%2C+0%2C+0.08715%7D%2C+%7B0%2C+1%2C+0%7D%2C+%7B-0.08715%2C+0%2C+0.99619%7D%7D+.+%7B%7B1%2C+0%2C+0%7D%2C+%7B0%2C+-0.41614%2C+-0.90929%7D%2C+%7B0%2C+0.90929%2C+-0.41614%7D%7D

The rotation matrix is then multiplied with a scaling matrix:

http://www.wolframalpha.com/input/?i=%7B%7B0.93611%2C+0.216794%2C+0.276916%7D%2C+%7B0.340717%2C+-0.363939%2C+-0.866855%7D%2C+%7B-0.08715%2C+0.905826%2C+-0.414555%7D%7D.%7B%7B0.9%2C+0%2C+0%7D%2C+%7B0%2C+1.1%2C+0%7D%2C+%7B0%2C+0%2C+0.9%7D%7D

I'm not applying a translation since I've problems already using just a rotation and scaling.

I've tried to use the following matlab function but I don't get the expected results for my test data: http://www.mathworks.com/matlabcentral/fileexchange/26186-absolute-orientation-horns-method

My second approach was to use the methods described in the following publication but the results also differ from the expected results: http://cis.jhu.edu/software/lddmm-similitude/umeyama.pdf

My last approach was to solve 3 sets of linear equations to get the transformation coefficients for the transformation matrix. This approach results in the expected values but I can't find any publication or whatsoever which describes (and therefore validates) this approach.

My first question is, why the first methods calculate a single scaling value instead of 3 values (for each of the coordinate axis). Do these methods assume that the points are scaled in general? My second question would be how to solve my problem and have a publication which validates my results.

Thanks in advance. Edit: Fixed the link to the paper.

[Matt J]

My second approach was to use the methods described in the following publication

The link to the "publication" is the same as the link to the first method.

Answer 1

[Matt J]

My first question is, why the first methods calculate a single scaling value instead of 3 values

Because that is the scope of the absolute orientation problem. If you know your transformation has isotropic scaling, a single scaling value is a desired assumption of the transformation model.

My second question would be how to solve my problem and have a publication which validates my results.

If you have the Image Processing Toolbox, use CP2TFORM with the 'affine' transform type. Publications are referenced in the References section of "doc cp2tform".

[Matthias]

Thanks for your answer.

Just to get this straight: The first two (and all similar - I found lots of publications which often result in approx. the same transformation matrix) methods don't work in my case because my transformation doesn't have isotropic scaling?

I'll have a look at CP2TFORM and the publications. From what I read so far, they're using another approach to solve this problem than simply solving 3 sets of linear equations (would be great if they'd used my working approach :D).

[Matt J]

Just to get this straight: The first two ... methods don't work in my case because my transformation doesn't have isotropic scaling?

That is the reason the first method is not applicable. You didn't provide a correct link to the second method (see my Comment above) so I can only speculate on that.

From what I read so far, they're using another approach to solve this problem than simply solving 3 sets of linear equations (would be great if they'd used my working approach :D).

There are reasons why your approach is not preferred. If you're using only 3 equations, derived from only three of the N data points, it means that you're using only a small subset of the data points available for the parameter estimation. This reduces estimation accuracy.

Also, by directly solving a set of linear equations, your method does not take advantage of the a priori knowledge that the transformation has only 9 degrees of freedom (or 6 degrees if translation is omitted). This also is sub-optimal for accuracy.

[Matthias]

Sorry, i missed your comment. The correct link is: http://cis.jhu.edu/software/lddmm-similitude/umeyama.pdf ("Least squares estimation of transformation parameters between to point patterns")

I probably was not clear enough when describing my method, sorry. I'm not using 3 equations, I'm using 3 sets of n equations (n is the number of corresponding points). Maybe the following pdf helps to describe my approach: http://root.phpcoding.org/equationapproach.pdf

After calculating the coefficients I'm using a SVD on the 3x3 matrix (m_11 to m_33) to separate the rotation from the scaling component.

I'd appreciate it if you could have a short look on my approach and tell me whether this could work or is complete nonsense.

[Matt J]

Your document doesn't show what you are doing with the SVD, so I can't easily evaluate it. From what I can tell, it doesn't look like nonsense, but it doesn't look like it is optimal either. It still looks like you are solving for 12 unknowns instead of only 9. A simple test would be to take the result of cp2tform for several sets of data with significant noise and compare its squared error to the squared error given by your approach.

[Matthias]

You're right, I do solve for 12 unknowns instead of only 9. I just understood the problem regarding this issue.

Regarding the SVD Step:

SVD(M) = L * S * R^T

S is a diagonal matrix so I'm using it as the scaling matrix. L and R are unitary matrices, therefore the product of these two matrices is unitary, so I'm simply using the product as the rotation matrix.

=> Scaling = S

=> Rotation = L * R^T

I continued my search for a paper regarding this topic but I still couldn't find anything which suits to my problem. The papers either deal with rigid transformations, affine transformations with isotropic scaling (just like the paper linked in my first post) or implement some registration algorithms which focus on finding corresponding points. I can't believe that I'm the first person facing this problem...

Using CP2TFORM to test my approach is a great idea but I just realized this isn't part of the "normal" matlab application and I'd have to buy it.

Unless you have any additional thoughts on this (which I'd love to hear) thanks for all of your help, replies and time :)

[Matt J]

Scaling = S Rotation = L * R^T

I don't see the logic behind this. You were looking for a decomposition M=Scaling*Rotation and you say you've decided to do this as

M=S*(L*R^T)

but since the SVD also says that M=L*S*R^T it means you're hypothesizing an equivalence

L*S*R^T = S*L*R^T

which would be satisfied if L*S=S*L, but I'm not sure why you think the order of the multiplication doesn't matter.

I can't believe that I'm the first person facing this problem...

Obviously, you aren't the first person by far because, as we've been saying, cp2tform already does what you want. You just lack the Toolbox containing it. However, you could still read the references in "doc cpt2form"

http://www.mathworks.com/help/images/ref/cp2tform.html

and see what you can implement yourself.

Aside from that, Google shows a lot of references to 3D affine point matching. This one looks promising

http://www.ncbi.nlm.nih.gov/pubmed/20635158

but i haven't read it.

[Matthias]

Regarding the SVD: I don't think that L*S equals S*L. I just tried to calculate the rotation and scaling component in this way ("trial and error") and the results matched with the expected results (aside from a permutation of the diagonal elements of the scaling matrix which I don't care about).

I read the references in doc cpt2form but they don't really match with my problem (they refer to other sections of the doc than the affine transformation). I also read the paper you found (<http://www.ncbi.nlm.nih.gov/pubmed/20635158> - thanks for the link) but it's using an affine transformation with isotropic scaling (basically it's almost the same as "Least squares estimation of transformation parameters between to point patterns").

I just found the following modified version of the ICP algorithm which seems to use an anisotropic scaling but I haven't figured out yet how to implement this method: http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=4379798 It does look promising to me...

Answer 2

[Matt J]

You've said you don't have the Image Processing Toolbox, but if you have the Optimization Toolbox, you could use LSQNONLIN to solve it.

If you don't have any toolboxes at all, there are some nonlinear solvers on the FEX that don't seem to require anything but core MATLAB, e.g.,

http://www.mathworks.com/matlabcentral/fileexchange/39564-lmfnlsq2

http://www.mathworks.com/matlabcentral/fileexchange/17534-lmfnlsq-solution-of-nonlinear-least-squares