title: Point cloud alignment using principal component analysis with CGAL
body: <p>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.</p>
  
  <p>I first compute the centroid of the point set, and translate all points such that the origin becomes the new centroid.</p>
  
  <p>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:</p>
  
  <pre><code>Plane plane;
  linear_least_squares_fitting_3(points.begin(), points.end(), 
      plane, CGAL::Dimension_tag&lt;0&gt;());
  
  auto dir1 = dir2vec(plane.base1().direction());
  auto dir2 = dir2vec(plane.base2().direction());
  auto normal = dir1 ^ dir2; // cross product
  normal.normalize(); dir1.normalize(); dir2.normalize();
  </code></pre>
  
  <p>The <code>dir2vec</code> function converts a <code>CGAL::Direction_3</code> object to an equivalent <code>osg::Vec3d</code> object (I am using the OpenSceneGraph graphics engine). Finally, I rotate everything to the unit axes using the following code:</p>
  
  <pre><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 = [&amp;](Vec3d const &amp;p) {
      return p * r1 * r2 * r3;
  };
  transform(osgPoints.begin(), osgPoints.end(), osgPoints.begin(), rotate);
  </code></pre>
  
  <p>Here, <code>osgPoints</code> is an <code>vector&lt;osg::Vec3d&gt;</code>. For testing purposes, I translate the centroid of the rotated points back to original location, so both point clouds don't overlap.</p>
  
  <pre><code>Vec3d center = point2vec(centroid);
  auto tocentroid = [&amp;](Vec3d const &amp;v) {
      return v + center;
  };
  transform(osgPoints.begin(), osgPoints.end(), osgPoints.begin(), tocentroid);
  </code></pre>
  
  <p>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 <a href="http://www.kmussche.nl/PCA_Align.png" rel="nofollow">this image</a>. 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 <code>linear_least_squares_fitting_3</code> gives slightly different answers, as one of the planes is rotated a little bit with respect to the other.</p>
  
  <p><a href="http://www.kmussche.nl/PCA_Align2.png" rel="nofollow">Here is another</a> 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.</p>
  
  <p>Does anybody know why this happens, and, how I can prevent it?</p>