2D projection of a Surface Mesh and get the edges

[maximecharriere]

Question

I have a 3D mesh that use the Surface Mesh structure (a scan of a head in the pic).
I need to calculate the projection of this 3D mesh onto a plane and get the edges of this projection (the red line in the pic).
At the moment I'm using the Polygon_mesh_slicer and cutting a slice of the mesh in the middle, but I would prefer to have a projection.

What would be your approach? What package and function would you use?

What I've found

I have found in the doc the Projection_traits_3 but I don't really understand how to use it.

MWE (using the Polygon_mesh_slicer)

using Scalar     = double;
using Kernel     = typename CGAL::Simple_cartesian<Scalar>;
using Point_3    = typename Kernel::Point_3;
using Mesh       = typename CGAL::Surface_mesh<Point_3>;
using Plane_3    = typename Kernel::Plane_3;
using Polyline_3 = typename std::vector<Point_3>;

/* Import scan */
Mesh scan;
std::ifstream fd(filename, std::ios::in | std::ios::binary);
CGAL::IO::read_PLY(fd, scan);
fd.close()

/* Slice head in the midle*/
CGAL::Polygon_mesh_slicer<Mesh, Kernel> slicer(scan);
std::vector<Polyline_3> polylines;
slicer(Plane_3(1,0,0,0), std::back_inserter(polylines));

/* Draw result */
typedef CGAL::Polygon_2<Kernel> Polygon_2;
typedef CGAL::Polygon_set_2<Kernel> Polygon_set_2;

Polygon_set_2 polygonSet;
for(const auto& polyline : polylines) {
  Polygon_2 polygon;
  for(const auto& point : polyline) {
    polygon.push_back(Point_2(point.x(), point.z()));
  }
  polygonSet.insert(polygon);
}
CGAL::draw(polygonSet);

[afabri]

The projection traits can be used with a constrained 2D triangulation. It exists for the xy plane (by just ignoring z of 3D points) but also for an arbitrary plane. The projection direction is the normal of the plane.

[sloriot]

The brute force-way to get that is to put all edges in a 2-CDT and extract the contour accessible from the infinite vertex. You can also filter faces that have a normal opposite to the projection direction if you want to improve things.

I know @efifogel already did that for a personal project , maybe using https://doc.cgal.org/latest/Envelope_3/index.html.

[maximecharriere]

Okay, I've spent the last 5 hours figuring out how Constrained_Delaunay_triangulation_2 and Projection_traits_3 work, and I think I've arrived at a solution.

using Kernel    = typename CGAL::Exact_predicates_inexact_constructions_kernel;
using GT        = typename CGAL::Projection_traits_yz_3<Kernel>;
using Point_3   = typename Kernel::Point_3;
using Mesh      = typename CGAL::Surface_mesh<Point_3>;
               
using Itag      = typename CGAL::Exact_predicates_tag;
using CDT       = typename CGAL::Constrained_Delaunay_triangulation_2<GT, CGAL::Default, Itag>;


int main()
{
  Mesh scan;
  std::ifstream fd("D:\\cube.ply", std::ios::in | std::ios::binary);
  CGAL::IO::read_PLY(fd, scan);
  fd.close();

  GT gt;
  CDT cdt(gt);
  cdt.insert(scan.points().begin(), scan.points().end());

  CDT::Vertex_circulator start = cdt.infinite_vertex()->incident_vertices();
  CDT::Vertex_circulator vertex = start;
  do
  {
    std::cout << vertex->point() << std::endl;
    vertex++;
  } while(vertex != start);
}

Last small question. Is it not possible to have a Ranges for the incident_vertices() function to use a C++11 for loop like that ?
The part with the Vertex_circulator in my code is not very beautiful.

for(const auto& vertex: cdt.infinite_vertex()->incident_vertices()) {
   std::cout << vertex->point() << std::endl;
}

[afabri]

You must work harder. Your code will list the vertices on the convex hull. That works for the cube as it is convex, but a head is not

You have to walk along the chain of edges marked as constraint. Note that you will also get edges "inside" the head.

[sloriot]

That's the first step of the mark_domain() here

[efifogel]

Hi, You can use fairly simple code to find the 2D silhouette of a polygonal surface. It is described in Section 2.3 of the book "CGAL Arrangements and Their Applications"; se, e.g., http://acg.cs.tau.ac.il/cgal-arrangement-book. All the code samples can be obtained from here: http://acg.cs.tau.ac.il/cgal-arrangement-book/source-code-data-and-miscellaneous-files . The main cpp is called polytope_projection.cpp <http://acg.cs.tau.ac.il/cgal-arrangement-book/examples/chapter-2/polytope_projection.cpp>. It depends on some header files that can be obtained from the same page. The example first computes the convex hull of the polygonal surface, but I think that this can be skipped just as well. Efi ____ _ ____ _ /_____/_) o /__________ __ // (____ ( ( ( (_/ (_/-(-'_(/ _/

…

On Wed, 6 Apr 2022 at 12:17, Sebastien Loriot ***@***.***> wrote: That's the first step of the mark_domain() here <https://doc.cgal.org/latest/Triangulation_2/Triangulation_2_2polygon_triangulation_8cpp-example.html> — Reply to this email directly, view it on GitHub <#6472 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABVBNOA6KPEFPZITOZRN62LVDVJDFANCNFSM5SS66ALQ> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[maximecharriere]

Thank you ! I will read this chapter and try your different solutions. I'll get back to you when I know more.

[maximecharriere]

Are you sure this solution works for a concave polygon?
In the book, chapter 2.3, it's written : "Given a convex polytope P, the program obtains the outline of the shadow of P cast on the xy-plane,..."

[efifogel]

