From f440031f5c9488fffd7be20ac2a6b33c448ed322 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Wed, 28 Oct 2020 08:51:09 +0100
Subject: [PATCH 01/19] Fix precision of Eigen diagonalize traits

---
 Solver_interface/include/CGAL/Eigen_diagonalize_traits.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Solver_interface/include/CGAL/Eigen_diagonalize_traits.h b/Solver_interface/include/CGAL/Eigen_diagonalize_traits.h
index 60aee088ce1..9c6281e8188 100644
--- a/Solver_interface/include/CGAL/Eigen_diagonalize_traits.h
+++ b/Solver_interface/include/CGAL/Eigen_diagonalize_traits.h
@@ -65,7 +65,7 @@ private:
     {
       for(std::size_t j=i; j<dim; ++j)
       {
-        m(i,j) = static_cast<float>(cov[(dim * i) + j - ((i * (i+1)) / 2)]);
+        m(i,j) = static_cast<FT>(cov[(dim * i) + j - ((i * (i+1)) / 2)]);
 
         if(i != j)
           m(j,i) = m(i,j);

From d1b9e872b2a75ed41548383aa5a2cb98ca6ef8c4 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Wed, 28 Oct 2020 09:10:07 +0100
Subject: [PATCH 02/19] Add assertions if measure is null

---
 .../include/CGAL/PCA_util_Eigen.h                    | 12 ++++++++++++
 .../CGAL/linear_least_squares_fitting_circles_2.h    |  4 ++++
 .../CGAL/linear_least_squares_fitting_rectangles_2.h |  2 ++
 .../CGAL/linear_least_squares_fitting_segments_2.h   |  4 +++-
 .../CGAL/linear_least_squares_fitting_triangles_2.h  |  2 ++
 5 files changed, 23 insertions(+), 1 deletion(-)

diff --git a/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h b/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h
index dfabb50a0b0..81b69cd1dc1 100644
--- a/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h
+++ b/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h
@@ -89,6 +89,8 @@ assemble_covariance_matrix_3(InputIterator first,
     mass += area;
   }
 
+  CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
+
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
   covariance[0] += mass * (-1.0 * c.x() * c.x());
@@ -175,6 +177,8 @@ assemble_covariance_matrix_3(InputIterator first,
     mass += volume;
   }
 
+  CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
+
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
   covariance[0] += mass * (- c.x() * c.x());
@@ -267,6 +271,8 @@ assemble_covariance_matrix_3(InputIterator first,
     mass += area;
   }
 
+  CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
+
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
   covariance[0] += mass * (-1.0 * c.x() * c.x());
@@ -350,6 +356,8 @@ assemble_covariance_matrix_3(InputIterator first,
     mass += volume;
   }
 
+  CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
+
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
   covariance[0] += mass * (-1.0 * c.x() * c.x());
@@ -433,6 +441,8 @@ assemble_covariance_matrix_3(InputIterator first,
     mass += area;
   }
 
+  CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
+
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
   covariance[0] += mass * (-1.0 * c.x() * c.x());
@@ -589,6 +599,8 @@ assemble_covariance_matrix_3(InputIterator first,
     mass += length;
   }
 
+  CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
+
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
   covariance[0] += mass * (-1.0 * c.x() * c.x());
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_circles_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_circles_2.h
index dc45b98497e..052354745e2 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_circles_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_circles_2.h
@@ -107,6 +107,8 @@ linear_least_squares_fitting_2(InputIterator first,
     mass += area;
   }
 
+  CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
+
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
   covariance[0] += mass * (-1.0 * c.x() * c.x());
@@ -211,6 +213,8 @@ linear_least_squares_fitting_2(InputIterator first,
     mass += length;
   }
 
+  CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
+
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
   covariance[0] += mass * (-1.0 * c.x() * c.x());
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h
index 7aad5b86554..237ee868889 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h
@@ -110,6 +110,8 @@ linear_least_squares_fitting_2(InputIterator first,
     mass += area;
   }
 
+  CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
+
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
   covariance[0] += mass * (-1.0 * c.x() * c.x());
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_2.h
index f290e741eef..34738280ccb 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_2.h
@@ -67,7 +67,7 @@ linear_least_squares_fitting_2(InputIterator first,
 
   // assemble 2nd order moment about the origin.
   FT temp[4] = {1.0, 0.5, 0.5, 1.0};
-  Matrix moment = (1.0/3.0) * init_matrix<FT>(2,temp);
+  Matrix moment = FT(1.0/3.0) * init_matrix<FT>(2,temp);
 
   for(InputIterator it = first;
       it != beyond;
@@ -99,6 +99,8 @@ linear_least_squares_fitting_2(InputIterator first,
     mass += length;
   }
 
+  CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
+
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
   covariance[0] += mass * (-1.0 * c.x() * c.x());
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
index 2cff2b8081e..9185911fa61 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
@@ -108,6 +108,8 @@ linear_least_squares_fitting_2(InputIterator first,
     mass += area;
   }
 
+  CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
+
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
   covariance[0] += mass * (-1.0 * c.x() * c.x());

From e58561aecabffe828efee2440a38a5770647214c Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Wed, 28 Oct 2020 10:39:36 +0100
Subject: [PATCH 03/19] Make Eigen_diagonalize_traits compatible with epeck

---
 .../include/CGAL/Eigen_diagonalize_traits.h      | 16 +++++++++++++---
 1 file changed, 13 insertions(+), 3 deletions(-)

diff --git a/Solver_interface/include/CGAL/Eigen_diagonalize_traits.h b/Solver_interface/include/CGAL/Eigen_diagonalize_traits.h
index 9c6281e8188..eb2bb042b2f 100644
--- a/Solver_interface/include/CGAL/Eigen_diagonalize_traits.h
+++ b/Solver_interface/include/CGAL/Eigen_diagonalize_traits.h
@@ -30,6 +30,15 @@
 
 namespace CGAL {
 
+namespace internal {
+
+template <typename FT>
+struct Restricted_FT { typedef double type; };
+template <>
+struct Restricted_FT<float> { typedef float type; };
+
+}
+
 /// \ingroup PkgSolverInterfaceRef
 ///
 /// The class `Eigen_diagonalize_traits` provides an interface to the
@@ -47,14 +56,15 @@ namespace CGAL {
 template <typename FT, unsigned int dim = 3>
 class Eigen_diagonalize_traits
 {
+  typedef typename internal::Restricted_FT<FT>::type RFT;
 public:
   typedef std::array<FT, dim>                  Vector;
   typedef std::array<FT, dim*dim>              Matrix;
   typedef std::array<FT, (dim * (dim+1) / 2)>  Covariance_matrix;
 
 private:
-  typedef Eigen::Matrix<FT, dim, dim>            EigenMatrix;
-  typedef Eigen::Matrix<FT, dim, 1>              EigenVector;
+  typedef Eigen::Matrix<RFT, dim, dim>            EigenMatrix;
+  typedef Eigen::Matrix<RFT, dim, 1>              EigenVector;
 
   /// Construct the covariance matrix
   static EigenMatrix construct_covariance_matrix(const Covariance_matrix& cov)
@@ -65,7 +75,7 @@ private:
     {
       for(std::size_t j=i; j<dim; ++j)
       {
-        m(i,j) = static_cast<FT>(cov[(dim * i) + j - ((i * (i+1)) / 2)]);
+        m(i,j) = static_cast<RFT>(CGAL::to_double(cov[(dim * i) + j - ((i * (i+1)) / 2)]));
 
         if(i != j)
           m(j,i) = m(i,j);

From 35a096011ef4bb6e81691d026328da9cb12999ac Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Wed, 28 Oct 2020 10:39:49 +0100
Subject: [PATCH 04/19] Make PCA compatible with Epeck

---
 .../include/CGAL/PCA_util_Eigen.h             | 33 +++++++++++--------
 .../linear_least_squares_fitting_circles_2.h  |  6 ++--
 .../linear_least_squares_fitting_segments_2.h |  2 +-
 ...linear_least_squares_fitting_triangles_2.h |  2 +-
 ...linear_least_squares_fitting_triangles_3.h |  3 +-
 .../include/CGAL/centroid.h                   | 17 ++++------
 6 files changed, 33 insertions(+), 30 deletions(-)

diff --git a/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h b/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h
index 81b69cd1dc1..ff1ddc9b025 100644
--- a/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h
+++ b/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h
@@ -22,6 +22,12 @@ namespace CGAL {
 
 namespace internal {
 
+template <typename FT>
+FT approximate_cbrt (const FT& x)
+{
+  return static_cast<FT>(std::cbrt (CGAL::to_double(x)));
+}
+
 // assemble covariance matrix from a triangle set
 template < typename InputIterator,
            typename K >
@@ -67,7 +73,7 @@ assemble_covariance_matrix_3(InputIterator first,
                       t[0].y(), t[1].y(), t[2].y(),
                       t[0].z(), t[1].z(), t[2].z();
 
-    FT area = std::sqrt(t.squared_area());
+    FT area = CGAL::approximate_sqrt(t.squared_area());
 
     // skip zero measure primitives
     if(area == (FT)0.0)
@@ -236,15 +242,15 @@ assemble_covariance_matrix_3(InputIterator first,
                    t[1].y()-y0, t[3].y()-y0, t[5].y()-y0,
                    t[1].z()-z0, t[3].z()-z0, t[5].z()-z0};
     Matrix transformation (delta);
-    FT area = std::pow(delta[0]*delta[0] + delta[3]*delta[3] +
-                  delta[6]*delta[6],1/3.0)*std::pow(delta[1]*delta[1] +
-                  delta[4]*delta[4] + delta[7]*delta[7],1/3.0)*2 +
-                  std::pow(delta[0]*delta[0] + delta[3]*delta[3] +
-                  delta[6]*delta[6],1/3.0)*std::pow(delta[2]*delta[2] +
-                  delta[5]*delta[5] + delta[8]*delta[8],1/3.0)*2 +
-                  std::pow(delta[1]*delta[1] + delta[4]*delta[4] +
-                  delta[7]*delta[7],1/3.0)*std::pow(delta[2]*delta[2] +
-                  delta[5]*delta[5] + delta[8]*delta[8],1/3.0)*2;
+    FT area = approximate_cbrt(delta[0]*delta[0] + delta[3]*delta[3] +
+                  delta[6]*delta[6])*approximate_cbrt(delta[1]*delta[1] +
+                  delta[4]*delta[4] + delta[7]*delta[7])*2 +
+                  approximate_cbrt(delta[0]*delta[0] + delta[3]*delta[3] +
+                  delta[6]*delta[6])*approximate_cbrt(delta[2]*delta[2] +
+                  delta[5]*delta[5] + delta[8]*delta[8])*2 +
+                  approximate_cbrt(delta[1]*delta[1] + delta[4]*delta[4] +
+                  delta[7]*delta[7])*approximate_cbrt(delta[2]*delta[2] +
+                  delta[5]*delta[5] + delta[8]*delta[8])*2;
 
     // skip zero measure primitives
     if(area == (FT)0.0)
@@ -324,7 +330,7 @@ assemble_covariance_matrix_3(InputIterator first,
     const Sphere& t = *it;
 
     // defined for convenience.
-    FT radius = std::sqrt(t.squared_radius());
+    FT radius = CGAL::approximate_sqrt(t.squared_radius());
     Matrix transformation;
     transformation << radius, 0.0, 0.0,
                       0.0, radius, 0.0,
@@ -409,7 +415,7 @@ assemble_covariance_matrix_3(InputIterator first,
 
     // defined for convenience.
     // FT example = CGAL::to_double(t[0].x());
-    FT radius = std::sqrt(t.squared_radius());
+    FT radius = CGAL::approximate_sqrt(t.squared_radius());
     Matrix transformation;
     transformation << radius, 0.0,    0.0,
                       0.0,    radius, 0.0,
@@ -576,8 +582,7 @@ assemble_covariance_matrix_3(InputIterator first,
     transformation << t[0].x(), t[1].x(), 0.0,
                       t[0].y(), t[1].y(), 0.0,
                       t[0].z(), t[1].z(), 1.0;
-    using std::sqrt;
-    FT length = sqrt(t.squared_length());
+    FT length = CGAL::approximate_sqrt(t.squared_length());
 
     // skip zero measure primitives
     if(length == (FT)0.0)
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_circles_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_circles_2.h
index 052354745e2..3b4f679ce3b 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_circles_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_circles_2.h
@@ -83,7 +83,7 @@ linear_least_squares_fitting_2(InputIterator first,
 
     // defined for convenience.
     // FT example = CGAL::to_double(t[0].x());
-    FT radius = std::sqrt(t.squared_radius());
+    FT radius = CGAL::approximate_sqrt(t.squared_radius());
     FT delta[4] = {radius, 0.0,
                    0.0, radius};
     Matrix transformation = init_matrix<FT>(2,delta);
@@ -189,7 +189,7 @@ linear_least_squares_fitting_2(InputIterator first,
 
     // defined for convenience.
     // FT example = CGAL::to_double(t[0].x());
-    FT radius = std::sqrt(t.squared_radius());
+    FT radius = CGAL::approximate_sqrt(t.squared_radius());
     FT delta[4] = {radius, 0.0,
                    0.0, radius};
     Matrix transformation = init_matrix<FT>(2,delta);
@@ -199,7 +199,7 @@ linear_least_squares_fitting_2(InputIterator first,
     // Find the 2nd order moment for the circle wrt to the origin by an affine transformation.
 
     // Transform the standard 2nd order moment using the transformation matrix
-    transformation = 0.5 * length * transformation * moment * LA::transpose(transformation);
+    transformation = FT(0.5) * length * transformation * moment * LA::transpose(transformation);
 
     // Translate the 2nd order moment to the center of the circle.
     FT x0 = t.center().x();
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_2.h
index 34738280ccb..77f7626f702 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_2.h
@@ -83,7 +83,7 @@ linear_least_squares_fitting_2(InputIterator first,
                    t[0].y(), t[1].y()};
     Matrix transformation = init_matrix<FT>(2,delta);
     using std::sqrt;
-    FT length = sqrt(t.squared_length());
+    FT length = CGAL::approximate_sqrt(t.squared_length());
     CGAL_assertion(length != 0.0);
 
     // Find the 2nd order moment for the segment wrt to the origin by an affine transformation.
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
index 9185911fa61..1fc81ebebd5 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
@@ -88,7 +88,7 @@ linear_least_squares_fitting_2(InputIterator first,
                    t[1].y() - y0, t[2].y() - y0};
 
     Matrix transformation = init_matrix<FT>(2,delta);
-    FT area = 0.5 * std::abs(LA::determinant(transformation));
+    FT area = 0.5 * CGAL::abs(LA::determinant(transformation));
     CGAL_assertion(area!=0);
 
     // Find the 2nd order moment for the triangle wrt to the origin by an affine transformation.
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_3.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_3.h
index 2337428e1a8..09fabe92dae 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_3.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_3.h
@@ -77,7 +77,7 @@ linear_least_squares_fitting_3(InputIterator first,
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Segment> segments;
+   std::list<Segment> segments;
   for(InputIterator it = first;
       it != beyond;
       it++)
@@ -88,6 +88,7 @@ linear_least_squares_fitting_3(InputIterator first,
     segments.push_back(Segment(t[2],t[0]));
   }
 
+
   // compute fitting plane
   return linear_least_squares_fitting_3(segments.begin(),segments.end(),plane,c,(Segment*)nullptr,k,tag,
                                         diagonalize_traits);
diff --git a/Principal_component_analysis_LGPL/include/CGAL/centroid.h b/Principal_component_analysis_LGPL/include/CGAL/centroid.h
index 63fae0354ce..b69025fcb29 100644
--- a/Principal_component_analysis_LGPL/include/CGAL/centroid.h
+++ b/Principal_component_analysis_LGPL/include/CGAL/centroid.h
@@ -120,9 +120,7 @@ centroid(InputIterator begin,
       it++)
   {
     const Segment& s = *it;
-    using std::abs;
-    using std::sqrt;
-    FT length = sqrt(abs(s.squared_length()));
+    FT length = CGAL::approximate_sqrt(CGAL::abs(s.squared_length()));
     Point c = K().construct_midpoint_2_object()(s[0],s[1]);
     v = v + length * (c - ORIGIN);
     sum_lengths += length;
@@ -216,7 +214,7 @@ centroid(InputIterator begin,
       it++)
   {
     const Triangle& triangle = *it;
-    FT unsigned_area = std::abs(triangle.area());
+    FT unsigned_area = CGAL::abs(triangle.area());
     Point c = K().construct_centroid_2_object()(triangle[0],triangle[1],triangle[2]);
     v = v + unsigned_area * (c - ORIGIN);
     sum_areas += unsigned_area;
@@ -252,7 +250,7 @@ centroid(InputIterator begin,
       it++)
   {
     const Circle& s = *it;
-    FT radius = std::sqrt(s.squared_radius());
+    FT radius = CGAL::approximate_sqrt(s.squared_radius());
     Point c = s.center();
     v = v + radius * (c - ORIGIN);
     sum_lengths += radius;
@@ -381,7 +379,7 @@ centroid(InputIterator begin,
       it++)
   {
     const Iso_rectangle& r = *it;
-    FT unsigned_area = std::abs(r.area());
+    FT unsigned_area = CGAL::abs(r.area());
     Point c = K().construct_centroid_2_object()(r[0],r[1],r[2],r[3]);
     v = v + unsigned_area * (c - ORIGIN);
     sum_areas += unsigned_area;
@@ -442,8 +440,7 @@ centroid(InputIterator begin,
       it++)
   {
     const Segment& s = *it;
-    using std::sqrt;
-    FT length = sqrt(s.squared_length());
+    FT length = CGAL::approximate_sqrt(s.squared_length());
     Point c = CGAL::midpoint(s.source(),s.target());
     // Point c = K().construct_midpoint_3_object()(s[0],s[1]);
     //Point c = Point((s[0][0] + s[1][0])/2.0, (s[0][1] + s[1][1])/2.0, (s[0][2] + s[1][2])/2.0);
@@ -539,7 +536,7 @@ centroid(InputIterator begin,
       it++)
   {
     const Triangle& triangle = *it;
-    FT unsigned_area = std::sqrt(triangle.squared_area());
+    FT unsigned_area = CGAL::approximate_sqrt(triangle.squared_area());
     Point c = K().construct_centroid_3_object()(triangle[0],triangle[1],triangle[2]);
     v = v + unsigned_area * (c - ORIGIN);
     sum_areas += unsigned_area;
@@ -608,7 +605,7 @@ centroid(InputIterator begin,
       it++)
   {
     const Sphere& sphere = *it;
-    FT unsigned_volume = sphere.squared_radius() * std::sqrt(sphere.squared_radius());
+    FT unsigned_volume = sphere.squared_radius() * CGAL::approximate_sqrt(sphere.squared_radius());
     Point c = sphere.center();
     v = v + unsigned_volume * (c - ORIGIN);
     sum_volumes += unsigned_volume;

From e657cb7316c79cebc49808eae04a017a9d7bfaad Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Wed, 28 Oct 2020 13:03:34 +0100
Subject: [PATCH 05/19] Avoid filling temp lists using a subiterator for
 subdimensions

---
 .../linear_least_squares_fitting_cuboids_3.h  | 127 ++++++---------
 ...inear_least_squares_fitting_rectangles_2.h |  43 ++---
 .../linear_least_squares_fitting_segments_2.h |  18 +--
 .../linear_least_squares_fitting_segments_3.h |  41 ++---
 ...inear_least_squares_fitting_tetrahedra_3.h | 149 ++++++------------
 ...linear_least_squares_fitting_triangles_2.h |  39 ++---
 ...linear_least_squares_fitting_triangles_3.h |  82 +++-------
 7 files changed, 163 insertions(+), 336 deletions(-)

diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_cuboids_3.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_cuboids_3.h
index 23d80a5839f..4b552111dd2 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_cuboids_3.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_cuboids_3.h
@@ -20,6 +20,7 @@
 #include <CGAL/PCA_util.h>
 #include <CGAL/linear_least_squares_fitting_points_3.h>
 #include <CGAL/linear_least_squares_fitting_segments_3.h>
+#include <CGAL/Subiterator.h>
 
 #include <list>
 #include <iterator>
@@ -112,29 +113,22 @@ linear_least_squares_fitting_3(InputIterator first,
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Segment> segments;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Iso_cuboid& t = *it;
-    segments.push_back(Segment(t[0],t[1]));
-    segments.push_back(Segment(t[1],t[2]));
-    segments.push_back(Segment(t[2],t[3]));
-    segments.push_back(Segment(t[3],t[0]));
-    segments.push_back(Segment(t[4],t[5]));
-    segments.push_back(Segment(t[5],t[6]));
-    segments.push_back(Segment(t[6],t[7]));
-    segments.push_back(Segment(t[7],t[4]));
-    segments.push_back(Segment(t[5],t[0]));
-    segments.push_back(Segment(t[1],t[6]));
-    segments.push_back(Segment(t[2],t[7]));
-    segments.push_back(Segment(t[3],t[4]));
-  }
+  auto converter = [](const Iso_cuboid& c, std::size_t idx) -> Segment
+    {
+      if (idx < 7)
+        return Segment (c[idx], c[idx+1]);
+      if (idx < 10)
+        return Segment (c[idx - 6], c[(idx-1)%8]);
+      if (idx == 10)
+        return Segment (c[5], c[0]);
+      CGAL_assertion (idx == 11);
+      return Segment (c[7], c[4]);
+    };
 
-  // compute fitting plane
-  return linear_least_squares_fitting_3(segments.begin(),segments.end(),plane,c,(Segment*)nullptr,k,tag,
-                                        diagonalize_traits);
+  return linear_least_squares_fitting_3
+    (make_subiterator<Segment, 12> (first, converter),
+     make_subiterator<Segment, 12> (beyond),
+     plane,c,(Segment*)nullptr,k,tag,diagonalize_traits);
 
 } // end linear_least_squares_fitting_cuboids_3
 
@@ -154,29 +148,16 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Point_3 Point;
   typedef typename K::Iso_cuboid_3 Iso_cuboid;
+  auto converter = [](const Iso_cuboid& c, std::size_t idx) -> Point { return c[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Point> points;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Iso_cuboid& t = *it;
-    points.push_back(t[0]);
-    points.push_back(t[1]);
-    points.push_back(t[2]);
-    points.push_back(t[3]);
-    points.push_back(t[4]);
-    points.push_back(t[5]);
-    points.push_back(t[6]);
-    points.push_back(t[7]);
-  }
-
-  // compute fitting plane
-  return linear_least_squares_fitting_3(points.begin(),points.end(),plane,c,(Point*)nullptr,k,tag,
-                                        diagonalize_traits);
+  return linear_least_squares_fitting_3
+    (make_subiterator<Point, 8> (first, converter),
+     make_subiterator<Point, 8> (beyond),
+     plane,c,(Point*)nullptr,k,tag,
+     diagonalize_traits);
 
 } // end linear_least_squares_fitting_cuboids_3
 
@@ -259,33 +240,26 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Segment_3 Segment;
   typedef typename K::Iso_cuboid_3 Iso_cuboid;
+  auto converter = [](const Iso_cuboid& c, std::size_t idx) -> Segment
+    {
+      if (idx < 7)
+        return Segment (c[idx], c[idx+1]);
+      if (idx < 10)
+        return Segment (c[idx - 6], c[(idx-1)%8]);
+      if (idx == 10)
+        return Segment (c[5], c[0]);
+      CGAL_assertion (idx == 11);
+      return Segment (c[7], c[4]);
+    };
+
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Segment> segments;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Iso_cuboid& t = *it;
-    segments.push_back(Segment(t[0],t[1]));
-    segments.push_back(Segment(t[1],t[2]));
-    segments.push_back(Segment(t[2],t[3]));
-    segments.push_back(Segment(t[3],t[0]));
-    segments.push_back(Segment(t[4],t[5]));
-    segments.push_back(Segment(t[5],t[6]));
-    segments.push_back(Segment(t[6],t[7]));
-    segments.push_back(Segment(t[7],t[4]));
-    segments.push_back(Segment(t[5],t[0]));
-    segments.push_back(Segment(t[1],t[6]));
-    segments.push_back(Segment(t[2],t[7]));
-    segments.push_back(Segment(t[3],t[4]));
-  }
-
-  // compute fitting line
-  return linear_least_squares_fitting_3(segments.begin(),segments.end(),line,c,(Segment*)nullptr,k,tag,
-                                        diagonalize_traits);
+  return linear_least_squares_fitting_3
+    (make_subiterator<Segment, 12> (first, converter),
+     make_subiterator<Segment, 12> (beyond),
+     line,c,(Segment*)nullptr,k,tag,diagonalize_traits);
 
 } // end linear_least_squares_fitting_cuboids_3
 
@@ -305,29 +279,16 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Point_3 Point;
   typedef typename K::Iso_cuboid_3 Iso_cuboid;
+  auto converter = [](const Iso_cuboid& c, std::size_t idx) -> Point { return c[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Point> points;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Iso_cuboid& t = *it;
-    points.push_back(t[0]);
-    points.push_back(t[1]);
-    points.push_back(t[2]);
-    points.push_back(t[3]);
-    points.push_back(t[4]);
-    points.push_back(t[5]);
-    points.push_back(t[6]);
-    points.push_back(t[7]);
-  }
-
-  // compute fitting line
-  return linear_least_squares_fitting_3(points.begin(),points.end(),line,c,(Point*)nullptr,k,tag,
-                                        diagonalize_traits);
+  return linear_least_squares_fitting_3
+    (make_subiterator<Point, 8> (first, converter),
+     make_subiterator<Point, 8> (beyond),
+     line,c,(Point*)nullptr,k,tag,
+     diagonalize_traits);
 
 } // end linear_least_squares_fitting_cuboids_3
 
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h
index 237ee868889..579287ce003 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h
@@ -19,6 +19,7 @@
 #include <CGAL/centroid.h>
 #include <CGAL/Linear_algebraCd.h>
 #include <CGAL/PCA_util.h>
+#include <CGAL/Subiterator.h>
 
 #include <iterator>
 #include <list>
@@ -157,25 +158,16 @@ linear_least_squares_fitting_2(InputIterator first,
   // types
   typedef typename K::Iso_rectangle_2 Iso_rectangle;
   typedef typename K::Segment_2         Segment_2;
+  auto converter = [](const Iso_rectangle& r, std::size_t idx) -> Segment_2 { return Segment_2(r[idx], r[(idx+1)%4]); };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Segment_2> segments;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Iso_rectangle& t = *it;
-    segments.push_back(Segment_2(t[0],t[1]));
-    segments.push_back(Segment_2(t[1],t[2]));
-    segments.push_back(Segment_2(t[2],t[3]));
-    segments.push_back(Segment_2(t[3],t[0]));
-  }
-
-  return linear_least_squares_fitting_2(segments.begin(),segments.end(),line,c,
-                                        (Segment_2*)nullptr, K(),tag,
-                                        diagonalize_traits);
+  return linear_least_squares_fitting_2
+    (make_subiterator<Segment_2, 4> (first, converter),
+     make_subiterator<Segment_2, 4> (beyond),
+     line,c,(Segment_2*)nullptr,K(),tag,
+     diagonalize_traits);
 
 } // end linear_least_squares_fitting_2 for rectangle set with 1D tag
 
@@ -196,25 +188,16 @@ linear_least_squares_fitting_2(InputIterator first,
   // types
   typedef typename K::Iso_rectangle_2 Iso_rectangle;
   typedef typename K::Point_2         Point_2;
+  auto converter = [](const Iso_rectangle& r, std::size_t idx) -> Point_2 { return r[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Point_2> points;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Iso_rectangle& t = *it;
-    points.push_back(Point_2(t[0]));
-    points.push_back(Point_2(t[1]));
-    points.push_back(Point_2(t[2]));
-    points.push_back(Point_2(t[3]));
-  }
-
-  return linear_least_squares_fitting_2(points.begin(),points.end(),line,c,
-                                        (Point_2*)nullptr, K(),tag,
-                                        diagonalize_traits);
+  return linear_least_squares_fitting_2
+    (make_subiterator<Point_2, 4> (first, converter),
+     make_subiterator<Point_2, 4> (beyond),
+     line,c,(Point_2*)nullptr,K(),tag,
+     diagonalize_traits);
 
 } // end linear_least_squares_fitting_2 for rectangle set with 0D tag
 
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_2.h
index 77f7626f702..aa3a3432fdf 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_2.h
@@ -19,6 +19,7 @@
 #include <CGAL/centroid.h>
 #include <CGAL/Linear_algebraCd.h>
 #include <CGAL/PCA_util.h>
+#include <CGAL/Subiterator.h>
 
 #include <iterator>
 #include <list>
@@ -146,21 +147,16 @@ linear_least_squares_fitting_2(InputIterator first,
   // types
   typedef typename K::Point_2  Point;
   typedef typename K::Segment_2 Segment;
+  auto converter = [](const Segment& s, std::size_t idx) -> Point { return s[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Point> points;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Segment& s = *it;
-    points.push_back(s[0]);
-    points.push_back(s[1]);
-  }
-  return linear_least_squares_fitting_2(points.begin(),points.end(),line,c,(Point*)nullptr,k,tag,
-                                        diagonalize_traits);
+  return linear_least_squares_fitting_2
+    (make_subiterator<Point, 2> (first, converter),
+     make_subiterator<Point, 2> (beyond),
+     line,c,(Point*)nullptr,k,tag,
+     diagonalize_traits);
 
 } // end linear_least_squares_fitting_2 for segment set with 1D tag
 
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_3.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_3.h
index 79a22b08cec..7fababd309b 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_3.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_3.h
@@ -19,6 +19,7 @@
 #include <CGAL/centroid.h>
 #include <CGAL/PCA_util.h>
 #include <CGAL/linear_least_squares_fitting_points_3.h>
+#include <CGAL/Subiterator.h>
 
 #include <list>
 #include <iterator>
@@ -74,24 +75,16 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Segment_3  Segment;
   typedef typename K::Point_3  Point;
+  auto converter = [](const Segment& s, std::size_t idx) -> Point { return s[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Point> points;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Segment& t = *it;
-    points.push_back(t[0]);
-    points.push_back(t[1]);
-  }
-
-  // compute fitting plane
-  return linear_least_squares_fitting_3(points.begin(),points.end(),plane,c,(Point*)nullptr,k,tag,
-                                        diagonalize_traits);
-
+  return linear_least_squares_fitting_3
+    (make_subiterator<Point, 2> (first, converter),
+     make_subiterator<Point, 2> (beyond),
+     plane,c,(Point*)nullptr,k,tag,
+     diagonalize_traits);
 } // end linear_least_squares_fitting_segments_3
 
 // fits a line to a 3D segment set
@@ -141,24 +134,16 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Segment_3  Segment;
   typedef typename K::Point_3  Point;
+  auto converter = [](const Segment& s, std::size_t idx) -> Point { return s[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Point> points;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Segment& t = *it;
-    points.push_back(t[0]);
-    points.push_back(t[1]);
-  }
-
-  // compute fitting plane
-  return linear_least_squares_fitting_3(points.begin(),points.end(),line,c,(Point*)nullptr,k,tag,
-                                        diagonalize_traits);
-
+  return linear_least_squares_fitting_3
+    (make_subiterator<Point, 2> (first, converter),
+     make_subiterator<Point, 2> (beyond),
+     line,c,(Point*)nullptr,k,tag,
+     diagonalize_traits);
 } // end linear_least_squares_fitting_segments_3
 
 } // end namespace internal
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_tetrahedra_3.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_tetrahedra_3.h
index 4875556d73e..b7f36c518d2 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_tetrahedra_3.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_tetrahedra_3.h
@@ -22,6 +22,7 @@
 #include <CGAL/linear_least_squares_fitting_points_3.h>
 #include <CGAL/linear_least_squares_fitting_segments_3.h>
 #include <CGAL/linear_least_squares_fitting_triangles_3.h>
+#include <CGAL/Subiterator.h>
 
 #include <iterator>
 
@@ -75,26 +76,17 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Tetrahedron_3    Tetrahedron;
   typedef typename K::Triangle_3 Triangle;
+  auto converter = [](const Tetrahedron& t, std::size_t idx) -> Triangle
+    { return Triangle(t[idx], t[(idx+1)%4], t[(idx+2)%4]); };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Triangle> triangles;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Tetrahedron& t = *it;
-    triangles.push_back(Triangle(t[0],t[1],t[2]));
-    triangles.push_back(Triangle(t[0],t[2],t[3]));
-    triangles.push_back(Triangle(t[0],t[3],t[1]));
-    triangles.push_back(Triangle(t[3],t[1],t[2]));
- }
-
-  // compute fitting plane
-  return linear_least_squares_fitting_3(triangles.begin(),triangles.end(),plane,c,(Triangle*)nullptr,k,tag,
-                                        diagonalize_traits);
-
+  return linear_least_squares_fitting_3
+    (make_subiterator<Triangle, 4> (first, converter),
+     make_subiterator<Triangle, 4> (beyond),
+     plane,c,(Triangle*)nullptr,k,tag,
+     diagonalize_traits);
 } // end linear_least_squares_fitting_tetrahedrons_3
 
 // fits a plane to a 3D tetrahedron set
@@ -113,28 +105,21 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Tetrahedron_3    Tetrahedron;
   typedef typename K::Segment_3 Segment;
+  auto converter = [](const Tetrahedron& t, std::size_t idx) -> Segment
+    {
+      if (idx < 4)
+        return Segment (t[idx], t[(idx+1)%4]);
+      return Segment (t[idx-4], t[idx-2]);
+    };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Segment> segments;
-
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Tetrahedron& t = *it;
-    segments.push_back(Segment(t[0],t[1]));
-    segments.push_back(Segment(t[1],t[2]));
-    segments.push_back(Segment(t[1],t[3]));
-    segments.push_back(Segment(t[2],t[3]));
-    segments.push_back(Segment(t[0],t[2]));
-    segments.push_back(Segment(t[0],t[3]));
- }
-
-  // compute fitting plane
-  return linear_least_squares_fitting_3(segments.begin(),segments.end(),plane,c,(Segment*)nullptr,k,tag,
-                                        diagonalize_traits);
+  return linear_least_squares_fitting_3
+    (make_subiterator<Segment, 6> (first, converter),
+     make_subiterator<Segment, 6> (beyond),
+     plane,c,(Segment*)nullptr,k,tag,
+     diagonalize_traits);
 
 } // end linear_least_squares_fitting_tetrahedrons_3
 
@@ -154,26 +139,16 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Tetrahedron_3    Tetrahedron;
   typedef typename K::Point_3 Point;
+  auto converter = [](const Tetrahedron& t, std::size_t idx) -> Point { return t[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Point> points;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Tetrahedron& t = *it;
-    points.push_back(t[0]);
-    points.push_back(t[1]);
-    points.push_back(t[2]);
-    points.push_back(t[3]);
- }
-
-  // compute fitting plane
-  return linear_least_squares_fitting_3(points.begin(),points.end(),plane,c,(Point*)nullptr,k,tag,
-                                        diagonalize_traits);
-
+  return linear_least_squares_fitting_3
+    (make_subiterator<Point, 4> (first, converter),
+     make_subiterator<Point, 4> (beyond),
+     plane,c,(Point*)nullptr,k,tag,
+     diagonalize_traits);
 } // end linear_least_squares_fitting_tetrahedrons_3
 
 // fits a line to a 3D tetrahedron set
@@ -223,26 +198,17 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Tetrahedron_3    Tetrahedron;
   typedef typename K::Triangle_3 Triangle;
+  auto converter = [](const Tetrahedron& t, std::size_t idx) -> Triangle
+    { return Triangle(t[idx], t[(idx+1)%4], t[(idx+2)%4]); };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Triangle> triangles;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Tetrahedron& t = *it;
-    triangles.push_back(Triangle(t[0],t[1],t[2]));
-    triangles.push_back(Triangle(t[0],t[2],t[3]));
-    triangles.push_back(Triangle(t[0],t[3],t[1]));
-    triangles.push_back(Triangle(t[3],t[1],t[2]));
- }
-
-  // compute fitting line
-  return linear_least_squares_fitting_3(triangles.begin(),triangles.end(),line,c,(Triangle*)nullptr,k,tag,
-                                        diagonalize_traits);
-
+  return linear_least_squares_fitting_3
+    (make_subiterator<Triangle, 4> (first, converter),
+     make_subiterator<Triangle, 4> (beyond),
+     line,c,(Triangle*)nullptr,k,tag,
+     diagonalize_traits);
 } // end linear_least_squares_fitting_tetrahedrons_3
 
 // fits a line to a 3D tetrahedron set
@@ -261,28 +227,21 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Tetrahedron_3    Tetrahedron;
   typedef typename K::Segment_3 Segment;
+  auto converter = [](const Tetrahedron& t, std::size_t idx) -> Segment
+    {
+      if (idx < 4)
+        return Segment (t[idx], t[(idx+1)%4]);
+      return Segment (t[idx-4], t[idx-2]);
+    };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Segment> segments;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Tetrahedron& t = *it;
-    segments.push_back(Segment(t[0],t[1]));
-    segments.push_back(Segment(t[1],t[2]));
-    segments.push_back(Segment(t[1],t[3]));
-    segments.push_back(Segment(t[2],t[3]));
-    segments.push_back(Segment(t[0],t[2]));
-    segments.push_back(Segment(t[0],t[3]));
- }
-
-  // compute fitting line
-  return linear_least_squares_fitting_3(segments.begin(),segments.end(),line,c,(Segment*)nullptr,k,tag,
-                                        diagonalize_traits);
-
+  return linear_least_squares_fitting_3
+    (make_subiterator<Segment, 6> (first, converter),
+     make_subiterator<Segment, 6> (beyond),
+     line,c,(Segment*)nullptr,k,tag,
+     diagonalize_traits);
 } // end linear_least_squares_fitting_tetrahedrons_3
 
 // fits a line to a 3D tetrahedron set
@@ -301,26 +260,16 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Tetrahedron_3    Tetrahedron;
   typedef typename K::Point_3 Point;
+  auto converter = [](const Tetrahedron& t, std::size_t idx) -> Point { return t[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Point> points;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Tetrahedron& t = *it;
-    points.push_back(t[0]);
-    points.push_back(t[1]);
-    points.push_back(t[2]);
-    points.push_back(t[3]);
- }
-
-  // compute fitting line
-  return linear_least_squares_fitting_3(points.begin(),points.end(),line,c,(Point*)nullptr,k,tag,
-                                        diagonalize_traits);
-
+  return linear_least_squares_fitting_3
+    (make_subiterator<Point, 4> (first, converter),
+     make_subiterator<Point, 4> (beyond),
+     line,c,(Point*)nullptr,k,tag,
+     diagonalize_traits);
 } // end linear_least_squares_fitting_tetrahedra_3
 
 } // end namespace internal
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
index 1fc81ebebd5..3cd5d651b24 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
@@ -159,23 +159,16 @@ linear_least_squares_fitting_2(InputIterator first,
   // types
   typedef typename Kernel::Triangle_2 Triangle;
   typedef typename Kernel::Segment_2  Segment;
+  auto converter = [](const Triangle& t, std::size_t idx) -> Segment { return Segment(t[idx], t[(idx+1)%3]); };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Segment> segments;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Triangle& t = *it;
-    segments.push_back(Segment(t[0],t[1]));
-    segments.push_back(Segment(t[1],t[2]));
-    segments.push_back(Segment(t[2],t[0]));
-  }
-
-  return linear_least_squares_fitting_2(segments.begin(),segments.end(),line,c,tag,Kernel(),
-                                        diagonalize_traits);
+  return linear_least_squares_fitting_2
+    (make_subiterator<Segment, 3> (first, converter),
+     make_subiterator<Segment, 3> (beyond),
+     line,c,(Segment*)nullptr,Kernel(),tag,
+     diagonalize_traits);
 
 } // end linear_least_squares_fitting_2 for triangle set with 1D tag
 
@@ -196,24 +189,16 @@ linear_least_squares_fitting_2(InputIterator first,
 
   typedef typename Kernel::Triangle_2 Triangle;
   typedef typename Kernel::Point_2 Point;
+  auto converter = [](const Triangle& t, std::size_t idx) -> Point { return t[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Point> points;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Triangle& t = *it;
-    points.push_back(Point(t[0]));
-    points.push_back(Point(t[1]));
-    points.push_back(Point(t[2]));
-  }
-
-  return linear_least_squares_fitting_2(points.begin(),points.end(),line,c,tag,Kernel(),
-                                        diagonalize_traits);
-
+  return linear_least_squares_fitting_2
+    (make_subiterator<Point, 3> (first, converter),
+     make_subiterator<Point, 3> (beyond),
+     line,c,(Point*)nullptr,Kernel(),tag,
+     diagonalize_traits);
 } // end linear_least_squares_fitting_2 for triangle set with 0D tag
 
 } // end namespace internal
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_3.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_3.h
index 09fabe92dae..3c403399cbb 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_3.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_3.h
@@ -18,6 +18,7 @@
 #include <CGAL/basic.h>
 #include <CGAL/centroid.h>
 #include <CGAL/PCA_util.h>
+#include <CGAL/Subiterator.h>
 
 #include <list>
 #include <iterator>
@@ -73,25 +74,16 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Triangle_3  Triangle;
   typedef typename K::Segment_3  Segment;
+  auto converter = [](const Triangle& t, std::size_t idx) -> Segment { return Segment(t[idx], t[(idx+1)%3]); };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-   std::list<Segment> segments;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Triangle& t = *it;
-    segments.push_back(Segment(t[0],t[1]));
-    segments.push_back(Segment(t[1],t[2]));
-    segments.push_back(Segment(t[2],t[0]));
-  }
-
-
-  // compute fitting plane
-  return linear_least_squares_fitting_3(segments.begin(),segments.end(),plane,c,(Segment*)nullptr,k,tag,
-                                        diagonalize_traits);
+  return linear_least_squares_fitting_3
+    (make_subiterator<Segment, 3> (first, converter),
+     make_subiterator<Segment, 3> (beyond),
+     plane,c,(Segment*)nullptr,k,tag,
+     diagonalize_traits);
 
 } // end linear_least_squares_fitting_triangles_3
 
@@ -111,23 +103,16 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Triangle_3  Triangle;
   typedef typename K::Point_3  Point;
+  auto converter = [](const Triangle& t, std::size_t idx) -> Point { return t[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
-  std::list<Point> points;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Triangle& t = *it;
-    points.push_back(t[0]);
-    points.push_back(t[1]);
-    points.push_back(t[2]);
-  }
 
-  // compute fitting plane
-  return linear_least_squares_fitting_3(points.begin(),points.end(),plane,c,(Point*)nullptr,k,tag,
-                                        diagonalize_traits);
+  return linear_least_squares_fitting_3
+    (make_subiterator<Point, 3> (first, converter),
+     make_subiterator<Point, 3> (beyond),
+     plane,c,(Point*)nullptr,k,tag,
+     diagonalize_traits);
 
 } // end linear_least_squares_fitting_triangles_3
 
@@ -178,25 +163,16 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Triangle_3  Triangle;
   typedef typename K::Segment_3  Segment;
+  auto converter = [](const Triangle& t, std::size_t idx) -> Segment { return Segment(t[idx], t[(idx+1)%3]); };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Segment> segments;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Triangle& t = *it;
-    segments.push_back(Segment(t[0],t[1]));
-    segments.push_back(Segment(t[1],t[2]));
-    segments.push_back(Segment(t[2],t[0]));
-  }
-
-  // compute fitting line
-  return linear_least_squares_fitting_3(segments.begin(),segments.end(),line,c,(Segment*)nullptr,k,tag,
-                                        diagonalize_traits);
-
+  return linear_least_squares_fitting_3
+    (make_subiterator<Segment, 3> (first, converter),
+     make_subiterator<Segment, 3> (beyond),
+     line,c,(Segment*)nullptr,k,tag,
+     diagonalize_traits);
 } // end linear_least_squares_fitting_triangles_3
 
 // fits a line to a 3D triangle set
@@ -215,24 +191,16 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Triangle_3  Triangle;
   typedef typename K::Point_3  Point;
+  auto converter = [](const Triangle& t, std::size_t idx) -> Point { return t[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
-  std::list<Point> points;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Triangle& t = *it;
-    points.push_back(t[0]);
-    points.push_back(t[1]);
-    points.push_back(t[2]);
-  }
-
-  // compute fitting line
-  return linear_least_squares_fitting_3(points.begin(),points.end(),line,c,(Point*)nullptr,k,tag,
-                                        diagonalize_traits);
 
+  return linear_least_squares_fitting_3
+    (make_subiterator<Point, 3> (first, converter),
+     make_subiterator<Point, 3> (beyond),
+     line,c,(Point*)nullptr,k,tag,
+     diagonalize_traits);
 } // end linear_least_squares_fitting_triangles_3
 
 } // end namespace internal

From 77b2c1ecc4c0265242a763193bdf1cce5f1505d1 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Wed, 28 Oct 2020 13:05:52 +0100
Subject: [PATCH 06/19] Add test for kernels

---
 .../test_kernels.cpp                          | 100 ++++++++++++++++++
 1 file changed, 100 insertions(+)
 create mode 100644 Principal_component_analysis/test/Principal_component_analysis/test_kernels.cpp

diff --git a/Principal_component_analysis/test/Principal_component_analysis/test_kernels.cpp b/Principal_component_analysis/test/Principal_component_analysis/test_kernels.cpp
new file mode 100644
index 00000000000..b96b98e298a
--- /dev/null
+++ b/Principal_component_analysis/test/Principal_component_analysis/test_kernels.cpp
@@ -0,0 +1,100 @@
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+#include <CGAL/linear_least_squares_fitting_2.h>
+#include <CGAL/linear_least_squares_fitting_3.h>
+
+// Generate default objects so that the test compiles *AND* runs fine
+template <typename K, typename O> O default_object(const O&) { return O(); }
+template <typename K> typename K::Point_2 default_object(const typename K::Point_2&)
+{ return typename K::Point_2(0,0); }
+template <typename K> typename K::Segment_2 default_object(const typename K::Segment_2&)
+{ return typename K::Segment_2(typename K::Point_2(0,0), typename K::Point_2(0,1)); }
+template <typename K> typename K::Circle_2 default_object(const typename K::Circle_2&)
+{ return typename K::Circle_2(typename K::Point_2(0,0), typename K::FT(1.0)); }
+template <typename K> typename K::Triangle_2 default_object(const typename K::Triangle_2&)
+{ return typename K::Triangle_2(typename K::Point_2(0,0), typename K::Point_2(0,1), typename K::Point_2(1,0)); }
+template <typename K> typename K::Iso_rectangle_2 default_object(const typename K::Iso_rectangle_2&)
+{ return typename K::Iso_rectangle_2(typename K::Point_2(0,0), typename K::Point_2(1,1)); }
+template <typename K> typename K::Point_3 default_object(const typename K::Point_3&)
+{ return typename K::Point_3(0,0,0); }
+template <typename K> typename K::Segment_3 default_object(const typename K::Segment_3&)
+{ return typename K::Segment_3(typename K::Point_3(0,0,0), typename K::Point_3(0,0,1)); }
+template <typename K> typename K::Sphere_3 default_object(const typename K::Sphere_3&)
+{ return typename K::Sphere_3(typename K::Point_3(0,0,0), typename K::FT(1.0)); }
+template <typename K> typename K::Triangle_3 default_object(const typename K::Triangle_3&)
+{ return typename K::Triangle_3(typename K::Point_3(0,0,0), typename K::Point_3(0,0,1), typename K::Point_3(0,1,0)); }
+template <typename K> typename K::Tetrahedron_3 default_object(const typename K::Tetrahedron_3&)
+{ return typename K::Tetrahedron_3(typename K::Point_3(0,0,0), typename K::Point_3(0,0,1),
+                                   typename K::Point_3(0,1,0), typename K::Point_3(1,0,0)); }
+template <typename K> typename K::Iso_cuboid_3 default_object(const typename K::Iso_cuboid_3&)
+{ return typename K::Iso_cuboid_3(typename K::Point_3(0,0,0), typename K::Point_3(1,1,1)); }
+
+template <typename Kernel, typename Object, int dim>
+void test_2d()
+{
+  std::array<Object, 1> dummy = { default_object<Kernel>(Object()) };
+  typename Kernel::Line_2 line;
+  typename Kernel::Point_2 centroid;
+  CGAL::linear_least_squares_fitting_2 (dummy.begin(), dummy.end(), line, centroid,
+                                        CGAL::Dimension_tag<dim>(), Kernel(),
+                                        CGAL::Eigen_diagonalize_traits<typename Kernel::FT, 2>());
+}
+
+template <typename Kernel, typename Object, int dim>
+void test_3d()
+{
+  std::array<Object, 1> dummy = { default_object<Kernel>(Object()) };
+  typename Kernel::Line_3 line;
+  typename Kernel::Plane_3 plane;
+  typename Kernel::Point_3 centroid;
+  CGAL::linear_least_squares_fitting_3 (dummy.begin(), dummy.end(), line, centroid,
+                                        CGAL::Dimension_tag<dim>(), Kernel(),
+                                        CGAL::Eigen_diagonalize_traits<typename Kernel::FT, 3>());
+  CGAL::linear_least_squares_fitting_3 (dummy.begin(), dummy.end(), plane, centroid,
+                                        CGAL::Dimension_tag<dim>(), Kernel(),
+                                        CGAL::Eigen_diagonalize_traits<typename Kernel::FT, 3>());
+}
+
+template <typename Kernel>
+void test_kernel()
+{
+  test_2d<Kernel, typename Kernel::Point_2, 0>();
+  test_2d<Kernel, typename Kernel::Segment_2, 1>();
+  test_2d<Kernel, typename Kernel::Segment_2, 0>();
+  test_2d<Kernel, typename Kernel::Circle_2, 2>();
+  test_2d<Kernel, typename Kernel::Circle_2, 1>();
+  test_2d<Kernel, typename Kernel::Triangle_2, 2>();
+  test_2d<Kernel, typename Kernel::Triangle_2, 1>();
+  test_2d<Kernel, typename Kernel::Triangle_2, 0>();
+  test_2d<Kernel, typename Kernel::Iso_rectangle_2, 2>();
+  test_2d<Kernel, typename Kernel::Iso_rectangle_2, 1>();
+  test_2d<Kernel, typename Kernel::Iso_rectangle_2, 0>();
+  test_3d<Kernel, typename Kernel::Point_3, 0>();
+  test_3d<Kernel, typename Kernel::Segment_3, 1>();
+  test_3d<Kernel, typename Kernel::Segment_3, 0>();
+  test_3d<Kernel, typename Kernel::Sphere_3, 3>();
+  test_3d<Kernel, typename Kernel::Sphere_3, 2>();
+  test_3d<Kernel, typename Kernel::Triangle_3, 2>();
+  test_3d<Kernel, typename Kernel::Triangle_3, 1>();
+  test_3d<Kernel, typename Kernel::Triangle_3, 0>();
+  test_3d<Kernel, typename Kernel::Tetrahedron_3, 3>();
+  test_3d<Kernel, typename Kernel::Tetrahedron_3, 2>();
+  test_3d<Kernel, typename Kernel::Tetrahedron_3, 1>();
+  test_3d<Kernel, typename Kernel::Tetrahedron_3, 0>();
+  test_3d<Kernel, typename Kernel::Iso_cuboid_3, 3>();
+  test_3d<Kernel, typename Kernel::Iso_cuboid_3, 2>();
+  test_3d<Kernel, typename Kernel::Iso_cuboid_3, 1>();
+  test_3d<Kernel, typename Kernel::Iso_cuboid_3, 0>();
+}
+
+int main()
+{
+  test_kernel<CGAL::Simple_cartesian<float> >();
+  test_kernel<CGAL::Simple_cartesian<double> >();
+  test_kernel<CGAL::Exact_predicates_inexact_constructions_kernel>();
+  test_kernel<CGAL::Exact_predicates_exact_constructions_kernel>();
+
+  return EXIT_SUCCESS;
+}

From 913c29bbca543f7da094a0bb6a098bf489ee3cc0 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Thu, 29 Oct 2020 09:24:36 +0100
Subject: [PATCH 07/19] Make PCA of iso cuboid more conform to paper formulas

---
 .../include/CGAL/PCA_util_Eigen.h             | 29 ++++++++++++-------
 1 file changed, 18 insertions(+), 11 deletions(-)

diff --git a/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h b/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h
index ff1ddc9b025..b095511ba60 100644
--- a/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h
+++ b/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h
@@ -149,18 +149,25 @@ assemble_covariance_matrix_3(InputIterator first,
 
     // defined for convenience.
     // FT example = CGAL::to_double(t[0].x());
-    FT x0 = t[0].x();
-    FT y0 = t[0].y();
-    FT z0 = t[0].z();
-    FT delta[9] = {t[1].x()-x0, t[3].x()-x0, t[5].x()-x0,
-                   t[1].y()-y0, t[3].y()-y0, t[5].y()-y0,
-                   t[1].z()-z0, t[3].z()-z0, t[5].z()-z0};
-    Matrix transformation (delta);
-    FT volume = t.volume();
+    FT x0 = t.xmin();
+    FT y0 = t.ymin();
+    FT z0 = t.zmin();
+
+    FT x1 = t.xmax();
+    FT y1 = t.ymax();
+    FT z1 = t.zmax();
+
+    Matrix transformation;
+    transformation << x1 - x0, 0      , 0      ,
+                      0      , y1 - y0, 0      ,
+                      0      , 0      , z1 - z0;
+
+    FT volume = (x1-x0) * (y1-y0) * (z1-z0);
 
     // skip zero measure primitives
     if(volume == (FT)0.0)
       continue;
+    CGAL_assertion(volume > 0.0);
 
     // Find the 2nd order moment for the cuboid wrt to the origin by an affine transformation.
 
@@ -168,9 +175,9 @@ assemble_covariance_matrix_3(InputIterator first,
     transformation = volume * transformation * moment * transformation.transpose();
 
     // Translate the 2nd order moment to the minimum corner (x0,y0,z0) of the cuboid.
-    FT xav0 = (delta[0] + delta[1] + delta[2])/4.0;
-    FT yav0 = (delta[3] + delta[4] + delta[5])/4.0;
-    FT zav0 = (delta[6] + delta[7] + delta[8])/4.0;
+    FT xav0 = (x1 - x0) / (4.0);
+    FT yav0 = (y1 - y0) / (4.0);
+    FT zav0 = (z1 - z0) / (4.0);
 
     // and add to covariance matrix
     covariance[0] += transformation(0,0) + volume * (2*x0*xav0 + x0*x0);

From 5e727465eb79439c261fdbc624429b4c36743d76 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Thu, 29 Oct 2020 14:15:38 +0100
Subject: [PATCH 08/19] Fix centroid computations for Iso_cuboid_3

---
 Principal_component_analysis_LGPL/include/CGAL/centroid.h | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Principal_component_analysis_LGPL/include/CGAL/centroid.h b/Principal_component_analysis_LGPL/include/CGAL/centroid.h
index b69025fcb29..5d6bee250d2 100644
--- a/Principal_component_analysis_LGPL/include/CGAL/centroid.h
+++ b/Principal_component_analysis_LGPL/include/CGAL/centroid.h
@@ -714,7 +714,7 @@ centroid(InputIterator begin,
   {
     const Iso_cuboid& cuboid = *it;
     FT unsigned_area = 2 * ((cuboid.xmax()-cuboid.xmin())*(cuboid.ymax()-cuboid.ymin()) + (cuboid.xmax()-cuboid.xmin())*(cuboid.zmax()-cuboid.zmin()) + (cuboid.ymax()-cuboid.ymin())*(cuboid.zmax()-cuboid.zmin()));
-    Point c = K().construct_centroid_3_object()(cuboid[0],cuboid[1],cuboid[3],cuboid[5]);
+    Point c = K().construct_midpoint_3_object()(cuboid[0],cuboid[7]);
     v = v + unsigned_area * (c - ORIGIN);
     sum_areas += unsigned_area;
   }
@@ -747,7 +747,7 @@ centroid(InputIterator begin,
   {
     const Iso_cuboid& cuboid = *it;
     FT unsigned_volume = cuboid.volume();
-    Point c = K().construct_centroid_3_object()(cuboid[0],cuboid[1],cuboid[3],cuboid[5]);
+    Point c = K().construct_midpoint_3_object()(cuboid[0],cuboid[7]);
     v = v + unsigned_volume * (c - ORIGIN);
     sum_volumes += unsigned_volume;
   }

From 5c39c3bd245894c075546905c94c4eec5360d78e Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Thu, 29 Oct 2020 15:11:03 +0100
Subject: [PATCH 09/19] Fix PCA for Iso_cuboid_3 volumic

---
 Principal_component_analysis/include/CGAL/PCA_util_Eigen.h | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h b/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h
index b095511ba60..516efed7c55 100644
--- a/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h
+++ b/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h
@@ -175,9 +175,9 @@ assemble_covariance_matrix_3(InputIterator first,
     transformation = volume * transformation * moment * transformation.transpose();
 
     // Translate the 2nd order moment to the minimum corner (x0,y0,z0) of the cuboid.
-    FT xav0 = (x1 - x0) / (4.0);
-    FT yav0 = (y1 - y0) / (4.0);
-    FT zav0 = (z1 - z0) / (4.0);
+    FT xav0 = (x1 - x0) / (2.0);
+    FT yav0 = (y1 - y0) / (2.0);
+    FT zav0 = (z1 - z0) / (2.0);
 
     // and add to covariance matrix
     covariance[0] += transformation(0,0) + volume * (2*x0*xav0 + x0*x0);

From c216426908d9a13292cbe336e7b132c153ccb6a4 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Thu, 29 Oct 2020 15:16:44 +0100
Subject: [PATCH 10/19] Add missing file

---
 .../include/CGAL/Subiterator.h                | 99 +++++++++++++++++++
 1 file changed, 99 insertions(+)
 create mode 100644 Principal_component_analysis/include/CGAL/Subiterator.h

diff --git a/Principal_component_analysis/include/CGAL/Subiterator.h b/Principal_component_analysis/include/CGAL/Subiterator.h
new file mode 100644
index 00000000000..fbadd1980c2
--- /dev/null
+++ b/Principal_component_analysis/include/CGAL/Subiterator.h
@@ -0,0 +1,99 @@
+// Copyright (c) 2020  GeometryFactory (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org)
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
+//
+// Author(s) : Simon Giraudot
+//
+
+#ifndef CGAL_PCA_SUBITERATOR_H
+#define CGAL_PCA_SUBITERATOR_H
+
+#include <boost/iterator/iterator_facade.hpp>
+
+namespace CGAL
+{
+
+template <typename InputIterator, typename ValueType, int Size>
+class Subiterator
+  : public boost::iterator_facade<Subiterator<InputIterator, ValueType, Size>,
+                                  ValueType,
+                                  std::input_iterator_tag>
+{
+public:
+  using Self = Subiterator<InputIterator, ValueType, Size>;
+  using Facade = boost::iterator_facade<Self, ValueType, std::input_iterator_tag>;
+
+  using Input_type = typename std::iterator_traits<InputIterator>::value_type;
+  using Output_type = ValueType;
+
+  using Converter = std::function<Output_type(const Input_type&, std::size_t)>;
+
+private:
+
+  Converter m_converter;
+  std::size_t m_next_index;
+  InputIterator m_base;
+  mutable Output_type m_current;
+
+public:
+
+  Subiterator() { }
+
+  Subiterator(InputIterator begin, const Converter& converter)
+    : m_converter(converter), m_next_index(0), m_base(begin)
+  { }
+
+  Subiterator(InputIterator end)
+    : m_next_index(0), m_base(end)
+  { }
+
+private:
+
+  friend class boost::iterator_core_access;
+
+  void increment()
+  {
+    ++ m_next_index;
+    if (m_next_index == Size)
+    {
+      ++ m_base;
+      m_next_index = 0;
+    }
+  }
+
+  bool equal(const Self& other) const
+  {
+    return this->m_base == other.m_base && this->m_next_index == other.m_next_index;
+  }
+
+  Output_type& dereference() const
+  {
+    m_current = m_converter (*m_base, m_next_index);
+    return const_cast<Output_type&>(m_current);
+  }
+};
+
+template <typename ValueType, int Size, typename InputIterator>
+Subiterator<InputIterator, ValueType, Size>
+make_subiterator (InputIterator begin,
+                  const typename Subiterator<InputIterator, ValueType, Size>::Converter& converter)
+{
+  return Subiterator<InputIterator, ValueType, Size>(begin, converter);
+}
+
+template <typename ValueType, int Size, typename InputIterator>
+Subiterator<InputIterator, ValueType, Size>
+make_subiterator (InputIterator end)
+{
+  return Subiterator<InputIterator, ValueType, Size>(end);
+}
+
+
+} // namespace CGAL
+
+#endif // CGAL_PCA_SUBITERATOR_H

From 2bbb9f24daf32cc6c4677fa280fa2b4f5418e03d Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Mon, 2 Nov 2020 15:30:32 +0100
Subject: [PATCH 11/19] Fix Iso_cuboid_3 version with tag<2> + decrease error
 limit in tests

---
 .../include/CGAL/PCA_util_Eigen.h             | 66 +++++++++++++------
 .../test_validity.cpp                         |  2 +-
 2 files changed, 47 insertions(+), 21 deletions(-)

diff --git a/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h b/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h
index 516efed7c55..af8e8fb634e 100644
--- a/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h
+++ b/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h
@@ -17,6 +17,7 @@
 #include <CGAL/Eigen_diagonalize_traits.h>
 #include <CGAL/Default_diagonalize_traits.h>
 #include <CGAL/Dimension.h>
+#include <CGAL/Subiterator.h>
 
 namespace CGAL {
 
@@ -215,6 +216,32 @@ assemble_covariance_matrix_3(InputIterator first,
                              const CGAL::Dimension_tag<2>&,
                              const Eigen_diagonalize_traits<typename K::FT, 3>&)
 {
+#if 1
+  typedef typename K::FT              FT;
+  typedef typename K::Iso_cuboid_3    Iso_cuboid;
+  typedef typename K::Triangle_3      Triangle;
+  auto converter = [](const Iso_cuboid& c, std::size_t idx) -> Triangle
+    {
+      // Decomposition of 6 faces of the cuboid into 12 triangles
+      static constexpr std::array<std::array<std::size_t, 3>, 12 > indices
+      = {{ { 0, 1, 2 }, { 0, 2, 3 }, { 2, 3, 4 }, { 2, 4, 7 },
+           { 3, 4, 5 }, { 3, 5, 0 }, { 4, 5, 6 }, { 4, 6, 7 },
+           { 5, 6, 1 }, { 5, 1, 0 }, { 6, 7, 2 }, { 6, 2, 1 } }};
+      return Triangle (c[indices[idx][0]], c[indices[idx][1]], c[indices[idx][2]]);
+    };
+
+  assemble_covariance_matrix_3
+    (make_subiterator<Triangle, 12> (first, converter),
+     make_subiterator<Triangle, 12> (beyond),
+     covariance, c, K(), (Triangle*)nullptr, CGAL::Dimension_tag<2>(),
+     Eigen_diagonalize_traits<FT, 3>());
+
+#else
+  // This variant uses the standard formulas but seems to be broken
+  // (line/plane estimated appear to be wrong). In the absence of a
+  // reliable fix so far, the above workaround applying PCA to a
+  // decomposition of the cuboid into triangles is used.
+
   typedef typename K::FT FT;
   typedef typename K::Iso_cuboid_3 Iso_cuboid;
   typedef typename Eigen::Matrix<FT, 3, 3> Matrix;
@@ -242,26 +269,25 @@ assemble_covariance_matrix_3(InputIterator first,
     const Iso_cuboid& t = *it;
 
     // defined for convenience.
-    FT x0 = t[0].x();
-    FT y0 = t[0].y();
-    FT z0 = t[0].z();
-    FT delta[9] = {t[1].x()-x0, t[3].x()-x0, t[5].x()-x0,
-                   t[1].y()-y0, t[3].y()-y0, t[5].y()-y0,
-                   t[1].z()-z0, t[3].z()-z0, t[5].z()-z0};
-    Matrix transformation (delta);
-    FT area = approximate_cbrt(delta[0]*delta[0] + delta[3]*delta[3] +
-                  delta[6]*delta[6])*approximate_cbrt(delta[1]*delta[1] +
-                  delta[4]*delta[4] + delta[7]*delta[7])*2 +
-                  approximate_cbrt(delta[0]*delta[0] + delta[3]*delta[3] +
-                  delta[6]*delta[6])*approximate_cbrt(delta[2]*delta[2] +
-                  delta[5]*delta[5] + delta[8]*delta[8])*2 +
-                  approximate_cbrt(delta[1]*delta[1] + delta[4]*delta[4] +
-                  delta[7]*delta[7])*approximate_cbrt(delta[2]*delta[2] +
-                  delta[5]*delta[5] + delta[8]*delta[8])*2;
+    FT x0 = t.xmin();
+    FT y0 = t.ymin();
+    FT z0 = t.zmin();
+
+    FT x1 = t.xmax();
+    FT y1 = t.ymax();
+    FT z1 = t.zmax();
+
+    Matrix transformation;
+    transformation << x1 - x0, 0      , 0      ,
+                      0      , y1 - y0, 0      ,
+                      0      , 0      , z1 - z0;
+
+    FT area = FT(2) * ((x1-x0)*(y1-y0) + (x1-x0)*(z1-z0) + (y1-y0)*(z1-z0));
 
     // skip zero measure primitives
     if(area == (FT)0.0)
       continue;
+    CGAL_assertion(area > 0.0);
 
     // Find the 2nd order moment for the cuboid wrt to the origin by an affine transformation.
 
@@ -269,9 +295,9 @@ assemble_covariance_matrix_3(InputIterator first,
     transformation = area * transformation * moment * transformation.transpose();
 
     // Translate the 2nd order moment to the minimum corner (x0,y0,z0) of the cuboid.
-    FT xav0 = (delta[0] + delta[1] + delta[2])/4.0;
-    FT yav0 = (delta[3] + delta[4] + delta[5])/4.0;
-    FT zav0 = (delta[6] + delta[7] + delta[8])/4.0;
+    FT xav0 = (x1 - x0) / (2.0);
+    FT yav0 = (y1 - y0) / (2.0);
+    FT zav0 = (z1 - z0) / (2.0);
 
     // and add to covariance matrix
     covariance[0] += transformation(0,0) + area * (2*x0*xav0 + x0*x0);
@@ -294,7 +320,7 @@ assemble_covariance_matrix_3(InputIterator first,
   covariance[3] += mass * (-1.0 * c.y() * c.y());
   covariance[4] += mass * (-1.0 * c.z() * c.y());
   covariance[5] += mass * (-1.0 * c.z() * c.z());
-
+#endif
 }
 
 // assemble covariance matrix from a sphere set
diff --git a/Principal_component_analysis/test/Principal_component_analysis/test_validity.cpp b/Principal_component_analysis/test/Principal_component_analysis/test_validity.cpp
index 26bddf76fe7..c9a4e42ee61 100644
--- a/Principal_component_analysis/test/Principal_component_analysis/test_validity.cpp
+++ b/Principal_component_analysis/test/Principal_component_analysis/test_validity.cpp
@@ -225,7 +225,7 @@ void assert_quality (const std::vector<Point_3>& points, const Fitted& fitted)
   std::cerr << "mean distance = " << mean_dist << std::endl;
 
   CGAL_assertion_code
-    (double limit = 1e-3 * std::sqrt (CGAL::squared_distance (points.front(), points.back())));
+    (double limit = 1e-5 * std::sqrt (CGAL::squared_distance (points.front(), points.back())));
   CGAL_assertion (mean_dist < limit);
 }
 

From 6f19f920b059eb912c44a78173dba6a730937454 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Wed, 4 Nov 2020 11:00:21 +0100
Subject: [PATCH 12/19] Add missing include

---
 Solver_interface/include/CGAL/Eigen_diagonalize_traits.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Solver_interface/include/CGAL/Eigen_diagonalize_traits.h b/Solver_interface/include/CGAL/Eigen_diagonalize_traits.h
index eb2bb042b2f..03ebf8fe353 100644
--- a/Solver_interface/include/CGAL/Eigen_diagonalize_traits.h
+++ b/Solver_interface/include/CGAL/Eigen_diagonalize_traits.h
@@ -26,7 +26,7 @@
 
 //#define DO_NOT_USE_EIGEN_COMPUTEDIRECT_FOR_DIAGONALIZATION
 
-#include <CGAL/array.h>
+#include <CGAL/number_utils.h>
 
 namespace CGAL {
 

From ce76bf08f7e933728633736c51f94d83de0c9c9f Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Wed, 3 Mar 2021 11:43:46 +0100
Subject: [PATCH 13/19] Add missing include

---
 Solver_interface/include/CGAL/Eigen_diagonalize_traits.h | 1 +
 1 file changed, 1 insertion(+)

diff --git a/Solver_interface/include/CGAL/Eigen_diagonalize_traits.h b/Solver_interface/include/CGAL/Eigen_diagonalize_traits.h
index 03ebf8fe353..38ee3c7c4c7 100644
--- a/Solver_interface/include/CGAL/Eigen_diagonalize_traits.h
+++ b/Solver_interface/include/CGAL/Eigen_diagonalize_traits.h
@@ -27,6 +27,7 @@
 //#define DO_NOT_USE_EIGEN_COMPUTEDIRECT_FOR_DIAGONALIZATION
 
 #include <CGAL/number_utils.h>
+#include <CGAL/double.h>
 
 namespace CGAL {
 

From a4ac6061cfe8cce378ecbf13bc52293329cf4915 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Thu, 25 Mar 2021 10:40:59 +0100
Subject: [PATCH 14/19] Fix conversion warnings

---
 .../include/CGAL/linear_least_squares_fitting_rectangles_2.h  | 4 ++--
 .../include/CGAL/linear_least_squares_fitting_triangles_2.h   | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h
index 579287ce003..8c37f5960bb 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h
@@ -68,8 +68,8 @@ linear_least_squares_fitting_2(InputIterator first,
   typename DiagonalizeTraits::Covariance_matrix covariance = {{ 0., 0., 0. }};
 
   // assemble 2nd order moment about the origin.
-  FT temp[4] = {1/3.0, 0.25,
-                0.25,  1/3.0};
+  FT temp[4] = {FT(1/3.0), FT(0.25),
+                FT(0.25),  FT(1/3.0)};
   Matrix moment = init_matrix<FT>(2,temp);
 
   for(InputIterator it = first;
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
index 3cd5d651b24..7d32f2bcdbd 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
@@ -68,8 +68,8 @@ linear_least_squares_fitting_2(InputIterator first,
   typename DiagonalizeTraits::Covariance_matrix covariance = {{ 0., 0., 0. }};
 
   // assemble the 2nd order moment about the origin.
-  FT temp[4] = {1/12.0, 1/24.0,
-                1/24.0, 1/12.0};
+  FT temp[4] = {FT(1/12.0), FT(1/24.0),
+                FT(1/24.0), FT(1/12.0)};
 
   Matrix moment = init_matrix<FT>(2,temp);
 

From 3efe2ec339de85014ac766dc6b36fe825adeeff3 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Wed, 31 Mar 2021 09:40:19 +0200
Subject: [PATCH 15/19] Fix warnings: CGAL objects use int parameters for
 operator[]

---
 .../include/CGAL/PCA_util_Eigen.h                    |  4 ++--
 .../include/CGAL/Subiterator.h                       |  4 ++--
 .../CGAL/linear_least_squares_fitting_cuboids_3.h    |  8 ++++----
 .../CGAL/linear_least_squares_fitting_rectangles_2.h |  4 ++--
 .../CGAL/linear_least_squares_fitting_segments_2.h   |  2 +-
 .../CGAL/linear_least_squares_fitting_segments_3.h   |  4 ++--
 .../CGAL/linear_least_squares_fitting_tetrahedra_3.h | 12 ++++++------
 .../CGAL/linear_least_squares_fitting_triangles_2.h  |  4 ++--
 .../CGAL/linear_least_squares_fitting_triangles_3.h  |  8 ++++----
 9 files changed, 25 insertions(+), 25 deletions(-)

diff --git a/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h b/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h
index af8e8fb634e..5cd1dbbd7d8 100644
--- a/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h
+++ b/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h
@@ -220,10 +220,10 @@ assemble_covariance_matrix_3(InputIterator first,
   typedef typename K::FT              FT;
   typedef typename K::Iso_cuboid_3    Iso_cuboid;
   typedef typename K::Triangle_3      Triangle;
-  auto converter = [](const Iso_cuboid& c, std::size_t idx) -> Triangle
+  auto converter = [](const Iso_cuboid& c, int idx) -> Triangle
     {
       // Decomposition of 6 faces of the cuboid into 12 triangles
-      static constexpr std::array<std::array<std::size_t, 3>, 12 > indices
+      static constexpr std::array<std::array<int, 3>, 12 > indices
       = {{ { 0, 1, 2 }, { 0, 2, 3 }, { 2, 3, 4 }, { 2, 4, 7 },
            { 3, 4, 5 }, { 3, 5, 0 }, { 4, 5, 6 }, { 4, 6, 7 },
            { 5, 6, 1 }, { 5, 1, 0 }, { 6, 7, 2 }, { 6, 2, 1 } }};
diff --git a/Principal_component_analysis/include/CGAL/Subiterator.h b/Principal_component_analysis/include/CGAL/Subiterator.h
index fbadd1980c2..186776918be 100644
--- a/Principal_component_analysis/include/CGAL/Subiterator.h
+++ b/Principal_component_analysis/include/CGAL/Subiterator.h
@@ -31,12 +31,12 @@ public:
   using Input_type = typename std::iterator_traits<InputIterator>::value_type;
   using Output_type = ValueType;
 
-  using Converter = std::function<Output_type(const Input_type&, std::size_t)>;
+  using Converter = std::function<Output_type(const Input_type&, int)>;
 
 private:
 
   Converter m_converter;
-  std::size_t m_next_index;
+  int m_next_index;
   InputIterator m_base;
   mutable Output_type m_current;
 
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_cuboids_3.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_cuboids_3.h
index 4b552111dd2..c0071304d6a 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_cuboids_3.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_cuboids_3.h
@@ -113,7 +113,7 @@ linear_least_squares_fitting_3(InputIterator first,
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  auto converter = [](const Iso_cuboid& c, std::size_t idx) -> Segment
+  auto converter = [](const Iso_cuboid& c, int idx) -> Segment
     {
       if (idx < 7)
         return Segment (c[idx], c[idx+1]);
@@ -148,7 +148,7 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Point_3 Point;
   typedef typename K::Iso_cuboid_3 Iso_cuboid;
-  auto converter = [](const Iso_cuboid& c, std::size_t idx) -> Point { return c[idx]; };
+  auto converter = [](const Iso_cuboid& c, int idx) -> Point { return c[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
@@ -240,7 +240,7 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Segment_3 Segment;
   typedef typename K::Iso_cuboid_3 Iso_cuboid;
-  auto converter = [](const Iso_cuboid& c, std::size_t idx) -> Segment
+  auto converter = [](const Iso_cuboid& c, int idx) -> Segment
     {
       if (idx < 7)
         return Segment (c[idx], c[idx+1]);
@@ -279,7 +279,7 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Point_3 Point;
   typedef typename K::Iso_cuboid_3 Iso_cuboid;
-  auto converter = [](const Iso_cuboid& c, std::size_t idx) -> Point { return c[idx]; };
+  auto converter = [](const Iso_cuboid& c, int idx) -> Point { return c[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h
index 8c37f5960bb..34a00b79519 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h
@@ -158,7 +158,7 @@ linear_least_squares_fitting_2(InputIterator first,
   // types
   typedef typename K::Iso_rectangle_2 Iso_rectangle;
   typedef typename K::Segment_2         Segment_2;
-  auto converter = [](const Iso_rectangle& r, std::size_t idx) -> Segment_2 { return Segment_2(r[idx], r[(idx+1)%4]); };
+  auto converter = [](const Iso_rectangle& r, int idx) -> Segment_2 { return Segment_2(r[idx], r[(idx+1)%4]); };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
@@ -188,7 +188,7 @@ linear_least_squares_fitting_2(InputIterator first,
   // types
   typedef typename K::Iso_rectangle_2 Iso_rectangle;
   typedef typename K::Point_2         Point_2;
-  auto converter = [](const Iso_rectangle& r, std::size_t idx) -> Point_2 { return r[idx]; };
+  auto converter = [](const Iso_rectangle& r, int idx) -> Point_2 { return r[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_2.h
index aa3a3432fdf..5407c3fe497 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_2.h
@@ -147,7 +147,7 @@ linear_least_squares_fitting_2(InputIterator first,
   // types
   typedef typename K::Point_2  Point;
   typedef typename K::Segment_2 Segment;
-  auto converter = [](const Segment& s, std::size_t idx) -> Point { return s[idx]; };
+  auto converter = [](const Segment& s, int idx) -> Point { return s[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_3.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_3.h
index 7fababd309b..b685d8d1f41 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_3.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_3.h
@@ -75,7 +75,7 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Segment_3  Segment;
   typedef typename K::Point_3  Point;
-  auto converter = [](const Segment& s, std::size_t idx) -> Point { return s[idx]; };
+  auto converter = [](const Segment& s, int idx) -> Point { return s[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
@@ -134,7 +134,7 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Segment_3  Segment;
   typedef typename K::Point_3  Point;
-  auto converter = [](const Segment& s, std::size_t idx) -> Point { return s[idx]; };
+  auto converter = [](const Segment& s, int idx) -> Point { return s[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_tetrahedra_3.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_tetrahedra_3.h
index b7f36c518d2..cfe029d5f0c 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_tetrahedra_3.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_tetrahedra_3.h
@@ -76,7 +76,7 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Tetrahedron_3    Tetrahedron;
   typedef typename K::Triangle_3 Triangle;
-  auto converter = [](const Tetrahedron& t, std::size_t idx) -> Triangle
+  auto converter = [](const Tetrahedron& t, int idx) -> Triangle
     { return Triangle(t[idx], t[(idx+1)%4], t[(idx+2)%4]); };
 
   // precondition: at least one element in the container.
@@ -105,7 +105,7 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Tetrahedron_3    Tetrahedron;
   typedef typename K::Segment_3 Segment;
-  auto converter = [](const Tetrahedron& t, std::size_t idx) -> Segment
+  auto converter = [](const Tetrahedron& t, int idx) -> Segment
     {
       if (idx < 4)
         return Segment (t[idx], t[(idx+1)%4]);
@@ -139,7 +139,7 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Tetrahedron_3    Tetrahedron;
   typedef typename K::Point_3 Point;
-  auto converter = [](const Tetrahedron& t, std::size_t idx) -> Point { return t[idx]; };
+  auto converter = [](const Tetrahedron& t, int idx) -> Point { return t[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
@@ -198,7 +198,7 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Tetrahedron_3    Tetrahedron;
   typedef typename K::Triangle_3 Triangle;
-  auto converter = [](const Tetrahedron& t, std::size_t idx) -> Triangle
+  auto converter = [](const Tetrahedron& t, int idx) -> Triangle
     { return Triangle(t[idx], t[(idx+1)%4], t[(idx+2)%4]); };
 
   // precondition: at least one element in the container.
@@ -227,7 +227,7 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Tetrahedron_3    Tetrahedron;
   typedef typename K::Segment_3 Segment;
-  auto converter = [](const Tetrahedron& t, std::size_t idx) -> Segment
+  auto converter = [](const Tetrahedron& t, int idx) -> Segment
     {
       if (idx < 4)
         return Segment (t[idx], t[(idx+1)%4]);
@@ -260,7 +260,7 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Tetrahedron_3    Tetrahedron;
   typedef typename K::Point_3 Point;
-  auto converter = [](const Tetrahedron& t, std::size_t idx) -> Point { return t[idx]; };
+  auto converter = [](const Tetrahedron& t, int idx) -> Point { return t[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
index 7d32f2bcdbd..91c5d2b6449 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
@@ -159,7 +159,7 @@ linear_least_squares_fitting_2(InputIterator first,
   // types
   typedef typename Kernel::Triangle_2 Triangle;
   typedef typename Kernel::Segment_2  Segment;
-  auto converter = [](const Triangle& t, std::size_t idx) -> Segment { return Segment(t[idx], t[(idx+1)%3]); };
+  auto converter = [](const Triangle& t, int idx) -> Segment { return Segment(t[idx], t[(idx+1)%3]); };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
@@ -189,7 +189,7 @@ linear_least_squares_fitting_2(InputIterator first,
 
   typedef typename Kernel::Triangle_2 Triangle;
   typedef typename Kernel::Point_2 Point;
-  auto converter = [](const Triangle& t, std::size_t idx) -> Point { return t[idx]; };
+  auto converter = [](const Triangle& t, int idx) -> Point { return t[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_3.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_3.h
index 3c403399cbb..a2a1e96bec4 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_3.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_3.h
@@ -74,7 +74,7 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Triangle_3  Triangle;
   typedef typename K::Segment_3  Segment;
-  auto converter = [](const Triangle& t, std::size_t idx) -> Segment { return Segment(t[idx], t[(idx+1)%3]); };
+  auto converter = [](const Triangle& t, int idx) -> Segment { return Segment(t[idx], t[(idx+1)%3]); };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
@@ -103,7 +103,7 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Triangle_3  Triangle;
   typedef typename K::Point_3  Point;
-  auto converter = [](const Triangle& t, std::size_t idx) -> Point { return t[idx]; };
+  auto converter = [](const Triangle& t, int idx) -> Point { return t[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
@@ -163,7 +163,7 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Triangle_3  Triangle;
   typedef typename K::Segment_3  Segment;
-  auto converter = [](const Triangle& t, std::size_t idx) -> Segment { return Segment(t[idx], t[(idx+1)%3]); };
+  auto converter = [](const Triangle& t, int idx) -> Segment { return Segment(t[idx], t[(idx+1)%3]); };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
@@ -191,7 +191,7 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Triangle_3  Triangle;
   typedef typename K::Point_3  Point;
-  auto converter = [](const Triangle& t, std::size_t idx) -> Point { return t[idx]; };
+  auto converter = [](const Triangle& t, int idx) -> Point { return t[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);

From 53c24ef300943d9d985ae84ec0380b2e5624f9dd Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Tue, 6 Apr 2021 08:30:06 +0200
Subject: [PATCH 16/19] Fix conversion warnings in test kernels

---
 .../test_kernels.cpp                          | 31 ++++++++++++-------
 1 file changed, 19 insertions(+), 12 deletions(-)

diff --git a/Principal_component_analysis/test/Principal_component_analysis/test_kernels.cpp b/Principal_component_analysis/test/Principal_component_analysis/test_kernels.cpp
index b96b98e298a..606f3abbc56 100644
--- a/Principal_component_analysis/test/Principal_component_analysis/test_kernels.cpp
+++ b/Principal_component_analysis/test/Principal_component_analysis/test_kernels.cpp
@@ -5,31 +5,38 @@
 #include <CGAL/linear_least_squares_fitting_2.h>
 #include <CGAL/linear_least_squares_fitting_3.h>
 
+template <typename K>
+typename K::Point_2 point_2 (int x, int y)
+{ return typename K::Point_2 (typename K::FT(x), typename K::FT(y)); }
+template <typename K>
+typename K::Point_3 point_3 (int x, int y, int z)
+{ return typename K::Point_3 (typename K::FT(x), typename K::FT(y), typename K::FT(z)); }
+
 // Generate default objects so that the test compiles *AND* runs fine
 template <typename K, typename O> O default_object(const O&) { return O(); }
 template <typename K> typename K::Point_2 default_object(const typename K::Point_2&)
-{ return typename K::Point_2(0,0); }
+{ return point_2<K>(0,0); }
 template <typename K> typename K::Segment_2 default_object(const typename K::Segment_2&)
-{ return typename K::Segment_2(typename K::Point_2(0,0), typename K::Point_2(0,1)); }
+{ return typename K::Segment_2(point_2<K>(0,0), point_2<K>(0,1)); }
 template <typename K> typename K::Circle_2 default_object(const typename K::Circle_2&)
-{ return typename K::Circle_2(typename K::Point_2(0,0), typename K::FT(1.0)); }
+{ return typename K::Circle_2(point_2<K>(0,0), typename K::FT(1.0)); }
 template <typename K> typename K::Triangle_2 default_object(const typename K::Triangle_2&)
-{ return typename K::Triangle_2(typename K::Point_2(0,0), typename K::Point_2(0,1), typename K::Point_2(1,0)); }
+{ return typename K::Triangle_2(point_2<K>(0,0), point_2<K>(0,1), point_2<K>(1,0)); }
 template <typename K> typename K::Iso_rectangle_2 default_object(const typename K::Iso_rectangle_2&)
-{ return typename K::Iso_rectangle_2(typename K::Point_2(0,0), typename K::Point_2(1,1)); }
+{ return typename K::Iso_rectangle_2(point_2<K>(0,0), point_2<K>(1,1)); }
 template <typename K> typename K::Point_3 default_object(const typename K::Point_3&)
-{ return typename K::Point_3(0,0,0); }
+{ return point_3<K>(0,0,0); }
 template <typename K> typename K::Segment_3 default_object(const typename K::Segment_3&)
-{ return typename K::Segment_3(typename K::Point_3(0,0,0), typename K::Point_3(0,0,1)); }
+{ return typename K::Segment_3(point_3<K>(0,0,0), point_3<K>(0,0,1)); }
 template <typename K> typename K::Sphere_3 default_object(const typename K::Sphere_3&)
-{ return typename K::Sphere_3(typename K::Point_3(0,0,0), typename K::FT(1.0)); }
+{ return typename K::Sphere_3(point_3<K>(0,0,0), typename K::FT(1.0)); }
 template <typename K> typename K::Triangle_3 default_object(const typename K::Triangle_3&)
-{ return typename K::Triangle_3(typename K::Point_3(0,0,0), typename K::Point_3(0,0,1), typename K::Point_3(0,1,0)); }
+{ return typename K::Triangle_3(point_3<K>(0,0,0), point_3<K>(0,0,1), point_3<K>(0,1,0)); }
 template <typename K> typename K::Tetrahedron_3 default_object(const typename K::Tetrahedron_3&)
-{ return typename K::Tetrahedron_3(typename K::Point_3(0,0,0), typename K::Point_3(0,0,1),
-                                   typename K::Point_3(0,1,0), typename K::Point_3(1,0,0)); }
+{ return typename K::Tetrahedron_3(point_3<K>(0,0,0), point_3<K>(0,0,1),
+                                   point_3<K>(0,1,0), point_3<K>(1,0,0)); }
 template <typename K> typename K::Iso_cuboid_3 default_object(const typename K::Iso_cuboid_3&)
-{ return typename K::Iso_cuboid_3(typename K::Point_3(0,0,0), typename K::Point_3(1,1,1)); }
+{ return typename K::Iso_cuboid_3(point_3<K>(0,0,0), point_3<K>(1,1,1)); }
 
 template <typename Kernel, typename Object, int dim>
 void test_2d()

From 7d716c0be7c6728494722b9f603720474d44d9b2 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Tue, 6 Apr 2021 08:38:56 +0200
Subject: [PATCH 17/19] Fix conversion warnings

---
 .../include/CGAL/PCA_util_Eigen.h             | 96 +++++++++----------
 .../linear_least_squares_fitting_circles_2.h  | 12 +--
 ...inear_least_squares_fitting_rectangles_2.h |  6 +-
 .../linear_least_squares_fitting_segments_2.h |  8 +-
 ...linear_least_squares_fitting_triangles_2.h |  6 +-
 5 files changed, 64 insertions(+), 64 deletions(-)

diff --git a/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h b/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h
index 5cd1dbbd7d8..e7fd18b8842 100644
--- a/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h
+++ b/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h
@@ -56,9 +56,9 @@ assemble_covariance_matrix_3(InputIterator first,
 
   // assemble 2nd order moment about the origin.
   Matrix moment;
-  moment << 1.0/12.0, 1.0/24.0, 1.0/24.0,
-            1.0/24.0, 1.0/12.0, 1.0/24.0,
-            1.0/24.0, 1.0/24.0, 1.0/12.0;
+  moment << FT(1.0/12.0), FT(1.0/24.0), FT(1.0/24.0),
+            FT(1.0/24.0), FT(1.0/12.0), FT(1.0/24.0),
+            FT(1.0/24.0), FT(1.0/24.0), FT(1.0/12.0);
 
   for(InputIterator it = first;
       it != beyond;
@@ -100,12 +100,12 @@ assemble_covariance_matrix_3(InputIterator first,
 
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
-  covariance[0] += mass * (-1.0 * c.x() * c.x());
-  covariance[1] += mass * (-1.0 * c.x() * c.y());
-  covariance[2] += mass * (-1.0 * c.z() * c.x());
-  covariance[3] += mass * (-1.0 * c.y() * c.y());
-  covariance[4] += mass * (-1.0 * c.z() * c.y());
-  covariance[5] += mass * (-1.0 * c.z() * c.z());
+  covariance[0] += -mass * (c.x() * c.x());
+  covariance[1] += -mass * (c.x() * c.y());
+  covariance[2] += -mass * (c.z() * c.x());
+  covariance[3] += -mass * (c.y() * c.y());
+  covariance[4] += -mass * (c.z() * c.y());
+  covariance[5] += -mass * (c.z() * c.z());
 
 }
 
@@ -176,9 +176,9 @@ assemble_covariance_matrix_3(InputIterator first,
     transformation = volume * transformation * moment * transformation.transpose();
 
     // Translate the 2nd order moment to the minimum corner (x0,y0,z0) of the cuboid.
-    FT xav0 = (x1 - x0) / (2.0);
-    FT yav0 = (y1 - y0) / (2.0);
-    FT zav0 = (z1 - z0) / (2.0);
+    FT xav0 = (x1 - x0) / FT(2.0);
+    FT yav0 = (y1 - y0) / FT(2.0);
+    FT zav0 = (z1 - z0) / FT(2.0);
 
     // and add to covariance matrix
     covariance[0] += transformation(0,0) + volume * (2*x0*xav0 + x0*x0);
@@ -314,12 +314,12 @@ assemble_covariance_matrix_3(InputIterator first,
 
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
-  covariance[0] += mass * (-1.0 * c.x() * c.x());
-  covariance[1] += mass * (-1.0 * c.x() * c.y());
-  covariance[2] += mass * (-1.0 * c.z() * c.x());
-  covariance[3] += mass * (-1.0 * c.y() * c.y());
-  covariance[4] += mass * (-1.0 * c.z() * c.y());
-  covariance[5] += mass * (-1.0 * c.z() * c.z());
+  covariance[0] += -mass * (c.x() * c.x());
+  covariance[1] += -mass * (c.x() * c.y());
+  covariance[2] += -mass * (c.z() * c.x());
+  covariance[3] += -mass * (c.y() * c.y());
+  covariance[4] += -mass * (c.z() * c.y());
+  covariance[5] += -mass * (c.z() * c.z());
 #endif
 }
 
@@ -350,9 +350,9 @@ assemble_covariance_matrix_3(InputIterator first,
 
   // assemble 2nd order moment about the origin.
   Matrix moment;
-  moment << 4.0/15.0, 0.0,      0.0,
-            0.0,      4.0/15.0, 0.0,
-            0.0,      0.0,      4.0/15.0;
+  moment << FT(4.0/15.0), FT(0.0),      FT(0.0),
+            FT(0.0),      FT(4.0/15.0), FT(0.0),
+            FT(0.0),      FT(0.0),      FT(4.0/15.0);
 
   for(InputIterator it = first;
       it != beyond;
@@ -399,12 +399,12 @@ assemble_covariance_matrix_3(InputIterator first,
 
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
-  covariance[0] += mass * (-1.0 * c.x() * c.x());
-  covariance[1] += mass * (-1.0 * c.x() * c.y());
-  covariance[2] += mass * (-1.0 * c.z() * c.x());
-  covariance[3] += mass * (-1.0 * c.y() * c.y());
-  covariance[4] += mass * (-1.0 * c.z() * c.y());
-  covariance[5] += mass * (-1.0 * c.z() * c.z());
+  covariance[0] += -mass * (c.x() * c.x());
+  covariance[1] += -mass * (c.x() * c.y());
+  covariance[2] += -mass * (c.z() * c.x());
+  covariance[3] += -mass * (c.y() * c.y());
+  covariance[4] += -mass * (c.z() * c.y());
+  covariance[5] += -mass * (c.z() * c.z());
 
 }
 // assemble covariance matrix from a sphere set
@@ -434,9 +434,9 @@ assemble_covariance_matrix_3(InputIterator first,
 
   // assemble 2nd order moment about the origin.
   Matrix moment;
-  moment << 4.0/3.0, 0.0,     0.0,
-            0.0,     4.0/3.0, 0.0,
-            0.0,     0.0,     4.0/3.0;
+  moment << FT(4.0/3.0), FT(0.0),     FT(0.0),
+            FT(0.0),     FT(4.0/3.0), FT(0.0),
+            FT(0.0),     FT(0.0),     FT(4.0/3.0);
 
   for(InputIterator it = first;
       it != beyond;
@@ -484,12 +484,12 @@ assemble_covariance_matrix_3(InputIterator first,
 
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
-  covariance[0] += mass * (-1.0 * c.x() * c.x());
-  covariance[1] += mass * (-1.0 * c.x() * c.y());
-  covariance[2] += mass * (-1.0 * c.z() * c.x());
-  covariance[3] += mass * (-1.0 * c.y() * c.y());
-  covariance[4] += mass * (-1.0 * c.z() * c.y());
-  covariance[5] += mass * (-1.0 * c.z() * c.z());
+  covariance[0] += -mass * (c.x() * c.x());
+  covariance[1] += -mass * (c.x() * c.y());
+  covariance[2] += -mass * (c.z() * c.x());
+  covariance[3] += -mass * (c.y() * c.y());
+  covariance[4] += -mass * (c.z() * c.y());
+  covariance[5] += -mass * (c.z() * c.z());
 
 }
 
@@ -520,9 +520,9 @@ assemble_covariance_matrix_3(InputIterator first,
   //     5
   // assemble 2nd order moment about the origin.
   Matrix moment;
-  moment << 1.0/60.0,  1.0/120.0, 1.0/120.0,
-            1.0/120.0, 1.0/60.0,  1.0/120.0,
-            1.0/120.0, 1.0/120.0, 1.0/60.0;
+  moment << FT(1.0/60.0),  FT(1.0/120.0), FT(1.0/120.0),
+            FT(1.0/120.0), FT(1.0/60.0),  FT(1.0/120.0),
+            FT(1.0/120.0), FT(1.0/120.0), FT(1.0/60.0);
 
   Matrix accum; // zero by default
   accum << 0, 0, 0, 0, 0, 0, 0, 0, 0;
@@ -597,9 +597,9 @@ assemble_covariance_matrix_3(InputIterator first,
 
   // assemble 2nd order moment about the origin.
   Matrix moment;
-  moment << 1.0/3.0, 0.5/3.0, 0.0,
-            0.5/3.0, 1.0/3.0, 0.0,
-            0.0,     0.0,     0.0;
+  moment << FT(1.0/3.0), FT(0.5/3.0), FT(0.0),
+            FT(0.5/3.0), FT(1.0/3.0), FT(0.0),
+            FT(0.0),     FT(0.0),     FT(0.0);
 
   for(InputIterator it = first;
       it != beyond;
@@ -641,12 +641,12 @@ assemble_covariance_matrix_3(InputIterator first,
 
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
-  covariance[0] += mass * (-1.0 * c.x() * c.x());
-  covariance[1] += mass * (-1.0 * c.x() * c.y());
-  covariance[2] += mass * (-1.0 * c.z() * c.x());
-  covariance[3] += mass * (-1.0 * c.y() * c.y());
-  covariance[4] += mass * (-1.0 * c.z() * c.y());
-  covariance[5] += mass * (-1.0 * c.z() * c.z());
+  covariance[0] += -mass * (c.x() * c.x());
+  covariance[1] += -mass * (c.x() * c.y());
+  covariance[2] += -mass * (c.z() * c.x());
+  covariance[3] += -mass * (c.y() * c.y());
+  covariance[4] += -mass * (c.z() * c.y());
+  covariance[5] += -mass * (c.z() * c.z());
 
 }
 
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_circles_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_circles_2.h
index 3b4f679ce3b..715f76cd0d0 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_circles_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_circles_2.h
@@ -111,9 +111,9 @@ linear_least_squares_fitting_2(InputIterator first,
 
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
-  covariance[0] += mass * (-1.0 * c.x() * c.x());
-  covariance[1] += mass * (-1.0 * c.x() * c.y());
-  covariance[2] += mass * (-1.0 * c.y() * c.y());
+  covariance[0] += -mass * (c.x() * c.x());
+  covariance[1] += -mass * (c.x() * c.y());
+  covariance[2] += -mass * (c.y() * c.y());
 
   // solve for eigenvalues and eigenvectors.
   // eigen values are sorted in ascending order,
@@ -217,9 +217,9 @@ linear_least_squares_fitting_2(InputIterator first,
 
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
-  covariance[0] += mass * (-1.0 * c.x() * c.x());
-  covariance[1] += mass * (-1.0 * c.x() * c.y());
-  covariance[2] += mass * (-1.0 * c.y() * c.y());
+  covariance[0] += -mass * (c.x() * c.x());
+  covariance[1] += -mass * (c.x() * c.y());
+  covariance[2] += -mass * (c.y() * c.y());
 
   // solve for eigenvalues and eigenvectors.
   // eigen values are sorted in ascending order,
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h
index 34a00b79519..c3d8c2602f7 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h
@@ -115,9 +115,9 @@ linear_least_squares_fitting_2(InputIterator first,
 
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
-  covariance[0] += mass * (-1.0 * c.x() * c.x());
-  covariance[1] += mass * (-1.0 * c.x() * c.y());
-  covariance[2] += mass * (-1.0 * c.y() * c.y());
+  covariance[0] += -mass * (c.x() * c.x());
+  covariance[1] += -mass * (c.x() * c.y());
+  covariance[2] += -mass * (c.y() * c.y());
 
   // solve for eigenvalues and eigenvectors.
   // eigen values are sorted in ascending order,
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_2.h
index 5407c3fe497..a6d581d97df 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_2.h
@@ -104,9 +104,9 @@ linear_least_squares_fitting_2(InputIterator first,
 
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
-  covariance[0] += mass * (-1.0 * c.x() * c.x());
-  covariance[1] += mass * (-1.0 * c.x() * c.y());
-  covariance[2] += mass * (-1.0 * c.y() * c.y());
+  covariance[0] += -mass * ( c.x() * c.x());
+  covariance[1] += -mass * (c.x() * c.y());
+  covariance[2] += -mass * (c.y() * c.y());
 
   // solve for eigenvalues and eigenvectors.
   // eigen values are sorted in ascending order,
@@ -128,7 +128,7 @@ linear_least_squares_fitting_2(InputIterator first,
     // isotropic case (infinite number of directions)
     // by default: assemble a line that goes through
     // the centroid and with a default horizontal vector.
-    line = Line(c, Vector(1.0, 0.0));
+    line = Line(c, Vector(FT(1), FT(0)));
     return (FT)0.0;
   }
 } // end linear_least_squares_fitting_2 for segment set with 1D tag
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
index 91c5d2b6449..15843aa40fb 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
@@ -112,9 +112,9 @@ linear_least_squares_fitting_2(InputIterator first,
 
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
-  covariance[0] += mass * (-1.0 * c.x() * c.x());
-  covariance[1] += mass * (-1.0 * c.x() * c.y());
-  covariance[2] += mass * (-1.0 * c.y() * c.y());
+  covariance[0] += -mass * (c.x() * c.x());
+  covariance[1] += -mass * (c.x() * c.y());
+  covariance[2] += -mass * (c.y() * c.y());
 
   //  std::cout<<"cov: "<<covariance[0]*covariance[2]<<" =? "<<covariance[1]*covariance[1]<<std::endl;
 

From c9a331ec192e37b4c949768d0f738f48ae7d3497 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Tue, 13 Apr 2021 10:18:06 +0200
Subject: [PATCH 18/19] Fix more conversion warnings

---
 .../include/CGAL/linear_least_squares_fitting_circles_2.h   | 2 +-
 .../CGAL/linear_least_squares_fitting_rectangles_2.h        | 4 ++--
 .../include/CGAL/linear_least_squares_fitting_triangles_2.h | 6 +++---
 3 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_circles_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_circles_2.h
index 715f76cd0d0..c0698acc85f 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_circles_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_circles_2.h
@@ -241,7 +241,7 @@ linear_least_squares_fitting_2(InputIterator first,
     // isotropic case (infinite number of directions)
     // by default: assemble a line that goes through
     // the centroid and with a default horizontal vector.
-    line = Line(c, Vector(1.0, 0.0));
+    line = Line(c, Vector(FT(1), FT(0)));
     return (FT)0.0;
   }
 } // end linear_least_squares_fitting_2 for circle set with 1D tag
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h
index c3d8c2602f7..7fec0fea638 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h
@@ -101,8 +101,8 @@ linear_least_squares_fitting_2(InputIterator first,
     transformation = area * transformation * moment * LA::transpose(transformation);
 
     // Translate the 2nd order moment to the center of the rectangle.
-    FT xav0 = (x1-x0)/2.0;
-    FT yav0 = (y2-y0)/2.0;
+    FT xav0 = (x1-x0)/FT(2);
+    FT yav0 = (y2-y0)/FT(2);
     // and add to covariance matrix
     covariance[0] += transformation[0][0] + area * (x0*xav0*2 + x0*x0);
     covariance[1] += transformation[0][1] + area * (x0*yav0 + xav0*y0 + x0*y0);
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
index 15843aa40fb..61857a079e8 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
@@ -88,7 +88,7 @@ linear_least_squares_fitting_2(InputIterator first,
                    t[1].y() - y0, t[2].y() - y0};
 
     Matrix transformation = init_matrix<FT>(2,delta);
-    FT area = 0.5 * CGAL::abs(LA::determinant(transformation));
+    FT area = FT(0.5) * CGAL::abs(LA::determinant(transformation));
     CGAL_assertion(area!=0);
 
     // Find the 2nd order moment for the triangle wrt to the origin by an affine transformation.
@@ -97,8 +97,8 @@ linear_least_squares_fitting_2(InputIterator first,
     transformation = 2 * area * transformation * moment * LA::transpose(transformation);
 
     // Translate the 2nd order moment to (x0,y0).
-    FT xav0 = (delta[0]+delta[1])/3.0;
-    FT yav0 = (delta[2]+delta[3])/3.0;
+    FT xav0 = (delta[0]+delta[1])/FT(3);
+    FT yav0 = (delta[2]+delta[3])/FT(3);
 
     // and add to the covariance matrix
     covariance[0] += transformation[0][0] + area * (x0*xav0*2 + x0*x0);

From fe3da4cc9fc4fc3bee6ca0c1c800c48b411a130d Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Wed, 14 Apr 2021 09:45:06 +0200
Subject: [PATCH 19/19] Fix conversion warnings AGAIN

---
 .../include/CGAL/linear_least_squares_fitting_circles_2.h       | 2 +-
 .../include/CGAL/linear_least_squares_fitting_rectangles_2.h    | 2 +-
 .../include/CGAL/linear_least_squares_fitting_triangles_2.h     | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_circles_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_circles_2.h
index c0698acc85f..d6abf3d1482 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_circles_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_circles_2.h
@@ -135,7 +135,7 @@ linear_least_squares_fitting_2(InputIterator first,
     // isotropic case (infinite number of directions)
     // by default: assemble a line that goes through
     // the centroid and with a default horizontal vector.
-    line = Line(c, Vector(1.0, 0.0));
+    line = Line(c, Vector(FT(1), FT(0)));
     return (FT)0.0;
   }
 
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h
index 7fec0fea638..3bec5e8d4d5 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h
@@ -139,7 +139,7 @@ linear_least_squares_fitting_2(InputIterator first,
     // isotropic case (infinite number of directions)
     // by default: assemble a line that goes through
     // the centroid and with a default horizontal vector.
-    line = Line(c, Vector(1.0, 0.0));
+    line = Line(c, Vector(FT(1), FT(0)));
     return (FT)0.0;
   }
 } // end linear_least_squares_fitting_2 for rectangle set with 2D tag
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
index 61857a079e8..ab96ad6d408 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
@@ -138,7 +138,7 @@ linear_least_squares_fitting_2(InputIterator first,
     // isotropic case (infinite number of directions)
     // by default: assemble a line that goes through
     // the centroid and with a default horizontal vector.
-    line = Line(c, Vector(1.0, 0.0));
+    line = Line(c, Vector(FT(1), FT(0)));
     return (FT)0.0;
   }
 } // end linear_least_squares_fitting_2 for triangle set with 2D tag