Merge pull request #363 from qnzhou/master

Triangulation

Former-commit-id: adb01e11ac5456a926287cc299ea375aac49d2ab

include/igl/copyleft/cgal/delaunay_triangulation.cpp

@@ -0,0 +1,24 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Qingnan Zhou <qnzhou@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 "delaunay_triangulation.h"
+#include "../../delaunay_triangulation.h"
+#include "orient2D.h"
+#include "incircle.h"
+
+template<
+  typename DerivedV,
+  typename DerivedF>
+IGL_INLINE void igl::copyleft::cgal::delaunay_triangulation(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    Eigen::PlainObjectBase<DerivedF>& F)
+{
+  typedef typename DerivedV::Scalar Scalar;
+  igl::delaunay_triangulation(V, orient2D<Scalar>, incircle<Scalar>, F);
+}
+

include/igl/copyleft/cgal/delaunay_triangulation.h

@@ -0,0 +1,47 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Qingan Zhou <qnzhou@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_COPYLEFT_CGAL_DELAUNAY_TRIANGULATION_H
+#define IGL_COPYLEFT_CGAL_DELAUNAY_TRIANGULATION_H
+
+#include "../../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+
+      // Given a set of points in 2D, return a Delaunay triangulation of these
+      // points.
+      //
+      // Inputs:
+      //   V  #V by 2 list of vertex positions
+      //
+      // Outputs:
+      //   F  #F by 3 of faces in Delaunay triangulation.
+      template<
+        typename DerivedV,
+        typename DerivedF
+        >
+      IGL_INLINE void delaunay_triangulation(
+          const Eigen::PlainObjectBase<DerivedV>& V,
+          Eigen::PlainObjectBase<DerivedF>& F);
+    }
+  }
+}
+
+
+
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "delaunay_triangulation.cpp"
+#endif
+#endif

include/igl/copyleft/cgal/incircle.cpp

@@ -0,0 +1,39 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Qingan Zhou <qnzhou@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 "incircle.h"
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+template<typename Scalar>
+IGL_INLINE short igl::copyleft::cgal::incircle(
+    const Scalar pa[2],
+    const Scalar pb[2],
+    const Scalar pc[2],
+    const Scalar pd[2])
+{
+  typedef CGAL::Exact_predicates_exact_constructions_kernel Epeck;
+  typedef CGAL::Exact_predicates_inexact_constructions_kernel Epick;
+  typedef typename std::conditional<std::is_same<Scalar, Epeck::FT>::value,
+          Epeck, Epick>::type Kernel;
+
+  switch(CGAL::side_of_oriented_circle(
+        typename Kernel::Point_2(pa[0], pa[1]),
+        typename Kernel::Point_2(pb[0], pb[1]),
+        typename Kernel::Point_2(pc[0], pc[1]),
+        typename Kernel::Point_2(pd[0], pd[1]))) {
+    case CGAL::ON_POSITIVE_SIDE:
+      return 1;
+    case CGAL::ON_NEGATIVE_SIDE:
+      return -1;
+    case CGAL::ON_ORIENTED_BOUNDARY:
+      return 0;
+    default:
+      throw "Invalid incircle result";
+  }
+}

include/igl/copyleft/cgal/incircle.h

@@ -0,0 +1,39 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Qingan Zhou <qnzhou@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_COPYLEFT_CGAL_INCIRCLE_H
+#define IGL_COPYLEFT_CGAL_INCIRCLE_H
+
+#include "../../igl_inline.h"
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // Inputs:
+      //   pa,pb,pc,pd  2D points.
+      // Output:
+      //   1 if pd is inside of the oriented circle formed by pa,pb,pc.
+      //   0 if pd is co-circular with pa,pb,pc.
+      //  -1 if pd is outside of the oriented circle formed by pa,pb,pc.
+      template <typename Scalar>
+      IGL_INLINE short incircle(
+          const Scalar pa[2],
+          const Scalar pb[2],
+          const Scalar pc[2],
+          const Scalar pd[2]);
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "incircle.cpp"
+#endif
+#endif

