From e703405f680125f96b2232a9c0df7827c6d7caf3 Mon Sep 17 00:00:00 2001 From: Jeff Emmett Date: Sun, 22 Feb 2026 00:10:52 +0000 Subject: [PATCH] Move conic trajectory visualization into rTrips as route planner Relocate the conic intersection calculator from standalone rConic module into rTrips at /routes sub-page. Adds enhanced visualization with ghost conic curves, distinct arc colors, overlap zones, crossing angle lines, pulsing intersection markers, arc labels, and SVG legend. - New math: sampleConicCurve() for ghost curve tracing, conicTangentAt() for tangent/crossing angle computation - Route planner served at /{space}/trips/routes with OSRM proxy - Removed standalone conic module registration from server/index.ts - Updated vite build config to build folk-route-planner under trips Co-Authored-By: Claude Opus 4.6 --- .../trips/components/folk-route-planner.ts | 531 ++++++++++++ modules/trips/components/route-planner.css | 7 + modules/trips/lib/conic-math.ts | 763 ++++++++++++++++++ modules/trips/lib/projection.ts | 71 ++ modules/trips/lib/types.ts | 64 ++ modules/trips/mod.ts | 33 + server/index.ts | 2 - vite.config.ts | 25 +- 8 files changed, 1481 insertions(+), 15 deletions(-) create mode 100644 modules/trips/components/folk-route-planner.ts create mode 100644 modules/trips/components/route-planner.css create mode 100644 modules/trips/lib/conic-math.ts create mode 100644 modules/trips/lib/projection.ts create mode 100644 modules/trips/lib/types.ts diff --git a/modules/trips/components/folk-route-planner.ts b/modules/trips/components/folk-route-planner.ts new file mode 100644 index 0000000..65fa8e7 --- /dev/null +++ b/modules/trips/components/folk-route-planner.ts @@ -0,0 +1,531 @@ +import type { GeoPoint, FittedRoute, RouteInput, IntersectionPoint, ConicArc } from "../lib/types"; +import { createProjection, type Projection } from "../lib/projection"; +import { fitRoute, findAllIntersections, bezierToSVGCubic, sampleConicCurve, conicTangentAt } from "../lib/conic-math"; + +// Monaco demo routes that cross each other +const DEMO_ROUTES: [RouteInput, RouteInput] = [ + { + label: "Route A", + start: { lng: 7.4267, lat: 43.7396 }, // Casino Monte-Carlo + end: { lng: 7.4163, lat: 43.7275 }, // Fontvieille + }, + { + label: "Route B", + start: { lng: 7.4317, lat: 43.7519 }, // Larvotto Beach area + end: { lng: 7.4197, lat: 43.7311 }, // Monaco-Ville / Palace + }, +]; + +// Distinct shades per arc segment +const COLORS_A = ["#00dcff", "#00b8d9", "#0099b3", "#007a8c"]; +const COLORS_B = ["#ffb400", "#e6a200", "#cc9000", "#b37e00"]; + +class FolkRoutePlanner extends HTMLElement { + private shadow: ShadowRoot; + private routeA: FittedRoute | null = null; + private routeB: FittedRoute | null = null; + private intersections: IntersectionPoint[] = []; + private projection: Projection | null = null; + private loading = false; + private error = ""; + private inputA: RouteInput = DEMO_ROUTES[0]; + private inputB: RouteInput = DEMO_ROUTES[1]; + + constructor() { + super(); + this.shadow = this.attachShadow({ mode: "open" }); + } + + connectedCallback() { + this.render(); + } + + private getApiBase(): string { + const path = window.location.pathname; + const parts = path.split("/").filter(Boolean); + if (parts.length >= 2 && parts[1] === "trips") { + return `/${parts[0]}/trips`; + } + return "/demo/trips"; + } + + private async fetchRoute(input: RouteInput): Promise { + const res = await fetch(`${this.getApiBase()}/api/route`, { + method: "POST", + headers: { "Content-Type": "application/json" }, + body: JSON.stringify({ start: input.start, end: input.end }), + }); + + if (!res.ok) { + const err = await res.json().catch(() => ({ error: "Unknown error" })); + throw new Error(err.error || `OSRM returned ${res.status}`); + } + + const data = await res.json(); + if (data.code !== "Ok" || !data.routes?.length) { + throw new Error("No route found between the given points"); + } + + const route = data.routes[0]; + const rawCoords: GeoPoint[] = route.geometry.coordinates.map( + ([lng, lat]: [number, number]) => ({ lng, lat }), + ); + + return { + rawCoords, + projectedPoints: [], // filled after projection + arcs: [], + distance: route.distance, + duration: route.duration, + }; + } + + private async compute() { + this.loading = true; + this.error = ""; + this.render(); + + try { + // Fetch both routes + const [rawA, rawB] = await Promise.all([ + this.fetchRoute(this.inputA), + this.fetchRoute(this.inputB), + ]); + + // Create unified projection from all coordinates + const allGeo = [...rawA.rawCoords, ...rawB.rawCoords]; + this.projection = createProjection(allGeo); + + // Project and fit route A + rawA.projectedPoints = rawA.rawCoords.map((g) => this.projection!.toSVG(g)); + rawA.arcs = fitRoute(rawA.projectedPoints); + this.routeA = rawA; + + // Project and fit route B + rawB.projectedPoints = rawB.rawCoords.map((g) => this.projection!.toSVG(g)); + rawB.arcs = fitRoute(rawB.projectedPoints); + this.routeB = rawB; + + // Find intersections + this.intersections = findAllIntersections(this.routeA.arcs, this.routeB.arcs); + } catch (e: any) { + this.error = e.message || String(e); + } + + this.loading = false; + this.render(); + } + + private readInputs() { + const q = (sel: string) => this.shadow.querySelector(sel) as HTMLInputElement | null; + const parse = (v: string | undefined, fallback: number) => { + const n = parseFloat(v || ""); + return isNaN(n) ? fallback : n; + }; + this.inputA = { + label: "Route A", + start: { + lat: parse(q("#a-start-lat")?.value, this.inputA.start.lat), + lng: parse(q("#a-start-lng")?.value, this.inputA.start.lng), + }, + end: { + lat: parse(q("#a-end-lat")?.value, this.inputA.end.lat), + lng: parse(q("#a-end-lng")?.value, this.inputA.end.lng), + }, + }; + this.inputB = { + label: "Route B", + start: { + lat: parse(q("#b-start-lat")?.value, this.inputB.start.lat), + lng: parse(q("#b-start-lng")?.value, this.inputB.start.lng), + }, + end: { + lat: parse(q("#b-end-lat")?.value, this.inputB.end.lat), + lng: parse(q("#b-end-lng")?.value, this.inputB.end.lng), + }, + }; + } + + private renderSVG(): string { + if (!this.routeA || !this.routeB || !this.projection) return ""; + + const vb = this.projection.viewBox; + + let svg = ``; + + // Defs: grid + pulse animation + svg += ` + + + + + `; + + // Layer 1: Grid (done above) + + // Layer 2: Overlap zone rectangles + svg += this.renderOverlapZones(this.routeA.arcs, this.routeB.arcs); + + // Layer 3: Ghost conic curves (dashed, low opacity) + svg += this.renderGhostCurves(this.routeA.arcs, "A"); + svg += this.renderGhostCurves(this.routeB.arcs, "B"); + + // Layer 4: Route waypoints (subtle) + svg += this.renderWaypoints(this.routeA.projectedPoints, "rgba(0,220,255,0.2)"); + svg += this.renderWaypoints(this.routeB.projectedPoints, "rgba(255,180,0,0.2)"); + + // Layer 5: Fitted arc paths (solid, distinct shades) + svg += this.renderArcs(this.routeA.arcs, COLORS_A, "routeA"); + svg += this.renderArcs(this.routeB.arcs, COLORS_B, "routeB"); + + // Layer 6: Arc labels + svg += this.renderArcLabels(this.routeA.arcs, "A", COLORS_A); + svg += this.renderArcLabels(this.routeB.arcs, "B", COLORS_B); + + // Layer 7: Intersection crossing lines + svg += this.renderCrossingLines(); + + // Layer 8: Pulsing intersection markers + labels + for (let i = 0; i < this.intersections.length; i++) { + const pt = this.intersections[i].point; + svg += ``; + svg += ``; + svg += `X${i + 1}`; + } + + // Layer 9: Start/end markers + svg += this.renderMarker(this.routeA.projectedPoints[0], "A0", "#00dcff"); + svg += this.renderMarker(this.routeA.projectedPoints[this.routeA.projectedPoints.length - 1], "A1", "#00dcff"); + svg += this.renderMarker(this.routeB.projectedPoints[0], "B0", "#ffb400"); + svg += this.renderMarker(this.routeB.projectedPoints[this.routeB.projectedPoints.length - 1], "B1", "#ffb400"); + + // Legend + svg += this.renderLegend(vb); + + svg += ``; + return svg; + } + + private renderOverlapZones(arcsA: ConicArc[], arcsB: ConicArc[]): string { + let svg = ""; + for (const aA of arcsA) { + const bbA = this.arcBounds(aA.points); + for (const aB of arcsB) { + const bbB = this.arcBounds(aB.points); + // Compute intersection of bounding boxes + const x1 = Math.max(bbA.minX, bbB.minX); + const y1 = Math.max(bbA.minY, bbB.minY); + const x2 = Math.min(bbA.maxX, bbB.maxX); + const y2 = Math.min(bbA.maxY, bbB.maxY); + if (x2 > x1 && y2 > y1) { + const pad = 8; + svg += ``; + } + } + } + return svg; + } + + private renderGhostCurves(arcs: ConicArc[], routeId: string): string { + const colors = routeId === "A" ? COLORS_A : COLORS_B; + let svg = ""; + for (let i = 0; i < arcs.length; i++) { + const arc = arcs[i]; + if (arc.points.length < 2 || arc.type === "degenerate") continue; + + const extend = 80; + const sampled = sampleConicCurve(arc.coeffs, arc.points, extend); + if (sampled.length < 2) continue; + + let d = `M ${sampled[0].x.toFixed(1)} ${sampled[0].y.toFixed(1)}`; + for (let j = 1; j < sampled.length; j++) { + d += ` L ${sampled[j].x.toFixed(1)} ${sampled[j].y.toFixed(1)}`; + } + const color = colors[i % colors.length]; + svg += ``; + } + return svg; + } + + private renderWaypoints(pts: { x: number; y: number }[], color: string): string { + return pts.map((p) => + `` + ).join(""); + } + + private renderArcs(arcs: ConicArc[], colors: string[], id: string): string { + let svg = ""; + for (let i = 0; i < arcs.length; i++) { + const arc = arcs[i]; + if (arc.beziers.length === 0) continue; + + const first = arc.beziers[0].p0; + let d = `M ${first.x.toFixed(1)} ${first.y.toFixed(1)} `; + for (const b of arc.beziers) { + d += bezierToSVGCubic(b) + " "; + } + const color = colors[i % colors.length]; + svg += ``; + } + return svg; + } + + private renderArcLabels(arcs: ConicArc[], prefix: string, colors: string[]): string { + let svg = ""; + for (let i = 0; i < arcs.length; i++) { + const arc = arcs[i]; + if (arc.points.length < 2) continue; + // Label at the midpoint of the arc + const midIdx = Math.floor(arc.points.length / 2); + const mp = arc.points[midIdx]; + const label = `${prefix}.${i + 1} ${arc.type}`; + const color = colors[i % colors.length]; + // Background for readability + svg += ``; + svg += `${label}`; + } + return svg; + } + + private renderCrossingLines(): string { + if (!this.routeA || !this.routeB) return ""; + let svg = ""; + const lineLen = 20; + + for (let i = 0; i < this.intersections.length; i++) { + const ix = this.intersections[i]; + const pt = ix.point; + const arcA = this.routeA.arcs[ix.arcIndexA]; + const arcB = this.routeB.arcs[ix.arcIndexB]; + if (!arcA || !arcB) continue; + + const tA = conicTangentAt(arcA.coeffs, pt); + const tB = conicTangentAt(arcB.coeffs, pt); + + // Draw tangent lines + const colorA = COLORS_A[ix.arcIndexA % COLORS_A.length]; + const colorB = COLORS_B[ix.arcIndexB % COLORS_B.length]; + svg += ``; + svg += ``; + + // Compute crossing angle + const dot = Math.abs(tA.dx * tB.dx + tA.dy * tB.dy); + const angle = Math.acos(Math.min(1, dot)) * (180 / Math.PI); + svg += `${angle.toFixed(1)}\u00B0`; + } + return svg; + } + + private renderMarker(pt: { x: number; y: number }, label: string, color: string): string { + return `` + + `${label}`; + } + + private renderLegend(viewBox: string): string { + // Parse viewBox to position legend in top-right + const parts = viewBox.split(" ").map(Number); + const vx = parts[0], vy = parts[1], vw = parts[2], vh = parts[3]; + const lx = vx + vw - 155; + const ly = vy + 10; + let svg = ``; + svg += ``; + // Ghost curve + svg += ``; + svg += `Ghost conic`; + // Fitted arc + svg += ``; + svg += `Fitted arc`; + // Overlap zone + svg += ``; + svg += `Overlap zone`; + // Intersection + svg += ``; + svg += `Intersection`; + svg += ``; + return svg; + } + + private arcBounds(pts: { x: number; y: number }[]): { minX: number; maxX: number; minY: number; maxY: number } { + let minX = Infinity, maxX = -Infinity, minY = Infinity, maxY = -Infinity; + for (const p of pts) { + if (p.x < minX) minX = p.x; + if (p.x > maxX) maxX = p.x; + if (p.y < minY) minY = p.y; + if (p.y > maxY) maxY = p.y; + } + return { minX, maxX, minY, maxY }; + } + + private renderInfo(): string { + if (!this.routeA || !this.routeB) return ""; + + const fmtDist = (m: number) => m >= 1000 ? `${(m / 1000).toFixed(1)}km` : `${Math.round(m)}m`; + const fmtTime = (s: number) => { + const min = Math.floor(s / 60); + return min > 0 ? `${min}m ${Math.round(s % 60)}s` : `${Math.round(s)}s`; + }; + + let html = `
`; + + // Summary row + html += `
+ Route A: ${this.routeA.arcs.length} arcs | ${fmtDist(this.routeA.distance)} | ${fmtTime(this.routeA.duration)} + Route B: ${this.routeB.arcs.length} arcs | ${fmtDist(this.routeB.distance)} | ${fmtTime(this.routeB.duration)} + ${this.intersections.length} intersection${this.intersections.length !== 1 ? "s" : ""} +
`; + + // Arc details + html += `

