songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #6864 from janetournois/Mesh_3-fix_weighted_images-GF 

Mesh 3 - fix construction of weights image for labeled image input
This commit is contained in:
Laurent Rineau 2022-10-13 15:14:11 +02:00
parent ae8101d615 9bb4fb4d15
commit adc5bd4677
2 changed files with 77 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Mesh_3/examples/Mesh_3/mesh_3D_weighted_image.cpp
View File

@ -39,7 +39,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);

84
Mesh_3/include/CGAL/Mesh_3/generate_label_weights.h
View File

@ -21,7 +21,7 @@
#include <itkImage.h>
#include <itkImageDuplicator.h>
#include <itkBinaryThresholdImageFilter.h>
#include <itkRecursiveGaussianImageFilter.h>
#include <itkDiscreteGaussianImageFilter.h>
#include <itkMaximumImageFilter.h>
#include <iostream>
@ -108,6 +108,40 @@ SIGN get_sign()
// SGN_UNKNOWN
}
#ifdef CGAL_MESH_3_WEIGHTED_IMAGES_DEBUG
template<typename Image_word_type>
void convert_itk_to_image_3(itk::Image<Image_word_type, 3>* 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_type>(), //image word kind WK_FIXED, WK_FLOAT, WK_UNKNOWN
internal::get_sign<Image_word_type>()); //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);
_writeImage(img, filename);
}
#endif
}//namespace internal
/// @cond INTERNAL
@ -141,9 +175,14 @@ CGAL::Image_3 generate_label_weights_with_known_word_type(const CGAL::Image_3& i
std::set<Image_word_type> 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_word_type>(image)
== internal::count_non_white_pixels<Image_word_type>(itk_img.GetPointer()));
#endif
using DuplicatorType = itk::ImageDuplicator<ImageType>;
using IndicatorFilter = itk::BinaryThresholdImageFilter<ImageType, WeightsType>;
using GaussianFilterType = itk::RecursiveGaussianImageFilter<WeightsType, WeightsType>;
using GaussianFilterType = itk::DiscreteGaussianImageFilter<WeightsType, WeightsType>;
using MaximumImageFilterType = itk::MaximumImageFilter<WeightsType>;
std::vector<typename ImageType::Pointer> indicators(labels.size());
@ -178,12 +217,27 @@ 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->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
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();
@ -197,13 +251,21 @@ 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::cout << "\tnon zero in max ("
<< label << ")\t= " << internal::count_non_white_pixels(blured_max.GetPointer()) << 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 ("
<< id << ")\t= " << internal::count_non_white_pixels(blured_max.GetPointer()) << std::endl;
#endif
}
//copy pixels to weights
std::copy(blured_max->GetBufferPointer(),
@ -235,7 +297,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`