function to determine if edge already exists nearby + test

include/igl/edge_exists_near.cpp

@@ -0,0 +1,102 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 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 "edge_exists_near.h"
+
+template <
+  typename DeriveduE,
+  typename DerivedEMAP,
+  typename uE2EType,
+  typename Index>
+IGL_INLINE bool igl::edge_exists_near(
+  const Eigen::MatrixBase<DeriveduE> & uE,
+  const Eigen::MatrixBase<DerivedEMAP> & EMAP,
+  const std::vector<std::vector< uE2EType> > & uE2E,
+  const Index & a,
+  const Index & b,
+  const Index & uei)
+{
+  std::vector<Index> face_queue;
+  face_queue.reserve(32);
+  std::vector<Index> pushed;
+  // 32 is faster than 8
+  pushed.reserve(32);
+  assert(a!=b);
+  // Not handling case where (a,b) is edge on face incident on uei
+  // since this can't happen for edge-flipping.
+  assert(a!=uE(uei,0));
+  assert(a!=uE(uei,1));
+  assert(b!=uE(uei,0));
+  assert(b!=uE(uei,1));
+  // starting with the (2) faces incident on e, consider all faces
+  // incident on edges containing either a or b.
+  //
+  // face_queue  Queue containing faces incident on exactly one of a/b
+  // Using a vector seems mildly faster
+  assert(EMAP.size()%3 == 0);
+  const Index num_faces = EMAP.size()/3;
+  const Index f1 = uE2E[uei][0]%num_faces;
+  const Index f2 = uE2E[uei][1]%num_faces;
+  // map is faster than unordered_map here, and vector + brute force
+  // is_member check is even faster
+  face_queue.push_back(f1);
+  pushed.push_back(f1);
+  face_queue.push_back(f2);
+  pushed.push_back(f2);
+  while(!face_queue.empty())
+  {
+    const Index f = face_queue.back();
+    face_queue.pop_back();
+    // consider each edge of this face
+    for(int c = 0;c<3;c++)
+    {
+      // Unique edge id
+      const Index uec = EMAP(c*num_faces+f);
+      const Index s = uE(uec,0);
+      const Index d = uE(uec,1);
+      const bool ona = s == a || d == a;
+      const bool onb = s == b || d == b;
+      // Is this the edge we're looking for?
+      if(ona && onb)
+      {
+        return true;
+      }
+      // not incident on either?
+      if(!ona && !onb)
+      {
+        continue;
+      }
+      // loop over all incident half-edges
+      for(const auto & he : uE2E[uec])
+      {
+        // face of this he
+        const Index fhe = he%num_faces;
+        bool already_pushed = false;
+        for(const auto & fp : pushed)
+        {
+          if(fp == fhe)
+          {
+            already_pushed = true;
+            break;
+          }
+        }
+        if(!already_pushed)
+        {
+          pushed.push_back(fhe);
+          face_queue.push_back(fhe);
+        }
+      }
+    }
+  }
+  return false;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template bool igl::edge_exists_near<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, int, int>(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, int const&, int const&, int const&);
+#endif 

include/igl/edge_exists_near.h

@@ -0,0 +1,47 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 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_EDGE_EXISTS_NEAR_H
+#define IGL_EDGE_EXISTS_NEAR_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+namespace igl
+{
+  // Does edge (a,b) exist in the edges of all faces incident on
+  // existing unique edge uei.
+  //
+  // Inputs:
+  //   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
+  //   E  #F*3 by 2 list of half-edges
+  //   a  1st end-point of query edge
+  //   b  2nd end-point of query edge
+  //   uei  index into uE/uE2E of unique edge
+  // Returns true if edge exists near uei.
+  //
+  // See also: unique_edge_map
+  template <
+    typename DeriveduE,
+    typename DerivedEMAP,
+    typename uE2EType,
+    typename Index>
+  IGL_INLINE bool edge_exists_near(
+    const Eigen::MatrixBase<DeriveduE> & uE,
+    const Eigen::MatrixBase<DerivedEMAP> & EMAP,
+    const std::vector<std::vector< uE2EType> > & uE2E,
+    const Index & a,
+    const Index & b,
+    const Index & uei);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "edge_exists_near.cpp"
+#endif
+#endif
+

tests/include/igl/edge_exists_near.cpp

@@ -0,0 +1,26 @@
+#include <test_common.h>
+#include <igl/edge_exists_near.h>
+#include <igl/unique_edge_map.h>
+
+TEST(edge_exists_near,tet)
+{
+  const Eigen::MatrixXi F = (Eigen::MatrixXi(4,3)<<
+     0,1,2,
+     0,2,3,
+     0,3,1,
+     1,3,2).finished();
+  Eigen::MatrixXi E,uE;
+  Eigen::VectorXi EMAP;
+  std::vector<std::vector<int> > uE2E;
+  igl::unique_edge_map(F,E,uE,EMAP,uE2E);
+  for(int uei = 0;uei<uE2E.size();uei++)
+  {
+    for(int i = 0;i<4;i++)
+    {
+      for(int j = 0;j<4;j++)
+      {
+        ASSERT_TRUE(igl::edge_exists_near(uE,EMAP,uE2E,i,j,uei));
+      }
+    }
+  }
+}