smooth subdomain separation from labelled images

[juanSilvaHenao]

Issue Details

Greetings,

I am using mesh_3D_weighted_image_with_detection_of_features.cpp to generate a 3D mesh from a segmented image with several subdomains. The code works without any compilation issues however, the generated meshes seem to have problems separating the different subdomains like it is shown in the image

I am using feature protection both inside the image and in bbox and also tried using weights. I also tried tuning the facet parameters (reducing facet size and distance) as well as edge length and cell size but nothing does the trick. Is there any way I can guarantee that there is no cell sharing between certain subdomains? Also, is there a way to force two adjacent subdomain to hace separate nodes in the boudary (e.g. for contact simulation)?

I send attached the image file I am trying to mesh as well as my source code, wich is a modified version of mesh_3D_weighted_image_with_detection_of_features.cpp

case_6.zip

Thanks in advance

Source Code

#define CGAL_MESH_3_WEIGHTED_IMAGES_DEBUG
#define CGAL_MESH_3_VERBOSE 1

#include <vtkNew.h>
#include <vtkImageData.h>
#include <vtkNIFTIImageReader.h>
#include <vtkImageReader.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/Labeled_mesh_domain_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/Image_3.h>
#include <CGAL/IO/read_vtk_image_data.h>
#include <CGAL/IO/output_to_vtu.h>
#include <boost/lexical_cast.hpp>
#include <boost/functional.hpp>
#include <filesystem>
#include "mesh_3D_gray_vtk_image.h"
#include <CGAL/Mesh_domain_with_polyline_features_3.h>
#include <CGAL/Mesh_3/Detect_features_in_image.h>
#include <CGAL/Mesh_3/generate_label_weights.h>
#include <CGAL/Mesh_constant_domain_field_3.h>

//typedef short Image_word_type;

// Domain

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;

#ifdef CGAL_CONCURRENT_MESH_3
typedef CGAL::Parallel_tag Concurrency_tag;
#else
typedef CGAL::Sequential_tag Concurrency_tag;
#endif

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

// Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
typedef CGAL::Mesh_constant_domain_field_3<Mesh_domain::R,
    Mesh_domain::Index> Sizing_field;


namespace params = CGAL::parameters;



namespace fs = std::filesystem;

int main(int argc, char* argv[])
{
    // Loads image
    if (argc < 5) {
        std::cerr << "Usage: " << argv[0] << " <nii in file> <vtu out file> <cell_size> <element_ratio>"
            << std::endl;
        return 1;
    }

    vtkImageData* vtk_image = nullptr;

    fs::path path(argv[1]);

    vtkNIFTIImageReader* reader = vtkNIFTIImageReader::New();
    reader->SetFileName(argv[1]);
    reader->Update();
    vtk_image = reader->GetOutput();
    vtk_image->Print(std::cerr);
    std::cout << ".nii image loaded";

    if (vtk_image == nullptr) {
        std::cout << "No image loaded" << std::endl;
        return 0;
    }
    CGAL::Image_3 image = CGAL::IO::read_vtk_image_data(vtk_image);
    if (image.image() == nullptr) {
        std::cerr << "could not create a CGAL::Image_3 from the vtk image\n";
        return 0;
    }


    /// [Domain creation]
    const float sigma = (std::max)(image.vx(), (std::max)(image.vy(), image.vz()));
    CGAL::Image_3 img_weights = CGAL::Mesh_3::generate_label_weights(image, sigma);

    Mesh_domain domain = Mesh_domain::create_labeled_image_mesh_domain(image,
        params::weights = std::ref(img_weights),
        params::value_outside(0.f),
        params::relative_error_bound(1e-10),
        params::features_detector = CGAL::Mesh_3::Detect_features_in_image());
  
    double cell_size = boost::lexical_cast<double>(argv[3]);
    double facet_angle =30;
    double facet_distance = 0.00001 * cell_size;
    double element_ratio = boost::lexical_cast<double>(argv[4]);;
    double min_size = 0.25 * cell_size;  // element_ratio* cell_size; (used for testing only)
    int volume_dimension = 3;



    // Mesh criteria
    Mesh_criteria criteria(params::facet_angle(facet_angle).
        non_manifold().
        facet_size(cell_size).
        facet_min_size(min_size).
        facet_distance(facet_distance).
        edge_size(cell_size).
        edge_min_size(min_size).
        cell_radius_edge_ratio(3).
        cell_size(cell_size).
        cell_min_size(min_size));


    // Meshing
    C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria,
        params::lloyd(params::time_limit(30)).
        no_perturb().
        no_exude());
    // Output
    std::ofstream file(argv[2]);
    CGAL::IO::output_to_vtu(file, c3t3, CGAL::IO::ASCII);

    return 0;

Environment

• Operating system (Windows/Mac/Linux, 32/64 bits): Windows 64 bits
• Compiler: visual studio
• Release or debug mode: release
• Specific flags used (if any):
• CGAL version: 5.6

[sloriot]

could you provide the parameters to pass to the command line?

[juanSilvaHenao]

Yes:
it takes the .nii image file, the vtu out file, cell size and element_ratio (this value doesn't matter as it is fixed in the code to 0.25*cell_size), so with the file provided it will be:

case_6.nii case_6.vtu 2 0.5

[janetournois]

Can you please tell us where this unexpected connection happens (coordinates of the points), so that we can spot it too in the input image and output mesh?
My guess for now is that you should use a smaller facet_distance criterion

[juanSilvaHenao]

the coordinates of one of the points are:
(72.6996, 50.5809, 111.179)
it is one of the points in the bottom right of the attached image. I am already using a very small facet distance (0.00002) but the problem is still present, for example in regions of the domain where the two subdomains are close together but there should be space between them, for example around the node with coordinates:
(35.474, 113.694, 12.2524)

[janetournois]

Hello,

I could generate a tet mesh from your input with almost the same meshing criteria, and the volumes are nicely separated :

Actually, facet_distance does not have to be that small (I used 0.01 and it could probably be higher).
However, the minimal sizes edge_min_size, facet_min_size and cell_min_size that you were using were too large. These criteria have been designed to prevent any other meshing criteria from triggering refinement of simplices smaller than that, hence the issues that you are facing.

[juanSilvaHenao]

Thank you! I tried once again removing the min_size parameters and it works with a facet distance of 0.08 :). I started using these parameters (min_sizes) to speed up the meshing process without realizing that it will induce this error.

One final question, is there a way to indicate that two adjacent subdomains should be touching but not sharing any nodes?

This can be useful to simulate contact between two bodies without meshing them separately (which I assume is also an option).

Thanks for your help!

[janetournois]

Mesh_3 meshes all subdomains simultaneously and ensures that connectivity between subdomains is valid, with subdomains separated by surface facets.
You can separate the subdomains a posteriori and duplicate surface vertices. Would it help?

[juanSilvaHenao]

I think that can be what I need. Is there any example or anywhere I can get information that I can use as a guide for doing this?
Sorry, I am fairly new to CGAL so that might be a very basic question....

[janetournois]

There isn't really an example of how to do that.

As an intermediate step, the function CGAL::facets_in_complex_3_to_triangle_mesh() can be used to export the surface mesh corresponding to a given surface patch index

[juanSilvaHenao]

Ok, i'll try starting with it. Thank yo so much for your help!