diff --git a/docs/algorithms.rst b/docs/algorithms.rst
index 8fae8e20fe4ac60bf8a2f712cf841eb89ce10252..f1779f9ba2c10e447f9780ccaf17e9facb017121 100644
--- a/docs/algorithms.rst
+++ b/docs/algorithms.rst
@@ -1,4 +1,65 @@
 Algorithms
 ==========
 
-TODO
+The headers located in ``polymesh/algorithms/*`` or included by ``polymesh/algorithms.hh`` contain useful algorithms operating on meshes.
+In contrast to the :doc:`properties` or the functional-inspired :doc:`smart-ranges`, these algorithms are often non-trivial and perform substantial transformations on the mesh.
+As most built-in geometry-related methods, the algorithms tend to be very generic, providing extensive customization points and using the :doc:`vector-math` wrapper to support many different math types.
+
+A reference of all algorithms can be found in :ref:`algorithms-ref`.
+
+Components
+----------
+
+This category of algorithms contains methods to compute connected components on meshes.
+
+::
+
+    #include <polymesh/algorithms/components.hh>
+
+    pm::Mesh m;
+    load(m, /* ... */);
+
+    // computes all vertex components of the mesh
+    int vcomp_count;
+    auto vcomp = pm::vertex_components(m, vcomp_count);
+    for (auto v : m.vertices())
+        std::cout << "vertex " << int(v) << " belongs to component " << vcomp[v] << " of " << vcomp_count << std::endl;
+
+    // special iterator to iterate over a single component
+    pm::face_handle f = /* ... */;
+    for (auto ff : pm::face_component(f))
+        std::cout << "face " << int(ff) << " is in the same connected component as face " << int(f) << std::endl;
+
+
+.. doxygenfunction:: polymesh::vertex_components
+
+.. doxygenfunction:: polymesh::face_components
+
+.. doxygenfunction:: polymesh::vertex_component
+
+.. doxygenfunction:: polymesh::face_component(face_handle)
+
+.. doxygenfunction:: polymesh::face_component(vertex_handle)
+
+
+Deduplicate
+-----------
+
+A small helper that merges vertices based on a user criterion.
+
+::
+
+    #include <polymesh/algorithms/deduplicate.hh>
+
+    pm::Mesh m;
+    auto pos = m.vertices().make_attribute<tg::pos3>();
+    load("some-file.stl", m, pos);
+
+    // merges all vertices with the same position
+    pm::deduplicate(m, pos);
+
+.. doxygenfunction:: polymesh::deduplicate
+
+TODO: preserve line breaks in doxygen
+
+TODO: decimate, dedup, delaunay, edge_split, subdivision, interpolation, iteration, normal-estimation, normalize, operations, optimizations, sampling, smoothing, stats, topology, tracing, triangulate
diff --git a/docs/reference.rst b/docs/reference.rst
index c1db0bdd80832036bceb646dc27dbbdeb896ef93..675a7940671feedbd09dff8b4ad3f00983cb9e1e 100644
--- a/docs/reference.rst
+++ b/docs/reference.rst
@@ -160,6 +160,8 @@ Geometric Properties
 
 .. doxygenfunction:: polymesh::bary_interpolate
 
+.. doxygenfunction:: polymesh::barycoords_of
+
 .. doxygenfunction:: polymesh::edge_length(edge_handle, vertex_attribute<Pos3> const&)
 
 .. doxygenfunction:: polymesh::edge_length(halfedge_handle, vertex_attribute<Pos3> const&)
@@ -196,6 +198,24 @@ Geometric Properties
 
 .. doxygenfunction:: polymesh::can_collapse_without_flips
 
+.. _algorithms-ref:
+
+Algorithms
+----------
+
+.. doxygenfunction:: polymesh::vertex_components
+
+.. doxygenfunction:: polymesh::face_components
+
+.. doxygenfunction:: polymesh::vertex_component
+
+.. doxygenfunction:: polymesh::face_component(face_handle)
+
+.. doxygenfunction:: polymesh::face_component(vertex_handle)
+
+.. doxygenfunction:: polymesh::deduplicate
+
+
 Low-Level API
 -------------
 
