Implement function that projects to a Bezier curve (#688)

* UI section for the projection function

* added bezier project impl

* Fix project function and add test for it

* Search method

* Re-use comptued distances

* Update comments

* rebase project changes

* clean up tests and library code

* use built-in functions and destructure syntax

* Remove redundant project implementation

* Fix typo, add lut size as parameter and add constant

* address comments

Co-authored-by: Thomas Cheng <contact.chengthomas@gmail.com>

bezier-rs/docs/interactive-docs/src/App.vue

@@ -21,7 +21,7 @@
 import { defineComponent, markRaw } from "vue";
 
 import { drawText, drawPoint, drawBezier, drawLine, getContextFromCanvas, drawBezierHelper, COLORS } from "@/utils/drawing";
-import { WasmBezierInstance } from "@/utils/types";
+import { Point, WasmBezierInstance } from "@/utils/types";
 
 import ExamplePane from "@/components/ExamplePane.vue";
 import SliderExample from "@/components/SliderExample.vue";
@@ -226,6 +226,16 @@ export default defineComponent({
 						],
 					},
 				},
+				{
+					name: "Project",
+					callback: (canvas: HTMLCanvasElement, bezier: WasmBezierInstance, options: Record<string, number>, mouseLocation: Point | null): void => {
+						if (mouseLocation != null) {
+							const context = getContextFromCanvas(canvas);
+							const closestPoint = JSON.parse(bezier.project(mouseLocation.x, mouseLocation.y));
+							drawLine(context, mouseLocation, closestPoint, COLORS.NON_INTERACTIVE.STROKE_1);
+						}
+					},
+				},
 			],
 		};
 	},

bezier-rs/docs/interactive-docs/src/components/BezierDrawing.ts

@@ -1,6 +1,6 @@
 import { WasmBezier } from "@/../wasm/pkg";
 import { COLORS, drawBezier, drawPoint, getContextFromCanvas, getPointSizeByIndex } from "@/utils/drawing";
-import { BezierCallback, BezierPoint, BezierStyleConfig, WasmBezierMutatorKey, WasmBezierInstance } from "@/utils/types";
+import { BezierCallback, BezierPoint, BezierStyleConfig, Point, WasmBezierMutatorKey, WasmBezierInstance } from "@/utils/types";
 
 // Offset to increase selectable range, used to make points easier to grab
 const FUDGE_FACTOR = 3;
@@ -81,10 +81,9 @@ class BezierDrawing {
 				selectedPoint.x = mx;
 				selectedPoint.y = my;
 				this.bezier[selectedPoint.mutator](selectedPoint.x, selectedPoint.y);
-				this.clearFigure();
 			}
 		}
