From 8c09f72039360505ffc169fe5d13a64c1f0415de Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Tue, 20 Sep 2022 15:54:54 +0200 Subject: [PATCH 1/6] add debug code --- .../CGAL/Mesh_3/generate_label_weights.h | 60 ++++++++++++++++++- 1 file changed, 58 insertions(+), 2 deletions(-) diff --git a/Mesh_3/include/CGAL/Mesh_3/generate_label_weights.h b/Mesh_3/include/CGAL/Mesh_3/generate_label_weights.h index 6f3d9d77bab..b77f28035e1 100644 --- a/Mesh_3/include/CGAL/Mesh_3/generate_label_weights.h +++ b/Mesh_3/include/CGAL/Mesh_3/generate_label_weights.h @@ -108,6 +108,39 @@ SIGN get_sign() // SGN_UNKNOWN } +template +void convert_itk_to_image_3(itk::Image* const itk_img, + const char* filename = "") +{ + auto t = itk_img->GetOrigin(); + auto v = itk_img->GetSpacing(); + auto region = itk_img->GetRequestedRegion(); + + _image* img + = _createImage(region.GetSize(0), region.GetSize(1), region.GetSize(2), + 1, //vectorial dimension + v[0], v[1], v[2], + sizeof(Image_word_type), //image word size in bytes + internal::get_wordkind(), //image word kind WK_FIXED, WK_FLOAT, WK_UNKNOWN + internal::get_sign()); //image word sign + Image_word_type* img_ptr = (Image_word_type*)(img->data); + + const int size = region.GetSize(0) * region.GetSize(1) * region.GetSize(2); + std::fill(img_ptr, + img_ptr + size, + Image_word_type(0)); + img->tx = t[0]; + img->ty = t[1]; + img->tz = t[2]; + + std::copy(itk_img->GetBufferPointer(), + itk_img->GetBufferPointer() + size, + img_ptr); + + if(filename != "") + _writeImage(img, filename); +} + }//namespace internal /// @cond INTERNAL @@ -140,6 +173,8 @@ CGAL::Image_3 generate_label_weights_with_known_word_type(const CGAL::Image_3& i typename ImageType::Pointer itk_img = ImageType::New(); std::set labels; internal::convert_image_3_to_itk(image, itk_img.GetPointer(), labels); + CGAL_assertion(internal::count_non_white_pixels(image) + == internal::count_non_white_pixels(itk_img.GetPointer())); using DuplicatorType = itk::ImageDuplicator; using IndicatorFilter = itk::BinaryThresholdImageFilter; @@ -178,12 +213,26 @@ CGAL::Image_3 generate_label_weights_with_known_word_type(const CGAL::Image_3& i indicator->SetUpperThreshold(label); indicator->Update(); +#ifdef CGAL_MESH_3_WEIGHTED_IMAGES_DEBUG + std::ostringstream oss; + oss << "indicator_" << id << ".inr.gz"; + std::cout << "filename = " << oss.str().c_str() << std::endl; + internal::convert_itk_to_image_3(indicator->GetOutput(), oss.str().c_str()); +#endif + //perform gaussian smoothing typename GaussianFilterType::Pointer smoother = GaussianFilterType::New(); smoother->SetInput(indicator->GetOutput()); smoother->SetSigma(sigma); smoother->Update(); +#ifdef CGAL_MESH_3_WEIGHTED_IMAGES_DEBUG + std::ostringstream oss1; + oss1 << "smooth_" << id << ".inr.gz"; + std::cout << "filename = " << oss1.str().c_str() << std::endl; + internal::convert_itk_to_image_3(smoother->GetOutput(), oss1.str().c_str()); +#endif + //take the max of smoothed indicator functions if (id == 0) blured_max = smoother->GetOutput(); @@ -199,9 +248,16 @@ CGAL::Image_3 generate_label_weights_with_known_word_type(const CGAL::Image_3& i id++; #ifdef CGAL_MESH_3_WEIGHTED_IMAGES_DEBUG - std::cout << "AFTER MAX (label = " << label << ") : " << std::endl; + std::ostringstream oss2; + oss2 << "max_" << "all" << ".inr.gz"; + std::cout << "filename = " << oss2.str().c_str() << std::endl; + internal::convert_itk_to_image_3(blured_max.GetPointer(), oss2.str().c_str()); +#endif + +#ifdef CGAL_MESH_3_WEIGHTED_IMAGES_DEBUG +// std::cout << "AFTER MAX (label = " << label << ") : " << std::endl; std::cout << "\tnon zero in max (" - << label << ")\t= " << internal::count_non_white_pixels(blured_max.GetPointer()) << std::endl; + << id << ")\t= " << internal::count_non_white_pixels(blured_max.GetPointer()) << std::endl; #endif } From af908bea25876689ed4a127f4fecf8e23cc93161 Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Thu, 22 Sep 2022 12:48:46 +0200 Subject: [PATCH 2/6] use ITK DiscreteGaussianImageFilter with variance = 1 voxel sigma parameter is difficult to use for end user --- Mesh_3/examples/Mesh_3/mesh_3D_weighted_image.cpp | 4 +++- Mesh_3/include/CGAL/Mesh_3/generate_label_weights.h | 13 +++++++++---- 2 files changed, 12 insertions(+), 5 deletions(-) diff --git a/Mesh_3/examples/Mesh_3/mesh_3D_weighted_image.cpp b/Mesh_3/examples/Mesh_3/mesh_3D_weighted_image.cpp index f6b072b3dcd..4fde2cf4831 100644 --- a/Mesh_3/examples/Mesh_3/mesh_3D_weighted_image.cpp +++ b/Mesh_3/examples/Mesh_3/mesh_3D_weighted_image.cpp @@ -1,3 +1,5 @@ +#define CGAL_MESH_3_WEIGHTED_IMAGES_DEBUG + #include #include @@ -39,7 +41,7 @@ int main(int argc, char* argv[]) /// [Loads image] /// [Domain creation] - const float sigma = 10.f; + 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); diff --git a/Mesh_3/include/CGAL/Mesh_3/generate_label_weights.h b/Mesh_3/include/CGAL/Mesh_3/generate_label_weights.h index b77f28035e1..dd1363f0e9b 100644 --- a/Mesh_3/include/CGAL/Mesh_3/generate_label_weights.h +++ b/Mesh_3/include/CGAL/Mesh_3/generate_label_weights.h @@ -21,7 +21,7 @@ #include #include #include -#include +#include #include #include @@ -173,12 +173,15 @@ CGAL::Image_3 generate_label_weights_with_known_word_type(const CGAL::Image_3& i typename ImageType::Pointer itk_img = ImageType::New(); std::set labels; internal::convert_image_3_to_itk(image, itk_img.GetPointer(), labels); + +#ifdef CGAL_MESH_3_WEIGHTED_IMAGES_DEBUG CGAL_assertion(internal::count_non_white_pixels(image) == internal::count_non_white_pixels(itk_img.GetPointer())); +#endif using DuplicatorType = itk::ImageDuplicator; using IndicatorFilter = itk::BinaryThresholdImageFilter; - using GaussianFilterType = itk::RecursiveGaussianImageFilter; + using GaussianFilterType = itk::DiscreteGaussianImageFilter; using MaximumImageFilterType = itk::MaximumImageFilter; std::vector indicators(labels.size()); @@ -222,8 +225,9 @@ CGAL::Image_3 generate_label_weights_with_known_word_type(const CGAL::Image_3& i //perform gaussian smoothing typename GaussianFilterType::Pointer smoother = GaussianFilterType::New(); + smoother->SetUseImageSpacing(true);//variance/std deviation is counted real world distances smoother->SetInput(indicator->GetOutput()); - smoother->SetSigma(sigma); + smoother->SetVariance(sigma*sigma); smoother->Update(); #ifdef CGAL_MESH_3_WEIGHTED_IMAGES_DEBUG @@ -246,6 +250,8 @@ CGAL::Image_3 generate_label_weights_with_known_word_type(const CGAL::Image_3& i } id++; + } + #ifdef CGAL_MESH_3_WEIGHTED_IMAGES_DEBUG std::ostringstream oss2; @@ -259,7 +265,6 @@ CGAL::Image_3 generate_label_weights_with_known_word_type(const CGAL::Image_3& i std::cout << "\tnon zero in max (" << id << ")\t= " << internal::count_non_white_pixels(blured_max.GetPointer()) << std::endl; #endif - } //copy pixels to weights std::copy(blured_max->GetBufferPointer(), From 0f4ac2bd04e82a6d42a6a5cd037ee97d7797ddf7 Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Thu, 22 Sep 2022 14:05:59 +0200 Subject: [PATCH 3/6] remove debug macro from example code --- Mesh_3/examples/Mesh_3/mesh_3D_weighted_image.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/Mesh_3/examples/Mesh_3/mesh_3D_weighted_image.cpp b/Mesh_3/examples/Mesh_3/mesh_3D_weighted_image.cpp index 4fde2cf4831..0a8f47a8dfa 100644 --- a/Mesh_3/examples/Mesh_3/mesh_3D_weighted_image.cpp +++ b/Mesh_3/examples/Mesh_3/mesh_3D_weighted_image.cpp @@ -1,5 +1,3 @@ -#define CGAL_MESH_3_WEIGHTED_IMAGES_DEBUG - #include #include From d6d2188cc7950fc5c9d08c60894a73053f373fd2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Wed, 5 Oct 2022 19:33:00 +0200 Subject: [PATCH 4/6] fix warning --- Mesh_3/include/CGAL/Mesh_3/generate_label_weights.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Mesh_3/include/CGAL/Mesh_3/generate_label_weights.h b/Mesh_3/include/CGAL/Mesh_3/generate_label_weights.h index dd1363f0e9b..cc8aeb61173 100644 --- a/Mesh_3/include/CGAL/Mesh_3/generate_label_weights.h +++ b/Mesh_3/include/CGAL/Mesh_3/generate_label_weights.h @@ -137,7 +137,7 @@ void convert_itk_to_image_3(itk::Image* const itk_img, itk_img->GetBufferPointer() + size, img_ptr); - if(filename != "") + if(std::strlen(filename)!=0) _writeImage(img, filename); } From c64a0d86484f66dc8844fb341c3c35ba9685c82c Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Thu, 6 Oct 2022 10:44:26 +0200 Subject: [PATCH 5/6] add more details about how to choose sigma --- Mesh_3/include/CGAL/Mesh_3/generate_label_weights.h | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/Mesh_3/include/CGAL/Mesh_3/generate_label_weights.h b/Mesh_3/include/CGAL/Mesh_3/generate_label_weights.h index cc8aeb61173..d34f849f35b 100644 --- a/Mesh_3/include/CGAL/Mesh_3/generate_label_weights.h +++ b/Mesh_3/include/CGAL/Mesh_3/generate_label_weights.h @@ -296,7 +296,13 @@ CGAL::Image_3 generate_label_weights_with_known_word_type(const CGAL::Image_3& i * * @param image the input labeled image from which the weights image is computed. * Both will then be used to construct a `Labeled_mesh_domain_3`. -* @param sigma the standard deviation parameter of the internal Gaussian filter +* @param sigma the standard deviation parameter of the internal Gaussian filter, +* measured in real-world distances. The size of a voxel (e.g. shortest length +* or longest length) usually is a good value for this parameter. +* Note that if `sigma` is too small, the "stair-effect" of meshing from +* a voxel image can appear. On the other side, if `sigma` is too large, +* thin volumes (basically one voxel thick) may be lost in the meshing process +* because the computed weights are too blurry. * * @returns a `CGAL::Image_3` of weights used to build a quality `Labeled_mesh_domain_3`, * with the same dimensions as `image` From 9bb4fb4d1516223af458a4ea7ce8d236bfb61df6 Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Thu, 6 Oct 2022 15:10:37 +0200 Subject: [PATCH 6/6] protect debug code with macro and remove default filename "" --- Mesh_3/include/CGAL/Mesh_3/generate_label_weights.h | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/Mesh_3/include/CGAL/Mesh_3/generate_label_weights.h b/Mesh_3/include/CGAL/Mesh_3/generate_label_weights.h index d34f849f35b..7343015f1da 100644 --- a/Mesh_3/include/CGAL/Mesh_3/generate_label_weights.h +++ b/Mesh_3/include/CGAL/Mesh_3/generate_label_weights.h @@ -108,9 +108,10 @@ SIGN get_sign() // SGN_UNKNOWN } +#ifdef CGAL_MESH_3_WEIGHTED_IMAGES_DEBUG template void convert_itk_to_image_3(itk::Image* const itk_img, - const char* filename = "") + const char* filename) { auto t = itk_img->GetOrigin(); auto v = itk_img->GetSpacing(); @@ -137,9 +138,9 @@ void convert_itk_to_image_3(itk::Image* const itk_img, itk_img->GetBufferPointer() + size, img_ptr); - if(std::strlen(filename)!=0) - _writeImage(img, filename); + _writeImage(img, filename); } +#endif }//namespace internal