diff --git a/CHANGELOG.md b/CHANGELOG.md
index f6d89a3..18969ab 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -15,6 +15,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Changed
 
+  - `subdivision`: pass-by-value → Eigen::Ref<const> (3 functions)
+  - `straight_skeleton_2`: const& → Eigen::Ref<const> (8 functions)
+  - `triangulation`: added const to all Eigen::Ref params (3 functions)
+  - `meshing`: removed extra & and added const (7 functions)
+
 ### Removed
 
 
diff --git a/src/meshing.cpp b/src/meshing.cpp
index 5d1794c..3ac44f4 100644
--- a/src/meshing.cpp
+++ b/src/meshing.cpp
@@ -39,8 +39,8 @@ typedef boost::graph_traits<Dual>::edge_descriptor     edge_descriptor;
 
 std::tuple<compas::RowMatrixXd, compas::RowMatrixXi>
 pmp_trimesh_remesh(
-    Eigen::Ref<compas::RowMatrixXd> vertices_a,
-    Eigen::Ref<compas::RowMatrixXi> faces_a,
+    Eigen::Ref<const compas::RowMatrixXd> vertices_a,
+    Eigen::Ref<const compas::RowMatrixXi> faces_a,
     double target_edge_length,
     unsigned int number_of_iterations,
     bool do_project)
@@ -132,13 +132,13 @@ std::vector<double> get_vertex_areas(compas::Mesh mesh, std::vector<double> face
 }
 
 std::tuple<
-compas::RowMatrixXd, 
+compas::RowMatrixXd,
 compas::RowMatrixXi,
