Spikes after make mesh

[TorgeW]

Hey,

We are using cgal to created meshes from labeled head images.

We are using the following code:

struct Sizing_field
{
    typedef ::FT FT;
    typedef Point Point_3;
    typedef Mesh_domain_img::Index Index; 
    double tx, ty, tz;
    double vx, vy, vz;
    std::size_t sx, sy, sz;
    float *sizing_field_image;

    FT operator()(const Point_3 &p, const int, const Index&) const
    {
        // I add 0.5 offset because CGAL assumes centered voxels
        const std::size_t i = (p.x() - tx)/vx + 0.5;
        const std::size_t j = (p.y() - ty)/vy + 0.5;
        const std::size_t k = (p.z() - tz)/vz + 0.5;
        if (i < 0 || j < 0 || k < 0 || i >= sx || j >= sy || k >= sz) {
            std::cerr << "trying to access sizing field out-of-bounds" << std::endl;
            return 1;
        }
	      FT val = sizing_field_image[i + sx*j + sx*sy*k];
        return val;
    };
};

int _mesh_image_sizing_field(
  char *fn_image, char *fn_out,
  float facet_angle, float *facet_size, float *facet_distance,
  float cell_radius_edge_ratio, float *cell_size,
  bool optimize, int num_threads
)
{
  /// Load image
  CGAL::Image_3 image;
  if(!image.read(fn_image)){
    std::cerr << "Error: Cannot read file " <<  fn_image << std::endl;
    return EXIT_FAILURE;
  }

  // Set max number of threads
  tbb::global_control global_limit(tbb::global_control::max_allowed_parallelism, num_threads);

  // Mesh domain
  Mesh_domain_img domain = Mesh_domain_img::create_labeled_image_mesh_domain(image, 1e-10);

  // Mesh criteria
  Sizing_field sizing_field_cell = {
     image.tx(), image.ty(), image.tz(),
     image.vx(), image.vy(), image.vz(),
     image.xdim(), image.ydim(), image.zdim(),
     cell_size
  };
  Sizing_field sizing_field_facet_distance = {
     image.tx(), image.ty(), image.tz(),
     image.vx(), image.vy(), image.vz(),
     image.xdim(), image.ydim(), image.zdim(),
     facet_distance
  };
  Sizing_field sizing_field_facet_size = {
     image.tx(), image.ty(), image.tz(),
     image.vx(), image.vy(), image.vz(),
     image.xdim(), image.ydim(), image.zdim(),
     facet_size
  };
  Facet_criteria_img facet_criteria(
          facet_angle,
          sizing_field_facet_size,
          sizing_field_facet_distance
  );
  Cell_criteria_img cell_criteria(
          cell_radius_edge_ratio,
          sizing_field_cell
  );
  Mesh_criteria_img criteria(facet_criteria, cell_criteria);
  
  // Mesh generation
  std::cout << "Began meshing \n";
  C3t3_img c3t3 = CGAL::make_mesh_3<C3t3_img>(domain, criteria, no_perturb(), no_exude());

  // Output
  std::ofstream medit_file(fn_out);
  c3t3.output_to_medit(medit_file);
  return EXIT_SUCCESS;
}

This is the result where you can clearly see the spikes.

We are currently relabeling / splitting the spikes to resolve this issue but maybe there is a better solution?

[lrineau]

That is not expected. I suggest to decrease the value of facet_distance to something significantly lower than the the width of a voxel (min of {vx, vy, vz}).

Would you mind sharing the 3D image? Maybe there is noise in it?

[TorgeW]

Thanks for the fast reply!

There is no noise in the 3D image, its a label image.
final_tissues.nii.gz

I will try your suggestion and come back after! :)

[TorgeW]

So decreasing the facet_distance indeed reduces the amount of spiking but it also opens up new issues:
There are almost 5 times as many tetrahedra and we can see an increased staircase effect on the surfaces (Especially the skin).
Keeping the element count low and receiving a somewhat smooth but also accurate mesh is of cause crucial.

[lrineau]

Can you please give use the parameter you use in _mesh_image_sizing_field:

  float facet_angle, float *facet_size, float *facet_distance,
  float cell_radius_edge_ratio, float *cell_size,

I would like to reproduce and experiment with your input data file.

