From 18fc6a8ec1ebc6a98ae02abcfae1e44293c41fb4 Mon Sep 17 00:00:00 2001 From: Arttu Miettinen Date: Fri, 3 Jan 2020 10:44:49 +0100 Subject: [PATCH] Fixes (cgalmesh part of) issue #7 by replacing mesh initialization in cgalmesh with the custom initialization by Laurent Irineau. --- cgalv2m.m | 10 +- ...tialise_triangulation_from_labeled_image.h | 219 ++++++++++++++++ tools/cgalmesh/mesh_3D_image.cpp | 234 +++++++++--------- 3 files changed, 339 insertions(+), 124 deletions(-) create mode 100644 tools/cgalmesh/initialise_triangulation_from_labeled_image.h diff --git a/cgalv2m.m b/cgalv2m.m index 3a37499..428571e 100644 --- a/cgalv2m.m +++ b/cgalv2m.m @@ -10,7 +10,7 @@ % author: Qianqian Fang (q.fang at neu.edu) % % input: -% vol: a volumetric binary image +% vol: a volumetric binary or labeled image % ix,iy,iz: subvolume selection indices in x,y,z directions % opt: parameters for CGAL mesher, if opt is a structure, then % opt.radbound: defines the maximum surface element size @@ -37,10 +37,10 @@ fprintf(1,'creating surface and tetrahedral mesh from a multi-domain volume ...\n'); -dtype=class(vol); -if(~(islogical(vol) || strcmp(dtype,'uint8'))) - error('cgalmesher can only handle uint8 volumes, you have to convert your image to unit8 first.'); -end +%dtype=class(vol); +%if(~(islogical(vol) || strcmp(dtype,'uint8'))) +% error('cgalmesher can only handle uint8 volumes, you have to convert your image to unit8 first.'); +%end if(~any(vol)) error('no labeled regions found in the input volume.'); diff --git a/tools/cgalmesh/initialise_triangulation_from_labeled_image.h b/tools/cgalmesh/initialise_triangulation_from_labeled_image.h new file mode 100644 index 0000000..2d9b65b --- /dev/null +++ b/tools/cgalmesh/initialise_triangulation_from_labeled_image.h @@ -0,0 +1,219 @@ +// Copyright (c) 2015,2016 GeometryFactory +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// +// $URL: https://github.com/CGAL/cgal/blob/releases/CGAL-5.0-I-190/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h $ +// $Id: initialize_triangulation_from_labeled_image.h 254d60f 2019-10-19T15:23:19+02:00 Sébastien Loriot +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// +// 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 +#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::Mesh_3::internal::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::Geom_traits Gt; + typedef typename Gt::FT FT; + typedef typename Tr::Weighted_point Weighted_point; + typedef typename Tr::Bare_point Bare_point; + typedef typename Tr::Segment Segment_3; + typedef typename Tr::Vertex_handle Vertex_handle; + typedef typename Tr::Cell_handle Cell_handle; + typedef typename Gt::Vector_3 Vector_3; + typedef MeshDomain Mesh_domain; + Tr& tr = c3t3.triangulation(); + typename Gt::Compare_weighted_squared_radius_3 cwsr = + tr.geom_traits().compare_weighted_squared_radius_3_object(); + typename Gt::Construct_point_3 cp = + tr.geom_traits().construct_point_3_object(); + typename Gt::Construct_weighted_point_3 cwp = + tr.geom_traits().construct_weighted_point_3_object(); + if (protect_features) { + init_tr_from_labeled_image_call_init_features + (c3t3, domain, criteria, + CGAL::Mesh_3::internal::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); + } + } + for (const Vector_3& v : directions) + { + const Bare_point test = it->first + v; + const typename Mesh_domain::Intersection intersect = + construct_intersection(Segment_3(it->first, test)); + if (std::get<2>(intersect) != 0) + { + const Bare_point& bpi = std::get<0>(intersect); + Weighted_point pi = cwp(bpi); + // This would cause trouble to optimizers + // check pi will not be hidden + typename Tr::Locate_type lt; + 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) + { + const Weighted_point& wp = tr.point(vit); + if (cwsr(wp, FT(0)) == CGAL::SMALLER) // 0 < wp's weight + conflict_vertices.push_back(vit); + } + } + bool pi_inside_protecting_sphere = false; + for (Vertex_handle cv : conflict_vertices) + { + if (tr.is_infinite(cv)) + continue; + const Weighted_point& cv_wp = tr.point(cv); + if (cwsr(cv_wp, FT(0)) == CGAL::EQUAL) // 0 == wp's weight + continue; + // if the (squared) distance between bpi and cv is smaller or equal than cv's weight + if (cwsr(cv_wp, -tr.min_squared_distance(bpi, cp(cv_wp))) != CGAL::LARGER) + { + pi_inside_protecting_sphere = true; + break; + } + } + if (pi_inside_protecting_sphere) + continue; + const typename Mesh_domain::Index index = std::get<1>(intersect); + 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); + } + // 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::Mesh_3::internal::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 \ No newline at end of file diff --git a/tools/cgalmesh/mesh_3D_image.cpp b/tools/cgalmesh/mesh_3D_image.cpp index d2e7349..9a8e6e4 100644 --- a/tools/cgalmesh/mesh_3D_image.cpp +++ b/tools/cgalmesh/mesh_3D_image.cpp @@ -8,141 +8,137 @@ #include #include -// start PV // #include #include #include -// end PV // -// Domain -struct K: public CGAL::Exact_predicates_inexact_constructions_kernel {}; +#include "initialise_triangulation_from_labeled_image.h" + +// Image typedef CGAL::Image_3 Image; -typedef CGAL::Labeled_image_mesh_domain_3 Mesh_domain; + +// Domain +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Labeled_mesh_domain_3 Mesh_domain; // Triangulation -typedef CGAL::Mesh_triangulation_3::type Tr; +typedef CGAL::Parallel_tag Concurrency_tag; +typedef CGAL::Mesh_triangulation_3::type Tr; typedef CGAL::Mesh_complex_3_in_triangulation_3 C3t3; // Mesh Criteria typedef CGAL::Mesh_criteria_3 Mesh_criteria; -typedef Mesh_criteria::Facet_criteria Facet_criteria; -typedef Mesh_criteria::Cell_criteria Cell_criteria; +typedef Mesh_criteria::Facet_criteria Facet_criteria; +typedef Mesh_criteria::Cell_criteria Cell_criteria; -// start PV // -typedef CGAL::Mesh_constant_domain_field_3 Sizing_field_cell; -// end PV // +typedef CGAL::Mesh_constant_domain_field_3 Sizing_field_cell; -void usage(char *exename){ - printf("usage:\n\t%s input.inr output.mesh \ +void usage(char *exename) { + printf( +"usage:\n\t%s input.inr output.mesh \ \n\ -example:\n\t%s input.inr output.mesh 30 6 4 3 8 123456789\n",exename,exename); +example:\n\t%s input.inr output.mesh 30 6 4 3 8 123456789\n", exename, exename); exit(1); } -int main(int argc,char *argv[]) + +int main(int argc, char *argv[]) { - // Loads image - -// start PV // - //float angle=30.f,ssize=6.f,approx=4.f,reratio=3.f,vsize=8.f; - float angle=30.f,ssize=6.f,approx=4.f,reratio=3.f,maxvol=0.f; -// end PV // - int labelid=0, lid; - - printf("Volume/Surface Mesh Generation Utility (Based on CGAL %s)\n\ -(modified for iso2mesh by Qianqian Fang and Peter Varga)\nhttp://iso2mesh.sf.net\n\n",CGAL_VERSION_STR); - - if(argc!=3&&argc!=8&&argc!=9){ - usage(argv[0]); - } - if(argc>=8){ - angle=atof(argv[3]); - ssize=atof(argv[4]); - approx=atof(argv[5]); - reratio=atof(argv[6]); -// start PV // - //vsize=atof(argv[7]); -// end PV // - } - if(argc==9 && atoi(argv[8])>0){ - printf("RNG seed=%d\n",atoi(argv[8])); - CGAL::Random rd(atoi(argv[8])); - CGAL::Random::State st; - rd.save_state(st); - CGAL::default_random.restore_state(st); - } - Image image; - image.read(argv[1]); - Mesh_domain domain(image); - - // Mesh criteria - Facet_criteria facet_criteria(angle, ssize, approx); // angle, size, approximation - -// start PV // - //Cell_criteria cell_criteria(reratio, vsize); // radius-edge ratio, size - - std::string vsize; - vsize=argv[7]; - - std::vector vsize_vect; - std::vector labels_vect; - - std::stringstream stream_vsize(vsize); - std::string word_vsize; - - while( std::getline(stream_vsize, word_vsize, ':') ){ - int len=sscanf(word_vsize.c_str(), "%d=%f", &lid, &maxvol); - if(len==2){ - labelid=lid; - if(maxvol<=0.f){ - std::cerr << "cell volume must be positive" << std::endl; - exit(-2); - } - labels_vect.push_back(labelid); - vsize_vect.push_back(maxvol); - }else{ - len=sscanf(word_vsize.c_str(), "%f", &maxvol); - if(len!=1 || maxvol<0 || word_vsize.find_first_of('=')!=std::string::npos){ - std::cerr << "invalid sizing field label '" << word_vsize << "', please check your command" << std::endl; - exit(-1); - } - labels_vect.push_back(++labelid); - vsize_vect.push_back(maxvol); - } - } - - float max_vsize = *max_element(vsize_vect.begin(),vsize_vect.end()); - - Sizing_field_cell vsize_cell(max_vsize); - int volume_dimension = 3; - - for (int i=0; i(domain, criteria); - - // Output - std::ofstream medit_file; - if(argc>=8) - medit_file.open(argv[2]); - else - medit_file.open("output.mesh"); - c3t3.output_to_medit(medit_file); - medit_file.close(); - - return 0; + //float angle=30.f,ssize=6.f,approx=4.f,reratio=3.f,vsize=8.f; + float angle = 30.f, ssize = 6.f, approx = 4.f, reratio = 3.f, maxvol = 0.f; + int labelid = 0, lid; + + printf("Volume/Surface Mesh Generation Utility (Based on CGAL %s)\n\ +(modified for iso2mesh by Qianqian Fang and Peter Varga)\nhttp://iso2mesh.sf.net\n\n", CGAL_VERSION_STR); + + if (argc != 3 && argc != 8 && argc != 9) { + usage(argv[0]); + } + if (argc >= 8) { + angle = (float)atof(argv[3]); + ssize = (float)atof(argv[4]); + approx = (float)atof(argv[5]); + reratio = (float)atof(argv[6]); + } + if (argc == 9 && atoi(argv[8]) > 0) { + printf("RNG seed=%d\n", atoi(argv[8])); + CGAL::Random rd(atoi(argv[8])); + CGAL::Random::State st; + rd.save_state(st); + CGAL::default_random.restore_state(st); + } + + Image image; + image.read(argv[1]); + Mesh_domain domain = Mesh_domain::create_labeled_image_mesh_domain(image); + + // Mesh criteria + Facet_criteria facet_criteria(angle, ssize, approx); // angle, size, approximation + + std::string vsize; + vsize = argv[7]; + + std::vector vsize_vect; + std::vector labels_vect; + + std::stringstream stream_vsize(vsize); + std::string word_vsize; + + while (std::getline(stream_vsize, word_vsize, ':')) { + int len = sscanf(word_vsize.c_str(), "%d=%f", &lid, &maxvol); + if (len == 2) { + labelid = lid; + if (maxvol <= 0.f) { + std::cerr << "cell volume must be positive" << std::endl; + exit(-2); + } + labels_vect.push_back(labelid); + vsize_vect.push_back(maxvol); + } + else { + len = sscanf(word_vsize.c_str(), "%f", &maxvol); + if (len != 1 || maxvol < 0 || word_vsize.find_first_of('=') != std::string::npos) { + std::cerr << "invalid sizing field label '" << word_vsize << "', please check your command" << std::endl; + exit(-1); + } + labels_vect.push_back(++labelid); + vsize_vect.push_back(maxvol); + } + } + + float max_vsize = (float)*max_element(vsize_vect.begin(), vsize_vect.end()); + + Sizing_field_cell vsize_cell(max_vsize); + int volume_dimension = 3; + + for (int i = 0; i < vsize_vect.size(); i++) + vsize_cell.set_size(vsize_vect[i], volume_dimension, + domain.index_from_subdomain_index(labels_vect[i])); + + std::cout << "Mesh sizes are (label=size) "; + for (int i = 0; i < vsize_vect.size(); i++) { + std::cout << "(" << labels_vect[i] << "=" << vsize_vect[i] << ") "; + } + std::cout << std::endl; + + Cell_criteria cell_criteria(reratio, vsize_cell); + + Mesh_criteria criteria(facet_criteria, cell_criteria); + + // Meshing + C3t3 c3t3; + CGAL_IMAGE_IO_CASE(image.image(), initialize_triangulation_from_labeled_image(c3t3, domain, image, criteria, (Word)0)) + + CGAL::refine_mesh_3(c3t3, domain, criteria); + + // Output + std::ofstream medit_file; + if (argc >= 8) + medit_file.open(argv[2]); + else + medit_file.open("output.mesh"); + c3t3.output_to_medit(medit_file); + medit_file.close(); + + return 0; }