added 205 laplacian

update pybind
added example of usage in matlab


Former-commit-id: b0f890d25be2027f6701157c3b42f38c20643014

.gitignore

@@ -83,3 +83,4 @@ Untitled.ipynb
 python/.ipynb_checkpoints
 python/py_igl/todo
 python/__pycache__
+iglhelpers.pyc

python/102_DrawMesh.m

@@ -0,0 +1,12 @@
+% Load a mesh in OFF format
+V = py.igl.eigen.MatrixXd();
+F = py.igl.eigen.MatrixXi();
+py.igl.readOFF('../tutorial/shared/beetle.off', V, F);
+
+V
+
+% Plot the mesh
+viewer = py.igl.viewer.Viewer();
+viewer.data.set_mesh(V, F);
+viewer.launch();
+

python/205_Laplacian.py

@@ -0,0 +1,109 @@
+from __future__ import print_function
+import igl
+from iglhelpers import *
+import math
+
+global V
+global U
+global F
+global L
+
+V = igl.eigen.MatrixXd()
+U = igl.eigen.MatrixXd()
+F = igl.eigen.MatrixXi()
+
+L = igl.eigen.SparseMatrixd()
+viewer = igl.viewer.Viewer()
+
+# Load a mesh in OFF format
+igl.readOFF("../tutorial/shared/cow.off", V, F)
+
+# Compute Laplace-Beltrami operator: #V by #V
+igl.cotmatrix(V,F,L)
+
+# Alternative construction of same Laplacian
+G = igl.eigen.SparseMatrixd()
+K = igl.eigen.SparseMatrixd()
+
+# Gradient/Divergence
+igl.grad(V,F,G);
+
+# Diagonal per-triangle "mass matrix"
+dblA = igl.eigen.MatrixXd()
+igl.doublearea(V,F,dblA)
+
+# Place areas along diagonal #dim times
+
+T = (dblA.replicate(3,1)*0.5).asDiagonal() * 1
+
+# Laplacian K built as discrete divergence of gradient or equivalently
+# discrete Dirichelet energy Hessian
+
+temp = -G.transpose()
+K = -G.transpose() * T * G
+print("|K-L|: ",(K-L).norm())
+
+def key_pressed(viewer, key, modifier):
+    global V
+    global U
+    global F
+    global L
+
+    if key == ord('r') or key == ord('R'):
+        U = V;
+    elif key == ord(' '):
+
+        # Recompute just mass matrix on each step
+        M = igl.eigen.SparseMatrixd()
+
+        igl.massmatrix(U,F,igl.MASSMATRIX_TYPE_BARYCENTRIC,M);
+
+        # Solve (M-delta*L) U = M*U
+        S = (M - 0.001*L)
+
+        solver = igl.eigen.SimplicialLLTsparse(S)
+
+        U = solver.solve(M*U)
+
+        # Compute centroid and subtract (also important for numerics)
+        dblA = igl.eigen.MatrixXd()
+        igl.doublearea(U,F,dblA)
+
+        print(dblA.sum())
+
+        area = 0.5*dblA.sum()
+        BC = igl.eigen.MatrixXd()
+        igl.barycenter(U,F,BC)
+        centroid = p2e(np.array([[0.0,0.0,0.0]]))
+
+        for i in range(0,BC.rows()):
+            centroid += 0.5*dblA[i,0]/area*BC.row(i)
+
+        U -= centroid.replicate(U.rows(),1)
+
+        # Normalize to unit surface area (important for numerics)
+        U = U / math.sqrt(area)
+    else:
+        return False
+
+    # Send new positions, update normals, recenter
+    viewer.data.set_vertices(U)
+    viewer.data.compute_normals()
+    viewer.core.align_camera_center(U,F)
+    return True
+
+# Use original normals as pseudo-colors
+N = igl.eigen.MatrixXd()
+igl.per_vertex_normals(V,F,N)
+C = N.rowwiseNormalized()*0.5+0.5;
+
+# Initialize smoothing with base mesh
+U = V
+viewer.data.set_mesh(U, F)
+viewer.data.set_colors(C)
+viewer.callback_key_pressed = key_pressed
+
+print("Press [space] to smooth.")
+print("Press [r] to reset.")
+
+viewer.launch()

python/CurvatureDirections.m

