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..920e8d3eebe --- /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 3c977493db1..70dd6bfdf1f 100644 --- a/Polyhedron/demo/Polyhedron/C3t3_type.h +++ b/Polyhedron/demo/Polyhedron/C3t3_type.h @@ -9,6 +9,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 @@ -50,8 +51,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..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 @@ -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,29 @@ 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 + + + + typedef Facet_extra_criterion + Extra_criterion; + + + // facet_criteria.add(new Extra_criterion(*p_domain)); + Mesh_criteria criteria(edge_criteria, facet_criteria, cell_criteria); + CGAL::Timer timer; timer.start(); Scene_c3t3_item* p_new_item = new Scene_c3t3_item(CGAL::make_mesh_3(*p_domain, diff --git a/Polyhedron/demo/Polyhedron/Scene_c3t3_item.cpp b/Polyhedron/demo/Polyhedron/Scene_c3t3_item.cpp index 57d25896936..b2921353fd0 100644 --- a/Polyhedron/demo/Polyhedron/Scene_c3t3_item.cpp +++ b/Polyhedron/demo/Polyhedron/Scene_c3t3_item.cpp @@ -691,8 +691,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()); @@ -1088,7 +1086,14 @@ void Scene_c3t3_item::compute_elements() //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(); @@ -1099,7 +1104,6 @@ void Scene_c3t3_item::compute_elements() 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(); - if(cell->subdomain_index() == 0) { QColor color = d->colors[cell->neighbor(index)->subdomain_index()];