Add curve closest point finding and distance calc

feat/remove-ga
Mark Tolmacs 1 month ago
parent 33c8c3006f
commit 22bf71e0e5
No known key found for this signature in database

@ -1,5 +1,10 @@
import "../utils/test-utils";
import { curve, curveIntersectLine } from "./curve";
import {
curve,
curveClosestPoint,
curveIntersectLine,
curvePointDistance,
} from "./curve";
import { line } from "./line";
import { pointFrom } from "./point";
@ -45,4 +50,34 @@ describe("Math curve", () => {
]);
});
});
describe("point closest to other", () => {
it("point can be found", () => {
const c = curve(
pointFrom(-50, -50),
pointFrom(10, -50),
pointFrom(10, 50),
pointFrom(50, 50),
);
const p = pointFrom(0, 0);
expect([curveClosestPoint(c, p)]).toCloselyEqualPoints([
[5.965462100367372, -3.04104878946646],
]);
});
});
describe("point shortest distance", () => {
it("can be determined", () => {
const c = curve(
pointFrom(-50, -50),
pointFrom(10, -50),
pointFrom(10, 50),
pointFrom(50, 50),
);
const p = pointFrom(0, 0);
expect(curvePointDistance(c, p)).toBeCloseTo(6.695873043213627);
});
});
});