@@ -0,0 +1,59 @@
+V = py.igl.eigen.MatrixXd();
+F = py.igl.eigen.MatrixXi();
+py.igl.read_triangle_mesh('../tutorial/shared/fertility.off', V, F);
+
+% Alternative discrete mean curvature
+HN = py.igl.eigen.MatrixXd();
+L = py.igl.eigen.SparseMatrixd();
+M = py.igl.eigen.SparseMatrixd();
+Minv = py.igl.eigen.SparseMatrixd();
+
+
+py.igl.cotmatrix(V,F,L)
+py.igl.massmatrix(V,F,py.igl.MASSMATRIX_TYPE_VORONOI,M)
+
+py.igl.invert_diag(M,Minv)
+
+% Laplace-Beltrami of position
+HN = -Minv*(L*V)
+
+% Extract magnitude as mean curvature
+H = HN.rowwiseNorm()
+
+% Compute curvature directions via quadric fitting
+PD1 = py.igl.eigen.MatrixXd()
+PD2 = py.igl.eigen.MatrixXd()
+
+PV1 = py.igl.eigen.MatrixXd()
+PV2 = py.igl.eigen.MatrixXd()
+
+py.igl.principal_curvature(V,F,PD1,PD2,PV1,PV2)
+
+% Mean curvature
+H = 0.5*(PV1+PV2)
+
+viewer = py.igl.viewer.Viewer()
+viewer.data.set_mesh(V, F)
+
+% Compute pseudocolor
+C = py.igl.eigen.MatrixXd()
+py.igl.parula(H,true,C)
+
+viewer.data.set_colors(C)
+
+% Average edge length for sizing
+avg = py.igl.avg_edge_length(V,F)
+
+% Draw a blue segment parallel to the minimal curvature direction
+red = py.iglhelpers.p2e(py.numpy.array([[0.8,0.2,0.2]]))
+blue = py.iglhelpers.p2e(py.numpy.array([[0.2,0.2,0.8]]))
+
+viewer.data.add_edges(V + PD1*avg, V - PD1*avg, blue)
+
+% Draw a red segment parallel to the maximal curvature direction
+viewer.data.add_edges(V + PD2*avg, V - PD2*avg, red)
+
+% Hide wireframe
+viewer.core.show_lines = false
+
+viewer.launch()

python/py_igl.cpp

@@ -19,6 +19,7 @@
 #include <igl/grad.h>
 #include <igl/avg_edge_length.h>
 #include <igl/barycenter.h>
+#include <igl/doublearea.h>
 
 void python_export_igl(py::module &m)
 {
@@ -39,5 +40,6 @@ void python_export_igl(py::module &m)
 #include "py_igl/py_grad.cpp"
 #include "py_igl/py_avg_edge_length.cpp"
 #include "py_igl/py_barycenter.cpp"
+#include "py_igl/py_doublearea.cpp"
 
 }

python/py_igl/py_doublearea.cpp

@@ -0,0 +1,54 @@
+m.def("doublearea", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  Eigen::MatrixXd& dblA
+)
+{
+  return igl::doublearea(V,F,dblA);
+}, __doc_igl_doublearea,
+py::arg("V"), py::arg("F"), py::arg("dblA"));
+
+m.def("doublearea", []
+(
+  const Eigen::MatrixXd& A,
+  const Eigen::MatrixXd& B,
+  const Eigen::MatrixXd& C,
+  Eigen::MatrixXd& D
+)
+{
+  return igl::doublearea(A,B,C,D);
+}, __doc_igl_doublearea,
+py::arg("A"), py::arg("B"), py::arg("C"), py::arg("D"));
+
+m.def("doublearea_single", []
+(
+  const Eigen::MatrixXd& A,
+  const Eigen::MatrixXd& B,
+  const Eigen::MatrixXd& C
+)
+{
+  return igl::doublearea_single(A,B,C);
+}, __doc_igl_doublearea,
+py::arg("A"), py::arg("B"), py::arg("C"));
+
+m.def("doublearea", []
+(
+  const Eigen::MatrixXd& l,
+  Eigen::MatrixXd& dblA
+)
+{
+  return igl::doublearea(l,dblA);
+}, __doc_igl_doublearea,
+py::arg("l"), py::arg("dblA"));
+
+m.def("doublearea_quad", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  Eigen::MatrixXd& dblA
+)
+{
+  return igl::doublearea_quad(V,F,dblA);
+}, __doc_igl_doublearea,
+py::arg("V"), py::arg("F"), py::arg("dblA"));

python/py_vector.cpp

