CMakeLists.txt

# Created by the script cgal_create_CMakeLists
# This is the CMake script for compiling a set of CGAL applications.

cmake_minimum_required(VERSION 3.1...3.13)

project( data )


# CGAL and its components
find_package( CGAL REQUIRED QUIET OPTIONAL_COMPONENTS  )

if ( NOT CGAL_FOUND )

  message(STATUS "This project requires the CGAL library, and will not be compiled.")
  return()  

endif()


# Boost and its components
find_package( Boost REQUIRED )
if ( NOT Boost_FOUND )

  message(STATUS "This project requires the Boost library, and will not be compiled.")

  return()  

endif()

# include for local directory

# include for local package


# Creating entries for all C++ files with "main" routine
# ##########################################################


create_single_source_cgal_program( "remesh_polyhedral_surface.cpp" )
create_single_source_cgal_program( "orient_polygon_soup.cpp" )

find_package( TBB REQUIRED )
include(CGAL_target_use_TBB)
CGAL_target_use_TBB(orient_polygon_soup)

orient_polygon_soup.cpp

#include <array>
#include <vector>
#include <iostream>
#include <fstream>
#include <utility>

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>

#include <CGAL/IO/OBJ_reader.h>
#include <CGAL/IO/OFF_reader.h>

#include <CGAL/Surface_mesh.h>

#include <CGAL/AABB_tree.h>
#include <CGAL/AABB_traits.h>
#include <CGAL/AABB_face_graph_triangle_primitive.h>

#include <CGAL/Polygon_mesh_processing/orient_polygon_soup.h>
#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
#include <CGAL/Polygon_mesh_processing/shape_predicates.h>
#include <CGAL/Polygon_mesh_processing/compute_normal.h>
#include <CGAL/Polygon_mesh_processing/connected_components.h>
#include <CGAL/Polygon_mesh_processing/orientation.h>
#include <CGAL/Polygon_mesh_processing/stitch_borders.h>

#ifdef CGAL_LINKED_WITH_TBB
#include <tbb/parallel_for.h>
#else
#endif // CGAL_LINKED_WITH_TBB
#include <boost/iterator/counting_iterator.hpp>
#include <boost/iterator/filter_iterator.hpp>

#if CGAL_VERSION_NR > 1050000000
template <class PointRange, class PolygonRange>
bool
duplicate_incompatible_edge_in_polygon_soup(PointRange& points,
                                            PolygonRange& polygons)
{
  std::size_t inital_nb_pts = points.size();
  typedef CGAL::Polygon_mesh_processing::internal::
    Polygon_soup_orienter<PointRange, PolygonRange> Orienter;

  Orienter orienter(points, polygons);
  orienter.fill_edge_map();
  // make edges to duplicate
  for(std::size_t i1=0;i1<points.size();++i1)
    for(const typename Orienter::Internal_map_type::value_type& i2_and_pids : orienter.edges[i1])
      if (i2_and_pids.second.size() > 1)
        orienter.set_edge_marked(i1,i2_and_pids.first,orienter.marked_edges);
  orienter.duplicate_singular_vertices();

  return inital_nb_pts==points.size();
}
#else
template <class PointRange, class PolygonRange>
bool
duplicate_incompatible_edge_in_polygon_soup(PointRange& points,
                                            PolygonRange& polygons)
{
  std::size_t inital_nb_pts = points.size();
  typedef CGAL::Polygon_mesh_processing::internal::
    Polygon_soup_orienter<PointRange, PolygonRange> Orienter;

  Orienter orienter(points, polygons);
  orienter.fill_edge_map();
  // make edges to duplicate
  for(typename Orienter::Edge_map_iterator eit = orienter.edges.begin(); eit!= orienter.edges.end(); ++eit)
    if (eit->second.size() > 1)
      orienter.set_edge_marked(eit->first.first,eit->first.second,orienter.marked_edges);
  orienter.duplicate_singular_vertices();

  return inital_nb_pts==points.size();
}
#endif

