diff --git a/editor/src/messages/tool/tool_messages/freehand_tool.rs b/editor/src/messages/tool/tool_messages/freehand_tool.rs
index c0d728182..3e0cae9b3 100644
--- a/editor/src/messages/tool/tool_messages/freehand_tool.rs
+++ b/editor/src/messages/tool/tool_messages/freehand_tool.rs
@@ -11,7 +11,7 @@ use graph_craft::document::NodeId;
 use graphene_core::uuid::generate_uuid;
 use graphene_core::vector::VectorModificationType;
 use graphene_core::Color;
-use graphene_std::vector::{PointId, SegmentId};
+use graphene_std::vector::{HandleId, PointId, SegmentId};
 
 use glam::DVec2;
 
@@ -175,13 +175,117 @@ impl ToolTransition for FreehandTool {
 
 #[derive(Clone, Debug, Default)]
 struct FreehandToolData {
-	extend_from_start: bool,
-	end_point: Option<(DVec2, PointId)>,
+	positions: Vec<DVec2>,
 	dragged: bool,
+	required_tangent: DVec2,
+	last_tangent: DVec2,
+
+	start: Option<(PointId, DVec2)>,
+	end: Option<(PointId, DVec2)>,
+	segment: Option<SegmentId>,
+	last_segment: Option<SegmentId>,
+
 	weight: f64,
 	layer: Option<LayerNodeIdentifier>,
 }
 
+impl FreehandToolData {
+	fn smooth(&mut self, responses: &mut VecDeque<Message>) {
+		const MAX_POSITIONS: usize = 16;
+
+		let Some(layer) = self.layer else { return };
+		let fit = if self.positions.len() < MAX_POSITIONS {
+			bezier_rs::Subpath::<PointId>::fit_cubic(&self.positions, 1, self.required_tangent, 0.)
+		} else {
+			None
+		};
+
+		let Some(bezier) = fit.and_then(|subpath| subpath.iter().next()) else {
+			// Start a new segment
+			if let Some(point) = self.end {
+				self.positions.clear();
+				self.positions.push(point.1);
+				self.start = Some(point);
+				self.end = None;
+				self.required_tangent = self.last_tangent;
+				self.last_segment = self.segment;
+				self.segment = None;
+			}
+			return;
+		};
+
+		let set_point = |point: &mut Option<(PointId, DVec2)>, position: DVec2, layer, responses: &mut VecDeque<Message>| {
+			if let Some((id, current)) = point {
+				let delta = position - *current;
+				*current += delta;
+				responses.add(GraphOperationMessage::Vector {
+					layer,
+					modification_type: VectorModificationType::ApplyPointDelta { point: *id, delta },
+				});
+				*id
+			} else {
+				let id = PointId::generate();
+				*point = Some((id, position));
+				responses.add(GraphOperationMessage::Vector {
+					layer,
+					modification_type: VectorModificationType::InsertPoint { id, position },
+				});
+				id
+			}
+		};
+
+		let start = set_point(&mut self.start, bezier.start, layer, responses);
+		let end = set_point(&mut self.end, bezier.end, layer, responses);
+		let points = [start, end];
+
+		let handles = [bezier.handle_start().map(|handle| handle - bezier.start), bezier.handle_end().map(|handle| handle - bezier.end)];
+		let segment = if let Some(segment) = self.segment {
+			responses.add(GraphOperationMessage::Vector {
+				layer,
+				modification_type: VectorModificationType::SetHandles { segment, handles },
+			});
+			segment
+		} else {
+			let id = SegmentId::generate();
+			self.segment = Some(id);
+			responses.add(GraphOperationMessage::Vector {
+				layer,
+				modification_type: VectorModificationType::InsertSegment { id, points, handles },
+			});
+			id
+		};
+
+		// Mirror handles if appropriate
+		if let Some(last_segment) = self.last_segment {
+			if last_segment != segment {
+				responses.add(GraphOperationMessage::Vector {
+					layer,
+					modification_type: VectorModificationType::SetG1Continuous {
+						handles: [HandleId::end(last_segment), HandleId::primary(segment)],
+						enabled: true,
+					},
+				});
+			}
+		}
+
+		self.last_tangent = match bezier.handles {
+			bezier_rs::BezierHandles::Linear => bezier.end - bezier.start,
+			bezier_rs::BezierHandles::Quadratic { handle } => bezier.end - handle,
+			bezier_rs::BezierHandles::Cubic { handle_end, .. } => bezier.end - handle_end,
+		}
+		.normalize_or_zero();
+	}
+
+	fn push(&mut self, document: &DocumentMessageHandler, layer: LayerNodeIdentifier, viewport: DVec2) {
+		let transform = document.metadata().transform_to_viewport(layer);
+		let position = transform.inverse().transform_point2(viewport);
+
+		if self.positions.last() != Some(&position) && position.is_finite() {
+			self.positions.push(position);
+		}
+	}
+}
+
 impl Fsm for FreehandToolFsmState {
 	type ToolData = FreehandToolData;
 	type ToolOptions = FreehandOptions;
@@ -207,48 +311,50 @@ impl Fsm for FreehandToolFsmState {
 			(FreehandToolFsmState::Ready, FreehandToolMessage::DragStart) => {
 				responses.add(DocumentMessage::StartTransaction);
 
-				tool_data.dragged = false;
-				tool_data.extend_from_start = false;
 				tool_data.weight = tool_options.line_weight;
+				tool_data.positions.clear();
+				tool_data.required_tangent = DVec2::ZERO;
+				tool_data.last_tangent = DVec2::ZERO;
+				tool_data.start = None;
+				tool_data.end = None;
+				tool_data.segment = None;
+				tool_data.last_segment = None;
+				tool_data.dragged = false;
 
-				// Extend an endpoint of the selected path
-				if let Some((layer, _, position)) = should_extend(document, input.mouse.position, crate::consts::SNAP_POINT_TOLERANCE) {
+				let layer = if let Some((layer, point, position)) = should_extend(document, input.mouse.position, crate::consts::SNAP_POINT_TOLERANCE) {
+					tool_data.positions.push(position);
+					tool_data.start = Some((point, position));
 					tool_data.layer = Some(layer);
+					layer
+				} else {
+					responses.add(DocumentMessage::DeselectAllLayers);
 
-					extend_path_with_next_segment(tool_data, position, responses);
+					let parent = document.new_layer_parent(true);
 
-					return FreehandToolFsmState::Drawing;
-				}
+					let nodes = {
+						let node_type = resolve_document_node_type("Path").expect("Path node does not exist");
+						let node = node_type.to_document_node_default_inputs([], Default::default());
 
-				responses.add(DocumentMessage::DeselectAllLayers);
+						HashMap::from([(NodeId(0), node)])
+					};
 
-				let parent = document.new_layer_parent(true);
+					let layer = graph_modification_utils::new_custom(NodeId(generate_uuid()), nodes, parent, responses);
+					tool_options.fill.apply_fill(layer, responses);
+					tool_options.stroke.apply_stroke(tool_data.weight, layer, responses);
 
-				let nodes = {
-					let node_type = resolve_document_node_type("Path").expect("Path node does not exist");
-					let node = node_type.to_document_node_default_inputs([], Default::default());
-
-					HashMap::from([(NodeId(0), node)])
+					tool_data.layer = Some(layer);
+					parent
 				};
-
-				let layer = graph_modification_utils::new_custom(NodeId(generate_uuid()), nodes, parent, responses);
-				tool_options.fill.apply_fill(layer, responses);
-				tool_options.stroke.apply_stroke(tool_data.weight, layer, responses);
-				tool_data.layer = Some(layer);
-
-				let transform = document.metadata().transform_to_viewport(layer);
-				let position = transform.inverse().transform_point2(input.mouse.position);
-
-				extend_path_with_next_segment(tool_data, position, responses);
+				tool_data.push(&document, layer, input.mouse.position);
+				tool_data.smooth(responses);
 
 				FreehandToolFsmState::Drawing
 			}
 			(FreehandToolFsmState::Drawing, FreehandToolMessage::PointerMove) => {
 				if let Some(layer) = tool_data.layer {
-					let transform = document.metadata().transform_to_viewport(layer);
-					let position = transform.inverse().transform_point2(input.mouse.position);
-
-					extend_path_with_next_segment(tool_data, position, responses);
+					tool_data.push(&document, layer, input.mouse.position);
+					tool_data.smooth(responses);
+					tool_data.dragged = true;
 				}
 
 				FreehandToolFsmState::Drawing
@@ -294,34 +400,3 @@ impl Fsm for FreehandToolFsmState {
 		responses.add(FrontendMessage::UpdateMouseCursor { cursor: MouseCursorIcon::Default });
 	}
 }
-
-fn extend_path_with_next_segment(tool_data: &mut FreehandToolData, position: DVec2, responses: &mut VecDeque<Message>) {
-	if !tool_data.end_point.map_or(true, |(last_pos, _)| position != last_pos) || !position.is_finite() {
-		return;
-	}
-
-	let Some(layer) = tool_data.layer else { return };
-
-	let id = PointId::generate();
-	responses.add(GraphOperationMessage::Vector {
-		layer,
-		modification_type: VectorModificationType::InsertPoint { id, position },
-	});
-
-	if let Some((_, previous_position)) = tool_data.end_point {
-		let next_id = SegmentId::generate();
-		let points = [previous_position, id];
-
-		responses.add(GraphOperationMessage::Vector {
-			layer,
-			modification_type: VectorModificationType::InsertSegment {
-				id: next_id,
-				points,
-				handles: [None, None],
-			},
-		});
-	}
-
-	tool_data.dragged = true;
-	tool_data.end_point = Some((position, id));
-}
diff --git a/libraries/bezier-rs/src/bezier_fit.rs b/libraries/bezier-rs/src/bezier_fit.rs
new file mode 100644
index 000000000..897890d1b
--- /dev/null
+++ b/libraries/bezier-rs/src/bezier_fit.rs
@@ -0,0 +1,552 @@
+/*
+* Modifications and Rust port copyright (C) 2024 by 0Hypercube.
+*
+* Original version by lib2geom: <https://gitlab.com/inkscape/lib2geom>
+*
+* The entirety of this file is specially licensed under MPL 1.1 terms:
+*
+*  Original code published in:
+*    An Algorithm for Automatically Fitting Digitized Curves
+*    by Philip J. Schneider
+*   "Graphics Gems", Academic Press, 1990
+*
+*  Authors:
+*    Philip J. Schneider
+*    Lauris Kaplinski <lauris@kaplinski.com>
+*    Peter Moulder <pmoulder@mail.csse.monash.edu.au>
+*
+*  Copyright (C) 1990 Philip J. Schneider
+*  Copyright (C) 2001 Lauris Kaplinski
+*  Copyright (C) 2001 Ximian, Inc.
+*  Copyright (C) 2003,2004 Monash University
+
+*  Original authors listed in the version control history of the following files:
+*  - https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier-utils.cpp
+*
+* This file is free software; you can redistribute it and/or modify it
+* either under the terms of the Mozilla Public License Version 1.1 (the
+* "MPL").
+*
+* The contents of this file are subject to the Mozilla Public License
+* Version 1.1 (the "License"); you may not use this file except in
+* compliance with the License. You may obtain a copy of the License at
+* https://www.mozilla.org/MPL/1.1/
+*
+* This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+* OF ANY KIND, either express or implied. See the MPL for the specific
+* language governing rights and limitations.
+*/
+
+use glam::DVec2;
+
+// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier-utils.cpp#L192
+fn bezier_fit_cubic_full(bezier: &mut [DVec2], data: &[DVec2], t_hat_1: DVec2, t_hat_2: DVec2, error: f64, max_beziers: usize) -> Option<usize> {
+	if data.len() < 2 {
+		return Some(0);
+	}
+	if data.len() == 2 {
+		// Fit 2 points
+		bezier[0] = data[0];
+		bezier[3] = data[data.len() - 1];
+		let dist = bezier[0].distance(bezier[3]) / 3.;
+		if dist.is_finite() {
+			bezier[1] = if t_hat_1 == DVec2::ZERO { (2. * bezier[0] + bezier[3]) / 3. } else { bezier[0] + dist * t_hat_1 };
+			bezier[2] = if t_hat_2 == DVec2::ZERO { (bezier[0] + 2. * bezier[3]) / 3. } else { bezier[3] + dist * t_hat_2 };
+		} else {
+			bezier[1] = bezier[0];
+			bezier[2] = bezier[3];
+		}
+		return Some(1);
+	}
+
+	// Parameterize points, and attempt to fit curve
+
+	let mut u = chord_length_parameterize(data);
+	if u[data.len() - 1] == 0. {
+		return Some(0);
+	}
+
+	generate_bezier(bezier, data, &u, t_hat_1, t_hat_2, error);
+	reparameterize(data, &mut u, bezier);
+
+	// Find max deviation of points to fitted curve.
+	let tolerance = (error + 1e-9).sqrt();
+	let (mut split_point, mut max_error_ratio) = compute_max_error_ratio(data, &u, bezier, tolerance, 0);
+
+	if max_error_ratio.abs() <= 1. {
+		return Some(1);
+	}
+
+	// If error not too large, then try some reparameterization and iteration.
+	if (0.0..=3.).contains(&max_error_ratio) {
+		let max_iterations = 4; // Max times to try iterating
+		for _ in 0..max_iterations {
+			generate_bezier(bezier, data, &u, t_hat_1, t_hat_2, error);
+			reparameterize(data, &mut u, bezier);
+			(split_point, max_error_ratio) = compute_max_error_ratio(data, &u, bezier, tolerance, split_point);
+			if (max_error_ratio.abs()) <= 1. {
+				return Some(1);
+			}
+		}
+	}
+	let is_corner = max_error_ratio < 0.;
+
+	if is_corner {
+		assert!(split_point < data.len());
+		if split_point == 0 {
+			if t_hat_1 == DVec2::ZERO {
+				// Got spike even with unconstrained initial tangent.
+				split_point += 1;
+			} else {
+				return bezier_fit_cubic_full(bezier, data, DVec2::ZERO, t_hat_2, error, max_beziers);
+			}
+		} else if split_point == data.len() - 1 {
+			if t_hat_2 == DVec2::ZERO {
+				// Got spike even with unconstrained final tangent.
+				split_point -= 1;
+			} else {
+				return bezier_fit_cubic_full(bezier, data, t_hat_1, DVec2::ZERO, error, max_beziers);
+			}
+		}
+	}
+
+	if 1 < max_beziers {
+		// Fitting failed -- split at max error point and fit recursively
+		let rec_max_beziers1 = max_beziers - 1;
+
+		let [rec_t_hat_1, rec_t_hat_2] = if is_corner {
+			if !(0 < split_point && split_point < data.len() - 1) {
+				return None;
+			}
+			[DVec2::ZERO; 2]
+		} else {
+			// Unit tangent vector at splitPoint.
+			let rec_t_hat_2 = darray_center_tangent(data, split_point, data.len());
+			[-rec_t_hat_2, rec_t_hat_2]
+		};
+		let nsegs1 = bezier_fit_cubic_full(bezier, &data[..split_point + 1], t_hat_1, rec_t_hat_2, error, rec_max_beziers1)?; //
+
+		assert!(nsegs1 != 0);
+		let rec_max_beziers2 = max_beziers - nsegs1;
+		let nsegs2 = bezier_fit_cubic_full(&mut bezier[nsegs1 * 4..], &data[split_point..], rec_t_hat_1, t_hat_2, error, rec_max_beziers2)?;
+		Some(nsegs1 + nsegs2)
+	} else {
+		None
+	}
+}
+
+// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier-utils.cpp#L538
+fn reparameterize(d: &[DVec2], u: &mut [f64], bez_curve: &[DVec2]) {
+	let len = u.len();
+	assert!(2 <= len);
+
+	let last = len - 1;
+	assert!(bez_curve[0] == d[0]);
+	assert!(bez_curve[3] == d[last]);
+	assert!(u[0] == 0.);
+	assert!(u[last] == 1.);
+	// Otherwise, consider including 0 and last in the below loop.
+
+	for i in 1..last {
+		u[i] = newton_raphson_root_find(bez_curve, d[i], u[i]);
+	}
+}
+
+// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier-utils.cpp#L567
+fn newton_raphson_root_find(q: &[DVec2], p: DVec2, u: f64) -> f64 {
+	assert!(0. <= u);
+	assert!(u <= 1.);
+
+	// Generate control vertices for Q'.
+	let mut q1 = [DVec2::ZERO; 3];
+	for i in 0..3 {
+		q1[i] = 3. * (q[i + 1] - q[i]);
+	}
+
+	// Generate control vertices for Q''.
+	let mut q2 = [DVec2::ZERO; 2];
+	for i in 0..2 {
+		q2[i] = 2. * (q1[i + 1] - q1[i]);
+	}
+
+	// Compute Q(u), Q'(u) and Q''(u).
+
+	let q_u = crate::Bezier::from_cubic_dvec2(q[0], q[1], q[2], q[3]).evaluate(crate::TValue::Parametric(u));
+	let q1_u = crate::Bezier::from_quadratic_dvec2(q1[0], q1[1], q1[2]).evaluate(crate::TValue::Parametric(u));
+	let q2_u = crate::Bezier::from_linear_dvec2(q2[0], q2[1]).evaluate(crate::TValue::Parametric(u));
+
+	// Compute f(u)/f'(u), where f is the derivative wrt u of distsq(u) = 0.5 * the square of the
+	// distance from P to Q(u).  Here we're using Newton-Raphson to find a stationary point in the
+	// distsq(u), hopefully corresponding to a local minimum in distsq (and hence a local minimum
+	// distance from P to Q(u)).
+	let diff = q_u - p;
+	let numerator = diff.dot(q1_u);
+	let denominator = q1_u.dot(q1_u) + diff.dot(q2_u);
+
+	let mut improved_u = if denominator > 0. {
+		// One iteration of Newton-Raphson:
+		// improved_u = u - f(u)/f'(u)
+		u - (numerator / denominator)
+	} else {
+		// Using Newton-Raphson would move in the wrong direction (towards a local maximum rather
+		// than local minimum), so we move an arbitrary amount in the right direction.
+		if numerator > 0. {
+			u * 0.98 - 0.001
+		} else if numerator < 0. {
+			// Deliberately asymmetrical, to reduce the chance of cycling.
+			0.031 + u * 0.98
+		} else {
+			u
+		}
+	};
+
+	if !improved_u.is_finite() {
+		improved_u = u;
+	} else {
+		improved_u = improved_u.clamp(0., 1.);
+	}
+
+	// Ensure that improved_u isn't actually worse.
+	{
+		let diff_lensq = diff.length_squared();
+		let mut proportion = 0.125;
+		loop {
+			let bezier_pt = crate::Bezier::from_cubic_dvec2(q[0], q[1], q[2], q[3]).evaluate(crate::TValue::Parametric(improved_u));
+			if (bezier_pt - p).length_squared() > diff_lensq {
+				if proportion > 1. {
+					// g_warning("found proportion %g", proportion);
+					improved_u = u;
+					break;
+				}
+				improved_u = (1. - proportion) * improved_u + proportion * u;
+			} else {
+				break;
+			}
+			proportion += 0.125;
+		}
+	}
+
+	improved_u
+}
+
+// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier-utils.cpp#L807
+fn darray_center_tangent(d: &[DVec2], center: usize, len: usize) -> DVec2 {
+	assert!(center != 0);
+	assert!(center < len - 1);
+
+	if d[center + 1] == d[center - 1] {
+		// Rotate 90 degrees in an arbitrary direction.
+		let diff = d[center] - d[center - 1];
+		diff.perp()
+	} else {
+		d[center - 1] - d[center + 1]
+	}
+	.normalize()
+}
+
+// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier-utils.cpp#L387
+fn estimate_lengths(bezier: &mut [DVec2], data: &[DVec2], u_prime: &[f64], t_hat_1: DVec2, t_hat_2: DVec2) {
+	let len = data.len();
+	let mut c = [[0.; 2]; 2]; // Matrix C
+	let mut x = [0.; 2]; // Matrix X
+
+	// First and last control points of the Bezier curve are positioned exactly at the first and last data points.
+	bezier[0] = data[0];
+	bezier[3] = data[len - 1];
+
+	for i in 0..len {
+		// Bezier control point coefficients.
+		let b0 = (1. - u_prime[i]) * (1. - u_prime[i]) * (1. - u_prime[i]);
+		let b1 = 3. * u_prime[i] * (1. - u_prime[i]) * (1. - u_prime[i]);
+		let b2 = 3. * u_prime[i] * u_prime[i] * (1. - u_prime[i]);
+		let b3 = u_prime[i] * u_prime[i] * u_prime[i];
+
+		// rhs for eqn
+		let a1 = b1 * t_hat_1;
+		let a2 = b2 * t_hat_2;
+
+		c[0][0] += a1.dot(a1);
+		c[0][1] += a1.dot(a2);
+		c[1][0] = c[0][1];
+		c[1][1] += a2.dot(a2);
+
+		// Additional offset to the data point from the predicted point if we were to set bezier[1] to bezier[0] and bezier[2] to bezier[3].
+		let shortfall = data[i] - ((b0 + b1) * bezier[0]) - ((b2 + b3) * bezier[3]);
+		x[0] += a1.dot(shortfall);
+		x[1] += a2.dot(shortfall);
+	}
+
+	// We've constructed a pair of equations in the form of a matrix product C * alpha = X. Now solve for alpha.
+
+	// Compute the determinants of C and X.
+	let det_c0_c1 = c[0][0] * c[1][1] - c[1][0] * c[0][1];
+	let [mut alpha_l, mut alpha_r] = if det_c0_c1 != 0. {
+		// Apparently Kramer's rule.
+		let det_c0_x = c[0][0] * x[1] - c[0][1] * x[0];
+		let det_x_c1 = x[0] * c[1][1] - x[1] * c[0][1];
+		[det_x_c1 / det_c0_c1, det_c0_x / det_c0_c1]
+	} else {
+		// The matrix is under-determined. Try requiring alpha_l == alpha_r.
+		//
+		// One way of implementing the constraint alpha_l == alpha_r is to treat them as the same
+		// variable in the equations.  We can do this by adding the columns of C to form a single
+		// column, to be multiplied by alpha to give the column vector X.
+		//
+		// We try each row in turn.
+		let c0 = c[0][0] + c[0][1];
+		if c0 != 0. {
+			[x[0] / c0; 2]
+		} else {
+			let c1 = c[1][0] + c[1][1];
+			if c1 != 0. {
+				[x[1] / c1; 2]
+			} else {
+				// Let the below code handle this.
+				[0.; 2]
+			}
+		}
+	};
+
+	// If alpha negative, use the Wu/Barsky heuristic (see text).
+	// (If alpha is 0, you get coincident control points that lead to divide by zero in any subsequent NewtonRaphsonRootFind() call.)
+	// TODO: Check whether this special-casing is necessary now that NewtonRaphsonRootFind handles non-positive denominator.
+	if alpha_l < 1.0e-6 || alpha_r < 1.0e-6 {
+		alpha_l = data[0].distance(data[len - 1]) / 3.;
+		alpha_r = alpha_l;
+	}
+
+	// Control points 1 and 2 are positioned an alpha distance out on the tangent vectors, left and right, respectively.
+	bezier[1] = alpha_l * t_hat_1 + bezier[0];
+	bezier[2] = alpha_r * t_hat_2 + bezier[3];
+}
+
+// ComputeLeftTangent, ComputeRightTangent, ComputeCenterTangent:
+// Approximate unit tangents at endpoints and "center" of digitized curve
+
+// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier-utils.cpp#L706
+fn darray_left_tangent_simple(d: &[DVec2]) -> DVec2 {
+	let len = d.len();
+	assert!(len >= 2);
+	assert!(d[0] != d[1]);
+	(d[1] - d[0]).normalize()
+}
+
+// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier-utils.cpp#L724
+fn darray_right_tangent_simple(d: &[DVec2]) -> DVec2 {
+	let len = d.len();
+	assert!(2 <= len);
+	let last = len - 1;
+	let prev = last - 1;
+	assert!(d[last] != d[prev]);
+	(d[prev] - d[last]).normalize()
+}
+
+// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier-utils.cpp#L745
+fn darray_left_tangent(d: &[DVec2], tolerance_sq: f64) -> DVec2 {
+	let len = d.len();
+	assert!(2 <= len);
+	assert!(0. <= tolerance_sq);
+	let mut i = 1;
+	loop {
+		let pi = d[i];
+		let t = pi - d[0];
+		let distsq = t.length_squared();
+		if tolerance_sq < distsq {
+			return (t).normalize();
+		}
+		i += 1;
+		if i == len {
+			return if distsq == 0. { darray_left_tangent_simple(d) } else { (t).normalize() };
+		}
+	}
+}
+
+// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier-utils.cpp#L776
+fn darray_right_tangent(d: &[DVec2], tolerance_sq: f64) -> DVec2 {
+	let len = d.len();
+	assert!(2 <= len);
+	assert!(0. <= tolerance_sq);
+	let last = len - 1;
+	let mut i = last - 1;
+	loop {
+		let pi = d[i];
+		let t = pi - d[last];
+		let dist_sq = t.length_squared();
+		if tolerance_sq < dist_sq {
+			return t.normalize_or_zero();
+		}
+		if i == 0 {
+			return if dist_sq == 0. { darray_right_tangent_simple(d) } else { (t).normalize_or_zero() };
+		}
+		i -= 1;
+	}
+}
+// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier-utils.cpp#L360
+fn generate_bezier(bezier: &mut [DVec2], data: &[DVec2], u: &[f64], t_hat_1: DVec2, t_hat_2: DVec2, tolerance_sq: f64) {
+	let est1 = t_hat_1 == DVec2::ZERO;
+	let est2 = t_hat_2 == DVec2::ZERO;
+	let mut est_t_hat_1 = if est1 { darray_left_tangent(data, tolerance_sq) } else { t_hat_1 };
+	let est_t_hat_2 = if est2 { darray_right_tangent(data, tolerance_sq) } else { t_hat_2 };
+	estimate_lengths(bezier, data, u, est_t_hat_1, est_t_hat_2);
+	// We find that darray_right_tangent tends to produce better results for our current freehand tool than full estimation.
+	if est1 {
+		estimate_bi(bezier, 1, data, u);
+		if bezier[1] != bezier[0] {
+			est_t_hat_1 = (bezier[1] - bezier[0]).normalize_or_zero();
+		}
+		estimate_lengths(bezier, data, u, est_t_hat_1, est_t_hat_2);
+	}
+}
+
+// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier-utils.cpp#L492
+fn estimate_bi(bezier: &mut [DVec2], ei: usize, data: &[DVec2], u: &[f64]) {
+	if !(1..=2).contains(&ei) {
+		return;
+	}
+	let oi = 3 - ei;
+	let mut num = [0., 0.];
+	let mut den = 0.;
+	for i in 0..data.len() {
+		let ui = u[i];
+		let b = [((1. - ui) * (1. - ui) * (1. - ui)), (3. * ui * (1. - ui) * (1. - ui)), (3. * ui * ui * (1. - ui)), (ui * ui * ui)];
+
+		for (d, num_d) in num.iter_mut().enumerate() {
+			*num_d += b[ei] * (b[0] * bezier[0][d] + b[oi] * bezier[oi][d] + b[3] * bezier[3][d] + -data[i][d]);
+		}
+		den -= b[ei] * b[ei];
+	}
+
+	if den != 0. {
+		for (d, num_d) in num.iter().enumerate() {
+			bezier[ei][d] = num_d / den;
+		}
+	} else {
+		bezier[ei] = (oi as f64 * bezier[0] + ei as f64 * bezier[3]) / 3.;
+	}
+}
+
+// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier-utils.cpp#L898
+fn compute_max_error_ratio(d: &[DVec2], u: &[f64], bezier_curve: &[DVec2], tolerance: f64, mut split_point: usize) -> (usize, f64) {
+	let len = u.len();
+	assert!(2 <= len);
+	let last = len - 1;
+	assert!(bezier_curve[0] == d[0]);
+	assert!(bezier_curve[3] == d[last]);
+	assert!(u[0] == 0.);
+	assert!(u[last] == 1.);
+	// I.e. assert that the error for the first & last points is zero.
+	// Otherwise we should include those points in the below loop.
+	// The assertion is also necessary to ensure 0 < splitPoint < last.
+
+	let mut max_distance_sq = 0.; // Maximum error
+	let mut max_hook_ratio = 0.;
+	let mut snap_end = 0;
+	let mut previous_point = bezier_curve[0];
+	for i in 1..=last {
+		let current_point = crate::Bezier::from_cubic_dvec2(bezier_curve[0], bezier_curve[1], bezier_curve[2], bezier_curve[3]).evaluate(crate::TValue::Parametric(u[i]));
+		let distsq = (current_point - d[i]).length_squared();
+		if distsq > max_distance_sq {
+			max_distance_sq = distsq;
+			split_point = i;
+		}
+		let hook_ratio = compute_hook(previous_point, current_point, 0.5 * (u[i - 1] + u[i]), bezier_curve, tolerance);
+		if max_hook_ratio < hook_ratio {
+			max_hook_ratio = hook_ratio;
+			snap_end = i;
+		}
+		previous_point = current_point;
+	}
+
+	let dist_ratio = (max_distance_sq).sqrt() / tolerance;
+	let error_ratio = if max_hook_ratio <= dist_ratio {
+		dist_ratio
+	} else {
+		assert!(0 < snap_end);
+		split_point = snap_end - 1;
+		-max_hook_ratio
+	};
+	assert!(error_ratio == 0. || ((split_point < last) && (split_point != 0 || error_ratio < 0.)));
+	(split_point, error_ratio)
+}
+
+// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier-utils.cpp#L969
+fn compute_hook(a: DVec2, b: DVec2, u: f64, bezier_curve: &[DVec2], tolerance: f64) -> f64 {
+	let point = crate::Bezier::from_cubic_dvec2(bezier_curve[0], bezier_curve[1], bezier_curve[2], bezier_curve[3]).evaluate(crate::TValue::Parametric(u));
+
+	let distance = ((a + b) / 2.).distance(point);
+	if distance < tolerance {
+		return 0.;
+	}
+	let allowed = a.distance(b) + tolerance;
+	distance / allowed
+}
+
+/// A value from 0..1 for each point in the path containing the total distance from the start along the path
+// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier-utils.cpp#L833
+fn chord_length_parameterize(d: &[DVec2]) -> Vec<f64> {
+	let len = d.len();
+	if len < 2 {
+		return Vec::new();
+	}
+
+	// First let u[i] equal the distance travelled along the path from d[0] to d[i]
+	let mut u = vec![0.; len];
+	for i in 1..len {
+		let dist = d[i].distance(d[i - 1]);
+		u[i] = u[i - 1] + dist;
+	}
+
+	// Scale to 0 ..= 1
+	let tot_len = u[len - 1];
+	if tot_len <= 0. {
+		return Vec::new();
+	}
+	if tot_len.is_finite() {
+		for u in u.iter_mut().skip(1) {
+			*u /= tot_len;
+		}
+	} else {
+		// fallback to even space
+		for (i, u) in u.iter_mut().enumerate().skip(1) {
+			*u = i as f64 / (len - 1) as f64;
+		}
+	}
+
+	// Ensure last is exactly 1.
+	u[len - 1] = 1.;
+	u
+}
+
+impl<PointId: crate::Identifier> crate::Subpath<PointId> {
+	pub fn fit_cubic(points: &[DVec2], max_segs: usize, tangent: DVec2, tolerance_sq: f64) -> Option<Self> {
+		let mut b = vec![DVec2::ZERO; 4 * max_segs];
+		let len = bezier_fit_cubic_full(&mut b, points, tangent, DVec2::ZERO, tolerance_sq, max_segs)?;
+		if len < 1 {
+			return None;
+		}
+		let beziers = (0..len)
+			.map(|i| crate::Bezier::from_cubic_dvec2(b[i * 4], b[i * 4 + 1], b[i * 4 + 2], b[i * 4 + 3]))
+			.collect::<Vec<_>>();
+		Some(Self::from_beziers(&beziers, false))
+	}
+}
+
+#[test]
+fn generate_bezier_test() {
+	let src_b = vec![DVec2::new(5., -3.), DVec2::new(8., 0.), DVec2::new(4., 2.), DVec2::new(3., 3.)];
+	let t = [
+		0., 0.001, 0.03, 0.05, 0.09, 0.13, 0.18, 0.25, 0.29, 0.33, 0.39, 0.44, 0.51, 0.57, 0.62, 0.69, 0.75, 0.81, 0.91, 0.93, 0.97, 0.98, 0.999, 1.,
+	];
+
+	let data = t
+		.iter()
+		.map(|&t| crate::Bezier::from_cubic_dvec2(src_b[0], src_b[1], src_b[2], src_b[3]).evaluate(crate::TValue::Parametric(t)))
+		.collect::<Vec<_>>();
+	let t_hat_1 = (src_b[1] - src_b[0]).normalize();
+	let t_hat_2 = (src_b[2] - src_b[3]).normalize();
+
+	let mut est_b = vec![DVec2::ZERO; 4];
+	generate_bezier(&mut est_b, &data, &t, t_hat_1, t_hat_2, 1.);
+
+	assert_eq!(src_b, est_b);
+}
diff --git a/libraries/bezier-rs/src/lib.rs b/libraries/bezier-rs/src/lib.rs
index 06faf15d8..0be07df48 100644
--- a/libraries/bezier-rs/src/lib.rs
+++ b/libraries/bezier-rs/src/lib.rs
@@ -4,6 +4,7 @@
 pub(crate) mod compare;
 
 mod bezier;
+mod bezier_fit;
 mod consts;
 mod poisson_disk;
 mod polynomial;
@@ -12,6 +13,7 @@ mod symmetrical_basis;
 mod utils;
 
 pub use bezier::*;
+pub use bezier_fit::*;
 pub use subpath::*;
 pub use symmetrical_basis::*;
 pub use utils::{Cap, Join, SubpathTValue, TValue, TValueType};
diff --git a/website/other/bezier-rs-demos/src/features/subpath-features.ts b/website/other/bezier-rs-demos/src/features/subpath-features.ts
index b67d62f1c..8ddc33d08 100644
--- a/website/other/bezier-rs-demos/src/features/subpath-features.ts
+++ b/website/other/bezier-rs-demos/src/features/subpath-features.ts
@@ -7,6 +7,27 @@ const subpathFeatures = {
 		name: "Constructor",
 		callback: (subpath: WasmSubpathInstance): string => subpath.to_svg(),
 	},
+	fit: {
+		name: "Fit Points",
+		callback: (subpath: WasmSubpathInstance, options: Record<string, number>, _: undefined): string => subpath.fit(options.max_segments, options.error),
+
+		inputOptions: [
+			{
+				variable: "error",
+				min: 0,
+				max: 25,
+				step: 1,
+				default: 10,
+			},
+			{
+				variable: "max_segments",
+				min: 1,
+				max: 3,
+				step: 1,
+				default: 2,
+			},
+		],
+	},
 	insert: {
 		name: "Insert",
 		callback: (subpath: WasmSubpathInstance, options: Record<string, number>, _: undefined): string => subpath.insert(options.t, SUBPATH_T_VALUE_VARIANTS[options.TVariant]),
diff --git a/website/other/bezier-rs-demos/wasm/src/subpath.rs b/website/other/bezier-rs-demos/wasm/src/subpath.rs
index ecf62d7cf..afd287af2 100644
--- a/website/other/bezier-rs-demos/wasm/src/subpath.rs
+++ b/website/other/bezier-rs-demos/wasm/src/subpath.rs
@@ -76,6 +76,29 @@ impl WasmSubpath {
 		subpath_svg
 	}
 
+	pub fn fit(&self, max_segments: f64, tolerance: f64) -> String {
+		let mut result = String::new();
+
+		let subpath = self.0.clone();
+		let points = subpath.manipulator_groups().iter().map(|group| group.anchor).collect::<Vec<_>>();
+
+		let tolerance_sq = 0.02 * tolerance * tolerance * (0.2 * tolerance - 2.).exp();
+		if let Some(subpath) = Subpath::<EmptyId>::fit_cubic(&points, (max_segments as usize).max(1), DVec2::ZERO, tolerance_sq) {
+			subpath.to_svg(
+				&mut result,
+				CURVE_ATTRIBUTES.to_string().replace(BLACK, RED),
+				ANCHOR_ATTRIBUTES.to_string().replace(BLACK, RED),
+				HANDLE_ATTRIBUTES.to_string().replace(GRAY, RED),
+				HANDLE_LINE_ATTRIBUTES.to_string().replace(GRAY, RED),
+			);
+		}
+		for point in points {
+			result += &draw_circle(point, 4., BLACK, 1.5, WHITE);
+		}
+
+		wrap_svg_tag(result)
+	}
+
 	fn to_filled_svg(&self) -> String {
 		let mut subpath_svg = String::new();
 		self.0.to_svg(