diff --git a/BGL/include/CGAL/boost/parameter.h b/BGL/include/CGAL/boost/parameter.h index 9c28ac96b7a2..14ab09b04748 100644 --- a/BGL/include/CGAL/boost/parameter.h +++ b/BGL/include/CGAL/boost/parameter.h @@ -110,6 +110,7 @@ BOOST_PARAMETER_NAME( (pointer_to_stop_atomic_boolean, tag ) pointer_to_stop_ato BOOST_PARAMETER_NAME( (function, tag ) function_) BOOST_PARAMETER_NAME( (bounding_object, tag ) bounding_object_) BOOST_PARAMETER_NAME( (relative_error_bound, tag ) relative_error_bound_) +BOOST_PARAMETER_NAME( (weights, tag) weights_) BOOST_PARAMETER_NAME( (p_rng, tag ) p_rng_) BOOST_PARAMETER_NAME( (null_subdomain_index, tag ) null_subdomain_index_) BOOST_PARAMETER_NAME( (construct_surface_patch_index, tag ) construct_surface_patch_index_) diff --git a/CGAL_ImageIO/include/CGAL/Image_3.h b/CGAL_ImageIO/include/CGAL/Image_3.h index 8ede8eefec57..e0ecd7ce77b9 100644 --- a/CGAL_ImageIO/include/CGAL/Image_3.h +++ b/CGAL_ImageIO/include/CGAL/Image_3.h @@ -102,7 +102,7 @@ class CGAL_IMAGEIO_EXPORT Image_3 // std::cerr << "Image_3::copy_constructor\n"; } - Image_3(_image* im, Own own_the_data = OWN_THE_DATA) + explicit Image_3(_image* im, Own own_the_data = OWN_THE_DATA) { private_read(im, own_the_data); } @@ -146,9 +146,9 @@ class CGAL_IMAGEIO_EXPORT Image_3 double vy() const { return image_ptr->vy; } double vz() const { return image_ptr->vz; } - double tx() const { return image_ptr->tx; } - double ty() const { return image_ptr->ty; } - double tz() const { return image_ptr->tz; } + float tx() const { return image_ptr->tx; } + float ty() const { return image_ptr->ty; } + float tz() const { return image_ptr->tz; } float value(const std::size_t i, const std::size_t j, @@ -157,6 +157,11 @@ class CGAL_IMAGEIO_EXPORT Image_3 return ::evaluate(image(),i,j,k); } + bool is_valid() const + { + return image_ptr.get() != nullptr; + } + public: bool read(const char* file) diff --git a/Documentation/doc/biblio/cgal_manual.bib b/Documentation/doc/biblio/cgal_manual.bib index 9ac6cdc2a186..1c49fe6ffd27 100644 --- a/Documentation/doc/biblio/cgal_manual.bib +++ b/Documentation/doc/biblio/cgal_manual.bib @@ -3256,6 +3256,12 @@ @article{cgal:Wwshap-eepec-20 publisher = {ACM} } +@article{stalling1998weighted, + title={Weighted labels for 3D image segmentation}, + author={Stalling, Detlev and Z{\"o}ckler, Malte and Sander, Oliver and Hege, Hans-Christian}, + year={1998} +} + % ---------------------------------------------------------------------------- % END OF BIBFILE % ---------------------------------------------------------------------------- diff --git a/Installation/cmake/modules/CGAL_ITK_support.cmake b/Installation/cmake/modules/CGAL_ITK_support.cmake new file mode 100644 index 000000000000..61ed67e2f2ba --- /dev/null +++ b/Installation/cmake/modules/CGAL_ITK_support.cmake @@ -0,0 +1,7 @@ +if(ITK_FOUND AND NOT TARGET CGAL::ITK_support) + add_library(CGAL::ITK_support INTERFACE IMPORTED) + set_target_properties(CGAL::ITK_support PROPERTIES + INTERFACE_COMPILE_DEFINITIONS "CGAL_USE_ITK" + INTERFACE_INCLUDE_DIRECTORIES "${ITK_INCLUDE_DIRS}" + INTERFACE_LINK_LIBRARIES "${ITK_LIBRARIES}") +endif() diff --git a/Installation/cmake/modules/list_of_whitelisted_headers.cmake b/Installation/cmake/modules/list_of_whitelisted_headers.cmake index 5aa71d9560b5..9cbef5d8a577 100644 --- a/Installation/cmake/modules/list_of_whitelisted_headers.cmake +++ b/Installation/cmake/modules/list_of_whitelisted_headers.cmake @@ -41,7 +41,7 @@ set(list_of_whitelisted_headers_txt [=[ CGAL/IO/Triangulation_geomview_ostream_3.h CGAL/IO/Triangulation_geomview_ostream_2.h CGAL/IO/Polyhedron_geomview_ostream.h - + CGAL/Mesh_3/generate_label_weights.h ]=]) separate_arguments(list_of_whitelisted_headers UNIX_COMMAND ${list_of_whitelisted_headers_txt}) diff --git a/Mesh_3/doc/Mesh_3/CGAL/Labeled_mesh_domain_3.h b/Mesh_3/doc/Mesh_3/CGAL/Labeled_mesh_domain_3.h index 1dbef3e1efb8..59808d448996 100644 --- a/Mesh_3/doc/Mesh_3/CGAL/Labeled_mesh_domain_3.h +++ b/Mesh_3/doc/Mesh_3/CGAL/Labeled_mesh_domain_3.h @@ -233,6 +233,12 @@ The parameters are optional unless otherwise specified.
  • `parameters::image` (mandatory) the input 3D image. Must be a `CGAL::Image_3` object. +
  • `parameters::weights` an input 3D image that provides +weights associated to each voxel (the word type is `unsigned char`, +and the voxels values are integers between 0 and 255). +The weights image can be generated with `CGAL::Mesh_3::generate_label_weights()`. +Its dimensions must be the same as the dimensions of `parameters::image`. +
  • `parameter::value_outside` the value attached to voxels outside of the domain to be meshed. Its default value is `0`. @@ -246,6 +252,11 @@ From the example (\ref Mesh_3/mesh_3D_image.cpp): \snippet Mesh_3/mesh_3D_image.cpp Domain creation +From the example (\ref Mesh_3/mesh_3D_weighted_image.cpp), +where the labeled image is used with a precomputed 3D image of weights : + +\snippet Mesh_3/mesh_3D_weighted_image.cpp Domain creation + */ template static diff --git a/Mesh_3/doc/Mesh_3/Doxyfile.in b/Mesh_3/doc/Mesh_3/Doxyfile.in index 5c8e03ccafd0..62ea73403e60 100644 --- a/Mesh_3/doc/Mesh_3/Doxyfile.in +++ b/Mesh_3/doc/Mesh_3/Doxyfile.in @@ -5,7 +5,8 @@ ALIASES += "cgalDescribePolylineType=A polyline is defined as a sequence of poin INPUT += \ ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Polyhedral_complex_mesh_domain_3.h \ - ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Mesh_domain_with_polyline_features_3.h + ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Mesh_domain_with_polyline_features_3.h \ + ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Mesh_3/generate_label_weights.h PROJECT_NAME = "CGAL ${CGAL_DOC_VERSION} - 3D Mesh Generation" HTML_EXTRA_FILES = ${CGAL_PACKAGE_DOC_DIR}/fig/implicit_domain_3.jpg \ diff --git a/Mesh_3/doc/Mesh_3/Mesh_3.txt b/Mesh_3/doc/Mesh_3/Mesh_3.txt index a8b806d783c7..08e82efd875c 100644 --- a/Mesh_3/doc/Mesh_3/Mesh_3.txt +++ b/Mesh_3/doc/Mesh_3/Mesh_3.txt @@ -706,6 +706,24 @@ The resulting mesh is shown in \cgalFigureRef{figureliver_3d_image_mesh}. Cut view of a 3D mesh produced from a segmented liver image. Code from subsection \ref Mesh_3_subsection_examples_3d_image generates this file. \cgalFigureEnd +\subsubsection Mesh_3DomainsFrom3DImagesWithWeights Domains From Segmented 3D Images, with Weights + +When a segmented image is given as input, the generated mesh surface sometimes sticks too closely +to the voxels surface, causing an aliasing effect. +A solution to generate a smooth and accurate output surface was described by Stalling et al in +\cgalCite{stalling1998weighted}. It consists in generating a second input image, made +of integer coefficients called *weights*, and use those weights to define smoother domain boundaries. +The 3D image of weights can be generated using `CGAL::Mesh_3::generate_weights()` as shown in +the following example. + +\cgalExample{Mesh_3/mesh_3D_weighted_image.cpp} + +\cgalFigureBegin{figure_weightedImage, weighted_images.jpg} +Surface of the output mesh generated with a very small `facet_distance` +without the weights (left, 25563 vertices) and with the weights (right, 19936 vertices). +\cgalFigureEnd + + \subsubsection Mesh_3DomainsFrom3DImagesWithCustomInitialization Domains From 3D Images, with a Custom Initialization The example \ref Mesh_3/mesh_3D_image_with_custom_initialization.cpp is a modification diff --git a/Mesh_3/doc/Mesh_3/PackageDescription.txt b/Mesh_3/doc/Mesh_3/PackageDescription.txt index f7a09d8f4738..240ca4bf2941 100644 --- a/Mesh_3/doc/Mesh_3/PackageDescription.txt +++ b/Mesh_3/doc/Mesh_3/PackageDescription.txt @@ -115,6 +115,7 @@ and their associated classes: - `CGAL::lloyd_optimize_mesh_3()` - `CGAL::odt_optimize_mesh_3()` - `CGAL::facets_in_complex_3_to_triangle_mesh()` +- `CGAL::Mesh_3::generate_label_weights()` \cgalCRPSection{CGAL::parameters Functions} diff --git a/Mesh_3/doc/Mesh_3/examples.txt b/Mesh_3/doc/Mesh_3/examples.txt index c1c5e4de1a81..aa44908504c1 100644 --- a/Mesh_3/doc/Mesh_3/examples.txt +++ b/Mesh_3/doc/Mesh_3/examples.txt @@ -4,6 +4,7 @@ \example Mesh_3/mesh_3D_gray_image_with_custom_initialization.cpp \example Mesh_3/mesh_3D_image_with_features.cpp \example Mesh_3/mesh_3D_image_with_custom_initialization.cpp +\example Mesh_3/mesh_3D_weighted_image.cpp \example Mesh_3/random_labeled_image.h \example CGAL/Mesh_3/initialize_triangulation_from_gray_image.h \example CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h diff --git a/Mesh_3/doc/Mesh_3/fig/weighted_images.jpg b/Mesh_3/doc/Mesh_3/fig/weighted_images.jpg new file mode 100644 index 000000000000..39a3dca2f9cf Binary files /dev/null and b/Mesh_3/doc/Mesh_3/fig/weighted_images.jpg differ diff --git a/Mesh_3/examples/Mesh_3/CMakeLists.txt b/Mesh_3/examples/Mesh_3/CMakeLists.txt index 3c1602a2e267..c1230c541dee 100644 --- a/Mesh_3/examples/Mesh_3/CMakeLists.txt +++ b/Mesh_3/examples/Mesh_3/CMakeLists.txt @@ -161,10 +161,21 @@ if(TARGET CGAL::CGAL_ImageIO) create_single_source_cgal_program("mesh_3D_image_variable_size.cpp") target_link_libraries(mesh_3D_image_variable_size PUBLIC CGAL::Eigen3_support) + + find_package(ITK QUIET COMPONENTS ITKCommon ITKThresholding ITKSmoothing ITKImageIntensity) + if(ITK_FOUND) + include(CGAL_ITK_support) + create_single_source_cgal_program("mesh_3D_weighted_image.cpp") + target_link_libraries(mesh_3D_weighted_image + PUBLIC CGAL::Eigen3_support CGAL::ITK_support) + else(ITK_FOUND) + message(STATUS "NOTICE: The examples that need ITK will not be compiled.") + endif(ITK_FOUND) + else() message( STATUS - "NOTICE: The examples mesh_3D_image.cpp, mesh_3D_image_variable_size.cpp, mesh_optimization_example.cpp and mesh_optimization_lloyd_example.cpp need CGAL_ImageIO to be configured with ZLIB support, and will not be compiled." + "NOTICE: The examples mesh_3D_image.cpp, mesh_3D_weighted_image.cpp, mesh_3D_image_variable_size.cpp, mesh_optimization_example.cpp and mesh_optimization_lloyd_example.cpp need CGAL_ImageIO to be configured with ZLIB support, and will not be compiled." ) endif() @@ -179,6 +190,7 @@ if(CGAL_ACTIVATE_CONCURRENT_MESH_3 AND TARGET CGAL::TBB_support) foreach( target mesh_3D_image + mesh_3D_weighted_image mesh_3D_image_variable_size mesh_3D_image_with_custom_initialization mesh_3D_gray_image_with_custom_initialization diff --git a/Mesh_3/examples/Mesh_3/mesh_3D_weighted_image.cpp b/Mesh_3/examples/Mesh_3/mesh_3D_weighted_image.cpp new file mode 100644 index 000000000000..0dc3e8e8be2e --- /dev/null +++ b/Mesh_3/examples/Mesh_3/mesh_3D_weighted_image.cpp @@ -0,0 +1,67 @@ +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +// Domain +using K = CGAL::Exact_predicates_inexact_constructions_kernel; +using Mesh_domain = CGAL::Labeled_mesh_domain_3; + +// Triangulation +using Tr = CGAL::Mesh_triangulation_3::type; +using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3; + +// Criteria +using Mesh_criteria = CGAL::Mesh_criteria_3; + +// To avoid verbose function and named parameters call +using namespace CGAL::parameters; + +int main(int argc, char* argv[]) +{ + /// [Loads image] + const char* fname = (argc > 1) ? argv[1] : "data/liver.inr.gz"; + CGAL::Image_3 image; + if(!image.read(fname)){ + std::cerr << "Error: Cannot read file " << fname << std::endl; + return EXIT_FAILURE; + } + /// [Loads image] + + /// [Domain creation] + const float sigma = 10.f; + CGAL::Image_3 img_weights = + CGAL::Mesh_3::generate_label_weights(image, sigma); + + Mesh_domain domain + = Mesh_domain::create_labeled_image_mesh_domain(image, + weights = img_weights, + relative_error_bound = 1e-6); + /// [Domain creation] + + // Mesh criteria + Mesh_criteria criteria(facet_angle=30, facet_size=6, facet_distance=0.5, + cell_radius_edge_ratio=3, cell_size=8); + + /// [Meshing] + C3t3 c3t3 = CGAL::make_mesh_3(domain, criteria); + /// [Meshing] + + // Output + std::ofstream medit_file("out.mesh"); + c3t3.output_to_medit(medit_file); + std::ofstream bin_file("out.binary.cgal", std::ios_base::binary); + CGAL::IO::save_binary_file(bin_file, c3t3); + + return 0; +} diff --git a/Mesh_3/include/CGAL/Labeled_mesh_domain_3.h b/Mesh_3/include/CGAL/Labeled_mesh_domain_3.h index e2b1af1e8994..1e18a4dc124c 100644 --- a/Mesh_3/include/CGAL/Labeled_mesh_domain_3.h +++ b/Mesh_3/include/CGAL/Labeled_mesh_domain_3.h @@ -39,6 +39,7 @@ // support for `CGAL::Image_3` #include #include +#include // support for implicit functions #include @@ -418,6 +419,7 @@ class Labeled_mesh_domain_3 : (optional (relative_error_bound_, (const FT&), FT(1e-3)) + (weights_, (const CGAL::Image_3&), CGAL::Image_3()) (value_outside_, *, 0) (p_rng_, (CGAL::Random*), (CGAL::Random*)(0)) (image_values_to_subdomain_indices_, *, @@ -429,18 +431,37 @@ class Labeled_mesh_domain_3 : ) { namespace p = CGAL::parameters; - return Labeled_mesh_domain_3 - (create_labeled_image_wrapper + if (weights_.is_valid()) + { + return Labeled_mesh_domain_3 + (create_weighted_labeled_image_wrapper (image_, + weights_, image_values_to_subdomain_indices_, value_outside_), - Mesh_3::internal::compute_bounding_box(image_), - p::relative_error_bound = relative_error_bound_, - p::p_rng = p_rng_, - p::null_subdomain_index = - create_null_subdomain_index(null_subdomain_index_), - p::construct_surface_patch_index = - create_construct_surface_patch_index(construct_surface_patch_index_)); + Mesh_3::internal::compute_bounding_box(image_), + p::relative_error_bound = relative_error_bound_, + p::p_rng = p_rng_, + p::null_subdomain_index = + create_null_subdomain_index(null_subdomain_index_), + p::construct_surface_patch_index = + create_construct_surface_patch_index(construct_surface_patch_index_)); + } + else + { + return Labeled_mesh_domain_3 + (create_labeled_image_wrapper + (image_, + image_values_to_subdomain_indices_, + value_outside_), + Mesh_3::internal::compute_bounding_box(image_), + p::relative_error_bound = relative_error_bound_, + p::p_rng = p_rng_, + p::null_subdomain_index = + create_null_subdomain_index(null_subdomain_index_), + p::construct_surface_patch_index = + create_construct_surface_patch_index(construct_surface_patch_index_)); + } } BOOST_PARAMETER_MEMBER_FUNCTION( @@ -866,6 +887,33 @@ class Labeled_mesh_domain_3 : transform_fct(value_outside)); } + template + static + Function + create_weighted_labeled_image_wrapper_with_know_word_type + (const CGAL::Image_3& image, + const CGAL::Image_3& weights, + const Functor& image_values_to_subdomain_indices, + const FT& value_outside) + { + using Mesh_3::internal::Create_labeled_image_values_to_subdomain_indices; + typedef Create_labeled_image_values_to_subdomain_indices C_i_v_t_s_i; + typedef typename C_i_v_t_s_i::type Image_values_to_subdomain_indices; + Image_values_to_subdomain_indices transform_fct = + C_i_v_t_s_i()(image_values_to_subdomain_indices); + + typedef Mesh_3::Image_plus_weights_to_labeled_function_wrapper< + Image_word_type, + int, //interpolation type + unsigned char, // Weights_type, + Subdomain_index> Wrapper; + return Wrapper(image, + weights, + transform_fct, + transform_fct(value_outside)); + } + template static Function @@ -885,6 +933,27 @@ class Labeled_mesh_domain_3 : return Function(); } + template + static + Function + create_weighted_labeled_image_wrapper(const CGAL::Image_3& image, + const CGAL::Image_3& weights, + const Functor& image_values_to_subdomain_indices, + const FT& value_outside) + { + CGAL_IMAGE_IO_CASE(image.image(), + return create_weighted_labeled_image_wrapper_with_know_word_type + (image, + weights, + image_values_to_subdomain_indices, + value_outside); + ); + CGAL_error_msg("This place should never be reached, because it would mean " + "the image word type is a type that is not handled by " + "CGAL_ImageIO."); + return Function(); + } + static Construct_surface_patch_index create_construct_surface_patch_index(const Null_functor&) { diff --git a/Mesh_3/include/CGAL/Mesh_3/Image_plus_weights_to_labeled_function_wrapper.h b/Mesh_3/include/CGAL/Mesh_3/Image_plus_weights_to_labeled_function_wrapper.h new file mode 100644 index 000000000000..e1d54932eb69 --- /dev/null +++ b/Mesh_3/include/CGAL/Mesh_3/Image_plus_weights_to_labeled_function_wrapper.h @@ -0,0 +1,141 @@ +// Copyright (c) 2016,2021 GeometryFactory Sarl (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// +// Author(s) : Laurent Rineau +// + +#ifndef CGAL_MESH_3_IMAGE_PLUS_WEIGHTS_TO_LABELED_FUNCTION_WRAPPER_H +#define CGAL_MESH_3_IMAGE_PLUS_WEIGHTS_TO_LABELED_FUNCTION_WRAPPER_H + +#include + +#include +#include + +namespace CGAL { + +namespace ImageIO { + +template +class Weighted_indicator_factory +{ + const Image_3& image; + const Image_3& weights; + +public: + + Weighted_indicator_factory(const Image_3& image, + const Image_3& weights) + : image(image) + , weights(weights) + {} + + class Indicator + { + const Weighted_indicator_factory& f; + const Image_word_type label; + public: + Indicator(const Weighted_indicator_factory& f, + const Image_word_type label) : f(f), label(label) {} + + double operator()(const Image_word_type& x) const + { + const std::ptrdiff_t offset = &x - (Image_word_type *)f.image.data(); + const Weights_type w = (std::max)( + Weights_type(128), ((Weights_type *)f.weights.data())[offset]); + return (x == label) ? w : (255 - w); + } + }; // end nested class Indicator + + Indicator indicator(const Image_word_type label) const + { + return Indicator(*this, label); + } +}; // end template class ImageIO::Weighted_indicator_factory + +} // end namespace CGAL::ImageIO + +namespace Mesh_3 { + + +/** + * @class Image_plus_weights_to_labeled_function_wrapper + * + * Wraps a pair of images into a labeled function which takes his values into + * N. Uses weighted trilinear interpolation. + */ +template +class Image_plus_weights_to_labeled_function_wrapper +{ +public: + typedef std::function + Image_values_to_labels; + + // Types + typedef Return_type return_type; + typedef Image_word_type word_type; + typedef CGAL::Image_3 Image_; + + /// Constructor + Image_plus_weights_to_labeled_function_wrapper + (const Image_& image, + const Image_& weights_image, + Image_values_to_labels transform = Identity(), + const Interpolation_type value_outside = Interpolation_type()) + : r_im_(image) + , r_weights_im_(weights_image) + , transform(transform) + , value_outside(value_outside) + , indicator_factory(image, weights_image) + { + } + + // Default copy constructor and assignment operator are ok + + /// Destructor + ~Image_plus_weights_to_labeled_function_wrapper() {} + + /** + * Returns an int corresponding to the label at point \c p + * @param p the input point + * @return the label at point \c p + */ + template + return_type operator()(const Point_3& p) const + { + return static_cast(transform( + r_im_.template labellized_trilinear_interpolation( + CGAL::to_double(p.x()), + CGAL::to_double(p.y()), + CGAL::to_double(p.z()), + value_outside, + indicator_factory))); + } + +private: + /// Labeled image to wrap + const Image_& r_im_; + const Image_& r_weights_im_; + const Image_values_to_labels transform; + const Interpolation_type value_outside; + CGAL::ImageIO::Weighted_indicator_factory indicator_factory; +}; // end class Image_plus_weights_to_labeled_function_wrapper + + +} // end namespace Mesh_3 + +} // end namespace CGAL + +#endif // CGAL_MESH_3_IMAGE_PLUS_WEIGHTS_TO_LABELED_FUNCTION_WRAPPER_H diff --git a/Mesh_3/include/CGAL/Mesh_3/generate_label_weights.h b/Mesh_3/include/CGAL/Mesh_3/generate_label_weights.h new file mode 100644 index 000000000000..6f3d9d77babf --- /dev/null +++ b/Mesh_3/include/CGAL/Mesh_3/generate_label_weights.h @@ -0,0 +1,259 @@ +// Copyright (c) 2021 GeometryFactory +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// +// Author(s) : Laurent Rineau, Jane Tournois + +#ifndef CGAL_MESH_3_GENERATE_LABEL_WEIGHTS_H +#define CGAL_MESH_3_GENERATE_LABEL_WEIGHTS_H + +#include + +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace CGAL { +namespace Mesh_3 { +namespace internal { + +template +void convert_image_3_to_itk(const CGAL::Image_3& image, + itk::Image* const itk_img, + LabelsSet& labels) +{ + const double spacing[3] = {image.vx(), image.vy(), image.vz()}; + itk_img->SetSpacing(spacing); + + const double origin[3] = {image.tx(), image.ty(), image.tz()}; + itk_img->SetOrigin(origin); + + using ImageType = itk::Image; + typename ImageType::IndexType corner = {{0, 0, 0 }}; + typename ImageType::SizeType size = {{image.xdim(), image.ydim(), image.zdim()}}; + typename ImageType::RegionType region(corner, size); + itk_img->SetRegions(region); + + itk_img->Allocate(); + + const Image_word_type* img_begin = static_cast(image.data()); + std::copy(img_begin, img_begin + image.size(), itk_img->GetBufferPointer()); + + labels.insert(img_begin, img_begin + image.size()); +} + +#ifdef CGAL_MESH_3_WEIGHTED_IMAGES_DEBUG +template +int count_non_white_pixels(const CGAL::Image_3& image) +{ + auto diff255 = [&](const Image_word_type p) + { + return p != 255; + }; + const Image_word_type* img_begin = static_cast(image.data()); + return std::count_if(img_begin, + img_begin + image.size(), + diff255); +} + +template +int count_non_white_pixels(itk::Image* itk_img) +{ + auto diff255 = [&](const Image_word_type p) + { + return p != 255; + }; + auto size = itk_img->GetLargestPossibleRegion().GetSize(); + return std::count_if(itk_img->GetBufferPointer(), + itk_img->GetBufferPointer() + size[0]*size[1]*size[2], + diff255); +} +#endif //CGAL_MESH_3_WEIGHTED_IMAGES_DEBUG + +template +WORD_KIND get_wordkind() +{ + if (std::is_floating_point::value) + return WK_FLOAT; + else + return WK_FIXED; +/** unknown (uninitialized) */ +// WK_UNKNOWN +} + +template +SIGN get_sign() +{ + if (std::is_signed::value) + return SGN_SIGNED; + else + return SGN_UNSIGNED; +/** unknown (uninitialized or floating point words) */ +// SGN_UNKNOWN +} + +}//namespace internal + +/// @cond INTERNAL +template +CGAL::Image_3 generate_label_weights_with_known_word_type(const CGAL::Image_3& image, + const float& sigma) +{ + typedef unsigned char Weights_type; //from 0 t 255 + const std::size_t img_size = image.size(); + + //create weights image + _image* weights + = _createImage(image.xdim(), image.ydim(), image.zdim(), + 1, //vectorial dimension + image.vx(), image.vy(), image.vz(), + sizeof(Weights_type), //image word size in bytes + internal::get_wordkind(), //image word kind WK_FIXED, WK_FLOAT, WK_UNKNOWN + internal::get_sign()); //image word sign + Weights_type* weights_ptr = (Weights_type*)(weights->data); + std::fill(weights_ptr, + weights_ptr + img_size, + Weights_type(0)); + weights->tx = image.tx(); + weights->ty = image.ty(); + weights->tz = image.tz(); + + //convert image to itkImage + using ImageType = itk::Image; + using WeightsType = itk::Image; + typename ImageType::Pointer itk_img = ImageType::New(); + std::set labels; + internal::convert_image_3_to_itk(image, itk_img.GetPointer(), labels); + + using DuplicatorType = itk::ImageDuplicator; + using IndicatorFilter = itk::BinaryThresholdImageFilter; + using GaussianFilterType = itk::RecursiveGaussianImageFilter; + using MaximumImageFilterType = itk::MaximumImageFilter; + + std::vector indicators(labels.size()); + typename DuplicatorType::Pointer duplicator = DuplicatorType::New(); + duplicator->SetInputImage(itk_img); + duplicator->Update(); + + for (std::size_t id = 0; id < labels.size(); ++id) + { + if (id > 0) + { + duplicator->SetInputImage(indicators[id - 1]); + duplicator->Update(); + } + indicators[id] = duplicator->GetOutput(); + } + + int id = 0; + typename WeightsType::Pointer blured_max = WeightsType::New(); + for (Image_word_type label : labels) + { +#ifdef CGAL_MESH_3_WEIGHTED_IMAGES_DEBUG + std::cout << "\nLABEL = " << label << std::endl; +#endif + + //compute "indicator image" for "label" + typename IndicatorFilter::Pointer indicator = IndicatorFilter::New(); + indicator->SetInput(indicators[id]); + indicator->SetOutsideValue(0); + indicator->SetInsideValue(255); + indicator->SetLowerThreshold(label); + indicator->SetUpperThreshold(label); + indicator->Update(); + + //perform gaussian smoothing + typename GaussianFilterType::Pointer smoother = GaussianFilterType::New(); + smoother->SetInput(indicator->GetOutput()); + smoother->SetSigma(sigma); + smoother->Update(); + + //take the max of smoothed indicator functions + if (id == 0) + blured_max = smoother->GetOutput(); + else + { + typename MaximumImageFilterType::Pointer maximumImageFilter = MaximumImageFilterType::New(); + maximumImageFilter->SetInput(0, blured_max); + maximumImageFilter->SetInput(1, smoother->GetOutput()); + maximumImageFilter->Update(); + blured_max = maximumImageFilter->GetOutput(); + } + + 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; +#endif + } + + //copy pixels to weights + std::copy(blured_max->GetBufferPointer(), + blured_max->GetBufferPointer() + img_size, + weights_ptr); + + CGAL::Image_3 weights_img(weights); + +#ifdef CGAL_MESH_3_WEIGHTED_IMAGES_DEBUG + std::cout << "non white in image \t= " + << internal::count_non_white_pixels(image) << std::endl; + std::cout << "non white in weights \t= " + << internal::count_non_white_pixels(weights_img) << std::endl; + std::cout << "non white in itkWeights \t= " + << internal::count_non_white_pixels(blured_max.GetPointer()) << std::endl; + _writeImage(weights, "weights-image.inr.gz"); +#endif + + return weights_img; +} +/// @endcond + +/*! +* Free function that generates a `CGAL::Image_3` of weights associated to each +* voxel of `image`, to make the output mesh surfaces smoother. +* The weights image is generated using the algorithm described by Stalling et al +* in \cgalCite{stalling1998weighted}. +* The [Insight toolkit](https://itk.org/) is needed to compile this function. +* +* @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 +* +* @returns a `CGAL::Image_3` of weights used to build a quality `Labeled_mesh_domain_3`, +* with the same dimensions as `image` +*/ + +CGAL::Image_3 generate_label_weights(const CGAL::Image_3& image, + const float& sigma) +{ + CGAL_IMAGE_IO_CASE(image.image(), + return generate_label_weights_with_known_word_type(image, sigma); + ); + CGAL_error_msg("This place should never be reached, because it would mean " + "the image word type is a type that is not handled by " + "CGAL_ImageIO."); + return CGAL::Image_3(); +} + +}//namespace Mesh_3 +}//namespace CGAL + +#endif // CGAL_MESH_3_GENERATE_LABEL_WEIGHTS_H diff --git a/Polyhedron/demo/Polyhedron/CMakeLists.txt b/Polyhedron/demo/Polyhedron/CMakeLists.txt index ade09e4d8369..a5da27ffd13b 100644 --- a/Polyhedron/demo/Polyhedron/CMakeLists.txt +++ b/Polyhedron/demo/Polyhedron/CMakeLists.txt @@ -268,6 +268,7 @@ if(CGAL_Qt5_FOUND AND Qt5_FOUND) endif() add_item(scene_tetrahedra_item Scene_tetrahedra_item.cpp) target_link_libraries(scene_tetrahedra_item PUBLIC scene_c3t3_item) + add_item(scene_transform_item Plugins/PCA/Scene_facegraph_transform_item.cpp) add_item(scene_edit_box_item Plugins/PCA/Scene_edit_box_item.cpp) add_item(scene_image_item Scene_image_item.cpp) diff --git a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/CMakeLists.txt b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/CMakeLists.txt index ac143f6dec03..89d05476c462 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/CMakeLists.txt +++ b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/CMakeLists.txt @@ -28,6 +28,12 @@ target_link_libraries( ${OPENGL_gl_LIBRARY}) target_include_directories(mesh_3_plugin PRIVATE include) +find_package(ITK QUIET COMPONENTS ITKCommon ITKThresholding ITKSmoothing ITKImageIntensity) +if(ITK_FOUND) + include(CGAL_ITK_support) + target_link_libraries(mesh_3_plugin PUBLIC CGAL::ITK_support) +endif(ITK_FOUND) + find_package(VTK QUIET COMPONENTS vtkImagingGeneral vtkIOImage NO_MODULE) if(VTK_FOUND) if(VTK_USE_FILE) diff --git a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_3_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_3_plugin.cpp index 5dcf518f0ec2..e2139eb1b14a 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_3_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_3_plugin.cpp @@ -36,6 +36,10 @@ auto make_not_null(T&& t) { #ifdef CGAL_MESH_3_DEMO_ACTIVATE_SEGMENTED_IMAGES #include "Scene_image_item.h" #include "Image_type.h" +#ifdef CGAL_USE_ITK +#include +#endif + #endif #include "Meshing_thread.h" @@ -590,6 +594,29 @@ void Mesh_3_plugin::mesh_3(const Mesh_type mesh_type, } } } + + //Labeled (weighted) image + connect(ui.useWeights_checkbox, SIGNAL(toggled(bool)), + ui.weightsSigma, SLOT(setEnabled(bool))); + connect(ui.useWeights_checkbox, SIGNAL(toggled(bool)), + ui.weightsSigma_label, SLOT(setEnabled(bool))); + ui.weightsSigma->setValue(1.); + bool input_is_labeled_img = (image_item != nullptr && !image_item->isGray()); + ui.labeledImgGroup->setVisible(input_is_labeled_img); + +#ifndef CGAL_USE_ITK + if (input_is_labeled_img) + { + ui.labeledImgGroup->setDisabled(true); + ui.labeledImgGroup->setToolTip( + QString("The use of weighted images is disabled " + "because the Insight Toolkit (ITK) is not available.")); + ui.useWeights_checkbox->setDisabled(true); + ui.weightsSigma_label->setDisabled(true); + ui.weightsSigma->setDisabled(true); + } +#endif + // ----------------------------------- // Get values // ----------------------------------- @@ -621,6 +648,9 @@ void Mesh_3_plugin::mesh_3(const Mesh_type mesh_type, const float iso_value = float(ui.iso_value_spinBox->value()); const float value_outside = float(ui.value_outside_spinBox->value()); const bool inside_is_less = ui.inside_is_less_checkBox->isChecked(); + const float sigma_weights = ui.useWeights_checkbox->isChecked() + ? ui.weightsSigma->value() : 0.f; + as_facegraph = (mesh_type == Mesh_type::SURFACE_ONLY) ? ui.facegraphCheckBox->isChecked() : false; @@ -692,6 +722,18 @@ void Mesh_3_plugin::mesh_3(const Mesh_type mesh_type, QMessageBox::critical(mw, tr(""), tr("ERROR: no data in selected item")); return; } +#ifdef CGAL_USE_ITK + if ( sigma_weights > 0 + && sigma_weights != image_item->sigma_weights()) + { + image_item->set_image_weights( + CGAL::Mesh_3::generate_label_weights(*pImage, sigma_weights), + sigma_weights); + } +#endif + const Image* pWeights = sigma_weights > 0 + ? image_item->image_weights() + : nullptr; Scene_polylines_item::Polylines_container plc; @@ -711,7 +753,8 @@ void Mesh_3_plugin::mesh_3(const Mesh_type mesh_type, image_item->isGray(), iso_value, value_outside, - inside_is_less); + inside_is_less, + pWeights); break; } default: diff --git a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_3_plugin_cgal_code.cpp b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_3_plugin_cgal_code.cpp index 4eee5f1aa9ff..f51c46022ed4 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_3_plugin_cgal_code.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_3_plugin_cgal_code.cpp @@ -203,7 +203,8 @@ Meshing_thread* cgal_code_mesh_3(const Image* pImage, bool is_gray, float iso_value, float value_outside, - bool inside_is_less) + bool inside_is_less, + const Image* pWeights) { if (nullptr == pImage) { return nullptr; } @@ -221,12 +222,50 @@ Meshing_thread* cgal_code_mesh_3(const Image* pImage, param.tet_shape = tet_shape; param.manifold = manifold; param.image_3_ptr = pImage; + param.weights_ptr = pWeights; Scene_c3t3_item* p_new_item = new Scene_c3t3_item(surface_only); + + QString tooltip = QString("\" With the following mesh parameters" + "
      " + "
    • Angle: %1
      • " + "
    • Edge size bound: %2
      • " + "
    • Facets size bound: %3
      • " + "
    • Approximation bound: %4
      • ") + .arg(facet_angle) + .arg(edge_size) + .arg(facet_sizing) + .arg(facet_approx); + if (!surface_only) + tooltip += QString("
    • Tetrahedra size bound: %1
      • ") + .arg(tet_sizing); + tooltip += QString("
    • Use Weighted Image: %1
      • ") + .arg(pWeights == nullptr ? "No" : "Yes"); + tooltip += "
    "; + + p_new_item->setProperty("toolTip", tooltip); + if(!is_gray) { namespace p = CGAL::parameters; - Image_mesh_domain* p_domain = new Image_mesh_domain + Image_mesh_domain* p_domain; +#ifdef CGAL_USE_ITK + if(nullptr != pWeights) + { + p_domain = new Image_mesh_domain + (Image_mesh_domain::create_labeled_image_mesh_domain + (p::image = *pImage, + p::weights = *pWeights, + p::relative_error_bound = 1e-6, + p::construct_surface_patch_index = + [](int i, int j) { return (i * 1000 + j); } + ) + ); + } + else +#endif + { + p_domain = new Image_mesh_domain (Image_mesh_domain::create_labeled_image_mesh_domain (p::image = *pImage, p::relative_error_bound = 1e-6, @@ -234,6 +273,7 @@ Meshing_thread* cgal_code_mesh_3(const Image* pImage, [](int i, int j) { return (i * 1000 + j); } ) ); + } if(protect_features && polylines.empty()){ std::vector > polylines_on_bbox; diff --git a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_3_plugin_cgal_code.h b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_3_plugin_cgal_code.h index 2e72edebc9c8..5a0e5cc0a1d1 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_3_plugin_cgal_code.h +++ b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_3_plugin_cgal_code.h @@ -63,5 +63,6 @@ Meshing_thread* cgal_code_mesh_3(const CGAL::Image_3* pImage, bool is_gray = false, float iso_value = 3.f, float value_outside = 0.f, - bool inside_is_less = true); + bool inside_is_less = true, + const CGAL::Image_3* pWeights = nullptr); #endif diff --git a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_function.h b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_function.h index addbc6b8ea7f..d5249e256b91 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_function.h +++ b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_function.h @@ -53,6 +53,7 @@ struct Mesh_parameters bool detect_connected_components; int manifold; const CGAL::Image_3* image_3_ptr; + const CGAL::Image_3* weights_ptr; bool use_sizing_field_with_aabb_tree; inline QStringList log() const; diff --git a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Meshing_dialog.ui b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Meshing_dialog.ui index cc9a9b224bdf..776adbc87661 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Meshing_dialog.ui +++ b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Meshing_dialog.ui @@ -9,8 +9,8 @@ 0 0 - 576 - 928 + 568 + 1049 @@ -342,10 +342,10 @@ - 25,00 + - 25,00 + 25.00 @@ -516,27 +516,69 @@ - - - Create as FaceGraph + + + Labeled Images + + + + + Qt::RightToLeft + + + Use Weighted Image + + + Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter + + + + + + + Qt::LeftToRight + + + + + + false + + + + + + + false + + + standard deviation parameter of the internal Gaussian filter + + + Sigma + + + Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter + + + + + + + false + + + + - - - Qt::Vertical - - - QSizePolicy::MinimumExpanding - - - - 20 - 0 - + + + Create as FaceGraph - + @@ -558,13 +600,31 @@ + detectComponents + protectEdges + protect + sharpEdgesAngle + edgeSizing + noEdgeSizing + approx noApprox + facetSizing noFacetSizing + facetAngle noAngle + tetSizing noTetSizing + tetShape noTetShape + value_outside_spinBox + iso_value_spinBox + inside_is_less_checkBox + useWeights_checkbox + weightsSigma + facegraphCheckBox checkBox manifoldCheckBox + facetTopology diff --git a/Polyhedron/demo/Polyhedron/Scene_image_item.cpp b/Polyhedron/demo/Polyhedron/Scene_image_item.cpp index 64f756237176..229aff8a7ffe 100644 --- a/Polyhedron/demo/Polyhedron/Scene_image_item.cpp +++ b/Polyhedron/demo/Polyhedron/Scene_image_item.cpp @@ -429,6 +429,8 @@ struct Scene_image_item_priv Scene_image_item_priv(int display_scale, bool hidden, Scene_image_item* parent) : m_initialized(false) , m_voxel_scale(display_scale) + , m_image_weights() + , m_sigma_weights(0.f) { item = parent; is_ogl_4_3 = Three::mainViewer()->isOpenGL_4_3(); @@ -450,6 +452,8 @@ struct Scene_image_item_priv bool is_ogl_4_3; Scene_image_item* item; internal::Vertex_buffer_helper* helper; + Image m_image_weights; + float m_sigma_weights; //#endif // SCENE_SEGMENTED_IMAGE_GL_BUFFERS_AVAILABLE }; @@ -496,6 +500,24 @@ Scene_image_item::compute_bbox() const m_image->image()->tz+(m_image->zdim()-1) * m_image->vz())); } +const Image* +Scene_image_item::image_weights() const +{ + return &d->m_image_weights; +} +void +Scene_image_item::set_image_weights(const Image& img, const float sigma) +{ + d->m_image_weights = img; + d->m_sigma_weights = sigma; +} +float +Scene_image_item::sigma_weights() const +{ + return d->m_sigma_weights; +} + + void Scene_image_item::draw(Viewer_interface* viewer) const { diff --git a/Polyhedron/demo/Polyhedron/Scene_image_item.h b/Polyhedron/demo/Polyhedron/Scene_image_item.h index 39630ae6d39a..f1470860c0b4 100644 --- a/Polyhedron/demo/Polyhedron/Scene_image_item.h +++ b/Polyhedron/demo/Polyhedron/Scene_image_item.h @@ -44,6 +44,11 @@ class SCENE_IMAGE_ITEM_EXPORT Scene_image_item const Image* image() const { return m_image; } bool isGray(); Image* m_image; + + const Image* image_weights() const; + void set_image_weights(const Image& img, const float sigma); + float sigma_weights() const; + void invalidateOpenGLBuffers(); void initializeBuffers(Viewer_interface *) const; void computeElements() const; diff --git a/Testsuite/test/collect_cgal_testresults_from_cmake b/Testsuite/test/collect_cgal_testresults_from_cmake index 54b5d502ce0b..f6aca34a87d1 100755 --- a/Testsuite/test/collect_cgal_testresults_from_cmake +++ b/Testsuite/test/collect_cgal_testresults_from_cmake @@ -60,7 +60,7 @@ print_testresult() RESULT="t" fi else - if grep -E -q 'NOTICE: .*(need|require|incompatible).*and.*will not be' CompilerOutput_$1 + if grep -E -q 'NOTICE: .*(need|require|incompatible).*will not be' CompilerOutput_$1 then RESULT="r" else