[TorgeW]

see other comment

[TorgeW]

final_tissues.nii.gz
facet_distance: distance_field.nii.gz
facet_size, cell_size: https://drive.google.com/drive/folders/1fNgNHeDCAnyHimkDlNt2gEmOQ1D4SXoP?usp=sharinghttps://drive.google.com/drive/folders/1fNgNHeDCAnyHimkDlNt2gEmOQ1D4SXoP?usp=sharing

facet_angle=30,
cell_radius_edge_ratio=3

If you need anything else, feel free to ask!

[TorgeW]

@lrineau any updates on this? We would really appreciate a solution that would not generate spikes but still generate a mesh with a course surface true to the specification in the facet_distance.

[janetournois]

Hello @TorgeW
I am working on your issue, and I noticed that the 3 images (final_tissues, distance_field and facet_size/cell_size) do not have the same bounding boxes. Is it what you expect?

[TorgeW]

If the bounding box is defined by having a box around the data, excluding the outside (0 for final_tissues, 2 for distance_field and 20 for cell_size) than a small difference is to be expected. The cell_size is smoothed so it will be slightly bigger and also has a 1 voxel thick shell on the outside of the image (an artifact from processing). The distance field might also not have the exact same bounding box as tissue_labeling_upsampled.nii.gz. All of this shouldn't really be an issue.

One thing, i made a slight mistake with this data so here is a valid dataset
( The label image was not 100% correct, but again this doesn't change the underlying problem)
https://drive.google.com/drive/folders/1fNgNHeDCAnyHimkDlNt2gEmOQ1D4SXoP?usp=sharinghttps://drive.google.com/drive/folders/1fNgNHeDCAnyHimkDlNt2gEmOQ1D4SXoP?usp=sharing

[janetournois]

There is something I still don't understand.
The 3 images do not have the same bounding boxes. Should size_field and distance_field be translated by some vector, to have the images aligned?

[TorgeW]

You mean the image size, yes.
I uploaded the uncropped label image, corrected it to the cropped version, I am sorry for the confusion.

Now all 3 images should have the same size.

[janetournois]

Hello @TorgeW

thank you for the new data.

I can reproduce the spikes in the regions describing a subdomain that is only one-voxel-thick at most. To mesh very thin regions like that and avoid spikes, the distance field should return a distance which is at most half of the thickness of the domain.

The 3 straightforward ways for separating surfaces are driven by the Mesh_criteria_3 :

• small approximation error to ensure separation by choosing a small enough distance bound to the input surface
• small facet size that will automatically reduce surface approximation error
• small cell size to ensure that there are at least two "layers" of tetrahedra in each mesh subdomain

All of them are local criteria, that can be defined as you did by a sizing field.

I hope this helps!

[TorgeW]

Hey @janetournois

Thanks for your answer.

Unfortunately spikes also occur in regions where describing a subdomain that is more than one voxel thick (see images), although they are more frequent in thin regions. They are usually occurring regions that are one tetrahedron thick but there are also spikes in regions thicker than one tetrahedron.

All 3 ways you are mentioning here will increase the surface triangle density (and therefore also the tetrahedron density) in the mesh, which is exactly what we want to keep as low as possible. Our current solution involves relabeling (based on the barycenter of the tetrahedron in the label image) and splitting of spikes which brings us pretty close to a good result with a course mesh but its a complex process that can miss spikes.

[janetournois]

Hello @TorgeW

please have a look at the pull-request #7754.
It introduces the possibility to use the manifold() criterion, that fixes the non-manifold edges and vertices on surfaces, together with the facet_min_size criterion, that prevents surface facets from over-refinement when another criterion (like manifoldness) cannot be achieved, for example because the input itself is not manifold.

If the output mesh is too dense for your needs, you can run CGAL::tetrahedral_isotropic_remeshing() that will keep the validity of surfaces and un-refine them where it is required.

This combination of parameters and options should solve your issues.

[TorgeW]

Thanks for your help! :)
We will check it out when its available! :)

[janetournois]

The PR #7754 has been merged, it is now available in the master branch, and will be published in the upcoming CGAL-6.0 version.

[janetournois]

I'm closing this discussion now.
Fell free to reopen it if your meshing issues are still not fixed.