include/igl/copyleft/cgal/insphere.cpp

@@ -0,0 +1,41 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Qingan Zhou <qnzhou@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 "insphere.h"
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+template<typename Scalar>
+IGL_INLINE short igl::copyleft::cgal::insphere(
+    const Scalar pa[3],
+    const Scalar pb[3],
+    const Scalar pc[3],
+    const Scalar pd[3],
+    const Scalar pe[3])
+{
+  typedef CGAL::Exact_predicates_exact_constructions_kernel Epeck;
+  typedef CGAL::Exact_predicates_inexact_constructions_kernel Epick;
+  typedef typename std::conditional<std::is_same<Scalar, Epeck::FT>::value,
+          Epeck, Epick>::type Kernel;
+
+  switch(CGAL::side_of_oriented_sphere(
+        typename Kernel::Point_3(pa[0], pa[1], pa[2]),
+        typename Kernel::Point_3(pb[0], pb[1], pb[2]),
+        typename Kernel::Point_3(pc[0], pc[1], pc[2]),
+        typename Kernel::Point_3(pd[0], pd[1], pd[2]),
+        typename Kernel::Point_3(pe[0], pe[1], pe[2]))) {
+    case CGAL::ON_POSITIVE_SIDE:
+      return 1;
+    case CGAL::ON_NEGATIVE_SIDE:
+      return -1;
+    case CGAL::ON_ORIENTED_BOUNDARY:
+      return 0;
+    default:
+      throw "Invalid incircle result";
+  }
+}

include/igl/copyleft/cgal/insphere.h

@@ -0,0 +1,40 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Qingan Zhou <qnzhou@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_COPYLEFT_CGAL_INSPHERE_H
+#define IGL_COPYLEFT_CGAL_INSPHERE_H
+
+#include "../../igl_inline.h"
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // Inputs:
+      //   pa,pb,pc,pd,pe  3D points.
+      // Output:
+      //   1 if pe is inside of the oriented sphere formed by pa,pb,pc,pd.
+      //   0 if pe is co-spherical with pa,pb,pc,pd.
+      //  -1 if pe is outside of the oriented sphere formed by pa,pb,pc,pd.
+      template <typename Scalar>
+      IGL_INLINE short insphere(
+          const Scalar pa[3],
+          const Scalar pb[3],
+          const Scalar pc[3],
+          const Scalar pd[3],
+          const Scalar pe[3]);
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "insphere.cpp"
+#endif
+#endif

include/igl/copyleft/cgal/lexicographic_triangulation.cpp

@@ -0,0 +1,24 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Alec Jacobson <alecjacobson@gmail.com>
+//                    Qingan Zhou <qnzhou@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 "lexicographic_triangulation.h"
+#include "../../lexicographic_triangulation.h"
+#include "orient2D.h"
+
+template<
+  typename DerivedP,
+  typename DerivedF
+  >
+IGL_INLINE void igl::copyleft::cgal::lexicographic_triangulation(
+    const Eigen::PlainObjectBase<DerivedP>& P,
+    Eigen::PlainObjectBase<DerivedF>& F)
+{
+  typedef typename DerivedP::Scalar Scalar;
+  igl::lexicographic_triangulation(P, orient2D<Scalar>, F);
+}

include/igl/copyleft/cgal/lexicographic_triangulation.h

