From fbfec59341dccb23cfe2cb6ea13adbfebfb716c2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Wed, 25 Mar 2020 15:25:59 +0100 Subject: [PATCH] Add post-processing: 2D optimization along OBB axes --- .../Concepts/OrientedBoundingBoxTraits.h | 10 +- .../Optimal_bounding_box/internal/evolution.h | 23 +- .../internal/optimize_2.h | 216 ++++++++++++++++++ .../oriented_bounding_box.h | 4 +- 4 files changed, 238 insertions(+), 15 deletions(-) create mode 100644 Optimal_bounding_box/include/CGAL/Optimal_bounding_box/internal/optimize_2.h diff --git a/Optimal_bounding_box/doc/Optimal_bounding_box/Concepts/OrientedBoundingBoxTraits.h b/Optimal_bounding_box/doc/Optimal_bounding_box/Concepts/OrientedBoundingBoxTraits.h index 9eb5aee1daa..f31d11252ec 100644 --- a/Optimal_bounding_box/doc/Optimal_bounding_box/Concepts/OrientedBoundingBoxTraits.h +++ b/Optimal_bounding_box/doc/Optimal_bounding_box/Concepts/OrientedBoundingBoxTraits.h @@ -6,7 +6,7 @@ The concept `OrientedBoundingBoxTraits_3` describes the requirements of the trai used in the function `CGAL::oriented_bounding_box()`, and in particular the need for a 3x3 matrix type. -\cgalRefines `ConvexHullTraits_3` +\cgalRefines `Kernel` \cgalHasModel `CGAL::Oriented_bounding_box_traits_3` @@ -21,18 +21,10 @@ public: /// and be compatible with the type `Point_3` typedef CGAL::Aff_transformation_3 Aff_transformation_3; - /// A construction object that must provide the function operator: - /// `CGAL::Bbox_3 operator()(const Point_3&)`, - /// which returns an axis-aligned bounding box that contains the point - typedef unspecified_type Construct_bbox_3; - /// A 3x3 matrix type; model of the concept `SvdTraits::Matrix` and which supports /// matrix-matrix and scalar-matrix multiplication, as well as matrix-matrix addition typedef unspecified_type Matrix; /// Returns the unitary matrix `Q` obtained in the QR-decomposition of the matrix `m` Matrix get_Q(const Matrix& m) const; - - /// Returns the 3D box construction functor - Construct_bbox_3 construct_bbox_3_object() const; }; diff --git a/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/internal/evolution.h b/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/internal/evolution.h index 72b62959043..d6649fbe28a 100644 --- a/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/internal/evolution.h +++ b/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/internal/evolution.h @@ -19,6 +19,7 @@ #include #include #include +#include #include #include @@ -157,15 +158,27 @@ public: std::cout << std::endl; #endif + // optimize the current best rotation by using the exact OBB 2D algorithm + // along the axes of the current best OBB + m_best_v = &(m_population.get_best_vertex()); + Matrix& best_m = m_best_v->matrix(); + +#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG + std::cout << "new best matrix: " << std::endl << best_m << std::endl; + std::cout << "fitness: " << m_best_v->fitness() << std::endl; +#endif + + optimize_along_OBB_axes(best_m, m_points, m_traits); + m_best_v->fitness() = compute_fitness(best_m, m_points, m_traits); + // stopping criteria - m_best_v = &(get_best_vertex(m_population)); - const FT new_fit_value = m_best_v->fitness_value(); + const FT new_fit_value = m_best_v->fitness(); const FT difference = new_fit_value - prev_fit_value; #ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG - const Matrix& best_m = m_best_v->matrix(); - std::cout << "new best matrix: " << std::endl << best_m << std::endl; - std::cout << "value difference with previous: " << difference << std::endl; + std::cout << "post 2D optimization matrix: " << std::endl << best_m << std::endl; + std::cout << "new fit value: " << new_fit_value << std::endl; + std::cout << "difference: " << difference << std::endl; #endif if(CGAL::abs(difference) < tolerance * new_fit_value) diff --git a/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/internal/optimize_2.h b/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/internal/optimize_2.h new file mode 100644 index 00000000000..74306a03c4d --- /dev/null +++ b/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/internal/optimize_2.h @@ -0,0 +1,216 @@ +// Copyright (c) 2020 GeometryFactory (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s) : Mael Rouxel-Labbé +// +#ifndef CGAL_OPTIMAL_BOUNDING_BOX_INTERNAL_OPTIMIZE_2_H +#define CGAL_OPTIMAL_BOUNDING_BOX_INTERNAL_OPTIMIZE_2_H + +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace CGAL { +namespace Optimal_bounding_box { +namespace internal { + +enum PROJECTION_DIRECTION +{ + ALONG_X = 0, + ALONG_Y, + ALONG_Z +}; + +// Now, we would like to do all of this with projection traits... Unfortunately, it's missing +// a couple of functors (Has_on_negative_side_2, for example), which are necessary in +// CGAL::Min_quadrilateral_default_traits_2. And while we know we have a generic case here, +// it's a bit tedious to get something well-defined in the generic case: for example, +// what should Has_on_negative_side_2 do if the Line_3 is orthogonal to the projection plane? +// +// So, easier to just bail out to real 2D... +template +std::pair +compute_2D_deviation(const PointRange& points, + const PROJECTION_DIRECTION dir, + const Traits& traits) +{ + typedef typename Traits::FT FT; + typedef typename Traits::Point_2 Point_2; + typedef typename Traits::Vector_2 Vector_2; + typedef typename Traits::Point_3 Point_3; + + std::vector points_2D; + points_2D.reserve(points.size()); + for(const Point_3& pt : points) + { + if(dir == ALONG_X) + points_2D.emplace_back(pt.y(), pt.z()); + else if(dir == ALONG_Y) + points_2D.emplace_back(pt.x(), pt.z()); + else if(dir == ALONG_Z) + points_2D.emplace_back(pt.x(), pt.y()); + } + + std::vector extreme_points; + ch_akl_toussaint(points_2D.begin(), points_2D.end(), std::back_inserter(extreme_points), traits); + + CGAL::Polygon_2 pol; + CGAL::Min_quadrilateral_default_traits_2 mrt; + CGAL::min_rectangle_2(extreme_points.begin(), extreme_points.end(), std::back_inserter(pol), mrt); + + CGAL_assertion(pol.size() == 4); + CGAL_assertion(pol.is_counterclockwise_oriented()); + + // Compute the angle between the angle necessary to rotate the rectangle onto the reference frame + auto bot_pos = pol.bottom_vertex(); + auto next_pos = bot_pos; + ++next_pos; + if(next_pos == pol.vertices_end()) + next_pos = pol.begin(); + + const Point_2& p = *bot_pos; + const Point_2& q = *next_pos; + + const Vector_2 pq = traits.construct_vector_2_object()(p, q); + double n = sqrt(to_double(traits.compute_squared_length_2_object()(pq))); + CGAL_assertion(n != 0.); + + const double dot = pq.x(); // that's the scalar product of PQ with V(1, 0) (Ox) + double cosine = dot / n; + + if(cosine > 1.) + cosine = 1.; + if(cosine < -1.) + cosine = -1.; + + double theta = std::acos(cosine); + if(theta > 0.25 * CGAL_PI) // @todo is there a point to this + theta = 0.5 * CGAL_PI - theta; + + return std::make_pair(pol.area(), theta); +} + +template +void optimize_along_OBB_axes(typename Traits::Matrix& rot, + const PointRange& points, + const Traits& traits) +{ + typedef typename Traits::FT FT; + typedef typename Traits::Point_3 Point; + typedef typename Traits::Matrix Matrix; + + CGAL_static_assertion((std::is_same::type, Point>::value)); + + std::vector rotated_points; + rotated_points.size(); + + FT xmin, ymin, zmin, xmax, ymax, zmax; + xmin = ymin = zmin = FT{std::numeric_limits::max()}; + xmax = ymax = zmax = FT{std::numeric_limits::lowest()}; + + for(const Point& pt : points) + { + const FT rx = rot(0, 0) * pt.x() + rot(0, 1) * pt.y() + rot(0, 2) * pt.z(); + const FT ry = rot(1, 0) * pt.x() + rot(1, 1) * pt.y() + rot(1, 2) * pt.z(); + const FT rz = rot(2, 0) * pt.x() + rot(2, 1) * pt.y() + rot(2, 2) * pt.z(); + + rotated_points.emplace_back(rx, ry, rz); + + xmin = (std::min)(xmin, rx); + ymin = (std::min)(ymin, ry); + zmin = (std::min)(zmin, rz); + xmax = (std::max)(xmax, rx); + ymax = (std::max)(ymax, ry); + zmax = (std::max)(zmax, rz); + } + + const FT lx = xmax - xmin; + const FT ly = ymax - ymin; + const FT lz = zmax - zmin; + + std::array angles; + std::array volumes; + + FT area_xy; + std::tie(area_xy, angles[0]) = compute_2D_deviation(rotated_points, ALONG_Z, traits); + volumes[0] = lz * area_xy; + + FT area_xz; + std::tie(area_xz, angles[1]) = compute_2D_deviation(rotated_points, ALONG_Y, traits); + volumes[1] = ly * area_xz; + + FT area_yz; + std::tie(area_yz, angles[2]) = compute_2D_deviation(rotated_points, ALONG_X, traits); + volumes[2] = lx * area_yz; + + auto it = std::min_element(volumes.begin(), volumes.end()); + typename std::iterator_traits::difference_type d = std::distance(volumes.begin(), it); + +#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG + std::cout << "volumes: " << volumes[0] << " " << volumes[1] << " " << volumes[2] << std::endl; + std::cout << "angles: " << angles[0] << " " << angles[1] << " " << angles[2] << std::endl; + std::cout << "min at " << d << std::endl; +#endif + + if(d == 0) // Along_Z + { + const double c = std::cos(angles[0]); + const double s = std::sin(angles[0]); + + Matrix opt; + opt.set(0, 0, c); opt.set(0, 1, s); opt.set(0, 2, 0); + opt.set(1, 0, -s); opt.set(1, 1, c); opt.set(1, 2, 0); + opt.set(2, 0, 0); opt.set(2, 1, 0); opt.set(2, 2, 1); + + rot = opt * rot; + } + else if(d == 1) // Along_Y + { + const double c = std::cos(angles[1]); + const double s = std::sin(angles[1]); + + Matrix opt; + opt.set(0, 0, c); opt.set(0, 1, 0); opt.set(0, 2, -s); + opt.set(1, 0, 0); opt.set(1, 1, 1); opt.set(1, 2, 0); + opt.set(2, 0, s); opt.set(2, 1, 0); opt.set(2, 2, c); + + rot = opt * rot; + } + else if(d == 2) // Along_X + { + const double c = std::cos(angles[2]); + const double s = std::sin(angles[2]); + + Matrix opt; + opt.set(0, 0, 1); opt.set(0, 1, 0); opt.set(0, 2, 0); + opt.set(1, 0, 0); opt.set(1, 1, c); opt.set(1, 2, s); + opt.set(2, 0, 0); opt.set(2, 1, -s); opt.set(2, 2, c); + + rot = opt * rot; + } + else + { + CGAL_assertion(false); + } +} + +} // namespace internal +} // namespace Optimal_bounding_box +} // namespace CGAL + +#endif // CGAL_OPTIMAL_BOUNDING_BOX_INTERNAL_OPTIMIZE_2_H diff --git a/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/oriented_bounding_box.h b/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/oriented_bounding_box.h index 7dafd269b66..5a4a8819b6f 100644 --- a/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/oriented_bounding_box.h +++ b/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/oriented_bounding_box.h @@ -26,6 +26,7 @@ #include #include #include +#include #include #include #include @@ -171,7 +172,8 @@ void construct_oriented_bounding_box(const PointRange& points, if(use_ch) // construct the convex hull to reduce the number of points { std::vector ch_points; - extreme_points_3(points, std::back_inserter(ch_points), traits); + CGAL::Convex_hull_traits_3 CH_traits; + extreme_points_3(points, std::back_inserter(ch_points), CH_traits); #ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG std::cout << ch_points.size() << " points on the convex hull" << std::endl;