diff --git a/CMakeLists.txt b/CMakeLists.txt
index bb824d5ab237a855ea8719e71a4f248fc031b6bb..4148520e7b8adaa51f7cbadae02f828cf22e8174 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -17,6 +17,7 @@ set(CMAKE_CXX_STANDARD_REQUIRED ON)
 find_package(Python REQUIRED COMPONENTS Interpreter Development.Module)
 find_package(pybind11 CONFIG REQUIRED)
 find_package(CGAL REQUIRED COMPONENTS ImageIO)
+find_package(Eigen3 3.3 REQUIRED NO_MODULE)
 
 # make cache variables for install destinations
 include(GNUInstallDirs)
diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt
index d560f0225c04a4e104a8b438b7456ff6d6bcdf66..2aed2d18fd05ddad6c473103b376761b75f3de52 100644
--- a/cpp/CMakeLists.txt
+++ b/cpp/CMakeLists.txt
@@ -12,5 +12,6 @@ add_library(${PROJECT_NAME} SHARED
     ${${PROJECT_NAME}_HEADERS}
     ${${PROJECT_NAME}_SOURCES}
 )
+target_link_libraries(${PROJECT_NAME} PUBLIC Eigen3::Eigen)
 
 install(TARGETS ${PROJECT_NAME} DESTINATION "${CMAKE_INSTALL_PREFIX}/${CMAKE_INSTALL_LIBDIR}")
\ No newline at end of file
diff --git a/cpp/LabelledTetrahedralMesh.cpp b/cpp/LabelledTetrahedralMesh.cpp
index 0e4d977734fffe55cb17a5488ca40b077626a829..e7ab4d2ec7cfaf467f77615084409c63a22a8ba2 100644
--- a/cpp/LabelledTetrahedralMesh.cpp
+++ b/cpp/LabelledTetrahedralMesh.cpp
@@ -1,12 +1,19 @@
 #include "LabelledTetrahedralMesh.h"
 
 
-std::vector< std::vector<double> > LabelledTetrahedralMesh::vertexPoints(void) const
+Eigen::ArrayX3d LabelledTetrahedralMesh::vertexPoints(void) const
 {
-    return this->_vertex_points;
+    Eigen::ArrayX3d vertex_points(this->_vertex_points.size(), 3);
+    for (int i=0; i<this->_vertex_points.size(); i++) {
+        for (int dim=0; dim<3; dim++) {
+            vertex_points(i, dim) = this->_vertex_points[i][dim];
+        }
+    }
+    return vertex_points;
 }
 
-void LabelledTetrahedralMesh::addPoint(const std::vector<double>& point)
+void LabelledTetrahedralMesh::addPoint(Eigen::Array3d point)
 {
-    this->_vertex_points.push_back(point);
+    std::vector<double> _point = std::vector<double> {point(0, 0), point(0, 1), point(0, 2)};
+    this->_vertex_points.push_back(_point);
 }
diff --git a/cpp/LabelledTetrahedralMesh.h b/cpp/LabelledTetrahedralMesh.h
index 16bf3cfbc7025ee1bb64596e4db347aee0d1540b..c6be404731067ed8f319b01d4d270bad909b7c33 100644
--- a/cpp/LabelledTetrahedralMesh.h
+++ b/cpp/LabelledTetrahedralMesh.h
@@ -1,6 +1,7 @@
 #pragma once
 
 #include <vector>
+#include <Eigen/Dense>
 
 class LabelledTetrahedralMesh
 {
@@ -9,8 +10,8 @@ public:
     ~LabelledTetrahedralMesh(void) = default;
 
 public:
-    std::vector< std::vector<double> > vertexPoints(void) const;
-    void addPoint(const std::vector<double>& point);
+    Eigen::ArrayX3d vertexPoints(void) const;
+    void addPoint(Eigen::Array3d point);
 
 private:
     std::vector< std::vector<double> > _vertex_points;
diff --git a/tests/test_tetrahedral_mesh.py b/tests/test_tetrahedral_mesh.py
index 1ef3268bed30bb8d7f82cfa2b0a91287cd307769..3ae1ffd962f7e0e0392a5e40869eebedc3fc25cb 100644
--- a/tests/test_tetrahedral_mesh.py
+++ b/tests/test_tetrahedral_mesh.py
@@ -13,4 +13,4 @@ class TestTetrahedralMesh(unittest.TestCase):
         mesh = LabelledTetrahedralMesh()
         mesh.add_point([0, 0, 0])
         points = mesh.vertex_points()
-        assert len(points) == 1
+        assert points.shape == (1, 3)
diff --git a/wrp/LabelledTetrahedralMeshWrapper.cpp b/wrp/LabelledTetrahedralMeshWrapper.cpp
index 79ce03db1162cf8c89865b9e5a67d84ef59cb16a..e07adf85fb2e390cb106ff81d5fb56ec37f44426 100644
--- a/wrp/LabelledTetrahedralMeshWrapper.cpp
+++ b/wrp/LabelledTetrahedralMeshWrapper.cpp
@@ -1,5 +1,6 @@
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h> // This header is needed to work with STL containers like std::vector
+#include <pybind11/eigen.h>
 
 #include "LabelledTetrahedralMesh.h"