@@ -0,0 +1,47 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Alec Jacobson <alecjacobson@gmail.com>
+//                    Qingan Zhou <qnzhou@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_COPYLEFT_CGAL_LEXICOGRAPHIC_TRIANGULATION_H
+#define IGL_COPYLEFT_CGAL_LEXICOGRAPHIC_TRIANGULATION_H
+#include "../../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+
+      // Given a set of points in 2D, return a lexicographic triangulation of these
+      // points.
+      //
+      // Inputs:
+      //   P  #P by 2 list of vertex positions
+      //
+      // Outputs:
+      //   F  #F by 3 of faces in lexicographic triangulation.
+      template<
+        typename DerivedP,
+        typename DerivedF
+        >
+      IGL_INLINE void lexicographic_triangulation(
+          const Eigen::PlainObjectBase<DerivedP>& P,
+          Eigen::PlainObjectBase<DerivedF>& F);
+    }
+  }
+}
+
+
+
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "lexicographic_triangulation.cpp"
+#endif
+#endif

include/igl/copyleft/cgal/orient2D.cpp

@@ -0,0 +1,37 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Qingan Zhou <qnzhou@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 "orient2D.h"
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+template<typename Scalar>
+IGL_INLINE short igl::copyleft::cgal::orient2D(
+    const Scalar pa[2],
+    const Scalar pb[2],
+    const Scalar pc[2])
+{
+  typedef CGAL::Exact_predicates_exact_constructions_kernel Epeck;
+  typedef CGAL::Exact_predicates_inexact_constructions_kernel Epick;
+  typedef typename std::conditional<std::is_same<Scalar, Epeck::FT>::value,
+          Epeck, Epick>::type Kernel;
+
+  switch(CGAL::orientation(
+        typename Kernel::Point_2(pa[0], pa[1]),
+        typename Kernel::Point_2(pb[0], pb[1]),
+        typename Kernel::Point_2(pc[0], pc[1]))) {
+    case CGAL::LEFT_TURN:
+      return 1;
+    case CGAL::RIGHT_TURN:
+      return -1;
+    case CGAL::COLLINEAR:
+      return 0;
+    default:
+      throw "Invalid orientation";
+  }
+}

include/igl/copyleft/cgal/orient2D.h

@@ -0,0 +1,38 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Qingan Zhou <qnzhou@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_COPYLEFT_CGAL_ORIENT_2D_H
+#define IGL_COPYLEFT_CGAL_ORIENT_2D_H
+
+#include "../../igl_inline.h"
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // Inputs:
+      //   pa,pb,pc   2D points.
+      // Output:
+      //   1 if pa,pb,pc are counterclockwise oriented.
+      //   0 if pa,pb,pc are counterclockwise oriented.
+      //  -1 if pa,pb,pc are clockwise oriented.
+      template <typename Scalar>
+      IGL_INLINE short orient2D(
+          const Scalar pa[2],
+          const Scalar pb[2],
+          const Scalar pc[2]);
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "orient2D.cpp"
+#endif
+#endif

include/igl/copyleft/cgal/orient3D.cpp

@@ -0,0 +1,39 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Qingan Zhou <qnzhou@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 "orient3D.h"
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+template<typename Scalar>
+IGL_INLINE short igl::copyleft::cgal::orient3D(
+    const Scalar pa[3],
+    const Scalar pb[3],
+    const Scalar pc[3],
+    const Scalar pd[3])
+{
+  typedef CGAL::Exact_predicates_exact_constructions_kernel Epeck;
+  typedef CGAL::Exact_predicates_inexact_constructions_kernel Epick;
+  typedef typename std::conditional<std::is_same<Scalar, Epeck::FT>::value,
+          Epeck, Epick>::type Kernel;
+
+  switch(CGAL::orientation(
+        typename Kernel::Point_3(pa[0], pa[1], pa[2]),
+        typename Kernel::Point_3(pb[0], pb[1], pb[2]),
+        typename Kernel::Point_3(pc[0], pc[1], pc[2]),
+        typename Kernel::Point_3(pd[0], pd[1], pd[2]))) {
+    case CGAL::POSITIVE:
+      return 1;
+    case CGAL::NEGATIVE:
+      return -1;
+    case CGAL::COPLANAR:
+      return 0;
+    default:
+      throw "Invalid orientation";
+  }
+}

