From 18bbaede12e9831192b167cdd95d5c895a01f8e7 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Wed, 9 Dec 2015 06:30:15 +0100 Subject: [PATCH 1/6] Integrate code for respecting polylines on an Image_3 --- .../internal/Mesh_3/Graph_manipulations.h | 98 +++++ .../CGAL/internal/Mesh_3/split_in_polylines.h | 212 ++++++++++ Mesh_3/include/CGAL/polylines_to_protect.h | 369 ++++++++++++++++++ Polyhedron/demo/Polyhedron/C3t3_type.h | 6 +- Polyhedron/demo/Polyhedron/CMakeLists.txt | 1 - .../Mesh_3_plugin/Facet_extra_criterion.h | 124 ++++++ .../Plugins/Mesh_3_plugin/Mesh_3_plugin.cpp | 10 +- .../Mesh_3_plugin/Mesh_3_plugin_cgal_code.cpp | 22 +- 8 files changed, 835 insertions(+), 7 deletions(-) create mode 100644 Mesh_3/include/CGAL/internal/Mesh_3/Graph_manipulations.h create mode 100644 Mesh_3/include/CGAL/internal/Mesh_3/split_in_polylines.h create mode 100644 Mesh_3/include/CGAL/polylines_to_protect.h create mode 100644 Polyhedron/demo/Polyhedron/Plugins/Mesh_3_plugin/Facet_extra_criterion.h diff --git a/Mesh_3/include/CGAL/internal/Mesh_3/Graph_manipulations.h b/Mesh_3/include/CGAL/internal/Mesh_3/Graph_manipulations.h new file mode 100644 index 00000000000..4e1b31788d2 --- /dev/null +++ b/Mesh_3/include/CGAL/internal/Mesh_3/Graph_manipulations.h @@ -0,0 +1,98 @@ +// Copyright (c) 2015 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_GRAPH_MANIPULATIONS +#define CGAL_INTERNAL_MESH_3_INTERNAL_GRAPH_MANIPULATIONS + +#include +// Assumes the point is a CGAL point. + +#include + +#include + +namespace CGAL { +namespace internal { +namespace Mesh_3 { + +template +struct Graph_manipulations +{ + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::edge_descriptor edge_descriptor; + + std::map p2v; + Graph& g; + + Graph_manipulations(Graph& g) : g(g) {} + + vertex_descriptor get_vertex(const Point_3& p) { + typename std::map::iterator + it = p2v.find(p); + if(it == p2v.end()){ + vertex_descriptor v0 = add_vertex(g); + p2v[p] = v0; + g[v0] = p; + return v0; + } else { + return it->second; + } + } + + vertex_descriptor split(const Point_3& a, const Point_3& b) { + const Point_3 mid = a < b ? midpoint(a, b) : midpoint(b, a); + vertex_descriptor vmid = get_vertex(mid); + typename std::map::iterator + it_a = p2v.find(a), + it_b = p2v.find(b); + if(it_a != p2v.end() && it_b != p2v.end()) { + vertex_descriptor va = it_a->second; + vertex_descriptor vb = it_b->second; + edge_descriptor edge; + bool b; + // test if the edge is already here, using add_edge + boost::tie(edge, b) = add_edge(va, vb, g); + remove_edge(edge, g); + if(!b) { + // The edge was already here. + try_add_edge(va, vmid); + try_add_edge(vb, vmid); + return vmid; + } + } + return vmid; + } + + bool try_add_edge(vertex_descriptor v1, vertex_descriptor v2) { + if(v1 != v2) { + edge_descriptor edge; + bool b; + boost::tie(edge, b) = add_edge(v1, v2, g); + return b; + } else + return false; + } +}; // struct template Graph_manipulations + +} // namespace Mesh_3 +} // namespace internal +} // namespace CGAL + +#endif CGAL_INTERNAL_MESH_3_INTERNAL_GRAPH_MANIPULATIONS diff --git a/Mesh_3/include/CGAL/internal/Mesh_3/split_in_polylines.h b/Mesh_3/include/CGAL/internal/Mesh_3/split_in_polylines.h new file mode 100644 index 00000000000..d7e88a1b421 --- /dev/null +++ b/Mesh_3/include/CGAL/internal/Mesh_3/split_in_polylines.h @@ -0,0 +1,212 @@ +// Copyright (c) 2012-2015 GeometryFactory Sarl (France) +// All rights reserved. +// 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) : Andreas Fabri, Laurent Rineau + +#ifndef CGAL_SPLIT_IN_POLYLINES_H +#define CGAL_SPLIT_IN_POLYLINES_H + +#include +#include +#include +#include +#include +#include + +namespace CGAL { +namespace internal { +namespace Mesh_3 { + +template +void dump_graph_edges(std::ostream& out, const Graph& g) +{ + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::edge_descriptor edge_descriptor; + + BOOST_FOREACH(edge_descriptor e, edges(g)) + { + vertex_descriptor s = source(e, g); + vertex_descriptor t = target(e, g); + out.precision(17); + out << "2 " << g[s] << " " << g[t] << "\n"; + } +} + +template +void dump_graph_edges(const char* filename, const Graph& g) +{ + std::ofstream out(filename); + dump_graph_edges(out, g); +} + +/// Splits a graph at vertices with degree higher than two. +/// The vertices are duplicated, and new incident edges created. +template +void split_in_polylines(Graph& G, Kernel) +{ + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::edge_descriptor edge_descriptor; + typedef typename boost::graph_traits::vertex_iterator vertex_iterator; + typedef typename boost::graph_traits::out_edge_iterator out_edge_iterator; + + vertex_iterator b,e; + boost::tie(b,e) = vertices(G); + std::vector V(b,e); + for(typename std::vector::iterator it = V.begin(); + it != V.end(); + ++it){ + vertex_descriptor v = *it; + bool split = false; + + if(out_degree(v,G) > 2) { + split = true; + } else if(out_degree(v, G) == 2) { + out_edge_iterator out_edge_it, out_edges_end; + boost::tie(out_edge_it, out_edges_end) = out_edges(v, G); + + vertex_descriptor v1 = target(*out_edge_it++, G); + vertex_descriptor v2 = target(*out_edge_it++, G); + CGAL_assertion(out_edge_it == out_edges_end); + + const typename Kernel::Point_3& p = G[v]; + const typename Kernel::Point_3& p1 = G[v1]; + const typename Kernel::Point_3& p2 = G[v2]; + + const typename Kernel::Vector_3 e1 = p1 - p; + const typename Kernel::Vector_3 e2 = p2 - p; + const typename Kernel::FT sc_prod = e1 * e2; + if( sc_prod >= 0 || // angle < 135 degrees (3*pi/4) + (sc_prod < 0 && + CGAL::square(sc_prod) < (e1 * e1) * (e2 * e2) / 2 ) ) + { + split = true; + } + } + + if(split) { + out_edge_iterator b,e; + boost::tie(b,e) = out_edges(v,G); + std::vector E(b,e); + for(unsigned int i = 1; i < E.size(); ++i){ + edge_descriptor e = E[i]; + vertex_descriptor w = target(e,G); + remove_edge(e,G); + vertex_descriptor vc = add_vertex(G); + G[vc] = G[v]; + add_edge(vc,w,G); + } + } + } + CGAL_assertion_code( + BOOST_FOREACH(vertex_descriptor v, vertices(G)) + { + typename boost::graph_traits::degree_size_type + n = out_degree(v, G); + + CGAL_assertion(n == 1 || n == 2); + } + BOOST_FOREACH(edge_descriptor e, edges(G)) + { + vertex_descriptor v = target(e,G); + vertex_descriptor w = source(e, G); + CGAL_assertion(v != w); + CGAL_assertion(G[v] != G[w]); + } + ) +} + + +template +void split_in_polylines(Graph& G, Polylines_container& polylines, Kernel k) +{ + typedef typename Polylines_container::value_type Polyline; + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::vertex_iterator vertex_iterator; + typedef typename boost::graph_traits::out_edge_iterator out_edge_iterator; + + // dump_graph_edges("edges.polylines.txt", G); + split_in_polylines(G, k); + std::set terminal; + vertex_iterator b,e; + for(boost::tie(b,e) = vertices(G); b!=e; ++b){ + if(degree(*b,G) == 1){ + terminal.insert(*b); + } + } + + while(! terminal.empty()){ + typename std::set::iterator it = terminal.begin(); + vertex_descriptor u = *it; + terminal.erase(u); + Polyline V; + polylines.push_back(V); + Polyline& polyline = polylines.back(); + polyline.push_back(G[u]); + + while(out_degree(u,G)!=0){ + CGAL_assertion(out_degree(u,G) == 1); + out_edge_iterator b,e; + boost::tie(b,e) = out_edges(u,G); + vertex_descriptor v = target(*b,G); + CGAL_assertion(G[v] != polyline.back()); + polyline.push_back(G[v]); + remove_edge(b,G); + u = v; + } + terminal.erase(u); + + if(polyline.back() == polyline.front()) + { + CGAL_assertion(polyline.size() > 3); + // Fake cycle. We intended that cycle to be split at polyline.front() + // Split the line in two, arbitrary. + std::size_t n = polyline.size() / 2; + Polyline new_line(polyline.begin() + n, + polyline.end()); + polyline.resize(n+1); + polylines.push_back(new_line); + } + } + // dump_graph_edges("only-cycle-edges.polylines.txt", G); + + std::size_t nb_cycles = 0; + // process cycles + while(num_edges(G) != 0) + { + vertex_descriptor u = source(*edges(G).first, G); + + Polyline V; + polylines.push_back(V); + Polyline& polyline = polylines.back(); + polyline.push_back(G[u]); + + ++nb_cycles; + + CGAL_assertion_code(bool first = true); + while(out_degree(u,G)!=0){ + CGAL_assertion(out_degree(u,G) == 1 || + (first && out_degree(u, G) == 2)); + out_edge_iterator b,e; + boost::tie(b,e) = out_edges(u,G); + vertex_descriptor v = target(*b,G); + polyline.push_back(G[v]); + remove_edge(b,G); + u = v; + CGAL_assertion_code(first = false); + } + } +} + +} // namespace Mesh_3 +} // namespace internal +} // namespace CGAL + +#endif // CGAL_SPLIT_IN_POLYLINES_H diff --git a/Mesh_3/include/CGAL/polylines_to_protect.h b/Mesh_3/include/CGAL/polylines_to_protect.h new file mode 100644 index 00000000000..c5f8717cb35 --- /dev/null +++ b/Mesh_3/include/CGAL/polylines_to_protect.h @@ -0,0 +1,369 @@ +// Copyright (c) 2015 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_POLYLINES_TO_PROTECT_H +#define CGAL_POLYLINES_TO_PROTECT_H + +#include +#include +#include +#include +#include +#include + +namespace CGAL { + + + template + void + polylines_to_protect(const CGAL::Image_3& cgal_image, + const double vx, const double vy, const double vz, + std::vector >& polylines) + { + typedef typename Kernel_traits