@ -1,4 +1,4 @@
import { pointFrom } from "./point";
import { pointDistance, pointFrom } from "./point";
import type { Curve, GlobalPoint, Line, LocalPoint } from "./types";
/**
@ -18,7 +18,11 @@ export function curve<Point extends GlobalPoint | LocalPoint>(
return [a, b, c, d] as Curve<Point>;
}
/*computes intersection between a cubic spline and a line segment*/
/**
* Computes intersection between a cubic spline and a line segment
*
* @href https://www.particleincell.com/2013/cubic-line-intersection/
*/
export function curveIntersectLine<Point extends GlobalPoint | LocalPoint>(
p: Curve<Point>,
l: Line<Point>,
@ -27,8 +31,18 @@ export function curveIntersectLine<Point extends GlobalPoint | LocalPoint>(
const B = l[0][0] - l[1][0]; //B=x1-x2
const C = l[0][0] * (l[0][1] - l[1][1]) + l[0][1] * (l[1][0] - l[0][0]); //C=x1*(y1-y2)+y1*(x2-x1)
const bx = bezierCoefficients(p[0][0], p[1][0], p[2][0], p[3][0]);
const by = bezierCoefficients(p[0][1], p[1][1], p[2][1], p[3][1]);
const bx = [
-p[0][0] + 3 * p[1][0] + -3 * p[2][0] + p[3][0],
3 * p[0][0] - 6 * p[1][0] + 3 * p[2][0],
-3 * p[0][0] + 3 * p[1][0],
p[0][0],
];
const by = [
-p[0][1] + 3 * p[1][1] + -3 * p[2][1] + p[3][1],
3 * p[0][1] - 6 * p[1][1] + 3 * p[2][1],
-3 * p[0][1] + 3 * p[1][1],
p[0][1],
];
const P: [number, number, number, number] = [
A * bx[0] + B * by[0] /*t^3*/,
@ -69,7 +83,9 @@ export function curveIntersectLine<Point extends GlobalPoint | LocalPoint>(
.filter((x) => x !== null);
}
/*based on http://mysite.verizon.net/res148h4j/javascript/script_exact_cubic.html#the%20source%20code*/
/*
* Based on http://mysite.verizon.net/res148h4j/javascript/script_exact_cubic.html#the%20source%20code
*/
function cubicRoots(P: [number, number, number, number]) {
const a = P[0];
const b = P[1];
@ -128,215 +144,96 @@ function cubicRoots(P: [number, number, number, number]) {
return t;
}
function bezierCoefficients(P0: number, P1: number, P2: number, P3: number) {
const Z = [];
Z[0] = -P0 + 3 * P1 + -3 * P2 + P3;
Z[1] = 3 * P0 - 6 * P1 + 3 * P2;
Z[2] = -3 * P0 + 3 * P1;
Z[3] = P0;
return Z;
}
// export const curveRotate = <Point extends LocalPoint | GlobalPoint>(
// curve: Curve<Point>,
// angle: Radians,
// origin: Point,
// ) => {
// return curve.map((p) => pointRotateRads(p, origin, angle));
// };
// /**
// *
// * @param pointsIn
// * @param curveTightness
// * @returns
// */
// export function curveToBezier<Point extends LocalPoint | GlobalPoint>(
// pointsIn: readonly Point[],
// curveTightness = 0,
// ): Point[] {
// const len = pointsIn.length;
// if (len < 3) {
// throw new Error("A curve must have at least three points.");
// }
// const out: Point[] = [];
// if (len === 3) {
// out.push(
// pointFrom(pointsIn[0][0], pointsIn[0][1]), // Points need to be cloned
// pointFrom(pointsIn[1][0], pointsIn[1][1]), // Points need to be cloned
// pointFrom(pointsIn[2][0], pointsIn[2][1]), // Points need to be cloned
// pointFrom(pointsIn[2][0], pointsIn[2][1]), // Points need to be cloned
// );
// } else {
// const points: Point[] = [];
// points.push(pointsIn[0], pointsIn[0]);
// for (let i = 1; i < pointsIn.length; i++) {
// points.push(pointsIn[i]);
// if (i === pointsIn.length - 1) {
// points.push(pointsIn[i]);
// }
// }
// const b: Point[] = [];
// const s = 1 - curveTightness;
// out.push(pointFrom(points[0][0], points[0][1]));
// for (let i = 1; i + 2 < points.length; i++) {
// const cachedVertArray = points[i];
// b[0] = pointFrom(cachedVertArray[0], cachedVertArray[1]);
// b[1] = pointFrom(
// cachedVertArray[0] + (s * points[i + 1][0] - s * points[i - 1][0]) / 6,
// cachedVertArray[1] + (s * points[i + 1][1] - s * points[i - 1][1]) / 6,
// );
// b[2] = pointFrom(
// points[i + 1][0] + (s * points[i][0] - s * points[i + 2][0]) / 6,
// points[i + 1][1] + (s * points[i][1] - s * points[i + 2][1]) / 6,
// );
// b[3] = pointFrom(points[i + 1][0], points[i + 1][1]);
// out.push(b[1], b[2], b[3]);
// }
// }
// return out;
// }
// /**
// *
// * @param t
// * @param controlPoints
// * @returns
// */
// export const cubicBezierPoint = <Point extends LocalPoint | GlobalPoint>(
// t: number,
// controlPoints: Curve<Point>,
// ): Point => {
// const [p0, p1, p2, p3] = controlPoints;
// const x =
// Math.pow(1 - t, 3) * p0[0] +
// 3 * Math.pow(1 - t, 2) * t * p1[0] +
// 3 * (1 - t) * Math.pow(t, 2) * p2[0] +
// Math.pow(t, 3) * p3[0];
// const y =
// Math.pow(1 - t, 3) * p0[1] +
// 3 * Math.pow(1 - t, 2) * t * p1[1] +
// 3 * (1 - t) * Math.pow(t, 2) * p2[1] +
// Math.pow(t, 3) * p3[1];
// return pointFrom(x, y);
// };
// /**
// *
// * @param point
// * @param controlPoints
// * @returns
// */
// export const cubicBezierDistance = <Point extends LocalPoint | GlobalPoint>(
// point: Point,
// controlPoints: Curve<Point>,
// ) => {
// // Calculate the closest point on the Bezier curve to the given point
// const t = findClosestParameter(point, controlPoints);
// // Calculate the coordinates of the closest point on the curve
// const [closestX, closestY] = cubicBezierPoint(t, controlPoints);
// // Calculate the distance between the given point and the closest point on the curve
// const distance = Math.sqrt(
// (point[0] - closestX) ** 2 + (point[1] - closestY) ** 2,
// );
// return distance;
// };
// const solveCubic = (a: number, b: number, c: number, d: number) => {
// // This function solves the cubic equation ax^3 + bx^2 + cx + d = 0
// const roots: number[] = [];
// const discriminant =
// 18 * a * b * c * d -
// 4 * Math.pow(b, 3) * d +
// Math.pow(b, 2) * Math.pow(c, 2) -
// 4 * a * Math.pow(c, 3) -
// 27 * Math.pow(a, 2) * Math.pow(d, 2);
// if (discriminant >= 0) {
// const C = Math.cbrt((discriminant + Math.sqrt(discriminant)) / 2);
// const D = Math.cbrt((discriminant - Math.sqrt(discriminant)) / 2);
// const root1 = (-b - C - D) / (3 * a);
// const root2 = (-b + (C + D) / 2) / (3 * a);
// const root3 = (-b + (C + D) / 2) / (3 * a);
// roots.push(root1, root2, root3);
// } else {
// const realPart = -b / (3 * a);
// const root1 =
// 2 * Math.sqrt(-b / (3 * a)) * Math.cos(Math.acos(realPart) / 3);
// const root2 =
// 2 *
// Math.sqrt(-b / (3 * a)) *
// Math.cos((Math.acos(realPart) + 2 * Math.PI) / 3);
// const root3 =
// 2 *
// Math.sqrt(-b / (3 * a)) *
// Math.cos((Math.acos(realPart) + 4 * Math.PI) / 3);
// roots.push(root1, root2, root3);
// }
// return roots;
// };
// const findClosestParameter = <Point extends LocalPoint | GlobalPoint>(
// point: Point,
// controlPoints: Curve<Point>,
// ) => {
// // This function finds the parameter t that minimizes the distance between the point
// // and any point on the cubic Bezier curve.
// const [p0, p1, p2, p3] = controlPoints;
// // Use the direct formula to find the parameter t
// const a = p3[0] - 3 * p2[0] + 3 * p1[0] - p0[0];
// const b = 3 * p2[0] - 6 * p1[0] + 3 * p0[0];
// const c = 3 * p1[0] - 3 * p0[0];
// const d = p0[0] - point[0];
// const rootsX = solveCubic(a, b, c, d);
// // Do the same for the y-coordinate
// const e = p3[1] - 3 * p2[1] + 3 * p1[1] - p0[1];
// const f = 3 * p2[1] - 6 * p1[1] + 3 * p0[1];
// const g = 3 * p1[1] - 3 * p0[1];
// const h = p0[1] - point[1];
// const rootsY = solveCubic(e, f, g, h);
// // Select the real root that is between 0 and 1 (inclusive)
// const validRootsX = rootsX.filter((root) => root >= 0 && root <= 1);
// const validRootsY = rootsY.filter((root) => root >= 0 && root <= 1);
/**
* Finds the closest point on the Bezier curve from another point
*
* @param x
* @param y
* @param P0
* @param P1
* @param P2
* @param P3
* @param tolerance
* @param maxIterations
* @returns
*/
export function curveClosestPoint<Point extends GlobalPoint | LocalPoint>(
c: Curve<Point>,
p: Point,
tolerance: number = 1e-6,
maxIterations: number = 100,
): Point {
const [P0, P1, P2, P3] = c;
let t = 0.5; // Initial guess for t
for (let i = 0; i < maxIterations; i++) {
const B = [
(1 - t) ** 3 * P0[0] +
3 * (1 - t) ** 2 * t * P1[0] +
3 * (1 - t) * t ** 2 * P2[0] +
t ** 3 * P3[0],
(1 - t) ** 3 * P0[1] +
3 * (1 - t) ** 2 * t * P1[1] +
3 * (1 - t) * t ** 2 * P2[1] +
t ** 3 * P3[1],
]; // Current point on the curve
const dB = [
3 * (1 - t) ** 2 * (P1[0] - P0[0]) +
6 * (1 - t) * t * (P2[0] - P1[0]) +
3 * t ** 2 * (P3[0] - P2[0]),
3 * (1 - t) ** 2 * (P1[1] - P0[1]) +
6 * (1 - t) * t * (P2[1] - P1[1]) +
3 * t ** 2 * (P3[1] - P2[1]),
]; // Derivative at t
// Compute f(t) and f'(t)
const f = (p[0] - B[0]) * dB[0] + (p[1] - B[1]) * dB[1];
const df =
(-1 * dB[0]) ** 2 -
dB[1] ** 2 +
(p[0] - B[0]) *
(-6 * (1 - t) * (P1[0] - P0[0]) +
6 * (1 - 2 * t) * (P2[0] - P1[0]) +
6 * t * (P3[0] - P2[0])) +
(p[1] - B[1]) *
(-6 * (1 - t) * (P1[1] - P0[1]) +
6 * (1 - 2 * t) * (P2[1] - P1[1]) +
6 * t * (P3[1] - P2[1]));
// Check for convergence
if (Math.abs(f) < tolerance) {
break;
}
// if (validRootsX.length === 0 || validRootsY.length === 0) {
// // No valid roots found, use the midpoint as a fallback
// return 0.5;
// }
// Update t using Newton-Raphson
t = t - f / df;
// // Choose the parameter t that minimizes the distance
// let minDistance = Infinity;
// let closestT = 0;
// Clamp t to [0, 1] to stay within the curve segment
t = Math.max(0, Math.min(1, t));
}
// for (const rootX of validRootsX) {
// for (const rootY of validRootsY) {
// const distance = Math.sqrt(
// (rootX - point[0]) ** 2 + (rootY - point[1]) ** 2,
// );
// if (distance < minDistance) {
// minDistance = distance;
// closestT = (rootX + rootY) / 2; // Use the average for a smoother result
// }
// }
// }
// Return the closest point on the curve
return pointFrom(
(1 - t) ** 3 * P0[0] +
3 * (1 - t) ** 2 * t * P1[0] +
3 * (1 - t) * t ** 2 * P2[0] +
t ** 3 * P3[0],
(1 - t) ** 3 * P0[1] +
3 * (1 - t) ** 2 * t * P1[1] +
3 * (1 - t) * t ** 2 * P2[1] +
t ** 3 * P3[1],
);
}
// return closestT;
// };
/**
* Determines the distance between a point and the closest point on the
* Bezier curve.
*
* @param c The curve to test
* @param p The point to measure from
*/
export function curvePointDistance<Point extends GlobalPoint | LocalPoint>(
c: Curve<Point>,
p: Point,
) {
return pointDistance(p, curveClosestPoint(c, p));
}

Loading…
Cancel
Save