merged branch that deals with 1D features on the sides of a cube

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 <CGAL/Kernel/global_functions_3.h>
+// Assumes the point is a CGAL point.
+
+#include <boost/graph/graph_traits.hpp>
+
+#include <map>
+
+namespace CGAL {
+namespace internal {
+namespace Mesh_3 {
+
+template <typename Graph, typename Point_3>
+struct Graph_manipulations
+{
+  typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex_descriptor;
+  typedef typename boost::graph_traits<Graph>::edge_descriptor edge_descriptor;
+
+  std::map<Point_3, vertex_descriptor> p2v;
+  Graph& g;
+
+  Graph_manipulations(Graph& g) : g(g) {}
+
+  vertex_descriptor get_vertex(const Point_3& p) {
+    typename std::map<Point_3, vertex_descriptor>::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<Point_3, vertex_descriptor>::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

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 <vector>
+#include <set>
+#include <iostream>
+#include <boost/foreach.hpp>
+#include <CGAL/number_utils.h>
+#include <boost/graph/graph_traits.hpp>
+
+namespace CGAL {
+namespace internal {
+namespace Mesh_3 {
+
+template <typename Graph>
+void dump_graph_edges(std::ostream& out, const Graph& g)
+{
+  typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex_descriptor;
+  typedef typename boost::graph_traits<Graph>::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 <typename Graph>
+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 <typename Graph, typename Kernel>
+void split_in_polylines(Graph& G, Kernel)
+{
+  typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex_descriptor;
+  typedef typename boost::graph_traits<Graph>::edge_descriptor edge_descriptor;
+  typedef typename boost::graph_traits<Graph>::vertex_iterator vertex_iterator;
+  typedef typename boost::graph_traits<Graph>::out_edge_iterator out_edge_iterator;
+
+  vertex_iterator b,e;
+  boost::tie(b,e) = vertices(G);
+  std::vector<vertex_descriptor> V(b,e);
+  for(typename std::vector<vertex_descriptor>::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<edge_descriptor> 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<Graph>::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 <typename Graph, typename Polylines_container, typename Kernel>
+void split_in_polylines(Graph& G, Polylines_container& polylines, Kernel k)
+{
+  typedef typename Polylines_container::value_type Polyline;
+  typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex_descriptor;
+  typedef typename boost::graph_traits<Graph>::vertex_iterator vertex_iterator;
+  typedef typename boost::graph_traits<Graph>::out_edge_iterator out_edge_iterator;
+
+  // dump_graph_edges("edges.polylines.txt", G);
+  split_in_polylines(G, k);
+  std::set<vertex_descriptor> 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<vertex_descriptor>::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

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 <vector>
+#include <map>
+#include <CGAL/Image_3.h>
+#include <CGAL/internal/Mesh_3/split_in_polylines.h>
+#include <CGAL/internal/Mesh_3/Graph_manipulations.h>
+#include <boost/graph/adjacency_list.hpp>
+
+namespace CGAL {
+
+
+  template <typename P>
+  void
+  polylines_to_protect(const CGAL::Image_3& cgal_image,
+                       const double vx, const double vy, const double vz,
+                       std::vector<std::vector<P> >& polylines)
+  {
+    typedef typename Kernel_traits<P>::Kernel K;
+    typedef P Point_3;
+    typedef boost::adjacency_list<boost::setS, boost::setS, boost::undirectedS, Point_3> Graph;
+    typedef boost::graph_traits<Graph>::vertex_descriptor vertex_descriptor;
+    typedef boost::graph_traits<Graph>::edge_iterator edge_iterator;
+
+    const int xdim = static_cast<int>(cgal_image.xdim());
+    const int ydim = static_cast<int>(cgal_image.ydim());
+    const int zdim = static_cast<int>(cgal_image.zdim());
+
+  const int image_dims[3] = { xdim, ydim, zdim };
+
+  Graph graph;
+  internal::Mesh_3::Graph_manipulations<Graph, Point_3> 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<int, 3> 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<Pixel, Point_3, Image_word_type> Enriched_pixel;
+
+          array<array<Enriched_pixel, 2>, 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<Image_word_type, int> 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<Point_3>& polyline, polylines) {
+      out << polyline.size();
+      BOOST_FOREACH(const Point_3& p, polyline ) {
+        out << " " << p;
+      }
+      out << std::endl;
+    }
+  }
+  }
+
+
+  template <typename P>
+  void
+  polylines_to_protect(const CGAL::Image_3& cgal_image,
+                       std::vector<std::vector<P> >& polylines)
+  {
+    polylines_to_protect(cgal_image, 
+                         cgal_image.vx(), cgal_image.vy(),cgal_image.vz(),
+                         polylines);
+  }
+
+
+} // namespace CGAL
+
+
+#endif // CGAL_POLYLINES_TO_PROTECT_H
+

Polyhedron/demo/Polyhedron/C3t3_type.h

@@ -9,6 +9,7 @@
 #ifdef CGAL_MESH_3_DEMO_ACTIVATE_SEGMENTED_IMAGES
 #include "Image_type.h"
 #include <CGAL/Labeled_image_mesh_domain_3.h>
+#include <CGAL/Mesh_domain_with_polyline_features_3.h>
 #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<int, int>
 
 #ifdef CGAL_MESH_3_DEMO_ACTIVATE_SEGMENTED_IMAGES
-typedef CGAL::Labeled_image_mesh_domain_3<Image,Kernel>           Image_domain;
+typedef CGAL::Labeled_image_mesh_domain_3<Image,Kernel>           Image_domain1;
-typedef CGAL::Polyhedron_demo_labeled_mesh_domain_3<Image_domain> Image_mesh_domain;
+typedef CGAL::Polyhedron_demo_labeled_mesh_domain_3<Image_domain1> Image_domain;
+typedef CGAL::Mesh_domain_with_polyline_features_3<Image_domain> Image_mesh_domain;
 #endif
 #ifdef CGAL_MESH_3_DEMO_ACTIVATE_IMPLICIT_FUNCTIONS
 typedef Wrapper<Kernel>                                              Function_wrapper;

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 )

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 <CGAL/Mesh_3/mesh_standard_criteria.h>
+
+
+template <typename Tr, typename Domain, typename Visitor_>
+class Facet_extra_criterion :
+  public CGAL::Mesh_3::Abstract_criterion<Tr, Visitor_>
+{
+private:
+  typedef typename Tr::Facet Facet;
+
+  typedef CGAL::Mesh_3::Abstract_criterion<Tr,Visitor_> Base;
+  typedef typename Base::Quality Quality;
+  typedef typename Base::Badness Badness;
+
+  typedef Facet_extra_criterion<Tr,Domain,Visitor_> 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

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

Polyhedron/demo/Polyhedron/Plugins/Mesh_3_plugin/Mesh_3_plugin_cgal_code.cpp

@@ -2,6 +2,7 @@
 
 #include <CGAL/Mesh_criteria_3.h>
 #include <CGAL/make_mesh_3.h>
+#include <CGAL/polylines_to_protect.h>
 #include <CGAL/Bbox_3.h>
 
 #include <Polyhedron_type.h>
@@ -11,6 +12,7 @@
 #include <CGAL/Three/Scene_item.h>
 
 #include "Mesh_function.h"
+#include "Facet_extra_criterion.h"
 
 #include <CGAL/Timer.h>
 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<std::vector<Point_3> > 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<Tr,
+                                Image_mesh_domain,
+                                Mesh_criteria::Facet_criteria::Visitor>
+    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<C3t3>(*p_domain,

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()];