include/igl/copyleft/cgal/orient3D.h

@@ -0,0 +1,39 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Qingan Zhou <qnzhou@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_COPYLEFT_CGAL_ORIENT_3D_H
+#define IGL_COPYLEFT_CGAL_ORIENT_3D_H
+
+#include "../../igl_inline.h"
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // Inputs:
+      //   pa,pb,pc,pd  3D points.
+      // Output:
+      //   1 if pa,pb,pc,pd forms a tet of positive volume.
+      //   0 if pa,pb,pc,pd are coplanar.
+      //  -1 if pa,pb,pc,pd forms a tet of negative volume.
+      template <typename Scalar>
+      IGL_INLINE short orient3D(
+          const Scalar pa[3],
+          const Scalar pb[3],
+          const Scalar pc[3],
+          const Scalar pd[3]);
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "orient3D.cpp"
+#endif
+#endif

include/igl/delaunay_triangulation.cpp

@@ -0,0 +1,83 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Qingnan Zhou <qnzhou@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 "delaunay_triangulation.h"
+#include "flip_edge.h"
+#include "lexicographic_triangulation.h"
+#include "unique_edge_map.h"
+
+#include <vector>
+#include <sstream>
+
+template<
+  typename DerivedV,
+  typename Orient2D,
+  typename InCircle,
+  typename DerivedF>
+IGL_INLINE void igl::delaunay_triangulation(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    Orient2D orient2D,
+    InCircle incircle,
+    Eigen::PlainObjectBase<DerivedF>& F)
+{
+  assert(V.cols() == 2);
+  typedef typename DerivedF::Scalar Index;
+  typedef typename DerivedV::Scalar Scalar;
+  igl::lexicographic_triangulation(V, orient2D, F);
+  const size_t num_faces = F.rows();
+  if (num_faces == 0) {
+    // Input points are degenerate.  No faces will be generated.
+    return;
+  }
+  assert(F.cols() == 3);
+
+  Eigen::MatrixXi E;
+  Eigen::MatrixXi uE;
+  Eigen::VectorXi EMAP;
+  std::vector<std::vector<Index> > uE2E;
+  igl::unique_edge_map(F, E, uE, EMAP, uE2E);
+
+  auto is_delaunay = [&V,&F,&uE2E,num_faces,&incircle](size_t uei) {
+    auto& half_edges = uE2E[uei];
+    if (half_edges.size() != 2) {
+      throw "Cannot flip non-manifold or boundary edge";
+    }
+
+    const size_t f1 = half_edges[0] % num_faces;
+    const size_t f2 = half_edges[1] % num_faces;
+    const size_t c1 = half_edges[0] / num_faces;
+    const size_t c2 = half_edges[1] / num_faces;
+    assert(c1 < 3);
+    assert(c2 < 3);
+    assert(f1 != f2);
+    const size_t v1 = F(f1, (c1+1)%3);
+    const size_t v2 = F(f1, (c1+2)%3);
+    const size_t v4 = F(f1, c1);
+    const size_t v3 = F(f2, c2);
+    const Scalar p1[] = {V(v1, 0), V(v1, 1)};
+    const Scalar p2[] = {V(v2, 0), V(v2, 1)};
+    const Scalar p3[] = {V(v3, 0), V(v3, 1)};
+    const Scalar p4[] = {V(v4, 0), V(v4, 1)};
+    auto orientation = incircle(p1, p2, p4, p3);
+    return orientation <= 0;
+  };
+
+  bool all_delaunay = false;
+  while(!all_delaunay) {
+    all_delaunay = true;
+    for (size_t i=0; i<uE2E.size(); i++) {
+      if (uE2E[i].size() == 2) {
+        if (!is_delaunay(i)) {
+          all_delaunay = false;
+          flip_edge(F, E, uE, EMAP, uE2E, i);
+        }
+      }
+    }
+  }
+}
+

