diff --git a/CMakeLists.txt b/CMakeLists.txt
index 3a2d3d8..59279d7 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -239,3 +239,4 @@ add_nanobind_module(_subdivision src/subdivision.cpp)
 add_nanobind_module(_polylines src/polylines.cpp)
 
 add_nanobind_module(_geodesics src/geodesics.cpp)
+add_nanobind_module(_isolines src/isolines.cpp)
diff --git a/docs/api/compas_cgal.geodesics.md b/docs/api/compas_cgal.geodesics.md
new file mode 100644
index 0000000..00c64ff
--- /dev/null
+++ b/docs/api/compas_cgal.geodesics.md
@@ -0,0 +1 @@
+# ::: compas_cgal.geodesics
diff --git a/docs/api/compas_cgal.isolines.md b/docs/api/compas_cgal.isolines.md
new file mode 100644
index 0000000..7d97915
--- /dev/null
+++ b/docs/api/compas_cgal.isolines.md
@@ -0,0 +1 @@
+# ::: compas_cgal.isolines
diff --git a/docs/api/compas_cgal.polylines.md b/docs/api/compas_cgal.polylines.md
new file mode 100644
index 0000000..ee3d158
--- /dev/null
+++ b/docs/api/compas_cgal.polylines.md
@@ -0,0 +1 @@
+# ::: compas_cgal.polylines
diff --git a/docs/examples/example_isolines.md b/docs/examples/example_isolines.md
new file mode 100644
index 0000000..64a6983
--- /dev/null
+++ b/docs/examples/example_isolines.md
@@ -0,0 +1,16 @@
+# Isolines from Scalar Fields
+
+This example demonstrates isoline extraction from arbitrary vertex scalar fields using COMPAS CGAL.
+
+The visualization shows an elephant mesh with isolines extracted from geodesic distances computed from five source points: the four feet and the snout.
+
+Key Features:
+
+* Generic isoline extraction from any vertex scalar field via ``isolines``
+* Support for explicit isovalues or automatic even spacing
+* Adaptive resampling for smoother output
+* Optional Laplacian smoothing
+
+```python
+---8<--- "docs/examples/example_isolines.py"
+```
diff --git a/docs/examples/example_isolines.py b/docs/examples/example_isolines.py
new file mode 100644
index 0000000..51dc525
--- /dev/null
+++ b/docs/examples/example_isolines.py
@@ -0,0 +1,72 @@
+from pathlib import Path
+
+import numpy as np
+from compas.colors import Color
+from compas.datastructures import Mesh
+from compas.geometry import Polyline
+from compas_viewer import Viewer
+
+from compas_cgal.geodesics import heat_geodesic_distances
+from compas_cgal.isolines import isolines
+from compas_cgal.subdivision import mesh_subdivide_loop
+
+# =============================================================================
+# Load mesh and subdivide
+# =============================================================================
+
+FILE = Path(__file__).parent.parent.parent / "data" / "elephant.off"
+mesh = Mesh.from_off(FILE)
+mesh.quads_to_triangles()
+
+V, F = mesh_subdivide_loop(mesh.to_vertices_and_faces(), k=2) # k=4 for proper smoothness
+mesh = Mesh.from_vertices_and_faces(V.tolist(), F.tolist())
+
+# =============================================================================
+# Find source vertices: 4 feet and snout
+# =============================================================================
+
+V = np.array(mesh.vertices_attributes("xyz"))
+
+# feet: find lowest Y vertex in each XZ quadrant
+y_min = V[:, 1].min()
+low_verts = np.where(V[:, 1] < y_min + 0.15)[0]
+
+feet = []
+for sx in [-1, 1]:  # back/front (X)
+    for sz in [-1, 1]:  # left/right (Z)
+        mask = (np.sign(V[low_verts, 0]) == sx) & (np.sign(V[low_verts, 2]) == sz)
+        candidates = low_verts[mask]
+        if len(candidates):
+            foot = candidates[np.argmin(V[candidates, 1])]
+            feet.append(int(foot))
+
+# snout: max X (trunk tip)
+snout = int(np.argmax(V[:, 0]))
+
+sources = feet + [snout]
+
+# =============================================================================
+# Compute scalar field and extract isolines
+# =============================================================================
+
+vf = mesh.to_vertices_and_faces()
+distances = heat_geodesic_distances(vf, sources)
+
+for key, d in zip(mesh.vertices(), distances):
+    mesh.vertex_attribute(key, "distance", d)
+
+polylines = isolines(mesh, "distance", n=300, smoothing=0, resample=False)
+
+# =============================================================================
+# Viz
+# =============================================================================
+
+viewer = Viewer()
+
+viewer.scene.add(mesh, show_lines=False)
+
+for pts in polylines:
+    points = [pts[i].tolist() for i in range(len(pts))]
+    viewer.scene.add(Polyline(points), linecolor=Color.red(), lineswidth=3)
+
+viewer.show()
diff --git a/mkdocs.yml b/mkdocs.yml
index 160b44c..01ef0e1 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -142,6 +142,7 @@ nav:
       - Booleans: examples/example_booleans.md
       - Geodesics: examples/example_geodesics.md
       - Intersections: examples/example_intersections.md