Route A Arcs

`; + for (let i = 0; i < this.routeA.arcs.length; i++) { + const arc = this.routeA.arcs[i]; + html += `
A.${i + 1}: ${arc.type} e=${arc.eccentricity.toFixed(3)}
`; + } + html += `
`; + + html += `

Route B Arcs

`; + for (let i = 0; i < this.routeB.arcs.length; i++) { + const arc = this.routeB.arcs[i]; + html += `
B.${i + 1}: ${arc.type} e=${arc.eccentricity.toFixed(3)}
`; + } + html += `
`; + + // Intersection details + if (this.intersections.length > 0 && this.projection) { + html += `

Intersections

`; + for (let i = 0; i < this.intersections.length; i++) { + const ix = this.intersections[i]; + const geo = this.projection.toGeo(ix.point); + html += `
X${i + 1}: ${geo.lat.toFixed(5)}, ${geo.lng.toFixed(5)} (A.${ix.arcIndexA + 1} / B.${ix.arcIndexB + 1})
`; + } + html += `
`; + } + + html += `
`; + return html; + } + + private render() { + this.shadow.innerHTML = ` + +
+
+

rTrips — Route Planner

+

Plan paths between destinations. Fit conic arcs and find where routes intersect.

+
+ +
+
+ + + + to + + +
+
+ + + + to + + +
+
+ + +
+
+ + ${this.error ? `
${this.error}
` : ""} + ${this.loading ? `
Fetching routes and computing conic fits...
` : ""} + + ${this.routeA && this.routeB ? ` +
${this.renderSVG()}
+ ${this.renderInfo()} +
+ Conic: Ax\u00B2+Bxy+Cy\u00B2+Dx+Ey+F=0 | Discriminant \u0394=B\u00B2\u22124AC | + Intersection via Sylvester resultant \u2192 degree-4 polynomial \u2192 companion matrix eigenvalues +
+ ` : ""} + + ${!this.routeA && !this.routeB && !this.loading && !this.error ? ` +
+ Enter start/end coordinates for two routes, or click Use Demo to load Monaco routes. +
+ ` : ""} +
`; + + // Attach event listeners + this.shadow.querySelector("#btn-compute")?.addEventListener("click", () => { + this.readInputs(); + this.compute(); + }); + + this.shadow.querySelector("#btn-demo")?.addEventListener("click", () => { + this.inputA = DEMO_ROUTES[0]; + this.inputB = DEMO_ROUTES[1]; + this.compute(); + }); + } + + static define() { + if (!customElements.get("folk-route-planner")) { + customElements.define("folk-route-planner", FolkRoutePlanner); + } + } +} + +FolkRoutePlanner.define(); + +export { FolkRoutePlanner }; diff --git a/modules/trips/components/route-planner.css b/modules/trips/components/route-planner.css new file mode 100644 index 0000000..4ebd6bd --- /dev/null +++ b/modules/trips/components/route-planner.css @@ -0,0 +1,7 @@ +folk-route-planner { + display: block; + min-height: 400px; + padding: 20px; + max-width: 1200px; + margin: 0 auto; +} diff --git a/modules/trips/lib/conic-math.ts b/modules/trips/lib/conic-math.ts new file mode 100644 index 0000000..dfd81f0 --- /dev/null +++ b/modules/trips/lib/conic-math.ts @@ -0,0 +1,763 @@ +import type { Point2D, ConicCoeffs, ConicType, ConicArc, RationalQuadBezier, IntersectionPoint } from "./types"; + +const EPS = 1e-10; + +// ──────────────────────────────────────────────────────── +// Conic fitting via null space of design matrix +// ──────────────────────────────────────────────────────── + +/** + * Fit a general conic Ax²+Bxy+Cy²+Dx+Ey+F=0 through a set of points. + * Uses SVD-like approach: find the eigenvector of M^T*M with smallest eigenvalue. + * Points are normalized to [-1,1] before fitting for numerical stability. + */ +export function fitConic(points: Point2D[]): ConicCoeffs { + if (points.length < 5) { + return { A: 0, B: 0, C: 0, D: 0, E: 0, F: 1 }; // degenerate + } + + // Normalize coordinates to [-1,1] + let mx = 0, my = 0; + for (const p of points) { mx += p.x; my += p.y; } + mx /= points.length; my /= points.length; + let sx = 0, sy = 0; + for (const p of points) { + sx = Math.max(sx, Math.abs(p.x - mx)); + sy = Math.max(sy, Math.abs(p.y - my)); + } + sx = sx || 1; sy = sy || 1; + + // Build design matrix D (n x 6): each row = [x², xy, y², x, y, 1] + const n = points.length; + const D: number[][] = []; + for (const p of points) { + const x = (p.x - mx) / sx; + const y = (p.y - my) / sy; + D.push([x * x, x * y, y * y, x, y, 1]); + } + + // Compute M = D^T * D (6x6) + const M = Array.from({ length: 6 }, () => new Float64Array(6)); + for (let i = 0; i < 6; i++) { + for (let j = i; j < 6; j++) { + let sum = 0; + for (let k = 0; k < n; k++) sum += D[k][i] * D[k][j]; + M[i][j] = sum; + M[j][i] = sum; + } + } + + // Find eigenvector with smallest eigenvalue via inverse power iteration + const coeffs = smallestEigenvector6(M); + + // Denormalize: transform coefficients back to original coordinate space + // Substituting x_norm = (x - mx)/sx, y_norm = (y - my)/sy into the conic equation + const [a, b, c, d, e, f] = coeffs; + const A = a / (sx * sx); + const B = b / (sx * sy); + const C = c / (sy * sy); + const D2 = d / sx - 2 * a * mx / (sx * sx) - b * my / (sx * sy); + const E2 = e / sy - 2 * c * my / (sy * sy) - b * mx / (sx * sy); + const F2 = a * mx * mx / (sx * sx) + b * mx * my / (sx * sy) + c * my * my / (sy * sy) + - d * mx / sx - e * my / sy + f; + + return { A, B, C, D: D2, E: E2, F: F2 }; +} + +/** + * Find the eigenvector corresponding to the smallest eigenvalue of a 6x6 symmetric matrix. + * Uses inverse power iteration with a small shift. + */ +function smallestEigenvector6(M: Float64Array[]): number[] { + // Estimate smallest eigenvalue with Gershgorin circle theorem lower bound + let minEig = Infinity; + for (let i = 0; i < 6; i++) { + let r = 0; + for (let j = 0; j < 6; j++) if (i !== j) r += Math.abs(M[i][j]); + minEig = Math.min(minEig, M[i][i] - r); + } + const shift = minEig - 0.001; + + // M_shifted = M - shift*I + const Ms = Array.from({ length: 6 }, (_, i) => { + const row = new Float64Array(6); + for (let j = 0; j < 6; j++) row[j] = M[i][j]; + row[i] -= shift; + return row; + }); + + // LU decomposition of Ms + const LU = Ms.map((r) => Float64Array.from(r)); + const piv = Array.from({ length: 6 }, (_, i) => i); + for (let k = 0; k < 6; k++) { + let maxVal = 0, maxRow = k; + for (let i = k; i < 6; i++) { + if (Math.abs(LU[i][k]) > maxVal) { maxVal = Math.abs(LU[i][k]); maxRow = i; } + } + if (maxRow !== k) { + [LU[k], LU[maxRow]] = [LU[maxRow], LU[k]]; + [piv[k], piv[maxRow]] = [piv[maxRow], piv[k]]; + } + if (Math.abs(LU[k][k]) < 1e-15) LU[k][k] = 1e-15; + for (let i = k + 1; i < 6; i++) { + LU[i][k] /= LU[k][k]; + for (let j = k + 1; j < 6; j++) LU[i][j] -= LU[i][k] * LU[k][j]; + } + } + + function solveLU(b: number[]): number[] { + const pb = piv.map((i) => b[i]); + // Forward substitution + for (let i = 1; i < 6; i++) { + for (let j = 0; j < i; j++) pb[i] -= LU[i][j] * pb[j]; + } + // Back substitution + for (let i = 5; i >= 0; i--) { + for (let j = i + 1; j < 6; j++) pb[i] -= LU[i][j] * pb[j]; + pb[i] /= LU[i][i]; + } + return pb; + } + + // Power iteration on (M - shift*I)^{-1} + let v = [1, 1, 1, 1, 1, 1]; + for (let iter = 0; iter < 30; iter++) { + v = solveLU(v); + let norm = 0; + for (const x of v) norm += x * x; + norm = Math.sqrt(norm) || 1; + for (let i = 0; i < 6; i++) v[i] /= norm; + } + return v; +} + +// ──────────────────────────────────────────────────────── +// Conic classification +// ──────────────────────────────────────────────────────── + +export function classifyConic(c: ConicCoeffs): { type: ConicType; eccentricity: number } { + const delta = c.B * c.B - 4 * c.A * c.C; + + // 3x3 determinant for degeneracy check + const det3 = c.A * (c.C * c.F - c.E * c.E / 4) + - c.B / 2 * (c.B * c.F / 2 - c.E * c.D / 4) + + c.D / 2 * (c.B * c.E / 4 - c.C * c.D / 2); + + if (Math.abs(det3) < EPS) return { type: "degenerate", eccentricity: 0 }; + + if (Math.abs(delta) < EPS) { + return { type: "parabola", eccentricity: 1 }; + } + + if (delta < 0) { + // Ellipse or circle + if (Math.abs(c.A - c.C) < EPS && Math.abs(c.B) < EPS) { + return { type: "circle", eccentricity: 0 }; + } + // Eccentricity from eigenvalues of [[A, B/2], [B/2, C]] + const tr = c.A + c.C; + const det = c.A * c.C - (c.B / 2) ** 2; + const disc = Math.sqrt(Math.max(0, tr * tr - 4 * det)); + const lam1 = (tr + disc) / 2; + const lam2 = (tr - disc) / 2; + const ratio = Math.min(Math.abs(lam1), Math.abs(lam2)) / Math.max(Math.abs(lam1), Math.abs(lam2)); + const e = Math.sqrt(Math.max(0, 1 - ratio)); + return { type: "ellipse", eccentricity: e }; + } + + // Hyperbola + const tr = c.A + c.C; + const det = c.A * c.C - (c.B / 2) ** 2; + const disc = Math.sqrt(Math.max(0, tr * tr - 4 * det)); + const lam1 = (tr + disc) / 2; + const lam2 = (tr - disc) / 2; + const absRatio = Math.abs(lam1 / (lam2 || 1)); + const e = Math.sqrt(1 + absRatio); + return { type: "hyperbola", eccentricity: Math.min(e, 10) }; // cap for display +} + +// ──────────────────────────────────────────────────────── +// Route segmentation +// ──────────────────────────────────────────────────────── + +/** + * Split route points into overlapping segments for conic fitting. + * Each segment has maxPerArc points with 2-point overlap. + */ +export function segmentRoute(points: Point2D[], maxPerArc = 10): Point2D[][] { + if (points.length <= maxPerArc) return [points]; + + const segments: Point2D[][] = []; + const step = maxPerArc - 2; // overlap of 2 + for (let i = 0; i < points.length - 4; i += step) { + const end = Math.min(i + maxPerArc, points.length); + segments.push(points.slice(i, end)); + if (end >= points.length) break; + } + return segments; +} + +// ──────────────────────────────────────────────────────── +// Conic arc → Bézier conversion for SVG rendering +// ──────────────────────────────────────────────────────── + +/** + * Convert a conic arc (defined by its fitted points) into piecewise rational quadratic Béziers. + * For each consecutive pair of points, computes a Bézier arc using the conic gradient for tangents. + */ +export function conicArcToBeziers(coeffs: ConicCoeffs, points: Point2D[]): RationalQuadBezier[] { + if (points.length < 2) return []; + + const beziers: RationalQuadBezier[] = []; + for (let i = 0; i < points.length - 1; i++) { + const p0 = points[i]; + const p2 = points[i + 1]; + + // Tangent at each point: gradient of Ax²+Bxy+Cy²+Dx+Ey+F + // ∇F = (2Ax+By+D, Bx+2Cy+E), tangent is perpendicular to gradient + const g0x = 2 * coeffs.A * p0.x + coeffs.B * p0.y + coeffs.D; + const g0y = coeffs.B * p0.x + 2 * coeffs.C * p0.y + coeffs.E; + const g2x = 2 * coeffs.A * p2.x + coeffs.B * p2.y + coeffs.D; + const g2y = coeffs.B * p2.x + 2 * coeffs.C * p2.y + coeffs.E; + + // Tangent = perpendicular to gradient: (-gy, gx) + const t0x = -g0y, t0y = g0x; + const t2x = -g2y, t2y = g2x; + + // Control point = intersection of tangent lines from p0 and p2 + const p1 = intersectLines(p0.x, p0.y, t0x, t0y, p2.x, p2.y, t2x, t2y); + + if (p1) { + // Weight from the conic — use the midpoint deviation + const mid = { x: (p0.x + p2.x) / 2, y: (p0.y + p2.y) / 2 }; + const d0 = Math.hypot(p1.x - mid.x, p1.y - mid.y); + const d1 = Math.hypot(p2.x - p0.x, p2.y - p0.y); + const w = d1 > EPS ? Math.min(Math.max(d1 / (2 * d0 + d1), 0.1), 5) : 1; + + beziers.push({ p0, p1, p2, w }); + } else { + // Tangent lines are parallel → straight segment + beziers.push({ + p0, p2, + p1: { x: (p0.x + p2.x) / 2, y: (p0.y + p2.y) / 2 }, + w: 1, + }); + } + } + return beziers; +} + +/** Intersect two lines: (px, py) + t*(dx, dy) and (qx, qy) + s*(ex, ey) */ +function intersectLines( + px: number, py: number, dx: number, dy: number, + qx: number, qy: number, ex: number, ey: number, +): Point2D | null { + const denom = dx * ey - dy * ex; + if (Math.abs(denom) < EPS) return null; + const t = ((qx - px) * ey - (qy - py) * ex) / denom; + return { x: px + t * dx, y: py + t * dy }; +} + +// ──────────────────────────────────────────────────────── +// Conic intersection via Sylvester resultant +// ──────────────────────────────────────────────────────── + +/** + * Find intersection points of two general conics. + * Uses the resultant method: eliminate x to get a degree-4 polynomial in y, + * then solve the quartic via companion matrix eigenvalues. + */ +export function findIntersections(c1: ConicCoeffs, c2: ConicCoeffs): Point2D[] { + // Each conic: A*x² + (B*y + D)*x + (C*y² + E*y + F) = 0 + // Treat as quadratics in x with y-dependent coefficients: + // a_i(y)*x² + b_i(y)*x + c_i(y) = 0 + // a₁ = A₁, b₁(y) = B₁y + D₁, c₁(y) = C₁y² + E₁y + F₁ + // Resultant eliminates x, giving poly in y of degree ≤ 4 + + // Build the resultant polynomial coefficients in y + // Res(f,g) = a₁²c₂² + a₂²c₁² - a₁a₂(b₁c₂+b₂c₁) + b₁²a₂c₂ + ... (Sylvester det) + // It's cleaner to use the 4x4 Sylvester matrix determinant directly. + + // Sylvester matrix for two quadratics in x: + // | a₁ b₁(y) c₁(y) 0 | + // | 0 a₁ b₁(y) c₁(y)| + // | a₂ b₂(y) c₂(y) 0 | + // | 0 a₂ b₂(y) c₂(y)| + // + // Each entry is a polynomial in y. The determinant is a polynomial in y of degree ≤ 4. + // We evaluate at 5 values of y and interpolate to find the coefficients. + + const yVals = [-2, -1, 0, 1, 2]; + const detVals = yVals.map((y) => { + const b1 = c1.B * y + c1.D; + const cc1 = c1.C * y * y + c1.E * y + c1.F; + const b2 = c2.B * y + c2.D; + const cc2 = c2.C * y * y + c2.E * y + c2.F; + + // 4x4 determinant + const m = [ + [c1.A, b1, cc1, 0], + [0, c1.A, b1, cc1], + [c2.A, b2, cc2, 0], + [0, c2.A, b2, cc2], + ]; + return det4x4(m); + }); + + // Interpolate: we have 5 points → degree-4 polynomial + const polyCoeffs = lagrangeInterpolate(yVals, detVals); + + // Find real roots of the degree-4 polynomial + const yRoots = findRealRoots(polyCoeffs); + + // For each y root, solve for x from the original conics + const results: Point2D[] = []; + for (const y of yRoots) { + const xSolutions = solveForX(c1, c2, y); + for (const x of xSolutions) { + // Verify solution satisfies both conics (within tolerance) + const v1 = evalConic(c1, x, y); + const v2 = evalConic(c2, x, y); + const maxCoeff = Math.max( + Math.abs(c1.A) + Math.abs(c1.C) + Math.abs(c1.F), + Math.abs(c2.A) + Math.abs(c2.C) + Math.abs(c2.F), + 1, + ); + if (Math.abs(v1) < maxCoeff * 0.1 && Math.abs(v2) < maxCoeff * 0.1) { + results.push({ x, y }); + } + } + } + + return deduplicatePoints(results); +} + +function evalConic(c: ConicCoeffs, x: number, y: number): number { + return c.A * x * x + c.B * x * y + c.C * y * y + c.D * x + c.E * y + c.F; +} + +function det4x4(m: number[][]): number { + const [a, b, c, d] = m; + return ( + a[0] * (b[1] * (c[2] * d[3] - c[3] * d[2]) - b[2] * (c[1] * d[3] - c[3] * d[1]) + b[3] * (c[1] * d[2] - c[2] * d[1])) + - a[1] * (b[0] * (c[2] * d[3] - c[3] * d[2]) - b[2] * (c[0] * d[3] - c[3] * d[0]) + b[3] * (c[0] * d[2] - c[2] * d[0])) + + a[2] * (b[0] * (c[1] * d[3] - c[3] * d[1]) - b[1] * (c[0] * d[3] - c[3] * d[0]) + b[3] * (c[0] * d[1] - c[1] * d[0])) + - a[3] * (b[0] * (c[1] * d[2] - c[2] * d[1]) - b[1] * (c[0] * d[2] - c[2] * d[0]) + b[2] * (c[0] * d[1] - c[1] * d[0])) + ); +} + +/** + * Lagrange interpolation: given (x_i, y_i) pairs, return polynomial coefficients [a0, a1, ..., an] + * such that p(x) = a0 + a1*x + a2*x² + ... + */ +function lagrangeInterpolate(xs: number[], ys: number[]): number[] { + const n = xs.length; + const coeffs = new Array(n).fill(0); + + for (let i = 0; i < n; i++) { + // Build the i-th Lagrange basis polynomial + const basis = [ys[i]]; + for (let j = 0; j < n; j++) { + if (j === i) continue; + const scale = 1 / (xs[i] - xs[j]); + const newBasis = new Array(basis.length + 1).fill(0); + for (let k = 0; k < basis.length; k++) { + newBasis[k] += basis[k] * (-xs[j]) * scale; + newBasis[k + 1] += basis[k] * scale; + } + basis.length = newBasis.length; + for (let k = 0; k < newBasis.length; k++) basis[k] = newBasis[k]; + } + for (let k = 0; k < basis.length && k < n; k++) coeffs[k] += basis[k]; + } + return coeffs; +} + +/** + * Find real roots of a polynomial a0 + a1*x + a2*x² + ... + an*x^n = 0 + * Uses companion matrix eigenvalues via QR iteration. + */ +function findRealRoots(coeffs: number[]): number[] { + // Trim trailing near-zero coefficients + while (coeffs.length > 1 && Math.abs(coeffs[coeffs.length - 1]) < EPS) { + coeffs.pop(); + } + const deg = coeffs.length - 1; + if (deg <= 0) return []; + if (deg === 1) return [-coeffs[0] / coeffs[1]]; + if (deg === 2) return solveQuadratic(coeffs[2], coeffs[1], coeffs[0]); + + // Build companion matrix + const n = deg; + const C = Array.from({ length: n }, () => new Float64Array(n)); + const lead = coeffs[n]; + for (let i = 0; i < n; i++) { + C[i][n - 1] = -coeffs[i] / lead; + } + for (let i = 1; i < n; i++) { + C[i][i - 1] = 1; + } + + // QR iteration to find eigenvalues + return qrEigenvalues(C); +} + +/** Solve ax²+bx+c=0 for real roots */ +function solveQuadratic(a: number, b: number, c: number): number[] { + if (Math.abs(a) < EPS) { + if (Math.abs(b) < EPS) return []; + return [-c / b]; + } + const disc = b * b - 4 * a * c; + if (disc < -EPS) return []; + if (disc < EPS) return [-b / (2 * a)]; + const sq = Math.sqrt(disc); + return [(-b + sq) / (2 * a), (-b - sq) / (2 * a)]; +} + +/** + * QR iteration for eigenvalues of an n×n matrix (n ≤ 4). + * Returns only real eigenvalues. + */ +function qrEigenvalues(A: Float64Array[]): number[] { + const n = A.length; + const H = A.map((r) => Float64Array.from(r)); + + // Hessenberg reduction first + for (let k = 0; k < n - 2; k++) { + // Find the pivot + let maxVal = 0, maxIdx = k + 1; + for (let i = k + 1; i < n; i++) { + if (Math.abs(H[i][k]) > maxVal) { maxVal = Math.abs(H[i][k]); maxIdx = i; } + } + if (maxVal < EPS) continue; + if (maxIdx !== k + 1) { + [H[k + 1], H[maxIdx]] = [H[maxIdx], H[k + 1]]; + for (let i = 0; i < n; i++) { + const tmp = H[i][k + 1]; H[i][k + 1] = H[i][maxIdx]; H[i][maxIdx] = tmp; + } + } + + for (let i = k + 2; i < n; i++) { + const factor = H[i][k] / H[k + 1][k]; + for (let j = k; j < n; j++) H[i][j] -= factor * H[k + 1][j]; + for (let j = 0; j < n; j++) H[j][k + 1] += factor * H[j][i]; + } + } + + // QR iteration with Wilkinson shift + const eigenvalues: number[] = []; + let m = n; + for (let maxIter = 0; maxIter < 200 && m > 0; maxIter++) { + if (m === 1) { + eigenvalues.push(H[0][0]); + break; + } + + // Check for convergence on sub-diagonal + if (Math.abs(H[m - 1][m - 2]) < EPS * (Math.abs(H[m - 1][m - 1]) + Math.abs(H[m - 2][m - 2]) + 1)) { + eigenvalues.push(H[m - 1][m - 1]); + m--; + continue; + } + + // Wilkinson shift: eigenvalue of bottom-right 2x2 closer to H[m-1][m-1] + const a11 = H[m - 2][m - 2], a12 = H[m - 2][m - 1]; + const a21 = H[m - 1][m - 2], a22 = H[m - 1][m - 1]; + const tr = a11 + a22; + const det = a11 * a22 - a12 * a21; + const disc = tr * tr - 4 * det; + let shift: number; + if (disc >= 0) { + const sq = Math.sqrt(disc); + const e1 = (tr + sq) / 2, e2 = (tr - sq) / 2; + shift = Math.abs(e1 - a22) < Math.abs(e2 - a22) ? e1 : e2; + } else { + shift = tr / 2; + } + + // QR step: (H - shift*I) = Q*R, then H = R*Q + shift*I + // Implemented via Givens rotations + for (let i = 0; i < m; i++) H[i][i] -= shift; + + const cs = new Float64Array(m - 1); + const sn = new Float64Array(m - 1); + for (let i = 0; i < m - 1; i++) { + const a = H[i][i], b = H[i + 1][i]; + const r = Math.hypot(a, b); + if (r < EPS) { cs[i] = 1; sn[i] = 0; continue; } + cs[i] = a / r; sn[i] = b / r; + // Apply rotation to rows i, i+1 + for (let j = 0; j < m; j++) { + const t1 = H[i][j], t2 = H[i + 1][j]; + H[i][j] = cs[i] * t1 + sn[i] * t2; + H[i + 1][j] = -sn[i] * t1 + cs[i] * t2; + } + } + // Apply Q^T from the right + for (let i = 0; i < m - 1; i++) { + for (let j = 0; j < m; j++) { + const t1 = H[j][i], t2 = H[j][i + 1]; + H[j][i] = cs[i] * t1 + sn[i] * t2; + H[j][i + 1] = -sn[i] * t1 + cs[i] * t2; + } + } + + for (let i = 0; i < m; i++) H[i][i] += shift; + } + + return eigenvalues.filter((v) => isFinite(v)); +} + +/** Given a y value, solve for x from both conics and return candidate x values */ +function solveForX(c1: ConicCoeffs, c2: ConicCoeffs, y: number): number[] { + // From conic 1: A₁x² + (B₁y+D₁)x + (C₁y²+E₁y+F₁) = 0 + const a1 = c1.A; + const b1 = c1.B * y + c1.D; + const cc1 = c1.C * y * y + c1.E * y + c1.F; + const roots1 = solveQuadratic(a1, b1, cc1); + + const a2 = c2.A; + const b2 = c2.B * y + c2.D; + const cc2 = c2.C * y * y + c2.E * y + c2.F; + const roots2 = solveQuadratic(a2, b2, cc2); + + // Return x values that appear in both solutions (within tolerance) + const results: number[] = []; + for (const x1 of roots1) { + for (const x2 of roots2) { + if (Math.abs(x1 - x2) < Math.max(Math.abs(x1), 1) * 0.01) { + results.push((x1 + x2) / 2); + } + } + } + // If no matching roots, return all roots from conic 1 (the resultant guarantees they exist) + if (results.length === 0) return roots1; + return results; +} + +function deduplicatePoints(pts: Point2D[], tol = 1): Point2D[] { + const out: Point2D[] = []; + for (const p of pts) { + if (!out.some((q) => Math.hypot(q.x - p.x, q.y - p.y) < tol)) { + out.push(p); + } + } + return out; +} + +// ──────────────────────────────────────────────────────── +// High-level: process a route into conic arcs +// ──────────────────────────────────────────────────────── + +export function fitRoute(points: Point2D[]): ConicArc[] { + const segments = segmentRoute(points); + return segments.map((seg) => { + const coeffs = fitConic(seg); + const { type, eccentricity } = classifyConic(coeffs); + const beziers = conicArcToBeziers(coeffs, seg); + return { coeffs, type, eccentricity, points: seg, beziers }; + }); +} + +/** + * Find all intersection points between arcs of two fitted routes. + * Filters to points that lie within the bounding box of both arcs. + */ +export function findAllIntersections(arcsA: ConicArc[], arcsB: ConicArc[]): IntersectionPoint[] { + const results: IntersectionPoint[] = []; + + for (let ia = 0; ia < arcsA.length; ia++) { + for (let ib = 0; ib < arcsB.length; ib++) { + // Quick bounding-box overlap check + const bbA = arcBounds(arcsA[ia].points); + const bbB = arcBounds(arcsB[ib].points); + if (!boxesOverlap(bbA, bbB)) continue; + + const pts = findIntersections(arcsA[ia].coeffs, arcsB[ib].coeffs); + for (const pt of pts) { + // Filter: point must be within expanded bounding boxes of both arcs + if (inBox(pt, bbA, 20) && inBox(pt, bbB, 20)) { + results.push({ point: pt, arcIndexA: ia, arcIndexB: ib }); + } + } + } + } + + // Deduplicate across arc pairs + const unique: IntersectionPoint[] = []; + for (const r of results) { + if (!unique.some((u) => Math.hypot(u.point.x - r.point.x, u.point.y - r.point.y) < 3)) { + unique.push(r); + } + } + return unique; +} + +interface BBox { minX: number; maxX: number; minY: number; maxY: number } + +function arcBounds(pts: Point2D[]): BBox { + let minX = Infinity, maxX = -Infinity, minY = Infinity, maxY = -Infinity; + for (const p of pts) { + if (p.x < minX) minX = p.x; + if (p.x > maxX) maxX = p.x; + if (p.y < minY) minY = p.y; + if (p.y > maxY) maxY = p.y; + } + return { minX, maxX, minY, maxY }; +} + +function boxesOverlap(a: BBox, b: BBox): boolean { + return a.minX <= b.maxX && a.maxX >= b.minX && a.minY <= b.maxY && a.maxY >= b.minY; +} + +function inBox(p: Point2D, bb: BBox, margin: number): boolean { + return p.x >= bb.minX - margin && p.x <= bb.maxX + margin + && p.y >= bb.minY - margin && p.y <= bb.maxY + margin; +} + +// ──────────────────────────────────────────────────────── +// Ghost conic curve sampling & tangent +// ──────────────────────────────────────────────────────── + +/** + * Return the tangent direction at a point on the conic (perpendicular to gradient). + * Gradient of F = (2Ax+By+D, Bx+2Cy+E), tangent = (-∂F/∂y, ∂F/∂x). + */ +export function conicTangentAt(coeffs: ConicCoeffs, pt: Point2D): { dx: number; dy: number } { + const gx = 2 * coeffs.A * pt.x + coeffs.B * pt.y + coeffs.D; + const gy = coeffs.B * pt.x + 2 * coeffs.C * pt.y + coeffs.E; + const len = Math.hypot(gx, gy) || 1; + return { dx: -gy / len, dy: gx / len }; +} + +/** + * Sample the full conic curve extending `extend` units beyond the arc endpoints. + * Marches along the tangent direction, projecting back onto the conic via Newton's method. + * Returns the extended point array (backward extension + original arc + forward extension). + */ +export function sampleConicCurve(coeffs: ConicCoeffs, arcPoints: Point2D[], extend: number): Point2D[] { + if (arcPoints.length < 2) return [...arcPoints]; + + const { type } = classifyConic(coeffs); + if (type === "degenerate") return [...arcPoints]; + + // Clamp extension to 50% of arc length, minimum 20 + const arcLen = arcPoints.reduce((sum, p, i) => + i === 0 ? 0 : sum + Math.hypot(p.x - arcPoints[i - 1].x, p.y - arcPoints[i - 1].y), 0); + const maxExt = Math.max(20, Math.min(extend, arcLen * 0.5)); + const step = Math.max(2, maxExt / 25); + + const backward = marchConic(coeffs, arcPoints[0], arcPoints[1], -step, maxExt); + const forward = marchConic(coeffs, arcPoints[arcPoints.length - 1], arcPoints[arcPoints.length - 2], -step, maxExt); + + // backward is in reverse order, forward is in forward order (away from arc end) + return [...backward.reverse(), ...arcPoints, ...forward]; +} + +/** + * March along a conic curve from `start` in the direction away from `prev`. + * At each step, advance along the tangent then project back onto the conic. + */ +function marchConic(coeffs: ConicCoeffs, start: Point2D, prev: Point2D, stepSize: number, maxDist: number): Point2D[] { + const points: Point2D[] = []; + let cur = { ...start }; + + // Initial direction: away from prev + let dirX = start.x - prev.x; + let dirY = start.y - prev.y; + let dLen = Math.hypot(dirX, dirY) || 1; + dirX /= dLen; dirY /= dLen; + + const absStep = Math.abs(stepSize); + let traveled = 0; + + for (let i = 0; i < 30 && traveled < maxDist; i++) { + // Advance along tangent direction + const t = conicTangentAt(coeffs, cur); + // Choose tangent direction consistent with our travel direction + const dot = t.dx * dirX + t.dy * dirY; + const sign = dot >= 0 ? 1 : -1; + let nx = cur.x + sign * t.dx * absStep; + let ny = cur.y + sign * t.dy * absStep; + + // Newton projection: push (nx,ny) back onto the conic + for (let iter = 0; iter < 8; iter++) { + const val = evalConicAt(coeffs, nx, ny); + const gx = 2 * coeffs.A * nx + coeffs.B * ny + coeffs.D; + const gy = coeffs.B * nx + 2 * coeffs.C * ny + coeffs.E; + const g2 = gx * gx + gy * gy; + if (g2 < 1e-20) break; + const lambda = val / g2; + nx -= lambda * gx; + ny -= lambda * gy; + if (Math.abs(val) < 1e-6) break; + } + + // Bail if projection diverged + if (!isFinite(nx) || !isFinite(ny)) break; + if (Math.hypot(nx - cur.x, ny - cur.y) > absStep * 3) break; + + dirX = nx - cur.x; dirY = ny - cur.y; + dLen = Math.hypot(dirX, dirY) || 1; + dirX /= dLen; dirY /= dLen; + + cur = { x: nx, y: ny }; + traveled += absStep; + points.push({ ...cur }); + } + + return points; +} + +function evalConicAt(c: ConicCoeffs, x: number, y: number): number { + return c.A * x * x + c.B * x * y + c.C * y * y + c.D * x + c.E * y + c.F; +} + +// ──────────────────────────────────────────────────────── +// Bézier → SVG path conversion +// ──────────────────────────────────────────────────────── + +/** + * Convert a rational quadratic Bézier to an SVG cubic Bézier path command. + * Degree elevation: rational quad → approximated cubic. + */ +export function bezierToSVGCubic(b: RationalQuadBezier): string { + // Approximate rational quadratic as cubic via degree elevation + // For w ≈ 1, this is nearly exact. For w far from 1, sample instead. + const w = b.w; + if (w < 0.3 || w > 3) { + // Sample the rational curve and emit as polyline + return sampleBezierPath(b, 16); + } + + const c1x = b.p0.x + (2 * w * (b.p1.x - b.p0.x)) / (1 + w) * (2 / 3); + const c1y = b.p0.y + (2 * w * (b.p1.y - b.p0.y)) / (1 + w) * (2 / 3); + const c2x = b.p2.x + (2 * w * (b.p1.x - b.p2.x)) / (1 + w) * (2 / 3); + const c2y = b.p2.y + (2 * w * (b.p1.y - b.p2.y)) / (1 + w) * (2 / 3); + + return `C ${c1x.toFixed(1)} ${c1y.toFixed(1)} ${c2x.toFixed(1)} ${c2y.toFixed(1)} ${b.p2.x.toFixed(1)} ${b.p2.y.toFixed(1)}`; +} + +/** Sample a rational quadratic Bézier and return SVG line-to commands */ +function sampleBezierPath(b: RationalQuadBezier, samples: number): string { + const parts: string[] = []; + for (let i = 1; i <= samples; i++) { + const t = i / samples; + const pt = evalRationalQuadBezier(b, t); + parts.push(`L ${pt.x.toFixed(1)} ${pt.y.toFixed(1)}`); + } + return parts.join(" "); +} + +function evalRationalQuadBezier(b: RationalQuadBezier, t: number): Point2D { + const omt = 1 - t; + const w0 = omt * omt; + const w1 = 2 * b.w * t * omt; + const w2 = t * t; + const denom = w0 + w1 + w2; + return { + x: (w0 * b.p0.x + w1 * b.p1.x + w2 * b.p2.x) / denom, + y: (w0 * b.p0.y + w1 * b.p1.y + w2 * b.p2.y) / denom, + }; +} diff --git a/modules/trips/lib/projection.ts b/modules/trips/lib/projection.ts new file mode 100644 index 0000000..46cad8a --- /dev/null +++ b/modules/trips/lib/projection.ts @@ -0,0 +1,71 @@ +import type { GeoPoint, Point2D } from "./types"; + +const DEG_TO_RAD = Math.PI / 180; +const METERS_PER_DEG = 111_320; // at equator + +export interface Projection { + toSVG(geo: GeoPoint): Point2D; + toGeo(pt: Point2D): GeoPoint; + viewBox: string; +} + +/** + * Create an equirectangular projection centered on the centroid of all points. + * Scales so the bounding box fits ~800 SVG units wide with padding. + */ +export function createProjection(allPoints: GeoPoint[], targetWidth = 800, padding = 40): Projection { + if (allPoints.length === 0) { + return { toSVG: () => ({ x: 0, y: 0 }), toGeo: () => ({ lng: 0, lat: 0 }), viewBox: "0 0 800 600" }; + } + + // Centroid + const refLat = allPoints.reduce((s, p) => s + p.lat, 0) / allPoints.length; + const refLng = allPoints.reduce((s, p) => s + p.lng, 0) / allPoints.length; + const cosLat = Math.cos(refLat * DEG_TO_RAD); + + // Project all points to find bounds (meters from centroid) + const projected = allPoints.map((p) => ({ + x: (p.lng - refLng) * cosLat * METERS_PER_DEG, + y: -(p.lat - refLat) * METERS_PER_DEG, // SVG y-down + })); + + let minX = Infinity, maxX = -Infinity, minY = Infinity, maxY = -Infinity; + for (const p of projected) { + if (p.x < minX) minX = p.x; + if (p.x > maxX) maxX = p.x; + if (p.y < minY) minY = p.y; + if (p.y > maxY) maxY = p.y; + } + + const rangeX = maxX - minX || 1; + const rangeY = maxY - minY || 1; + const scale = (targetWidth - 2 * padding) / Math.max(rangeX, rangeY); + const cx = (minX + maxX) / 2; + const cy = (minY + maxY) / 2; + const svgW = rangeX * scale + 2 * padding; + const svgH = rangeY * scale + 2 * padding; + + function toSVG(geo: GeoPoint): Point2D { + const mx = (geo.lng - refLng) * cosLat * METERS_PER_DEG; + const my = -(geo.lat - refLat) * METERS_PER_DEG; + return { + x: (mx - cx) * scale + svgW / 2, + y: (my - cy) * scale + svgH / 2, + }; + } + + function toGeo(pt: Point2D): GeoPoint { + const mx = (pt.x - svgW / 2) / scale + cx; + const my = (pt.y - svgH / 2) / scale + cy; + return { + lng: mx / (cosLat * METERS_PER_DEG) + refLng, + lat: -(my / METERS_PER_DEG) + refLat, + }; + } + + return { + toSVG, + toGeo, + viewBox: `0 0 ${Math.round(svgW)} ${Math.round(svgH)}`, + }; +} diff --git a/modules/trips/lib/types.ts b/modules/trips/lib/types.ts new file mode 100644 index 0000000..91c2cc2 --- /dev/null +++ b/modules/trips/lib/types.ts @@ -0,0 +1,64 @@ +/** A point in projected SVG coordinates */ +export interface Point2D { + x: number; + y: number; +} + +/** A geographic coordinate from OSRM */ +export interface GeoPoint { + lng: number; + lat: number; +} + +/** General conic coefficients: Ax² + Bxy + Cy² + Dx + Ey + F = 0 */ +export interface ConicCoeffs { + A: number; + B: number; + C: number; + D: number; + E: number; + F: number; +} + +/** Classified conic type derived from the discriminant B²-4AC */ +export type ConicType = "ellipse" | "parabola" | "hyperbola" | "circle" | "degenerate"; + +/** Rational quadratic Bézier: exact conic arc from p0 to p2 via control p1, weight w */ +export interface RationalQuadBezier { + p0: Point2D; + p1: Point2D; + p2: Point2D; + w: number; +} + +/** A fitted conic arc segment */ +export interface ConicArc { + coeffs: ConicCoeffs; + type: ConicType; + eccentricity: number; + points: Point2D[]; + beziers: RationalQuadBezier[]; +} + +/** An intersection point between two conic arcs */ +export interface IntersectionPoint { + point: Point2D; + arcIndexA: number; + arcIndexB: number; +} + +/** A fully processed route: OSRM geometry → projected → conic-fitted */ +export interface FittedRoute { + rawCoords: GeoPoint[]; + projectedPoints: Point2D[]; + arcs: ConicArc[]; + distance: number; + duration: number; +} + +/** Input for a single route query */ +export interface RouteInput { + start: GeoPoint; + end: GeoPoint; + label: string; +} diff --git a/modules/trips/mod.ts b/modules/trips/mod.ts index 62c0bea..6a854a8 100644 --- a/modules/trips/mod.ts +++ b/modules/trips/mod.ts @@ -14,6 +14,8 @@ import { getModuleInfoList } from "../../shared/module"; import type { RSpaceModule } from "../../shared/module"; import { verifyEncryptIDToken, extractToken } from "@encryptid/sdk/server"; +const OSRM_URL = process.env.OSRM_URL || "http://osrm-backend:5000"; + const routes = new Hono(); // ── DB initialization ── @@ -217,6 +219,37 @@ routes.patch("/api/packing/:id", async (c) => { return c.json(rows[0]); }); +// ── OSRM proxy for route planner ── +routes.post("/api/route", async (c) => { + const body = await c.req.json<{ start: { lng: number; lat: number }; end: { lng: number; lat: number } }>(); + const { start, end } = body; + if (!start || !end) return c.json({ error: "start and end are required" }, 400); + + const url = `${OSRM_URL}/route/v1/driving/${start.lng},${start.lat};${end.lng},${end.lat}?geometries=geojson&overview=full`; + try { + const res = await fetch(url); + const data = await res.json(); + return c.json(data, res.status as any); + } catch (e) { + return c.json({ error: "OSRM backend unavailable", detail: String(e) }, 502); + } +}); + +// ── Route planner page ── +routes.get("/routes", (c) => { + const space = c.req.param("space") || "demo"; + return c.html(renderShell({ + title: `${space} — Route Planner | rTrips`, + moduleId: "trips", + spaceSlug: space, + modules: getModuleInfoList(), + theme: "light", + styles: ``, + body: ``, + scripts: ``, + })); +}); + // ── Page route ── routes.get("/", (c) => { const space = c.req.param("space") || "demo"; diff --git a/server/index.ts b/server/index.ts index c4c8b03..cdd9026 100644 --- a/server/index.ts +++ b/server/index.ts @@ -59,7 +59,6 @@ import { networkModule } from "../modules/network/mod"; import { tubeModule } from "../modules/tube/mod"; import { inboxModule } from "../modules/inbox/mod"; import { dataModule } from "../modules/data/mod"; -import { conicModule } from "../modules/conic/mod"; import { splatModule } from "../modules/splat/mod"; import { spaces } from "./spaces"; import { renderShell } from "./shell"; @@ -86,7 +85,6 @@ registerModule(networkModule); registerModule(tubeModule); registerModule(inboxModule); registerModule(dataModule); -registerModule(conicModule); registerModule(splatModule); // ── Config ── diff --git a/vite.config.ts b/vite.config.ts index 896441c..9db0300 100644 --- a/vite.config.ts +++ b/vite.config.ts @@ -626,38 +626,37 @@ export default defineConfig({ resolve(__dirname, "dist/modules/data/data.css"), ); - // Build conic module component + // Build route planner component (part of trips module) await build({ configFile: false, - root: resolve(__dirname, "modules/conic/components"), + root: resolve(__dirname, "modules/trips/components"), resolve: { alias: { - "../lib/types": resolve(__dirname, "modules/conic/lib/types.ts"), - "../lib/conic-math": resolve(__dirname, "modules/conic/lib/conic-math.ts"), - "../lib/projection": resolve(__dirname, "modules/conic/lib/projection.ts"), + "../lib/types": resolve(__dirname, "modules/trips/lib/types.ts"), + "../lib/conic-math": resolve(__dirname, "modules/trips/lib/conic-math.ts"), + "../lib/projection": resolve(__dirname, "modules/trips/lib/projection.ts"), }, }, build: { emptyOutDir: false, - outDir: resolve(__dirname, "dist/modules/conic"), + outDir: resolve(__dirname, "dist/modules/trips"), lib: { - entry: resolve(__dirname, "modules/conic/components/folk-conic-viewer.ts"), + entry: resolve(__dirname, "modules/trips/components/folk-route-planner.ts"), formats: ["es"], - fileName: () => "folk-conic-viewer.js", + fileName: () => "folk-route-planner.js", }, rollupOptions: { output: { - entryFileNames: "folk-conic-viewer.js", + entryFileNames: "folk-route-planner.js", }, }, }, }); - // Copy conic CSS - mkdirSync(resolve(__dirname, "dist/modules/conic"), { recursive: true }); + // Copy route planner CSS copyFileSync( - resolve(__dirname, "modules/conic/components/conic.css"), - resolve(__dirname, "dist/modules/conic/conic.css"), + resolve(__dirname, "modules/trips/components/route-planner.css"), + resolve(__dirname, "dist/modules/trips/route-planner.css"), ); // Build splat module component