Transposed rotation matrix

[ghost]

For rigid registration in function transform_point_cloud() the rotation matrix is used incorrectly:

def transform_point_cloud(self, Y=None):
        if Y is None:
            self.TY = self.s * np.dot(self.Y, self.R) + self.t
            return
        else:
            return self.s * np.dot(Y, self.R) + self.t

In these expressions the rotation matrix should be used transposed:

def transform_point_cloud(self, Y=None):
        if Y is None:
            self.TY = self.s * np.dot(self.Y, self.R.T) + self.t
            return
        else:
            return self.s * np.dot(Y, self.R.T) + self.t

Therefore, as it is now, the solution iterates around transposed rotation matrix, and the resulting rotation matrix is also transposed.

[ghost]

Transformation with rotation matrix R x v in numpy can be nicely written using the matmul function or its shortcut @ like this:

res = (R @ v.T).T
# or
res = v @ R.T

Looking at this, the error in the code becomes obvious.

[siavashk]

Thanks. Looking at the paper and the code, you seem to be right. Can you submit a pull request to the development branch?

[siavashk]

Is this still an issue?

[ghost]

Nice way of solving issues! Now it looks like the code has no issues :)

P.s: I gave up since the code is written as an attempt to directly translate MATLAB to numpy and it is hard to read it. Simple transposition in the transform_point_cloud() would not solve the issue.

[gattia]

@VladimirEI what else is the issue beyond the transpose? Seems like this package has a fair amount of traction, so fixing any known bugs would be beneficial to many.

[siavashk]

I am not sure if there is an issue to begin with. At the time when I was writing the initial commit for this package, I found it more intuitive to work with the transposed matrix than the original. For example, i.e. it is more intuitive to transform a point cloud using np.matmul(R, Y) as opposed to np.matmul(np.transpose(Y), np.transpose(R)) in the original paper.

By the time @VladimirEI brought the issue up, I had already forgotten about it. If this is an actual issue, I would welcome a pull request.

[ghost]

Maybe there is actually no issue, probably it depends on how the result rotation matrix is used. As I remember, when I tried to use the code it resulted in transposed rotation matrix and wrong translation vector, so I did not use the code after all. What I am sure is the proper way to use rotation matrix of size [3x3] on vectors are of size [Nx3] in usual Numpy is like so:

res = (R @ v.T).T
# or
res = v @ R.T

So I thought this was the reason for my wrong results. Trying to fix was not easy as the code is using [3xN] size format for vectors, as in MATLAB, not the one adopted in Numpy.

[gattia]

Sorry all, I started this up and then didn't respond.

@siavashk, I think the point that @VladimirEI is making is embedded in your response. In transform_point_cloud() it is written as np.dot(self.Y, self.R), where self.Y are the points and self.R is the rotation matrix. This is being written backwards compared to what most people would think about when applying a rotation matrix (the way I expect it to be written is how you wrote it in your first example np.matmul(R, Y).

Conventionally, we expect a DxD rotation matrix R and a DxM list of points Y where D is the number of dimensions and M is the number of points in the dataset. Note that Y may be initialized as DxM or MxD, but in either case, when written as above np.matmul(R,Y) the Y input should be shaped DxM (regardless of whether needs to be transposed to get there or not).

When written the way it is in pycpd, we have np.matmul(Y,R) which, can also work if we transpose things (as you pointed out). So, in this case, Y should be aMxD matrix (i.e., should be Y.transpose from my example) and R should also be transposed... it will be the same shape, but if the same R matrix as in the first example is used it must be transposed to get the same result from the transformation.

E.g.,:

import numpy as np
R = np.identity(3)
R[0,1] = 0.93; R[2,1] = 0.5
Y = np.ones((3,100))

YT_conventional = np.matmul(R, Y)
YT_conventional_backwards = np.matmul(Y.T, R.T)
YT_conventional_backwards_dont_transpose_R = np.matmul(Y.T, R)

print(YT_conventional[:,0])
print(YT_conventional_backwards[0,:])
print(YT_conventional_backwards_dont_transpose_R[0,:])

In the above, the first two give the same answer, but the last one (that doesnt transpose R) is different.

Therefore, while everything works in pycpd as is becuase it uses things consistently, unless I am missing a final transpose of R before returning it to the user, then pycpd is learning/converging on R.T instead of learning a conventional R. This really isn't a problem if you are relying on pycpd to do everything, and are taking the outputted transformed point cloud from pycpd. Or you are applying that transform using functions written by pycpd. However, if you take the outputted transformation matrix and try to use it in a conventional way ( np.matmul(R,Y), or np.matmul(Y.T, R.T) ) then you will get funky results... and this will likely be confusing to users, especially if it isn't documented somewhere.

TLDR: I think that @VladimirEI makes a point that should be considered because I beleive that the returned R is actually transposed compared to a conventional rotation matrix and may be/likely is confusing.