diff --git a/.gitignore b/.gitignore
index e3cb392a1..a98275229 100644
--- a/.gitignore
+++ b/.gitignore
@@ -75,3 +75,8 @@ optional/build
 lib
 tutorial/XXX_test/CMakeLists.txt
 tutorial/XXX_test/main.cpp
+python/build
+python/.idea
+untitled
+Untitled.ipynb
+python/.ipynb_checkpoints
diff --git a/python/101_FileIO.py b/python/101_FileIO.py
new file mode 100644
index 000000000..277a46b9d
--- /dev/null
+++ b/python/101_FileIO.py
@@ -0,0 +1,15 @@
+import sys, os
+import numpy as np
+import igl
+
+# Load a mesh in OFF format
+V = igl.eigen.MatrixXd()
+F = igl.eigen.MatrixXi()
+igl.readOFF("../tutorial/shared/cube.off", V, F)
+
+# Print the vertices and faces matrices
+print("Vertices: \n", V, sep='')
+print("Faces: \n", F, sep='')
+
+# Save the mesh in OBJ format
+igl.writeOBJ("cube.obj",V,F)
diff --git a/python/102_DrawMesh.py b/python/102_DrawMesh.py
new file mode 100755
index 000000000..e0b363e0d
--- /dev/null
+++ b/python/102_DrawMesh.py
@@ -0,0 +1,22 @@
+import sys, os
+import numpy as np
+import igl
+
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+import matplotlib.tri as mtri
+
+# Load a mesh in OFF format
+V = igl.eigen.MatrixXd()
+F = igl.eigen.MatrixXi()
+igl.readOFF("../tutorial/shared/beetle.off", V, F)
+
+# Convert the mesh to numpy matrices (without copying it)
+Vn = np.array(V, copy=False)
+Fn = np.array(F, copy=False)
+
+# Plot using matplotlib
+fig = plt.figure()
+ax = fig.add_subplot(1, 1, 1, projection='3d')
+ax.plot_trisurf(Vn[:,0], Vn[:,1], Vn[:,2], triangles=Fn, cmap=plt.cm.Spectral)
+plt.show()
diff --git a/python/201_Normals.py b/python/201_Normals.py
new file mode 100755
index 000000000..ce673bdad
--- /dev/null
+++ b/python/201_Normals.py
@@ -0,0 +1,24 @@
+
+
+
+# Load a mesh in OFF format
+V = igl.eigen.MatrixXd()
+F = igl.eigen.MatrixXi()
+igl.readOFF("../shared/fandisk.off", V, F);
+
+# Compute per-face normals
+N_faces = igl.eigen.MatrixXd()
+igl::per_face_normals(V,F,N_faces);
+print("igl::per_face_normals: \n", N_faces, sep='')
+
+# Compute per-vertex normals
+N_vertices = igl.eigen.MatrixXd()
+igl::per_vertex_normals(V,F,N_vertices);
+print("igl::per_face_normals: \n", N_vertices, sep='')
+
+# Compute per-corner normals, |dihedral angle| > 20 degrees --> crease
+N_corners = igl.eigen.MatrixXd()
+igl::per_corner_normals(V,F,20,N_corners);
+print("igl::per_face_normals: \n", N_corners, sep='')
+
+}
diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt
new file mode 100644
index 000000000..6341a33e4
--- /dev/null
+++ b/python/CMakeLists.txt
@@ -0,0 +1,118 @@
+# CMakeLists.txt -- Build system for the pybind11 examples
+#
+# Copyright (c) 2015 Wenzel Jakob <wenzel@inf.ethz.ch>
+#
+# All rights reserved. Use of this source code is governed by a
+# BSD-style license that can be found in the LICENSE file.
+
+cmake_minimum_required(VERSION 2.8)
+
+project(pybind)
+
+if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
+  message(STATUS "Setting build type to 'MinSizeRel' as none was specified.")
+  set(CMAKE_BUILD_TYPE MinSizeRel CACHE STRING "Choose the type of build." FORCE)
+  set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS "Debug" "Release"
+    "MinSizeRel" "RelWithDebInfo")
+endif()
+
+if (true)
+SET(PYTHON_LIBRARIES "/usr/local/Cellar/python3/3.4.3_2/Frameworks/Python.framework/Versions/3.4/lib/libpython3.4m.dylib")
+SET(PYTHON_INCLUDE_DIR "/usr/local/Cellar/python3/3.4.3_2/Frameworks/Python.framework/Versions/3.4/include/python3.4m")
+
+#set(Python_ADDITIONAL_VERSIONS 3.4 3.5 3.6)
+find_package(PythonLibs REQUIRED)
+find_package(PythonInterp REQUIRED)
+
+else()
+
+SET(PYTHON_LIBRARIES "/usr/local/Cellar/python3/3.4.3_2/Frameworks/Python.framework/Versions/3.4/lib/libpython3.4m.dylib")
+SET(PYTHON_INCLUDE_DIR "/usr/local/Cellar/python3/3.4.3_2/Frameworks/Python.framework/Versions/3.4/include/python3.4m")
+
+set(Python_ADDITIONAL_VERSIONS 3.4 3.5 3.6)
+find_package(PythonLibs 3 REQUIRED)
+find_package(PythonInterp 3 REQUIRED)
+
+endif()
+
+string(TOUPPER "${CMAKE_BUILD_TYPE}" U_CMAKE_BUILD_TYPE)
+if (UNIX)
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11")
+  if (NOT ${U_CMAKE_BUILD_TYPE} MATCHES DEBUG)
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fvisibility=hidden -flto")
+  endif()
+endif()
+
+# Compile with compiler warnings turned on
+if(MSVC)
+  if(CMAKE_CXX_FLAGS MATCHES "/W[0-4]")
+    string(REGEX REPLACE "/W[0-4]" "/W4" CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS}")
+  else()
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /W4")
+  endif()
+else()
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra")
+endif()
+
+include_directories(${PYTHON_INCLUDE_DIR} include)
+
+## include pybing
+include_directories(${PROJECT_SOURCE_DIR}/../external/pybind11/include)
+
+## include eigen
+include_directories(${PROJECT_SOURCE_DIR}/../external/nanogui/ext/eigen)
+
+## include libigl
+include_directories(${PROJECT_SOURCE_DIR}/../include)
+
+add_library(igl SHARED
+  python.cpp
+  py_vector.cpp
+  py_igl.cpp
+  py_doc.cpp
+)
+
+set_target_properties(igl PROPERTIES PREFIX "")
+set_target_properties(igl PROPERTIES LIBRARY_OUTPUT_DIRECTORY ${PROJECT_SOURCE_DIR})
+
+if (WIN32)
+  if (MSVC)
+    # Enforce size-based optimization and link time code generation on MSVC (~30% smaller binaries in experiments)
+    set_target_properties(igl PROPERTIES COMPILE_FLAGS "/Os /GL")
+    set_target_properties(igl PROPERTIES LINK_FLAGS "/LTCG")
+  endif()
+
+  # .PYD file extension on Windows
+  set_target_properties(igl PROPERTIES SUFFIX ".pyd")
+
+  # Link against the Python shared library
+  target_link_libraries(igl ${PYTHON_LIBRARY})
+elseif (UNIX)
+  # It's quite common to have multiple copies of the same Python version
+  # installed on one's system. E.g.: one copy from the OS and another copy
+  # that's statically linked into an application like Blender or Maya.
+  # If we link our plugin library against the OS Python here and import it
+  # into Blender or Maya later on, this will cause segfaults when multiple
+  # conflicting Python instances are active at the same time.
+
+  # Windows does not seem to be affected by this issue. The solution for Linux
+  # and Mac OS is simple: we just don't link against the Python library. The
+  # resulting shared library will have missing symbols, but that's perfectly
+  # fine -- they will be resolved at import time.
+
+  # .SO file extension on Linux/Mac OS
+  set_target_properties(igl PROPERTIES SUFFIX ".so")
+
+  # Strip unnecessary sections of the binary on Linux/Mac OS
+  if(APPLE)
+    set_target_properties(igl PROPERTIES MACOSX_RPATH ".")
+    set_target_properties(igl PROPERTIES LINK_FLAGS "-undefined dynamic_lookup -dead_strip")
+    if (NOT ${U_CMAKE_BUILD_TYPE} MATCHES DEBUG)
+      add_custom_command(TARGET igl POST_BUILD COMMAND strip -u -r ${CMAKE_CURRENT_BINARY_DIR}/igl.so)
+    endif()
+  else()
+    if (NOT ${U_CMAKE_BUILD_TYPE} MATCHES DEBUG)
+      add_custom_command(TARGET igl POST_BUILD COMMAND strip ${CMAKE_CURRENT_BINARY_DIR}/igl.so)
+    endif()
+  endif()
+endif()
diff --git a/python/py_igl.cpp b/python/py_igl.cpp
new file mode 100644
index 000000000..06a9d0eb1
--- /dev/null
+++ b/python/py_igl.cpp
@@ -0,0 +1,98 @@
+#include <Eigen/Dense>
+
+#include "python.h"
+#include <igl/readOFF.h>
+#include <igl/writeOBJ.h>
+#include <igl/per_face_normals.h>
+
+void python_export_igl(py::module &m) {
+
+// readOFF.h
+
+  m.def("readOFF", []
+  (
+    const std::string str,
+    Eigen::MatrixXd& V,
+    Eigen::MatrixXi& F
+  )
+  {
+    return igl::readOFF(str,V,F);
+  }, __doc_readOFF,
+  py::arg("str"), py::arg("V"), py::arg("F"));
+
+  m.def("readOFF", []
+  (
+    const std::string str,
+    Eigen::MatrixXd& V,
+    Eigen::MatrixXi& F,
+    Eigen::MatrixXd& N
+  )
+  {
+    return igl::readOFF(str,V,F,N);
+  }, __doc_readOFF,
+  py::arg("str"), py::arg("V"), py::arg("F"), py::arg("N"));
+
+// writeOBJ.h
+m.def("writeOBJ", []
+(
+  const std::string str,
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  const Eigen::MatrixXd& CN,
+  const Eigen::MatrixXi& FN,
+  const Eigen::MatrixXd& TC,
+  const Eigen::MatrixXi& FTC
+)
+{
+  return igl::writeOBJ(str,V,F,CN,FN,TC,FTC);
+}, __doc_writeOBJ,
+py::arg("str"), py::arg("V"), py::arg("F"), py::arg("CN"), py::arg("FN"), py::arg("TC"), py::arg("FTC"));
+
+m.def("writeOBJ", []
+(
+  const std::string str,
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F
+)
+{
+  return igl::writeOBJ(str,V,F);
+}, __doc_writeOBJ,
+py::arg("str"), py::arg("V"), py::arg("F"));
+
+// per_face_normals
+
+m.def("per_face_normals", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  const Eigen::VectorXd& Z,
+  Eigen::MatrixXd& N
+)
+{
+  return igl::per_face_normals(V,F,Z,N);
+}, __doc_per_face_normals,
+py::arg("V"), py::arg("F"), py::arg("Z"), py::arg("N"));
+
+m.def("per_face_normals", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  Eigen::MatrixXd& N
+)
+{
+  return igl::per_face_normals(V,F,N);
+}, __doc_per_face_normals,
+py::arg("V"), py::arg("F"), py::arg("N"));
+
+m.def("per_face_normals_stable", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  Eigen::MatrixXd& N
+)
+{
+  return igl::per_face_normals_stable(V,F,N);
+}, __doc_per_face_normals,
+py::arg("V"), py::arg("F"), py::arg("N"));
+
+}
diff --git a/python/py_vector.cpp b/python/py_vector.cpp
new file mode 100644
index 000000000..f5bbafc2f
--- /dev/null
+++ b/python/py_vector.cpp
@@ -0,0 +1,380 @@
+#include <Eigen/Dense>
+
+#include "python.h"
+
+template <typename Type> void init_fixed_from_buffer_3(Type &v, py::buffer &b) {
+    typedef typename Type::Scalar Scalar;
+
+    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 && info.strides[0] == sizeof(Scalar)) ||
+          (info.ndim == 2 &&
+              ((info.shape[0] == 1 && info.strides[0] == sizeof(Scalar) &&
+                info.shape[1] == 3) ||
+               (info.shape[1] == 1 && info.strides[1] == sizeof(Scalar) &&
+                info.shape[0] == 3)))))
+        throw std::runtime_error("Incompatible buffer dimension!");
+
+    memcpy(v.data(), info.ptr, sizeof(Scalar) * 3);
+}
+
+/// Creates Python bindings for an Eigen order-1 tensor of size 3 (i.e. a vector/normal/point)
+template <typename Type>
+py::class_<Type> bind_eigen_1_3(py::module &m, const char *name,
+                                py::object parent = py::object()) {
+    typedef typename Type::Scalar Scalar;
+
+    py::class_<Type> vector(m, name, parent);
+    vector
+        /* Constructors */
+        .def(py::init<>())
+        .def(py::init<Scalar>())
+        .def(py::init<Scalar, Scalar, Scalar>())
+        .def("__init__", [](Type &v, const std::vector<Scalar> &v2) {
+            if (v2.size() != 3)
+                throw std::runtime_error("Incompatible size!");
+            memcpy(v.data(), &v2[0], sizeof(Scalar) * 3);
+        })
+        .def("__init__", [](Type &v, py::buffer b) {
+            init_fixed_from_buffer_3(v, b);
+        })
+
+        /* Initialization */
+        .def("setConstant", [](Type &m, Scalar value) { m.setConstant(value); })
+        .def("setZero", [](Type &m) { m.setZero(); })
+
+        /* Arithmetic operators (def_cast forcefully casts the result back to a
+           Matrix 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 * Scalar())
+        .def_cast(py::self / Scalar())
+        .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)
+
+        /* Python protocol implementations */
+        .def("__len__", [](const Type &) { return (int) 3; })
+        .def("__repr__", [](const Type &v) {
+            std::ostringstream oss;
+            oss << v;
+            return oss.str();
+        })
+        .def("__getitem__", [](const Type &c, int i) {
+            if (i < 0 || i >= 3)
+                throw py::index_error();
+            return c[i];
+         })
+        .def("__setitem__", [](Type &c, int i, Scalar v) {
+             if (i < 0 || i >= 3)
+                 throw py::index_error();
+            c[i] = 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(),
+                1, { (size_t) 3 },
+                { sizeof(Scalar) }
+            );
+        });
+    return vector;
+}
+
+/// Creates Python bindings for a dynamic Eigen order-1 tensor (i.e. a vector)
+template <typename Type>
+py::class_<Type> bind_eigen_1(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> vector(m, name, parent);
+
+    vector
+        /* Constructors */
+        .def(py::init<>())
+        .def(py::init<size_t>())
+        .def("__init__", [](Type &v, const std::vector<Scalar> &v2) {
+            new (&v) Type(v2.size());
+            memcpy(v.data(), &v2[0], sizeof(Scalar) * v2.size());
+        })
+        .def("__init__", [](Type &v, py::buffer b) {
+            py::buffer_info info = b.request();
+            if (info.format != py::format_descriptor<Scalar>::value()) {
+                throw std::runtime_error("Incompatible buffer format!");
+            } else if (info.ndim == 1 && info.strides[0] == sizeof(Scalar)) {
+                new (&v) Type(info.shape[0]);
+                memcpy(v.data(), info.ptr, sizeof(Scalar) * info.shape[0]);
+            } else if (info.ndim == 2 && ((info.shape[0] == 1 && info.strides[0] == sizeof(Scalar))
+                                       || (info.shape[1] == 1 && info.strides[1] == sizeof(Scalar)))) {
+                new (&v) Type(info.shape[0] * info.shape[1]);
+                memcpy(v.data(), info.ptr, sizeof(Scalar) * info.shape[0] * info.shape[1]);
+            } else {
+                throw std::runtime_error("Incompatible buffer dimension!");
+            }
+        })
+
+        /* Size query functions */
+        .def("size", [](const Type &m) { return m.size(); })
+        .def("cols", &Type::cols)
+        .def("rows", &Type::rows)
+
+        /* Initialization */
+        .def("setZero", [](Type &m) { m.setZero(); })
+        .def("setConstant", [](Type &m, Scalar value) { m.setConstant(value); })
+
+        /* Resizing */
+        .def("resize", [](Type &m, size_t s0) { m.resize(s0); })
+        .def("resizeLike", [](Type &m, const Type &m2) { m.resizeLike(m2); })
+        .def("conservativeResize", [](Type &m, size_t s0) { m.conservativeResize(s0); })
+
+        /* 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 * Scalar())
+        .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)
+
+        /* Python protocol implementations */
+        .def("__repr__", [](const Type &v) {
+            std::ostringstream oss;
+            oss << v.transpose();
+            return oss.str();
+        })
+        .def("__getitem__", [](const Type &m, size_t i) {
+            if (i >= (size_t) m.size())
+                throw py::index_error();
+            return m[i];
+         })
+        .def("__setitem__", [](Type &m, size_t i, Scalar v) {
+            if (i >= (size_t) m.size())
+                throw py::index_error();
+            m[i] = 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(),
+                1,                       /* Number of dimensions */
+                { (size_t) m.size() },   /* Buffer dimensions */
+                { sizeof(Scalar) }       /* Strides (in bytes) for each index */
+            );
+         })
+
+        /* Static initializers */
+        .def_static("Zero", [](size_t n) { return Type(Type::Zero(n)); })
+        .def_static("Ones", [](size_t n) { return Type(Type::Ones(n)); })
+        .def_static("Constant", [](size_t n, Scalar value) { return Type(Type::Constant(n, value)); });
+    return vector;
+}
+
+/// Creates Python bindings for a dynamic Eigen order-2 tensor (i.e. a matrix)
+template <typename Type>
+py::class_<Type> bind_eigen_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", &Type::cols)
+        .def("rows", &Type::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_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");
+
+    bind_eigen_1<Eigen::VectorXd> (me, "VectorXd");
+    bind_eigen_1<Eigen::VectorXi> (me, "VectorXi");
+    bind_eigen_2<Eigen::MatrixXd> (me, "MatrixXd");
+    bind_eigen_2<Eigen::MatrixXi> (me, "MatrixXi");
+
+
+    /* Bindings for <vector.h> */
+    auto vector3 = bind_eigen_1_3<Eigen::Vector3d>(me, "Vector3d");
+    vector3
+        .def("norm", [](const Eigen::Vector3d &v) { return v.norm(); })
+        .def("squaredNorm", [](const Eigen::Vector3d &v) { return v.squaredNorm(); })
+        .def("normalize", [](Eigen::Vector3d &v) { v.normalize(); })
+        .def("normalized", [](const Eigen::Vector3d &v) -> Eigen::Vector3d { return v.normalized(); })
+        .def("dot", [](const Eigen::Vector3d &v1, const Eigen::Vector3d &v2) { return v1.dot(v2); })
+        .def("cross", [](const Eigen::Vector3d &v1, const Eigen::Vector3d &v2) -> Eigen::Vector3d { return v1.cross(v2); })
+        .def_property("x", [](const Eigen::Vector3d &v) -> double { return v.x(); },
+                           [](Eigen::Vector3d &v, double x) { v.x() = x; }, "X coordinate")
+        .def_property("y", [](const Eigen::Vector3d &v) -> double { return v.y(); },
+                           [](Eigen::Vector3d &v, double y) { v.y() = y; }, "Y coordinate")
+        .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>();
+}
diff --git a/python/python.cpp b/python/python.cpp
new file mode 100644
index 000000000..83543ce2d
--- /dev/null
+++ b/python/python.cpp
@@ -0,0 +1,17 @@
+#include "python.h"
+#include <sstream>
+#include <string>
+#include <fstream>
+
+extern void python_export_vector(py::module &);
+extern void python_export_igl(py::module &);
+
+PYTHON_PLUGIN(igl) {
+    py::init_threading();
+    py::module m("igl", "Python wrappers for libigl");
+
+    python_export_vector(m);
+    python_export_igl(m);
+
+    return m.ptr();
+}
diff --git a/python/python.h b/python/python.h
new file mode 100644
index 000000000..fb412ba6f
--- /dev/null
+++ b/python/python.h
@@ -0,0 +1,11 @@
+#pragma once
+
+#include <pybind/pybind.h>
+#include <pybind/operators.h>
+#include <pybind/complex.h>
+#include <pybind/numpy.h>
+#include <pybind/stl.h>
+
+#include "py_doc.h"
+
+namespace py = pybind;