-compas::RowMatrixXd, 
+compas::RowMatrixXd,
 std::vector<std::vector<int>>>
 pmp_trimesh_remesh_dual(
-    Eigen::Ref<compas::RowMatrixXd> vertices_a,
-    Eigen::Ref<compas::RowMatrixXi> faces_a,
+    Eigen::Ref<const compas::RowMatrixXd> vertices_a,
+    Eigen::Ref<const compas::RowMatrixXi> faces_a,
     const std::vector<int>& fixed_vertices,
     double length_factor,
     unsigned int number_of_iterations,
@@ -693,8 +693,8 @@ pmp_trimesh_remesh_dual(
 }
 
 void pmp_pull(
-    Eigen::Ref<compas::RowMatrixXd> vertices_a,
-    Eigen::Ref<compas::RowMatrixXi> faces_a,
+    Eigen::Ref<const compas::RowMatrixXd> vertices_a,
+    Eigen::Ref<const compas::RowMatrixXi> faces_a,
     Eigen::Ref<compas::RowMatrixXd> vertices_b,
     Eigen::Ref<compas::RowMatrixXd> normals_b)
 {
@@ -823,8 +823,8 @@ void pmp_pull(
 
 
 void pmp_project(
-    Eigen::Ref<compas::RowMatrixXd> vertices_a,
-    Eigen::Ref<compas::RowMatrixXi> faces_a,
+    Eigen::Ref<const compas::RowMatrixXd> vertices_a,
+    Eigen::Ref<const compas::RowMatrixXi> faces_a,
     Eigen::Ref<compas::RowMatrixXd> vertices_b)
 {
     /////////////////////////////////////////////////////////////////////////////////
diff --git a/src/meshing.h b/src/meshing.h
index 6378ec5..1dafebd 100644
--- a/src/meshing.h
+++ b/src/meshing.h
@@ -20,8 +20,8 @@ namespace compas {
      * @return Surface area of the mesh
      */
     double pmp_area(
-        Eigen::Ref<const RowMatrixXd>& vertices,
-        Eigen::Ref<const RowMatrixXi>& faces);
+        Eigen::Ref<const RowMatrixXd> vertices,
+        Eigen::Ref<const RowMatrixXi> faces);
 
     /**
      * @brief Calculate the volume enclosed by a mesh
@@ -30,8 +30,8 @@ namespace compas {
      * @return Volume enclosed by the mesh
      */
     double pmp_volume(
-        Eigen::Ref<const RowMatrixXd>& vertices,
-        Eigen::Ref<const RowMatrixXi>& faces);
+        Eigen::Ref<const RowMatrixXd> vertices,
+        Eigen::Ref<const RowMatrixXi> faces);
 
     /**
      * @brief Calculate the centroid of a mesh
@@ -40,8 +40,8 @@ namespace compas {
      * @return Coordinates of the mesh centroid
      */
     std::vector<double> pmp_centroid(
-        Eigen::Ref<const RowMatrixXd>& vertices,
-        Eigen::Ref<const RowMatrixXi>& faces);
+        Eigen::Ref<const RowMatrixXd> vertices,
+        Eigen::Ref<const RowMatrixXi> faces);
 
 } // namespace compas
 
@@ -59,8 +59,8 @@ namespace compas {
  */
 std::tuple<compas::RowMatrixXd, compas::RowMatrixXi>
 pmp_trimesh_remesh(
-    Eigen::Ref<compas::RowMatrixXd> vertices_a,
-    Eigen::Ref<compas::RowMatrixXi> faces_a,
+    Eigen::Ref<const compas::RowMatrixXd> vertices_a,
+    Eigen::Ref<const compas::RowMatrixXi> faces_a,
     double target_edge_length,
     unsigned int number_of_iterations = 10,
     bool do_project = true);
@@ -79,13 +79,13 @@ pmp_trimesh_remesh(
  */
 
 std::tuple<
-compas::RowMatrixXd, 
+compas::RowMatrixXd,
 compas::RowMatrixXi,
-compas::RowMatrixXd, 
-std::vector<std::vector<int>>> 
+compas::RowMatrixXd,
+std::vector<std::vector<int>>>
 pmp_trimesh_remesh_dual(
-    Eigen::Ref<compas::RowMatrixXd> vertices_a,
-    Eigen::Ref<compas::RowMatrixXi> faces_a,
+    Eigen::Ref<const compas::RowMatrixXd> vertices_a,
+    Eigen::Ref<const compas::RowMatrixXi> faces_a,
     const std::vector<int>& fixed_vertices,
     double length_factor=1.0,
     unsigned int number_of_iterations=10,
@@ -103,8 +103,8 @@ pmp_trimesh_remesh_dual(
  */
 
 void pmp_pull(
-    Eigen::Ref<compas::RowMatrixXd> vertices_a,
-    Eigen::Ref<compas::RowMatrixXi> faces_a,
+    Eigen::Ref<const compas::RowMatrixXd> vertices_a,
+    Eigen::Ref<const compas::RowMatrixXi> faces_a,
     Eigen::Ref<compas::RowMatrixXd> vertices_b,
     Eigen::Ref<compas::RowMatrixXd> normals_b);
 
@@ -116,6 +116,6 @@ void pmp_pull(
  * @param vertices_b Matrix of vertex positions as Nx3 matrix in row-major order (float64)
  */
 void pmp_project(
-    Eigen::Ref<compas::RowMatrixXd> vertices_a,
-    Eigen::Ref<compas::RowMatrixXi> faces_a,
+    Eigen::Ref<const compas::RowMatrixXd> vertices_a,
+    Eigen::Ref<const compas::RowMatrixXi> faces_a,
     Eigen::Ref<compas::RowMatrixXd> vertices_b);
diff --git a/src/straight_skeleton_2.cpp b/src/straight_skeleton_2.cpp
index 0819cc7..882f949 100644
--- a/src/straight_skeleton_2.cpp
+++ b/src/straight_skeleton_2.cpp
@@ -126,7 +126,7 @@ data_to_polygon_with_holes(
 
 std::tuple<compas::RowMatrixXd, std::vector<int>, compas::RowMatrixXi, std::vector<int>>
 pmp_create_interior_straight_skeleton(
-    const compas::RowMatrixXd& vertices)
+    Eigen::Ref<const compas::RowMatrixXd> vertices)
 {
     Polygon_2 polygon = data_to_polygon(vertices);
     SsPtr skeleton = CGAL::create_interior_straight_skeleton_2(polygon.vertices_begin(), polygon.vertices_end());
@@ -135,7 +135,7 @@ pmp_create_interior_straight_skeleton(
 
 std::tuple<compas::RowMatrixXd, std::vector<int>, compas::RowMatrixXi, std::vector<int>>
 pmp_create_interior_straight_skeleton_with_holes(
-    const compas::RowMatrixXd& boundary_vertices,
+    Eigen::Ref<const compas::RowMatrixXd> boundary_vertices,
     const std::vector<compas::RowMatrixXd>& hole_vertices)
 {
     Polygon_with_holes polygon = data_to_polygon_with_holes(boundary_vertices, hole_vertices);
@@ -145,7 +145,7 @@ pmp_create_interior_straight_skeleton_with_holes(
 
 std::vector<compas::RowMatrixXd>
 pmp_create_offset_polygons_2_inner(
-    const compas::RowMatrixXd& vertices,
+    Eigen::Ref<const compas::RowMatrixXd> vertices,
     double& offset_distance)
 {
     Polygon_2 polygon = data_to_polygon(vertices);
@@ -161,7 +161,7 @@ pmp_create_offset_polygons_2_inner(
 
 std::vector<std::vector<compas::RowMatrixXd>>
 pmp_create_offset_polygons_2_inner_with_holes(
-    const compas::RowMatrixXd& boundary_vertices,
+    Eigen::Ref<const compas::RowMatrixXd> boundary_vertices,
     const std::vector<compas::RowMatrixXd>& hole_vertices,
     double& offset_distance)
 {
@@ -188,7 +188,7 @@ pmp_create_offset_polygons_2_inner_with_holes(
 
 std::vector<compas::RowMatrixXd>
 pmp_create_offset_polygons_2_outer(
-    const compas::RowMatrixXd& vertices,
+    Eigen::Ref<const compas::RowMatrixXd> vertices,
     double& offset_distance)
 {
     Polygon_2 polygon = data_to_polygon(vertices);
@@ -204,7 +204,7 @@ pmp_create_offset_polygons_2_outer(
 
 std::vector<std::vector<compas::RowMatrixXd>>
 pmp_create_offset_polygons_2_outer_with_holes(
-    const compas::RowMatrixXd& boundary_vertices,
+    Eigen::Ref<const compas::RowMatrixXd> boundary_vertices,
     const std::vector<compas::RowMatrixXd>& hole_vertices,
     double& offset_distance)
 {
@@ -226,9 +226,9 @@ pmp_create_offset_polygons_2_outer_with_holes(
 
 std::vector<compas::RowMatrixXd>
 pmp_create_weighted_offset_polygons_2_inner(
-    const compas::RowMatrixXd& vertices,
+    Eigen::Ref<const compas::RowMatrixXd> vertices,
     double offset_distance,
-    const compas::RowMatrixXd& edge_weights)
+    Eigen::Ref<const compas::RowMatrixXd> edge_weights)
 {
     if (edge_weights.rows() != vertices.rows()) {
         throw std::invalid_argument("Number of weights must match number of polygon vertices");
@@ -264,9 +264,9 @@ pmp_create_weighted_offset_polygons_2_inner(
 
 std::vector<compas::RowMatrixXd>
 pmp_create_weighted_offset_polygons_2_outer(
-    const compas::RowMatrixXd& vertices,
+    Eigen::Ref<const compas::RowMatrixXd> vertices,
     double offset_distance,
-    const compas::RowMatrixXd& edge_weights)
+    Eigen::Ref<const compas::RowMatrixXd> edge_weights)
 {
     if (edge_weights.rows() != vertices.rows()) {
         throw std::invalid_argument("Number of weights must match number of polygon vertices");
diff --git a/src/straight_skeleton_2.h b/src/straight_skeleton_2.h
index 748093d..a6fcfa8 100644
--- a/src/straight_skeleton_2.h
+++ b/src/straight_skeleton_2.h
@@ -22,9 +22,9 @@
  *         - Matrix of skeleton edges as vertex pairs (Kx2, int32)
  *         - Vector of source edge indices from input polygon for each skeleton edge
  */
-std::tuple<compas::RowMatrixXd, std::vector<int>, compas::RowMatrixXi, std::vector<int>> 
+std::tuple<compas::RowMatrixXd, std::vector<int>, compas::RowMatrixXi, std::vector<int>>
 pmp_create_interior_straight_skeleton(
-    const compas::RowMatrixXd& vertices
+    Eigen::Ref<const compas::RowMatrixXd> vertices
 );
 
 /**
@@ -38,9 +38,9 @@ pmp_create_interior_straight_skeleton(
  *         - Matrix of skeleton edges as vertex pairs (Qx2, int32)
  *         - Vector of source edge indices from input polygon for each skeleton edge
  */
-std::tuple<compas::RowMatrixXd, std::vector<int>, compas::RowMatrixXi, std::vector<int>> 
+std::tuple<compas::RowMatrixXd, std::vector<int>, compas::RowMatrixXi, std::vector<int>>
 pmp_create_interior_straight_skeleton_with_holes(
-    const compas::RowMatrixXd& boundary_vertices,
+    Eigen::Ref<const compas::RowMatrixXd> boundary_vertices,
     const std::vector<compas::RowMatrixXd>& hole_vertices
 );
 
@@ -53,7 +53,7 @@ pmp_create_interior_straight_skeleton_with_holes(
  */
 std::vector<compas::RowMatrixXd>
 pmp_create_offset_polygons_2_inner(
-    const compas::RowMatrixXd& vertices,
+    Eigen::Ref<const compas::RowMatrixXd> vertices,
     double& offset_distance
 );
 
@@ -68,7 +68,7 @@ pmp_create_offset_polygons_2_inner(
  */
 std::vector<std::vector<compas::RowMatrixXd>>
 pmp_create_offset_polygons_2_inner_with_holes(
-    const compas::RowMatrixXd& boundary_vertices,
+    Eigen::Ref<const compas::RowMatrixXd> boundary_vertices,
     const std::vector<compas::RowMatrixXd>& hole_vertices,
     double& offset_distance
 );
@@ -88,7 +88,7 @@ pmp_create_offset_polygons_2_inner_with_holes(
  */
 std::vector<compas::RowMatrixXd>
 pmp_create_offset_polygons_2_outer(
-    const compas::RowMatrixXd& vertices,
+    Eigen::Ref<const compas::RowMatrixXd> vertices,
     double& offset_distance
 );
 
@@ -103,7 +103,7 @@ pmp_create_offset_polygons_2_outer(
  */
 std::vector<std::vector<compas::RowMatrixXd>>
 pmp_create_offset_polygons_2_outer_with_holes(
-    const compas::RowMatrixXd& boundary_vertices,
+    Eigen::Ref<const compas::RowMatrixXd> boundary_vertices,
     const std::vector<compas::RowMatrixXd>& hole_vertices,
     double& offset_distance
 );
@@ -120,9 +120,9 @@ pmp_create_offset_polygons_2_outer_with_holes(
  */
 std::vector<compas::RowMatrixXd>
 pmp_create_weighted_offset_polygons_2_inner(
-    const compas::RowMatrixXd& vertices,
+    Eigen::Ref<const compas::RowMatrixXd> vertices,
     double offset_distance,
-    const compas::RowMatrixXd& edge_weights
+    Eigen::Ref<const compas::RowMatrixXd> edge_weights
 );
 
 /**
@@ -137,9 +137,9 @@ pmp_create_weighted_offset_polygons_2_inner(
  */
 std::vector<compas::RowMatrixXd>
 pmp_create_weighted_offset_polygons_2_outer(
-    const compas::RowMatrixXd& vertices,
+    Eigen::Ref<const compas::RowMatrixXd> vertices,
     double offset_distance,
-    const compas::RowMatrixXd& edge_weights
+    Eigen::Ref<const compas::RowMatrixXd> edge_weights
 );
 
 void init_straight_skeleton_2(nb::module_& m);
\ No newline at end of file
diff --git a/src/subdivision.cpp b/src/subdivision.cpp
index b712206..fba52b4 100644
--- a/src/subdivision.cpp
+++ b/src/subdivision.cpp
@@ -2,8 +2,8 @@
 
 std::tuple<compas::RowMatrixXd, compas::RowMatrixXi>
 subd_catmullclark(
-    compas::RowMatrixXd vertices,
-    std::vector<std::vector<int>> faces,
+    Eigen::Ref<const compas::RowMatrixXd> vertices,
+    const std::vector<std::vector<int>>& faces,
     unsigned int num_iterations)
 {
     compas::Mesh mesh = compas::ngon_from_vertices_and_faces(vertices, faces);
@@ -14,8 +14,8 @@ subd_catmullclark(
 
 std::tuple<compas::RowMatrixXd, compas::RowMatrixXi>
 subd_loop(
-    compas::RowMatrixXd vertices,
-    compas::RowMatrixXi faces,
+    Eigen::Ref<const compas::RowMatrixXd> vertices,
+    Eigen::Ref<const compas::RowMatrixXi> faces,
     unsigned int num_iterations)
 {
     compas::Mesh mesh = compas::mesh_from_vertices_and_faces(vertices, faces);
@@ -26,8 +26,8 @@ subd_loop(
 
 std::tuple<compas::RowMatrixXd, compas::RowMatrixXi>
 subd_sqrt3(
-    compas::RowMatrixXd vertices,
-    compas::RowMatrixXi faces,
+    Eigen::Ref<const compas::RowMatrixXd> vertices,
+    Eigen::Ref<const compas::RowMatrixXi> faces,
     unsigned int num_iterations)
 {
     compas::Mesh mesh = compas::mesh_from_vertices_and_faces(vertices, faces);
diff --git a/src/subdivision.h b/src/subdivision.h
index 6babb1b..0abd9af 100644
--- a/src/subdivision.h
+++ b/src/subdivision.h
@@ -18,8 +18,8 @@
  */
 std::tuple<compas::RowMatrixXd, compas::RowMatrixXi>
 subd_catmullclark(
-    compas::RowMatrixXd vertices,
-    std::vector<std::vector<int>> faces,
+    Eigen::Ref<const compas::RowMatrixXd> vertices,
+    const std::vector<std::vector<int>>& faces,
     unsigned int num_iterations);
 
 /**
@@ -35,8 +35,8 @@ subd_catmullclark(
  */
 std::tuple<compas::RowMatrixXd, compas::RowMatrixXi>
 subd_loop(
-    compas::RowMatrixXd vertices,
-    compas::RowMatrixXi faces,
+    Eigen::Ref<const compas::RowMatrixXd> vertices,
+    Eigen::Ref<const compas::RowMatrixXi> faces,
     unsigned int num_iterations);
 
 /**
@@ -52,6 +52,6 @@ subd_loop(
  */
 std::tuple<compas::RowMatrixXd, compas::RowMatrixXi>
 subd_sqrt3(
-    compas::RowMatrixXd vertices,
-    compas::RowMatrixXi faces,
+    Eigen::Ref<const compas::RowMatrixXd> vertices,
+    Eigen::Ref<const compas::RowMatrixXi> faces,
     unsigned int num_iterations);
\ No newline at end of file
diff --git a/src/triangulation.cpp b/src/triangulation.cpp
index b4b5f55..084b9d1 100644
--- a/src/triangulation.cpp
+++ b/src/triangulation.cpp
@@ -30,7 +30,7 @@ using cdt3 = CGAL::Constrained_Delaunay_triangulation_2<compas::Kernel, tds3, it
 using criteria = CGAL::Delaunay_mesh_size_criteria_2<cdt3>;
 
 compas::RowMatrixXi
-pmp_delaunay_triangulation(Eigen::Ref<compas::RowMatrixXd> vertices)
+pmp_delaunay_triangulation(Eigen::Ref<const compas::RowMatrixXd> vertices)
 {
     dt triangulation;
     std::vector<std::pair<dt::Point, unsigned>> points(vertices.rows());
@@ -130,8 +130,8 @@ void mark_domains(cdt2 &triangulation)
 
 std::tuple<compas::RowMatrixXd, compas::RowMatrixXi>
 pmp_constrained_delaunay_triangulation(
-    Eigen::Ref<compas::RowMatrixXd> boundary_vertices,
-    Eigen::Ref<compas::RowMatrixXd> internal_vertices,
+    Eigen::Ref<const compas::RowMatrixXd> boundary_vertices,
+    Eigen::Ref<const compas::RowMatrixXd> internal_vertices,
     const std::vector<compas::RowMatrixXd> & holes,
     const std::vector<compas::RowMatrixXd> & curves,
     bool is_conforming
@@ -249,8 +249,8 @@ pmp_constrained_delaunay_triangulation(
 
 std::tuple<compas::RowMatrixXd, compas::RowMatrixXi>
 pmp_refined_delaunay_mesh(
-    Eigen::Ref<compas::RowMatrixXd> boundary_vertices,
-    Eigen::Ref<compas::RowMatrixXd> internal_vertices,
+    Eigen::Ref<const compas::RowMatrixXd> boundary_vertices,
+    Eigen::Ref<const compas::RowMatrixXd> internal_vertices,
     const std::vector<compas::RowMatrixXd> & holes,
     const std::vector<compas::RowMatrixXd> & curves,
     double min_angle,
diff --git a/src/triangulation.h b/src/triangulation.h
index 46b046d..1bc0b68 100644
--- a/src/triangulation.h
+++ b/src/triangulation.h
@@ -24,7 +24,7 @@
  */
 compas::RowMatrixXi
 pmp_delaunay_triangulation(
-    Eigen::Ref<compas::RowMatrixXd> vertices);
+    Eigen::Ref<const compas::RowMatrixXd> vertices);
 
 
 /**
@@ -39,8 +39,8 @@ pmp_delaunay_triangulation(
  */
 std::tuple<compas::RowMatrixXd, compas::RowMatrixXi>
 pmp_constrained_delaunay_triangulation(
-    Eigen::Ref<compas::RowMatrixXd> boundary_vertices,
-    Eigen::Ref<compas::RowMatrixXd> internal_vertices,
+    Eigen::Ref<const compas::RowMatrixXd> boundary_vertices,
+    Eigen::Ref<const compas::RowMatrixXd> internal_vertices,
     const std::vector<compas::RowMatrixXd> & holes,
     const std::vector<compas::RowMatrixXd> & curves,
     bool is_conforming = false
@@ -60,8 +60,8 @@ pmp_constrained_delaunay_triangulation(
  */
 std::tuple<compas::RowMatrixXd, compas::RowMatrixXi>
 pmp_refined_delaunay_mesh(
-    Eigen::Ref<compas::RowMatrixXd> boundary_vertices,
-    Eigen::Ref<compas::RowMatrixXd> internal_vertices,
+    Eigen::Ref<const compas::RowMatrixXd> boundary_vertices,
+    Eigen::Ref<const compas::RowMatrixXd> internal_vertices,
     const std::vector<compas::RowMatrixXd> & holes,
     const std::vector<compas::RowMatrixXd> & curves,
     double min_angle = 0.0,