Test cgal's fair() function using a cube mesh with a hole

main.cpp

/* adapted from cgal's example
** The purpose is to test fair() and expose issues,
** e.g. the faired region contains bad vertices (NaN)
*/

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polygon_mesh_processing/fair.h>
#include <CGAL/Polygon_mesh_processing/random_perturbation.h>
#include <CGAL/Polygon_mesh_processing/refine.h>
#include <CGAL/Polygon_mesh_processing/triangulate_hole.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/utility.h>
#include <fstream>
#include <map>

#include <glog/logging.h>

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_3 CGALPoint3;
typedef K::Vector_3 CGALVec3;
typedef CGAL::Surface_mesh<CGALPoint3> CGALMesh;
typedef CGALMesh::Vertex_index VHandle;
typedef CGALMesh::Face_index FHandle;
typedef CGALMesh::Edge_index EHandle;
typedef CGALMesh::Halfedge_index HEHandle;

bool isColinear(
    const CGALPoint3& p0,
    const CGALPoint3& p1,
    const CGALPoint3& p2,
    double eps /*  = 1e-9 */) {
  CGALVec3 v1 = p0 - p1;
  CGALVec3 v2 = p0 - p2;
  CGALVec3 v3 = p1 - p2;
  if (v1.squared_length() < eps || v2.squared_length() < eps || v3.squared_length() < eps) {
    return true;
  }

  v2 /= CGAL::sqrt(v2.squared_length());
  v1 /= CGAL::sqrt(v1.squared_length());
  auto projLen = std::abs(v2 * v1);
  if (projLen > 1.0 - eps) {
    return true;
  }
  return false;
}

bool checkThinFaces(const CGALMesh& mesh) {
  bool has_thinTris = false;
  for (const auto fi_cgal : mesh.faces()) {
    std::vector<CGALPoint3> f;
    for (const auto vidOfFace : CGAL::vertices_around_face(mesh.halfedge(fi_cgal), mesh)) {
      f.push_back(mesh.point(vidOfFace));
    }
    if (isColinear(f[0], f[1], f[2], 1e-9)) {
      LOG(WARNING) << "Thin triangles detected! DO NOT apply fairing to this hole.";
      has_thinTris = true;
    }
  }
  return has_thinTris;
}

bool edgeExistsInMesh(VHandle v1, VHandle v2, const CGALMesh& m) {
  return m.halfedge(v1, v2) != m.null_halfedge() || m.halfedge(v2, v1) != m.null_halfedge();
}

