signed distance and isosurface meshing

Former-commit-id: 9164537bc8da630155f50387952fd3cfe166f295

include/igl/cgal/complex_to_mesh.cpp

@@ -0,0 +1,150 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "complex_to_mesh.h"
+
+#include <igl/centroid.h>
+#include <igl/remove_unreferenced.h>
+
+#include <CGAL/Surface_mesh_default_triangulation_3.h>
+#include <set>
+#include <stack>
+
+template <typename Tr, typename DerivedV, typename DerivedF>
+IGL_INLINE bool igl::complex_to_mesh(
+  const CGAL::Complex_2_in_triangulation_3<Tr> & c2t3,
+  Eigen::PlainObjectBase<DerivedV> & V, 
+  Eigen::PlainObjectBase<DerivedF> & F)
+{
+  using namespace igl;
+  using namespace Eigen;
+  // CGAL/IO/Complex_2_in_triangulation_3_file_writer.h
+  using CGAL::Surface_mesher::number_of_facets_on_surface;
+
+  typedef typename CGAL::Complex_2_in_triangulation_3<Tr> C2t3;
+  typedef typename Tr::Finite_facets_iterator Finite_facets_iterator;
+  typedef typename Tr::Finite_vertices_iterator Finite_vertices_iterator;
+  typedef typename Tr::Facet Facet;
+  typedef typename Tr::Edge Edge;
+  typedef typename Tr::Vertex_handle Vertex_handle;
+
+  // Header.
+  const Tr& tr = c2t3.triangulation();
+
+  bool success = true;
+
+  const int n = tr.number_of_vertices();
+  const int m = c2t3.number_of_facets();
+
+  assert(m == number_of_facets_on_surface(tr));
+  
+  // Finite vertices coordinates.
+  std::map<Vertex_handle, int> v2i;
+  V.resize(n,3);
+  {
+    int v = 0;
+    for(Finite_vertices_iterator vit = tr.finite_vertices_begin();
+        vit != tr.finite_vertices_end();
+        ++vit)
+    {
+      V(v,0) = vit->point().x(); 
+      V(v,1) = vit->point().y(); 
+      V(v,2) = vit->point().z(); 
+      v2i[vit] = v++;
+    }
+  }
+
+  {
+    Finite_facets_iterator fit = tr.finite_facets_begin();
+    std::set<Facet> oriented_set;
+    std::stack<Facet> stack;
+
+    while ((int)oriented_set.size() != m) 
+    {
+      while ( fit->first->is_facet_on_surface(fit->second) == false ||
+        oriented_set.find(*fit) != oriented_set.end() ||
+        oriented_set.find(c2t3.opposite_facet(*fit)) !=
+        oriented_set.end() ) 
+      {
+        ++fit;
+      }
+      oriented_set.insert(*fit);
+      stack.push(*fit);
+      while(! stack.empty() )
+      {
+        Facet f = stack.top();
+        stack.pop();
+        for(int ih = 0 ; ih < 3 ; ++ih)
+        {
+          const int i1  = tr.vertex_triple_index(f.second, tr. cw(ih));
+          const int i2  = tr.vertex_triple_index(f.second, tr.ccw(ih));
+
+          const typename C2t3::Face_status face_status
+            = c2t3.face_status(Edge(f.first, i1, i2));
+          if(face_status == C2t3::REGULAR) 
+          {
+            Facet fn = c2t3.neighbor(f, ih);
+            if (oriented_set.find(fn) == oriented_set.end())
+            {
+              if(oriented_set.find(c2t3.opposite_facet(fn)) == oriented_set.end())
+              {
+                oriented_set.insert(fn);
+                stack.push(fn);
+              }else {
+                success = false; // non-orientable
+              }
+            }
+          }else if(face_status != C2t3::BOUNDARY)
+          {
+            success = false; // non manifold, thus non-orientable
+          }
+        } // end "for each neighbor of f"
+      } // end "stack non empty"
+    } // end "oriented_set not full"
+
+    F.resize(m,3);
+    int f = 0;
+    for(typename std::set<Facet>::const_iterator fit = 
+        oriented_set.begin();
+        fit != oriented_set.end();
+        ++fit)
+    {
+      const typename Tr::Cell_handle cell = fit->first;
+      const int& index = fit->second;
+      const int index1 = v2i[cell->vertex(tr.vertex_triple_index(index, 0))];
+      const int index2 = v2i[cell->vertex(tr.vertex_triple_index(index, 1))];
+      const int index3 = v2i[cell->vertex(tr.vertex_triple_index(index, 2))];
+      // This order is flipped
+      F(f,0) = index1;
+      F(f,1) = index2;
+      F(f,2) = index3;
+      f++;
+    }
+    assert(f == m);
+  } // end if(facets must be oriented)
+
+  // This CGAL code seems to randomly assign the global orientation
+  // Flip based on the signed volume.
+  Eigen::Vector3d cen;
+  double vol;
+  igl::centroid(V,F,cen,vol);
+  if(vol < 0)
+  {
+    // Flip
+    F = F.rowwise().reverse().eval();
+  }
+
+  // CGAL code somehow can end up with unreferenced vertices
+  VectorXi _1;
+  remove_unreferenced( MatrixXd(V), MatrixXi(F), V,F,_1);
+
+  return success;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+template bool igl::complex_to_mesh<CGAL::Delaunay_triangulation_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_data_structure_3<CGAL::Surface_mesh_vertex_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_vertex_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_ds_vertex_base_3<void> > >, CGAL::Triangulation_cell_base_with_circumcenter_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Surface_mesh_cell_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_cell_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_ds_cell_base_3<void> > > > >, CGAL::Default>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(CGAL::Complex_2_in_triangulation_3<CGAL::Delaunay_triangulation_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_data_structure_3<CGAL::Surface_mesh_vertex_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_vertex_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_ds_vertex_base_3<void> > >, CGAL::Triangulation_cell_base_with_circumcenter_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Surface_mesh_cell_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_cell_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_ds_cell_base_3<void> > > > >, CGAL::Default>, void> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+#endif

include/igl/cgal/complex_to_mesh.h

@@ -0,0 +1,41 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_COMPLEX_TO_MESH_H
+#define IGL_COMPLEX_TO_MESH_H
+#include "../igl_inline.h"
+
+#include <Eigen/Dense>
+#include <CGAL/Complex_2_in_triangulation_3.h>
+
+namespace igl 
+{
+  // Templates:
+  //   Tr  CGAL triangulation type, e.g.
+  //     CGAL::Surface_mesh_default_triangulation_3
+  // Inputs
+  //   c2t3  2-complex (surface) living in a 3d triangulation (e.g. result of
+  //     CGAL::make_surface_mesh)
+  // Outputs:
+  //   V  #V by 3 list of vertex positions
+  //   F  #F by 3 list of triangle indices
+  // Returns true iff conversion was successful, failure can ok if CGAL code
+  // can't figure out ordering.
+  //
+  template <typename Tr, typename DerivedV, typename DerivedF>
+  IGL_INLINE bool complex_to_mesh(
+    const CGAL::Complex_2_in_triangulation_3<Tr> & c2t3,
+    Eigen::PlainObjectBase<DerivedV> & V, 
+    Eigen::PlainObjectBase<DerivedF> & F);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "complex_to_mesh.cpp"
+#endif
+
+#endif
+

include/igl/cgal/signed_distance.cpp

@@ -0,0 +1,125 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "signed_distance.h"
+template <typename Kernel>
+IGL_INLINE typename Kernel::FT igl::signed_distance_pseudonormal(
+  const CGAL::AABB_tree<
+    CGAL::AABB_traits<Kernel, 
+      CGAL::AABB_triangle_primitive<Kernel, 
+        typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+      >
+    >
+  > & tree,
+  const std::vector<CGAL::Triangle_3<Kernel> > & T,
+  const Eigen::MatrixXi & F,
+  const Eigen::MatrixXd & FN,
+  const Eigen::MatrixXd & VN,
+  const Eigen::MatrixXd & EN,
+  const Eigen::VectorXi & EMAP,
+  const typename Kernel::Point_3 & q)
+{
+  using namespace Eigen;
+  using namespace std;
+  typedef typename Kernel::FT FT;
+  typedef typename Kernel::Point_3 Point_3;
+  typedef typename CGAL::Triangle_3<Kernel> Triangle_3;
+  typedef typename std::vector<Triangle_3>::iterator Iterator;
+  typedef typename CGAL::AABB_triangle_primitive<Kernel, Iterator> Primitive;
+  typedef typename CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
+  typedef typename CGAL::AABB_tree<AABB_triangle_traits> Tree;
+  typedef typename Tree::Point_and_primitive_id Point_and_primitive_id;
+
+  Point_and_primitive_id pp = tree.closest_point_and_primitive(q);
+  Point_3 & p = pp.first;
+  const auto & qp = q-p;
+  const FT sqrd = qp.squared_length();
+  Vector3d v(qp.x(),qp.y(),qp.z());
+  const int f = pp.second - T.begin();
+  const Triangle_3 & t = *pp.second;
+  // barycentric coordinates
+  const auto & area = [&p,&t](const int i, const int j)->FT
+  {
+    return sqrt(Triangle_3(p,t.vertex(i),t.vertex(j)).squared_area());
+  };
+  Vector3d b(area(1,2),area(2,0),area(0,1));
+  b /= b.sum();
+  // Determine which normal to use
+  Vector3d n;
+  const double epsilon = 1e-12;
+  const int type = (b.array()<=epsilon).cast<int>().sum();
+  switch(type)
+  {
+    case 2:
+      // Find vertex
+      for(int c = 0;c<3;c++)
+      {
+        if(b(c)>epsilon)
+        {
+          n = VN.row(F(f,c));
+          break;
+        }
+      }
+      break;
+    case 1:
+      // Find edge
+      for(int c = 0;c<3;c++)
+      {
+        if(b(c)<=epsilon)
+        {
+          n = EN.row(EMAP(F.rows()*c+f));
+          break;
+        }
+      }
+      break;
+    default:
+      assert(false && "all barycentric coords zero.");
+    case 0:
+      n = FN.row(f);
+      break;
+  }
+  const double s = (v.dot(n) >= 0 ? 1. : -1.);
+  return s*sqrt(sqrd);
+}
+
+template <typename Kernel>
+IGL_INLINE typename Kernel::FT igl::signed_distance_winding_number(
+  const CGAL::AABB_tree<
+    CGAL::AABB_traits<Kernel, 
+      CGAL::AABB_triangle_primitive<Kernel, 
+        typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+      >
+    >
+  > & tree,
+  const igl::WindingNumberAABB<Eigen::Vector3d> & hier,
+  const typename Kernel::Point_3 & q)
+{
+  typedef typename Kernel::FT FT;
+  typedef typename Kernel::Point_3 Point_3;
+  typedef typename CGAL::Triangle_3<Kernel> Triangle_3;
+  typedef typename std::vector<Triangle_3>::iterator Iterator;
+  typedef typename CGAL::AABB_triangle_primitive<Kernel, Iterator> Primitive;
+  typedef typename CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
+  typedef typename CGAL::AABB_tree<AABB_triangle_traits> Tree;
+  typedef typename Tree::Point_and_primitive_id Point_and_primitive_id;
+  using namespace Eigen;
+  using namespace std;
+  using namespace igl;
+  Point_and_primitive_id pp = tree.closest_point_and_primitive(q);
+  Point_3 & p = pp.first;
+  const auto & qp = q-p;
+  const Vector3d eq(q.x(),q.y(),q.z()); 
+  const FT sqrd = qp.squared_length();
+  const double w = hier.winding_number(eq);
+  const FT s = 1.-2.*w;
+  return s*sqrt(sqrd);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+template CGAL::Epick::FT igl::signed_distance_winding_number<CGAL::Epick>(CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Epick, CGAL::AABB_triangle_primitive<CGAL::Epick, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >::iterator, CGAL::Boolean_tag<false> > > > const&, igl::WindingNumberAABB<Eigen::Matrix<double, 3, 1, 0, 3, 1> > const&, CGAL::Epick::Point_3 const&);
+template CGAL::Epick::FT igl::signed_distance_pseudonormal<CGAL::Epick>(CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Epick, CGAL::AABB_triangle_primitive<CGAL::Epick, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >::iterator, CGAL::Boolean_tag<false> > > > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, CGAL::Epick::Point_3 const&);
+#endif

include/igl/cgal/signed_distance.h

@@ -0,0 +1,79 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_SIGNED_DISTANCE_H
+#define IGL_SIGNED_DISTANCE_H
+#include <igl/igl_inline.h>
+#include <igl/WindingNumberAABB.h>
+#include <Eigen/Core>
+#include <vector>
+#include "CGAL_includes.hpp"
+namespace igl
+{
+  enum SignedDistanceType
+  {
+    SIGNED_DISTANCE_TYPE_PSEUDONORMAL   = 0,
+    SIGNED_DISTANCE_TYPE_WINDING_NUMBER = 1,
+    SIGNED_DISTANCE_TYPE_DEFAULT        = 2,
+    NUM_SIGNED_DISTANCE_TYPE            = 3
+  };
+  // Computes signed distance to mesh
+  //
+  // Templates:
+  //   Kernal  CGAL computation and construction kernel (e.g.
+  //     CGAL::Simple_cartesian<double>)
+  // Inputs:
+  //   tree  AABB acceleration tree (see point_mesh_squared_distance.h)
+  //   T  #F list of CGAL triangles (see point_mesh_squared_distance.h)
+  //   F  #F by 3 list of triangle indices
+  //   FN  #F by 3 list of triangle normals 
+  //   VN  #V by 3 list of vertex normals (ANGLE WEIGHTING)
+  //   EN  #E by 3 list of edge normals (UNIFORM WEIGHTING)
+  //   EMAP  #F*3 mapping edges in F to E
+  //   q  Query point
+  // Returns signed distance to mesh
+  //
+  template <typename Kernel>
+  IGL_INLINE typename Kernel::FT signed_distance_pseudonormal(
+    const CGAL::AABB_tree<
+      CGAL::AABB_traits<Kernel, 
+        CGAL::AABB_triangle_primitive<Kernel, 
+          typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+        >
+      >
+    > & tree,
+    const std::vector<CGAL::Triangle_3<Kernel> > & T,
+    const Eigen::MatrixXi & F,
+    const Eigen::MatrixXd & FN,
+    const Eigen::MatrixXd & VN,
+    const Eigen::MatrixXd & EN,
+    const Eigen::VectorXi & EMAP,
+    const typename Kernel::Point_3 & q);
+  // Inputs:
+  //   tree  AABB acceleration tree (see point_mesh_squared_distance.h)
+  //   hier  Winding number evaluation hierarchy
+  //   q  Query point
+  // Returns signed distance to mesh
+  template <typename Kernel>
+  IGL_INLINE typename Kernel::FT signed_distance_winding_number(
+    const CGAL::AABB_tree<
+      CGAL::AABB_traits<Kernel, 
+        CGAL::AABB_triangle_primitive<Kernel, 
+          typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+        >
+      >
+    > & tree,
+    const igl::WindingNumberAABB<Eigen::Vector3d> & hier,
+    const typename Kernel::Point_3 & q);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "signed_distance.cpp"
+#endif
+
+#endif
+

include/igl/cgal/signed_distance_isosurface.cpp

@@ -0,0 +1,136 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "signed_distance_isosurface.h"
+#include "point_mesh_squared_distance.h"
+#include "complex_to_mesh.h"
+#include "signed_distance.h"
+
+#include <igl/per_face_normals.h>
+#include <igl/per_edge_normals.h>
+#include <igl/per_vertex_normals.h>
+#include <igl/centroid.h>
+#include <igl/WindingNumberAABB.h>
+#include <igl/matlab_format.h>
+#include <igl/remove_unreferenced.h>
+
+#include <CGAL/Surface_mesh_default_triangulation_3.h>
+#include <CGAL/Complex_2_in_triangulation_3.h>
+#include <CGAL/make_surface_mesh.h>
+#include <CGAL/Implicit_surface_3.h>
+#include <CGAL/Polyhedron_3.h>
+#include <CGAL/IO/output_surface_facets_to_polyhedron.h>
+// Axis-aligned bounding box tree for tet tri intersection
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_triangle_primitive.h>
+#include <vector>
+
+IGL_INLINE bool igl::signed_distance_isosurface(
+  const Eigen::MatrixXd & IV,
+  const Eigen::MatrixXi & IF,
+  const double level,
+  const double angle_bound,
+  const double radius_bound,
+  const double distance_bound,
+  const SignedDistanceType sign_type,
+  Eigen::MatrixXd & V,
+  Eigen::MatrixXi & F)
+{
+  using namespace std;
+
+  // default triangulation for Surface_mesher
+  typedef CGAL::Surface_mesh_default_triangulation_3 Tr;
+  // c2t3
+  typedef CGAL::Complex_2_in_triangulation_3<Tr> C2t3;
+  typedef Tr::Geom_traits GT;//Kernel
+  typedef GT::Sphere_3 Sphere_3;
+  typedef GT::Point_3 Point_3;
+  typedef GT::FT FT;
+  typedef std::function<FT (Point_3)> Function;
+  typedef CGAL::Implicit_surface_3<GT, Function> Surface_3;
+  typedef CGAL::Polyhedron_3<GT> Polyhedron;
+  typedef GT::Kernel Kernel;
+  typedef CGAL::Triangle_3<Kernel> Triangle_3; 
+  typedef typename std::vector<Triangle_3>::iterator Iterator;
+  typedef CGAL::AABB_triangle_primitive<Kernel, Iterator> Primitive;
+  typedef CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
+  typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
+  typedef typename Tree::Point_and_primitive_id Point_and_primitive_id;
+
+  Eigen::MatrixXd FN,VN,EN;
+  Eigen::MatrixXi E;
+  Eigen::VectorXi EMAP;
+  // "Signed Distance Computation Using the Angle Weighted Pseudonormal"
+  // [Bærentzen & Aanæs 2005]
+  per_face_normals(IV,IF,FN);
+  per_vertex_normals(IV,IF,PER_VERTEX_NORMALS_WEIGHTING_TYPE_ANGLE,VN);
+  per_edge_normals(IV,IF,PER_EDGE_NORMALS_WEIGHTING_TYPE_UNIFORM,EN,E,EMAP);
+  // Prepare distance computation
+  Tree tree;
+  vector<Triangle_3 > T;
+  point_mesh_squared_distance_precompute(IV,IF,tree,T);
+  WindingNumberAABB<Eigen::Vector3d> hier(IV,IF);
+  hier.grow();
+  Tr tr;            // 3D-Delaunay triangulation
+  C2t3 c2t3 (tr);   // 2D-complex in 3D-Delaunay triangulation
+  // defining the surface
+  const auto & IVmax = IV.colwise().maxCoeff();
+  const auto & IVmin = IV.colwise().minCoeff();
+  const double bbd = (IVmax-IVmin).norm();
+  const double r = bbd/2.;
+  const auto & IVmid = 0.5*(IVmax + IVmin);
+  // Supposedly the implict needs to evaluate to <0 at cmid...
+  // http://doc.cgal.org/latest/Surface_mesher/classCGAL_1_1Implicit__surface__3.html
+  Point_3 cmid(IVmid(0),IVmid(1),IVmid(2));
+  Function fun;
+  switch(sign_type)
+  {
+    default:
+      assert(false && "Unknown SignedDistanceType");
+    case SIGNED_DISTANCE_TYPE_DEFAULT:
+    case SIGNED_DISTANCE_TYPE_WINDING_NUMBER:
+      fun = 
+        [&tree,&hier,&level](const Point_3 & q) -> FT
+        {
+          return signed_distance_winding_number(tree,hier,q)-level;
+        };
+      break;
+    case SIGNED_DISTANCE_TYPE_PSEUDONORMAL:
+      fun = [&tree,&T,&IF,&FN,&VN,&EN,&EMAP,&level](const Point_3 & q) -> FT
+        {
+          return 
+            igl::signed_distance_pseudonormal(tree,T,IF,FN,VN,EN,EMAP,q) - 
+            level;
+        };
+      break;
+  }
+    //[&tree,&hier,&T,&IF,&FN,&VN,&EN,&EMAP,&level](const Point_3 & q) ->FT
+    //{
+    //  const FT p = signed_distance_pseudonormal(tree,T,IF,FN,VN,EN,EMAP,q);
+    //  const FT w = signed_distance_winding_number(tree,hier,q);
+    //  if(w*p < 0 && (fabs(w) > 0.1 || fabs(p) > 0.1))
+    //  {
+    //    cout<<"q=["<<q.x()<<","<<q.y()<<","<<q.z()<<"];"<<endl;
+    //    cout<<matlab_format(n.transpose().eval(),"n")<<endl;
+    //    cout<<matlab_format(b.transpose().eval(),"b")<<endl;
+    //    cout<<"Sig difference: "<<type<<endl;
+    //    cout<<"w: "<<w<<endl;
+    //    cout<<"p: "<<p<<endl;
+    //    exit(1);
+    //  }
+    //  return w;
+    //},
+  Sphere_3 bounding_sphere(cmid, (r+level)*(r+level));
+  Surface_3 surface(fun,bounding_sphere);
+  CGAL::Surface_mesh_default_criteria_3<Tr> 
+    criteria(angle_bound,radius_bound,distance_bound);
+  // meshing surface
+  CGAL::make_surface_mesh(c2t3, surface, criteria, CGAL::Manifold_tag());
+  // complex to (V,F)
+  return igl::complex_to_mesh(c2t3,V,F);
+}

include/igl/cgal/signed_distance_isosurface.h

@@ -0,0 +1,44 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_SIGNED_DISTANCE_ISOSURFACE_H
+#define IGL_SIGNED_DISTANCE_ISOSURFACE_H
+#include <igl/igl_inline.h>
+#include "signed_distance.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Inputs:
+  //   IV  #IV by 3 list of input mesh vertex positions
+  //   IF  #IF by 3 list of input triangle indices
+  //   level  iso-level to contour (e.g. 0)
+  //   angle_bound  lower bound on triangle angles (mesh quality) (e.g. 28)
+  //   radius_bound  upper bound on triangle size (mesh density?) (e.g. 0.02)
+  //   distance_bound  cgal mysterious parameter (mesh density?) (e.g. 0.01)
+  //   sign_type  method for computing distance _sign_ (see signed_distance.h)
+  // Outputs:
+  //   V  #V by 3 list of input mesh vertex positions
+  //   F  #F by 3 list of input triangle indices
+  //  
+  IGL_INLINE bool signed_distance_isosurface(
+    const Eigen::MatrixXd & IV,
+    const Eigen::MatrixXi & IF,
+    const double level,
+    const double angle_bound,
+    const double radius_bound,
+    const double distance_bound,
+    const SignedDistanceType sign_type,
+    Eigen::MatrixXd & V,
+    Eigen::MatrixXi & F);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "signed_distance_isosurface.cpp"
+#endif
+
+#endif
+