Mesh_3::Add handling of *.nii files in one example

Mesh_3/doc/Mesh_3/Mesh_3.txt

@@ -693,7 +693,7 @@ View of a remeshed surface. (Left) input mesh (Right) output mesh. Code from sub
 
 The following example produces a 3D mesh for a domain whose boundary surface
 is the isosurface associated to an isovalue inside the input gray-level
-3D image. In the distribution you can also find the example \ref Mesh_3/mesh_3D_gray_vtk_image.cpp which can deal with DICOM files as input.
+3D image. In the distribution you can also find the example \ref Mesh_3/mesh_3D_gray_vtk_image.cpp which can deal with `*.nii` as well as `DICOM` files as input.
 
 \cgalExample{Mesh_3/mesh_3D_gray_image.cpp}
 

Mesh_3/examples/Mesh_3/CMakeLists.txt

@@ -121,7 +121,7 @@ target_link_libraries(mesh_polyhedral_complex_sm PUBLIC CGAL::Eigen3_support)
 if(TARGET CGAL::CGAL_ImageIO)
   if(VTK_FOUND AND ("${VTK_VERSION_MAJOR}" GREATER "5" OR VTK_VERSION
                                                           VERSION_GREATER 5))
-    add_executable(mesh_3D_gray_vtk_image mesh_3D_gray_vtk_image.cpp)
+    create_single_source_cgal_program("mesh_3D_gray_vtk_image.cpp")
     target_link_libraries(
       mesh_3D_gray_vtk_image
       PUBLIC CGAL::Eigen3_support CGAL::CGAL CGAL::CGAL_ImageIO

Mesh_3/examples/Mesh_3/mesh_3D_gray_vtk_image.cpp

@@ -1,6 +1,8 @@
 
+#include <vtkNew.h>
 #include <vtkImageData.h>
 #include <vtkDICOMImageReader.h>
+#include <vtkNIFTIImageReader.h>
 #include <vtkImageReader.h>
 #include <vtkImageGaussianSmooth.h>
 #include <vtkDemandDrivenPipeline.h>
@@ -19,6 +21,8 @@
 #include <boost/lexical_cast.hpp>
 #include <boost/functional.hpp>
 
+#include <filesystem>
+
 typedef short Image_word_type;
 
 // Domain
@@ -43,36 +47,59 @@ class Less {
   }
 };
 
+namespace fs = std::filesystem;
+
 int main(int argc, char* argv[])
 {
   // Loads image
-  if(argc == 1){
-    std::cerr << "Usage:  " << argv[0] << " <directory with dicom data> iso_level=1  facet_size=1  facet_distance=0.1  cell_size=1\n";
-    return 0;
-  }
 
+  // Usage: mesh_3D_gray_vtk_image <nii file or dicom directory> iso_level=1  facet_size=1  facet_distance=0.1  cell_size=1
+
+  const std::string fname = (argc>1)?argv[1]:CGAL::data_file_path("images/squircle.nii");
+
+  vtkImageData* vtk_image = nullptr;
   Image_word_type iso = (argc>2)? boost::lexical_cast<Image_word_type>(argv[2]): 1;
   double fs = (argc>3)? boost::lexical_cast<double>(argv[3]): 1;
   double fd = (argc>4)? boost::lexical_cast<double>(argv[4]): 0.1;
   double cs = (argc>5)? boost::lexical_cast<double>(argv[5]): 1;
 
-  vtkDICOMImageReader*dicom_reader = vtkDICOMImageReader::New();
-  dicom_reader->SetDirectoryName(argv[1]);
-
-  vtkDemandDrivenPipeline*executive =
-    vtkDemandDrivenPipeline::SafeDownCast(dicom_reader->GetExecutive());
-  if (executive)
-    {
-      executive->SetReleaseDataFlag(0, 0); // where 0 is the port index
+  fs::path path(fname);
+
+  if(fs::is_regular_file(path)){
+    std::cout << "regular file" << std::endl;
+    if (path.has_extension()){
+      fs::path stem = path.stem();
+      if ((path.extension() == ".nii") || (stem.has_extension() && (stem.extension() == ".nii") && (path.extension() == ".gz"))) {
+        vtkNIFTIImageReader* reader = vtkNIFTIImageReader::New();
+        reader->SetFileName(fname.c_str());
+        reader->Update();
+        vtk_image = reader->GetOutput();
+        vtk_image->Print(std::cerr);
+      }
     }
-
-  vtkImageGaussianSmooth* smoother = vtkImageGaussianSmooth::New();
-  smoother->SetStandardDeviations(1., 1., 1.);
-  smoother->SetInputConnection(dicom_reader->GetOutputPort());
-  smoother->Update();
-  vtkImageData* vtk_image = smoother->GetOutput();
-  vtk_image->Print(std::cerr);
-
+  }
+  else if (fs::is_directory(path)) {
+    vtkDICOMImageReader* dicom_reader = vtkDICOMImageReader::New();
+    dicom_reader->SetDirectoryName(argv[1]);
+
+    vtkDemandDrivenPipeline* executive =
+      vtkDemandDrivenPipeline::SafeDownCast(dicom_reader->GetExecutive());
+    if (executive)
+      {
+        executive->SetReleaseDataFlag(0, 0); // where 0 is the port index
+      }
+
+    vtkImageGaussianSmooth* smoother = vtkImageGaussianSmooth::New();
+    smoother->SetStandardDeviations(1., 1., 1.);
+    smoother->SetInputConnection(dicom_reader->GetOutputPort());
+    smoother->Update();
+    vtk_image = smoother->GetOutput();
+    vtk_image->Print(std::cerr);
+  }
+  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";
@@ -84,7 +111,7 @@ int main(int argc, char* argv[])
   Mesh_domain domain = Mesh_domain::create_gray_image_mesh_domain
     (image,
      params::image_values_to_subdomain_indices(Less(iso)).
-             value_outside(0));
+             value_outside(iso+1));
   /// [Domain creation]
 
   // Mesh criteria