include/igl/delaunay_triangulation.h

@@ -0,0 +1,52 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Qingan Zhou <qnzhou@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_DELAUNAY_TRIANGULATION_H
+#define IGL_DELAUNAY_TRIANGULATION_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  // Given a set of points in 2D, return a Delaunay triangulation of these
+  // points.
+  //
+  // Inputs:
+  //   V  #V by 2 list of vertex positions
+  //   orient2D  A functor such that orient2D(pa, pb, pc) returns
+  //               1 if pa,pb,pc forms a conterclockwise triangle.
+  //              -1 if pa,pb,pc forms a clockwise triangle.
+  //               0 if pa,pb,pc are collinear.
+  //              where the argument pa,pb,pc are of type Scalar[2].
+  //   incircle  A functor such that incircle(pa, pb, pc, pd) returns
+  //               1 if pd is on the positive size of circumcirle of (pa,pb,pc)
+  //              -1 if pd is on the positive size of circumcirle of (pa,pb,pc)
+  //               0 if pd is cocircular with pa, pb, pc.
+  // Outputs:
+  //   F  #F by 3 of faces in Delaunay triangulation.
+  template<
+    typename DerivedV,
+    typename Orient2D,
+    typename InCircle,
+    typename DerivedF
+    >
+  IGL_INLINE void delaunay_triangulation(
+      const Eigen::PlainObjectBase<DerivedV>& V,
+      Orient2D orient2D,
+      InCircle incircle,
+      Eigen::PlainObjectBase<DerivedF>& F);
+}
+
+
+
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "delaunay_triangulation.cpp"
+#endif
+#endif

include/igl/flip_edge.cpp