int main(int argc, char* argv[]) {
  const char* filename = (argc > 1) ? argv[1] : "data/blobby.off";
  std::ifstream input(filename);
  CGALMesh mesh;
  if (!input) {
    std::cerr << "Not a valid input file." << std::endl;
    return 1;
  }
  if (!(input >> mesh) || mesh.is_empty() || !CGAL::is_triangle_mesh(mesh)) {
    std::cerr << "Cannot produce a valid CGAL mesh!" << std::endl;
    return 2;
  }

  /* find a hole boundary */
  HEHandle border_hei = CGALMesh::null_halfedge();
  for (HEHandle hei : mesh.halfedges()) {
    if (mesh.is_border(hei)) {
      border_hei = hei;
      break;
    }
  }

  if (border_hei == CGALMesh::null_halfedge()) {
    std::cerr << "Cannot find any border edge!" << std::endl;
    return 3;
  }

  HEHandle next_border_hei = border_hei;
  std::vector<VHandle> hole_vhdls;
  std::vector<CGALPoint3> hole_pts;
  std::vector<CGALPoint3> hole_3rd_pts;
  std::vector<FHandle> hole_faces;
  do {
    auto hole_vhdl = mesh.target(next_border_hei);
    hole_vhdls.push_back(hole_vhdl);
    hole_pts.push_back(mesh.point(hole_vhdl));

    HEHandle hei_border_face = mesh.opposite(next_border_hei);
    hole_faces.push_back(mesh.face(hei_border_face));
    VHandle otherVi_in_f = mesh.target(mesh.next(hei_border_face));
    hole_3rd_pts.push_back(mesh.point(otherVi_in_f));

    next_border_hei = mesh.next(next_border_hei);
  } while (next_border_hei != border_hei);
  CHECK_EQ(hole_3rd_pts.size(), hole_pts.size());

  // patch the hole
  std::vector<CGAL::Triple<int, int, int>> patch;
  CGAL::Polygon_mesh_processing::triangulate_hole_polyline(
      hole_pts, hole_3rd_pts, std::back_inserter(patch));
  LOG(INFO) << "Done: hole triangulated.";

  // two ways to apply the patch:
  // 1. lift the patch along normals and build a "scaffold surface" between lifted patch
  //    and the hole boundary. This avoids generating non-manifold edges as the patch faces
  //    might contain edges that already exist in the original mesh
  // 2. directly apply the patch, which could risk introducing non-manifold edges.
  //    In that case error will be reported
  std::vector<FHandle> fhdls_to_refine;
  std::vector<VHandle> vhdls_to_fair;
  bool perform_patch_lifting = true;
  if (perform_patch_lifting) {
    // determine lift length along some normal dir
    const auto f1 = patch[patch.size() / 2];
    const auto edge_len = CGAL::sqrt((hole_pts[f1.first] - hole_pts[f1.second]).squared_length());
    double lift_len = edge_len * 0.1;

    // lift hole points & add to surfacemesh
    std::vector<VHandle> lift_vhdls;
    for (int i = 0; i < hole_pts.size(); ++i) {
      std::vector<VHandle> f_vhdls;
      for (const auto vh : mesh.vertices_around_face(mesh.halfedge(hole_faces[i]))) {
        f_vhdls.push_back(vh);
      }

      const auto lift_vec =
          CGAL::unit_normal(
              mesh.point(f_vhdls[0]), mesh.point(f_vhdls[1]), mesh.point(f_vhdls[2])) *
          lift_len;

      CGALPoint3 p = hole_pts[i];
      p = p + lift_vec;
      auto hdl = mesh.add_vertex(p);
      CHECK(hdl != mesh.null_vertex()) << "Cannot add lifted hole vertex!";
      lift_vhdls.push_back(hdl);
    }
    LOG(INFO) << "Done: hole points lifted.";

    // add a strip of triangles that connect original hole boundary with lifted hole boundary
    auto n_holePts = hole_vhdls.size();
    for (auto i = 0; i < n_holePts; ++i) {
      FHandle f_hdl;
      f_hdl = mesh.add_face(hole_vhdls[i], hole_vhdls[(i + 1) % n_holePts], lift_vhdls[i]);
      CHECK(f_hdl != mesh.null_face()) << "Cannot add 1st tri of the quad!";
      f_hdl = mesh.add_face(
          lift_vhdls[i], hole_vhdls[(i + 1) % n_holePts], lift_vhdls[(i + 1) % n_holePts]);
      CHECK(f_hdl != mesh.null_face()) << "Cannot add 2nd tri of the quad!";
    }
    LOG(INFO) << "Done: scaffold faces added.";

    // add the lifted hole patch faces
    for (const auto& pf : patch) {
      auto hdl = mesh.add_face(lift_vhdls[pf.first], lift_vhdls[pf.second], lift_vhdls[pf.third]);
      CHECK(hdl != mesh.null_face()) << "Cannot add lifted patch face!";
      fhdls_to_refine.push_back(hdl);
    }

    vhdls_to_fair = lift_vhdls;
  } else {
    // add the hole patch faces
    for (const auto& pf : patch) {
      auto v1 = hole_vhdls[pf.first];
      auto v2 = hole_vhdls[pf.second];
      auto v3 = hole_vhdls[pf.third];
      CHECK(!edgeExistsInMesh(v1, v2, mesh))
          << "Invalid! An edge of the patch face is already existing in the mesh with hole!";
      CHECK(!edgeExistsInMesh(v1, v3, mesh))
          << "Invalid! An edge of the patch face is already existing in the mesh with hole!";
      CHECK(!edgeExistsInMesh(v2, v3, mesh))
          << "Invalid! An edge of the patch face is already existing in the mesh with hole!";
      auto hdl = mesh.add_face(hole_vhdls[pf.first], hole_vhdls[pf.second], hole_vhdls[pf.third]);
      CHECK(hdl != mesh.null_face()) << "Cannot add hole patch face!";
      fhdls_to_refine.push_back(hdl);
    }
  }

  // check for thin triangles
  bool has_thinTris = checkThinFaces(mesh);
  if (has_thinTris) {
    LOG(WARNING) << "Mesh contains very thin triangles. This could cause bad fairing results!";
  }

  // if there are thin faces, apply perturbation to only vertices-to-fair
  if (has_thinTris) {
    const auto f1 = patch[patch.size() / 2];
    const auto edge_len = CGAL::sqrt((hole_pts[f1.first] - hole_pts[f1.second]).squared_length());
    CGAL::Polygon_mesh_processing::random_perturbation(
        vhdls_to_fair, mesh, edge_len * 0.1, CGAL::parameters::do_project(false));

    // check again for thin triangles
    bool still_has_thinTris = checkThinFaces(mesh);
    if (still_has_thinTris) {
      LOG(WARNING) << "Mesh still contains very thin triangles after perturbation!!!";
    } else {
      LOG(INFO) << "Mesh is free of thin triangles after perturbation~~~";
    }
  }

  // refine lifted patch
  std::vector<FHandle> new_faces_refined;
  std::vector<VHandle> new_vts_refined;
  CGAL::Polygon_mesh_processing::refine(
      mesh,
      fhdls_to_refine,
      std::back_inserter(new_faces_refined),
      std::back_inserter(new_vts_refined));

  // perform fairing on the lifted patch
  vhdls_to_fair.insert(vhdls_to_fair.end(), new_vts_refined.begin(), new_vts_refined.end());
  bool fair_suc = CGAL::Polygon_mesh_processing::fair(
      mesh, vhdls_to_fair, CGAL::parameters::fairing_continuity(2));
  if (!fair_suc) {
    LOG(WARNING) << "Failed to perform fairing on the final patch!";
  }

  std::ofstream faired_off("cgal_faired.off");
  faired_off << mesh;
  faired_off.close();

  return 0;
}