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 <tbb/parallel_for.h>

#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>

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

K::Vector_3
get_face_normal(Polyhedron::Face_index f, const Polyhedron& tm)
{
  Polyhedron::Halfedge_index h = tm.halfedge(f);
  const K::Point_3& p1 = tm.point(source(h, tm));
  const K::Point_3& p2 = tm.point(target(h, tm));
  const K::Point_3& p3 = tm.point(target(next(h, tm), tm));

  return CGAL::cross_product(p2-p1, p3-p1);
}

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;
  }
  // remove degenerate faces
  std::vector<Polyhedron::Face_index> faces_to_remove;
  for (Polyhedron::Face_index f : poly_ref.faces())
    if( PMP::is_degenerate_triangle_face(f, poly_ref))
      faces_to_remove.push_back(f);
  for (Polyhedron::Face_index f : faces_to_remove)
    CGAL::Euler::remove_face(halfedge(f,poly_ref), poly_ref);
  faces_to_remove.clear();

  // now orient the faces
  typedef CGAL::AABB_face_graph_triangle_primitive<Polyhedron> Primitive;
  typedef CGAL::AABB_traits<K, Primitive> Tree_traits;
  CGAL::AABB_tree<Tree_traits> tree(poly_ref.faces().begin(),
                                    poly_ref.faces().end(),
                                    poly_ref);
  tree.build();
  tree.accelerate_distance_queries();
  auto process_facet =
    [&points_in, &tree, &poly_ref](auto& facet) mutable {
      static std::atomic<int> counter;
      const auto& p0 = points_in[facet[0]];
      const auto& p1 = points_in[facet[1]];
      const auto& p2 = points_in[facet[2]];
      const K::Point_3 mid = CGAL::centroid(p0, p1, p2);
      std::pair<K::Point_3, Polyhedron::Face_index> pt_and_f =
        tree.closest_point_and_primitive(mid);
      auto face_ref_normal = get_face_normal(pt_and_f.second, poly_ref);
      if(face_ref_normal * cross_product(p1-p0, p2-p0) < 0) {
        std::swap(facet[1], facet[2]);
      }
      if((++counter % 100) == 0) {
        std::cerr << ".";
      }
    };
  
  auto process_facets =
    [&process_facet, &faces_in](const int& index) {
      process_facet(faces_in[index]);
    };
  tbb::parallel_for(0, int(faces_in.size()), process_facets);
  // std::for_each(faces_in.begin(), faces_in.end(), process_facet);
  std::cerr << "\n";

  std::ofstream soup_out(argv[3]);
  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;
  // std::cerr << std::boolalpha << PMP::is_polygon_soup_a_polygon_mesh(faces_in)
  //           << '\n';
  // PMP::polygon_soup_to_polygon_mesh(points_in, faces_in, poly_out);
  // std::ofstream out(argv[3]);
  // out.precision(17);
  // out << poly_out;
}