added basic wrapper for sparse matrices, working on tutorial 203

Former-commit-id: 797a95c4fb4a3de6b0a64c53a13e6b31fef7d959

python/203_CurvatureDirections.py

@@ -1,56 +1,57 @@
 import igl
 
-
-
-
-
-  std::string filename = "../shared/fertility.off";
-  if(argc>1)
-  {
-    filename = argv[1];
-  }
-  // Load a mesh in OFF format
-  igl::read_triangle_mesh(filename, V, F);
-
-  // Alternative discrete mean curvature
-  MatrixXd HN;
-  SparseMatrix<double> L,M,Minv;
-  igl::cotmatrix(V,F,L);
-  igl::massmatrix(V,F,igl::MASSMATRIX_TYPE_VORONOI,M);
-  igl::invert_diag(M,Minv);
-  // Laplace-Beltrami of position
-  HN = -Minv*(L*V);
-  // Extract magnitude as mean curvature
-  VectorXd H = HN.rowwise().norm();
-
-  // Compute curvature directions via quadric fitting
-  MatrixXd PD1,PD2;
-  VectorXd PV1,PV2;
-  igl::principal_curvature(V,F,PD1,PD2,PV1,PV2);
-  // mean curvature
-  H = 0.5*(PV1+PV2);
-
-  igl::viewer::Viewer viewer;
-  viewer.data.set_mesh(V, F);
-
-
-  // Compute pseudocolor
-  MatrixXd C;
-  igl::parula(H,true,C);
-  viewer.data.set_colors(C);
-
-  // Average edge length for sizing
-  const double avg = igl::avg_edge_length(V,F);
-
-  // Draw a blue segment parallel to the minimal curvature direction
-  const RowVector3d red(0.8,0.2,0.2),blue(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();
-}
+V = igl.eigen.MatrixXd();
+F = igl.eigen.MatrixXi();
+igl.read_triangle_mesh("../tutorial/shared/fertility.off", V, F);
+
+# Alternative discrete mean curvature
+HN = igl.eigen.MatrixXd();
+L = igl.eigen.SparseMatrixd();
+M = igl.eigen.SparseMatrixd();
+Minv = igl.eigen.SparseMatrixd();
+
+
+igl.cotmatrix(V,F,L);
+igl.massmatrix(V,F,igl.MASSMATRIX_TYPE_VORONOI,M);
+
+igl.invert_diag(M,Minv);
+
+# Laplace-Beltrami of position
+temp = igl.eigen.SparseMatrixd();
+temp = L*V
+HN = -Minv*(L*V);
+
+# Extract magnitude as mean curvature
+#   VectorXd H = HN.rowwise().norm();
+#
+#   // Compute curvature directions via quadric fitting
+#   MatrixXd PD1,PD2;
+#   VectorXd PV1,PV2;
+#   igl::principal_curvature(V,F,PD1,PD2,PV1,PV2);
+#   // mean curvature
+#   H = 0.5*(PV1+PV2);
+#
+#   igl::viewer::Viewer viewer;
+#   viewer.data.set_mesh(V, F);
+#
+#
+#   // Compute pseudocolor
+#   MatrixXd C;
+#   igl::parula(H,true,C);
+#   viewer.data.set_colors(C);
+#
+#   // Average edge length for sizing
+#   const double avg = igl::avg_edge_length(V,F);
+#
+#   // Draw a blue segment parallel to the minimal curvature direction
+#   const RowVector3d red(0.8,0.2,0.2),blue(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

@@ -9,6 +9,10 @@
 #include <igl/per_vertex_normals.h>
 #include <igl/gaussian_curvature.h>
 #include <igl/jet.h>
+#include <igl/read_triangle_mesh.h>
+#include <igl/cotmatrix.h>
+#include <igl/massmatrix.h>
+#include <igl/invert_diag.h>
 
 void python_export_igl(py::module &m)
 {
@@ -19,4 +23,8 @@ void python_export_igl(py::module &m)
 #include "py_igl/py_per_vertex_normals.cpp"
 #include "py_igl/py_gaussian_curvature.cpp"
 #include "py_igl/py_jet.cpp"
+#include "py_igl/py_read_triangle_mesh.cpp"
+#include "py_igl/py_cotmatrix.cpp"
+#include "py_igl/py_massmatrix.cpp"
+#include "py_igl/py_invert_diag.cpp"
 }

