songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity
cgal/Straight_skeleton_extrusion_2/include/CGAL/extrude_skeleton.h

1142 lines
43 KiB
C++
Raw Blame History

// Copyright (c) 2023 GeometryFactory
// 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_SLS_EXTRUSION_EXTRUDE_SKELETON_H
#define CGAL_SLS_EXTRUSION_EXTRUDE_SKELETON_H
#include <CGAL/license/Straight_skeleton_extrusion_2.h>
#ifdef CGAL_SLS_DEBUG_DRAW
#include <CGAL/draw_straight_skeleton_2.h>
#include <CGAL/draw_polygon_2.h>
#include <CGAL/draw_polygon_with_holes_2.h>
#include <CGAL/draw_triangulation_2.h>
#endif
#include <CGAL/create_weighted_straight_skeleton_2.h>
#include <CGAL/create_weighted_straight_skeleton_from_polygon_with_holes_2.h>
#include <CGAL/create_weighted_offset_polygons_2.h>
#include <CGAL/create_weighted_offset_polygons_from_polygon_with_holes_2.h>
#include <CGAL/Straight_skeleton_2/Straight_skeleton_aux.h>
#include <CGAL/General_polygon_with_holes_2.h>
#include <CGAL/Polygon_mesh_processing/border.h>
#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
#include <CGAL/Polygon_mesh_processing/repair_polygon_soup.h>
#include <CGAL/Polygon_mesh_processing/stitch_borders.h>
#ifdef CGAL_SLS_OUTPUT_FILES
#include <CGAL/IO/polygon_soup_io.h>
#endif
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Projection_traits_3.h>
#include <CGAL/Triangulation_vertex_base_with_info_2.h>
#include <CGAL/mark_domain_in_triangulation.h>
#include <CGAL/Cartesian/function_objects.h>
#include <CGAL/enum.h>
#include <CGAL/IO/helpers.h>
#include <CGAL/Kernel_traits.h>
#include <CGAL/Polygon_2.h>
#include <boost/algorithm/clamp.hpp>
#include <boost/range/value_type.hpp>
#include <algorithm>
#include <iostream>
#include <map>
#include <type_traits>
#include <unordered_map>
#include <vector>
#include <optional>
#include <memory>
namespace CGAL {
namespace Straight_skeleton_extrusion {
namespace internal {
////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////
// Below is to handle vertical slabs.
enum class Slope
{
UNKNOWN = 0,
INWARD,
OUTWARD,
VERTICAL
};
template <typename FT>
inline constexpr FT default_extrusion_height()
{
return (std::numeric_limits<double>::max)();
}
// @todo Maybe this postprocessing is not really necessary? Do users really care if the point
// is not perfectly above the input contour edge (it generally cannot be anyway if the kernel
// is not exact except for some specific cases)?
#define CGAL_SLS_SNAP_TO_VERTICAL_SLABS
#ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS
// The purpose of this visitor is to snap back almost-vertical (see preprocessing_weights()) edges
// to actual vertical slabs.
template <typename HDS, typename GeomTraits>
typename GeomTraits::Point_2
snap_point_to_contour_halfedge_plane(const typename GeomTraits::Point_2& op,
typename HDS::Halfedge_const_handle ch)
{
using FT = typename GeomTraits::FT;
using Segment_2 = typename GeomTraits::Segment_2;
using Line_2 = typename GeomTraits::Line_2;
using Vertex_const_handle = typename HDS::Vertex_const_handle;
const Vertex_const_handle sv = ch->opposite()->vertex();
const Vertex_const_handle tv = ch->vertex();
if(sv->point().x() == tv->point().x())
{
// vertical edge
// std::cout << "vertical edge, snapping " << op << " to " << sv->point().x() << " " << op.y() << std::endl;
return { sv->point().x(), op.y() };
}
else if(sv->point().y() == tv->point().y())
{
// horizontal edge
// std::cout << "horizontal edge, snapping " << op << " to " << op.x() << " " << sv->point().y() << std::endl;
return { op.x(), sv->point().y() };
}
else
{
// Project orthogonally onto the halfedge
// @todo should the projection be along the direction of the other offset edge sharing this point?
Segment_2 s { sv->point(), tv->point() };
std::optional<Line_2> line = CGAL_SS_i::compute_normalized_line_coeffC2(s);
CGAL_assertion(bool(line)); // otherwise the skeleton would have failed already
FT px, py;
CGAL::line_project_pointC2(line->a(),line->b(),line->c(), op.x(),op.y(), px,py);
// std::cout << "snapping " << op << " to " << px << " " << py << std::endl;
return { px, py };
}
}
template <typename HDS, typename GeomTraits>
void snap_skeleton_vertex(typename HDS::Halfedge_const_handle hds_h,
typename HDS::Halfedge_const_handle contour_h,
std::map<typename GeomTraits::Point_2,
typename GeomTraits::Point_2>& snapped_positions)
{
typename HDS::Vertex_const_handle hds_tv = hds_h->vertex();
// this re-applies snapping towards contour_h even if the point was already snapped towards another contour
auto insert_result = snapped_positions.emplace(hds_tv->point(), hds_tv->point());
insert_result.first->second = snap_point_to_contour_halfedge_plane<HDS, GeomTraits>(insert_result.first->second, contour_h);
// std::cout << "snap_skeleton_vertex(V" << hds_tv->id() << " pt: " << hds_h->vertex()->point() << ")"
// << " to " << insert_result.first->second << std::endl;
}
template <typename GeomTraits, typename PointRange>
void apply_snapping(PointRange& points,
const std::map<typename GeomTraits::Point_2,
typename GeomTraits::Point_2>& snapped_positions)
{
using Point_3 = typename GeomTraits::Point_3;
for(Point_3& p3 : points)
{
auto it = snapped_positions.find({ p3.x(), p3.y() });
if(it != snapped_positions.end())
p3 = Point_3{it->second.x(), it->second.y(), p3.z()};
}
}
#endif // CGAL_SLS_SNAP_TO_VERTICAL_SLABS
template <typename StraightSkeleton_2,
typename OffsetBuilderTraits,
typename GeomTraits>
class Skeleton_offset_correspondence_builder_visitor
: public CGAL::Default_polygon_offset_builder_2_visitor<OffsetBuilderTraits, StraightSkeleton_2>
{
using Base = CGAL::Default_polygon_offset_builder_2_visitor<OffsetBuilderTraits, StraightSkeleton_2>;
using FT = typename GeomTraits::FT;
using Point_2 = typename GeomTraits::Point_2;
using SS_Halfedge_const_handle = typename StraightSkeleton_2::Halfedge_const_handle;
using HDS = typename StraightSkeleton_2::Base;
using HDS_Halfedge_const_handle = typename HDS::Halfedge_const_handle;
public:
Skeleton_offset_correspondence_builder_visitor(const StraightSkeleton_2& ss,
std::unordered_map<HDS_Halfedge_const_handle, Point_2>& offset_points
#ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS
, const FT vertical_weight
, std::map<Point_2, Point_2>& snapped_positions
#endif
)
: m_ss(ss)
, m_offset_points(offset_points)
#ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS
, m_vertical_weight(vertical_weight)
, m_snapped_positions(snapped_positions)
#endif
{ }
public:
void on_offset_contour_started() const
{
// std::cout << "~~ new contour ~~" << std::endl;
}
// can't modify the position yet because we need arrange_polygons() to still work properly
//
// @fixme on paper one could create a polygon thin-enough w.r.t. the max weight value
// such thatthere is a skeleton vertex that wants to be snapped to two different sides...
void on_offset_point(const Point_2& op,
SS_Halfedge_const_handle hook) const
{
CGAL_precondition(hook->is_bisector());
HDS_Halfedge_const_handle canonical_hook = (hook < hook->opposite()) ? hook : hook->opposite();
m_offset_points[canonical_hook] = op;
#ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS
CGAL_precondition(m_snapped_positions.count(op) == 0);
SS_Halfedge_const_handle contour_h1 = hook->defining_contour_edge();
CGAL_assertion(contour_h1->opposite()->is_border());
SS_Halfedge_const_handle contour_h2 = hook->opposite()->defining_contour_edge();
CGAL_assertion(contour_h2->opposite()->is_border());
const bool is_h1_vertical = (contour_h1->weight() == m_vertical_weight);
const bool is_h2_vertical = (contour_h2->weight() == m_vertical_weight);
// if both are vertical, it's the common vertex (which has to exist)
if(is_h1_vertical && is_h2_vertical)
{
CGAL_assertion(contour_h1->vertex() == contour_h2->opposite()->vertex() ||
contour_h2->vertex() == contour_h1->opposite()->vertex());
if(contour_h1->vertex() == contour_h2->opposite()->vertex())
m_snapped_positions[op] = contour_h1->vertex()->point();
else
m_snapped_positions[op] = contour_h2->vertex()->point();
}
else if(is_h1_vertical)
{
m_snapped_positions[op] = snap_point_to_contour_halfedge_plane<HDS, GeomTraits>(op, contour_h1);
}
else if(is_h2_vertical)
{
m_snapped_positions[op] = snap_point_to_contour_halfedge_plane<HDS, GeomTraits>(op, contour_h2);
}
#endif
}
private:
const StraightSkeleton_2& m_ss;
std::unordered_map<HDS_Halfedge_const_handle, Point_2>& m_offset_points;
#ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS
const FT m_vertical_weight;
std::map<Point_2, Point_2>& m_snapped_positions;
#endif
};
////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////
template <typename GeomTraits>
class Extrusion_builder
{
using Geom_traits = GeomTraits;
using FT = typename Geom_traits::FT;
using Point_2 = typename Geom_traits::Point_2;
using Segment_2 = typename Geom_traits::Segment_2;
using Line_2 = typename Geom_traits::Line_2;
using Point_3 = typename Geom_traits::Point_3;
using Vector_3 = typename Geom_traits::Vector_3;
using Polygon_2 = CGAL::Polygon_2<Geom_traits>;
using Polygon_with_holes_2 = CGAL::Polygon_with_holes_2<Geom_traits>;
using Offset_polygons = std::vector<std::shared_ptr<Polygon_2> >;
using Offset_polygons_with_holes = std::vector<std::shared_ptr<Polygon_with_holes_2> >;
using Straight_skeleton_2 = CGAL::Straight_skeleton_2<Geom_traits>;
using Straight_skeleton_2_ptr = std::shared_ptr<Straight_skeleton_2>;
using SS_Vertex_const_handle = typename Straight_skeleton_2::Vertex_const_handle;
using SS_Halfedge_const_handle = typename Straight_skeleton_2::Halfedge_const_handle;
using HDS = typename Straight_skeleton_2::Base;
using HDS_Vertex_const_handle = typename HDS::Vertex_const_handle;
using HDS_Halfedge_const_handle = typename HDS::Halfedge_const_handle;
using HDS_Face_handle = typename HDS::Face_handle;
// Standard CDT2 for the horizontal (z constant) faces
using Vb = CGAL::Triangulation_vertex_base_with_info_2<std::size_t, Geom_traits>;
using Vbb = CGAL::Triangulation_vertex_base_2<Geom_traits, Vb>;
using Fb = CGAL::Constrained_triangulation_face_base_2<Geom_traits>;
using TDS = CGAL::Triangulation_data_structure_2<Vb,Fb>;
using Itag = CGAL::No_constraint_intersection_requiring_constructions_tag;
using CDT = CGAL::Constrained_Delaunay_triangulation_2<Geom_traits, TDS, Itag>;
using CDT_Vertex_handle = typename CDT::Vertex_handle;
using CDT_Face_handle = typename CDT::Face_handle;
// Projection CDT2 for the lateral faces
using PK = CGAL::Projection_traits_3<Geom_traits>;
using PVbb = CGAL::Triangulation_vertex_base_with_info_2<std::size_t, PK>;
using PVb = CGAL::Triangulation_vertex_base_2<PK, PVbb>;
using PFb = CGAL::Constrained_triangulation_face_base_2<PK>;
using PTDS = CGAL::Triangulation_data_structure_2<PVb,PFb>;
using PCDT = CGAL::Constrained_Delaunay_triangulation_2<PK, PTDS, Itag>;
using PCDT_Vertex_handle = typename PCDT::Vertex_handle;
using PCDT_Face_handle = typename PCDT::Face_handle;
using Offset_builder_traits = CGAL::Polygon_offset_builder_traits_2<Geom_traits>;
using Visitor = Skeleton_offset_correspondence_builder_visitor<Straight_skeleton_2, Offset_builder_traits, Geom_traits>;
using Offset_builder = CGAL::Polygon_offset_builder_2<Straight_skeleton_2,
Offset_builder_traits,
Polygon_2,
Visitor>;
private:
Geom_traits m_gt;
public:
Extrusion_builder(const Geom_traits& gt) : m_gt(gt) { }
////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////
template <typename PolygonWithHoles,
typename PointRange, typename FaceRange>
void construct_horizontal_faces(const PolygonWithHoles& p,
const FT altitude,
PointRange& points,
FaceRange& faces,
const bool invert_faces = false)
{
#ifdef CGAL_SLS_DEBUG_DRAW
CGAL::draw(p);
#endif
CDT cdt;
try
{
cdt.insert_constraint(p.outer_boundary().begin(), p.outer_boundary().end(), true /*close*/);
for(auto h_it=p.holes_begin(); h_it!=p.holes_end(); ++h_it)
cdt.insert_constraint(h_it->begin(), h_it->end(), true /*close*/);
}
catch(const typename CDT::Intersection_of_constraints_exception&)
{
std::cerr << "Warning: Failed to triangulate horizontal face" << std::endl;
return;
}
std::size_t id = points.size(); // point ID offset (previous faces inserted their points)
for(CDT_Vertex_handle vh : cdt.finite_vertex_handles())
{
points.emplace_back(cdt.point(vh).x(), cdt.point(vh).y(), altitude);
vh->info() = id++;
}
#ifdef CGAL_SLS_DEBUG_DRAW
// CGAL::draw(cdt);
#endif
std::unordered_map<CDT_Face_handle, bool> in_domain_map;
boost::associative_property_map< std::unordered_map<CDT_Face_handle, bool> > in_domain(in_domain_map);
CGAL::mark_domain_in_triangulation(cdt, in_domain);
for(CDT_Face_handle f : cdt.finite_face_handles())
{
if(!get(in_domain, f))
continue;
// invert faces for the z=0 plane (bottom face)
if(invert_faces)
faces.push_back({f->vertex(0)->info(), f->vertex(2)->info(), f->vertex(1)->info()});
else
faces.push_back({f->vertex(0)->info(), f->vertex(1)->info(), f->vertex(2)->info()});
}
}
template <typename PointRange, typename FaceRange>
void construct_horizontal_faces(const Offset_polygons_with_holes& p_ptrs,
const FT altitude,
PointRange& points,
FaceRange& faces)
{
for(const auto& p_ptr : p_ptrs)
construct_horizontal_faces(*p_ptr, altitude, points, faces);
}
template <typename SLSFacePoints, typename PointRange, typename FaceRange>
void triangulate_skeleton_face(SLSFacePoints& face_points,
const bool invert_faces,
PointRange& points,
FaceRange& faces)
{
CGAL_precondition(face_points.size() >= 3);
// shift once to ensure that face_points[0] and face_points[1] are at z=0 and thus the normal is correct
std::rotate(face_points.rbegin(), face_points.rbegin() + 1, face_points.rend());
CGAL_assertion(face_points[0][2] == 0 && face_points[1][2] == 0);
const Vector_3 n = CGAL::cross_product(face_points[1] - face_points[0], face_points[2] - face_points[0]);
PK traits(n);
PCDT pcdt(traits);
try
{
pcdt.insert_constraint(face_points.begin(), face_points.end(), true /*close*/);
}
catch(const typename PCDT::Intersection_of_constraints_exception&)
{
std::cerr << "Warning: Failed to triangulate skeleton face" << std::endl;
return;
}
std::size_t id = points.size(); // point ID offset (previous faces inserted their points);
for(PCDT_Vertex_handle vh : pcdt.finite_vertex_handles())
{
points.push_back(pcdt.point(vh));
vh->info() = id++;
}
#ifdef CGAL_SLS_DEBUG_DRAW
// CGAL::draw(pcdt);
#endif
std::unordered_map<PCDT_Face_handle, bool> in_domain_map;
boost::associative_property_map< std::unordered_map<PCDT_Face_handle, bool> > in_domain(in_domain_map);
CGAL::mark_domain_in_triangulation(pcdt, in_domain);
for(PCDT_Face_handle f : pcdt.finite_face_handles())
{
if(!get(in_domain, f))
continue;
// invert faces for exterior skeletons
if(invert_faces)
faces.push_back({f->vertex(0)->info(), f->vertex(2)->info(), f->vertex(1)->info()});
else
faces.push_back({f->vertex(0)->info(), f->vertex(1)->info(), f->vertex(2)->info()});
}
}
// This version is for default height, so just gather all the full faces
void construct_lateral_faces(const Straight_skeleton_2& ss,
std::vector<Point_3>& points,
std::vector<std::vector<std::size_t> >& faces,
#ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS
const FT& vertical_weight,
std::map<Point_2, Point_2>& snapped_positions,
#endif
const bool ignore_frame_faces = false,
const bool invert_faces = false)
{
const HDS& hds = static_cast<const HDS&>(ss);
std::size_t fc = 0;
for(const HDS_Face_handle hds_f : CGAL::faces(hds))
{
std::vector<Point_3> face_points;
// If they exist (exterior skeleton), the first four faces of the SLS correspond
// to the outer frame, and should be ignored.
if(ignore_frame_faces && fc++ < 4)
continue;
HDS_Halfedge_const_handle hds_h = hds_f->halfedge(), done = hds_h;
#ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS
HDS_Halfedge_const_handle contour_h = hds_h->defining_contour_edge();
CGAL_assertion(hds_h == contour_h);
const bool is_vertical = (contour_h->weight() == vertical_weight);
#endif
do
{
HDS_Vertex_const_handle hds_tv = hds_h->vertex();
#ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS
// this computes the snapped position but does not change the geometry of the skeleton
if(is_vertical && !hds_tv->is_contour())
snap_skeleton_vertex<HDS, Geom_traits>(hds_h, contour_h, snapped_positions);
#endif
face_points.emplace_back(hds_tv->point().x(), hds_tv->point().y(), hds_tv->time());
hds_h = hds_h->next();
}
while(hds_h != done);
if(face_points.size() < 3)
{
std::cerr << "Warning: sm_vs has size 1 or 2: offset crossing face at a single point?" << std::endl;
continue;
}
triangulate_skeleton_face(face_points, invert_faces, points, faces);
}
}
void construct_lateral_faces(const Straight_skeleton_2& ss,
const Offset_builder& offset_builder,
const FT height,
std::vector<Point_3>& points,
std::vector<std::vector<std::size_t> >& faces,
const std::unordered_map<HDS_Halfedge_const_handle, Point_2>& offset_points,
#ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS
const FT& vertical_weight,
std::map<Point_2, Point_2>& snapped_positions,
#endif
const bool ignore_frame_faces = false,
const bool invert_faces = false)
{
CGAL_precondition(height != default_extrusion_height<FT>());
const bool extrude_upwards = is_positive(height);
const FT abs_height = CGAL::abs(height);
const HDS& hds = static_cast<const HDS&>(ss);
std::size_t fc = 0;
for(const HDS_Face_handle hds_f : CGAL::faces(hds))
{
// If they exist (exterior skeleton), the first four faces of the SLS correspond
// to the outer frame, and should be ignored.
if(ignore_frame_faces && fc++ < 4)
continue;
std::vector<Point_3> face_points;
HDS_Halfedge_const_handle hds_h = hds_f->halfedge(), done = hds_h;
#ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS
HDS_Halfedge_const_handle contour_h = hds_h->defining_contour_edge();
CGAL_assertion(hds_h == contour_h);
const bool is_vertical = (contour_h->weight() == vertical_weight);
#endif
do
{
HDS_Vertex_const_handle hds_sv = hds_h->opposite()->vertex();
HDS_Vertex_const_handle hds_tv = hds_h->vertex();
// Compare_offset_against_event_time compares offset to node->time(),
// so when the offset is greater or equal than the node->time(), the node is the face
auto compare_time_to_offset = [&](HDS_Vertex_const_handle node) -> CGAL::Comparison_result
{
if(node->is_contour())
return CGAL::LARGER; // offset > 0 and contour nodes' time is 0
else
return offset_builder.Compare_offset_against_event_time(abs_height, node);
};
const CGAL::Comparison_result sc = compare_time_to_offset(hds_sv);
const CGAL::Comparison_result tc = compare_time_to_offset(hds_tv);
// if the offset is crossing at the source, it will be added when seen as a target
// from the previous halfedge
if(sc != tc && sc != CGAL::EQUAL && tc != CGAL::EQUAL)
{
// std::cout << "sc != tc" << std::endl;
CGAL_assertion(sc != CGAL::EQUAL && tc != CGAL::EQUAL);
HDS_Halfedge_const_handle hds_off_h = hds_h;
if(hds_h->slope() == CGAL::NEGATIVE) // ensure same geometric point on both sides
hds_off_h = hds_off_h->opposite();
// The offset point must already been computed in the offset builder visitor
auto off_p = offset_points.find(hds_off_h);
CGAL_assertion(off_p != offset_points.end());
face_points.emplace_back(off_p->second.x(), off_p->second.y(), height);
}
if(tc != CGAL::SMALLER)
{
#ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS
if(is_vertical && !hds_tv->is_contour())
snap_skeleton_vertex<HDS, Geom_traits>(hds_h, contour_h, snapped_positions);
#endif
const Point_2& off_p = hds_tv->point();
// ->time() could be an approximation of the true time. If we are here, the target's time
// is smaller than the height, but we still sanitize it in case of numerical errors in ->time().
const FT time = boost::algorithm::clamp(hds_tv->time(), FT(0), abs_height);
face_points.emplace_back(off_p.x(), off_p.y(), extrude_upwards ? time : - time);
}
hds_h = hds_h->next();
}
while(hds_h != done);
if(face_points.size() < 3)
{
std::cerr << "Warning: sm_vs has size 1 or 2: offset crossing face at a single point?" << std::endl;
continue;
}
triangulate_skeleton_face(face_points, invert_faces, points, faces);
}
}
public:
// this is roughly "CGAL::create_interior_weighted_skeleton_and_offset_polygons_with_holes_2()",
// but we want to know the intermediate straight skeleton to build the lateral faces of the 3D mesh
template <typename PolygonWithHoles,
typename PointRange, typename FaceRange>
bool inward_construction(const PolygonWithHoles& pwh,
const std::vector<std::vector<FT> >& speeds,
const FT vertical_weight,
const FT height,
PointRange& points,
FaceRange& faces)
{
CGAL_precondition(!is_zero(height));
const FT abs_height = abs(height);
// bottom face (z=0)
construct_horizontal_faces(pwh, 0 /*height*/, points, faces, true /*invert faces*/);
// Avoid recomputing offset points multiple times
std::unordered_map<HDS_Halfedge_const_handle, Point_2> offset_points;
#ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS
// This is to deal with vertical slabs: we temporarily give a non-vertical slope to be able
// to construct the SLS, and we then snap back the position to the vertical planes.
// Note that points in non-vertical faces are also changed a bit
std::map<Point_2, Point_2> snapped_positions;
#endif
Straight_skeleton_2_ptr ss_ptr;
bool is_default_extrusion_height = (height == default_extrusion_height<FT>());
// @partial_wsls_pwh interior SLS of weighted polygons with holes can have skeleton faces with holes
// The current postprocessing is in the function EnforceSimpleConnectedness, but it has not yet
// been made compatible with partial skeletons.
if(is_default_extrusion_height || pwh.number_of_holes() != 0)
{
ss_ptr = CGAL::create_interior_weighted_straight_skeleton_2(
CGAL_SS_i::vertices_begin(pwh.outer_boundary()),
CGAL_SS_i::vertices_end(pwh.outer_boundary()),
pwh.holes_begin(), pwh.holes_end(),
std::begin(speeds[0]), std::end(speeds[0]),
std::next(std::begin(speeds)), std::end(speeds),
m_gt);
}
else
{
ss_ptr = CGAL_SS_i::create_partial_interior_weighted_straight_skeleton_2(
abs_height,
CGAL_SS_i::vertices_begin(pwh.outer_boundary()),
CGAL_SS_i::vertices_end(pwh.outer_boundary()),
pwh.holes_begin(), pwh.holes_end(),
std::begin(speeds[0]), std::end(speeds[0]),
std::next(std::begin(speeds)), std::end(speeds),
m_gt);
}
if(!ss_ptr)
{
std::cerr << "Error: encountered an error during skeleton construction" << std::endl;
return false;
}
#ifdef CGAL_SLS_DEBUG_DRAW
// print_straight_skeleton(*ss_ptr);
CGAL::draw(*ss_ptr);
#endif
if(is_default_extrusion_height)
{
#ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS
construct_lateral_faces(*ss_ptr, points, faces, vertical_weight, snapped_positions);
#else
construct_lateral_faces(*ss_ptr, points, faces);
#endif
}
else // height != default_extrusion_height<FT>()
{
#ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS
Visitor visitor(*ss_ptr, offset_points, vertical_weight, snapped_positions);
#else
Visitor visitor(*ss_ptr, offset_points);
#endif
Offset_builder ob(*ss_ptr, Offset_builder_traits(), visitor);
Offset_polygons raw_output;
ob.construct_offset_contours(abs_height, std::back_inserter(raw_output));
Offset_polygons_with_holes output = CGAL::arrange_offset_polygons_2<Polygon_with_holes_2>(raw_output);
construct_horizontal_faces(output, height, points, faces);
#ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS
construct_lateral_faces(*ss_ptr, ob, height, points, faces, offset_points, vertical_weight, snapped_positions);
#else
construct_lateral_faces(*ss_ptr, ob, height, points, faces, offset_points);
#endif
}
// if the height is negative (extruding downwards), we need to invert all the faces
if(is_negative(height))
{
for(auto& f : faces)
std::swap(f[0], f[1]);
}
#ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS
apply_snapping<Geom_traits>(points, snapped_positions);
#endif
return true;
}
template <typename PolygonWithHoles,
typename PointRange, typename FaceRange>
bool outward_construction(const PolygonWithHoles& pwh,
const std::vector<std::vector<FT> >& speeds,
const FT vertical_weight,
const FT height,
PointRange& points,
FaceRange& faces)
{
CGAL_precondition(!is_zero(height));
CGAL_precondition(height != default_extrusion_height<FT>()); // was checked before, this is just a reminder
const FT abs_height = abs(height);
// bottom face (z=0)
construct_horizontal_faces(pwh, 0 /*height*/, points, faces, true /*invert faces*/);
// Avoid recomputing offset points multiple times
std::unordered_map<HDS_Halfedge_const_handle, Point_2> offset_points;
#ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS
// This is to deal with vertical slabs: we temporarily give a non-vertical slope to be able
// to construct the SLS, and we then snap back the position to the vertical planes.
// Note that points in non-vertical faces are also changed a bit
std::map<Point_2, Point_2> snapped_positions;
#endif
Offset_polygons raw_output; // accumulates for both the outer boundary and the holes
// the exterior of a polygon with holes is the exterior of its outer boundary,
// and the interior of its inverted holes
//
// Start with the outer boundary
{
Straight_skeleton_2_ptr ss_ptr = CGAL_SS_i::create_partial_exterior_weighted_straight_skeleton_2(
abs_height,
CGAL_SS_i::vertices_begin(pwh.outer_boundary()),
CGAL_SS_i::vertices_end(pwh.outer_boundary()),
std::begin(speeds[0]), std::end(speeds[0]),
m_gt);
if(!ss_ptr)
{
std::cerr << "Error: encountered an error during outer skeleton construction" << std::endl;
return false;
}
#ifdef CGAL_SLS_DEBUG_DRAW
// print_straight_skeleton(*ss_ptr);
CGAL::draw(*ss_ptr);
#endif
#ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS
Visitor visitor(*ss_ptr, offset_points, vertical_weight, snapped_positions);
#else
Visitor visitor(*ss_ptr, offset_points);
#endif
Offset_builder ob(*ss_ptr, Offset_builder_traits(), visitor);
ob.construct_offset_contours(abs_height, std::back_inserter(raw_output));
// Manually filter the offset of the outer frame
std::swap(raw_output[0], raw_output.back());
raw_output.pop_back();
#ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS
construct_lateral_faces(*ss_ptr, ob, height, points, faces, offset_points,
vertical_weight, snapped_positions,
true /*ignore frame faces*/, true /*invert faces*/);
#else
construct_lateral_faces(*ss_ptr, ob, height, points, faces, offset_points,
true /*ignore frame faces*/, true /*invert faces*/);
#endif
}
// now, deal with the holes
std::size_t hole_id = 0;
for(auto hit=pwh.holes_begin(); hit!=pwh.holes_end(); ++hit, ++hole_id)
{
Polygon_2 hole = *hit; // intentional copy
hole.reverse_orientation();
// this is roughly "CGAL::create_exterior_weighted_skeleton_and_offset_polygons_with_holes_2()",
// but we want to know the intermediate straight skeleton to build the lateral faces of the 3D mesh
std::vector<Polygon_2> no_holes;
std::vector<std::vector<FT> > no_speeds;
Straight_skeleton_2_ptr ss_ptr = CGAL_SS_i::create_partial_interior_weighted_straight_skeleton_2(
abs_height,
CGAL_SS_i::vertices_begin(hole),
CGAL_SS_i::vertices_end(hole),
std::begin(no_holes), std::end(no_holes),
std::begin(speeds[1 + hole_id]), std::end(speeds[1 + hole_id]),
std::begin(no_speeds), std::end(no_speeds),
m_gt);
if(!ss_ptr)
{
std::cerr << "Error: encountered an error during skeleton construction" << std::endl;
return EXIT_FAILURE;
}
#ifdef CGAL_SLS_DEBUG_DRAW
// print_straight_skeleton(*ss_ptr);
CGAL::draw(*ss_ptr);
#endif
#ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS
Visitor visitor(*ss_ptr, offset_points, vertical_weight, snapped_positions);
#else
Visitor visitor(*ss_ptr, offset_points);
#endif
Offset_builder ob(*ss_ptr, Offset_builder_traits(), visitor);
ob.construct_offset_contours(abs_height, std::back_inserter(raw_output));
#ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS
construct_lateral_faces(*ss_ptr, ob, height, points, faces, offset_points,
vertical_weight, snapped_positions,
false /*no outer frame*/, true /*invert faces*/);
#else
construct_lateral_faces(*ss_ptr, ob, height, points, faces, offset_points,
false /*no outer frame*/, true /*invert faces*/);
#endif
}
// - the exterior offset of the outer boundary is built by creating an extra frame and turning
// the outer boundary into a hole. Hence, it needs to be reversed back to proper orientation
// - the exterior offset of the holes is built by reversing the holes and computing an internal
// skeleton. Hence, the result also needs to be reversed.
for(std::shared_ptr<Polygon_2> ptr : raw_output)
ptr->reverse_orientation();
Offset_polygons_with_holes output = CGAL::arrange_offset_polygons_2<Polygon_with_holes_2>(raw_output);
construct_horizontal_faces(output, height, points, faces);
// if the height is negative (extruding downwards), we need to invert all the faces
if(is_negative(height))
{
for(auto& f : faces)
std::swap(f[0], f[1]);
}
#ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS
apply_snapping<Geom_traits>(points, snapped_positions);
#endif
return true;
}
};
////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////
// convert angles (in degrees) to weights, and handle vertical angles
template <typename FT, typename AngleRange>
void convert_angles(AngleRange& angles)
{
CGAL_precondition(!angles.empty());
auto angle_to_weight = [](const FT& angle) -> FT
{
CGAL_precondition(0 < angle && angle < 180);
// @todo should this be an epsilon around 90°? As theta goes to 90°, tan(theta) goes to infinity
// and thus we could get numerical issues (overflows) if the kernel is not exact
if(angle == 90)
return 0;
else
return std::tan(CGAL::to_double(angle * CGAL_PI / 180));
};
for(auto& contour_angles : angles)
for(FT& angle : contour_angles)
angle = angle_to_weight(angle);
}
// handle vertical angles (inf speed)
// returns
// - whether the weights are positive or negative (inward / outward)
// - whether the input weights are valid
// - the weight of vertical slabs
template <typename FT, typename WeightRange>
std::tuple<Slope, bool, FT>
preprocess_weights(WeightRange& weights)
{
CGAL_precondition(!weights.empty());
Slope slope = Slope::UNKNOWN;
FT max_value = 0; // non-inf, maximum absolute value
for(auto& contour_weights : weights)
{
for(FT& w : contour_weights)
{
// '0' means a vertical slab, aka 90° angle (see preprocess_angles())
if(w == 0)
continue;
// determine whether weights indicate all inward or all outward
if(slope == Slope::UNKNOWN)
{
// w is neither 0 nor inf here
slope = (w > 0) ? Slope::INWARD : Slope::OUTWARD;
}
else if(slope == Slope::INWARD && w < 0)
{
std::cerr << "Error: mixing positive and negative weights is not yet supported" << std::endl;
return {Slope::UNKNOWN, false, FT(-1)};
}
else if(slope == Slope::OUTWARD && w > 0)
{
std::cerr << "Error: mixing positive and negative weights is not yet supported" << std::endl;
return {Slope::UNKNOWN, false, FT(-1)};
}
// if we are going outwards, it is just an interior skeleton with opposite weights
w = CGAL::abs(w);
if(w > max_value)
max_value = w;
}
}
if(slope == Slope::UNKNOWN)
{
std::cerr << "Warning: all edges vertical?" << std::endl;
slope = Slope::VERTICAL;
}
// Take a weight which is a large % of the max value to ensure there's no ambiguity
//
// Since the max value might not be very close to 90°, take the max between of the large-% weight
// and the weight corresponding to an angle of 89.9999999°
const FT weight_of_89d9999999 = 572957787.3425436; // tan(89.9999999°)
const FT scaled_max = (std::max)(weight_of_89d9999999, FT(1e3) * max_value);
for(auto& contour_weights : weights)
{
for(FT& w : contour_weights)
{
if(w == FT(0))
w = scaled_max;
}
}
return {slope, true, scaled_max};
}
template <typename PolygonWithHoles,
typename WeightRange,
typename PolygonMesh,
typename NamedParameters>
bool extrude_skeleton(const PolygonWithHoles& pwh,
WeightRange& weights,
PolygonMesh& out,
const NamedParameters& np)
{
namespace PMP = ::CGAL::Polygon_mesh_processing;
using Polygon_2 = typename PolygonWithHoles::Polygon_2;
using Point = typename boost::range_value<Polygon_2>::type;
using Default_kernel = typename Kernel_traits<Point>::type;
using Geom_traits = typename internal_np::Lookup_named_param_def<internal_np::geom_traits_t,
NamedParameters,
Default_kernel>::type;
using FT = typename Geom_traits::FT;
using Point_3 = typename Geom_traits::Point_3;
using parameters::choose_parameter;
using parameters::get_parameter;
const bool verbose = choose_parameter(get_parameter(np, internal_np::verbose), false);
Geom_traits gt = choose_parameter<Geom_traits>(get_parameter(np, CGAL::internal_np::geom_traits));
const FT height = choose_parameter(get_parameter(np, internal_np::maximum_height), default_extrusion_height<FT>());
Slope slope;
bool valid_input;
FT vertical_weight;
std::tie(slope, valid_input, vertical_weight) = preprocess_weights<FT>(weights);
if(!valid_input)
{
if(verbose)
std::cerr << "Error: invalid input weights" << std::endl;
return false;
}
if(verbose)
{
switch(slope)
{
case Slope::UNKNOWN: std::cout << "Slope is UNKNOWN??" << std::endl; break;
case Slope::INWARD: std::cout << "Slope is INWARD" << std::endl; break;
case Slope::OUTWARD: std::cout << "Slope is OUTWARD" << std::endl; break;
case Slope::VERTICAL: std::cout << "Slope is VERTICAL" << std::endl; break;
}
}
// End of preprocessing, start the actual skeleton computation
if(slope != Slope::INWARD && height == default_extrusion_height<FT>())
{
if(verbose)
std::cerr << "Error: height must be specified when using an outward (or vertical) slope" << std::endl;
return false;
}
// build a soup, to be converted to a mesh afterwards
std::vector<Point_3> points;
std::vector<std::vector<std::size_t> > faces;
points.reserve(2 * pwh.outer_boundary().size()); // just a reasonable guess
faces.reserve(2 * pwh.outer_boundary().size() + 2*pwh.number_of_holes());
Extrusion_builder<Geom_traits> builder(gt);
bool res;
if(slope != Slope::OUTWARD) // INWARD or VERTICAL
res = builder.inward_construction(pwh, weights, vertical_weight, height, points, faces);
else
res = builder.outward_construction(pwh, weights, vertical_weight, height, points, faces);
if(!res)
return false;
#ifdef CGAL_SLS_OUTPUT_FILES
// This soup provides one connected component per edge of the input polygon
CGAL::IO::write_polygon_soup("extruded_skeleton_soup.off", points, faces, CGAL::parameters::stream_precision(17));
#endif
// Convert the triangle soup to a triangle mesh
PMP::merge_duplicate_points_in_polygon_soup(points, faces);
if(!PMP::is_polygon_soup_a_polygon_mesh(faces))
PMP::orient_polygon_soup(points, faces);
CGAL_assertion(PMP::is_polygon_soup_a_polygon_mesh(faces));
PMP::polygon_soup_to_polygon_mesh(points, faces, out);
CGAL_warning(is_valid_polygon_mesh(out) && is_closed(out));
return true;
}
} // namespace internal
} // namespace Straight_skeleton_extrusion
// Documented in the Straight_skeleton_2 package
template <typename PolygonWithHoles,
typename PolygonMesh,
typename NamedParameters = parameters::Default_named_parameters>
bool extrude_skeleton(const PolygonWithHoles& pwh,
PolygonMesh& out,
const NamedParameters& np = parameters::default_values(),
std::enable_if_t<CGAL_SS_i::has_Hole_const_iterator<PolygonWithHoles>::value>* = nullptr)
{
namespace SSEI = Straight_skeleton_extrusion::internal;
using Polygon_2 = typename PolygonWithHoles::General_polygon_2;
using K = typename Kernel_traits<typename boost::range_value<Polygon_2>::type>::type;
using FT = typename K::FT;
using Default_speed_type = std::vector<std::vector<FT> >;
using parameters::choose_parameter;
using parameters::is_default_parameter;
using parameters::get_parameter_reference;
const bool has_weights = !(is_default_parameter<NamedParameters, internal_np::weights_param_t>::value);
const bool has_angles = !(is_default_parameter<NamedParameters, internal_np::angles_param_t>::value);
Default_speed_type def_speeds; // never used, but needed for compilation
if(has_weights)
{
// not using "Angles" here is on purpose, I want to copy the range but the NP is ref only
auto weights = choose_parameter(get_parameter_reference(np, internal_np::weights_param), def_speeds);
return SSEI::extrude_skeleton(pwh, weights, out, np);
}
else if(has_angles)
{
// not using "Weights" here is on purpose, I want to copy the range but the NP is ref only
auto angles = choose_parameter(get_parameter_reference(np, internal_np::angles_param), def_speeds);
SSEI::convert_angles<FT>(angles);
return SSEI::extrude_skeleton(pwh, angles, out, np);
}
else // neither angles nor weights were provided
{
std::vector<std::vector<FT> > uniform_weights;
uniform_weights.reserve(pwh.number_of_holes() + 1);
uniform_weights.push_back(std::vector<FT>(pwh.outer_boundary().size(), FT(1)));
for(const auto& hole : pwh.holes())
uniform_weights.push_back(std::vector<FT>(hole.size(), FT(1)));
return SSEI::extrude_skeleton(pwh, uniform_weights, out, np);
}
}
////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////
template <typename Polygon, // without holes
typename PolygonMesh,
typename NamedParameters = parameters::Default_named_parameters>
bool extrude_skeleton(const Polygon& p,
PolygonMesh& out,
const NamedParameters& np = parameters::default_values(),
std::enable_if_t<!CGAL_SS_i::has_Hole_const_iterator<Polygon>::value>* = nullptr)
{
using Polygon_with_holes = CGAL::General_polygon_with_holes_2<Polygon>;
return extrude_skeleton(Polygon_with_holes(p), out, np);
}
} // namespace CGAL
#endif // CGAL_SLS_EXTRUSION_EXTRUDE_SKELETON_H
Reference in New Issue View Git Blame Copy Permalink