diff --git a/.gitignore b/.gitignore index 4c4769bf1ab..1bb95887866 100644 --- a/.gitignore +++ b/.gitignore @@ -373,6 +373,7 @@ Mesh_3/examples/Mesh_3/*.mesh.* Mesh_3/examples/Mesh_3/*.off Mesh_3/examples/Mesh_3/*.png Mesh_3/examples/Mesh_3/.*.deps +Mesh_3/examples/Mesh_3/random-image.inr Mesh_3/examples/Mesh_3/Makefile Mesh_3/examples/Mesh_3/applications Mesh_3/examples/Mesh_3/cgal_test_with_cmake diff --git a/CGAL_ImageIO/include/CGAL/ImageIO.h b/CGAL_ImageIO/include/CGAL/ImageIO.h index c296b99aa39..70b4ee109e5 100644 --- a/CGAL_ImageIO/include/CGAL/ImageIO.h +++ b/CGAL_ImageIO/include/CGAL/ImageIO.h @@ -614,7 +614,7 @@ static_evaluate(const _image* image, template inline -Word +Word& static_evaluate(const _image* image, const std::size_t i, const std::size_t j, @@ -623,6 +623,15 @@ static_evaluate(const _image* image, return ((Word*)image->data)[(k * image->ydim + j) * image->xdim + i]; } +template +inline +Word& +static_evaluate(const _image* image, + const std::size_t i) +{ + return ((Word*)image->data)[i]; +} + } // end namespace IMAGEIO } // end namespace CGAL diff --git a/Mesh_3/doc/Mesh_3/Doxyfile.in b/Mesh_3/doc/Mesh_3/Doxyfile.in index ed15cf7f4eb..fb3d336a944 100644 --- a/Mesh_3/doc/Mesh_3/Doxyfile.in +++ b/Mesh_3/doc/Mesh_3/Doxyfile.in @@ -6,4 +6,8 @@ HTML_EXTRA_FILES = ${CGAL_PACKAGE_DOC_DIR}/fig/implicit_domain_3.jpg ${CGAL_PACKAGE_DOC_DIR}/fig/implicit_domain_5.jpg \ ${CGAL_PACKAGE_DOC_DIR}/fig/no-protection.png \ ${CGAL_PACKAGE_DOC_DIR}/fig/protection-box.png \ - ${CGAL_PACKAGE_DOC_DIR}/fig/protection-all.png + ${CGAL_PACKAGE_DOC_DIR}/fig/protection-all.png \ + ${CGAL_PACKAGE_DOC_DIR}/fig/no-custom-init.png \ + ${CGAL_PACKAGE_DOC_DIR}/fig/with-custom-init.png + +EXAMPLE_PATH += ${CGAL_PACKAGE_INCLUDE_DIR} diff --git a/Mesh_3/doc/Mesh_3/Mesh_3.txt b/Mesh_3/doc/Mesh_3/Mesh_3.txt index 7c375f4cd52..90f27f8e693 100644 --- a/Mesh_3/doc/Mesh_3/Mesh_3.txt +++ b/Mesh_3/doc/Mesh_3/Mesh_3.txt @@ -669,6 +669,91 @@ 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_3DomainsFromSegmented3DImagesWithCustomInitialization Domains From Segmented 3D Images, with a Custom Initialization + +The example \ref Mesh_3/mesh_3D_image_with_custom_initialization.cpp is a modification +of \ref Mesh_3/mesh_3D_image.cpp. The goal of that example is to show how +the default initialization of the triangulation, using random rays, can be +replaced by a new implementation. In this case, the initialization detects +all connected components in the 3D segmented image, and inserts points in +the triangulation for each connected component. + +For the meshing, in the previous example (\ref Mesh_3/mesh_3D_image.cpp), we called `make_mesh_3()` as follows. + +\snippet Mesh_3/mesh_3D_image.cpp Meshing + +In the example \ref Mesh_3/mesh_3D_image_with_custom_initialization.cpp, +that call is replaced by: + -# the creation of an empty `%c3t3` object, + -# a call to a non-documented function + `initialize_triangulation_from_labeled_image()` that inserts points in + the triangulation, + -# then the call to `refine_mesh_3()`. + +\snippet Mesh_3/mesh_3D_image_with_custom_initialization.cpp Meshing + +The code of the function `initialize_triangulation_from_labeled_image()` is +in the non-documented header \ref +CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h\. As it is +undocumented and may be removed or modified at any time, if you wish to +use it then you should copy-paste it to your user code. The code of that +function is rather complicated. The following lines show how to insert new +points in the `%c3t3` object, with the calls to +`MeshVertexBase_3::set_dimension()` and +`MeshVertexBase_3::set_index()`. + +\snippet CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h insert initial points + +The value of `index` must be consistent with the possible values of +`Mesh_domain::Index`. In \ref +CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h, it is +constructed using the API of the mesh domain, as follows. First the functor +`construct_intersect` is created + +\dontinclude CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h +\skip Construct_intersection construct_intersection = +\until construct_intersection_object +then the `%Mesh_domain::Intersection` object (a `%tuple` with three +elements) is constructed using a call to the functor `construct_intersection` +\skip Intersection intersect = +\until construct_intersection +and eventually `%index` is the element \#1 of `%intersect`. +\skipline get<1>(intersect) + +The result of the custom initialization can be seen in +\cgalFigureRef{mesh3custominitimage3D}. The generated 3D image contains a +big sphere at the center, and 50 smaller spheres, generated +randomly. Without the custom initialization, only the biggest component +(the sphere at the center) was initialized and meshed. With the custom +initialization, the initial `%c3t3` object contains points on all connected +components, and all spheres are meshed. + +\cgalFigureAnchor{mesh3custominitimage3D} +
+ + + + + +
+
+\cgalFigureCaptionBegin{mesh3custominitimage3D} +Left: the mesh without the custom initialization, only the big sphere at +the center is meshed + +Right: the mesh generated after the initialization of all connected components +\cgalFigureCaptionEnd + + +Note that the example \ref +Mesh_3/mesh_3D_image_with_custom_initialization.cpp also shows how to +create a 3D image using the undocumented API of CGAL_ImageIO. + +\snippet Mesh_3/mesh_3D_image_with_custom_initialization.cpp Create the image + +The code of the function `%random_labeled_image()` is in the header file \ref +Mesh_3/random_labeled_image.h\. + \subsection Mesh_3UsingVariableSizingField Using Variable Sizing Field \subsubsection Mesh_3SizingFieldasanAnalyticalFunction Sizing Field as an Analytical Function diff --git a/Mesh_3/doc/Mesh_3/examples.txt b/Mesh_3/doc/Mesh_3/examples.txt index 682b8406c87..3e4dfaaf201 100644 --- a/Mesh_3/doc/Mesh_3/examples.txt +++ b/Mesh_3/doc/Mesh_3/examples.txt @@ -2,6 +2,9 @@ \example Mesh_3/implicit_functions.cpp \example Mesh_3/mesh_3D_image.cpp \example Mesh_3/mesh_3D_image_with_features.cpp +\example Mesh_3/mesh_3D_image_with_custom_initialization.cpp +\example Mesh_3/random_labeled_image.h +\example CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h \example Mesh_3/mesh_3D_image_variable_size.cpp \example Mesh_3/mesh_implicit_domains.cpp \example Mesh_3/mesh_implicit_domains_2.cpp diff --git a/Mesh_3/doc/Mesh_3/fig/no-custom-init.png b/Mesh_3/doc/Mesh_3/fig/no-custom-init.png new file mode 100644 index 00000000000..41c313519fe Binary files /dev/null and b/Mesh_3/doc/Mesh_3/fig/no-custom-init.png differ diff --git a/Mesh_3/doc/Mesh_3/fig/with-custom-init.png b/Mesh_3/doc/Mesh_3/fig/with-custom-init.png new file mode 100644 index 00000000000..1ef18a62e1e Binary files /dev/null and b/Mesh_3/doc/Mesh_3/fig/with-custom-init.png differ diff --git a/Mesh_3/examples/Mesh_3/CMakeLists.txt b/Mesh_3/examples/Mesh_3/CMakeLists.txt index 8a1535a53b9..1ea85bed952 100644 --- a/Mesh_3/examples/Mesh_3/CMakeLists.txt +++ b/Mesh_3/examples/Mesh_3/CMakeLists.txt @@ -86,6 +86,7 @@ if ( CGAL_FOUND ) create_single_source_cgal_program( "mesh_optimization_example.cpp" ) create_single_source_cgal_program( "mesh_optimization_lloyd_example.cpp" ) create_single_source_cgal_program( "mesh_3D_image.cpp" ) + create_single_source_cgal_program( "mesh_3D_image_with_custom_initialization.cpp" ) create_single_source_cgal_program( "mesh_3D_image_variable_size.cpp" ) 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." ) diff --git a/Mesh_3/examples/Mesh_3/mesh_3D_image.cpp b/Mesh_3/examples/Mesh_3/mesh_3D_image.cpp index 6893f5d9594..8fa9d544102 100644 --- a/Mesh_3/examples/Mesh_3/mesh_3D_image.cpp +++ b/Mesh_3/examples/Mesh_3/mesh_3D_image.cpp @@ -1,4 +1,3 @@ - #include #include @@ -32,13 +31,14 @@ using namespace CGAL::parameters; int main(int argc, char* argv[]) { + /// [Loads image] const char* fname = (argc>1)?argv[1]:"data/liver.inr.gz"; - // Loads image CGAL::Image_3 image; if(!image.read(fname)){ std::cerr << "Error: Cannot read file " << fname << std::endl; return EXIT_FAILURE; } + /// [Loads image] // Domain Mesh_domain domain(image); @@ -47,8 +47,9 @@ int main(int argc, char* argv[]) Mesh_criteria criteria(facet_angle=30, facet_size=6, facet_distance=4, cell_radius_edge_ratio=3, cell_size=8); - // Meshing + /// [Meshing] C3t3 c3t3 = CGAL::make_mesh_3(domain, criteria); + /// [Meshing] // Output std::ofstream medit_file("out.mesh"); diff --git a/Mesh_3/examples/Mesh_3/mesh_3D_image_with_custom_initialization.cpp b/Mesh_3/examples/Mesh_3/mesh_3D_image_with_custom_initialization.cpp new file mode 100644 index 00000000000..9197c4239bf --- /dev/null +++ b/Mesh_3/examples/Mesh_3/mesh_3D_image_with_custom_initialization.cpp @@ -0,0 +1,64 @@ +#include "random_labeled_image.h" +#include + +#include +#include +#include + +#include + +#include +#include +#include + +#include + +// Domain +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Labeled_image_mesh_domain_3 Mesh_domain; + +#ifdef CGAL_CONCURRENT_MESH_3 +typedef CGAL::Parallel_tag Concurrency_tag; +#else +typedef CGAL::Sequential_tag Concurrency_tag; +#endif + +// Triangulation +typedef CGAL::Mesh_triangulation_3::type Tr; + +typedef CGAL::Mesh_complex_3_in_triangulation_3 C3t3; + +// Criteria +typedef CGAL::Mesh_criteria_3 Mesh_criteria; + +// To avoid verbose function and named parameters call +using namespace CGAL::parameters; + +int main() +{ + /// [Create the image] + CGAL::Image_3 image = random_labeled_image(); + /// [Create the image] + + // Domain + Mesh_domain domain(image); + + // Mesh criteria + Mesh_criteria criteria(facet_angle=30, facet_size=3, facet_distance=1, + cell_radius_edge_ratio=3, cell_size=3); + + /// [Meshing] + C3t3 c3t3; + initialize_triangulation_from_labeled_image(c3t3, + domain, + image, + criteria, + (unsigned char)0); + CGAL::refine_mesh_3(c3t3, domain, criteria); + /// [Meshing] + + // Output + CGAL::dump_c3t3(c3t3, "out"); + + return 0; +} diff --git a/Mesh_3/examples/Mesh_3/random_labeled_image.h b/Mesh_3/examples/Mesh_3/random_labeled_image.h new file mode 100644 index 00000000000..8d03d462d49 --- /dev/null +++ b/Mesh_3/examples/Mesh_3/random_labeled_image.h @@ -0,0 +1,53 @@ +#include +#include +#include + +CGAL::Image_3 random_labeled_image() +{ + const int dim = 400; + const unsigned char number_of_spheres = 50; + const int max_radius_of_spheres = 10; + const int radius_of_big_sphere = 80; + _image* image = _createImage(dim, dim, dim, 1, + 1.f, 1.f, 1.f, 1, + WK_FIXED, SGN_UNSIGNED); + unsigned char* ptr = (unsigned char*)(image->data); + std::fill(ptr, ptr+dim*dim*dim, '\0'); + + std::ptrdiff_t center = dim / 2; + CGAL::Random rand(0); + for(unsigned char n = number_of_spheres; n > 0 ; --n) { + std::size_t i, j, k; + do { + i = rand.uniform_smallint(1 + max_radius_of_spheres, + dim-2 - max_radius_of_spheres); + j = rand.uniform_smallint(1 + max_radius_of_spheres, + dim-2 - max_radius_of_spheres); + k = rand.uniform_smallint(1 + max_radius_of_spheres, + dim-2 - max_radius_of_spheres); + } while ( ( CGAL::square(double(center) - double(i)) + + CGAL::square(double(center) - double(j)) + + CGAL::square(double(center) - double(k)) ) + < + CGAL::square(double(radius_of_big_sphere) + 4 * max_radius_of_spheres) ); + std::ptrdiff_t radius = max_radius_of_spheres; + if(n==1) { + i = j = k = center; + radius = radius_of_big_sphere; + } + for(std::ptrdiff_t ii = - radius; ii <= radius; ++ii) + { + for(std::ptrdiff_t jj = - radius; jj <= radius; ++jj) + { + for(std::ptrdiff_t kk = - radius; kk <= radius; ++kk) + { + if(ii*ii + jj*jj + kk*kk > radius * radius) continue; + using CGAL::IMAGEIO::static_evaluate; + static_evaluate(image, i+ii, j+jj, k+kk) = n; + } + } + } + } + _writeImage(image, "random-image.inr"); + return CGAL::Image_3(image); +} diff --git a/Mesh_3/include/CGAL/Labeled_image_mesh_domain_3.h b/Mesh_3/include/CGAL/Labeled_image_mesh_domain_3.h index 982b3fa023f..f7a62f59562 100644 --- a/Mesh_3/include/CGAL/Labeled_image_mesh_domain_3.h +++ b/Mesh_3/include/CGAL/Labeled_image_mesh_domain_3.h @@ -32,6 +32,7 @@ #include #include #include +#include namespace CGAL { diff --git a/Mesh_3/include/CGAL/Labeled_mesh_domain_3.h b/Mesh_3/include/CGAL/Labeled_mesh_domain_3.h index af50ca3245f..655bad3442f 100644 --- a/Mesh_3/include/CGAL/Labeled_mesh_domain_3.h +++ b/Mesh_3/include/CGAL/Labeled_mesh_domain_3.h @@ -39,6 +39,10 @@ #include #include +#include + +#include + namespace CGAL { struct Null_subdomain_index { diff --git a/Mesh_3/include/CGAL/Mesh_3/Dump_c3t3.h b/Mesh_3/include/CGAL/Mesh_3/Dump_c3t3.h index aba5b56c898..85e830cd6a0 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Dump_c3t3.h +++ b/Mesh_3/include/CGAL/Mesh_3/Dump_c3t3.h @@ -32,9 +32,13 @@ namespace CGAL { template ::value && - is_streamable::value && - is_streamable::value && - is_streamable::value + is_streamable::value + && + (is_streamable::value || + Output_rep::is_specialized) + && + (is_streamable::value || + Output_rep::is_specialized) > struct Dump_c3t3 { void dump_c3t3(const C3t3& c3t3, std::string prefix) const { @@ -75,13 +79,17 @@ struct Dump_c3t3 { << "\n"; } - if(!is_streamable::value) { + if(!is_streamable::value && + !CGAL::Output_rep::is_specialized) + { std::cerr << " - C3t3::Surface_patch_index is not streamable\n"; std::cerr << " " << typeid(typename C3t3::Surface_patch_index).name() << "\n"; } - if(!is_streamable::value) { + if(!is_streamable::value && + !CGAL::Output_rep::is_specialized) + { std::cerr << " - C3t3::Subdomain_index is not streamable\n"; std::cerr << " " << typeid(typename C3t3::Subdomain_index).name() diff --git a/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h b/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h new file mode 100644 index 00000000000..db729c8ec6d --- /dev/null +++ b/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h @@ -0,0 +1,237 @@ +// Copyright (c) 2015,2016 GeometryFactory +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s) : Laurent Rineau + +#ifndef CGAL_MESH_3_INITIALIZE_TRIANGULATION_FROM_LABELED_IMAGE_H +#define CGAL_MESH_3_INITIALIZE_TRIANGULATION_FROM_LABELED_IMAGE_H + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +template +struct Get_point +{ + const double vx, vy, vz; + Get_point(const CGAL::Image_3* image) + : vx(image->vx()) + , vy(image->vy()) + , vz(image->vz()) + {} + + Point operator()(const std::size_t i, + const std::size_t j, + const std::size_t k) const + { + return Point(double(i) * vx, double(j) * vy, double(k) * vz); + } +}; +template +void init_tr_from_labeled_image_call_init_features(C3T3&, + const MeshDomain&, + const MeshCriteria&, + CGAL::Tag_false) +{ +} +template +void init_tr_from_labeled_image_call_init_features(C3T3& c3t3, + const MeshDomain& domain, + const MeshCriteria& criteria, + CGAL::Tag_true) +{ + CGAL::internal::Mesh_3::init_c3t3_with_features(c3t3, + domain, + criteria); + std::cout << c3t3.triangulation().number_of_vertices() + << " initial points on 1D-features" << std::endl; +} + + +template +void initialize_triangulation_from_labeled_image(C3T3& c3t3, + const MeshDomain& domain, + const CGAL::Image_3& image, + const MeshCriteria& criteria, + Image_word_type, + bool protect_features = false + ) +{ + typedef typename C3T3::Triangulation Tr; + typedef typename Tr::Point Weighted_point; + typedef typename Weighted_point::Point Point_3; + typedef typename Tr::Segment Segment_3; + typedef typename Tr::Geom_traits::Vector_3 Vector_3; + typedef typename Tr::Vertex_handle Vertex_handle; + typedef typename Tr::Cell_handle Cell_handle; + + typedef Point_3 Point; + typedef MeshDomain Mesh_domain; + + Tr& tr = c3t3.triangulation(); + + if(protect_features) { + init_tr_from_labeled_image_call_init_features + (c3t3, domain, criteria, + CGAL::internal::Mesh_3::Has_features()); + } + + const double max_v = (std::max)((std::max)(image.vx(), + image.vy()), + image.vz()); + + typedef std::vector > Seeds; + Seeds seeds; + Get_point get_point(&image); + std::cout << "Searching for connected components..." << std::endl; + CGAL::Identity no_transformation; + search_for_connected_components_in_labeled_image(image, + std::back_inserter(seeds), + CGAL::Emptyset_iterator(), + no_transformation, + get_point, + Image_word_type()); + std::cout << " " << seeds.size() << " components were found." << std::endl; + std::cout << "Construct initial points..." << std::endl; + for(typename Seeds::const_iterator it = seeds.begin(), end = seeds.end(); + it != end; ++it) + { + const double radius = double(it->second + 1)* max_v; + CGAL::Random_points_on_sphere_3 points_on_sphere_3(radius); + typename Mesh_domain::Construct_intersection construct_intersection = + domain.construct_intersection_object(); + + std::vector directions; + if(it->second < 2) { + // shoot in six directions + directions.push_back(Vector_3(-radius, 0, 0)); + directions.push_back(Vector_3(+radius, 0, 0)); + directions.push_back(Vector_3(0, -radius, 0)); + directions.push_back(Vector_3(0, +radius, 0)); + directions.push_back(Vector_3(0, 0, -radius)); + directions.push_back(Vector_3(0, 0, +radius)); + } else { + for(int i = 0; i < 20; ++i) + { + // shoot 20 random directions + directions.push_back(*points_on_sphere_3++ - CGAL::ORIGIN); + } + } + + BOOST_FOREACH(const Vector_3& v, directions) + { + const Point test = it->first + v; + const typename Mesh_domain::Intersection intersect = + construct_intersection(Segment_3(it->first, test)); + if (CGAL::cpp11::get<2>(intersect) != 0) + { + Point_3 pi = CGAL::cpp11::get<0>(intersect); + + // This would cause trouble to optimizers + // check pi will not be hidden + typename Tr::Locate_type lt; + Cell_handle c; + int li, lj; + Cell_handle pi_cell = tr.locate(pi, lt, li, lj); + if(lt != Tr::OUTSIDE_AFFINE_HULL) { + switch (tr.dimension()) + { //skip dimension 0 + case 1: + if (tr.side_of_power_segment(pi_cell, pi, true) != CGAL::ON_BOUNDED_SIDE) + continue; + break; + case 2: + if (tr.side_of_power_circle(pi_cell, 3, pi, true) != CGAL::ON_BOUNDED_SIDE) + continue; + break; + case 3: + if (tr.side_of_power_sphere(pi_cell, pi, true) != CGAL::ON_BOUNDED_SIDE) + continue; + } + } + + //check pi is not inside a protecting ball + std::vector conflict_vertices; + if (tr.dimension() == 3) + { + tr.vertices_on_conflict_zone_boundary(pi, pi_cell + , std::back_inserter(conflict_vertices)); + } + else + { + for (typename Tr::Finite_vertices_iterator vit = tr.finite_vertices_begin(); + vit != tr.finite_vertices_end(); ++vit) + { + if (vit->point().weight() > 0.) + conflict_vertices.push_back(vit); + } + } + bool pi_inside_protecting_sphere = false; + BOOST_FOREACH(Vertex_handle cv, conflict_vertices) + { + if (cv->point().weight() == 0.) + continue; + if (CGAL::compare_squared_distance(pi, cv->point(), cv->point().weight()) + != CGAL::LARGER) + { + pi_inside_protecting_sphere = true; + break; + } + } + if (pi_inside_protecting_sphere) + continue; + const typename Mesh_domain::Index index = CGAL::cpp11::get<1>(intersect); + + /// The following lines show how to insert initial points in the + /// `c3t3` object. [insert initial points] + Vertex_handle v = tr.insert(pi); + + // `v` could be null if `pi` is hidden by other vertices of `tr`. + CGAL_assertion(v != Vertex_handle()); + + c3t3.set_dimension(v, 2); // by construction, points are on surface + c3t3.set_index(v, index); + /// [insert initial points] + } + // else + // { + // std::cerr << + // boost::format("Error. Segment (%1%, %2%) does not intersect the surface!\n") + // % it->first % test; + // } + } + } + std::cout << " " << tr.number_of_vertices() << " initial points." << std::endl; + if ( c3t3.triangulation().dimension() != 3 ) + { + std::cout << " not enough points: triangulation.dimension() == " + << c3t3.triangulation().dimension() << std::endl; + CGAL::internal::Mesh_3::init_c3t3(c3t3, domain, criteria, 20); + std::cout << " -> " << tr.number_of_vertices() << " initial points." << std::endl; + } +} + +#endif // CGAL_MESH_3_INITIALIZE_TRIANGULATION_FROM_LABELED_IMAGE_H diff --git a/Mesh_3/include/CGAL/Mesh_3/search_for_connected_components_in_labeled_image.h b/Mesh_3/include/CGAL/Mesh_3/search_for_connected_components_in_labeled_image.h new file mode 100644 index 00000000000..205d0e8af1e --- /dev/null +++ b/Mesh_3/include/CGAL/Mesh_3/search_for_connected_components_in_labeled_image.h @@ -0,0 +1,249 @@ +// Copyright (c) 2015,2016 GeometryFactory +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s) : Laurent Rineau + +#ifndef CGAL_MESH_3_SEARCH_FOR_CONNECTED_COMPONENTS_IN_LABELED_IMAGE_H +#define CGAL_MESH_3_SEARCH_FOR_CONNECTED_COMPONENTS_IN_LABELED_IMAGE_H + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#ifdef CGAL_MESH_3_SEARCH_FOR_CONNECTED_COMPONENTS_IN_LABELED_IMAGE_VERBOSE +# include +# include +#endif // CGAL_MESH_3_SEARCH_FOR_CONNECTED_COMPONENTS_IN_LABELED_IMAGE_VERBOSE +template +void +search_for_connected_components_in_labeled_image(const CGAL::Image_3& image, + PointsOutputIterator it, + DomainsOutputIterator dom_it, + TransformOperator transform, + Construct_point point, + Image_word_type) +{ + const std::size_t nx = image.xdim(); + const std::size_t ny = image.ydim(); + const std::size_t nz = image.zdim(); + const std::size_t size = nx * ny * nz; + + typedef boost::uint16_t uint; + + if(nx > 65535 || ny > 65535 || nz > 65535) + { + CGAL_error_msg("The dimensions of the image must be lower than 2^16"); + } + + typedef typename TransformOperator::result_type Label; + + std::vector visited(size, false); + std::vector second_pass(size, false); + typedef boost::tuple Indices; + typedef std::queue > Indices_queue; + typedef std::vector Border_vector; + +#ifdef CGAL_MESH_3_SEARCH_FOR_CONNECTED_COMPONENTS_IN_LABELED_IMAGE_VERBOSE + int number_of_connected_components = 0; +#endif // CGAL_MESH_3_SEARCH_FOR_CONNECTED_COMPONENTS_IN_LABELED_IMAGE_VERBOSE + std::size_t voxel_index = 0; + for(uint k=0; k(image.image(), + voxel_index)); + *dom_it++ = current_label; + if(current_label == Label()) { + visited[voxel_index] = true; + second_pass[voxel_index] = true; + ++voxel_index; + continue; + } + +#ifdef CGAL_MESH_3_SEARCH_FOR_CONNECTED_COMPONENTS_IN_LABELED_IMAGE_VERBOSE + // if we reach here, (i, j, k) is a new connected component + ++number_of_connected_components; + std::cerr << boost::format("Found new connected component (#%5%) " + "at voxel (%1%, %2%, %3%), value=%4%, volume id=%6%\n") + % i % j % k + % (long)static_evaluate(image.image(), i, j, k) + % number_of_connected_components + % (int)current_label; +#endif // CGAL_MESH_3_SEARCH_FOR_CONNECTED_COMPONENTS_IN_LABELED_IMAGE_VERBOSE + + int nb_voxels = 0; + + Indices_queue queue; + Indices indices(i, j ,k, 0); + queue.push(indices); + + Border_vector border; + + /* + * First pass is a BFS to retrieve all the connected component, and + * its border. + * Second pass is a BFS initialized with all voxel of the border. + * The last voxel of that BFS is used as the seed. + */ + int pass = 1; // pass will be equal to 2 in second pass + + Indices bbox_min = indices; + Indices bbox_max = indices; + + while(!queue.empty()) // walk through the connected component + { + Indices indices = queue.front(); + queue.pop(); + + // warning: those indices i, j and k are local to the while loop + const uint i = boost::get<0>(indices); + const uint j = boost::get<1>(indices); + const uint k = boost::get<2>(indices); + const uint depth = boost::get<3>(indices); + + const size_t offset = i + nx * (j + ny * k); + const int m = (visited[offset] ? 1 : 0) + (second_pass[offset] ? 2 : 0); + if(m < pass) + { + if(pass == 1 ) + { + visited[offset] = true; + second_pass[offset] = false; + ++nb_voxels; + boost::get<0>(bbox_min) = (std::min)(i, boost::get<0>(bbox_min)); + boost::get<0>(bbox_max) = (std::max)(i, boost::get<0>(bbox_max)); + boost::get<1>(bbox_min) = (std::min)(j, boost::get<1>(bbox_min)); + boost::get<1>(bbox_max) = (std::max)(j, boost::get<1>(bbox_max)); + boost::get<2>(bbox_min) = (std::min)(k, boost::get<2>(bbox_min)); + boost::get<2>(bbox_max) = (std::max)(k, boost::get<2>(bbox_max)); + } else + { + CGAL_assertion(pass == 2); + visited[offset] = false; + second_pass[offset] = true; + } + + static const int neighbors_offset[6][3] = { { +1, 0, 0 }, + { -1, 0, 0 }, + { 0, +1, 0 }, + { 0, -1, 0 }, + { 0, 0, +1 }, + { 0, 0, -1 } }; + bool voxel_is_on_border = false; + + // Visit neighbors. + // (i_n, j_n, k_n) are indices of neighbors. + for(int n = 0; n < 6; ++n) + { + const ptrdiff_t i_n = i + neighbors_offset[n][0]; + const ptrdiff_t j_n = j + neighbors_offset[n][1]; + const ptrdiff_t k_n = k + neighbors_offset[n][2]; + if(i_n < 0 || i_n >= static_cast(nx) || + j_n < 0 || j_n >= static_cast(ny) || + k_n < 0 || k_n >= static_cast(nz)) + { + voxel_is_on_border = true; + continue; + } + else + { + const std::size_t offset_n = i_n + nx * (j_n + k_n * ny); + if(transform(static_evaluate(image.image(), + offset_n)) + == current_label) + { + const int m_n = (visited[offset_n] ? 1 : 0) + + (second_pass[offset_n] ? 2 : 0); + if(m_n < pass) { + Indices indices(uint(i_n), uint(j_n), uint(k_n), uint(depth+1)); + queue.push(indices); + } + } + else + voxel_is_on_border = true; + } + } // end for neighbors + + if(pass == 1 && voxel_is_on_border) + border.push_back(indices); + } // end if voxel not already visited + + if(queue.empty()) { + if(pass == 1) + { // End of first pass. Begin second pass with the voxels of + // the border. + for(typename Border_vector::const_iterator + border_it = border.begin(), border_end = border.end(); + border_it != border_end; ++border_it) + queue.push(*border_it); + pass = 2; + } + else // end of second pass, return the last visited voxel + { +// if(nb_voxels >= 100) + { + *it++ = std::make_pair(point(i, j, k), + depth+1); +#if CGAL_MESH_3_SEARCH_FOR_CONNECTED_COMPONENTS_IN_LABELED_IMAGE_VERBOSE > 1 + std::cerr << boost::format("Found seed %5%, which is voxel " + "(%1%, %2%, %3%), value=%4%\n") + % i % j % k + % (long)static_evaluate(image.image(), i, j, k) + % point(i, j, k); +#endif // CGAL_MESH_3_SEARCH_FOR_CONNECTED_COMPONENTS_IN_LABELED_IMAGE_VERBOSE>1 + } + } + } // end if queue.empty() + } // end while !queue.empty() (with local indices i, j, k) +#ifdef CGAL_MESH_3_SEARCH_FOR_CONNECTED_COMPONENTS_IN_LABELED_IMAGE_VERBOSE + std::cerr + << boost::format("There was %1% voxel(s) in that component.\n" + "The bounding box is (%2% %3% %4%, %5% %6% %7%).\n" + "%8% voxel(s) on border\n") + % nb_voxels + % boost::get<0>(bbox_min) % boost::get<1>(bbox_min) + % boost::get<2>(bbox_min) + % boost::get<0>(bbox_max) % boost::get<1>(bbox_max) + % boost::get<2>(bbox_max) + % border.size(); +#endif // CGAL_MESH_3_SEARCH_FOR_CONNECTED_COMPONENTS_IN_LABELED_IMAGE_VERBOSE + + ++voxel_index; + } // end for i,j,k +} // end function search_for_connected_components_in_labeled_image() + +#endif // CGAL_MESH_3_SEARCH_FOR_CONNECTED_COMPONENTS_IN_LABELED_IMAGE_H diff --git a/Mesh_3/include/CGAL/internal/Mesh_3/Handle_IO_for_pair_of_int.h b/Mesh_3/include/CGAL/internal/Mesh_3/Handle_IO_for_pair_of_int.h new file mode 100644 index 00000000000..777c47a6449 --- /dev/null +++ b/Mesh_3/include/CGAL/internal/Mesh_3/Handle_IO_for_pair_of_int.h @@ -0,0 +1,84 @@ +// Copyright (c) 2016 GeometryFactory +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s) : Laurent Rineau + +#ifndef CGAL_INTERNAL_MESH_3_INTERNAL_HANDLE_IO_FOR_PAIR_OF_INT_H +#define CGAL_INTERNAL_MESH_3_INTERNAL_HANDLE_IO_FOR_PAIR_OF_INT_H + +#include +#include +#include +#include + +namespace CGAL { +template <> +struct Get_io_signature > { + std::string operator()() const + { + return std::string("std::pair"); + } +}; // end Get_io_signature > + +inline std::ostream& operator<<(std::ostream& out, const std::pair& id) { + return out << id.first << " " << id.second; +} +inline std::istream& operator>>(std::istream& in, std::pair& id) { + return in >> id.first >> id.second; +} + +template <> +class Output_rep > : public IO_rep_is_specialized { + typedef std::pair T; + const T& t; +public: + //! initialize with a const reference to \a t. + Output_rep( const T& tt) : t(tt) {} + //! perform the output, calls \c operator\<\< by default. + std::ostream& operator()( std::ostream& out) const { + if(is_ascii(out)) { + out << t.first << " " << t.second; + } else { + CGAL::write(out, t.first); + CGAL::write(out, t.second); + } + return out; + } +}; + +template <> +class Input_rep > : public IO_rep_is_specialized { + typedef std::pair T; + T& t; +public: + //! initialize with a const reference to \a t. + Input_rep( T& tt) : t(tt) {} + //! perform the output, calls \c operator\<\< by default. + std::istream& operator()( std::istream& in) const { + if(is_ascii(in)) { + in >> t.first >> t.second; + } else { + CGAL::read(in, t.first); + CGAL::read(in, t.second); + } + return in; + } +}; +} // end namespace CGAL + +#endif // CGAL_INTERNAL_MESH_3_INTERNAL_HANDLE_IO_FOR_PAIR_OF_INT_H diff --git a/Mesh_3/include/CGAL/internal/Mesh_3/get_index.h b/Mesh_3/include/CGAL/internal/Mesh_3/get_index.h index 66eb3026449..25a58ff10a9 100644 --- a/Mesh_3/include/CGAL/internal/Mesh_3/get_index.h +++ b/Mesh_3/include/CGAL/internal/Mesh_3/get_index.h @@ -29,6 +29,7 @@ #include #include +#include namespace CGAL { namespace internal { @@ -87,13 +88,13 @@ struct Write_mesh_domain_index { switch(dimension) { case 0: { const Ci& ci = get_index(index); - if(is_ascii(os)) os << ci; + if(is_ascii(os)) os << oformat(ci); else CGAL::write(os, ci); } break; case 1: { const Si& si = get_index(index); - if(is_ascii(os)) os << si; + if(is_ascii(os)) os << oformat(si); else CGAL::write(os, si); } break; @@ -122,7 +123,7 @@ struct Read_mesh_domain_index { break; default: {// 3 typename MT::Subdomain_index di; - if(is_ascii(is)) is >> di; + if(is_ascii(is)) is >> iformat(di); else CGAL::read(is, di); return di; } @@ -146,13 +147,13 @@ struct Write_mesh_domain_index { switch(dimension) { case 2: { const Spi& spi = get_index(index); - if(is_ascii(os)) os << spi; + if(is_ascii(os)) os << oformat(spi); else CGAL::write(os, spi); } break; default: {// 3 const Di& di = get_index(index); - if(is_ascii(os)) os << di; + if(is_ascii(os)) os << oformat(di); else CGAL::write(os, di); } break; diff --git a/Mesh_3/include/CGAL/make_mesh_3.h b/Mesh_3/include/CGAL/make_mesh_3.h index 4e08d340c12..673066e7b65 100644 --- a/Mesh_3/include/CGAL/make_mesh_3.h +++ b/Mesh_3/include/CGAL/make_mesh_3.h @@ -257,7 +257,15 @@ struct C3t3_initializer < C3T3, MD, MC, true, CGAL::Tag_true > bool with_features, const int nb_initial_points = -1) { - if ( with_features ) { init_c3t3_with_features(c3t3,domain,criteria); } + if ( with_features ) { + init_c3t3_with_features(c3t3,domain,criteria); + + // If c3t3 initialization is not sufficient (may happen if there is only + // a planar curve as feature for example), add some surface points + if ( c3t3.triangulation().dimension() != 3 ) { + init_c3t3(c3t3, domain, criteria, nb_initial_points); + } + } else { init_c3t3(c3t3,domain,criteria,nb_initial_points); } } }; @@ -430,15 +438,8 @@ void make_mesh_3_impl(C3T3& c3t3, with_features, mesh_options.number_of_initial_points); - // If c3t3 initialization is not sufficient (may happen if there is only - // a planar curve as feature for example), add some surface points - if ( c3t3.triangulation().dimension() != 3 ) - { - internal::Mesh_3::init_c3t3(c3t3, domain, criteria, - mesh_options.number_of_initial_points); - } CGAL_assertion( c3t3.triangulation().dimension() == 3 ); - + // Build mesher and launch refinement process // Don't reset c3t3 as we just created it refine_mesh_3(c3t3, domain, criteria, diff --git a/Mesh_3/test/Mesh_3/test_meshing_3D_image.cpp b/Mesh_3/test/Mesh_3/test_meshing_3D_image.cpp index 9d8338d92ed..81b9d189cd1 100644 --- a/Mesh_3/test/Mesh_3/test_meshing_3D_image.cpp +++ b/Mesh_3/test/Mesh_3/test_meshing_3D_image.cpp @@ -46,8 +46,6 @@ public: typedef typename Mesh_criteria::Facet_criteria Facet_criteria; typedef typename Mesh_criteria::Cell_criteria Cell_criteria; - CGAL_USE_TYPE(typename Mesh_domain::Surface_patch_index); - //------------------------------------------------------- // Data generation //------------------------------------------------------- @@ -71,6 +69,10 @@ public: // Verify this->verify_c3t3_volume(c3t3, 1772330*0.95, 1772330*1.05); this->verify(c3t3,domain,criteria, Bissection_tag()); + + typedef typename Mesh_domain::Surface_patch_index Patch_id; + CGAL_static_assertion(CGAL::Output_rep::is_specialized); + CGAL_USE_TYPE(Patch_id); } }; diff --git a/Number_types/include/CGAL/CORE_BigRat.h b/Number_types/include/CGAL/CORE_BigRat.h index 0be849bec7e..6afc12f1fd7 100644 --- a/Number_types/include/CGAL/CORE_BigRat.h +++ b/Number_types/include/CGAL/CORE_BigRat.h @@ -155,7 +155,7 @@ public: }; template -class Output_rep< ::CORE::BigRat, F> { +class Output_rep< ::CORE::BigRat, F> : public IO_rep_is_specialized { const ::CORE::BigRat& t; public: //! initialize with a const reference to \a t. @@ -192,7 +192,9 @@ struct Needs_parens_as_product< ::CORE::BigRat >{ }; template <> -class Output_rep< ::CORE::BigRat, Parens_as_product_tag > { +class Output_rep< ::CORE::BigRat, Parens_as_product_tag > + : public IO_rep_is_specialized +{ const ::CORE::BigRat& t; public: // Constructor diff --git a/Number_types/include/CGAL/leda_rational.h b/Number_types/include/CGAL/leda_rational.h index f242f7637bd..903e3a14545 100644 --- a/Number_types/include/CGAL/leda_rational.h +++ b/Number_types/include/CGAL/leda_rational.h @@ -209,7 +209,7 @@ public: }; template -class Output_rep< leda_rational, F> { +class Output_rep< leda_rational, F> : public IO_rep_is_specialized { const leda_rational& t; public: //! initialize with a const reference to \a t. @@ -246,7 +246,9 @@ struct Needs_parens_as_product< leda_rational >{ }; template <> -class Output_rep< leda_rational, Parens_as_product_tag > { +class Output_rep< leda_rational, Parens_as_product_tag > + : public IO_rep_is_specialized +{ const leda_rational& t; public: // Constructor diff --git a/Number_types/include/CGAL/leda_real.h b/Number_types/include/CGAL/leda_real.h index b3503a86c07..b25aa39354b 100644 --- a/Number_types/include/CGAL/leda_real.h +++ b/Number_types/include/CGAL/leda_real.h @@ -212,7 +212,7 @@ template <> class Real_embeddable_traits< leda_real > template <> -class Output_rep< ::leda::real > { +class Output_rep< ::leda::real > : public IO_rep_is_specialized { const ::leda::real& t; public: //! initialize with a const reference to \a t. @@ -226,7 +226,9 @@ public: }; template <> -class Output_rep< ::leda::real, CGAL::Parens_as_product_tag > { +class Output_rep< ::leda::real, CGAL::Parens_as_product_tag > + : public IO_rep_is_specialized +{ const ::leda::real& t; public: //! initialize with a const reference to \a t. diff --git a/Number_types/test/Number_types/ioformat.cpp b/Number_types/test/Number_types/ioformat.cpp index c8f431898f7..cc270975a85 100644 --- a/Number_types/test/Number_types/ioformat.cpp +++ b/Number_types/test/Number_types/ioformat.cpp @@ -107,6 +107,7 @@ int main() // CORE #ifdef CGAL_USE_CORE + CGAL_static_assertion(CGAL::Output_rep::is_specialized == true); //bug in io for CORE. test_it("CORE::BigInt"); test_it("CORE::BigRat"); @@ -116,6 +117,7 @@ int main() // LEDA based NTs #ifdef CGAL_USE_LEDA + CGAL_static_assertion(CGAL::Output_rep::is_specialized == true); test_it("leda_integer"); test_it("leda_rational"); test_it("leda_bigfloat"); diff --git a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Image_res_dialog.ui b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Image_res_dialog.ui index 975d0d2962b..fe7a01cb767 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Image_res_dialog.ui +++ b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Image_res_dialog.ui @@ -6,8 +6,8 @@ 0 0 - 448 - 282 + 493 + 241 @@ -17,7 +17,7 @@ - Image Drawing Precision + Load a 3D image @@ -36,14 +36,17 @@ - Please choose the image type + Please choose the image &type + + + imageType - GroupBox + Drawing settings for a segment image @@ -77,11 +80,14 @@ 1 - Please choose the image drawing precision + Please choose the image drawing &precision true + + precisionList + @@ -97,7 +103,7 @@ 20 - 40 + 0 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 45f03ec43d4..12a32835ba4 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_3_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_3_plugin.cpp @@ -188,7 +188,7 @@ void Mesh_3_plugin::mesh_3(const bool surface_only) else if (NULL != image_item) { item = image_item; - features_protection_available = (polylines_item != NULL); + features_protection_available = (polylines_item != NULL) || !image_item->isGray(); bool fit_wrdtp = true; std::size_t img_wdim = image_item->image()->image()->wdim; @@ -298,6 +298,7 @@ void Mesh_3_plugin::mesh_3(const bool surface_only) ui.protect->setEnabled(features_protection_available); ui.protect->setChecked(features_protection_available); + ui.initializationGroup->setVisible(image_item != NULL && !image_item->isGray()); ui.grayImgGroup->setVisible(image_item != NULL && image_item->isGray()); if (poly_item != NULL) ui.volumeGroup->setVisible(!surface_only && poly_item->polyhedron()->is_closed()); @@ -321,6 +322,7 @@ void Mesh_3_plugin::mesh_3(const bool surface_only) const double tet_sizing = !ui.noTetSizing->isChecked() ? 0 : ui.tetSizing->value(); const double edge_size = !ui.noEdgeSizing->isChecked() ? DBL_MAX : ui.edgeSizing->value(); const bool protect_features = ui.protect->isChecked(); + const bool detect_connected_components = ui.detectComponents->isChecked(); const int manifold = ui.manifoldCheckBox->isChecked() ? 1 : 0; const float iso_value = ui.iso_value_spinBox->value(); const float value_outside = ui.value_outside_spinBox->value(); @@ -402,6 +404,7 @@ void Mesh_3_plugin::mesh_3(const bool surface_only) protect_features, manifold, scene, + detect_connected_components, image_item->isGray(), iso_value, value_outside, 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 9691ec013a8..66257245bf9 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 @@ -81,7 +81,8 @@ Meshing_thread* cgal_code_mesh_3(const Polyhedron* pMesh, param.manifold = manifold; param.protect_features = protect_features; - typedef ::Mesh_function Mesh_function; + typedef ::Mesh_function Mesh_function; Mesh_function* p_mesh_function = new Mesh_function(p_new_item->c3t3(), p_domain, param); return new Meshing_thread(p_mesh_function, p_new_item); @@ -124,7 +125,8 @@ Meshing_thread* cgal_code_mesh_3(const Implicit_function_interface* pfunction, param.edge_sizing = edge_size; param.manifold = manifold; - typedef ::Mesh_function Mesh_function; + typedef ::Mesh_function Mesh_function; Mesh_function* p_mesh_function = new Mesh_function(p_new_item->c3t3(), p_domain, param); return new Meshing_thread(p_mesh_function, p_new_item); @@ -150,6 +152,7 @@ Meshing_thread* cgal_code_mesh_3(const Image* pImage, bool protect_features, const int manifold, CGAL::Three::Scene_interface* scene, + bool detect_connected_components, bool is_gray, float iso_value, float value_outside, @@ -161,6 +164,7 @@ Meshing_thread* cgal_code_mesh_3(const Image* pImage, Mesh_parameters param; param.protect_features = protect_features; + param.detect_connected_components = detect_connected_components; param.facet_angle = facet_angle; param.facet_sizing = facet_sizing; param.facet_approx = facet_approx; @@ -168,6 +172,7 @@ Meshing_thread* cgal_code_mesh_3(const Image* pImage, param.edge_sizing = edge_size; param.tet_shape = tet_shape; param.manifold = manifold; + param.image_3_ptr = pImage; CGAL::Timer timer; Scene_c3t3_item* p_new_item = new Scene_c3t3_item; p_new_item->setScene(scene); @@ -187,7 +192,8 @@ Meshing_thread* cgal_code_mesh_3(const Image* pImage, p_domain->add_features(polylines.begin(), polylines.end()); } timer.start(); - typedef ::Mesh_function Mesh_function; + typedef ::Mesh_function Mesh_function; Mesh_function* p_mesh_function = new Mesh_function(p_new_item->c3t3(), p_domain, param); return new Meshing_thread(p_mesh_function, p_new_item); @@ -214,7 +220,8 @@ Meshing_thread* cgal_code_mesh_3(const Image* pImage, p_domain->add_features(polylines.begin(), polylines.end()); } timer.start(); - typedef ::Mesh_function Mesh_function; + typedef ::Mesh_function Mesh_function; Mesh_function* p_mesh_function = new Mesh_function(p_new_item->c3t3(), p_domain, param); return new Meshing_thread(p_mesh_function, p_new_item); 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 a18743d2902..59faf1b7bfd 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 @@ -56,6 +56,7 @@ Meshing_thread* cgal_code_mesh_3(const CGAL::Image_3* pImage, bool protect_features, const int manifold, CGAL::Three::Scene_interface* scene, + bool detect_connected_components, bool is_gray = false, float iso_value = 3.f, float value_outside = 0.f, diff --git a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_function.h b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_function.h index a6b4a6b3660..25537a96741 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_function.h +++ b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_function.h @@ -35,12 +35,17 @@ #include #include #include +#include #include "C3t3_type.h" #include "Meshing_thread.h" #include // for C3t3_initializer #include +namespace CGAL { + class Image_3; +} + struct Mesh_parameters { double facet_angle; @@ -51,13 +56,15 @@ struct Mesh_parameters double tet_sizing; double edge_sizing; bool protect_features; + bool detect_connected_components; int manifold; + const CGAL::Image_3* image_3_ptr; inline QStringList log() const; }; -template < typename Domain_ > +template < typename Domain_, typename Image_tag > class Mesh_function : public Mesh_function_interface { @@ -93,6 +100,9 @@ private: typedef CGAL::Mesh_3::Mesher_3 Mesher; + void initialize(const Mesh_criteria& criteria, CGAL::Tag_true); + void initialize(const Mesh_criteria& criteria, CGAL::Tag_false); + private: C3t3& c3t3_; Domain* domain_; @@ -121,6 +131,8 @@ log() const << QString("facet approx error: %1").arg(facet_approx) << QString("tet shape (radius-edge): %1").arg(tet_shape) << QString("tet max size: %1").arg(tet_sizing) + << QString("detect connected components: %1") + .arg(detect_connected_components) << QString("protect features: %1").arg(protect_features); } @@ -128,8 +140,8 @@ log() const // ----------------------------------- // Class Mesh_function // ----------------------------------- -template < typename D_ > -Mesh_function:: +template < typename D_, typename Tag > +Mesh_function:: Mesh_function(C3t3& c3t3, Domain* domain, const Mesh_parameters& p) : c3t3_(c3t3) , domain_(domain) @@ -146,8 +158,8 @@ Mesh_function(C3t3& c3t3, Domain* domain, const Mesh_parameters& p) } -template < typename D_ > -Mesh_function:: +template < typename D_, typename Tag > +Mesh_function:: ~Mesh_function() { delete domain_; @@ -162,9 +174,47 @@ CGAL::Mesh_facet_topology topology(int manifold) { CGAL::FACET_VERTICES_ON_SAME_SURFACE_PATCH); } -template < typename D_ > +template < typename D_, typename Tag > void -Mesh_function:: +Mesh_function:: +initialize(const Mesh_criteria& criteria, CGAL::Tag_true) // for an image +{ + if(p_.detect_connected_components) { + initialize_triangulation_from_labeled_image(c3t3_ + , *domain_ + , *p_.image_3_ptr + , criteria + , typename D_::Image_word_type() + , p_.protect_features); + } else { + initialize(criteria, CGAL::Tag_false()); + } +} + +template < typename D_, typename Tag > +void +Mesh_function:: +initialize(const Mesh_criteria& criteria, CGAL::Tag_false) // for the other domain types +{ + // Initialization of the mesh, either with the protection of sharp + // features, or with the initial points (or both). + // If `detect_connected_components==true`, the initialization is + // already done. + CGAL::internal::Mesh_3::C3t3_initializer< + C3t3, + Domain, + Mesh_criteria, + CGAL::internal::Mesh_3::has_Has_features::value >() + (c3t3_, + *domain_, + criteria, + p_.protect_features); +} + + +template < typename D_, typename Tag > +void +Mesh_function:: launch() { #ifdef CGAL_MESH_3_INITIAL_POINTS_NO_RANDOM_SHOOTING @@ -180,17 +230,7 @@ launch() Cell_criteria(p_.tet_shape, p_.tet_sizing)); - // Initialization of the mesh, either with the protection of sharp - // features, or with the initial points (or both). - CGAL::internal::Mesh_3::C3t3_initializer< - C3t3, - Domain, - Mesh_criteria, - CGAL::internal::Mesh_3::has_Has_features::value >() - (c3t3_, - *domain_, - criteria, - p_.protect_features); + initialize(criteria, CGAL::Boolean_tag()); // Build mesher and launch refinement process mesher_ = new Mesher(c3t3_, *domain_, criteria); @@ -219,27 +259,27 @@ launch() } -template < typename D_ > +template < typename D_, typename Tag > void -Mesh_function:: +Mesh_function:: stop() { continue_ = false; } -template < typename D_ > +template < typename D_, typename Tag > QStringList -Mesh_function:: +Mesh_function:: parameters_log() const { return p_.log(); } -template < typename D_ > +template < typename D_, typename Tag > QString -Mesh_function:: +Mesh_function:: status(double time_period) const { QString result; diff --git a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Meshing_dialog.ui b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Meshing_dialog.ui index d70ca06d79f..e9bc4d64665 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Meshing_dialog.ui +++ b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Meshing_dialog.ui @@ -10,7 +10,7 @@ 0 0 473 - 637 + 662 @@ -47,6 +47,38 @@ + + + + Mesh initialization + + + + + + + + + + + + + Qt::LeftToRight + + + &Detect all connected components + + + Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter + + + detectComponents + + + + + + @@ -88,9 +120,6 @@ - - Qt::RightToLeft - @@ -102,7 +131,10 @@ - Protect sharp edges + &Protect sharp edges + + + protect @@ -507,8 +539,8 @@ accept() - 397 - 544 + 403 + 655 157 @@ -523,8 +555,8 @@ reject() - 397 - 544 + 403 + 655 286 @@ -539,12 +571,12 @@ setVisible(bool) - 69 - 283 + 83 + 339 - 72 - 297 + 78 + 199 diff --git a/Stream_support/doc/Stream_support/CGAL/IO/io.h b/Stream_support/doc/Stream_support/CGAL/IO/io.h index d4de8ed80db..06b6ae5c1a0 100644 --- a/Stream_support/doc/Stream_support/CGAL/IO/io.h +++ b/Stream_support/doc/Stream_support/CGAL/IO/io.h @@ -148,14 +148,21 @@ Specializations of `Output_rep` should provide the following features: template< class F > struct Output_rep< Some_type, F > { -Output_rep( const Some_type& t ); -std::ostream& operator()( std::ostream& out ) const; + static const bool is_specialized = true; + Output_rep( const Some_type& t ); + std::ostream& operator()( std::ostream& out ) const; }; \endcode You can also specialize for a formatting tag `F`. +The constant `is_specialized` can be tested by meta-programming tools to +verify that a given type can be used with `oformat()`. Its value has to be +`true` in a specialization of `Output_rep`. When there is no specialization +for a type, the class template `Output_rep` defines `is_specialized` to the +default value `false`. + */ template< typename T, typename F > class Output_rep { diff --git a/Stream_support/include/CGAL/IO/io.h b/Stream_support/include/CGAL/IO/io.h index 9f7be26865b..e270de17db8 100644 --- a/Stream_support/include/CGAL/IO/io.h +++ b/Stream_support/include/CGAL/IO/io.h @@ -56,8 +56,27 @@ public: enum Mode {ASCII = 0, PRETTY, BINARY}; }; +template +struct IO_rep_is_specialized_aux +{ + static const bool is_specialized = true; +}; +template< class Dummy > +const bool IO_rep_is_specialized_aux::is_specialized; + +template +struct IO_rep_is_not_specialized_aux +{ + static const bool is_specialized = false; +}; +template< class Dummy > +const bool IO_rep_is_not_specialized_aux::is_specialized; + +typedef IO_rep_is_specialized_aux IO_rep_is_specialized; +typedef IO_rep_is_not_specialized_aux IO_rep_is_not_specialized; + template -class Output_rep { +class Output_rep : public IO_rep_is_not_specialized { const T& t; public: //! initialize with a const reference to \a t. @@ -95,7 +114,7 @@ oformat( const T& t, F) { return Output_rep(t); } * for external types not supporting our stream IO format. */ template -class Input_rep { +class Input_rep : public IO_rep_is_not_specialized { T& t; public: //! initialize with a reference to \a t. @@ -107,7 +126,7 @@ public: #if CGAL_FORCE_IFORMAT_DOUBLE || \ ( ( _MSC_VER > 1600 ) && (! defined( CGAL_NO_IFORMAT_DOUBLE )) ) template <> -class Input_rep { +class Input_rep : public IO_rep_is_specialized { double& t; public: //! initialize with a reference to \a t. diff --git a/Stream_support/test/Stream_support/test_ioformat.cpp b/Stream_support/test/Stream_support/test_ioformat.cpp index ec669bc1844..e1a45fd6d5b 100644 --- a/Stream_support/test/Stream_support/test_ioformat.cpp +++ b/Stream_support/test/Stream_support/test_ioformat.cpp @@ -59,6 +59,8 @@ void test_io(const NT& x){ } int main() { + CGAL_static_assertion(CGAL::Output_rep::is_specialized == false); + CGAL_static_assertion(CGAL::Input_rep::is_specialized == false); std::cout << "test_io: short "<< std::endl; test_io(12);