per edge normals

Former-commit-id: 866592823d07391ef040e759bbaaf46daf8971dc

include/igl/per_edge_normals.cpp

@@ -0,0 +1,54 @@
+#include "all_edges.h"
+#include "doublearea.h"
+#include "per_edge_normals.h"
+#include "per_face_normals.h"
+#include "unique_simplices.h"
+#include <vector>
+
+template <
+  typename DerivedV, 
+  typename DerivedF, 
+  typename DerivedN,
+  typename DerivedE,
+  typename DerivedEMAP>
+IGL_INLINE void igl::per_edge_normals(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedN> & N,
+  Eigen::PlainObjectBase<DerivedE> & E,
+  Eigen::PlainObjectBase<DerivedEMAP> & EMAP)
+{
+  using namespace Eigen;
+  using namespace std;
+  assert(F.cols() == 3 && "Faces must be triangles");
+  // number of faces
+  const int m = F.rows();
+  // All occurances of directed edges
+  MatrixXi allE;
+  all_edges(F,allE);
+  // Find unique undirected edges and mapping
+  VectorXi _;
+  unique_simplices(allE,E,_,EMAP);
+  // now sort(allE,2) == E(EMAP,:), that is, if EMAP(i) = j, then E.row(j) is
+  // the undirected edge corresponding to the directed edge allE.row(i).
+  MatrixXd FN;
+  per_face_normals(V,F,FN);
+
+  VectorXd dblA;
+  doublearea(V,F,dblA);
+  N.setConstant(E.rows(),3,0);
+  for(int f = 0;f<m;f++)
+  {
+    for(int c = 0;c<3;c++)
+    {
+      N.row(EMAP(f+c*m)) += dblA(f) * FN.row(f);
+    }
+  }
+  N.rowwise().normalize();
+  
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instanciation
+template void igl::per_edge_normals<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::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/per_edge_normals.h

@@ -0,0 +1,42 @@
+// 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_PER_EDGE_NORMALS_H
+#define IGL_PER_EDGE_NORMALS_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Compute face normals via vertex position list, face list
+  // Inputs:
+  //   V  #V by 3 eigen Matrix of mesh vertex 3D positions
+  //   F  #F by 3 eigen Matrix of face (triangle) indices
+  // Output:
+  //   N  #2 by 3 matrix of mesh edge 3D normals per row
+  //   E  #E by 2 matrix of edge indices per row
+  //   EMAP  #E by 1 matrix of indices from all edges to E
+  //
+  template <
+    typename DerivedV, 
+    typename DerivedF, 
+    typename DerivedN,
+    typename DerivedE,
+    typename DerivedEMAP>
+  IGL_INLINE void per_edge_normals(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    Eigen::PlainObjectBase<DerivedN> & N,
+    Eigen::PlainObjectBase<DerivedE> & E,
+    Eigen::PlainObjectBase<DerivedEMAP> & EMAP);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "per_edge_normals.cpp"
+#endif
+
+#endif
+

include/igl/unique_simplices.cpp

@@ -9,9 +9,16 @@
 #include "sort.h"
 #include "unique.h"
 
+template <
+  typename DerivedF,
+  typename DerivedFF,
+  typename DerivedIA,
+  typename DerivedIC>
 IGL_INLINE void igl::unique_simplices(
-  const Eigen::MatrixXi & F,
-  Eigen::MatrixXi & FF)
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedFF>& FF,
+  Eigen::PlainObjectBase<DerivedIA>& IA,
+  Eigen::PlainObjectBase<DerivedIC>& IC)
 {
   using namespace Eigen;
   using namespace igl;
@@ -19,7 +26,6 @@ IGL_INLINE void igl::unique_simplices(
   MatrixXi sortF, unusedI;
   igl::sort(F,2,1,sortF,unusedI);
   // Find unique faces
-  VectorXi IA,IC;
   MatrixXi C;
   igl::unique_rows(sortF,C,IA,IC);
   FF.resize(IA.size(),F.cols());
@@ -29,3 +35,19 @@ IGL_INLINE void igl::unique_simplices(
     FF.row(i) = F.row(IA(i));
   }
 }
+
+template <
+  typename DerivedF,
+  typename DerivedFF>
+IGL_INLINE void igl::unique_simplices(
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedFF>& FF)
+{
+  Eigen::VectorXi IA,IC;
+  return unique_simplices(F,FF,IA,IC);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instanciations
+template void igl::unique_simplices<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/unique_simplices.h

@@ -11,15 +11,30 @@
 #include <Eigen/Dense>
 namespace igl
 {
-  // Find *combinatorially* unique simplices in F
+  // Find *combinatorially* unique simplices in F.  **Order independent**
   //
   // Inputs:
   //   F  #F by simplex-size list of simplices
   // Outputs:
   //   FF  #FF by simplex-size list of unique simplices in F
+  //   IA  #FF index vector so that FF == sort(F(IA,:),2);
+  //   IC  #F index vector so that sort(F,2) == FF(IC,:);
+  template <
+    typename DerivedF,
+    typename DerivedFF,
+    typename DerivedIA,
+    typename DerivedIC>
   IGL_INLINE void unique_simplices(
-    const Eigen::MatrixXi & F,
-    Eigen::MatrixXi & FF);
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    Eigen::PlainObjectBase<DerivedFF>& FF,
+    Eigen::PlainObjectBase<DerivedIA>& IA,
+    Eigen::PlainObjectBase<DerivedIC>& IC);
+  template <
+    typename DerivedF,
+    typename DerivedFF>
+  IGL_INLINE void unique_simplices(
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    Eigen::PlainObjectBase<DerivedFF>& FF);
 }
 
 #ifndef IGL_STATIC_LIBRARY

tutorial/205_Laplacian/main.cpp

@@ -1,5 +1,6 @@
 #include <igl/readOFF.h>
 #include <igl/readDMAT.h>
+#include <igl/barycenter.h>
 #include <igl/cotmatrix.h>
 #include <igl/massmatrix.h>
 #include <igl/grad.h>
@@ -59,8 +60,20 @@ int main(int argc, char *argv[])
         Eigen::SimplicialLLT<Eigen::SparseMatrix<double > > solver(S);
         assert(solver.info() == Eigen::Success);
         U = solver.solve(M*U).eval();
+        // Compute centroid and subtract (also important for numerics)
+        VectorXd dblA;
+        igl::doublearea(U,F,dblA);
+        double area = 0.5*dblA.sum();
+        MatrixXd BC;
+        igl::barycenter(U,F,BC);
+        RowVector3d centroid(0,0,0);
+        for(int i = 0;i<BC.rows();i++)
+        {
+          centroid += 0.5*dblA(i)/area*BC.row(i);
+        }
+        U.rowwise() -= centroid;
         // Normalize to unit surface area (important for numerics)
-        U.array() /= sqrt(M.diagonal().array().sum());
+        U.array() /= sqrt(area);
         break;
       }
       default: