matrix to transform set of point to another set of points

I would like to start by saying that I'm not an expert mathematician, neither an expert in English and that's why I have so much problems at present.

What I'm trying to do is to generate 2 matrix R and T that fit the better 2 list of points 3D (I believe it would have the form B = A * R + T, with A and B matrices generated from my list of points.

I've found a lot of answer on this subject, but they are all time different, and there are only little details that change and I don't understand which method I should use. Rigid transform? Orthogonal matrix and shift? ... (My pointcloud should be only scaled, rotated and translated, no need for shearing or things like that)

To explain a little more, I have 2 points cloud; one is generated by photogrammetry, and the other is the "real point cloud". Because of photogrammetry, my reconstructed point cloud have no orientation, nor scale; and I would like to superpose them.
My idea was to pick 4 points (since this is 3D) in each pointclouds that represents the same points. Then find the transformation matrix using these points, and then, use this matrix on all my points.

But I don't find any ways to do it, and no method I tried works actually...

Can someone explain me easily what would be the steps to generate the matrix R and T?

• I want to nitpick on many things here, but until I have the time and motivation to actually nitpick, I'll give you some pointers. What you're trying to do is commonly called "point cloud registration". Automatic methods to achieve that are still a research topic. Since you can afford to manually pick some correspondences between the point clouds, it comes down to the model you want to use. Scaling simple to figure out in your scenario, and combination of rotation/translation is what many people call a "rigid transform", because you just move around an object rigidly without deformation. – N.Bach Sep 10 '18 at 19:20
• From there, the best rigid transform that matches the point clouds depend on what you want to call "best". The simplest method would use a mean least square error function. Assuming perfect match, any 3 pairs of points are enough to compute the matrices you want. I'll try to find a link explaining the computation later. – N.Bach Sep 10 '18 at 19:26
• Thanks for these precisions. I will take a look at what you said! I believe the mean least square error function will be enough in my case :) – Raph Schim Sep 11 '18 at 9:50

Let's agree on some notations so we know what we're talking about. Let $P=\{p_0,\ldots,p_n\}$ and $Q=\{q_0,\ldots,q_n\}$ be two points cloud representing the same scene/object, and assume they only differ by scale/rotation/translation. For simplicity, I'll just suppose that $p_i$ should be sent to point $q_i$ through the transform that you're looking for. In the more general context of point cloud registration, a common preliminary step consists in figuring out the parts of $P$ and $Q$ that overlaps, and which point of $P$ should be sent to which point of $Q$, but in your case we'll just assume you have that information available for enough pairs of points. I'll also assume you have perfect match between $P$ and $Q$.

Method

Say you want to transform $P$ into $Q$ through a rotation $R$, translation $T$, and scaling $\alpha$. Specifically $q_i = \alpha (R p_i+T)$.

If you get rid of $\alpha$, you'll fall back to the rigid transform case, so let's just do that. Note that in order to be a rotation matrix, $R$ has to be orthogonal, that is $R\times R^t = R^t\times R = \mathbf 1$. So in order to apply rigid transform techniques, you have to get rid of the scaling factor first.

Scaling

Given two pairs of correspondences $(p_i,q_i)$ and $(p_j,q_j)$, the scaling factor $\alpha$ must satisfy $\|q_i-q_j\| = \alpha \times \|p_i-p_j\|$. In theory you can compute $\alpha$ from any pair of correspondences, but in practice it's preferable to choose $p_i$ ($q_i$) far away from $p_j$ ($q_j$).

Rigid transform

Some more notations. Let the MLS error function be $$E^2(R,T) = \frac 1n\sum_{i=1}^n\left\| Rp_i+T - q_i \right\|^2$$ The "simple" method consists in finding $R$ and $T$ that minimize the error function $E$. Do note that the definition of $E$ depends on the correspondences $(p_i,q_i)$ you use. At the very minimum, you need three pairs of correspondences, and I'm pretty sure these three correspondences must satisfy some geometric constraint in order to yield a solution... probably something about being non-collinear, but if in doubt, just use more points. Least squares approaches usually benefit from having more samples, as long as there are no outliers.

Now as for the computation itself... if you don't mind quaternions, you can always have a look at Besl and McKay's 1992 publication, titled "Method for registration of 3-D shapes". All the formulas should be somewhere around page 5, though you may want to read what comes before that to understand the notations.

If you prefer to stick with matrices, you can just derive the error function. Apparently wikipedia calls that matrix calculus. As a freebie, here is what the derivative looks like for $\left\|Rp_i+T -q_i\right\|^2$: \begin{align*} \frac{\partial\left\|Rp_i+T -q_i\right\|^2}{\partial T} &= 2\left( T^t +(Rp_i)^t -q_i \right) \\ \frac{\partial\left\|Rp_i+T -q_i\right\|^2}{\partial R} &= 2\left( p_iT^t - p_iq_i^t \right) \end{align*} After summing, start by solving the derivative wrt $T$, you should find that $T$ is the translation between the two centers of mass after rotation was applied. Injecting that into the derivative wrt $R$, you should find that $R$ is the cross-covariance of $P$ and $Q$ multiplied by the inverse of some matrix.