python/py_igl/py_cotmatrix.cpp

@@ -0,0 +1,10 @@
+m.def("cotmatrix", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  Eigen::SparseMatrix<double>& L
+)
+{
+  return igl::cotmatrix(V,F,L);
+}, __doc_igl_cotmatrix,
+py::arg("V"), py::arg("F"), py::arg("L"));

python/py_igl/py_invert_diag.cpp

@@ -0,0 +1,9 @@
+m.def("invert_diag", []
+(
+  const Eigen::SparseMatrix<double>& X,
+  Eigen::SparseMatrix<double>& Y
+)
+{
+  return igl::invert_diag(X,Y);
+}, __doc_igl_invert_diag,
+py::arg("X"), py::arg("Y"));

python/py_igl/py_massmatrix.cpp

@@ -0,0 +1,19 @@
+py::enum_<igl::MassMatrixType>(m, "MassMatrixType")
+    .value("MASSMATRIX_TYPE_BARYCENTRIC", igl::MASSMATRIX_TYPE_BARYCENTRIC)
+    .value("MASSMATRIX_TYPE_VORONOI", igl::MASSMATRIX_TYPE_VORONOI)
+    .value("MASSMATRIX_TYPE_FULL", igl::MASSMATRIX_TYPE_FULL)
+    .value("MASSMATRIX_TYPE_DEFAULT", igl::MASSMATRIX_TYPE_DEFAULT)
+    .value("NUM_MASSMATRIX_TYPE", igl::NUM_MASSMATRIX_TYPE)
+    .export_values();
+
+m.def("massmatrix", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  const igl::MassMatrixType type,
+  Eigen::SparseMatrix<double>& M
+)
+{
+  return igl::massmatrix(V,F,type,M);
+}, __doc_igl_massmatrix,
+py::arg("V"), py::arg("F"), py::arg("type"), py::arg("M"));

python/py_igl/py_read_triangle_mesh.cpp

@@ -0,0 +1,36 @@
+m.def("read_triangle_mesh", []
+(
+  const std::string str,
+  std::vector<std::vector<double> >& V,
+  std::vector<std::vector<int> >& F
+)
+{
+  return igl::read_triangle_mesh(str,V,F);
+}, __doc_igl_read_triangle_mesh,
+py::arg("str"), py::arg("V"), py::arg("F"));
+
+m.def("read_triangle_mesh", []
+(
+  const std::string str,
+  Eigen::MatrixXd& V,
+  Eigen::MatrixXi& F
+)
+{
+  return igl::read_triangle_mesh(str,V,F);
+}, __doc_igl_read_triangle_mesh,
+py::arg("str"), py::arg("V"), py::arg("F"));
+
+m.def("read_triangle_mesh", []
+(
+  const std::string str,
+  Eigen::MatrixXd& V,
+  Eigen::MatrixXi& F,
+  std::string & dir,
+  std::string & base,
+  std::string & ext,
+  std::string & name
+)
+{
+  return igl::read_triangle_mesh(str,V,F,dir,base,ext,name);
+}, __doc_igl_read_triangle_mesh,
+py::arg("str"), py::arg("V"), py::arg("F"), py::arg("dir"), py::arg("base"), py::arg("ext"), py::arg("name"));

python/py_vector.cpp

@@ -1,4 +1,5 @@
 #include <Eigen/Dense>
+#include <Eigen/Sparse>
 
 #include "python.h"
 
@@ -339,18 +340,166 @@ py::class_<Type> bind_eigen_2(py::module &m, const char *name,
     return matrix;
 }
 
