mirror of https://github.com/CGAL/cgal
Add post-processing: 2D optimization along OBB axes
This commit is contained in:
Mael Rouxel-Labbé
parent
0c2b32ea48
commit
fbfec59341
|
|
@ -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<K> 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;
|
||||
};
|
||||
|
|
|
|||
|
|
@ -19,6 +19,7 @@
|
|||
#include <CGAL/Optimal_bounding_box/internal/fitness_function.h>
|
||||
#include <CGAL/Optimal_bounding_box/internal/helper.h>
|
||||
#include <CGAL/Optimal_bounding_box/internal/nelder_mead_functions.h>
|
||||
#include <CGAL/Optimal_bounding_box/internal/optimize_2.h>
|
||||
#include <CGAL/Optimal_bounding_box/internal/population.h>
|
||||
|
||||
#include <CGAL/Random.h>
|
||||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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 <CGAL/license/Optimal_bounding_box.h>
|
||||
|
||||
#include <CGAL/ch_akl_toussaint.h>
|
||||
#include <CGAL/min_quadrilateral_2.h>
|
||||
#include <CGAL/Polygon_2.h>
|
||||
#include <CGAL/number_type_config.h>
|
||||
|
||||
#include <iostream>
|
||||
#include <iterator>
|
||||
#include <limits>
|
||||
#include <utility>
|
||||
#include <vector>
|
||||
|
||||
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 <typename Traits, typename PointRange>
|
||||
std::pair<typename Traits::FT, typename Traits::FT>
|
||||
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<Point_2> 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<Point_2> extreme_points;
|
||||
ch_akl_toussaint(points_2D.begin(), points_2D.end(), std::back_inserter(extreme_points), traits);
|
||||
|
||||
CGAL::Polygon_2<Traits> pol;
|
||||
CGAL::Min_quadrilateral_default_traits_2<Traits> 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 <typename PointRange, typename Traits>
|
||||
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<typename boost::range_value<PointRange>::type, Point>::value));
|
||||
|
||||
std::vector<Point> rotated_points;
|
||||
rotated_points.size();
|
||||
|
||||
FT xmin, ymin, zmin, xmax, ymax, zmax;
|
||||
xmin = ymin = zmin = FT{std::numeric_limits<double>::max()};
|
||||
xmax = ymax = zmax = FT{std::numeric_limits<double>::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<FT, 3> angles;
|
||||
std::array<FT, 3> 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<decltype(it)>::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
|
||||
|
|
@ -26,6 +26,7 @@
|
|||
#include <CGAL/Bbox_3.h>
|
||||
#include <CGAL/boost/graph/helpers.h>
|
||||
#include <CGAL/convex_hull_3.h>
|
||||
#include <CGAL/Convex_hull_traits_3.h>
|
||||
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
|
||||
#include <CGAL/Iso_cuboid_3.h>
|
||||
#include <CGAL/Iterator_range.h>
|
||||
|
|
@ -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<Point> ch_points;
|
||||
extreme_points_3(points, std::back_inserter(ch_points), traits);
|
||||
CGAL::Convex_hull_traits_3<Traits> 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;
|
||||
|
|
|
|||
Loading…
Reference in New Issue