template <class Concurrency_tag, class K, class Point_3, class Triangle_face, class TriangleMesh>
void
orient_triangle_soup_with_reference_triangle_mesh(
  const TriangleMesh& tm_ref,
  std::vector<Point_3>& points,
  std::vector<Triangle_face>& triangles)
{
  namespace PMP = CGAL::Polygon_mesh_processing;

  typedef boost::graph_traits<TriangleMesh> GrT;
  typedef typename GrT::face_descriptor face_descriptor;

  typedef std::function<bool(face_descriptor)> Face_predicate;
  Face_predicate is_not_deg =
    [&tm_ref](face_descriptor f)
    {
      return !PMP::is_degenerate_triangle_face(f, tm_ref);
    };

  // build a tree filtereing degenerate faces
  typedef CGAL::AABB_face_graph_triangle_primitive<TriangleMesh> Primitive;
  typedef CGAL::AABB_traits<K, Primitive> Tree_traits;

  boost::filter_iterator<Face_predicate, typename GrT::face_iterator>
    begin(is_not_deg, faces(tm_ref).begin(), faces(tm_ref).end()),
    end(is_not_deg, faces(tm_ref).end(), faces(tm_ref).end());

  CGAL::AABB_tree<Tree_traits> tree(begin, end, tm_ref);

  // now orient the faces
  tree.build();
  tree.accelerate_distance_queries();
  auto process_facet =
    [&points, &tree, &tm_ref, &triangles](std::size_t fid) {
      const auto& p0 = points[triangles[fid][0]];
      const auto& p1 = points[triangles[fid][1]];
      const auto& p2 = points[triangles[fid][2]];
      const typename K::Point_3 mid = CGAL::centroid(p0, p1, p2);
      std::pair<typename K::Point_3, face_descriptor> pt_and_f =
        tree.closest_point_and_primitive(mid);
      auto face_ref_normal = PMP::compute_face_normal(pt_and_f.second, tm_ref);
      if(face_ref_normal * cross_product(p1-p0, p2-p0) < 0) {
        std::swap(triangles[fid][1], triangles[fid][2]);
      }
    };

#if !defined(CGAL_LINKED_WITH_TBB)
  CGAL_static_assertion_msg (!(boost::is_convertible<Concurrency_tag, CGAL::Parallel_tag>::value),
                             "Parallel_tag is enabled but TBB is unavailable.");
#else
  if (boost::is_convertible<Concurrency_tag,CGAL::Parallel_tag>::value)
    tbb::parallel_for(std::size_t(0), triangles.size(), std::size_t(1), process_facet);
  else
#endif
    std::for_each(
      boost::counting_iterator<std::size_t> (0),
      boost::counting_iterator<std::size_t> (triangles.size()),
      process_facet);
}

template <class TriangleMesh>
void merge_reversible_connected_components(TriangleMesh& tm)
{
  namespace PMP = CGAL::Polygon_mesh_processing;

  typedef boost::graph_traits<TriangleMesh> GrT;
  typedef typename GrT::face_descriptor face_descriptor;
  typedef typename GrT::halfedge_descriptor halfedge_descriptor;

  typedef typename boost::property_map<TriangleMesh, CGAL::vertex_point_t >::type Vpm;
  typedef typename boost::property_traits<Vpm>::value_type Point_3;
  Vpm vpm = get(CGAL::vertex_point, tm);

  typedef typename boost::property_map<TriangleMesh, CGAL::dynamic_face_property_t<std::size_t> >::type Fidmap;
  Fidmap fim = get(CGAL::dynamic_face_property_t<std::size_t>(), tm);
  std::size_t i=0;
  for (face_descriptor f : faces(tm))
    put(fim, f, i++);

  Fidmap f_cc_ids = get(CGAL::dynamic_face_property_t<std::size_t>(), tm);
  std::size_t nb_cc = PMP::connected_components(tm, f_cc_ids, CGAL::parameters::face_index_map(fim));

  std::vector<std::size_t> nb_faces_per_cc(nb_cc, 0);
  for (face_descriptor f : faces(tm))
    nb_faces_per_cc[ get(f_cc_ids, f) ]+=1;

  std::map< std::pair<Point_3, Point_3>, std::vector<halfedge_descriptor> > border_hedges_map;
  std::vector<halfedge_descriptor> border_hedges;
  typedef typename boost::property_map<TriangleMesh, CGAL::dynamic_halfedge_property_t<std::size_t> >::type Hidmap;
  Hidmap him = get(CGAL::dynamic_halfedge_property_t<std::size_t>(), tm);
  const std::size_t DV(-1);

  // fill endpoints -> hedges
  for (halfedge_descriptor h : halfedges(tm))
  {
    if ( CGAL::is_border(h, tm) )
    {
      put(him, h, DV);
      border_hedges_map[std::make_pair(get(vpm, source(h, tm)), get(vpm, target(h, tm)))].push_back(h);
      border_hedges.push_back(h);
    }
  }

  // set the id of boundary cc
  std::size_t id=0;
  for(halfedge_descriptor h : border_hedges)
  {
    if (get(him,h) == DV)
    {
      for (halfedge_descriptor hh : CGAL::halfedges_around_face(h, tm))
      {
        put(him, hh, id);
      }
      ++id;
    }
  }

  // check boundary fully matching
  std::vector<bool> border_cycle_handled(id, false);
  std::set<std::size_t> ccs_to_reverse;
  for ( const auto& p : border_hedges_map )
  {
    const std::vector<halfedge_descriptor>& hedges = p.second;
    if ( hedges.size()==2 && !border_cycle_handled[ get(him, hedges[0]) ]
                          && !border_cycle_handled[ get(him, hedges[1]) ] )
    {
      border_cycle_handled[ get(him, hedges[0]) ] = true;
      border_cycle_handled[ get(him, hedges[1]) ] = true;

      halfedge_descriptor h2=hedges[1];
      bool did_break=false;
      for(halfedge_descriptor h1 : CGAL::halfedges_around_face(hedges[0], tm))
      {
        if (get(vpm, target(h1,tm)) != get(vpm, target(h2,tm)))
        {
          did_break=true;
          break;
        }
        h2=next(h2,tm);
      }
      if (!did_break && h2==hedges[1])
      {
        std::size_t cc_id1 = get(f_cc_ids, face(opposite(hedges[0], tm), tm)),
                    cc_id2 = get(f_cc_ids, face(opposite(hedges[1], tm), tm));

        if (cc_id1!=cc_id2)
        {
          if (ccs_to_reverse.count(cc_id1)==0 && ccs_to_reverse.count(cc_id2)==0)
            ccs_to_reverse.insert( nb_faces_per_cc[cc_id1] < nb_faces_per_cc[cc_id2] ? cc_id1 : cc_id2 );
        }
      }
    }
  }

  // reverse ccs and stitches boundaries
  std::vector<face_descriptor> faces_to_reverse;
  for (face_descriptor f : faces(tm))
    if ( ccs_to_reverse.count( get(f_cc_ids, f) ) != 0 )
      faces_to_reverse.push_back(f);

  if ( !faces_to_reverse.empty() )
  {
    PMP::reverse_face_orientations(faces_to_reverse, tm);
    PMP::stitch_borders(tm);
  }
}


