gaussian curvature and internal angles

Former-commit-id: 0e9cd73d246fe2e401e56ece44b331cfb49ad874

include/igl/gaussian_curvature.cpp

@@ -0,0 +1,55 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2013 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 "gaussian_curvature.h"
+#include "internal_angles.h"
+#include "PI.h"
+#include <iostream>
+template <typename DerivedV, typename DerivedF, typename DerivedK>
+IGL_INLINE void igl::gaussian_curvature(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedK> & K)
+{
+  using namespace Eigen;
+  using namespace std;
+  // internal corner angles
+  Matrix<
+    typename DerivedV::Scalar,
+    DerivedF::RowsAtCompileTime,
+    DerivedF::ColsAtCompileTime> A;
+  internal_angles(V,F,A);
+  K.resize(V.rows(),1);
+  K.setConstant(V.rows(),1,2.*PI);
+  assert(A.rows() == F.rows());
+  assert(A.cols() == F.cols());
+  assert(K.rows() == V.rows());
+  assert(F.maxCoeff() < V.rows());
+  assert(K.cols() == 1);
+  const int Frows = F.rows();
+  //K_G(x_i) = (2π - ∑θj)
+//#ifndef IGL_GAUSSIAN_CURVATURE_OMP_MIN_VALUE
+//#  define IGL_GAUSSIAN_CURVATURE_OMP_MIN_VALUE 1000
+//#endif
+//#pragma omp parallel for if (Frows>IGL_GAUSSIAN_CURVATURE_OMP_MIN_VALUE)
+  for(int f = 0;f<Frows;f++)
+  {
+    // throw normal at each corner
+    for(int j = 0; j < 3;j++)
+    {
+      // Q: Does this need to be critical?
+      // H: I think so, sadly. Maybe there's a way to use reduction
+//#pragma omp critical
+      K(F(f,j),0) -=  A(f,j);
+    }
+  }
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+template void igl::gaussian_curvature<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+#endif

include/igl/gaussian_curvature.h

@@ -0,0 +1,32 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2013 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_GAUSSIAN_CURVATURE_H
+#define IGL_GAUSSIAN_CURVATURE_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Compute discrete gaussian curvature (angle deficit)
+  // Inputs:
+  //   V  #V by 3 eigen Matrix of mesh vertex 3D positions
+  //   F  #F by 3 eigen Matrix of face (triangle) indices
+  // Output:
+  //   K  #V by 1 eigen Matrix of discrete gaussian curvature values
+  template <typename DerivedV, typename DerivedF, typename DerivedK>
+  IGL_INLINE void gaussian_curvature(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    Eigen::PlainObjectBase<DerivedK> & K);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "gaussian_curvature.cpp"
+#endif
+
+#endif
+

include/igl/internal_angles.cpp

@@ -0,0 +1,48 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2013 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 "internal_angles.h"
+#include "edge_lengths.h"
+
+template <typename DerivedV, typename DerivedF, typename DerivedK>
+IGL_INLINE void igl::internal_angles(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedK> & K)
+{
+  using namespace Eigen;
+  // Edge lengths
+  Matrix<
+    typename DerivedV::Scalar,
+    DerivedF::RowsAtCompileTime,
+    DerivedF::ColsAtCompileTime> L;
+  edge_lengths(V,F,L);
+
+  assert(F.cols() == 3 && "F should contain triangles");
+  return internal_angles(L,K);
+}
+
+template <typename DerivedL, typename DerivedK>
+IGL_INLINE void igl::internal_angles(
+  const Eigen::PlainObjectBase<DerivedL>& L,
+  Eigen::PlainObjectBase<DerivedK> & K)
+{
+  assert(L.cols() == 3 && "Edge-lengths should come from triangles");
+  K.resize(L.rows(),L.cols());
+  for(int d = 0;d<3;d++)
+  {
+    const auto & s1 = L.col(d).array();
+    const auto & s2 = L.col((d+1)%3).array();
+    const auto & s3 = L.col((d+2)%3).array();
+    K.col(d) = ((s3.square() + s2.square() - s1.square())/(2.*s3*s2)).acos();
+  }
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+template void igl::internal_angles<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+#endif

include/igl/internal_angles.h

@@ -0,0 +1,40 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2013 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_INTERNAL_ANGLES_H
+#define IGL_INTERNAL_ANGLES_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Compute internal angles for a triangle mesh
+  // Inputs:
+  //   V  #V by 3 eigen Matrix of mesh vertex 3D positions
+  //   F  #F by 3 eigen Matrix of face (triangle) indices
+  // Output:
+  //   K  #F by 3 eigen Matrix of internal angles
+  //     for triangles, columns correspond to edges [1,2],[2,0],[0,1]
+  template <typename DerivedV, typename DerivedF, typename DerivedK>
+  IGL_INLINE void internal_angles(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    Eigen::PlainObjectBase<DerivedK> & K);
+  // Inputs:
+  //   L  #F by 3 list of edge lengths
+  template <typename DerivedL, typename DerivedK>
+  IGL_INLINE void internal_angles(
+    const Eigen::PlainObjectBase<DerivedL>& L,
+    Eigen::PlainObjectBase<DerivedK> & K);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "internal_angles.cpp"
+#endif
+
+#endif
+
+

include/igl/per_face_normals.cpp

@@ -22,8 +22,8 @@ IGL_INLINE void igl::per_face_normals(
 #pragma omp parallel for
   for(int i = 0; i < Frows;i++)
   {
-    Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v1 = V.row(F(i,1)) - V.row(F(i,0));
-    Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v2 = V.row(F(i,2)) - V.row(F(i,0));
+    const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v1 = V.row(F(i,1)) - V.row(F(i,0));
+    const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v2 = V.row(F(i,2)) - V.row(F(i,0));
     N.row(i) = v1.cross(v2);//.normalized();
     typename DerivedV::Scalar r = N.row(i).norm();
     if(r == 0)

include/igl/per_vertex_normals.cpp

@@ -32,19 +32,26 @@ IGL_INLINE void igl::per_vertex_normals(
   N.setZero(V.rows(),3);
 
   // loop over faces
-  int Frows = F.rows();
-#pragma omp parallel for
+  const int Frows = F.rows();
+//// Minimum number of iterms per openmp thread
+//#ifndef IGL_OMP_MIN_VALUE
+//#  define IGL_OMP_MIN_VALUE 1000
+//#endif
+//#pragma omp parallel for if (Frows>IGL_OMP_MIN_VALUE)
   for(int i = 0; i < Frows;i++)
   {
     // throw normal at each corner
     for(int j = 0; j < 3;j++)
     {
+      // Does this need to be critical?
+//#pragma omp critical
       N.row(F(i,j)) += FN.row(i);
     }
   }
   // normalize each row
   igl::normalize_row_lengths(N,N);
 }
+
 #ifndef IGL_HEADER_ONLY
 // Explicit template specialization
 // generated by autoexplicit.sh

include/igl/per_vertex_normals.h

@@ -26,7 +26,8 @@ namespace igl
     const Eigen::PlainObjectBase<DerivedV>& V,
     const Eigen::PlainObjectBase<DerivedF>& F,
     Eigen::PlainObjectBase<DerivedV> & N);
-
+  // Inputs:
+  //   FN  #F by 3 matrix of face (triangle) normals
   template <typename DerivedV, typename DerivedF>
   IGL_INLINE void per_vertex_normals(
     const Eigen::PlainObjectBase<DerivedV>& V,

include/igl/principal_curvature.cpp

@@ -296,10 +296,6 @@ public:
   
 };
 
-
-
-
-
 class comparer
 {
 public:

include/igl/principal_curvature.h

@@ -47,13 +47,12 @@ namespace igl
 
 template <typename DerivedV, typename DerivedF>
 IGL_INLINE void principal_curvature(
-                                     const Eigen::PlainObjectBase<DerivedV>& V,
-                                     const Eigen::PlainObjectBase<DerivedF>& F,
-                                     Eigen::PlainObjectBase<DerivedV>& PD1,
-                                     Eigen::PlainObjectBase<DerivedV>& PD2,
-                                     unsigned radius = 5,
-                                     bool useKring = true
-                                     );
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedV>& PD1,
+  Eigen::PlainObjectBase<DerivedV>& PD2,
+  unsigned radius = 5,
+  bool useKring = true);
 }
 
 

include/igl/viewer/TODOs.txt

@@ -1,4 +1,5 @@
 - rewrite in libigl style
+- separate various class into their own .h/.cpp pairs
 - remove use of double underscores (http://stackoverflow.com/a/224420/148668)
 - document inputs and outputs to all functions
 - document all member fields