Point cloud alignment using principal component analysis with CGAL

Question

I have sets of randomly sampled points on the surface of 3D objects. I want to be able to compute the similarity between two different objects. To make that work, I first have to make sure that the sample points of both objects I want to compare do have the same rotation and scale. I thought I could do this by orienting the principal component axes along the x/y/z axes, and scaling such that the longest principal component does have unit length.

I first compute the centroid of the point set, and translate all points such that the origin becomes the new centroid.

I do the principal component analysis using the CGAL linear_least_squares_fitting_3 function, which gives the best fitting plane through the points. I compute the normal of this plane by taking the cross product of both base vectors:

Plane plane;
linear_least_squares_fitting_3(points.begin(), points.end(), 
    plane, CGAL::Dimension_tag<0>());

auto dir1 = dir2vec(plane.base1().direction());
auto dir2 = dir2vec(plane.base2().direction());
auto normal = dir1 ^ dir2; // cross product
normal.normalize(); dir1.normalize(); dir2.normalize();

The dir2vec function converts a CGAL::Direction_3 object to an equivalent osg::Vec3d object (I am using the OpenSceneGraph graphics engine). Finally, I rotate everything to the unit axes using the following code:

Matrixd r1, r2, r3;
r1.makeRotate(normal, Vec3d(1,0,0));
r2.makeRotate(dir1 * r1, Vec3d(0,1,0));
r3.makeRotate(dir2 * r1 * r2, Vec3d(0,0,1));
auto rotate = [&](Vec3d const &p) {
    return p * r1 * r2 * r3;
};
transform(osgPoints.begin(), osgPoints.end(), osgPoints.begin(), rotate);

Here, osgPoints is an vector<osg::Vec3d>. For testing purposes, I translate the centroid of the rotated points back to original location, so both point clouds don't overlap.

Vec3d center = point2vec(centroid);
auto tocentroid = [&](Vec3d const &v) {
    return v + center;
};
transform(osgPoints.begin(), osgPoints.end(), osgPoints.begin(), tocentroid);

To test it, I use two copies of the same point set, however one is transformed (rotated and translated). The above code should undo the rotations, however the results are not what I did expect: See this image. The red lines indicate the base vectors of the best fitting planes and their normal. It looks like that the results of both calls to linear_least_squares_fitting_3 gives slightly different answers, as one of the planes is rotated a little bit with respect to the other.

Here is another image where both objects are positioned with their centroid in the origin. It is now clearly visible that the normals and base vectors fall together, but the points do not.

Does anybody know why this happens, and, how I can prevent it?

Answer 1

Fitting a plane to a set of points leaves one degree of freedom unconstrained. The plane is free to spin around its normal and the fit is equal. I don't know anything about CGAL, but I wouldn't be surprised to discover that they just find a convenient plane when finding the fit (probably a nearest projection from the original axes of the space).

If you did real PCA on the point cloud, I don't think you'd have that problem. Alternatively, perhaps you could rescale (stretch) your data along the normal discovered by the fitting algorithm and then find another fit. If you stretch the data out sufficiently, then the first plane found shouldn't be as good a fit as some orthogonal plane.

Answer 2

It indeed seemed that CGAL does not compute all principal components, as JCooper suggested. I switched to the ALGLIB library to do the PCA and now it works.