+      - Isolines: examples/example_isolines.md
       - Meshing: examples/example_meshing.md
       - Meshing Dual: examples/example_meshing_dual.md
       - Projection: examples/example_projection.md
@@ -162,9 +163,12 @@ nav:
       - Triangulation: examples/example_triangulation.md
   - API Reference:
       - compas_cgal.booleans: api/compas_cgal.booleans.md
+      - compas_cgal.geodesics: api/compas_cgal.geodesics.md
       - compas_cgal.intersections: api/compas_cgal.intersections.md
+      - compas_cgal.isolines: api/compas_cgal.isolines.md
       - compas_cgal.measure: api/compas_cgal.measure.md
       - compas_cgal.meshing: api/compas_cgal.meshing.md
+      - compas_cgal.polylines: api/compas_cgal.polylines.md
       - compas_cgal.projection: api/compas_cgal.projection.md
       - compas_cgal.reconstruction: api/compas_cgal.reconstruction.md
       - compas_cgal.skeletonization: api/compas_cgal.skeletonization.md
diff --git a/requirements-dev.txt b/requirements-dev.txt
index 8fdb1bc..3383202 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -2,6 +2,7 @@ bump-my-version
 cibuildwheel
 compas_invocations2
 invoke >=0.14
+-r requirements-docs.txt
 nanobind >=1.3.2
 pytest >=7.0.0
 pytest-cov >=4.0.0