Hi, The general framework works for any set of polygons. You need to (i) project the polygons onto a 2D plane, (ii) insert the projected edges into an arrangement, and (iii) remove the internal edges. If the input is a mesh, you can mark the edges, then traverse all facets; for each facet, traverse all boundary edges; if the edge hasn't been projected yet, mark it and project it. If the mesh is a 2D manifold (watertight), filter out any facet whose normal and projected direction form an obtuse angle. Efi ____ _ ____ _ /_____/_) o /__________ __ // (____ ( ( ( (_/ (_/-(-'_(/ _/

…

On Wed, 6 Apr 2022 at 17:20, Charrière Maxime ***@***.***> wrote: Are you sure this solution works for a concave polygon? In the book, chapter 2.3, it's written : *"Given a convex polytope P, the program obtains the outline of the shadow of P cast on the xy-plane,..."* — Reply to this email directly, view it on GitHub <#6472 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABVBNODIDZEVZ26RFHO4VW3VDWMTBANCNFSM5SS66ALQ> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[maximecharriere]

Hello,
After working a bit harder, I came up with this solution. It's based on the Triangulation_2/polygon_triangulation.cpp example.
To make the calculation faster, I could still filter the faces according to the angle between their normal and the direction of the projection, as you suggested.

Wouldn't it be a good thing to include a function in CGAL to get the borders (shadow from an infinit source of light) of any class based on the BGL concepts such as Surface_mesh , Polyhedron and the 2D triangulation .
You give the direction of the projection and you get the edges on the border.
Or at least publish a example dedicated on this feature.

My code

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Projection_traits_3.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Triangulation_face_base_with_info_2.h>
#include <CGAL/Point_set_3.h>
#include <CGAL/draw_point_set_3.h>


struct FaceInfo2
{
  FaceInfo2() {}
  int nesting_level;
  bool in_domain() {
    return nesting_level % 2 == 1;
  }
};


using Kernel              = typename CGAL::Exact_predicates_inexact_constructions_kernel;
using Point_3             = typename Kernel::Point_3;
using Vector_3            = typename Kernel::Vector_3;
using Mesh                = typename CGAL::Surface_mesh<Point_3>;
using Point_set_3         = typename CGAL::Point_set_3<Point_3>;

using Projection_traits_3 = typename CGAL::Projection_traits_3<Kernel>;
using Vb                  = typename CGAL::Triangulation_vertex_base_2<Projection_traits_3>;
using Fbb                 = typename CGAL::Triangulation_face_base_with_info_2<FaceInfo2, Projection_traits_3>;
using Fb                  = typename CGAL::Constrained_triangulation_face_base_2<Projection_traits_3, Fbb>;
using TDS                 = typename CGAL::Triangulation_data_structure_2<Vb, Fb>;
using Itag                = typename CGAL::Exact_predicates_tag;
using CDT                 = typename CGAL::Constrained_Delaunay_triangulation_2<Projection_traits_3, TDS, Itag>;
using Face_handle         = typename CDT::Face_handle;


// Based on the example "Triangulation_2/polygon_triangulation.cpp"
void getBorders(const CDT& cdt, std::list<CDT::Edge>& borders)
{
  for(CDT::Face_handle f : cdt.all_face_handles()) {
    f->info().nesting_level = -1;
  }

  std::list<Face_handle> queue;
  queue.push_back(cdt.infinite_face());
  while(!queue.empty()) {
    Face_handle fh = queue.front();
    queue.pop_front();
    if(fh->info().nesting_level == -1) {
      fh->info().nesting_level = 0;
      for(int i = 0; i < 3; i++) {
        CDT::Edge e(fh, i);
        Face_handle n = fh->neighbor(i);
        if(n->info().nesting_level == -1) {
          if(cdt.is_constrained(e)) borders.push_back(e);
          else queue.push_back(n);
        }
      }
    }
  }
}

int main(const int ac, const char* const av[])
{
  // Import 3D mesh
  Mesh mesh;
  std::ifstream fd("D:\\scan.ply", std::ios::in | std::ios::binary);
  CGAL::IO::read_PLY(fd, mesh);
  fd.close();

  // Construct the 2D triangulation of the 3D mesh projected in a given direction
  Vector_3 direction(1, 0, 0); //slower than using Projection_traits_yz_3
  Projection_traits_3 gt(direction);
  CDT cdt(gt);
  for(const auto& edge : mesh.edges()) {
    const Point_3& start = mesh.point(mesh.vertex(edge, 0));
    const Point_3& end = mesh.point(mesh.vertex(edge, 1));
    cdt.insert_constraint(start, end);
  }

  // Get the externals borders of the triangulation
  std::list<CDT::Edge> borders;
  getBorders(cdt, borders);


#ifdef CGAL_USE_BASIC_VIEWER

  // Draw result
  Point_set_3 ps;
  for(const auto& border : borders) {
    const CDT::Face& face = *border.first;
    const int i = border.second;
    const Point_3& start = face.vertex(face.cw(i))->point();
    const Point_3& end = face.vertex(face.ccw(i))->point();
    ps.insert(start);
    ps.insert(end);
  }
  CGAL::draw(ps);

#else

  // Print result
  for(const auto& border : borders) {
    const CDT::Face& face = *border.first;
    const int i = border.second;
    const Point_3& start = face.vertex(face.cw(i))->point();
    const Point_3& end = face.vertex(face.ccw(i))->point();
    std::cout << start << "  -->  " << end << std::endl;
  }

#endif // CGAL_USE_BASIC_VIEWER

  return EXIT_SUCCESS;
}

[sloriot]

Nice! Why don't you have a closed curve as output?

[afabri]

because he draws a point set.

[maximecharriere]

I tried to generate a closed Polygon_2, but the list of edges does not form a chain of edges, and I haven't found a simple solution to order it. So I displayed only the points and not a closed curve.

[sloriot]

Using the projection traits you can even see where are the edges contributing to the silhouette.

[maximecharriere]

How did you generate this result? I would be interested to see your code, if possible.

[sloriot]

I dumped the polyline in a file and opened it in the polyhedron demo.

[afabri]

Concerning the question why this is not in CGAL. Because it is too brute force. This may sound unfriendly or arrogant, but that is not my intention. If it was in CGAL it would be something that preprocesses the mesh and had an internal data structure in order to generate a silhouette in less than linear time for a given orientation.

[maximecharriere]

Okay, I understand. This is a high-level feature, and it is not the purpose of CGAL to provide this type of function. CGAL provides data structures and basic mathematical functions, and this kind of problem must be programmed directly by the user of the library.

[afabri]

Note also that you must not insert a constraint which in the projection is just a point.

[efifogel]

Hi,

As mentioned above, the book "CGAL Arrangements and Their Applications" (see, e.g., http://acg.cs.tau.ac.il/cgal-arrangement-book) touches this topic.
In particular look for
(i) Section 2.3 "Application: Obtaining Silhouettes of Polytopes", and
(ii) Section 8.4 "Application: Obtaining Silhouettes of Polyhedra"

All the code samples can be obtained from here.
The main files are called polytope_projection.cpp and polyhedron_projection.cpp, respectively.

The files depends on some header files that can be obtained from the same page.

Efi