Add swept extrusion

Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt

@@ -59,6 +59,7 @@ create_single_source_cgal_program("isotropic_remeshing_with_allow_move.cpp")
 create_single_source_cgal_program("triangle_mesh_autorefinement.cpp")
 create_single_source_cgal_program("soup_autorefinement.cpp")
 create_single_source_cgal_program("snap_polygon_soup.cpp")
+create_single_source_cgal_program("create_bended_cylinder.cpp")
 
 find_package(Eigen3 QUIET)
 include(CGAL_Eigen3_support)

Polygon_mesh_processing/examples/Polygon_mesh_processing/create_bended_cylinder.cpp

@@ -0,0 +1,55 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Polygon_mesh_processing/extrude.h>
+#include <CGAL/Polygon_mesh_processing/triangulate_hole.h>
+#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
+
+#include <fstream>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef K::Point_3 Point_3;
+typedef CGAL::Surface_mesh<Point_3> Mesh;
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+struct Face : public std::array<int,3>
+{
+  Face(int i, int j, int k)
+  {
+    (*this)[0] = i;
+    (*this)[1] = j;
+    (*this)[2] = k;
+  }
+};
+
+int main()
+{
+  //create a mesh of a disk
+  std::vector<Point_3> circle;
+  for (int i=0; i<25; ++i)
+    circle.push_back(Point_3(std::cos(2*CGAL_PI/25 * i), std::sin(2*CGAL_PI/25 * i), 0.));
+  std::vector<Face> triangles;
+  PMP::triangulate_hole_polyline(circle, std::back_inserter(triangles), CGAL::parameters::use_2d_constrained_delaunay_triangulation(true));
+  Mesh mesh;
+  PMP::polygon_soup_to_polygon_mesh(circle, triangles, mesh);
+
+  std::vector<Point_3> guide;
+
+  for (int i=0; i<120; ++i)
+    guide.push_back(Point_3(10*std::cos(2*CGAL_PI/25 * i), double(10*i)/24., 10*std::sin(2*CGAL_PI/25 * i)));
+
+  std::ofstream debug("guide.polylines.txt");
+  debug << guide.size();
+  for (auto p : guide)
+    debug << " " << p;
+  debug << std::endl;
+
+  Mesh out;
+  PMP::reverse_face_orientations(mesh);
+  std::ofstream("in.off") << mesh;
+
+  PMP::sweep_extrude(mesh, guide, out);
+
+  std::ofstream("swept_volume.off") << out;
+
+}

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/extrude.h

@@ -18,6 +18,7 @@
 
 
 #include <CGAL/Polygon_mesh_processing/orientation.h>
+#include <CGAL/Polygon_mesh_processing/transform.h>
 #include <CGAL/boost/graph/named_params_helper.h>
 #include <CGAL/Named_function_parameters.h>
 #include <CGAL/boost/graph/copy_face_graph.h>
@@ -298,5 +299,264 @@ void extrude_mesh(const InputMesh& input,
   extrude_mesh(input, output, bot,top, np_in, np_out);
 }
 