::Kernel K; + typedef P Point_3; + typedef boost::adjacency_list Graph; + typedef boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef boost::graph_traits::edge_iterator edge_iterator; + + const int xdim = static_cast(cgal_image.xdim()); + const int ydim = static_cast(cgal_image.ydim()); + const int zdim = static_cast(cgal_image.zdim()); + + const int image_dims[3] = { xdim, ydim, zdim }; + + Graph graph; + internal::Mesh_3::Graph_manipulations g_manip(graph); + + std::size_t + case4 = 0, // 4 colors + case211 = 0, case121 = 0, // 3 colors + case31 = 0, case22 = 0, // 2 colors + case1 = 0; // 1 color + + for(int axis = 0; axis < 3; ++axis) + { + for(int i = 0; i < xdim; i+= (axis == 0 ? (xdim-1) : 1 ) ) + for(int j = 0; j < ydim; j+= (axis == 1 ? (ydim-1) : 1 ) ) + for(int k = 0; k < zdim; k+= (axis == 2 ? (zdim-1) : 1 ) ) + { + + using CGAL::cpp11::array; + using CGAL::cpp11::tuple; + + typedef array Pixel; + + Pixel pix00 = {{i , j , k }}, + pix10 = pix00, pix01 = pix00, pix11 = pix00; + + const int axis_xx = (axis + 1) % 3; + const int axis_yy = (axis + 2) % 3; + + ++pix10[axis_xx]; + ++pix11[axis_xx]; + ++pix01[axis_yy]; + ++pix11[axis_yy]; + if(pix11[0] >= xdim || pix11[1] >= ydim || pix11[2] >= zdim) { + // we have gone too far + continue; + } + typedef unsigned char Image_word_type; + typedef tuple Enriched_pixel; + + array, 2> square = + {{ {{ Enriched_pixel(pix00, Point_3(), Image_word_type()), + Enriched_pixel(pix01, Point_3(), Image_word_type()) }}, + {{ Enriched_pixel(pix10, Point_3(), Image_word_type()), + Enriched_pixel(pix11, Point_3(), Image_word_type()) }} }}; + + std::map pixel_values_set; + for(int ii = 0; ii < 2; ++ii) + for(int jj = 0; jj < 2; ++jj) + { + const Pixel& pixel = get<0>(square[ii][jj]); + double x = pixel[0] * cgal_image.vx(); + double y = pixel[1] * cgal_image.vy(); + double z = pixel[2] * cgal_image.vz(); + get<1>(square[ii][jj]) = Point_3(x, y, z); + get<2>(square[ii][jj]) = cgal_image.value(pixel[0], + pixel[1], + pixel[2]); + ++pixel_values_set[get<2>(square[ii][jj])]; + } + const Point_3& p00 = get<1>(square[0][0]); + const Point_3& p10 = get<1>(square[1][0]); + const Point_3& p01 = get<1>(square[0][1]); + const Point_3& p11 = get<1>(square[1][1]); + + // + // Protect the edges of the cube + // + if(pix00[axis_xx] == 0) { + g_manip.try_add_edge(g_manip.get_vertex(p00), + g_manip.get_vertex(p01)); + } + if(pix11[axis_xx] == image_dims[axis_xx]-1) { + g_manip.try_add_edge(g_manip.get_vertex(p10), + g_manip.get_vertex(p11)); + } + if(pix00[axis_yy] == 0) { + g_manip.try_add_edge(g_manip.get_vertex(p00), + g_manip.get_vertex(p10)); + } + if(pix11[axis_yy] == image_dims[axis_yy]-1) { + g_manip.try_add_edge(g_manip.get_vertex(p01), + g_manip.get_vertex(p11)); + } + + // + // Protect lines inside the square + // + switch(pixel_values_set.size()) { + case 4: { + assert(get<2>(square[0][0]) != get<2>(square[0][1])); + assert(get<2>(square[0][0]) != get<2>(square[1][0])); + assert(get<2>(square[0][0]) != get<2>(square[1][1])); + assert(get<2>(square[1][0]) != get<2>(square[1][1])); + assert(get<2>(square[0][1]) != get<2>(square[1][1])); + assert(get<2>(square[0][1]) != get<2>(square[1][0])); +case_4: + // case 4 or case 2-2 + ++case4; + vertex_descriptor left = g_manip.split(p00, p01); + vertex_descriptor right = g_manip.split(p10, p11); + vertex_descriptor top = g_manip.split(p01, p11); + vertex_descriptor bottom = g_manip.split(p00, p10); + vertex_descriptor vmid = g_manip.get_vertex(midpoint(p00, p11)); + g_manip.try_add_edge(left , vmid); + g_manip.try_add_edge(right , vmid); + g_manip.try_add_edge(top , vmid); + g_manip.try_add_edge(bottom , vmid); + } + break; + case 3: { + if(get<2>(square[0][0]) == get<2>(square[1][1])) { + // Diagonal, but the wrong one. + // Vertical swap + std::swap(square[0][1], square[0][0]); + std::swap(square[1][1], square[1][0]); + } + if(get<2>(square[0][1]) == get<2>(square[1][0])) { + // diagonal case 1-2-1 + assert(get<2>(square[0][1]) == get<2>(square[1][0])); + assert(get<2>(square[1][1]) != get<2>(square[0][0])); + assert(get<2>(square[0][1]) != get<2>(square[0][0])); + assert(get<2>(square[0][1]) != get<2>(square[1][1])); + ++case121; + vertex_descriptor left = g_manip.split(p00, p01); + vertex_descriptor right = g_manip.split(p10, p11); + vertex_descriptor top = g_manip.split(p01, p11); + vertex_descriptor bottom = g_manip.split(p00, p10); + Point_3 inter_left = p00 + (1./3.) * (p11 - p00); + Point_3 inter_right = p00 + (2./3.) * (p11 - p00); + vertex_descriptor v_int_left = g_manip.get_vertex(inter_left); + vertex_descriptor v_int_right = g_manip.get_vertex(inter_right); + g_manip.try_add_edge(left, v_int_left); + g_manip.try_add_edge(v_int_left, bottom); + g_manip.try_add_edge(top, v_int_right); + g_manip.try_add_edge(v_int_right, right); + } else { + // case 2-1-1 + if(get<2>(square[0][0]) == get<2>(square[1][0])) { + // Diagonal swap + std::swap(square[0][1], square[1][0]); + } else + if(get<2>(square[0][1]) == get<2>(square[1][1])) { + // The other diagonal swap + std::swap(square[0][0], square[1][1]); + } else + if(get<2>(square[1][0]) == get<2>(square[1][1])) { + // Vertical swap + std::swap(square[0][0], square[1][0]); + std::swap(square[0][1], square[1][1]); + } + assert(get<2>(square[0][0]) == get<2>(square[0][1])); + assert(get<2>(square[0][0]) != get<2>(square[1][0])); + assert(get<2>(square[0][0]) != get<2>(square[1][1])); + ++case211; + Point_3 midleft = midpoint(p00, p01); + Point_3 midright = midpoint(p10, p11); + Point_3 inter = midleft + (2./3)*(midright-midleft); + vertex_descriptor v_inter = g_manip.get_vertex(inter); + vertex_descriptor right = g_manip.split(p10, p11); + vertex_descriptor top = g_manip.split(p01, p11); + vertex_descriptor bottom = g_manip.split(p00, p10); + g_manip.try_add_edge(top, v_inter); + g_manip.try_add_edge(bottom, v_inter); + g_manip.try_add_edge(right, v_inter); + } + } + break; + case 2: { + if(pixel_values_set.begin()->second == + pixel_values_set.rbegin()->second) + { + // Case of two colors with two pixels each. + + if(get<2>(square[0][0])==get<2>(square[1][0])) { + // case 2-2, diagonal swap + std::swap(square[0][1], square[1][0]); + assert(get<2>(square[0][0])==get<2>(square[0][1])); + } + if(get<2>(square[1][0])==get<2>(square[1][1])) { + // case 2-2, vertical swap + std::swap(square[0][1], square[1][1]); + std::swap(square[0][0], square[1][0]); + assert(get<2>(square[0][0])==get<2>(square[0][1])); + } + if(get<2>(square[0][1])==get<2>(square[1][1])) { + // case 2-2, diagonal swap + std::swap(square[0][0], square[1][1]); + assert(get<2>(square[0][0])==get<2>(square[0][1])); + } + + if(get<2>(square[0][0])==get<2>(square[0][1])) { + // vertical case 2-2 + ++case22; + assert(get<2>(square[1][0])==get<2>(square[1][1])); + assert(get<2>(square[1][0])!=get<2>(square[0][1])); + vertex_descriptor top = g_manip.split(p01, p11); + vertex_descriptor bottom = g_manip.split(p00, p10); + g_manip.try_add_edge(top, bottom); + } else { + // Else diagonal case case 2-2 + // Same as the case with 4 colors + assert(get<2>(square[0][0])==get<2>(square[1][1])); + assert(get<2>(square[1][0])==get<2>(square[0][1])); + assert(get<2>(square[0][0])!=get<2>(square[0][1])); + goto case_4; + } + } + else { + // case of two colors with one pixel green and three red + Image_word_type value_alone; + if(pixel_values_set.begin()->second == 1) { + value_alone = pixel_values_set.begin()->first; + } else { + assert(pixel_values_set.begin()->second == 3); + assert(pixel_values_set.rbegin()->second ==1); + value_alone = pixel_values_set.rbegin()->first; + } + if(get<2>(square[0][1]) == value_alone) { + // central symmetry + std::swap(square[0][1], square[1][0]); + std::swap(square[0][0], square[1][1]); + assert(get<2>(square[1][0]) == value_alone); + } + if(get<2>(square[1][1]) == value_alone) { + // vertical swap + std::swap(square[0][0], square[0][1]); + std::swap(square[1][0], square[1][1]); + assert(get<2>(square[1][0]) == value_alone); + } + if(get<2>(square[0][0]) == value_alone) { + // horizontal swap + std::swap(square[0][1], square[1][1]); + std::swap(square[0][0], square[1][0]); + assert(get<2>(square[1][0]) == value_alone); + } + ++case31; + assert(get<2>(square[1][0]) == value_alone); + assert(get<2>(square[1][0]) != get<2>(square[0][0])); + assert(get<2>(square[1][1]) == get<2>(square[0][0])); + assert(get<2>(square[0][1]) == get<2>(square[0][0])); + vertex_descriptor bottom = g_manip.split(p00, p10); + vertex_descriptor old = bottom; + + Point_3 inter; + + vertex_descriptor v_int; + for(double x = 0.55; x < 1.; x+= 0.05) + { + inter = p00 + + x * (p10 - p00) // x + + (1.-1./(2.*x)) * (p01 - p00); // y + v_int = g_manip.get_vertex(inter); + g_manip.try_add_edge(old, v_int); + old = v_int; + } + + vertex_descriptor right = g_manip.split(p10, p11); + g_manip.try_add_edge(v_int, right); + } + } + break; + default: // case 1 + ++case1; + // nothing to do + break; + } + } + } + // std::cerr << "case 4: " << case4 << std::endl; + // std::cerr << "case 2-1-1: " << case211 << std::endl; + // std::cerr << "case 1-2-1: " << case121 << std::endl; + // std::cerr << "case 3-1: " << case31 << std::endl; + // std::cerr << "case 2-2: " << case22 << std::endl; + // std::cerr << "case 1: " << case1 << std::endl; + + const std::ptrdiff_t nb_facets = + case4 + case211 + case121 + case31 + case22 + case1; + const std::ptrdiff_t expected_nb_facets = + 2*((xdim-1)*(ydim-1) + (ydim-1)*(zdim-1) + (xdim-1)*(zdim-1)); + + // std::cerr << "nb of facets: " << nb_facets << std::endl + // << " expected nb of facets: " << expected_nb_facets << std::endl; + + CGAL_assertion(nb_facets == expected_nb_facets); + CGAL_USE(nb_facets); CGAL_USE(expected_nb_facets); + + + internal::Mesh_3::split_in_polylines(graph, polylines, K()); + { + std::ofstream out("polylines.cgal.txt"); + out.precision(17); + + BOOST_FOREACH(const std::vector& polyline, polylines) { + out << polyline.size(); + BOOST_FOREACH(const Point_3& p, polyline ) { + out << " " << p; + } + out << std::endl; + } + } + } + + + template + void + polylines_to_protect(const CGAL::Image_3& cgal_image, + std::vector >& polylines) + { + polylines_to_protect(cgal_image, + cgal_image.vx(), cgal_image.vy(),cgal_image.vz(), + polylines); + } + + +} // namespace CGAL + + +#endif // CGAL_POLYLINES_TO_PROTECT_H + diff --git a/Polyhedron/demo/Polyhedron/C3t3_type.h b/Polyhedron/demo/Polyhedron/C3t3_type.h index e8be6c02511..beb15c389f0 100644 --- a/Polyhedron/demo/Polyhedron/C3t3_type.h +++ b/Polyhedron/demo/Polyhedron/C3t3_type.h @@ -5,6 +5,7 @@ #ifdef CGAL_MESH_3_DEMO_ACTIVATE_SEGMENTED_IMAGES #include "Image_type.h" #include +#include #include "Polyhedron_demo_mesh_3_labeled_mesh_domain_3.h" #endif #ifdef CGAL_MESH_3_DEMO_ACTIVATE_IMPLICIT_FUNCTIONS @@ -46,8 +47,9 @@ typedef CGAL::Polyhedral_mesh_domain_with_features_3< // The last `Tag_true` says the Patch_id type will be int, and not pair #ifdef CGAL_MESH_3_DEMO_ACTIVATE_SEGMENTED_IMAGES -typedef CGAL::Labeled_image_mesh_domain_3 Image_domain; -typedef CGAL::Polyhedron_demo_labeled_mesh_domain_3 Image_mesh_domain; +typedef CGAL::Labeled_image_mesh_domain_3 Image_domain1; +typedef CGAL::Polyhedron_demo_labeled_mesh_domain_3 Image_domain; +typedef CGAL::Mesh_domain_with_polyline_features_3 Image_mesh_domain; #endif #ifdef CGAL_MESH_3_DEMO_ACTIVATE_IMPLICIT_FUNCTIONS typedef Wrapper Function_wrapper; diff --git a/Polyhedron/demo/Polyhedron/CMakeLists.txt b/Polyhedron/demo/Polyhedron/CMakeLists.txt index 6a963e3b909..59b73060321 100644 --- a/Polyhedron/demo/Polyhedron/CMakeLists.txt +++ b/Polyhedron/demo/Polyhedron/CMakeLists.txt @@ -22,7 +22,6 @@ set(CMAKE_LIBRARY_OUTPUT_DIRECTORY "${CMAKE_RUNTIME_OUTPUT_DIRECTORY}") # Include this package's headers first include_directories( BEFORE ./ ./include ../../include ./CGAL_demo ) - list(INSERT CMAKE_MODULE_PATH 0 "${CMAKE_CURRENT_SOURCE_DIR}") add_subdirectory( implicit_functions ) diff --git a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3_plugin/Facet_extra_criterion.h b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3_plugin/Facet_extra_criterion.h new file mode 100644 index 00000000000..c50321c4bd9 --- /dev/null +++ b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3_plugin/Facet_extra_criterion.h @@ -0,0 +1,124 @@ +// Copyright (c) 2015 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_FACET_EXTRA_CRITERION_H +#define CGAL_FACET_EXTRA_CRITERION_H + +#include + + +template +class Facet_extra_criterion : + public CGAL::Mesh_3::Abstract_criterion +{ +private: + typedef typename Tr::Facet Facet; + + typedef CGAL::Mesh_3::Abstract_criterion Base; + typedef typename Base::Quality Quality; + typedef typename Base::Badness Badness; + + typedef Facet_extra_criterion Self; + + const Domain& domain; + +public: + /// Constructor + Facet_extra_criterion(const Domain& domain) : domain(domain) {} + /// Destructor + ~Facet_extra_criterion() {} + +protected: + virtual void do_accept(Visitor_& v) const + { + v.visit(*this); + } + + virtual Self* do_clone() const + { + // Call copy ctor on this + return new Self(*this); + } + + + virtual Badness do_is_bad (const Facet& f) const + { + typedef typename Tr::Vertex_handle Vertex_handle; + typedef typename Tr::Cell_handle Cell_handle; + typedef typename Domain::Surface_patch_index Surface_patch_index; + + const Cell_handle& ch = f.first; + const int& i = f.second; + + const Vertex_handle& v1 = ch->vertex((i+1)&3); + const Vertex_handle& v2 = ch->vertex((i+2)&3); + const Vertex_handle& v3 = ch->vertex((i+3)&3); + + Surface_patch_index index = Surface_patch_index(); + bool is_index_initialized = false; + + std::cerr << typeid(index).name() << std::endl; + if ( v1->in_dimension() == 2 ) + { + index = domain.surface_patch_index(v1->index()); + if(index != 0) { // (index.first != 0 && index.second != 0) + //index = Surface_patch_index(std::make_pair(1,1)); + } + is_index_initialized = true; + } + + if ( v2->in_dimension() == 2 ) + { + Surface_patch_index index2 = domain.surface_patch_index(v2->index()); + if(index2) { // (index2.first != 0 && index2.second != 0){ + //index2 = Surface_patch_index(1,1); + } + if ( is_index_initialized ) + { + if ( !(index2 == index) ) + { + return Badness(Quality(1)); + } + } + else + { + index = index2; + is_index_initialized = true; + } + } + + if ( v3->in_dimension() == 2 ) + { + Surface_patch_index index3 = domain.surface_patch_index(v3->index()); + if(index3 != 0) { // (index3.first != 0 && index3.second != 0) + // index3 = Surface_patch_index(1,1); + } + if ( is_index_initialized && !(index3 == index) ) + { + return Badness(Quality(1)); + } + } + + return Badness(); + } + +}; // end class Facet_extra_criterion + + +#endif // CGAL_FACET_EXTRA_CRITERION_H diff --git a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3_plugin/Mesh_3_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3_plugin/Mesh_3_plugin.cpp index bbc0fc75a45..72a6cf8de43 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3_plugin/Mesh_3_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3_plugin/Mesh_3_plugin.cpp @@ -55,6 +55,7 @@ Scene_item* cgal_code_mesh_3(const Image* pImage, const double facet_approx, const double tet_sizing, const double tet_shape, + const bool protect_features, CGAL::Three::Scene_interface* scene); #endif @@ -134,7 +135,7 @@ void Polyhedron_demo_mesh_3_plugin::mesh_3() Scene_item* item = NULL; bool features_protection_available = false; - if (NULL != poly_item) + if(NULL != poly_item) { item = poly_item; features_protection_available = true; @@ -143,7 +144,11 @@ void Polyhedron_demo_mesh_3_plugin::mesh_3() else if (NULL != function_item) { item = function_item; } #endif #ifdef CGAL_MESH_3_DEMO_ACTIVATE_SEGMENTED_IMAGES - else if (NULL != image_item) { item = image_item; } + else if (NULL != image_item) + { + item = image_item; + features_protection_available = true; + } #endif if (NULL == item) @@ -284,6 +289,7 @@ void Polyhedron_demo_mesh_3_plugin::mesh_3() approx, tet_sizing, radius_edge, + protect_features, scene); } #endif diff --git a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3_plugin/Mesh_3_plugin_cgal_code.cpp b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3_plugin/Mesh_3_plugin_cgal_code.cpp index 784b0471231..a7afa95c1d6 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3_plugin/Mesh_3_plugin_cgal_code.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3_plugin/Mesh_3_plugin_cgal_code.cpp @@ -2,6 +2,7 @@ #include #include +#include #include #include @@ -11,6 +12,7 @@ #include #include "Mesh_function.h" +#include "Facet_extra_criterion.h" #include using namespace CGAL::Three; @@ -41,7 +43,6 @@ Scene_item* cgal_code_mesh_3(const Polyhedron* pMesh, Facet_criteria facet_criteria(angle, facet_sizing, approx); // angle, size, approximation Cell_criteria cell_criteria(tet_shape, tet_sizing); // radius-edge ratio, size Mesh_criteria criteria(edge_criteria, facet_criteria, cell_criteria); - CGAL::Timer timer; timer.start(); std::cerr << "Meshing file \"" << qPrintable(filename) << "\"\n"; @@ -52,6 +53,7 @@ Scene_item* cgal_code_mesh_3(const Polyhedron* pMesh, std::cerr << "Build AABB tree..."; // Create domain Polyhedral_mesh_domain domain(*pMesh); + if(protect_features) { domain.detect_features(); } @@ -126,7 +128,7 @@ Scene_item* cgal_code_mesh_3(const Implicit_function_interface* pfunction, CGAL::Timer timer; p_new_item->set_scene(scene); - std::cerr << "done (" << timer.time() << " ms, " + std::cerr << "done (" << timer.time() << " sec, " << p_new_item->c3t3().triangulation().number_of_vertices() << " vertices)" << std::endl; @@ -154,6 +156,7 @@ Scene_item* cgal_code_mesh_3(const Image* pImage, const double facet_approx, const double tet_sizing, const double tet_shape, + const bool protect_features, CGAL::Three::Scene_interface* scene) { @@ -162,11 +165,26 @@ Scene_item* cgal_code_mesh_3(const Image* pImage, Image_mesh_domain* p_domain = new Image_mesh_domain(*pImage, 1e-6); + if(protect_features){ + std::vector > polylines; + CGAL::polylines_to_protect(*pImage, polylines); + p_domain->add_features(polylines.begin(), polylines.end()); + } + // Set mesh criteria Edge_criteria edge_criteria(facet_sizing); Facet_criteria facet_criteria(facet_angle, facet_sizing, facet_approx); // angle, size, approximation Cell_criteria cell_criteria(tet_shape, tet_sizing); // radius-edge ratio, size Mesh_criteria criteria(edge_criteria, facet_criteria, cell_criteria); + + + typedef Facet_extra_criterion + Extra_criterion; + + facet_criteria.add(new Extra_criterion(*p_domain)); + CGAL::Timer timer; timer.start(); Scene_c3t3_item* p_new_item = new Scene_c3t3_item(CGAL::make_mesh_3(*p_domain, From 056ea83fa1e25baa7b28be9112f9fd005d91cfe8 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Wed, 9 Dec 2015 10:04:20 +0100 Subject: [PATCH 2/6] move criteria --- .../Plugins/Mesh_3_plugin/Mesh_3_plugin_cgal_code.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3_plugin/Mesh_3_plugin_cgal_code.cpp b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3_plugin/Mesh_3_plugin_cgal_code.cpp index a7afa95c1d6..f5dd6daa68c 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3_plugin/Mesh_3_plugin_cgal_code.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3_plugin/Mesh_3_plugin_cgal_code.cpp @@ -175,7 +175,7 @@ Scene_item* cgal_code_mesh_3(const Image* pImage, Edge_criteria edge_criteria(facet_sizing); Facet_criteria facet_criteria(facet_angle, facet_sizing, facet_approx); // angle, size, approximation Cell_criteria cell_criteria(tet_shape, tet_sizing); // radius-edge ratio, size - Mesh_criteria criteria(edge_criteria, facet_criteria, cell_criteria); + typedef Facet_extra_criterion(*p_domain, From b97ccb832bc5f30baef5bbdf915a9ddc848dba1b Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Wed, 9 Dec 2015 10:28:55 +0100 Subject: [PATCH 3/6] capitalize file name --- Mesh_3/include/CGAL/polylines_to_protect.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Mesh_3/include/CGAL/polylines_to_protect.h b/Mesh_3/include/CGAL/polylines_to_protect.h index c5f8717cb35..920e8d3eebe 100644 --- a/Mesh_3/include/CGAL/polylines_to_protect.h +++ b/Mesh_3/include/CGAL/polylines_to_protect.h @@ -25,7 +25,7 @@ #include #include #include -#include +#include #include namespace CGAL { From c8bd7850618f3c2a921a65b154156311a97b8349 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Wed, 9 Dec 2015 12:35:18 +0100 Subject: [PATCH 4/6] do not use the facet criteria --- .../Plugins/Mesh_3_plugin/Mesh_3_plugin_cgal_code.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3_plugin/Mesh_3_plugin_cgal_code.cpp b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3_plugin/Mesh_3_plugin_cgal_code.cpp index f5dd6daa68c..9e46f202cc2 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3_plugin/Mesh_3_plugin_cgal_code.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3_plugin/Mesh_3_plugin_cgal_code.cpp @@ -183,7 +183,8 @@ Scene_item* cgal_code_mesh_3(const Image* pImage, Mesh_criteria::Facet_criteria::Visitor> Extra_criterion; - facet_criteria.add(new Extra_criterion(*p_domain)); + + // facet_criteria.add(new Extra_criterion(*p_domain)); Mesh_criteria criteria(edge_criteria, facet_criteria, cell_criteria); From fefa699a3d935d0fa2a7d733d640334c1c1ad319 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Wed, 9 Dec 2015 12:35:49 +0100 Subject: [PATCH 5/6] no need to compute a normal --- Polyhedron/demo/Polyhedron/Scene_c3t3_item.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/Polyhedron/demo/Polyhedron/Scene_c3t3_item.cpp b/Polyhedron/demo/Polyhedron/Scene_c3t3_item.cpp index ba901acf1d4..b83ad2e022a 100644 --- a/Polyhedron/demo/Polyhedron/Scene_c3t3_item.cpp +++ b/Polyhedron/demo/Polyhedron/Scene_c3t3_item.cpp @@ -662,8 +662,6 @@ void Scene_c3t3_item::draw_triangle_edges(const Kernel::Point_3& pa, const Kernel::Point_3& pc)const { #undef darker - Kernel::Vector_3 n = cross_product(pb - pa, pc - pa); - n = n / CGAL::sqrt(n*n); positions_lines.push_back(pa.x()); positions_lines.push_back(pa.y()); positions_lines.push_back(pa.z()); From 46eec2d9f054b9fcfeb5abf6c103e2661f0266c3 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Wed, 9 Dec 2015 13:43:11 +0100 Subject: [PATCH 6/6] do not compute sidedness twice --- .../demo/Polyhedron/Scene_c3t3_item.cpp | 37 ++++++++----------- 1 file changed, 15 insertions(+), 22 deletions(-) diff --git a/Polyhedron/demo/Polyhedron/Scene_c3t3_item.cpp b/Polyhedron/demo/Polyhedron/Scene_c3t3_item.cpp index b83ad2e022a..470059d1836 100644 --- a/Polyhedron/demo/Polyhedron/Scene_c3t3_item.cpp +++ b/Polyhedron/demo/Polyhedron/Scene_c3t3_item.cpp @@ -959,16 +959,16 @@ void Scene_c3t3_item::compute_elements() const if (isEmpty()) return; - const Kernel::Plane_3& plane = this->plane(); - GLdouble clip_plane[4]; - clip_plane[0] = -plane.a(); - clip_plane[1] = -plane.b(); - clip_plane[2] = -plane.c(); - clip_plane[3] = -plane.d(); - //The facets - { + { + const Kernel::Plane_3& plane = this->plane(); + GLdouble clip_plane[4]; + clip_plane[0] = plane.a(); + clip_plane[1] = plane.b(); + clip_plane[2] = plane.c(); + clip_plane[3] = plane.d(); + for (C3t3::Facet_iterator fit = c3t3().facets_begin(), end = c3t3().facets_end(); @@ -979,21 +979,14 @@ void Scene_c3t3_item::compute_elements() const const Kernel::Point_3& pa = cell->vertex((index + 1) & 3)->point(); const Kernel::Point_3& pb = cell->vertex((index + 2) & 3)->point(); const Kernel::Point_3& pc = cell->vertex((index + 3) & 3)->point(); - typedef Kernel::Oriented_side Side; - using CGAL::ON_ORIENTED_BOUNDARY; - const Side sa = plane.oriented_side(pa); - const Side sb = plane.oriented_side(pb); - const Side sc = plane.oriented_side(pc); - bool is_shown = false; - if (pa.x() * clip_plane[0] + pa.y() * clip_plane[1] + pa.z() * clip_plane[2] + clip_plane[3] > 0 - && pb.x() * clip_plane[0] + pb.y() * clip_plane[1] + pb.z() * clip_plane[2] + clip_plane[3] > 0 - && pc.x() * clip_plane[0] + pc.y() * clip_plane[1] + pc.z() * clip_plane[2] + clip_plane[3] > 0) - is_shown = true; - if (is_shown && sa != ON_ORIENTED_BOUNDARY && - sb != ON_ORIENTED_BOUNDARY && - sc != ON_ORIENTED_BOUNDARY && - sb == sa && sc == sa) + bool is_shown = + pa.x() * clip_plane[0] + pa.y() * clip_plane[1] + pa.z() * clip_plane[2] + clip_plane[3] < 0 + && pb.x() * clip_plane[0] + pb.y() * clip_plane[1] + pb.z() * clip_plane[2] + clip_plane[3] < 0 + && pc.x() * clip_plane[0] + pc.y() * clip_plane[1] + pc.z() * clip_plane[2] + clip_plane[3] < 0; + + + if (is_shown) { if(cell->subdomain_index() == 0) {