@@ -0,0 +1,148 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Qingan Zhou <qnzhou@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 "flip_edge.h"
+
+template <
+  typename DerivedF,
+  typename DerivedE,
+  typename DeriveduE,
+  typename DerivedEMAP,
+  typename uE2EType>
+IGL_INLINE void igl::flip_edge(
+  Eigen::PlainObjectBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedE> & E,
+  Eigen::PlainObjectBase<DeriveduE> & uE,
+  Eigen::PlainObjectBase<DerivedEMAP> & EMAP,
+  std::vector<std::vector<uE2EType> > & uE2E,
+  const size_t uei)
+{
+  typedef typename DerivedF::Scalar Index;
+  const size_t num_faces = F.rows();
+  assert(F.cols() == 3);
+  //          v1                 v1
+  //          /|\                / \
+  //         / | \              /f1 \
+  //     v3 /f2|f1\ v4  =>  v3 /_____\ v4
+  //        \  |  /            \ f2  /
+  //         \ | /              \   /
+  //          \|/                \ /
+  //          v2                 v2
+  auto& half_edges = uE2E[uei];
+  if (half_edges.size() != 2) {
+    throw "Cannot flip non-manifold or boundary edge";
+  }
+
+  const size_t f1 = half_edges[0] % num_faces;
+  const size_t f2 = half_edges[1] % num_faces;
+  const size_t c1 = half_edges[0] / num_faces;
+  const size_t c2 = half_edges[1] / num_faces;
+  assert(c1 < 3);
+  assert(c2 < 3);
+
+  assert(f1 != f2);
+  const size_t v1 = F(f1, (c1+1)%3);
+  const size_t v2 = F(f1, (c1+2)%3);
+  const size_t v4 = F(f1, c1);
+  const size_t v3 = F(f2, c2);
+  assert(F(f2, (c2+2)%3) == v1);
+  assert(F(f2, (c2+1)%3) == v2);
+
+  const size_t e_12 = half_edges[0];
+  const size_t e_24 = f1 + ((c1 + 1) % 3) * num_faces;
+  const size_t e_41 = f1 + ((c1 + 2) % 3) * num_faces;
+  const size_t e_21 = half_edges[1];
+  const size_t e_13 = f2 + ((c2 + 1) % 3) * num_faces;
+  const size_t e_32 = f2 + ((c2 + 2) % 3) * num_faces;
+  assert(E(e_12, 0) == v1);
+  assert(E(e_12, 1) == v2);
+  assert(E(e_24, 0) == v2);
+  assert(E(e_24, 1) == v4);
+  assert(E(e_41, 0) == v4);
+  assert(E(e_41, 1) == v1);
+  assert(E(e_21, 0) == v2);
+  assert(E(e_21, 1) == v1);
+  assert(E(e_13, 0) == v1);
+  assert(E(e_13, 1) == v3);
+  assert(E(e_32, 0) == v3);
+  assert(E(e_32, 1) == v2);
+
+  const size_t ue_24 = EMAP[e_24];
+  const size_t ue_41 = EMAP[e_41];
+  const size_t ue_13 = EMAP[e_13];
+  const size_t ue_32 = EMAP[e_32];
+
+  F(f1, 0) = v1;
+  F(f1, 1) = v3;
+  F(f1, 2) = v4;
+  F(f2, 0) = v2;
+  F(f2, 1) = v4;
+  F(f2, 2) = v3;
+
+  uE(uei, 0) = v3;
+  uE(uei, 1) = v4;
+
+  const size_t new_e_34 = f1;
+  const size_t new_e_41 = f1 + num_faces;
+  const size_t new_e_13 = f1 + num_faces*2;
+  const size_t new_e_43 = f2;
+  const size_t new_e_32 = f2 + num_faces;
+  const size_t new_e_24 = f2 + num_faces*2;
+
+  E(new_e_34, 0) = v3;
+  E(new_e_34, 1) = v4;
+  E(new_e_41, 0) = v4;
+  E(new_e_41, 1) = v1;
+  E(new_e_13, 0) = v1;
+  E(new_e_13, 1) = v3;
+  E(new_e_43, 0) = v4;
+  E(new_e_43, 1) = v3;
+  E(new_e_32, 0) = v3;
+  E(new_e_32, 1) = v2;
+  E(new_e_24, 0) = v2;
+  E(new_e_24, 1) = v4;
+
+  EMAP[new_e_34] = uei;
+  EMAP[new_e_43] = uei;
+  EMAP[new_e_41] = ue_41;
+  EMAP[new_e_13] = ue_13;
+  EMAP[new_e_32] = ue_32;
+  EMAP[new_e_24] = ue_24;
+
+  auto replace = [](std::vector<Index>& array, Index old_v, Index new_v) {
+    std::replace(array.begin(), array.end(), old_v, new_v);
+  };
+  replace(uE2E[uei], e_12, new_e_34);
+  replace(uE2E[uei], e_21, new_e_43);
+  replace(uE2E[ue_13], e_13, new_e_13);
+  replace(uE2E[ue_32], e_32, new_e_32);
+  replace(uE2E[ue_24], e_24, new_e_24);
+  replace(uE2E[ue_41], e_41, new_e_41);
+
+#ifndef NDEBUG
+  auto sanity_check = [&](size_t ue) {
+    const auto& adj_faces = uE2E[ue];
+    if (adj_faces.size() == 2) {
+      const size_t first_f  = adj_faces[0] % num_faces;
+      const size_t first_c  = adj_faces[0] / num_faces;
+      const size_t second_f = adj_faces[1] % num_faces;
+      const size_t second_c = adj_faces[1] / num_faces;
+      const size_t vertex_0 = F(first_f, (first_c+1) % 3);
+      const size_t vertex_1 = F(first_f, (first_c+2) % 3);
+      assert(vertex_0 == F(second_f, (second_c+2) % 3));
+      assert(vertex_1 == F(second_f, (second_c+1) % 3));
+    }
+  };
+
+  sanity_check(uei);
+  sanity_check(ue_13);
+  sanity_check(ue_32);
+  sanity_check(ue_24);
+  sanity_check(ue_41);
+#endif
+}

