diff --git a/crates/fj-kernel/src/algorithms/approx/curve.rs b/crates/fj-kernel/src/algorithms/approx/curve.rs
index 9610e278d..708bec430 100644
--- a/crates/fj-kernel/src/algorithms/approx/curve.rs
+++ b/crates/fj-kernel/src/algorithms/approx/curve.rs
@@ -15,28 +15,31 @@ use fj_math::{Circle, Point, Scalar};
 
 use crate::objects::{Curve, CurveKind, GlobalCurve, Vertex};
 
-use super::{Approx, Tolerance};
+use super::{Approx, ApproxPoint, Tolerance};
 
 impl Approx for (&Curve, RangeOnCurve) {
-    type Approximation = Vec<(Point<2>, Point<3>)>;
+    type Approximation = CurveApprox;
 
     fn approx(self, tolerance: Tolerance) -> Self::Approximation {
         let (curve, range) = self;
 
-        (curve.global_form(), range)
+        let points = (curve.global_form(), range)
             .approx(tolerance)
             .into_iter()
-            .map(|(point_curve, point_global)| {
+            .map(|point| {
                 let point_surface =
-                    curve.kind().point_from_curve_coords(point_curve);
-                (point_surface, point_global)
+                    curve.kind().point_from_curve_coords(point.local_form);
+                ApproxPoint::new(point_surface, point.global_form)
+                    .with_source((*curve, point.local_form))
             })
-            .collect()
+            .collect();
+
+        CurveApprox { points }
     }
 }
 
 impl Approx for (&GlobalCurve, RangeOnCurve) {
-    type Approximation = Vec<(Point<1>, Point<3>)>;
+    type Approximation = Vec<ApproxPoint<1>>;
 
     fn approx(self, tolerance: Tolerance) -> Self::Approximation {
         let (curve, range) = self;
@@ -62,7 +65,7 @@ fn approx_circle(
     circle: &Circle<3>,
     range: impl Into<RangeOnCurve>,
     tolerance: Tolerance,
-    points: &mut Vec<(Point<1>, Point<3>)>,
+    points: &mut Vec<ApproxPoint<1>>,
 ) {
     let radius = circle.a().magnitude();
     let range = range.into();
@@ -82,7 +85,7 @@ fn approx_circle(
         let point_curve = Point::from([angle]);
         let point_global = circle.point_from_circle_coords(point_curve);
 
-        points.push((point_curve, point_global));
+        points.push(ApproxPoint::new(point_curve, point_global));
     }
 }
 
@@ -137,6 +140,13 @@ impl RangeOnCurve {
     }
 }
 
+/// An approximation of a [`Curve`]
+#[derive(Debug, Eq, PartialEq, Hash, Ord, PartialOrd)]
+pub struct CurveApprox {
+    /// The points that approximate the curve
+    pub points: Vec<ApproxPoint<2>>,
+}
+
 #[cfg(test)]
 mod tests {
     use fj_math::Scalar;
diff --git a/crates/fj-kernel/src/algorithms/approx/cycle.rs b/crates/fj-kernel/src/algorithms/approx/cycle.rs
index 29f2ea94b..9610202d7 100644
--- a/crates/fj-kernel/src/algorithms/approx/cycle.rs
+++ b/crates/fj-kernel/src/algorithms/approx/cycle.rs
@@ -2,52 +2,55 @@
 //!
 //! See [`CycleApprox`].
 
-use fj_math::{Point, Segment};
+use fj_math::Segment;
 
 use crate::objects::Cycle;
 
-use super::{Approx, Tolerance};
+use super::{edge::EdgeApprox, Approx, ApproxPoint, Tolerance};
 
 impl Approx for &Cycle {
     type Approximation = CycleApprox;
 
     fn approx(self, tolerance: Tolerance) -> Self::Approximation {
-        let mut points = Vec::new();
-
-        for edge in self.edges() {
-            let edge_points = edge.approx(tolerance);
-            points.extend(edge_points);
-        }
-
-        if let Some(&point) = points.first() {
-            points.push(point);
-        }
-
-        CycleApprox { points }
+        let edges = self.edges().map(|edge| edge.approx(tolerance)).collect();
+        CycleApprox { edges }
     }
 }
 
 /// An approximation of a [`Cycle`]
-#[derive(Debug, Eq, PartialEq, Hash)]
+#[derive(Debug, Eq, PartialEq, Hash, Ord, PartialOrd)]
 pub struct CycleApprox {
-    /// The points that approximate the cycle
-    pub points: Vec<(Point<2>, Point<3>)>,
+    /// The approximated edges that make up the approximated cycle
+    pub edges: Vec<EdgeApprox>,
 }
 
 impl CycleApprox {
+    /// Compute the points that approximate the cycle
+    pub fn points(&self) -> Vec<ApproxPoint<2>> {
+        let mut points = Vec::new();
+
+        for edge_approx in &self.edges {
+            points.extend(edge_approx.points());
+        }
+
+        if let Some(point) = points.first() {
+            points.push(point.clone());
+        }
+
+        points
+    }
+
     /// Construct the segments that approximate the cycle
     pub fn segments(&self) -> Vec<Segment<3>> {
         let mut segments = Vec::new();
 
-        for segment in self.points.windows(2) {
+        for segment in self.points().windows(2) {
             // This can't panic, as we passed `2` to `windows`. Can be cleaned
             // up, once `array_windows` is stable.
-            let segment = [segment[0], segment[1]];
+            let segment = [&segment[0], &segment[1]];
 
-            segments.push(Segment::from(segment.map(|point| {
-                let (_, point_global) = point;
-                point_global
-            })));
+            segments
+                .push(Segment::from(segment.map(|point| point.global_form)));
         }
 
         segments
diff --git a/crates/fj-kernel/src/algorithms/approx/edge.rs b/crates/fj-kernel/src/algorithms/approx/edge.rs
index 3cefdc3fc..e5de04683 100644
--- a/crates/fj-kernel/src/algorithms/approx/edge.rs
+++ b/crates/fj-kernel/src/algorithms/approx/edge.rs
@@ -9,10 +9,13 @@ use fj_math::{Point, Scalar};
 
 use crate::objects::{Edge, GlobalVertex, SurfaceVertex, Vertex};
 
-use super::{curve::RangeOnCurve, Approx};
+use super::{
+    curve::{CurveApprox, RangeOnCurve},
+    Approx, ApproxPoint,
+};
 
 impl Approx for &Edge {
-    type Approximation = Vec<(Point<2>, Point<3>)>;
+    type Approximation = EdgeApprox;
 
     fn approx(self, tolerance: super::Tolerance) -> Self::Approximation {
         // The range is only used for circles right now.
@@ -71,14 +74,36 @@ impl Approx for &Edge {
 
         let range = RangeOnCurve { boundary };
 
-        let mut points = (self.curve(), range).approx(tolerance);
-        points.insert(
-            0,
-            (
-                range.start().surface_form().position(),
-                range.start().global_form().position(),
-            ),
+        let first = ApproxPoint::new(
+            range.start().surface_form().position(),
+            range.start().global_form().position(),
         );
+        let curve_approx = (self.curve(), range).approx(tolerance);
+
+        EdgeApprox {
+            first,
+            curve_approx,
+        }
+    }
+}
+
+/// An approximation of an [`Edge`]
+#[derive(Debug, Eq, PartialEq, Hash, Ord, PartialOrd)]
+pub struct EdgeApprox {
+    /// The point that approximates the first vertex of the curve
+    pub first: ApproxPoint<2>,
+
+    /// The approximation of the edge's curve
+    pub curve_approx: CurveApprox,
+}
+
+impl EdgeApprox {
+    /// Compute the points that approximate the edge
+    pub fn points(&self) -> Vec<ApproxPoint<2>> {
+        let mut points = Vec::new();
+
+        points.push(self.first.clone());
+        points.extend(self.curve_approx.points.clone());
 
         points
     }
diff --git a/crates/fj-kernel/src/algorithms/approx/face.rs b/crates/fj-kernel/src/algorithms/approx/face.rs
index 55b7b0586..d3f9819e1 100644
--- a/crates/fj-kernel/src/algorithms/approx/face.rs
+++ b/crates/fj-kernel/src/algorithms/approx/face.rs
@@ -2,13 +2,13 @@
 //!
 //! See [`FaceApprox`].
 
-use std::collections::HashSet;
+use std::collections::BTreeSet;
 
-use fj_math::Point;
+use fj_interop::mesh::Color;
 
 use crate::objects::Face;
 
-use super::{cycle::CycleApprox, Approx, Tolerance};
+use super::{cycle::CycleApprox, Approx, ApproxPoint, Tolerance};
 
 impl Approx for &Face {
     type Approximation = FaceApprox;
@@ -27,20 +27,15 @@ impl Approx for &Face {
         // would need to provide its own approximation, as the edges that bound
         // it have nothing to do with its curvature.
 
-        let mut points = HashSet::new();
         let mut exteriors = Vec::new();
-        let mut interiors = HashSet::new();
+        let mut interiors = BTreeSet::new();
 
         for cycle in self.exteriors() {
             let cycle = cycle.approx(tolerance);
-
-            points.extend(cycle.points.iter().copied());
             exteriors.push(cycle);
         }
         for cycle in self.interiors() {
             let cycle = cycle.approx(tolerance);
-
-            points.extend(cycle.points.iter().copied());
             interiors.insert(cycle);
         }
 
@@ -57,84 +52,37 @@ impl Approx for &Face {
         );
 
         FaceApprox {
-            points,
             exterior,
             interiors,
+            color: self.color(),
         }
     }
 }
 
 /// An approximation of a [`Face`]
-#[derive(Debug, Eq, PartialEq)]
+#[derive(Debug, Eq, PartialEq, Ord, PartialOrd)]
 pub struct FaceApprox {
-    /// All points that make up the approximation
-    ///
-    /// These could be actual vertices from the model, points that approximate
-    /// an edge, or points that approximate a face.
-    pub points: HashSet<(Point<2>, Point<3>)>,
-
     /// Approximation of the exterior cycle
     pub exterior: CycleApprox,
 
     /// Approximations of the interior cycles
-    pub interiors: HashSet<CycleApprox>,
+    pub interiors: BTreeSet<CycleApprox>,
+
+    /// The color of the approximated face
+    pub color: Color,
 }
 
-#[cfg(test)]
-mod tests {
-    use fj_math::{Point, Scalar};
-    use map_macro::set;
-
-    use crate::{
-        algorithms::approx::Approx,
-        objects::{Face, Surface},
-    };
-
-    use super::{CycleApprox, FaceApprox, Tolerance};
-
-    #[test]
-    fn for_face_closed() -> anyhow::Result<()> {
-        // Test a closed face, i.e. one that is completely encircled by edges.
-
-        let tolerance = Tolerance::from_scalar(Scalar::ONE)?;
-
-        let a = Point::from([0., 0.]);
-        let b = Point::from([3., 0.]);
-        let c = Point::from([3., 3.]);
-        let d = Point::from([0., 3.]);
-
-        let e = Point::from([1., 1.]);
-        let f = Point::from([2., 1.]);
-        let g = Point::from([2., 2.]);
-        let h = Point::from([1., 2.]);
-
-        let surface = Surface::xy_plane();
-        let face = Face::build(surface)
-            .polygon_from_points([a, b, c, d])
-            .with_hole([e, f, g, h]);
-
-        let a = (a, a.to_xyz());
-        let b = (b, b.to_xyz());
-        let c = (c, c.to_xyz());
-        let d = (d, d.to_xyz());
-        let e = (e, e.to_xyz());
-        let f = (f, f.to_xyz());
-        let g = (g, g.to_xyz());
-        let h = (h, h.to_xyz());
-
-        let approx = face.approx(tolerance);
-        let expected = FaceApprox {
-            points: set![a, b, c, d, e, f, g, h],
-            exterior: CycleApprox {
-                points: vec![a, b, c, d, a],
-            },
-            interiors: set![CycleApprox {
-                points: vec![e, f, g, h, e],
-            }],
-        };
-
-        assert_eq!(approx, expected);
-
-        Ok(())
+impl FaceApprox {
+    /// Compute all points that make up the approximation
+    pub fn points(&self) -> BTreeSet<ApproxPoint<2>> {
+        let mut points = BTreeSet::new();
+
+        points.extend(self.exterior.points());
+
+        for cycle_approx in &self.interiors {
+            points.extend(cycle_approx.points());
+        }
+
+        points
     }
 }
diff --git a/crates/fj-kernel/src/algorithms/approx/mod.rs b/crates/fj-kernel/src/algorithms/approx/mod.rs
index cf17774af..fc6bafcb1 100644
--- a/crates/fj-kernel/src/algorithms/approx/mod.rs
+++ b/crates/fj-kernel/src/algorithms/approx/mod.rs
@@ -6,6 +6,18 @@ pub mod edge;
 pub mod face;
 pub mod tolerance;
 
+use std::{
+    any::Any,
+    cmp::Ordering,
+    fmt::Debug,
+    hash::{Hash, Hasher},
+    rc::Rc,
+};
+
+use fj_math::Point;
+
+use crate::objects::Curve;
+
 pub use self::tolerance::{InvalidTolerance, Tolerance};
 
 /// Approximate an object
@@ -19,3 +31,72 @@ pub trait Approx {
     /// the actual object.
     fn approx(self, tolerance: Tolerance) -> Self::Approximation;
 }
+
+/// A point from an approximation, with local and global forms
+#[derive(Debug, Clone)]
+pub struct ApproxPoint<const D: usize> {
+    /// The local form of the point
+    pub local_form: Point<D>,
+
+    /// The global form of the points
+    pub global_form: Point<3>,
+
+    /// The optional source of the point
+    pub source: Option<Rc<dyn Source>>,
+}
+
+impl<const D: usize> ApproxPoint<D> {
+    /// Create an instance of `ApproxPoint`, without a source
+    pub fn new(local_form: Point<D>, global_form: Point<3>) -> Self {
+        Self {
+            local_form,
+            global_form,
+            source: None,
+        }
+    }
+
+    /// Attach a source to the point
+    pub fn with_source(self, source: impl Source) -> Self {
+        Self {
+            source: Some(Rc::new(source)),
+            ..self
+        }
+    }
+}
+
+impl<const D: usize> Eq for ApproxPoint<D> {}
+
+impl<const D: usize> PartialEq for ApproxPoint<D> {
+    fn eq(&self, other: &Self) -> bool {
+        self.local_form == other.local_form
+            && self.global_form == other.global_form
+    }
+}
+
+impl<const D: usize> Hash for ApproxPoint<D> {
+    fn hash<H: Hasher>(&self, state: &mut H) {
+        self.local_form.hash(state);
+        self.global_form.hash(state);
+    }
+}
+
+impl<const D: usize> Ord for ApproxPoint<D> {
+    fn cmp(&self, other: &Self) -> Ordering {
+        match self.local_form.cmp(&other.local_form) {
+            Ordering::Equal => {}
+            ord => return ord,
+        }
+        self.global_form.cmp(&other.global_form)
+    }
+}
+
+impl<const D: usize> PartialOrd for ApproxPoint<D> {
+    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
+        Some(self.cmp(other))
+    }
+}
+
+/// The source of an [`ApproxPoint`]
+pub trait Source: Any + Debug {}
+
+impl Source for (Curve, Point<1>) {}
diff --git a/crates/fj-kernel/src/algorithms/triangulate/mod.rs b/crates/fj-kernel/src/algorithms/triangulate/mod.rs
index 1d69b84ec..7927b9711 100644
--- a/crates/fj-kernel/src/algorithms/triangulate/mod.rs
+++ b/crates/fj-kernel/src/algorithms/triangulate/mod.rs
@@ -73,26 +73,23 @@ impl Triangulate for Face {
         let approx = self.approx(tolerance.into());
 
         let points: Vec<_> = approx
-            .points
+            .points()
             .into_iter()
-            .map(|(point_surface, point_global)| TriangulationPoint {
-                point_surface,
-                point_global,
+            .map(|point| TriangulationPoint {
+                point_surface: point.local_form,
+                point_global: point.global_form,
             })
             .collect();
         let face_as_polygon = Polygon::new(*surface)
             .with_exterior(
                 approx
                     .exterior
-                    .points
+                    .points()
                     .into_iter()
-                    .map(|(point_surface, _)| point_surface),
+                    .map(|point| point.local_form),
             )
             .with_interiors(approx.interiors.into_iter().map(|interior| {
-                interior
-                    .points
-                    .into_iter()
-                    .map(|(point_surface, _)| point_surface)
+                interior.points().into_iter().map(|point| point.local_form)
             }));
 
         let mut triangles = delaunay::triangulate(points);