Error when merging triangle soup to Polyhedron_3 mesh

[gy114]

Hi, I collected the triangles using marching tetrahedra algorithm, searching for the triangular pieces using the unique_cell_iterator in a P3DT3.

Here's the triangle soup:

After I call
CGAL::Polygon_mesh_processing::merge_duplicate_points_in_polygon_soup(points, triangles);
CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points, triangles, isoSurface);

I got this result:

It seems that the kernel I'm using (Epick) is removing some points that are actually not duplicates.
Is there an easy fix for this? Specifically,

(1) Can I manually set the tolerance? How should I measure such tolerance?

(2) Should I switch to Epeck?

Thanks.

[sloriot]

Try calling PMP::orient_polygon_soup() before polygon_soup_to_polygon_mesh().

[gy114]

Thanks. I'll try that.

[gy114]

I tried orienting triangles before merging to a polygon mesh. It worked. Here's the outcome:

[MaelRL]

Did you have an issue with using the function facets_in_complex_3_to_triangle_mesh() mentioned in the other discussion (#7904)?

[gy114]

Thanks. It actually worked perfectly for an orthogonal cuboid. I then tried to find a way to modify the mesh after its created, such that it could have better quality for a non-orthogonal lattice. I could not access the details of how parameters such as params::facet_angle(30) are evaluated. After transforming (x,y,z) -> (a1x + a2y + a3*z), the geometry such as facet_angle and facet_size etc. should be changed. Is there an easy way to modify how these criteria are evaluated during the mesh generation?

P.S. Here's the code that worked, and the image is the visualization of the *.off output file.

int testStandardPeriodicMeshGenerationAndWriteToVTP() {
    // Define the cube that bounds the periodic domain
    Iso_cuboid domain_bounds(P3_Point(-0.5, -0.5, -0.5), P3_Point(0.5, 0.5, 0.5));
    // Create domain using the periodic function directly
    //Periodic_mesh_domain domain = Periodic_mesh_domain(periodic_function, domain_bounds);
    Periodic_mesh_domain domain = Periodic_mesh_domain::create_implicit_mesh_domain(periodic_function, domain_bounds);
    Periodic_mesh_criteria criteria(params::facet_angle(30)
                                            .facet_size(0.05 )
                                            .facet_distance(0.025 )
                                            .cell_radius_edge_ratio(2.)
                                            .cell_size(0.05));
    C3t3 c3t3_bis = CGAL::make_periodic_3_mesh_3<C3t3>(domain, criteria,  params::no_odt().
                                                                             no_lloyd().
                                                                             no_perturb().
                                                                             no_exude());
    // Output
    //Periodic_mesh_domain::create
    std::ofstream outfile1("outfile1.off");
    //c3t3_bis.output_facets_in_complex_to_off(outfile);
    c3t3_bis.output_facets_in_complex_to_off(outfile1);
    std::cout << "Mesh has " << c3t3_bis.number_of_cells() << " cells." << std::endl;
    outfile1.close();
    Surface_mesh tmesh;
    CGAL::facets_in_complex_3_to_triangle_mesh(c3t3_bis, tmesh);
    write_off(tmesh,"outputSurface.off");
    write_to_vtp_SurfaceMesh("outputSurface.vtp", tmesh);
    return 0;
}

[MaelRL]

You could try something like:

• extract the C3T3 surface
• Scale it
• Use PMP::isotropic_remeshing and 1. Constrain boundary edges 2. Provide a projection functor (see the named parameters) that projects onto your final isosurface

[gy114]

Thanks. I'll try that.
I will post my results here.

[gy114]

Hi,

I was running the following code and it seemed to went to an infinite loop. I did call "split_long_edges". See down below.
Also, how does the projector functor work? I could not find an example on the repository. Should it take a point and return a point? If I keep the default value, does it mean that the projection is already to the scaled surface? Thanks.

#include <vtkSmartPointer.h>
#include <vtkUnstructuredGrid.h>
#include <vtkDoubleArray.h>
#include <vtkPoints.h>
#include <vtkCellData.h>
#include <vtkPointData.h>
#include <vtkPolyData.h>
#include <vtkPolygon.h>
#include <vtkXMLPolyDataWriter.h>


#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Periodic_3_triangulation_ds_vertex_base_3.h>
#include <CGAL/Periodic_3_triangulation_ds_cell_base_3.h>
#include <CGAL/Triangulation_vertex_base_with_info_3.h>
#include <CGAL/Triangulation_cell_base_with_info_3.h>
#include <CGAL/Periodic_3_Delaunay_triangulation_3.h>
#include <CGAL/Triangulation_data_structure_3.h>
#include <CGAL/make_periodic_3_mesh_3.h>
#include <CGAL/Periodic_3_mesh_triangulation_3.h>
#include <CGAL/Labeled_mesh_domain_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/facets_in_complex_3_to_triangle_mesh.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/Surface_mesh.h>

#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
#include <CGAL/Polygon_mesh_processing/repair_polygon_soup.h>
#include <CGAL/Polygon_mesh_processing/remesh.h>



#include <CGAL/number_type_config.h> // CGAL_PI

#include <cmath>
#include <random>
#include <iostream>
#include <unordered_set>
#include <armadillo>




typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Periodic_3_Delaunay_triangulation_traits_3<K> P3_Traits;
typedef CGAL::Periodic_3_triangulation_ds_vertex_base_3<> P3_VbDS;
typedef CGAL::Triangulation_vertex_base_3<P3_Traits, P3_VbDS >       P3_Vb;
typedef CGAL::Periodic_3_triangulation_ds_cell_base_3<> P3_CbDS;
typedef CGAL::Triangulation_cell_base_3<P3_Traits, P3_CbDS >       P3_Cb;


typedef CGAL::Triangulation_vertex_base_with_info_3<double, P3_Traits, P3_Vb > VbInfo; // Test to attach data to vertices.
typedef CGAL::Triangulation_cell_base_with_info_3<double, P3_Traits, P3_Cb > CbInfo; // Test to attach data to cells.
typedef CGAL::Triangulation_data_structure_3<VbInfo, CbInfo >                 TDS;
typedef CGAL::Periodic_3_Delaunay_triangulation_3<P3_Traits, TDS >               P3_Delaunay;
typedef P3_Delaunay::Point P3_Point;
typedef P3_Delaunay::Iso_cuboid Iso_cuboid;
typedef CGAL::Surface_mesh<P3_Point> Surface_mesh;


typedef CGAL::Labeled_mesh_domain_3<K>                      Periodic_mesh_domain;
// Triangulation
typedef CGAL::Periodic_3_mesh_triangulation_3<Periodic_mesh_domain>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr>                       C3t3;
// Criteria
typedef CGAL::Mesh_criteria_3<Tr>                           Periodic_mesh_criteria;
// To avoid verbose function and named parameters call
namespace params = CGAL::parameters;



void write_to_vtp_SurfaceMesh(const std::string& filename, const Surface_mesh& mesh) {
    vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
    vtkSmartPointer<vtkPolyData> polyData = vtkSmartPointer<vtkPolyData>::New();
    vtkSmartPointer<vtkCellArray> polygons = vtkSmartPointer<vtkCellArray>::New();

    // Transfer vertices to VTK points
    for (auto v : mesh.vertices()) {
        const P3_Point& p = mesh.point(v);
        points->InsertNextPoint(p.x(), p.y(), p.z());
    }

    // Transfer faces to VTK polygons
    for (auto f : mesh.faces()) {
        vtkSmartPointer<vtkPolygon> polygon = vtkSmartPointer<vtkPolygon>::New();
        polygon->GetPointIds()->SetNumberOfIds(mesh.degree(f));
        // Add vertex indices to the polygon
        size_t i = 0;
        for (auto v : vertices_around_face(mesh.halfedge(f), mesh)) {
            polygon->GetPointIds()->SetId(i++, v.idx());
        }
        polygons->InsertNextCell(polygon);
    }

    polyData->SetPoints(points);
    polyData->SetPolys(polygons);

    // Write to VTK file
    vtkSmartPointer<vtkXMLPolyDataWriter> writer = vtkSmartPointer<vtkXMLPolyDataWriter>::New();
    writer->SetFileName(filename.c_str());
    writer->SetInputData(polyData);
    //writer->SetDataModeToBinary();
    writer->SetDataModeToAscii();
    writer->Write();
}

double periodic_function(const P3_Point& p) {
    return -std::cos(2 * CGAL_PI * p[0]) - std::cos(2 * CGAL_PI * p[1]) - std::cos(2 * CGAL_PI * p[2]);
}


struct halfedge2edge
{
    halfedge2edge(const Surface_mesh& m, std::vector<boost::graph_traits<Surface_mesh>::edge_descriptor>& edges)
            : m_mesh(m), m_edges(edges)
    {}
    void operator()(const boost::graph_traits<Surface_mesh>::halfedge_descriptor& h) const
    {
        m_edges.push_back(edge(h, m_mesh));
    }
    const Surface_mesh& m_mesh;
    std::vector<boost::graph_traits<Surface_mesh>::edge_descriptor>& m_edges;
};


int testStandardPeriodicMeshGenerationAndWriteToVTP() {
    // Define the cube that bounds the periodic domain
    Iso_cuboid domain_bounds(P3_Point(-0.5, -0.5, -0.5), P3_Point(0.5, 0.5, 0.5));

    // Create domain using the periodic function directly
    Periodic_mesh_domain domain = Periodic_mesh_domain::create_implicit_mesh_domain(periodic_function, domain_bounds);
    Periodic_mesh_criteria criteria(params::facet_angle(30)
                                            .facet_size(0.05 )
                                            .facet_distance(0.025 )
                                            .cell_radius_edge_ratio(2.)
                                            .cell_size(0.05));
    C3t3 c3t3_bis = CGAL::make_periodic_3_mesh_3<C3t3>(domain, criteria,  params::no_odt().
            no_lloyd().
            no_perturb().
            no_exude());
    // Output
    std::ofstream outfile1("outfile1.off");
    c3t3_bis.output_facets_in_complex_to_off(outfile1);
    std::cout << "Mesh has " << c3t3_bis.number_of_cells() << " cells." << std::endl;
    outfile1.close();
    Surface_mesh tmesh;
    CGAL::facets_in_complex_3_to_triangle_mesh(c3t3_bis, tmesh);
    //We now scale the old points to the non-orthogonal unit cell.
    arma::mat B = {{1.5064,  -0.8697,   0},
                   {1.5064,   0.8697,   0},
                   {0,        0,   0.9398}};

    for (auto v : tmesh.vertices())  {
        P3_Point& p = tmesh.point(v);
        arma::vec armaP = {p.x(), p.y(), p.z()};
        arma::vec scaledArmaP = (B.t())*armaP;
        p = P3_Point(scaledArmaP(0), scaledArmaP(1), scaledArmaP(2));
    }
    double target_edge_length = 0.04;
    std::cout << "Split border...";
    std::vector<boost::graph_traits<Surface_mesh>::edge_descriptor> border;
    CGAL::Polygon_mesh_processing::border_halfedges(faces(tmesh),tmesh, boost::make_function_output_iterator(halfedge2edge(tmesh, border)));
    CGAL::Polygon_mesh_processing::split_long_edges(border, target_edge_length, tmesh);
    write_to_vtp_SurfaceMesh("outputSurfaceBorderSplit.vtp", tmesh);
    std::cout << "done." << std::endl;


    std::cout << "Start remeshing... "
              << " (" << num_faces(tmesh) << " faces)...";
    CGAL::Polygon_mesh_processing::isotropic_remeshing(faces(tmesh), target_edge_length, tmesh, params::number_of_iterations(3)
                                                                                                        .protect_constraints(true));
    std::cout << "done." << std::endl;
    write_to_vtp_SurfaceMesh("outputSurface.vtp", tmesh);
    return 0;
}

int main() {
    testStandardPeriodicMeshGenerationAndWriteToVTP();
    return 0;
}

This is the how "outputSurfaceBorderSplit.vtp" looks.