diff --git a/src/polymesh/Mesh.hh b/src/polymesh/Mesh.hh
index 188f4b9a17f6375a7bba5e08ce059623b881bdb9..a24298eae05702aca18b15da0e912522f71a6b6d 100644
--- a/src/polymesh/Mesh.hh
+++ b/src/polymesh/Mesh.hh
@@ -236,7 +236,7 @@ private:
 };
 } // namespace polymesh
 
-/// ======== IMPLEMENTATIONS ========
+// ======== IMPLEMENTATIONS ========
 
 #include "impl/impl_attributes.hh"
 #include "impl/impl_cursors.hh"
diff --git a/src/polymesh/algorithms/barycoords.hh b/src/polymesh/algorithms/barycoords.hh
deleted file mode 100644
index fcd8114257be85ebffb901dce1b8ff1855af4520..0000000000000000000000000000000000000000
--- a/src/polymesh/algorithms/barycoords.hh
+++ /dev/null
@@ -1,45 +0,0 @@
-#pragma once
-
-#include "../Mesh.hh"
-#include "../fields.hh"
-
-namespace polymesh
-{
-/// calculates the barycentric coordinates of a given point p within a face f
-/// NOTE: asserts that f is triangular
-/// NOTE: also works for other points in the same plane as f
-template <class Pos3, class Scalar = typename field3<Pos3>::scalar_t>
-std::array<Scalar, 3> barycoords(face_handle f, vertex_attribute<Pos3> const& positions, Pos3 p);
-
-/// ======== IMPLEMENTATION ========
-
-template <class Pos3, class Scalar>
-std::array<Scalar, 3> barycoords(face_handle f, vertex_attribute<Pos3> const& positions, Pos3 p)
-{
-    POLYMESH_ASSERT(f.vertices().size() == 3 && "only supports triangles");
-
-    auto ps = f.vertices().to_array<3>(positions);
-
-    auto e10 = ps[1] - ps[0];
-    auto e21 = ps[2] - ps[1];
-
-    auto n = field3<Pos3>::cross(e10, e21);
-
-    auto signed_area = [&](Pos3 const& v0, Pos3 const& v1, Pos3 const& v2) {
-        auto d1 = v1 - v0;
-        auto d2 = v2 - v0;
-
-        auto a = field3<Pos3>::cross(d1, d2);
-
-        return field3<Pos3>::dot(a, n);
-    };
-
-    auto a = signed_area(ps[0], ps[1], ps[2]);
-    auto a0 = signed_area(p, ps[1], ps[2]);
-    auto a1 = signed_area(p, ps[2], ps[0]);
-    auto a2 = signed_area(p, ps[0], ps[1]);
-
-    auto inv_a = Scalar(1) / a;
-    return {a0 * inv_a, a1 * inv_a, a2 * inv_a};
-}
-}
diff --git a/src/polymesh/algorithms/components.hh b/src/polymesh/algorithms/components.hh
index d964c52a28b1c700aad95f2c280af322103d4235..03aa1c23deb9e79d033cac18952cc5f3b9c9cd74 100644
--- a/src/polymesh/algorithms/components.hh
+++ b/src/polymesh/algorithms/components.hh
@@ -23,9 +23,11 @@ detail::bfs_range<vertex_tag> vertex_component(vertex_handle v);
 
 /// returns a range that iterates over all connected faces in BFS order
 detail::bfs_range<face_tag> face_component(face_handle f);
+
+/// returns a range that iterates over all connected faces in BFS order
 detail::bfs_range<face_tag> face_component(vertex_handle v);
 
-/// ======== IMPLEMENTATION ========
+// ======== IMPLEMENTATION ========
 
 inline detail::bfs_range<vertex_tag> vertex_component(vertex_handle v) { return {v}; }
 inline detail::bfs_range<face_tag> face_component(face_handle f) { return {f}; }
diff --git a/src/polymesh/algorithms/decimate.hh b/src/polymesh/algorithms/decimate.hh
index 1f5069f8ceecb7fc558a531047387f45830575a3..b908a30e503a865dce5e4ecca5976e203b5fd4e8 100644
--- a/src/polymesh/algorithms/decimate.hh
+++ b/src/polymesh/algorithms/decimate.hh
@@ -76,7 +76,7 @@ struct decimate_config
     }
 };
 
