diff --git a/src/distanceField/cpt.ts b/src/distanceField/cpt.ts index 0c52980..627701c 100644 --- a/src/distanceField/cpt.ts +++ b/src/distanceField/cpt.ts @@ -1,6 +1,6 @@ import type { Vector2 } from '../utils/Vector2.ts'; -import { findHullParabolas, transpose } from './utils.ts'; +/** Adapted from Felzenszwalb, P. F., & Huttenlocher, D. P. (2012). Distance Transforms of Sampled Functions. Theory of Computing, 8(1), 415–428. */ export function computeCPT( sedt: Float32Array[], cpt: Vector2[][], @@ -8,40 +8,62 @@ export function computeCPT( ycoords: Float32Array[] ): Vector2[][] { const length = sedt.length; + const tempArray = new Float32Array(length); + + // Pre-allocate hull arrays + const hullVertices: Vector2[] = []; + const hullIntersections: Vector2[] = []; for (let row = 0; row < length; row++) { - horizontalPass(sedt[row], xcoords[row]); + horizontalPass(sedt[row], xcoords[row], hullVertices, hullIntersections); } - transpose(sedt); + for (let i = 0; i < length; i++) { + for (let j = i + 1; j < length; j++) { + tempArray[0] = sedt[i][j]; + sedt[i][j] = sedt[j][i]; + sedt[j][i] = tempArray[0]; + } + } for (let row = 0; row < length; row++) { - horizontalPass(sedt[row], ycoords[row]); + horizontalPass(sedt[row], ycoords[row], hullVertices, hullIntersections); } - for (let col = 0; col < length; col++) { - for (let row = 0; row < length; row++) { - const y = ycoords[col][row]; - const x = xcoords[y][col]; - cpt[row][col] = { x, y }; - } + const len = length * length; + for (let i = 0; i < len; i++) { + const row = i % length; + const col = (i / length) | 0; + const y = ycoords[col][row]; + const x = xcoords[y][col]; + cpt[row][col].x = x; + cpt[row][col].y = y; } return cpt; } -function horizontalPass(singleRow: Float32Array, indices: Float32Array) { - const hullVertices: Vector2[] = []; - const hullIntersections: Vector2[] = []; +function horizontalPass( + singleRow: Float32Array, + indices: Float32Array, + hullVertices: Vector2[], + hullIntersections: Vector2[] +) { + // Clear hull arrays before use + hullVertices.length = 0; + hullIntersections.length = 0; + findHullParabolas(singleRow, hullVertices, hullIntersections); marchParabolas(singleRow, hullVertices, hullIntersections, indices); } function marchParabolas(row: Float32Array, verts: Vector2[], intersections: Vector2[], indices: Float32Array) { let k = 0; + const n = row.length; + const numVerts = verts.length; - for (let i = 0; i < row.length; i++) { - while (intersections[k + 1].x < i) { + for (let i = 0; i < n; i++) { + while (k < numVerts - 1 && intersections[k + 1].x < i) { k++; } const dx = i - verts[k].x; @@ -49,3 +71,38 @@ function marchParabolas(row: Float32Array, verts: Vector2[], intersections: Vect indices[i] = verts[k].x; } } + +function findHullParabolas(row: Float32Array, verts: Vector2[], intersections: Vector2[]) { + let k = 0; + + verts[k] = { x: 0, y: row[0] }; + intersections[k] = { x: -Infinity, y: 0 }; + intersections[k + 1] = { x: Infinity, y: 0 }; + + const n = row.length; + + for (let i = 1; i < n; i++) { + const s: Vector2 = { x: 0, y: 0 }; + const qx = i; + const qy = row[i]; + let p = verts[k]; + + // Calculate intersection + s.x = (qy + qx * qx - (p.y + p.x * p.x)) / (2 * (qx - p.x)); + + while (k > 0 && s.x <= intersections[k].x) { + k--; + p = verts[k]; + s.x = (qy + qx * qx - (p.y + p.x * p.x)) / (2 * (qx - p.x)); + } + + k++; + verts[k] = { x: qx, y: qy }; + intersections[k] = { x: s.x, y: 0 }; + intersections[k + 1] = { x: Infinity, y: 0 }; + } + + // Adjust the length of verts and intersections arrays + verts.length = k + 1; + intersections.length = k + 2; +} diff --git a/src/distanceField/fields.ts b/src/distanceField/fields.ts index 827de62..536f247 100644 --- a/src/distanceField/fields.ts +++ b/src/distanceField/fields.ts @@ -147,47 +147,6 @@ export class Fields { } } - private renderer( - pixelRenderer: ( - distance: number, - closestX: number, - closestY: number, - shapeColor: number, - row: number, - col: number - ) => { r: number; g: number; b: number } - ): ImageData { - const imageData = new ImageData(this.resolution, this.resolution); - const data = imageData.data; - const resolution = this.resolution; - - // Pre-cache arrays to avoid repeated property access - const edt = this.edt; - const cpt = this.cpt; - const colorField = this.colorField; - - // Process pixels in a single loop - for (let i = 0; i < resolution * resolution; i++) { - const row = i % resolution; - const col = (i / resolution) | 0; // Faster integer division - const index = i * 4; - - const distance = edt[row][col]; - const { x: closestX, y: closestY } = cpt[row][col]; - const shapeColor = colorField[closestX][closestY]; - - const color = pixelRenderer(distance, closestX, closestY, shapeColor, row, col); - - // Direct array access is faster than property access - data[index] = color.r; - data[index + 1] = color.g; - data[index + 2] = color.b; - data[index + 3] = 255; - } - - return imageData; - } - public generateImageData(): ImageData { const imageData = new ImageData(this.resolution, this.resolution); const data = imageData.data; diff --git a/src/distanceField/utils.ts b/src/distanceField/utils.ts deleted file mode 100644 index aca089f..0000000 --- a/src/distanceField/utils.ts +++ /dev/null @@ -1,41 +0,0 @@ -import type { Vector2 } from '../utils/Vector2.ts'; - -export function intersectParabolas(p: Vector2, q: Vector2): Vector2 { - const x = (q.y + q.x * q.x - (p.y + p.x * p.x)) / (2 * q.x - 2 * p.x); - return { x, y: 0 }; -} - -export function transpose(matrix: Float32Array[]) { - for (let i = 0; i < matrix.length; i++) { - for (let j = i + 1; j < matrix[i].length; j++) { - const temp = matrix[i][j]; - matrix[i][j] = matrix[j][i]; - matrix[j][i] = temp; - } - } -} - -export function findHullParabolas(row: Float32Array, verts: Vector2[], intersections: Vector2[]) { - let k = 0; - - verts[0] = { x: 0, y: row[0] }; - intersections[0] = { x: -Infinity, y: 0 }; - intersections[1] = { x: Infinity, y: 0 }; - - for (let i = 1; i < row.length; i++) { - const q: Vector2 = { x: i, y: row[i] }; - let p = verts[k]; - let s = intersectParabolas(p, q); - - while (s.x <= intersections[k].x) { - k--; - p = verts[k]; - s = intersectParabolas(p, q); - } - - k++; - verts[k] = q; - intersections[k] = s; - intersections[k + 1] = { x: Infinity, y: 0 }; - } -}