diff --git a/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h b/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h index 66494ba125a..dfabb50a0b0 100644 --- a/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h +++ b/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h @@ -144,9 +144,9 @@ assemble_covariance_matrix_3(InputIterator first, FT x0 = t[0].x(); FT y0 = t[0].y(); FT z0 = t[0].z(); - FT delta[9] = {t[0].x(), t[1].x(), t[2].x(), - t[0].y(), t[1].y(), t[2].y(), - t[0].z(), t[1].z(), t[2].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(); @@ -548,9 +548,9 @@ assemble_covariance_matrix_3(InputIterator first, // assemble 2nd order moment about the origin. Matrix moment; - moment << 1.0, 0.5, 0.0, - 0.5, 1.0, 0.0, - 0.0, 0.0, 0.0; + 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; for(InputIterator it = first; it != 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 67d23e04d64..7aad5b86554 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 @@ -171,7 +171,8 @@ linear_least_squares_fitting_2(InputIterator first, segments.push_back(Segment_2(t[3],t[0])); } - return linear_least_squares_fitting_2(segments.begin(),segments.end(),line,c,K(),tag, + return linear_least_squares_fitting_2(segments.begin(),segments.end(),line,c, + (Segment_2*)nullptr, K(),tag, diagonalize_traits); } // end linear_least_squares_fitting_2 for rectangle set with 1D tag @@ -209,7 +210,8 @@ linear_least_squares_fitting_2(InputIterator first, points.push_back(Point_2(t[3])); } - return linear_least_squares_fitting_2(points.begin(),points.end(),line,c,K(),tag, + return linear_least_squares_fitting_2(points.begin(),points.end(),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 8bf70443120..f290e741eef 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 @@ -157,7 +157,7 @@ linear_least_squares_fitting_2(InputIterator first, points.push_back(s[0]); points.push_back(s[1]); } - return linear_least_squares_fitting_2(points.begin(),points.end(),line,c,k,(Point*)nullptr,tag, + return linear_least_squares_fitting_2(points.begin(),points.end(),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/test/Principal_component_analysis/test_validity.cpp b/Principal_component_analysis/test/Principal_component_analysis/test_validity.cpp new file mode 100644 index 00000000000..26bddf76fe7 --- /dev/null +++ b/Principal_component_analysis/test/Principal_component_analysis/test_validity.cpp @@ -0,0 +1,378 @@ +/* + This test checks that the outputs of all variants of PCA are + correct. + + To do so: + - N Point_3 objects are generated on a random line + - N*N Point_3 objects are generated on a random plane + + For each variant of the PCA, we generate objects "centered" on the + reference points: + - Sphere_3 centered on each point with random radius + - Equilateral Triangle_3 centered on each point with random length + - Point_2 using first 2 coordinates of Point_3 + - etc. + + PCA is then computed for Line_2 in 2D or for Line_3 and Plane_3 in + 3D. Then, the mean distance between the original points (or 2D + variants if 2D case) and the estimated support is computed: if it + higher than 0.01% of the diameter of the point cloud, an assertion + is raised. + + This only tests an obvious case where the PCA is computed on a set + of exactly aligned objects, but it should help avoiding the + introduction of bugs. + */ + +#include +#include +#include +#include + +#include + +typedef CGAL::Simple_cartesian Kernel; +typedef Kernel::Point_3 Point_3; +typedef Kernel::Segment_3 Segment_3; +typedef Kernel::Vector_3 Vector_3; +typedef Kernel::Sphere_3 Sphere_3; +typedef Kernel::Triangle_3 Triangle_3; +typedef Kernel::Tetrahedron_3 Tetrahedron_3; +typedef Kernel::Iso_cuboid_3 Iso_cuboid_3; +typedef Kernel::Plane_3 Plane_3; +typedef Kernel::Line_3 Line_3; + +typedef Kernel::Point_2 Point_2; +typedef Kernel::Vector_2 Vector_2; +typedef Kernel::Circle_2 Circle_2; +typedef Kernel::Iso_rectangle_2 Iso_rectangle_2; +typedef Kernel::Segment_2 Segment_2; +typedef Kernel::Triangle_2 Triangle_2; +typedef Kernel::Line_2 Line_2; + +CGAL::Random rnd; + +/* + Functions used to generate objects "centered" on points (objects + whose centers of mass are on the points). + */ + +void generate_objects_centered_on_points (const std::vector& source, + std::vector& target) +{ + target.reserve (source.size()); + std::copy (source.begin(), source.end(), std::back_inserter (target)); +} + +void generate_objects_centered_on_points (const std::vector& source, + std::vector& target) +{ + target.reserve (source.size()); + for (std::size_t i = 0; i < source.size(); ++ i) + target.push_back (Point_2 (source[i].x(), source[i].y())); +} + +void generate_objects_centered_on_points (const std::vector& source, + std::vector& target) +{ + Vector_3 diff (rnd.get_double(0, 0.01), + rnd.get_double(0, 0.01), + rnd.get_double(0, 0.01)); + + target.reserve (source.size()); + for (std::size_t i = 0; i < source.size(); ++ i) + target.push_back (Segment_3 (source[i] + diff, source[i] - diff)); +} + +void generate_objects_centered_on_points (const std::vector& source, + std::vector& target) +{ + Vector_2 diff (rnd.get_double(0, 0.01), + rnd.get_double(0, 0.01)); + + target.reserve (source.size()); + for (std::size_t i = 0; i < source.size(); ++ i) + { + Point_2 p (source[i].x(), source[i].y()); + target.push_back (Segment_2 (p + diff, p - diff)); + } +} + +void generate_objects_centered_on_points (const std::vector& source, + std::vector& target) +{ + double radius = rnd.get_double(0, 0.01); + target.reserve (source.size()); + for (std::size_t i = 0; i < source.size(); ++ i) + target.push_back (Sphere_3 (source[i], radius * radius)); +} + +void generate_objects_centered_on_points (const std::vector& source, + std::vector& target) +{ + double radius = rnd.get_double(0, 0.01); + target.reserve (source.size()); + for (std::size_t i = 0; i < source.size(); ++ i) + { + Point_2 p (source[i].x(), source[i].y()); + target.push_back (Circle_2 (p, radius * radius)); + } +} + +void generate_objects_centered_on_points (const std::vector& source, + std::vector& target) +{ + // General equilateral triangles + double length = rnd.get_double(0, 0.01); + double height = length * std::sqrt(3.) / 2.; + + target.reserve (source.size()); + for (std::size_t i = 0; i < source.size(); ++ i) + target.push_back + (Triangle_3 (source[i] + Vector_3 (0, 0, 2. * height / 3.), + source[i] + Vector_3 (0, -length/2., -height/3.), + source[i] + Vector_3 (0, length/2., -height/3.))); +} + +void generate_objects_centered_on_points (const std::vector& source, + std::vector& target) +{ + // General equilateral triangles + double length = rnd.get_double(0, 0.01); + double height = length * std::sqrt(3.) / 2.; + + target.reserve (source.size()); + for (std::size_t i = 0; i < source.size(); ++ i) + { + Point_2 p (source[i].x(), source[i].y()); + target.push_back + (Triangle_2 (p + Vector_2 (0, 2. * height / 3.), + p + Vector_2 (-length/2., -height/3.), + p + Vector_2 (length/2., -height/3.))); + } +} + +void generate_objects_centered_on_points (const std::vector& source, + std::vector& target) +{ + // General regular tetrahedra + double length = rnd.get_double(0, 0.01); + double height = length * std::sqrt(2. / 3.); + + double height_tri = length * std::sqrt(3.) / 2.; + + target.reserve (source.size()); + for (std::size_t i = 0; i < source.size(); ++ i) + target.push_back + (Tetrahedron_3 (source[i] + Vector_3 (0, 0, 3. * height / 4.), + source[i] + Vector_3 (0, 2. * height_tri / 3., -height/4.), + source[i] + Vector_3 (-length/2., -height_tri / 3., -height/4.), + source[i] + Vector_3 (length/2., -height_tri / 3., -height/4.))); +} + +void generate_objects_centered_on_points (const std::vector& source, + std::vector& target) +{ + Vector_3 diff (rnd.get_double(0, 0.01), + rnd.get_double(0, 0.01), + rnd.get_double(0, 0.01)); + + target.reserve (source.size()); + for (std::size_t i = 0; i < source.size(); ++ i) + target.push_back (Iso_cuboid_3 (source[i] + diff, source[i] - diff)); +} + +void generate_objects_centered_on_points (const std::vector& source, + std::vector& target) +{ + Vector_2 diff (rnd.get_double(0, 0.01), + rnd.get_double(0, 0.01)); + + target.reserve (source.size()); + for (std::size_t i = 0; i < source.size(); ++ i) + { + Point_2 p (source[i].x(), source[i].y()); + target.push_back (Iso_rectangle_2 (p + diff, p - diff)); + } +} + +/* + Compute distances between original points to Line_3/Plane_3/Line_2 + */ + +double distance (const Point_3& p, const Line_3& l) +{ return std::sqrt (CGAL::squared_distance (p, l)); } +double distance (const Point_3& p, const Plane_3& pl) +{ return std::sqrt (CGAL::squared_distance (p, pl)); } +double distance (const Point_3& p, const Line_2& l) +{ return std::sqrt (CGAL::squared_distance (Point_2 (p.x(), p.y()), l)); } + +/* + Test quality of fit and raise assertion if too low. + */ + +template +void assert_quality (const std::vector& points, const Fitted& fitted) +{ + double mean_dist = 0; + for (std::size_t i = 0; i < points.size(); ++ i) + { + double dist = distance (points[i], fitted); + mean_dist += dist; + } + mean_dist /= points.size(); + + std::cerr << "mean distance = " << mean_dist << std::endl; + + CGAL_assertion_code + (double limit = 1e-3 * std::sqrt (CGAL::squared_distance (points.front(), points.back()))); + CGAL_assertion (mean_dist < limit); +} + +/* + Reference test functions, both in 2D and 3D. + */ + +template +void test (const std::vector& objects, + const std::vector& points, + const CGAL::Dimension_tag<2>& /* ambient dimension */) +{ + std::cerr << " CGAL::Dimension_tag<" << dim << ">: "; + Fitted fitted; + CGAL::linear_least_squares_fitting_2 + (objects.begin(), objects.end(), fitted, CGAL::Dimension_tag()); + assert_quality (points, fitted); +} + +template +void test (const std::vector& objects, + const std::vector& points, + const CGAL::Dimension_tag<3>& /* ambient dimension */) +{ + std::cerr << " CGAL::Dimension_tag<" << dim << ">: "; + Fitted fitted; + CGAL::linear_least_squares_fitting_3 + (objects.begin(), objects.end(), fitted, CGAL::Dimension_tag()); + assert_quality (points, fitted); +} + +template +void test (const std::vector& points) +{ + std::vector objects; + generate_objects_centered_on_points (points, objects); + test + (objects, points, typename CGAL::Ambient_dimension::type()); +} + + + + +int main() +{ + std::cerr << "Validity test with seed " << rnd.get_seed() << std::endl; + + Point_3 origin (rnd.get_double(), rnd.get_double(), rnd.get_double()); + Vector_3 base1 (rnd.get_double(), rnd.get_double(), rnd.get_double()); + Vector_3 base2 (rnd.get_double(), rnd.get_double(), rnd.get_double()); + + std::cerr << "Origin = " << origin << std::endl + << "Base1 = " << base1 << std::endl + << "Base2 = " << base2 << std::endl; + + std::size_t nb_points_on_line = 100; + std::vector points_on_line; + points_on_line.reserve (nb_points_on_line); + for (std::size_t i = 0; i < nb_points_on_line; ++ i) + points_on_line.push_back (Point_3 (origin + double(i) * base1)); + + std::vector points_on_plane; + points_on_plane.reserve (nb_points_on_line * nb_points_on_line); + for (std::size_t i = 0; i < nb_points_on_line; ++ i) + for (std::size_t j = 0; j < nb_points_on_line; ++ j) + points_on_plane.push_back (Point_3 (origin + double(i) * base1 + double(j) * base2)); + + std::cerr << std::endl << "=== 2D ===" << std::endl << std::endl; + + std::cerr << "[Testing line fitting on Point_2 objects]" << std::endl; + test (points_on_line); + + std::cerr << "[Testing line fitting on Segment_2 objects]" << std::endl; + test (points_on_line); + test (points_on_line); + + std::cerr << "[Testing line fitting on Circle_2 objects]" << std::endl; + test (points_on_line); + test (points_on_line); + + std::cerr << "[Testing line fitting on Triangle_2 objects]" << std::endl; + test (points_on_line); + test (points_on_line); + test (points_on_line); + + std::cerr << "[Testing line fitting on Iso_rectangle_2 objects]" << std::endl; + test (points_on_line); + test (points_on_line); + test (points_on_line); + + std::cerr << std::endl << "=== 3D ===" << std::endl << std::endl; + + std::cerr << "[Testing line fitting on Point_3 objects]" << std::endl; + test (points_on_line); + + std::cerr << "[Testing plane fitting on Point_3 objects]" << std::endl; + test (points_on_plane); + + std::cerr << "[Testing line fitting on Segment_3 objects]" << std::endl; + test (points_on_line); + test (points_on_line); + + std::cerr << "[Testing plane fitting on Segment_3 objects]" << std::endl; + test (points_on_plane); + test (points_on_plane); + + std::cerr << "[Testing line fitting on Sphere_3 objects]" << std::endl; + test (points_on_line); + test (points_on_line); + + std::cerr << "[Testing plane fitting on Sphere_3 objects]" << std::endl; + test (points_on_plane); + test (points_on_plane); + + std::cerr << "[Testing line fitting on Triangle_3 objects]" << std::endl; + test (points_on_line); + test (points_on_line); + test (points_on_line); + + std::cerr << "[Testing plane fitting on Triangle_3 objects]" << std::endl; + test (points_on_plane); + test (points_on_plane); + test (points_on_plane); + + std::cerr << "[Testing line fitting on Tetrahedron_3 objects]" << std::endl; + test (points_on_line); + test (points_on_line); + test (points_on_line); + test (points_on_line); + + std::cerr << "[Testing plane fitting on Tetrahedron_3 objects]" << std::endl; + test (points_on_plane); + test (points_on_plane); + test (points_on_plane); + test (points_on_plane); + + std::cerr << "[Testing line fitting on Iso_cuboid_3 objects]" << std::endl; + test (points_on_line); + test (points_on_line); + test (points_on_line); + test (points_on_line); + + std::cerr << "[Testing plane fitting on Iso_cuboid_3 objects]" << std::endl; + test (points_on_plane); + test (points_on_plane); + test (points_on_plane); + test (points_on_plane); + + return EXIT_SUCCESS; +}