-/*
+/**
  * Error-function-based incremental decimation
  *
  * Initial per-vertex errors must be provided in vertex_errors
@@ -97,6 +97,7 @@ void decimate(pm::Mesh& m, //
               pm::vertex_attribute<ErrorF>& errors,
               ConfigT const& config);
 
+/// calls decimate with a default configuration that decimates until a target vertex count is reached
 template <class Pos3, class ErrorF>
 void decimate_down_to(pm::Mesh& m, //
                       pm::vertex_attribute<Pos3>& pos,
@@ -106,6 +107,7 @@ void decimate_down_to(pm::Mesh& m, //
     return decimate(m, pos, errors, decimate_config<Pos3, ErrorF>::down_to(target_vertex_cnt));
 }
 
+/// calls decimate with a default configuration that decimates until a target error value is reached
 template <class Pos3, class ErrorF, class ErrorValueT>
 void decimate_up_to_error(pm::Mesh& m, //
                           pm::vertex_attribute<Pos3>& pos,
diff --git a/src/polymesh/algorithms/deduplicate.hh b/src/polymesh/algorithms/deduplicate.hh
index bc455c956c0b9a880330ce7ed12205fa8ea44130..e6cda0a46ab39de578dc8db3cd69f8b1842f5647 100644
--- a/src/polymesh/algorithms/deduplicate.hh
+++ b/src/polymesh/algorithms/deduplicate.hh
@@ -23,10 +23,12 @@ namespace polymesh
 /// CAUTION: currently only works on faces and will remove isolated vertices/edges
 ///
 /// returns number of removed vertices (-1 if deduplication failed (e.g. due to non-manifoldness))
+///
+/// TODO: use a function_ref and implement this in a .cc
 template <class KeyF>
 int deduplicate(Mesh& m, KeyF&& kf);
 
-/// ======== IMPLEMENTATION ========
+// ======== IMPLEMENTATION ========
 
 template <class KeyF>
 int deduplicate(Mesh& m, KeyF&& kf)
diff --git a/src/polymesh/algorithms/interpolation.hh b/src/polymesh/algorithms/interpolation.hh
index 06579ebf21fe4099efc0846d55ab050e21822a0a..29f393f6c44067562a27ba0c36753ccd8b3e5d20 100644
--- a/src/polymesh/algorithms/interpolation.hh
+++ b/src/polymesh/algorithms/interpolation.hh
@@ -49,7 +49,7 @@ T interpolate(handle_t h, attr_t<T> const& attr, W const* ws, int wcnt);
 template <class T, class WeightFuncT, class handle_t, template <class> class attr_t, class Enabled = decltype(std::declval<WeightFuncT>()(0, typename primitive<typename attr_t<T>::tag_t>::handle{}))>
 T interpolate(handle_t h, attr_t<T> const& attr, WeightFuncT&& wh);
 
-/// ======== IMPLEMENTATION ========
+// ======== IMPLEMENTATION ========
 
 namespace detail
 {
diff --git a/src/polymesh/algorithms/operations.hh b/src/polymesh/algorithms/operations.hh
index 38cacb354b20bca92eaf6dbf12c29a625d1ae7eb..37bdc56d2067658ec164b364f1c82fd414a303a3 100644
--- a/src/polymesh/algorithms/operations.hh
+++ b/src/polymesh/algorithms/operations.hh
@@ -13,7 +13,7 @@ void remove_faces(Mesh& m);
 /// NOTE: does NOT compactify!
 void remove_edges_and_faces(Mesh& m);
 
-/// ======== IMPLEMENTATION ========
+// ======== IMPLEMENTATION ========
 
 inline void remove_faces(Mesh& m)
 {
diff --git a/src/polymesh/algorithms/stats.hh b/src/polymesh/algorithms/stats.hh
index 9117c9a87a7648d955ceeee722ec932dd0280043..2507fa48d361f0e6b59910283c69c619dbf30f00 100644
--- a/src/polymesh/algorithms/stats.hh
+++ b/src/polymesh/algorithms/stats.hh
@@ -18,7 +18,7 @@ namespace polymesh
 template <class Vec3 = void>
 void print_stats(std::ostream& out, Mesh const& m, vertex_attribute<Vec3> const* position = nullptr);
 
-/// ======== IMPLEMENTATION ========
+// ======== IMPLEMENTATION ========
 template <class Vec3>
 void print_stats(std::ostream& out, Mesh const& m, vertex_attribute<Vec3> const* position)
 {
diff --git a/src/polymesh/algorithms/subdivision/sqrt3.hh b/src/polymesh/algorithms/subdivision/sqrt3.hh
index 920f271b60d0c6fa8704567fbc6f09d893d42859..2985f2cfe5869a5a6a3d4e2e2dd02c6c554d4c0e 100644
--- a/src/polymesh/algorithms/subdivision/sqrt3.hh
+++ b/src/polymesh/algorithms/subdivision/sqrt3.hh
@@ -9,7 +9,7 @@ namespace polymesh
 template <class VertexF>
 void subdivide_sqrt3(Mesh& m, VertexF&& vf);
 
-/// ======== IMPLEMENTATION ========
+// ======== IMPLEMENTATION ========
 
 template <class VertexF>
 void subdivide_sqrt3(Mesh& m, VertexF&& vf)
diff --git a/src/polymesh/algorithms/topology.hh b/src/polymesh/algorithms/topology.hh
index 16e308b467c6e0e8b6ec07a1d5741a4ef4ff4317..9be358bafd6f3b6ea410bdebc38240e92b5373e9 100644
--- a/src/polymesh/algorithms/topology.hh
+++ b/src/polymesh/algorithms/topology.hh
@@ -10,7 +10,7 @@ namespace polymesh
 /// (i.e. the last face that would be visited in a BFS)
 face_handle farthest_face(face_handle f);
 
-/// ======== IMPLEMENTATION ========
+// ======== IMPLEMENTATION ========
 
 inline face_handle farthest_face(face_handle f)
 {
diff --git a/src/polymesh/algorithms/tracing.hh b/src/polymesh/algorithms/tracing.hh
index af8a606e9f7ae040815ab6de0b752323169945de..d757a17af723e8927a6763e2d8da7a5370080c15 100644
--- a/src/polymesh/algorithms/tracing.hh
+++ b/src/polymesh/algorithms/tracing.hh
@@ -16,7 +16,7 @@ namespace polymesh
 template <class Scalar = float, class EdgeLengthF>
 std::pair<halfedge_handle, float> trace_step(halfedge_handle h, EdgeLengthF&& edge_length, Scalar x, Scalar d1_sqr, Scalar d2_sqr);
 
-/// ======== IMPLEMENTATION ========
+// ======== IMPLEMENTATION ========
 
 template <class Scalar, class EdgeLengthF>
 std::pair<halfedge_handle, float> trace_step(halfedge_handle h, EdgeLengthF&& edge_length, Scalar x, Scalar d1_sqr, Scalar d2_sqr)
diff --git a/src/polymesh/attributes/flags.hh b/src/polymesh/attributes/flags.hh
index e1fa07e90bf18c88fdeb8b99f223e4ccf78e0e2a..ea87558e307f20d953b9b73471dd4c082415ccf3 100644
--- a/src/polymesh/attributes/flags.hh
+++ b/src/polymesh/attributes/flags.hh
@@ -106,7 +106,7 @@ private:
     attribute<enum_t> entries;
 };
 
-/// ======== IMPLEMENTATION ========
+// ======== IMPLEMENTATION ========
 
 template <class enum_t, class mesh_ptr, class tag, class iterator>
 flags<enum_t, tag> make_flags(smart_collection<mesh_ptr, tag, iterator> const& c, enum_t initial_value)
diff --git a/src/polymesh/attributes/partitioning.hh b/src/polymesh/attributes/partitioning.hh
index 8b03ec6a5569535fda2add59f5b1febd53fa77bc..e38947261f3f6aefbe2fe7693ed7ca99a985abee 100644
--- a/src/polymesh/attributes/partitioning.hh
+++ b/src/polymesh/attributes/partitioning.hh
@@ -84,7 +84,7 @@ private:
     int partitions;
 };
 
-/// ======== IMPLEMENTATION ========
+// ======== IMPLEMENTATION ========
 
 template <class mesh_ptr, class tag, class iterator>
 partitioning<tag> make_partitioning(smart_collection<mesh_ptr, tag, iterator> const& c)
diff --git a/src/polymesh/detail/permutation.hh b/src/polymesh/detail/permutation.hh
index 954cf6c9214bb10ef64d656bd982899f60e2c4d2..ae880594843e5cd5e801bfa0c193e3a14ede0baa 100644
--- a/src/polymesh/detail/permutation.hh
+++ b/src/polymesh/detail/permutation.hh
@@ -19,7 +19,7 @@ bool is_valid_permutation(std::vector<int> const& p);
 /// p[curr_idx] = new_idx
 std::vector<std::pair<int, int>> transpositions_of(std::vector<int> const& p);
 
-/// ======== IMPLEMENTATION ========
+// ======== IMPLEMENTATION ========
 
 inline bool is_valid_permutation(std::vector<int> const& p)
 {
diff --git a/src/polymesh/objects/cone.hh b/src/polymesh/objects/cone.hh
index 9d4e76495d0605114f677777b4500258b047aa7c..4e17284dc8d70d644be19a851a8af275dbc5d1e0 100644
--- a/src/polymesh/objects/cone.hh
+++ b/src/polymesh/objects/cone.hh
@@ -13,7 +13,7 @@ namespace polymesh::objects
 template <class coneF>
 auto add_cone(Mesh& m, coneF&& qf, int segments, bool closed = true) -> decltype(qf(vertex_handle{}, float{}, float{}), vertex_handle{});
 
-/// ======== IMPLEMENTATION ========
+// ======== IMPLEMENTATION ========
 
 template <class coneF>
 auto add_cone(Mesh& m, coneF&& qf, int segments, bool closed) -> decltype(qf(vertex_handle{}, float{}, float{}), vertex_handle{})
diff --git a/src/polymesh/objects/cube.hh b/src/polymesh/objects/cube.hh
index c82b615accd8e3464178192a85bbd460e3f76558..11cc36c2e5bd1196f34f53806774a88b33f3a81d 100644
--- a/src/polymesh/objects/cube.hh
+++ b/src/polymesh/objects/cube.hh
@@ -18,7 +18,7 @@ auto add_cube(Mesh& m, CubeF&& cf) -> decltype(cf(vertex_handle{}, int{}, int{},
 template <class Pos3>
 auto add_cube(Mesh& m, vertex_attribute<Pos3>& pos) -> vertex_handle;
 
-/// ======== IMPLEMENTATION ========
+// ======== IMPLEMENTATION ========
 
 template <class Pos3>
 auto add_cube(Mesh& m, vertex_attribute<Pos3>& pos) -> vertex_handle
diff --git a/src/polymesh/objects/cylinder.hh b/src/polymesh/objects/cylinder.hh
index 7ef8c05aab09d248000514cc68a6f58bebd0d5b8..7e576b6cb9a717c3f0184d84deb1d45fd037e168 100644
--- a/src/polymesh/objects/cylinder.hh
+++ b/src/polymesh/objects/cylinder.hh
@@ -13,7 +13,7 @@ namespace polymesh::objects
 template <class CylinderF>
 auto add_cylinder(Mesh& m, CylinderF&& qf, int segments, bool closed = true) -> decltype(qf(vertex_handle{}, float{}, float{}), vertex_handle{});
 
-/// ======== IMPLEMENTATION ========
+// ======== IMPLEMENTATION ========
 
 template <class CylinderF>
 auto add_cylinder(Mesh& m, CylinderF&& qf, int segments, bool closed) -> decltype(qf(vertex_handle{}, float{}, float{}), vertex_handle{})
diff --git a/src/polymesh/objects/quad.hh b/src/polymesh/objects/quad.hh
index 7b381d64b3415e0da20a5c4c5be5be335133ea4d..229eeab9a0f4cea066c7c9ee5e9caa658e3c01b4 100644
--- a/src/polymesh/objects/quad.hh
+++ b/src/polymesh/objects/quad.hh
@@ -14,7 +14,7 @@ namespace objects
 template <class QuadF>
 auto add_quad(Mesh& m, QuadF&& qf, int w = 1, int h = 1) -> decltype(qf(vertex_handle{}, float{}, float{}), vertex_handle{});
 
-/// ======== IMPLEMENTATION ========
+// ======== IMPLEMENTATION ========
 
 template <class QuadF>
 auto add_quad(Mesh& m, QuadF&& qf, int w, int h) -> decltype(qf(vertex_handle{}, float{}, float{}), vertex_handle{})
diff --git a/src/polymesh/objects/uv_sphere.hh b/src/polymesh/objects/uv_sphere.hh
index 894870db31d38a193bb5af2cf78229920bf8a8ef..a44b2d9e2b3e47ae57d43a40a4392201a5118eb2 100644
--- a/src/polymesh/objects/uv_sphere.hh
+++ b/src/polymesh/objects/uv_sphere.hh
@@ -13,7 +13,7 @@ namespace polymesh::objects
 template <class SphereF>
 auto add_uv_sphere(Mesh& m, SphereF&& qf, int cnt_longitude, int cnt_latitude) -> decltype(qf(vertex_handle{}, float{}, float{}), vertex_handle{});
 
-/// ======== IMPLEMENTATION ========
+// ======== IMPLEMENTATION ========
 
 template <class SphereF>
 auto add_uv_sphere(Mesh& m, SphereF&& qf, int cnt_longitude, int cnt_latitude) -> decltype(qf(vertex_handle{}, float{}, float{}), vertex_handle{})
diff --git a/src/polymesh/properties.hh b/src/polymesh/properties.hh
index 8edc8a809aa8233d604b88cea59630f27945c794..c1090fb4bec23e86c2cc5effa9dfde214e6e775d 100644
--- a/src/polymesh/properties.hh
+++ b/src/polymesh/properties.hh
@@ -143,6 +143,12 @@ typename field3<Pos3>::vec_t triangle_normal_unorm(face_handle f, vertex_attribu
 template <class Pos3>
 Pos3 bary_interpolate(face_handle f, Pos3 bary, vertex_attribute<Pos3> const& position);
 
+/// calculates the barycentric coordinates of a given point p within a face f
+/// NOTE: asserts that f is triangular
+/// NOTE: also works for other points in the same plane as f
+template <class Pos3, class Scalar = typename field3<Pos3>::scalar_t>
+Pos3 barycoords_of(face_handle f, vertex_attribute<Pos3> const& positions, Pos3 p);
+
 /// returns the length of an edge
 template <class Pos3, class Scalar = typename field3<Pos3>::scalar_t>
 Scalar edge_length(edge_handle e, vertex_attribute<Pos3> const& position);
@@ -229,7 +235,7 @@ bool is_delaunay(edge_handle e, vertex_attribute<Pos3> const& position);
 template <class Pos3>
 bool can_collapse_without_flips(halfedge_handle h, Pos3 new_pos, vertex_attribute<Pos3> const& position);
 
-/// ======== IMPLEMENTATION ========
+// ======== IMPLEMENTATION ========
 
 inline bool is_boundary(vertex_handle v) { return v.is_boundary(); }
 inline bool is_vertex_boundary(vertex_handle v) { return v.is_boundary(); }
@@ -536,6 +542,36 @@ Pos3 bary_interpolate(face_handle f, Pos3 bary, vertex_attribute<Pos3> const& po
     return (v0 * bary[0] + v1 * bary[1] + v2 * bary[2]) / field3<Pos3>::scalar(1);
 }
 
+template <class Pos3, class Scalar>
+Pos3 barycoords_of(face_handle f, vertex_attribute<Pos3> const& positions, Pos3 p)
+{
+    POLYMESH_ASSERT(f.vertices().size() == 3 && "only supports triangles");
+
+    auto ps = f.vertices().to_array<3>(positions);
+
+    auto e10 = ps[1] - ps[0];
+    auto e21 = ps[2] - ps[1];
+
+    auto n = field3<Pos3>::cross(e10, e21);
+
+    auto signed_area = [&](Pos3 const& v0, Pos3 const& v1, Pos3 const& v2) {
+        auto d1 = v1 - v0;
+        auto d2 = v2 - v0;
+
+        auto a = field3<Pos3>::cross(d1, d2);
+
+        return field3<Pos3>::dot(a, n);
+    };
+
+    auto a = signed_area(ps[0], ps[1], ps[2]);
+    auto a0 = signed_area(p, ps[1], ps[2]);
+    auto a1 = signed_area(p, ps[2], ps[0]);
+    auto a2 = signed_area(p, ps[0], ps[1]);
+
+    auto inv_a = Scalar(1) / a;
+    return field3<Pos3>::make_pos(a0 * inv_a, a1 * inv_a, a2 * inv_a);
+}
+
 template <class Pos3, class Scalar>
 vertex_attribute<Scalar> vertex_voronoi_areas(vertex_attribute<Pos3> const& position)
 {