+/// Creates Python bindings for a dynamic Eigen sparse order-2 tensor (i.e. a matrix)
+template <typename Type>
+py::class_<Type> bind_eigen_sparse_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>())
+        // .def("__init__", [](Type &m, Scalar f) {
+        //     new (&m) Type(1, 1);
+        //     m(0, 0) = f;
+        // })
+        // .def("__init__", [](Type &m, py::buffer b) {
+        //     py::buffer_info info = b.request();
+        //     if (info.format != py::format_descriptor<Scalar>::value())
+        //         throw std::runtime_error("Incompatible buffer format!");
+        //     if (info.ndim == 1) {
+        //         new (&m) Type(info.shape[0], 1);
+        //         memcpy(m.data(), info.ptr, sizeof(Scalar) * m.size());
+        //     } else if (info.ndim == 2) {
+        //         if (info.strides[0] == sizeof(Scalar)) {
+        //             new (&m) Type(info.shape[0], info.shape[1]);
+        //             memcpy(m.data(), info.ptr, sizeof(Scalar) * m.size());
+        //         } else {
+        //             new (&m) Type(info.shape[1], info.shape[0]);
+        //             memcpy(m.data(), info.ptr, sizeof(Scalar) * m.size());
+        //             m.transposeInPlace();
+        //         }
+        //     } else {
+        //         throw std::runtime_error("Incompatible buffer dimension!");
+        //     }
+        // })
+
+        /* 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(); })
+        // .def("setConstant", [](Type &m, Scalar value) { m.setConstant(value); })
+
+        /* Resizing */
+        // .def("resize", [](Type &m, size_t s0, size_t s1) { m.resize(s0, s1); })
+        // .def("resizeLike", [](Type &m, const Type &m2) { m.resizeLike(m2); })
+        // .def("conservativeResize", [](Type &m, size_t s0, size_t s1) { m.conservativeResize(s0, s1); })
+
+        /* Component-wise operations */
+        // .def("cwiseAbs", &Type::cwiseAbs)
+        // .def("cwiseAbs2", &Type::cwiseAbs2)
+        // .def("cwiseSqrt", &Type::cwiseSqrt)
+        // .def("cwiseInverse", &Type::cwiseInverse)
+        // .def("cwiseMin", [](const Type &m1, const Type &m2) -> Type { return m1.cwiseMin(m2); })
+        // .def("cwiseMax", [](const Type &m1, const Type &m2) -> Type { return m1.cwiseMax(m2); })
+        // .def("cwiseMin", [](const Type &m1, Scalar s) -> Type { return m1.cwiseMin(s); })
+        // .def("cwiseMax", [](const Type &m1, Scalar s) -> Type { return m1.cwiseMax(s); })
+        // .def("cwiseProduct", [](const Type &m1, const Type &m2) -> Type { return m1.cwiseProduct(m2); })
+        // .def("cwiseQuotient", [](const Type &m1, const Type &m2) -> Type { return m1.cwiseQuotient(m2); })
+
+        /* 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(py::self * Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic>())
+        .def_cast(py::self / Scalar())
+
+        /* Arithmetic in-place operators */
+        // .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(py::self /= Scalar())
+
+        /* Comparison operators */
+        // .def(py::self == py::self)
+        // .def(py::self != py::self)
+
+        // .def("transposeInPlace", [](Type &m) { m.transposeInPlace(); })
+        // /* Other transformations */
+        // .def("transpose", [](Type &m) -> Type { return m.transpose(); })
+
+        /* Python protocol implementations */
+        .def("__repr__", [](const Type &v) {
+            std::ostringstream oss;
+            oss << v;
+            return oss.str();
+        })
+        // .def("__getitem__", [](const Type &m, std::pair<size_t, size_t> i) {
+        //     if (i.first >= (size_t) m.rows() || i.second >= (size_t) m.cols())
+        //         throw py::index_error();
+        //     return m(i.first, i.second);
+        //  })
+        // .def("__setitem__", [](Type &m, std::pair<size_t, size_t> i, Scalar v) {
+        //     if (i.first >= (size_t) m.rows() || i.second >= (size_t) m.cols())
+        //         throw py::index_error();
+        //     m(i.first, i.second) = v;
+        //  })
+
+        // /* Buffer access for interacting with NumPy */
+        // .def_buffer([](Type &m) -> py::buffer_info {
+        //     return py::buffer_info(
+        //         m.data(),                /* Pointer to buffer */
+        //         sizeof(Scalar),          /* Size of one scalar */
+        //         /* Python struct-style format descriptor */
+        //         py::format_descriptor<Scalar>::value(),
+        //         2,                       /* Number of dimensions */
+        //         { (size_t) m.rows(),     /* Buffer dimensions */
+        //           (size_t) m.cols() },
+        //         { sizeof(Scalar),        /* Strides (in bytes) for each index */
+        //           sizeof(Scalar) * m.rows() }
+        //     );
+        //  })
+
+        /* Static initializers */
+        // .def_static("Zero", [](size_t n, size_t m) { return Type(Type::Zero(n, m)); })
+        // .def_static("Ones", [](size_t n, size_t m) { return Type(Type::Ones(n, m)); })
+        // .def_static("Constant", [](size_t n, size_t m, Scalar value) { return Type(Type::Constant(n, m, value)); })
+        // .def_static("Identity", [](size_t n, size_t m) { return Type(Type::Identity(n, m)); })
+        ;
+    return matrix;
+}
+
+
 void python_export_vector(py::module &m) {
 
   py::module me = m.def_submodule(
     "eigen", "Wrappers for Eigen types");
 
+    /* Bindings for VectorXd */
     bind_eigen_1<Eigen::VectorXd> (me, "VectorXd");
+    py::implicitly_convertible<py::buffer, Eigen::VectorXd>();
+    py::implicitly_convertible<double, Eigen::VectorXd>();
+
+    /* Bindings for VectorXi */
     bind_eigen_1<Eigen::VectorXi> (me, "VectorXi");
+    py::implicitly_convertible<py::buffer, Eigen::VectorXi>();
+    py::implicitly_convertible<double, Eigen::VectorXi>();
+
+    /* Bindings for MatrixXd */
     bind_eigen_2<Eigen::MatrixXd> (me, "MatrixXd");
-    bind_eigen_2<Eigen::MatrixXi> (me, "MatrixXi");
+    py::implicitly_convertible<py::buffer, Eigen::MatrixXd>();
+    py::implicitly_convertible<double, Eigen::MatrixXd>();
 
+    /* Bindings for MatrixXi */
+    bind_eigen_2<Eigen::MatrixXi> (me, "MatrixXi");
+    py::implicitly_convertible<py::buffer, Eigen::MatrixXi>();
+    py::implicitly_convertible<double, Eigen::MatrixXi>();
 
-    /* Bindings for <vector.h> */
+    /* Bindings for Vector3d */
     auto vector3 = bind_eigen_1_3<Eigen::Vector3d>(me, "Vector3d");
     vector3
         .def("norm", [](const Eigen::Vector3d &v) { return v.norm(); })
@@ -366,15 +515,19 @@ void python_export_vector(py::module &m) {
         .def_property("z", [](const Eigen::Vector3d &v) -> double { return v.z(); },
                            [](Eigen::Vector3d &v, double z) { v.z() = z; }, "Z coordinate");
 
-    py::implicitly_convertible<py::buffer, Eigen::VectorXd>();
-    py::implicitly_convertible<py::buffer, Eigen::MatrixXd>();
-    py::implicitly_convertible<py::buffer, Eigen::VectorXi>();
-    py::implicitly_convertible<py::buffer, Eigen::MatrixXi>();
     py::implicitly_convertible<py::buffer, Eigen::Vector3d>();
-
-    py::implicitly_convertible<double, Eigen::VectorXd>();
-    py::implicitly_convertible<double, Eigen::MatrixXd>();
-    py::implicitly_convertible<double, Eigen::VectorXi>();
-    py::implicitly_convertible<double, Eigen::MatrixXi>();
     py::implicitly_convertible<double, Eigen::Vector3d>();
+
+    /* Bindings for SparseMatrix<double> */
+    bind_eigen_sparse_2< Eigen::SparseMatrix<double> > (me, "SparseMatrixd");
+
+    /* Bindings for SparseMatrix<int> */
+    bind_eigen_sparse_2< Eigen::SparseMatrix<int> > (me, "SparseMatrixi");
+
+
+
+
+
+
+
 }