+/// \cond SKIP_IN_MANUAL
+namespace sweep_impl {
+
+template <class Kernel>
+void get_rotation(typename Kernel::FT angle,
+                  typename Kernel::Vector_3 axis,
+                  typename Kernel::Aff_transformation_3& rot)
+{
+  //create matrix of the rotation
+  const typename Kernel::RT c = cos(angle),
+                            s = sin(angle),
+                            x(axis.x()), y(axis.y()), z(axis.z());
+  typename Kernel::RT matrix[12] =
+  {
+    x*x*(1-c)+c, x*y*(1-c)-z*s, x*z*(1-c)+y*s, 0,
+    x*y*(1-c)+z*s, y*y*(1-c)+c, y*z*(1-c)-x*s, 0,
+    x*z*(1-c)-y*s, y*z*(1-c)+x*s, z*z*(1-c)+c, 0
+  };
+  rot = typename Kernel::Aff_transformation_3(matrix[0],matrix[1],matrix[2],
+                                              matrix[3],matrix[4],matrix[5],
+                                              matrix[6],matrix[7],matrix[8],
+                                              matrix[9],matrix[10],matrix[11]);
+}
+
+template<class Kernel>
+bool get_transfo(const typename Kernel::Point_3& p1,
+                 const typename Kernel::Point_3& p2,
+                 const typename Kernel::Point_3& p3,
+                 typename Kernel::Aff_transformation_3& rot,
+                 typename Kernel::FT& angle,
+                 bool split_angle = false)
+{
+  if(CGAL::collinear(p1,p2,p3)){
+    angle = 0;
+    return false;
+  }
+  //find the axis of the rotation:
+  typename Kernel::Vector_3 vec1(p1,p2), vec2(p2,p3);
+  typename Kernel::Vector_3 axis = CGAL::cross_product(vec1, vec2);
+  axis /= CGAL::approximate_sqrt(axis.squared_length());
+  //find the angle of the rotation:
+  angle = 180+CGAL::approximate_dihedral_angle(
+        p2, p2+axis,
+        p1, p3);
+  angle = split_angle ? angle * CGAL_PI/360.0
+                      : angle * CGAL_PI/180;
+
+
+
+  get_rotation<Kernel>(angle, axis, rot);
+  return true;
+}
+
+} // end of sweep_impl namespace
+
+/// \endcond
+
+/**
+ * \ingroup PMP_meshing_grp
+ *
+ * fills `output` with a closed mesh bounding the volume swept by `input`
+ * when moved along the path defined by the polyline `guide`.
+ * Note that initial position of `guide` does not influence
+ * the final result as it is translated to start on `initial`.
+ * `output` is oriented so that the faces corresponding to `input`
+ * in `output` have the same orientation.
+ *
+ * \attention `output` may be self intersecting.
+ *
+ * @tparam InputMesh a model of the concept `FaceListGraph`
+ * @tparam OutputMesh a model of the concept `FaceListGraph` and `MutableFaceGraph`
+ * @tparam Polyline a model of `SinglePassRange` containing points from the same CGAL kernel as the point of the vertex point map used for `OutputMesh`.
+ * @tparam NamedParameters1 a sequence of \ref bgl_namedparameters "Named Parameters" for `InputMesh`
+ * @tparam NamedParameters2 a sequence of \ref bgl_namedparameters "Named Parameters" for `OutputMesh`
+ *
+ * @param input an open surface mesh to extrude.
+ * @param output a surface mesh that will contain the result of the extrusion.
+ * @param v the vector defining the direction of the extrusion
+ * @param np_in an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+ *
+ * \cgalNamedParamsBegin
+ *   \cgalParamNBegin{vertex_point_map}
+ *     \cgalParamDescription{a property map associating points to the vertices of `input`}
+ *     \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits<InputMesh>::%vertex_descriptor`
+ *                    as key type and `%Point_3` as value type}
+ *     \cgalParamDefault{`boost::get(CGAL::vertex_point, input)`}
+ *     \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t`
+ *                     should be available for the vertices of `input`.}
+ *   \cgalParamNEnd
+ *   \cgalParamNBegin{use_rotations}
+ *     \cgalParamDescription{If `false`, `input` will be only translated along `guide`.
+ *                           If `true`, in addition to the translation applied to each layer,
+ *                           a rotation will also be applied using the angle between
+ *                           two consecutive segments along the path.}
+ *     \cgalParamType{`bool`}
+ *     \cgalParamDefault{`true`}
+ *   \cgalParamNEnd
+ * \cgalNamedParamsEnd
+ *
+ * @param np_out an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+ *
+ * \cgalNamedParamsBegin
+ *   \cgalParamNBegin{vertex_point_map}
+ *     \cgalParamDescription{a property map associating points to the vertices of `output`}
+ *     \cgalParamType{a class model of `ReadWritePropertyMap` with `boost::graph_traits<OutputMesh>::%vertex_descriptor`
+ *                    as key type and `%Point_3` as value type}
+ *     \cgalParamDefault{`boost::get(CGAL::vertex_point, output)`}
+ *     \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t`
+ *                     should be available for the vertices of `output`.}
+ *   \cgalParamNEnd
+ * \cgalNamedParamsEnd
+ */
+template <class InputMesh,
+          class OutputMesh,
+          class Polyline,
+          class CGAL_NP_TEMPLATE_PARAMETERS_1,
+          class CGAL_NP_TEMPLATE_PARAMETERS_2>
+void sweep_extrude(/* const */ InputMesh& input,
+                   const Polyline& guide,
+                   OutputMesh& output,
+                   const CGAL_NP_CLASS_1& np_in = parameters::default_values(),
+                   const CGAL_NP_CLASS_2& np_out = parameters::default_values())
+{
+  using CGAL::parameters::choose_parameter;
+  using CGAL::parameters::get_parameter;
+
+  typedef typename boost::graph_traits<OutputMesh>::halfedge_descriptor halfedge_descriptor;
+  typedef typename boost::graph_traits<OutputMesh>::vertex_descriptor vertex_descriptor;
+  typedef typename GetGeomTraits<InputMesh, CGAL_NP_CLASS_1>::type GT;
+  // GT gt = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
+  typedef typename GetVertexPointMap < OutputMesh, CGAL_NP_CLASS_2>::type VPM_out;
+  VPM_out vpm_out = choose_parameter(get_parameter(np_out, internal_np::vertex_point),
+                                     get_property_map(vertex_point, output));
+
+  bool use_rotational_term =
+    choose_parameter(get_parameter(np_in, internal_np::use_rotations), true);
+
+  typedef typename GT::Point_3 Point_3;
+  typedef typename GT::Vector_3 Vector_3;
+  typedef typename GT::Aff_transformation_3 Aff_transformation_3;
+
+  std::vector< std::vector< halfedge_descriptor > > previous_border_cycles, last_border_cycles;
+  copy_face_graph(input, output);
+  reverse_face_orientations(input);
+// NOTE: starting from here input will be used to hold the last layer and we will update it
+// accordingly all along the patch to be in the correct position for the last step.
+
+// translate guide points
+  std::vector<Point_3> guide_points(guide.begin(), guide.end());
+  std::size_t gs = guide.size();
+
+  Vector_3 initial_trans(guide_points[0], get(vpm_out, *output.vertices().begin()));
+  for (std::size_t i=0; i<gs; ++i)
+    guide_points[i] = guide_points[i] + initial_trans;
+
+// collect halfedges per boundary cycle of the first layer (using input that is reversed on purpose)
+  std::size_t nb_input_hedges = num_halfedges(input);
+  std::vector<bool> handled(nb_input_hedges,false);
+  for(halfedge_descriptor h : input.halfedges())
+    if (is_border(h, input) && !handled[h])
+    {
+      previous_border_cycles.resize(previous_border_cycles.size()+1);
+      last_border_cycles.resize(last_border_cycles.size()+1);
+      do{
+        previous_border_cycles.back().push_back(h);
+        last_border_cycles.back().push_back(h);
+        handled[h]=true;
+        h=input.next(h);
+      }
+      while(previous_border_cycles.back().front()!=h);
+    }
+  if (previous_border_cycles.empty()) return;
+
+  // create the middle parts
+  std::size_t nb_strips = guide.size()-1;
+  for (std::size_t s=1; s<nb_strips; ++s)
+  {
+    Vector_3 dir = guide_points[s] - guide_points[s-1];
+
+    // transformations
+    Aff_transformation_3 trans(CGAL::TRANSLATION, dir);
+    Aff_transformation_3 to_origin(CGAL::TRANSLATION, guide_points[s]-CGAL::ORIGIN);
+    Aff_transformation_3 rot;
+    transform(trans, input); // update for last layer
+
+    double angle = 0;
+    if (use_rotational_term)
+    {
+      sweep_impl::get_transfo<GT>(guide_points[s-1], guide_points[s], guide_points[s+1], rot, angle, false);
+      // update for last layer
+      transform(to_origin.inverse(), input);
+      transform(rot, input);
+      transform(to_origin, input);
+    }
+
+    for (std::size_t ci=0;ci<previous_border_cycles.size();++ci)
+    {
+      std::size_t nbh = previous_border_cycles[ci].size();
+      std::vector<halfedge_descriptor> current_hedges(nbh);
+
+      for(std::size_t i=0; i<nbh; ++i)
+      {
+        current_hedges[i] = output.add_edge();
+        vertex_descriptor v = output.add_vertex();
+        output.set_target( current_hedges[i], v );
+        output.set_halfedge(v, current_hedges[i]);
+        output.point(v) = output.point( output.source(previous_border_cycles[ci][i]) );
+        // transform the points
+        output.point(v)=output.point(v).transform(trans);
+        if (use_rotational_term)
+        {
+          output.point(v)=output.point(v).transform(to_origin.inverse());
+          output.point(v)=output.point(v).transform(rot);
+          output.point(v)=output.point(v).transform(to_origin);
+        }
+      }
+
+      for(std::size_t i=0; i<nbh; ++i)
+      {
+        std::size_t j = (i+1)%nbh;
+        vertex_descriptor v = output.target(current_hedges[i]);
+        output.set_target(output.opposite(current_hedges[j]), v);
+        output.set_next(current_hedges[i], current_hedges[j]);
+        output.set_next(output.opposite(current_hedges[j]),
+                        output.opposite(current_hedges[i]));
+      }
+
+      extrude_impl::create_strip(previous_border_cycles[ci], current_hedges, output);
+
+      for(std::size_t i=0; i<nbh; ++i)
+      {
+        previous_border_cycles[ci][i]=output.opposite(current_hedges[i]);
+      }
+    }
+  }
+
+  // handle the last strip
+  CGAL_assertion(!output.has_garbage());
+  std::size_t h_shift = output.num_halfedges();
+
+  //translate the last layer
+  Vector_3 dir = guide_points[gs-1] - guide_points[gs-2];
+  for (vertex_descriptor v : input.vertices())
+    input.point( v ) = input.point( v ) + dir;
+
+  output.join(input);
+
+  for (std::size_t i=0; i<last_border_cycles.size(); ++i)
+  {
+    std::size_t nbh=last_border_cycles[i].size();
+    for (std::size_t k=0; k<nbh; ++k)
+    {
+      last_border_cycles[i][k] = halfedge_descriptor( last_border_cycles[i][k].idx() + h_shift );
+      CGAL_assertion( is_border(last_border_cycles[i][k], output) );
+    }
+    extrude_impl::create_strip(previous_border_cycles[i], last_border_cycles[i], output);
+  }
+}
+
 }} //end CGAL::PMP
 #endif //CGAL_POLYGON_MESH_PROCESSING_EXTRUDE_H

STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h

@@ -238,6 +238,7 @@ CGAL_add_named_parameter(scan_angle_t, scan_angle_map, scan_angle_map)
 CGAL_add_named_parameter(scanline_id_t, scanline_id_map, scanline_id_map)
 CGAL_add_named_parameter(min_points_per_cell_t, min_points_per_cell, min_points_per_cell)
 CGAL_add_named_parameter(scalar_t, scalar_map, scalar_map)
+CGAL_add_named_parameter(use_rotations_t, use_rotations, use_rotations)
 
 // List of named parameters used in Surface_mesh_approximation package
 CGAL_add_named_parameter(verbose_level_t, verbose_level, verbose_level)