turning number in 2D

include/igl/turning_number.cpp

@@ -0,0 +1,36 @@
+#include "turning_number.h"
+#include "PI.h"
+
+
+template <typename DerivedV>
+IGL_INLINE typename DerivedV::Scalar igl::turning_number(
+  const Eigen::MatrixBase<DerivedV> & V)
+{
+  typedef typename DerivedV::Scalar Scalar;
+  constexpr Scalar TWO_PI = 2.0 * igl::PI;
+    
+  const int n = V.rows();
+  Scalar total_angle = 0.0;
+
+  for (int i = 0; i < n; i++)
+  {
+    Eigen::Matrix<Scalar, 1, 2> current = V.row(i).head(2);
+    Eigen::Matrix<Scalar, 1, 2> next = V.row((i + 1) % n).head(2);
+    
+    Eigen::Matrix<Scalar, 1, 2> d1 = next - current;
+    Eigen::Matrix<Scalar, 1, 2> d2 = V.row((i + 2) % n).head(2) - next;
+
+    const Scalar angle = atan2(
+      d1(0)*d2(1) - d1(1)*d2(0), 
+      d1(0)*d2(0) + d1(1)*d2(1));
+    total_angle += angle;
+  }
+
+  return total_angle / TWO_PI;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+template double igl::turning_number<Eigen::Matrix<double, -1, 2, 0, -1, 2> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> > const&);
+template double igl::turning_number<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&);
+#endif

include/igl/turning_number.h

@@ -0,0 +1,22 @@
+#ifndef IGL_TURNING_NUMBER_H
+#define IGL_TURNING_NUMBER_H
+#include <Eigen/Core>
+#include "igl_inline.h"
+
+namespace igl
+{
+  /// Compute the turning number of a closed curve in the plane.
+  ///
+  /// @param[in] V  #V by 2 list of vertex positions
+  /// @return tn  turning number
+  ///
+  template <typename DerivedV>
+  IGL_INLINE typename DerivedV::Scalar turning_number(
+    const Eigen::MatrixBase<DerivedV> & V);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "turning_number.cpp"
+#endif
+
+#endif

tests/include/igl/turning_number.cpp

@@ -0,0 +1,38 @@
+#include <test_common.h>
+#include <igl/turning_number.h>
+#include <igl/matlab_format.h>
+#include <iostream>
+
+TEST_CASE("turning_number: pentagon", "[igl]")
+{
+  Eigen::MatrixXd V(5, 2);
+  V << 0, 0    ,
+       1, 0    ,
+       1, 1    ,
+       0.5, 1.5,
+       0, 1    ;
+  const double t_ccw = igl::turning_number(V);
+  // Reverse rows of V
+  V = V.colwise().reverse().eval();
+  const double t_cw = igl::turning_number(V);
+  REQUIRE( abs(t_ccw - 1) < 1e-16 );
+  REQUIRE( abs(t_cw - -1) < 1e-16 );
+}
+
+TEST_CASE("turning_number: heptagram", "[igl]")
+{
+  Eigen::MatrixXd V(7,2);
+  for (int i = 0; i < 7; i++)
+  {
+      // Multiply by 3 to get the {7/3} skipping pattern
+      double theta = 2.0 * M_PI * (3 * i) / 7;
+      V(i, 0) = std::cos(theta);
+      V(i, 1) = std::sin(theta);
+  }
+  const double t_ccw = igl::turning_number(V);
+  // Reverse rows of V
+  V = V.colwise().reverse().eval();
+  const double t_cw = igl::turning_number(V);
+  REQUIRE( abs(t_ccw - 3) < 1e-16 );
+  REQUIRE( abs(t_cw - -3) < 1e-16 );
+}