diff --git a/src/compas_cgal/__init__.py b/src/compas_cgal/__init__.py
index fb37ef4..4e5c6ce 100644
--- a/src/compas_cgal/__init__.py
+++ b/src/compas_cgal/__init__.py
@@ -18,6 +18,7 @@
 __all_plugins__ = [
     "compas_cgal.booleans",
     "compas_cgal.intersections",
+    "compas_cgal.isolines",
     "compas_cgal.meshing",
     "compas_cgal.measure",
     "compas_cgal.reconstruction",
diff --git a/src/compas_cgal/isolines.py b/src/compas_cgal/isolines.py
new file mode 100644
index 0000000..94467ce
--- /dev/null
+++ b/src/compas_cgal/isolines.py
@@ -0,0 +1,169 @@
+"""Isoline extraction from vertex scalar fields using CGAL."""
+
+from typing import List
+
+import numpy as np
+from numpy.typing import NDArray
+
+from compas.datastructures import Mesh
+from compas_cgal._isolines import isolines as _isolines
+from compas_cgal.types import PolylinesNumpy
+
+__all__ = ["isolines"]
+
+
+def _smooth_polyline(pts: NDArray, iterations: int = 1) -> NDArray:
+    """Apply Laplacian smoothing to a polyline, keeping endpoints fixed."""
+    if len(pts) < 3:
+        return pts
+    smoothed = pts.copy()
+    for _ in range(iterations):
+        new_pts = smoothed.copy()
+        for i in range(1, len(pts) - 1):
+            new_pts[i] = (smoothed[i - 1] + smoothed[i + 1]) / 2
+        smoothed = new_pts
+    return smoothed
+
+
+def _segment_lengths(pts: NDArray) -> NDArray:
+    """Compute segment lengths for a polyline."""
+    return np.linalg.norm(np.diff(pts, axis=0), axis=1)
+
+
+def _resample_polyline_adaptive(pts: NDArray, threshold: float = 2.0) -> NDArray:
+    """Resample only segments significantly longer than median.
+
+    Parameters
+    ----------
+    pts : NDArray
+        Polyline points (Nx3).
+    threshold : float
+        Segments longer than threshold * median_length get subdivided.
+
+    Returns
+    -------
+    NDArray
+        Resampled polyline.
+    """
+    if len(pts) < 3:
+        return pts
+
+    lengths = _segment_lengths(pts)
+    median_len = np.median(lengths)
+
+    if median_len == 0:
+        return pts
+
+    result = [pts[0]]
+    for i in range(len(pts) - 1):
+        p0, p1 = pts[i], pts[i + 1]
+        seg_len = lengths[i]
+
+        # subdivide long segments proportionally
+        if seg_len > threshold * median_len:
+            n_subdivs = int(np.ceil(seg_len / median_len))
+            for j in range(1, n_subdivs):
+                t = j / n_subdivs
+                result.append(p0 * (1 - t) + p1 * t)
+
+        result.append(p1)
+
+    return np.array(result)
+
+
+def _resample_polyline(pts: NDArray, factor: int) -> NDArray:
+    """Resample a polyline by inserting interpolated points between each pair."""
+    if len(pts) < 2 or factor < 2:
+        return pts
+    result = [pts[0]]
+    for i in range(len(pts) - 1):
+        p0, p1 = pts[i], pts[i + 1]
+        for j in range(1, factor):
+            t = j / factor
+            result.append(p0 * (1 - t) + p1 * t)
+        result.append(p1)
+    return np.array(result)
+
+
+def isolines(
+    mesh: Mesh,
+    scalars: str,
+    isovalues: List[float] | None = None,
+    n: int | None = None,
+    resample: int | bool = True,
+    smoothing: int = 0,
+) -> PolylinesNumpy:
+    """Extract isoline polylines from vertex scalar field.
+
+    Uses CGAL's refine_mesh_at_isolevel to extract isolines from a scalar
+    field defined at mesh vertices.
+
+    Parameters
+    ----------
+    mesh : :class:`compas.datastructures.Mesh`
+        A triangulated mesh.
+    scalars : str
+        Name of the vertex attribute containing scalar values.
+    isovalues : List[float], optional
+        Explicit isovalue thresholds for isoline extraction.
+    n : int, optional
+        Number of evenly spaced isovalues between scalar min and max.
+        The isovalues will exclude the endpoints.
+    resample : int or bool, optional
+        Polyline resampling mode. If True (default), adaptively resample
+        segments longer than 2x median length. If int > 1, uniformly
+        subdivide each segment into that many parts. If False or 1, disable.
+    smoothing : int, optional
+        Number of Laplacian smoothing iterations to apply to polylines.
+        Default is 0 (no smoothing).
+
+    Returns
+    -------
+    :attr:`compas_cgal.types.PolylinesNumpy`
+        List of polyline segments as Nx3 arrays of points.
+
+    Raises
+    ------
+    ValueError
+        If neither or both of `isovalues` and `n` are provided.
+
+    Examples
+    --------
+    >>> from compas.datastructures import Mesh
+    >>> from compas.geometry import Sphere
+    >>> from compas_cgal.geodesics import heat_geodesic_distances
+    >>> from compas_cgal.isolines import isolines
+    >>> sphere = Sphere(1.0)
+    >>> mesh = Mesh.from_shape(sphere, u=32, v=32)
+    >>> mesh.quads_to_triangles()
+    >>> vf = mesh.to_vertices_and_faces()
+    >>> distances = heat_geodesic_distances(vf, [0])
+    >>> for key, d in zip(mesh.vertices(), distances):
+    ...     mesh.vertex_attribute(key, "distance", d)
+    >>> polylines = isolines(mesh, "distance", n=5)
+
+    """
+    V = np.asarray(mesh.vertices_attributes("xyz"), dtype=np.float64, order="C")
+    F = np.asarray([mesh.face_vertices(f) for f in mesh.faces()], dtype=np.int32, order="C")
+    scalar_values = np.asarray(mesh.vertices_attribute(scalars), dtype=np.float64, order="C").reshape(-1, 1)
+
+    if isovalues is None and n is None:
+        raise ValueError("Either 'isovalues' or 'n' must be provided.")
+    if isovalues is not None and n is not None:
+        raise ValueError("Provide exactly one of 'isovalues' or 'n' (not both).")
+
+    if n is not None:
+        smin, smax = float(scalar_values.min()), float(scalar_values.max())
+        isovalues = np.linspace(smin, smax, n + 2)[1:-1].tolist()
+
+    polylines = list(_isolines(V, F, scalar_values, isovalues, 0))
+
+    if resample is True:
+        polylines = [_resample_polyline_adaptive(pl) for pl in polylines]
+    elif isinstance(resample, int) and resample > 1:
+        polylines = [_resample_polyline(pl, resample) for pl in polylines]
+
+    if smoothing > 0:
+        polylines = [_smooth_polyline(pl, smoothing) for pl in polylines]
+
+    return polylines
diff --git a/src/isolines.cpp b/src/isolines.cpp
new file mode 100644
index 0000000..6a2d674
--- /dev/null
+++ b/src/isolines.cpp
@@ -0,0 +1,175 @@
+#include "compas.h"
+
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Polygon_mesh_processing/refine_mesh_at_isolevel.h>
+#include <CGAL/boost/graph/split_graph_into_polylines.h>
+
+using vertex_descriptor = boost::graph_traits<compas::Mesh>::vertex_descriptor;
+using edge_descriptor = boost::graph_traits<compas::Mesh>::edge_descriptor;
+
+
+// Extract a single isoline with optional local refinement
+std::vector<std::vector<compas::Kernel::Point_3>>
+extract_single_isoline(
+    Eigen::Ref<const compas::RowMatrixXd> vertices,
+    Eigen::Ref<const compas::RowMatrixXi> faces,
+    Eigen::Ref<const compas::RowMatrixXd> scalars,
+    double isovalue,
+    double iso_spacing,
+    int refine)
+{
+    typedef boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS,
+        vertex_descriptor> IsolineGraph;
+    typedef boost::graph_traits<IsolineGraph>::vertex_descriptor ig_vertex;
+
+    // Fresh mesh for this isoline
+    compas::Mesh mesh = compas::mesh_from_vertices_and_faces(vertices, faces);
+    std::size_t num_verts_before = mesh.number_of_vertices();
+
+    // Create and populate scalar property map
+    auto scalar_pmap = mesh.add_property_map<vertex_descriptor, double>("v:scalar", 0.0).first;
+    for (vertex_descriptor vd : mesh.vertices()) {
+        scalar_pmap[vd] = scalars(vd.idx(), 0);
+    }
+
+    // Edge constraint map (only used for tracking, not refinement passes)
+    auto ecm = mesh.add_property_map<edge_descriptor, bool>("e:is_constrained", false).first;
+
+    // Pre-refine with values very close to the target to subdivide crossing edges
+    if (refine > 0 && iso_spacing > 0) {
+        // Use small offsets (fraction of spacing) to split edges that span the isovalue
+        double base_offset = iso_spacing * 0.1;  // 10% of spacing
+        for (int i = 1; i <= refine; ++i) {
+            double offset = base_offset * i / refine;
+
+            // Refine below - new vertices get the isovalue they were inserted at
+            CGAL::Polygon_mesh_processing::refine_mesh_at_isolevel(
+                mesh, scalar_pmap, isovalue - offset,
+                CGAL::parameters::edge_is_constrained_map(ecm));
+
+            // Update scalar values for new vertices
+            for (std::size_t idx = num_verts_before; idx < mesh.number_of_vertices(); ++idx) {
+                scalar_pmap[vertex_descriptor(idx)] = isovalue - offset;
+            }
+            num_verts_before = mesh.number_of_vertices();
+
+            // Refine above
+            CGAL::Polygon_mesh_processing::refine_mesh_at_isolevel(
+                mesh, scalar_pmap, isovalue + offset,
+                CGAL::parameters::edge_is_constrained_map(ecm));
+
+            for (std::size_t idx = num_verts_before; idx < mesh.number_of_vertices(); ++idx) {
+                scalar_pmap[vertex_descriptor(idx)] = isovalue + offset;
+            }
+            num_verts_before = mesh.number_of_vertices();
+        }
+
+        // Clear constraints from pre-refinement
+        for (auto e : mesh.edges()) {
+            put(ecm, e, false);
+        }
+    }
+
+    // Extract the actual isoline
+    CGAL::Polygon_mesh_processing::refine_mesh_at_isolevel(
+        mesh, scalar_pmap, isovalue,
+        CGAL::parameters::edge_is_constrained_map(ecm));
+
+    // Build graph from constrained edges
+    IsolineGraph iso_graph;
+    std::map<vertex_descriptor, ig_vertex> v_map;
+
+    for (auto e : mesh.edges()) {
+        if (get(ecm, e)) {
+            auto h = mesh.halfedge(e);
+            auto src = mesh.source(h);
+            auto tgt = mesh.target(h);
+
+            if (v_map.find(src) == v_map.end()) {
+                v_map[src] = boost::add_vertex(src, iso_graph);
+            }
+            if (v_map.find(tgt) == v_map.end()) {
+                v_map[tgt] = boost::add_vertex(tgt, iso_graph);
+            }
+            boost::add_edge(v_map[src], v_map[tgt], iso_graph);
+        }
+    }
+
+    // Visitor to collect polylines
+    struct PolylineVisitor {
+        const compas::Mesh& mesh;
+        std::vector<std::vector<compas::Kernel::Point_3>>& polylines;
+        std::vector<compas::Kernel::Point_3> current;
+
+        PolylineVisitor(const compas::Mesh& m, std::vector<std::vector<compas::Kernel::Point_3>>& p)
+            : mesh(m), polylines(p) {}
+
+        void start_new_polyline() { current.clear(); }
+        void add_node(ig_vertex v) {
+            vertex_descriptor vd = boost::get(boost::vertex_bundle, *graph_ptr)[v];
+            current.push_back(mesh.point(vd));
+        }
+        void end_polyline() {
+            if (current.size() >= 2) polylines.push_back(current);
+        }
+        const IsolineGraph* graph_ptr;
+    };
+
+    std::vector<std::vector<compas::Kernel::Point_3>> polyline_points;
+    PolylineVisitor visitor(mesh, polyline_points);
+    visitor.graph_ptr = &iso_graph;
+
+    CGAL::split_graph_into_polylines(iso_graph, visitor, CGAL::internal::IsTerminalDefault());
+
+    return polyline_points;
+}
+
+
+std::vector<compas::RowMatrixXd>
+pmp_isolines(
+    Eigen::Ref<const compas::RowMatrixXd> vertices,
+    Eigen::Ref<const compas::RowMatrixXi> faces,
+    Eigen::Ref<const compas::RowMatrixXd> scalars,
+    const std::vector<double>& isovalues,
+    int refine)
+{
+    // Compute spacing between isovalues for refinement offsets
+    double iso_spacing = 0.0;
+    if (isovalues.size() > 1) {
+        iso_spacing = std::abs(isovalues[1] - isovalues[0]);
+    } else {
+        double smin = scalars.minCoeff();
+        double smax = scalars.maxCoeff();
+        iso_spacing = (smax - smin) / 100.0;
+    }
+
+    std::vector<compas::RowMatrixXd> all_polylines;
+
+    // Process each isovalue independently
+    for (double isovalue : isovalues) {
+        auto polyline_points = extract_single_isoline(
+            vertices, faces, scalars, isovalue, iso_spacing, refine);
+
+        // Convert to output format
+        for (const auto& pl : polyline_points) {
+            compas::RowMatrixXd pts(pl.size(), 3);
+            for (std::size_t i = 0; i < pl.size(); ++i) {
+                pts(i, 0) = CGAL::to_double(pl[i].x());
+                pts(i, 1) = CGAL::to_double(pl[i].y());
+                pts(i, 2) = CGAL::to_double(pl[i].z());
+            }
+            all_polylines.push_back(std::move(pts));
+        }
+    }
+
+    return all_polylines;
+}
+
+
+NB_MODULE(_isolines, m) {
+    m.def(
+        "isolines",
+        &pmp_isolines,
+        "Extract isoline polylines from vertex scalar field.",
+        "vertices"_a, "faces"_a, "scalars"_a, "isovalues"_a, "refine"_a = 0);
+}