-		this.updateBezier();
+		this.updateBezier({ x: mx, y: my });
 	}
 
 	mouseDownHandler(evt: MouseEvent): void {
@@ -102,13 +101,13 @@ class BezierDrawing {
 
 	deselectPointHandler(): void {
 		if (this.dragIndex !== undefined) {
-			this.clearFigure();
 			this.dragIndex = null;
 			this.updateBezier();
 		}
 	}
 
-	updateBezier(options: Record<string, number> = {}): void {
+	updateBezier(mouseLocation?: Point, options: Record<string, number> = {}): void {
+		this.clearFigure();
 		if (Object.values(options).length !== 0) {
 			this.options = options;
 		}
@@ -146,7 +145,7 @@ class BezierDrawing {
 			// Draw the point that the curve was drawn through
 			drawPoint(this.ctx, this.points[1], getPointSizeByIndex(1, this.points.length), this.dragIndex === 1 ? COLORS.INTERACTIVE.SELECTED : COLORS.INTERACTIVE.STROKE_1);
 		}
-		this.callback(this.canvas, this.bezier, this.options);
+		this.callback(this.canvas, this.bezier, this.options, mouseLocation);
 	}
 
 	getCanvas(): HTMLCanvasElement {

bezier-rs/docs/interactive-docs/src/components/Example.vue

@@ -46,7 +46,7 @@ export default defineComponent({
 		options: {
 			deep: true,
 			handler() {
-				this.bezierDrawing.updateBezier(this.options);
+				this.bezierDrawing.updateBezier(undefined, this.options);
 			},
 		},
 	},

bezier-rs/docs/interactive-docs/src/utils/types.ts

@@ -4,7 +4,7 @@ export type WasmBezierInstance = InstanceType<WasmRawInstance["WasmBezier"]>;
 export type WasmBezierKey = keyof WasmBezierInstance;
 export type WasmBezierMutatorKey = "set_start" | "set_handle_start" | "set_handle_end" | "set_end";
 
-export type BezierCallback = (canvas: HTMLCanvasElement, bezier: WasmBezierInstance, options: Record<string, number>) => void;
+export type BezierCallback = (canvas: HTMLCanvasElement, bezier: WasmBezierInstance, options: Record<string, number>, mouseLocation?: Point) => void;
 
 export type SliderOption = {
 	min: number;

bezier-rs/docs/interactive-docs/wasm/src/lib.rs

@@ -98,4 +98,8 @@ impl WasmBezier {
 	pub fn trim(&self, t1: f64, t2: f64) -> WasmBezier {
 		WasmBezier(self.0.trim(t1, t2))
 	}
+
+	pub fn project(&self, x: f64, y: f64) -> JsValue {
+		vec_to_point(&self.0.project(DVec2::new(x, y), 20, 1e-4, 3, 10))
+	}
 }

bezier-rs/lib/src/lib.rs

@@ -331,6 +331,93 @@ impl Bezier {
 		};
 		bezier_starting_at_t1.split(adjusted_t2)[t2_split_side]
 	}
+
+	/// Returns the closest point on the curve to the provided point.
+	/// Uses a searching algorithm akin to binary search that can be customized using the following parameters:
+	/// - `lut_size` - Size of the lookup table for the initial passthrough
+	/// - `convergence_epsilon` - Difference used between floating point numbers to be considered as equal
+	/// - `convergence_limit` - Controls the number of iterations needed to consider that minimum distance to have converged
+	/// - `iteration_limit` - Controls the maximum total number of iterations to be used
+	pub fn project(&self, point: DVec2, lut_size: i32, convergence_epsilon: f64, convergence_limit: i32, iteration_limit: i32) -> DVec2 {
+		// First find the closest point from the results of a lookup table
+		let lut = self.compute_lookup_table(Some(lut_size));
+		let (minimum_position, minimum_distance) = utils::get_closest_point_in_lut(&lut, point);
+
+		// Get the t values to the left and right of the closest result in the lookup table
+		let mut left_t = (0.max(minimum_position - 1) as f64) / lut_size as f64;
+		let mut right_t = (lut_size.min(minimum_position + 1)) as f64 / lut_size as f64;
+
+		// Perform a finer search by finding closest t from 5 points between [left_t, right_t] inclusive
+		// Choose new left_t and right_t for a smaller range around the closest t and repeat the process
+		let mut final_t = left_t;
+		let mut distance;
+
+		// Increment minimum_distance to ensure that the distance < minimum_distance comparison will be true for at least one iteration
+		let mut new_minimum_distance = minimum_distance + 1.;
+		// Maintain the previous distance to identify convergence
+		let mut previous_distance;
+		// Counter to limit the number of iterations
+		let mut iteration_count = 0;
+		// Counter to identify how many iterations have had a similar result. Used for convergence test
+		let mut convergence_count = 0;
+		// Store calculated distances to minimize unnecessary recomputations
+		const NUM_DISTANCES: usize = 5;
+		let mut distances: [f64; NUM_DISTANCES] = [
+			point.distance(lut[0.max(minimum_position - 1) as usize]),
+			0.,
+			0.,
+			0.,
+			point.distance(lut[lut_size.min(minimum_position + 1) as usize]),
+		];
+
+		while left_t <= right_t && convergence_count < convergence_limit && iteration_count < iteration_limit {
+			previous_distance = new_minimum_distance;
+			let step = (right_t - left_t) / ((NUM_DISTANCES - 1) as f64);
+			let mut iterator_t = left_t;
+			let mut target_index = 0;
+			// Iterate through first 4 points and will handle the right most point later
+			for (step_index, table_distance) in distances.iter_mut().enumerate().take(4) {
+				// Use previously computed distance for the left most point, and compute new values for the others
+				if step_index == 0 {
+					distance = *table_distance;
+				} else {
+					distance = point.distance(self.compute(iterator_t));
+					*table_distance = distance;
+				}
+				if distance < new_minimum_distance {
+					new_minimum_distance = distance;
+					target_index = step_index;
+					final_t = iterator_t
+				}
+				iterator_t += step;
+			}
+			// Check right most edge separately since step may not perfectly add up to it (floating point errors)
+			if distances[NUM_DISTANCES - 1] < new_minimum_distance {
+				new_minimum_distance = distances[NUM_DISTANCES - 1];
+				final_t = right_t;
+			}
+
+			// Update left_t and right_t to be the t values (final_t +/- step), while handling the edges (i.e. if final_t is 0, left_t will be 0 instead of -step)
+			// Ensure that the t values never exceed the [0, 1] range
+			left_t = (final_t - step).max(0.);
+			right_t = (final_t + step).min(1.);
+
+			// Re-use the corresponding computed distances (target_index is the index corresponding to final_t)
+			// Since target_index is a u_size, can't subtract one if it is zero
+			distances[0] = distances[if target_index == 0 { 0 } else { target_index - 1 }];
+			distances[NUM_DISTANCES - 1] = distances[(target_index + 1).min(NUM_DISTANCES - 1)];
+
+			iteration_count += 1;
+			// update count for consecutive iterations of similar minimum distances
+			if previous_distance - new_minimum_distance < convergence_epsilon {
+				convergence_count += 1;
+			} else {
+				convergence_count = 0;
+			}
+		}
+
+		self.compute(final_t)
+	}
 }
 
 #[cfg(test)]
@@ -339,7 +426,7 @@ mod tests {
 	use glam::DVec2;
 
 	fn compare_points(p1: DVec2, p2: DVec2) -> bool {
-		DVec2::new(0.001, 0.001).cmpge(p1 - p2).all()
+		p1.abs_diff_eq(p2, 0.001)
 	}
 
 	#[test]
@@ -373,4 +460,14 @@ mod tests {
 		let bezier3 = Bezier::cubic_through_points(p1, p2, p3, 0., 91.7);
 		assert!(compare_points(bezier3.compute(0.), p2));
 	}
+
+	#[test]
+	fn project() {
+		let bezier1 = Bezier::from_cubic_coordinates(4., 4., 23., 45., 10., 30., 56., 90.);
+		assert!(bezier1.project(DVec2::new(100., 100.), 20, 0.0001, 3, 10) == DVec2::new(56., 90.));
+		assert!(bezier1.project(DVec2::new(0., 0.), 20, 0.0001, 3, 10) == DVec2::new(4., 4.));
+
+		let bezier2 = Bezier::from_quadratic_coordinates(0., 0., 0., 100., 100., 100.);
+		assert!(bezier2.project(DVec2::new(100., 0.), 20, 0.0001, 3, 10) == DVec2::new(0., 0.));
+	}
 }

bezier-rs/lib/src/utils.rs

@@ -29,3 +29,11 @@ pub fn compute_abc_for_cubic_through_points(start_point: DVec2, point_on_curve:
 
 	compute_abc_through_points(start_point, point_on_curve, end_point, t_cubed, cubed_one_minus_t)
 }
+
+pub fn get_closest_point_in_lut(lut: &[DVec2], point: DVec2) -> (i32, f64) {
+	lut.iter()
+		.enumerate()
+		.map(|(i, p)| (i as i32, point.distance(*p)))
+		.min_by(|x, y| (&(x.1)).partial_cmp(&(y.1)).unwrap())
+		.unwrap()
+}