@@ -74,6 +74,12 @@ py::class_<Type> bind_eigen_2(py::module &m, const char *name,
 
         .def("mean", [](const Type &m) {return m.mean();})
 
+        .def("sum", [](const Type &m) {return m.sum();})
+        .def("prod", [](const Type &m) {return m.prod();})
+        .def("trace", [](const Type &m) {return m.trace();})
+        .def("norm", [](const Type &m) {return m.norm();})
+        .def("squaredNorm", [](const Type &m) {return m.squaredNorm();})
+
         /* Component-wise operations */
         .def("cwiseAbs", &Type::cwiseAbs)
         .def("cwiseAbs2", &Type::cwiseAbs2)
@@ -91,6 +97,7 @@ py::class_<Type> bind_eigen_2(py::module &m, const char *name,
         .def("rowwiseProd", [](const Type &m) {return Type(m.rowwise().prod());} )
         .def("rowwiseMean", [](const Type &m) {return Type(m.rowwise().mean());} )
         .def("rowwiseNorm", [](const Type &m) {return Type(m.rowwise().norm());} )
+        .def("rowwiseNormalized", [](const Type &m) {return Type(m.rowwise().normalized());} )
         .def("rowwiseMinCoeff", [](const Type &m) {return Type(m.rowwise().minCoeff());} )
         .def("rowwiseMaxCoeff", [](const Type &m) {return Type(m.rowwise().maxCoeff());} )
 
@@ -101,6 +108,9 @@ py::class_<Type> bind_eigen_2(py::module &m, const char *name,
         .def("colwiseMinCoeff", [](const Type &m) {return Type(m.colwise().minCoeff());} )
         .def("colwiseMaxCoeff", [](const Type &m) {return Type(m.colwise().maxCoeff());} )
 
+        .def("replicate", [](const Type &m, const int& r, const int& c) {return Type(m.replicate(r,c));} )
+        .def("asDiagonal", [](const Type &m) {return Eigen::DiagonalMatrix<Scalar,Eigen::Dynamic,Eigen::Dynamic>(m.asDiagonal());} )
+
         /* Arithmetic operators (def_cast forcefully casts the result back to a
            Type to avoid type issues with Eigen's crazy expression templates) */
         .def_cast(-py::self)
@@ -258,9 +268,11 @@ py::class_<Type> bind_eigen_sparse_2(py::module &m, const char *name,
         .def("rows", [](const Type &m) { return m.rows(); })
 
         /* Initialization */
-        // .def("setZero", [](Type &m) { m.setZero(); })
-        // .def("setIdentity", [](Type &m) { m.setIdentity(); })
-        // .def("setConstant", [](Type &m, Scalar value) { m.setConstant(value); })
+        .def("setZero", [](Type &m) { m.setZero(); })
+        .def("setIdentity", [](Type &m) { m.setIdentity(); })
+
+        .def("transpose", [](Type &m) { return Type(m.transpose()); })
+        .def("norm", [](Type &m) { return m.norm(); })
 
         /* Resizing */
         // .def("resize", [](Type &m, size_t s0, size_t s1) { m.resize(s0, s1); })
@@ -286,7 +298,19 @@ py::class_<Type> bind_eigen_sparse_2(py::module &m, const char *name,
         .def_cast(py::self - py::self)
         .def_cast(py::self * py::self)
         .def_cast(py::self * Scalar())
+        .def_cast(Scalar() * py::self)
         // Special case, sparse * dense produces a dense matrix
+
+        // .def("__mul__", []
+        // (const Type &a, const Scalar& b)
+        // {
+        //   return Type(a * b);
+        // })
+        // .def("__rmul__", [](const Type& a, const Scalar& b)
+        // {
+        //   return Type(b * a);
+        // })
+
         .def("__mul__", []
         (const Type &a, const Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic>& b)
         {
@@ -297,6 +321,16 @@ py::class_<Type> bind_eigen_sparse_2(py::module &m, const char *name,
           return Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic>(b * a);
         })
 
+        .def("__mul__", []
+        (const Type &a, const Eigen::DiagonalMatrix<Scalar,Eigen::Dynamic,Eigen::Dynamic>& b)
+        {
+          return Type(a * b);
+        })
+        .def("__rmul__", [](const Type& a, const Eigen::DiagonalMatrix<Scalar,Eigen::Dynamic,Eigen::Dynamic>& b)
+        {
+          return Type(b * a);
+        })
+
         //.def(py::self * Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic>())
 //        .def_cast(py::self / Scalar())
 
@@ -384,6 +418,82 @@ py::class_<Type> bind_eigen_sparse_2(py::module &m, const char *name,
     return matrix;
 }
 
+/// Creates Python bindings for a diagonal Eigen sparse order-2 tensor (i.e. a matrix)
+template <typename Type>
+py::class_<Type> bind_eigen_diagonal_2(py::module &m, const char *name,
+                                py::object parent = py::object()) {
+    typedef typename Type::Scalar Scalar;
+
+    /* Many Eigen functions are templated and can't easily be referenced using
+       a function pointer, thus a big portion of the binding code below
+       instantiates Eigen code using small anonymous wrapper functions */
+    py::class_<Type> matrix(m, name, parent);
+
+    matrix
+        /* Constructors */
+        .def(py::init<>())
+        //.def(py::init<size_t, size_t>())
+
+        /* Size query functions */
+        .def("size", [](const Type &m) { return m.size(); })
+        .def("cols", [](const Type &m) { return m.cols(); })
+        .def("rows", [](const Type &m) { return m.rows(); })
+
+        /* Initialization */
+        .def("setZero", [](Type &m) { m.setZero(); })
+        .def("setIdentity", [](Type &m) { m.setIdentity(); })
+
+        /* Arithmetic operators (def_cast forcefully casts the result back to a
+           Type to avoid type issues with Eigen's crazy expression templates) */
+        // .def_cast(-py::self)
+        // .def_cast(py::self + py::self)
+        // .def_cast(py::self - py::self)
+        // .def_cast(py::self * py::self)
+        .def_cast(py::self * Scalar())
+        .def_cast(Scalar() * py::self)
+
+        // // Special case, sparse * dense produces a dense matrix
+        // .def("__mul__", []
+        // (const Type &a, const Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic>& b)
+        // {
+        //   return Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic>(a * b);
+        // })
+        // .def("__rmul__", [](const Type& a, const Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic>& b)
+        // {
+        //   return Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic>(b * a);
+        // })
+
+        .def("__mul__", []
+        (const Type &a, const Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic>& b)
+        {
+          return Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic>(a * b);
+        })
+        .def("__rmul__", [](const Type& a, const Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic>& b)
+        {
+          return Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic>(b * a);
+        })
+
+        .def("__mul__", []
+        (const Type &a, const Eigen::SparseMatrix<Scalar>& b)
+        {
+          return Eigen::SparseMatrix<Scalar>(a * b);
+        })
+        .def("__rmul__", [](const Type& a, const Eigen::SparseMatrix<Scalar>& b)
+        {
+          return Eigen::SparseMatrix<Scalar>(b * a);
+        })
+
+        /* Python protocol implementations */
+        .def("__repr__", [](const Type &/*v*/) {
+            std::ostringstream oss;
+            oss << "Eigen is not able to plot Diagonal Matrices";
+            return oss.str();
+        })
+
+        ;
+    return matrix;
+}
+
 
 void python_export_vector(py::module &m) {
 
@@ -440,7 +550,29 @@ void python_export_vector(py::module &m) {
     /* Bindings for SparseMatrix<int> */
     bind_eigen_sparse_2< Eigen::SparseMatrix<int> > (me, "SparseMatrixi");
 
+    /* Bindings for DiagonalMatrix<double> */
+    bind_eigen_diagonal_2< Eigen::DiagonalMatrix<double,Eigen::Dynamic,Eigen::Dynamic> > (me, "DiagonalMatrixd");
+
+    /* Bindings for DiagonalMatrix<int> */
+    bind_eigen_diagonal_2< Eigen::DiagonalMatrix<int,Eigen::Dynamic,Eigen::Dynamic> > (me, "DiagonalMatrixi");
+
+    /* Bindings for SimplicialLLT*/
+
+    py::class_<Eigen::SimplicialLLT<Eigen::SparseMatrix<double > >> simpliciallltsparse(me, "SimplicialLLTsparse");
+
+    simpliciallltsparse
+    .def(py::init<>())
+    .def(py::init<Eigen::SparseMatrix<double>>())
+    .def("info",[](const Eigen::SimplicialLLT<Eigen::SparseMatrix<double > >& s)
+    {
+       if (s.info() == Eigen::Success)
+          return "Success";
+       else
+          return "Numerical Issue";
+    })
+    .def("solve",[](const Eigen::SimplicialLLT<Eigen::SparseMatrix<double > >& s, const Eigen::MatrixXd& rhs) { return Eigen::MatrixXd(s.solve(rhs)); })
 
+    ;