Mesh and remesh

mesh_and_remesh.cpp

#define CGAL_MESH_3_VERBOSE 1
#define CGAL_TETRAHEDRAL_REMESHING_VERBOSE

#include "mesh_from_labeled_3d_array_feature_detect.hpp"

#include <cassert>
#include <vector>
#include <iostream>

#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>

#include <CGAL/make_mesh_3.h>
#include <CGAL/Image_3.h>

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Mesh_domain_with_polyline_features_3.h>
#include <CGAL/Labeled_mesh_domain_3.h>

#include <CGAL/Mesh_3/Detect_features_in_image.h>

#include <CGAL/tetrahedral_remeshing.h>
#include <CGAL/Tetrahedral_remeshing/Remeshing_triangulation_3.h>
#include <CGAL/Adaptive_remeshing_sizing_field.h>

#include <CGAL/IO/File_medit.h>

#include <fstream>
#include <map>
#include <string>

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Tetrahedral_remeshing::Remeshing_triangulation_3<K> Remeshing_triangulation;

// Utility to read key=value parameter file
std::map<std::string, std::string> read_remeshing_parameters(const std::string& filename) {
  std::ifstream file(filename);
  std::map<std::string, std::string> params;
  std::string line;
  while (std::getline(file, line)) {
    // Skip comments and empty lines
    if (line.empty() || line[0] == '#') continue;
    std::istringstream iss(line);
    std::string key;
    if (std::getline(iss, key, '=')) {
      std::string value;
      if (std::getline(iss, value)) {
        key.erase(0, key.find_first_not_of(" \t"));
        key.erase(key.find_last_not_of(" \t") + 1);
        value.erase(0, value.find_first_not_of(" \t"));
        value.erase(value.find_last_not_of(" \t") + 1);
        params[key] = value;
      }
    }
  }
  return params;
}


typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Labeled_mesh_domain_3<K> Image_domain;
typedef CGAL::Mesh_domain_with_polyline_features_3<Image_domain> Mesh_domain;

// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;

// Mesh Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
#include <map>

std::map<std::string, std::string> read_parameters(const std::string& filename) {
    std::ifstream file(filename);
    std::map<std::string, std::string> params;
    std::string line;

    while (std::getline(file, line)) {
        std::istringstream iss(line);
        std::string key;
        if (std::getline(iss, key, '=')) {
            std::string value;
            if (std::getline(iss, value)) {
                key.erase(0, key.find_first_not_of(" \t"));
                key.erase(key.find_last_not_of(" \t") + 1);
                value.erase(0, value.find_first_not_of(" \t"));
                value.erase(value.find_last_not_of(" \t") + 1);
                params[key] = value;
            }
        }
    }

    return params;
}

bool to_bool(const std::string& s) {
    return s == "1" || s == "true" || s == "True";
}


C3t3
generate_from_inr_feature_detect(
    const std::string & inr_filename,
    const std::string & outfile,
    const bool lloyd,
    const bool odt,
    const bool perturb,
    const bool exude,
    const double max_edge_size_at_feature_edges,
    const double min_facet_angle,
    const double max_radius_surface_delaunay_ball,
    const double max_facet_distance,
    const double max_circumradius_edge_ratio,
    const double max_cell_circumradius,
    const double exude_time_limit,
    const double exude_sliver_bound,
    const bool verbose,
    const int seed
    )
{
  CGAL::get_default_random() = CGAL::Random(seed);

  CGAL::Image_3 image;
  const bool success = image.read(inr_filename.c_str());
  if (!success) {
    throw "Could not read image file";
  }

  Mesh_domain cgal_domain =
    Mesh_domain::create_labeled_image_mesh_domain(image, CGAL::parameters::features_detector = CGAL::Mesh_3::Detect_features_in_image());

  CGAL::Bbox_3 bbox = cgal_domain.bbox();
  double diag = CGAL::sqrt(CGAL::square(bbox.xmax() - bbox.xmin()) +
                           CGAL::square(bbox.ymax() - bbox.ymin()) +
                           CGAL::square(bbox.zmax() - bbox.zmin()));

  double sizing;
  if (max_edge_size_at_feature_edges < 0.001) {
    sizing = diag * 0.05;
  } else {
    sizing = max_edge_size_at_feature_edges;
  }

  Mesh_criteria criteria(
      CGAL::parameters::edge_size=sizing,
      CGAL::parameters::facet_angle=min_facet_angle,
      CGAL::parameters::facet_size=max_radius_surface_delaunay_ball,
      CGAL::parameters::facet_distance=max_facet_distance,
      CGAL::parameters::facet_topology = CGAL::FACET_VERTICES_ON_SAME_SURFACE_PATCH,
      CGAL::parameters::cell_radius_edge_ratio=max_circumradius_edge_ratio,
      CGAL::parameters::cell_size=max_cell_circumradius
      );

  // Mesh generation
  if (!verbose) {
    // suppress output
    std::cerr.setstate(std::ios_base::failbit);
  }
  return CGAL::make_mesh_3<C3t3>(
      cgal_domain,
      criteria,
      lloyd ? CGAL::parameters::lloyd(CGAL::parameters::time_limit(0)) : CGAL::parameters::no_lloyd(),
      odt ? CGAL::parameters::odt(CGAL::parameters::time_limit(0)) : CGAL::parameters::no_odt(),
      perturb ? CGAL::parameters::perturb() : CGAL::parameters::no_perturb(),
      exude ?
        CGAL::parameters::exude(
          CGAL::parameters::time_limit = exude_time_limit,
          CGAL::parameters::sliver_bound = exude_sliver_bound
        ) :
        CGAL::parameters::no_exude()
      );
}