using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Polyhedron = CGAL::Surface_mesh<K::Point_3>;

#if defined(CGAL_LINKED_WITH_TBB)
#define TAG CGAL::Parallel_tag
#else
#define TAG CGAL::Sequential_tag
#endif

int main(int argc, char** argv)
{
  if(argc != 4) {
    std::cerr << argv[0]
              << "ref_OBJ input_OFF output_OFF\n";
    return EXIT_FAILURE;
  }
  // Load OBJ
  std::vector<K::Point_3> points_ref;
  std::vector<std::vector<std::size_t> > faces_ref;

  {
    std::ifstream in(argv[1]);
    if(!in || !CGAL::read_OBJ(in,points_ref,faces_ref))
    {
      return 1;
    }
  }
  // Load OFF
  std::vector<K::Point_3> points_in;
  std::vector<std::array<std::size_t, 3> > faces_in;
  {
    std::ifstream in(argv[2]);
    if(!in || !CGAL::read_OFF(in,points_in,faces_in))
    {
      return 1;
    }
  }

  Polyhedron poly_ref;
  namespace PMP = CGAL::Polygon_mesh_processing;
  PMP::orient_polygon_soup(points_ref,faces_ref);
  PMP::polygon_soup_to_polygon_mesh(points_ref, faces_ref, poly_ref);
  if (!CGAL::is_triangle_mesh(poly_ref))
  {
    std::cerr << "Input geometry is not triangulated." << std::endl;
    return EXIT_FAILURE;
  }

  orient_triangle_soup_with_reference_triangle_mesh<TAG, K>(poly_ref, points_in, faces_in);
  duplicate_incompatible_edge_in_polygon_soup(points_in, faces_in);

  std::ofstream soup_out("soup_out.off");
  soup_out << "OFF\n" << points_in.size() << " " << faces_in.size() << " 0\n";
  std::copy(points_in.begin(), points_in.end(), std::ostream_iterator<K::Point_3>(soup_out,"\n"));
  for (const auto& a : faces_in)
    soup_out << "3 " << a[0] << " " << a[1]<< " " << a[2] << "\n";

  Polyhedron poly_out;
  PMP::polygon_soup_to_polygon_mesh(points_in, faces_in, poly_out);
  std::ofstream("before_merge.off") << std::setprecision(17) << poly_out;

  merge_reversible_connected_components(poly_out);
  std::ofstream(argv[3]) << std::setprecision(17) << poly_out;

  return 0;
}