CGAL users discussion list

[Simon Giraudot <>]

Hi everyone,

After using the PCA function of CGAL on points 
(linear_least_squares_fitting_3), I noticed some weird behaviours: I 
attached a very simple source code example.

It performs a PCA to fit a plane to a set of 6 points (which happen to 
be coplanar). It seems to work well, the fitting score is 1.0 which 
makes sense as the points are coplanar (and if we display the plane, it 
indeed seems to fit the 6 points correctly).

But the plane we obtain is defined in a really weird way, the bases 
base1() and base2() have very very small norms (coordinates in e-30) and 
the orthogonal vector has a norm close to 1 but also with very small 
coordinates on the other 2 axis (whereas these coordinates should be 0 
considering the point set is exactly on a Z-plane).

If we call the method point() of the plane, we get a point with a huge X 
coordinate (in e28) and that is not on the plane (method has_on() 
returns false).

Then, we can notice that the method to_2d() also gives points with huge 
coordinates, which then gives points at (0,0,z) when we call the 
function to_3d().

My guess is that it all comes from the overly-small bases of the plane, 
but I don't know if this "bad" plane is due to the PCA or to the plane 
generation itself.

Notice that if we normalise the orthogonal vector and rebuilt the plane, 
it does not change anything (as the orthogonal vector was already close 
to unit).

The only way I found to correct the problem is to normalise the 2 bases 
and build a new plan using 3 points: centroid, centroid+base1() and 
centroid+base2(). This way, you can notice that the plane is 
well-defined and that the methods to_2d() and to_3d() give coherent results.

Is that behaviour known? It does not seem to occur often, the PCA still 
works fine on many data sets (even degenerate ones), so I have no idea 
why this particular point set produces this weird plane (that, from a 
PCA-fitting point of view, is still correct!).

Thanks!

--
Simon Giraudot
http://www-sop.inria.fr/members/Simon.Giraudot/index.html
3rd year PhD student

Titane team
Bureau Y310 Bâtiment Byron
INRIA Sophia Antipolis Méditerranée
2004 Route des Lucioles BP93
06902 SOPHIA ANTIPOLIS Cedex
Office: +33 (0)4 92 38 50 34
Mobile: +33 (0)6 07 77 13 59

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/linear_least_squares_fitting_3.h>
#include <vector>

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::FT FT;
typedef K::Line_3 Line;
typedef K::Plane_3 Plane;
typedef K::Point_3 Point;
typedef K::Vector_3 Vector;


int main (int argc, char** argv)
{
  // Define point set that creates the failure

  std::vector<Point> points;
  points.push_back (Point (0.30933860414961944,  0.20552380634418574, -0.40852958488336311));
  points.push_back (Point (0.2910569052260889,   0.23068922877300088, -0.40852958488336311));
  points.push_back (Point (0.30046051348904051,  0.21829778368242508, -0.40852958488336311));
  points.push_back (Point (0.30933860414961944, -0.20552380634418574, -0.40852958488336311));
  points.push_back (Point (0.30046051348904051, -0.21829778368242508, -0.40852958488336311));
  points.push_back (Point (0.2910569052260889,  -0.23068922877300088, -0.40852958488336311));

  // Compute plane + Centroid by PCA
  Plane plane;
  Point centroid;
  double score = linear_least_squares_fitting_3
    (points.begin (),points.end (), plane, centroid, CGAL::Dimension_tag<0> ());

  // Display infos
  std::cerr << "Fitting = " << score << std::endl;
  std::cerr << "Centroid = " << centroid << std::endl;
  std::cerr << "Plane   = " << plane << std::endl;
  std::cerr << "Normal  = " << plane.orthogonal_vector () << std::endl;
  std::cerr << "Base 1  = " << plane.base1 () << std::endl;
  std::cerr << "Base 2  = " << plane.base2 () << std::endl;
  std::cerr << "Point   = " << plane.point () << std::endl;
  if (!(plane.has_on (plane.point ())))
      std::cerr << " (point is not on the plane)" << std::endl;





  // Display coordinates + 2D projection + 3D coordinates of projection
  // (should be the same as original coordinates but ARE NOT)
  for (int i = 0; i < points.size (); ++ i)
    std::cerr << std::endl << "Point = " << points[i] << std::endl
	      << "To 2D = " << plane.to_2d (points[i]) << std::endl
	      << "To 3D = " << plane.to_3d (plane.to_2d (points[i])) << std::endl;

  Vector v = plane.orthogonal_vector ();
  v = v / std::sqrt (v * v);
  
  std::cerr << std::endl
	    << "Normalised normal = " << v << std::endl << std::endl;

  // Reminder: original plane
  std::cerr << "Original plane   = " << plane << std::endl;

  // Plane after normalising normal vector (does not work, same result,
  // vector was already close to unit)
  Plane p1 (centroid, v);
  std::cerr << "Normalised plane = " << p1 << std::endl;

  // THIS WORKS: normalising bases and constructing 3 points
  Vector v1 = plane.base1 ();
  Vector v2 = plane.base2 ();
  v1 = v1 / std::sqrt (v1 * v1);
  v2 = v2 / std::sqrt (v2 * v2);

  Point a = centroid + v1;
  Point b = centroid + v2;

  Plane p2 (centroid, a, b);
  std::cerr << "3-points plane   = " << p2 << std::endl;

  for (int i = 0; i < points.size (); ++ i)
    std::cerr << std::endl << "Point = " << points[i] << std::endl
	      << "To 2D = " << p2.to_2d(points[i]) << std::endl
	      << "To 3D = " << p2.to_3d (p2.to_2d (points[i])) << std::endl;

  return 0;
}