int main(int argc, char* argv[])
{
    // std::string mode = "image";
    // const char *modes = &mode[0];
    if (argc < 3) {
        std::cerr << "Usage: " << argv[0] << " input.txt input_remesh.txt \n";
        return 1;
    }

    std::map<std::string, std::string> p = read_parameters(argv[1]);

    const std::string inr_filename = p["inr_filename"];
    const std::string outfile = p["outfile"];
    const bool lloyd = to_bool(p["lloyd"]);
    const bool odt = to_bool(p["odt"]);
    const bool perturb = to_bool(p["perturb"]);
    const bool exude = to_bool(p["exude"]);
    const float sig = std::stof(p["sig"]);
    const double relative_error_bound_weights = std::stod(p["relative_error_bound_weights"]);
    const double max_edge_size_at_feature_edges = std::stod(p["max_edge_size_at_feature_edges"]);
    const double min_facet_angle = std::stod(p["min_facet_angle"]);
    const double max_radius_surface_delaunay_ball = std::stod(p["max_radius_surface_delaunay_ball"]);
    const double max_facet_distance = std::stod(p["max_facet_distance"]);
    const double max_circumradius_edge_ratio = std::stod(p["max_circumradius_edge_ratio"]);
    const double max_cell_circumradius = std::stod(p["max_cell_circumradius"]);
    const double exude_time_limit = std::stod(p["exude_time_limit"]);
    const double exude_sliver_bound = std::stod(p["exude_sliver_bound"]);
    const bool verbose = to_bool(p["verbose"]);
    const int seed = std::stoi(p["seed"]);

    C3t3 c3t3 = generate_from_inr_feature_detect(inr_filename,
                        outfile,
                        lloyd,
                        odt,
                        perturb,
                        exude,
                        max_edge_size_at_feature_edges,
                        min_facet_angle,
                        max_radius_surface_delaunay_ball,
                        max_facet_distance,
                        max_circumradius_edge_ratio,
                        max_cell_circumradius,
                        exude_time_limit,
                        exude_sliver_bound,
                        verbose,
                        seed);

    assert(c3t3.triangulation().is_valid(true));

    // Output
    std::ofstream medit_file(outfile);
    c3t3.output_to_medit(medit_file);
    medit_file.close();

    CGAL::dump_c3t3(c3t3, "c3t3_dump");


    // REMESHING

    std::map<std::string, std::string> rp = read_remeshing_parameters(argv[1]);

    const std::string input_filename = rp.count("input_mesh") ? rp["input_mesh"] : "data/sphere.mesh";
    const std::string output_filename = rp.count("output_mesh") ? rp["output_mesh"] : "after_remeshing.mesh";
    const double target_edge_length = rp.count("target_edge_length") ? std::stod(rp["target_edge_length"]) : 0.1;
    const int num_iterations = rp.count("num_iterations") ? std::stoi(rp["num_iterations"]) : 5;
    const bool use_adaptive = rp.count("adaptive") ? (rp["adaptive"] == "1" || rp["adaptive"] == "true") : true;

    auto tr = CGAL::convert_to_triangulation_3(std::move(c3t3));

    if (!tr.is_valid(true))
    {
      std::cerr << "Invalid input mesh.\n";
      return EXIT_FAILURE;
    }

    std::cout << "Remeshing..." << std::endl;

    if (use_adaptive) {
      CGAL::tetrahedral_isotropic_remeshing(
        tr,
        CGAL::create_adaptive_remeshing_sizing_field(tr),
        CGAL::parameters::number_of_iterations(num_iterations)
      );
    }
    else {
      CGAL::tetrahedral_isotropic_remeshing(
        tr,
        target_edge_length,
        CGAL::parameters::number_of_iterations(num_iterations)
      );
    }

    std::ofstream os(output_filename);
    if (!os) {
      std::cerr << "Error opening output file: " << output_filename << "\n";
      return EXIT_FAILURE;
    }

    CGAL::IO::write_MEDIT(os, tr);

    std::cout << "Output written to " << output_filename << std::endl;


    return 0;
}