include/igl/flip_edge.h

@@ -0,0 +1,52 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Qingan Zhou <qnzhou@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_FLIP_EDGE_H
+#define IGL_FLIP_EDGE_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+
+namespace igl
+{
+  // Flip an edge in a triangle mesh.  The edge specified by uei must have
+  // exactly **two** adjacent faces.  Violation will result in exception.
+  // Another warning: edge flipping could convert manifold mesh into
+  // non-manifold.
+  //
+  // Inputs:
+  //   F    #F by 3 list of triangles.
+  //   E    #F*3 by 2 list of all of directed edges
+  //   uE   #uE by 2 list of unique undirected edges
+  //   EMAP #F*3 list of indices into uE, mapping each directed edge to unique
+  //        undirected edge
+  //   uE2E #uE list of lists of indices into E of coexisting edges
+  //   ue   index into uE the edge to be flipped.
+  //
+  // Output:
+  //   Updated F, E, uE, EMAP and uE2E.
+  template <
+    typename DerivedF,
+    typename DerivedE,
+    typename DeriveduE,
+    typename DerivedEMAP,
+    typename uE2EType>
+  IGL_INLINE void flip_edge(
+    Eigen::PlainObjectBase<DerivedF> & F,
+    Eigen::PlainObjectBase<DerivedE> & E,
+    Eigen::PlainObjectBase<DeriveduE> & uE,
+    Eigen::PlainObjectBase<DerivedEMAP> & EMAP,
+    std::vector<std::vector<uE2EType> > & uE2E,
+    const size_t uei);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "flip_edge.cpp"
+#endif
+#endif

include/igl/lexicographic_triangulation.cpp

@@ -0,0 +1,128 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Alec Jacobson <alecjacobson@gmail.com>
+//                    Qingan Zhou <qnzhou@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 "lexicographic_triangulation.h"
+#include "sortrows.h"
+
+#include <vector>
+#include <list>
+
+template<
+  typename DerivedP,
+  typename Orient2D,
+  typename DerivedF
+  >
+IGL_INLINE void igl::lexicographic_triangulation(
+    const Eigen::PlainObjectBase<DerivedP>& P,
+    Orient2D orient2D,
+    Eigen::PlainObjectBase<DerivedF>& F)
+{
+  typedef typename DerivedP::Scalar Scalar;
+  const size_t num_pts = P.rows();
+  if (num_pts < 3) {
+    throw "At least 3 points are required for triangulation!";
+  }
+
+  // Sort points in lexicographic order.
+  DerivedP ordered_P;
+  Eigen::VectorXi order;
+  igl::sortrows(P, true, ordered_P, order);
+
+  std::vector<Eigen::Vector3i> faces;
+  std::list<size_t> boundary;
+  const Scalar p0[] = {ordered_P(0, 0), ordered_P(0, 1)};
+  const Scalar p1[] = {ordered_P(1, 0), ordered_P(1, 1)};
+  for (size_t i=2; i<num_pts; i++) {
+    const Scalar curr_p[] = {ordered_P(i, 0), ordered_P(i, 1)};
+    if (faces.size() == 0) {
+      // All points processed so far are collinear.
+      // Check if the current point is collinear with every points before it.
+      auto orientation = orient2D(p0, p1, curr_p);
+      if (orientation != 0) {
+        // Add a fan of triangles eminating from curr_p.
+        if (orientation > 0) {
+          for (size_t j=0; j<=i-2; j++) {
+            faces.push_back({order[j], order[j+1], order[i]});
+          }
+        } else if (orientation < 0) {
+          for (size_t j=0; j<=i-2; j++) {
+            faces.push_back({order[j+1], order[j], order[i]});
+          }
+        }
+        // Initialize current boundary.
+        boundary.insert(boundary.end(), order.data(), order.data()+i+1);
+        if (orientation < 0) {
+          boundary.reverse();
+        }
+      }
+    } else {
+      const size_t bd_size = boundary.size();
+      assert(bd_size >= 3);
+      std::vector<short> orientations;
+      for (auto itr=boundary.begin(); itr!=boundary.end(); itr++) {
+        auto next_itr = std::next(itr, 1);
+        if (next_itr == boundary.end()) {
+          next_itr = boundary.begin();
+        }
+        const Scalar bd_p0[] = {P(*itr, 0), P(*itr, 1)};
+        const Scalar bd_p1[] = {P(*next_itr, 0), P(*next_itr, 1)};
+        auto orientation = orient2D(bd_p0, bd_p1, curr_p);
+        if (orientation < 0) {
+          faces.push_back({*next_itr, *itr, order[i]});
+        }
+        orientations.push_back(orientation);
+      }
+
+      auto left_itr = boundary.begin();
+      auto right_itr = boundary.begin();
+      auto curr_itr = boundary.begin();
+      for (size_t j=0; j<bd_size; j++, curr_itr++) {
+        size_t prev = (j+bd_size-1) % bd_size;
+        if (orientations[j] >= 0 && orientations[prev] < 0) {
+          right_itr = curr_itr;
+        } else if (orientations[j] < 0 && orientations[prev] >= 0) {
+          left_itr = curr_itr;
+        }
+      }
+      assert(left_itr != right_itr);
+
+      for (auto itr=left_itr; itr!=right_itr; itr++) {
+        if (itr == boundary.end()) itr = boundary.begin();
+        if (itr == right_itr) break;
+        if (itr == left_itr || itr == right_itr) continue;
+        itr = boundary.erase(itr);
+        if (itr == boundary.begin()) {
+            itr = boundary.end();
+        } else {
+            itr--;
+        }
+      }
+
+      if (right_itr == boundary.begin()) {
+        assert(std::next(left_itr, 1) == boundary.end());
+        boundary.insert(boundary.end(), order[i]);
+      } else {
+        assert(std::next(left_itr, 1) == right_itr);
+        boundary.insert(right_itr, order[i]);
+      }
+    }
+  }
+
+  const size_t num_faces = faces.size();
+  if (num_faces == 0) {
+      // All input points are collinear.
+      // Do nothing here.
+  } else {
+      F.resize(num_faces, 3);
+      for (size_t i=0; i<num_faces; i++) {
+          F.row(i) = faces[i];
+      }
+  }
+}
+

include/igl/lexicographic_triangulation.h

@@ -0,0 +1,47 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Alec Jacobson <alecjacobson@gmail.com>
+//                    Qingan Zhou <qnzhou@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_LEXICOGRAPHIC_TRIANGULATION_H
+#define IGL_LEXICOGRAPHIC_TRIANGULATION_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  // Given a set of points in 2D, return a lexicographic triangulation of these
+  // points.
+  //
+  // Inputs:
+  //   P  #P by 2 list of vertex positions
+  //   orient2D  A functor such that orient2D(pa, pb, pc) returns
+  //               1 if pa,pb,pc forms a conterclockwise triangle.
+  //              -1 if pa,pb,pc forms a clockwise triangle.
+  //               0 if pa,pb,pc are collinear.
+  //              where the argument pa,pb,pc are of type Scalar[2].
+  //
+  // Outputs:
+  //   F  #F by 3 of faces in lexicographic triangulation.
+  template<
+    typename DerivedP,
+    typename Orient2D,
+    typename DerivedF
+    >
+  IGL_INLINE void lexicographic_triangulation(
+      const Eigen::PlainObjectBase<DerivedP>& P,
+      Orient2D orient2D,
+      Eigen::PlainObjectBase<DerivedF>& F);
+}
+
+
+
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "lexicographic_triangulation.cpp"
+#endif
+#endif