diff --git a/packages/math/src/geom.ts b/packages/math/src/geom.ts index c24a62f..fa9d028 100644 --- a/packages/math/src/geom.ts +++ b/packages/math/src/geom.ts @@ -4,7 +4,9 @@ */ import { polygonHull as d3_polygonHull, polygonCentroid as d3_polygonCentroid } from 'd3-polygon'; +import { ANGLE_EPSILON, HALF_PI } from './constants'; import { Extent } from './Extent'; +import { numWrap } from './number'; import { Vec2, vecLength } from './vector'; @@ -258,69 +260,130 @@ export function geomPolygonIntersectsPolygon(outer: Vec2[], inner: Vec2[], check } -/** Smallest Surrounding Rectangle. An Object containing `poly` and `angle` properties. */ -export interface SSR { - /** the smallest surrounding rectangle polygon */ +/** Surrounding Rectangle. An Object containing `poly` and `angle` properties. */ +export interface SurroundingRectangle { + /** the surrounding rectangle polygon */ poly: Vec2[]; /** angle offset from x axis */ angle: number; } -/** Return the Smallest Surrounding Rectangle for a given set of points - * @remarks - * http://gis.stackexchange.com/questions/22895/finding-minimum-area-rectangle-for-given-points - * http://gis.stackexchange.com/questions/3739/generalisation-strategies-for-building-outlines/3756#3756 +/** + * A generalized function for computing a surrounding rectangle for a given Array of points. + * The caller supplies `getScore` and `isBetter` functions that define the heuristic which + * determines what kind of rectangle to gather ("smallest" or "longest"). * @param points - * @returns rectangle if exists, null otherwise - * @example - * +-- p1 ------ p3 - * | | - * p0 ------ p2 --+ - * const points = [[0, -1], [5, 1], [10, -1], [15, 1]]; - * const ssr = geomGetSmallestSurroundingRectangle(points); - * // ssr.poly == [[0, -1], [0, 1], [15, 1], [15, -1], [0, -1]] - * // ssr.angle == 0 + * @param getScore + * @param isBetter + * @param initialBestScore + * @returns A SurroundingRectangle (polygon and angle) or `null` if the given points did not produce a valid hull. */ -export function geomGetSmallestSurroundingRectangle(points: Vec2[]): SSR | null { +export function getSurroundingRectangle( + points: Vec2[], + getScore: (extent: Extent) => number, + isBetter: (score: number, bestScore: number) => boolean, + initialBestScore: number +): SurroundingRectangle | null { const hull: Vec2[] | null = d3_polygonHull(points); if (!hull) return null; const centroid: Vec2 = d3_polygonCentroid(hull); - let minArea = Infinity; - let ssrExtent = new Extent(); - let ssrAngle = 0; - let c1: Vec2 = hull[0]; - - for (let i = 0; i <= hull.length - 1; i++) { + const centroidX = centroid[0]; + const centroidY = centroid[1]; + const checkedAngles = new Set(); + let bestScore = initialBestScore; + let bestExtent = new Extent(); + let bestAngle = 0; + + for (let i = 0; i < hull.length; i++) { + const c1: Vec2 = hull[i]; const c2: Vec2 = i === hull.length - 1 ? hull[0] : hull[i + 1]; - const angle = Math.atan2(c2[1] - c1[1], c2[0] - c1[0]); - const poly: Vec2[] = geomRotatePoints(hull, -angle, centroid); - const extent: Extent = poly.reduce((acc: Extent, point: Vec2) => { + const wrapped = numWrap(Math.atan2(c2[1] - c1[1], c2[0] - c1[0]), 0, HALF_PI); // normalize angle between 0..π/2 + const angle = (wrapped < ANGLE_EPSILON || (HALF_PI - wrapped) < ANGLE_EPSILON) ? 0 : wrapped; + const angleKey = Math.round(angle / ANGLE_EPSILON); + if (checkedAngles.has(angleKey)) continue; + checkedAngles.add(angleKey); + + const sin = Math.sin(-angle); + const cos = Math.cos(-angle); + const extent = new Extent(); + for (const point of hull) { + const radialX = point[0] - centroidX; + const radialY = point[1] - centroidY; + const x = radialX * cos - radialY * sin + centroidX; + const y = radialX * sin + radialY * cos + centroidY; + // update Extent min/max in-place for speed - acc.min[0] = Math.min(acc.min[0], point[0]); - acc.min[1] = Math.min(acc.min[1], point[1]); - acc.max[0] = Math.max(acc.max[0], point[0]); - acc.max[1] = Math.max(acc.max[1], point[1]); - return acc; - }, new Extent()); - - const area = extent.area(); - if (area < minArea) { - minArea = area; - ssrExtent = extent; - ssrAngle = angle; + extent.min[0] = Math.min(extent.min[0], x); + extent.min[1] = Math.min(extent.min[1], y); + extent.max[0] = Math.max(extent.max[0], x); + extent.max[1] = Math.max(extent.max[1], y); + } + + const score = getScore(extent); + if (isBetter(score, bestScore)) { + bestScore = score; + bestExtent = extent; + bestAngle = angle; } - c1 = c2; } return { - poly: geomRotatePoints(ssrExtent.polygon(), ssrAngle, centroid), - angle: ssrAngle + poly: geomRotatePoints(bestExtent.polygon(), bestAngle, centroid), + angle: bestAngle }; } +/** Return the Smallest Surrounding Rectangle for a given Array of points + * @remarks + * http://gis.stackexchange.com/questions/22895/finding-minimum-area-rectangle-for-given-points + * http://gis.stackexchange.com/questions/3739/generalisation-strategies-for-building-outlines/3756#3756 + * @param points + * @returns The smallest `SurroundingRectangle` by area, or `null` if the given points did not produce a valid hull. + * @example + * p5 --- p4 + * | | + * | p3 ------ p2 + * | | + * p0 ------------- p1 + * const footprint = [[0, 0], [8, 0], [8, 2], [3, 2], [3, 6], [0, 6], [0, 0]]; + * const points = geomRotatePoints(footprint, Math.PI / 6, [0, 0]); + * const ssr = geomGetSmallestSurroundingRectangle(points); + * // ssr.angle == Math.PI / 6 + */ +export function geomGetSmallestSurroundingRectangle(points: Vec2[]): SurroundingRectangle | null { + const getScore = (extent: Extent) => extent.area(); + const isBest = (score: number, best: number) => score < best; + return getSurroundingRectangle(points, getScore, isBest, Infinity); +} + + +/** Return the Longest Surrounding Rectangle for a given Array of points + * @remarks + * Loops over the convex hull edges, rotates to each edge angle, and chooses the + * rectangle with the maximum side length. + * @param points + * @returns The longest `SurroundingRectangle` by side length, or `null` if the given points did not produce a valid hull. + * @example + * p5 --- p4 + * | | + * | p3 ------ p2 + * | | + * p0 ------------- p1 + * const footprint = [[0, 0], [8, 0], [8, 2], [3, 2], [3, 6], [0, 6], [0, 0]]; + * const points = geomRotatePoints(footprint, Math.PI / 6, [0, 0]); + * const lsr = geomGetLongestSurroundingRectangle(points); + * // lsr.angle ~= 1.4196541601696426 + */ +export function geomGetLongestSurroundingRectangle(points: Vec2[]): SurroundingRectangle | null { + const getScore = (extent: Extent) => Math.max(Math.abs(extent.max[0] - extent.min[0]), Math.abs(extent.max[1] - extent.min[1])); + const isBest = (score: number, best: number) => score > best; + return getSurroundingRectangle(points, getScore, isBest, -Infinity); +} + + /** Return the length of the given path * @param path * @returns length diff --git a/packages/math/test/geom.test.js b/packages/math/test/geom.test.js index 5cd2eaa..f4f6559 100644 --- a/packages/math/test/geom.test.js +++ b/packages/math/test/geom.test.js @@ -9,6 +9,39 @@ assert.closeTo = function(a, b, epsilon = 1e-9) { } }; +function rectangleSideLengths(rectangle) { + const poly = rectangle.poly; + const sideA = Math.hypot(poly[1][0] - poly[0][0], poly[1][1] - poly[0][1]); + const sideB = Math.hypot(poly[2][0] - poly[1][0], poly[2][1] - poly[1][1]); + return sideA > sideB ? [sideA, sideB] : [sideB, sideA]; +} + +function rectangleArea(rectangle) { + const [longSide, shortSide] = rectangleSideLengths(rectangle); + return longSide * shortSide; +} + +function assertRectangleContainsPoints(rectangle, points) { + const origin = rectangle.poly[0]; + const edgeA = [rectangle.poly[1][0] - origin[0], rectangle.poly[1][1] - origin[1]]; + const edgeB = [rectangle.poly[3][0] - origin[0], rectangle.poly[3][1] - origin[1]]; + const edgeALengthSquare = edgeA[0] * edgeA[0] + edgeA[1] * edgeA[1]; + const edgeBLengthSquare = edgeB[0] * edgeB[0] + edgeB[1] * edgeB[1]; + + for (const point of points) { + const pointVector = [point[0] - origin[0], point[1] - origin[1]]; + const alongA = (pointVector[0] * edgeA[0] + pointVector[1] * edgeA[1]) / edgeALengthSquare; + const alongB = (pointVector[0] * edgeB[0] + pointVector[1] * edgeB[1]) / edgeBLengthSquare; + assert.ok(alongA >= -1e-9 && alongA <= 1 + 1e-9); + assert.ok(alongB >= -1e-9 && alongB <= 1 + 1e-9); + } +} + +function offAxisLShapedBuildingPoints() { + const footprint = [[0, 0], [8, 0], [8, 2], [3, 2], [3, 6], [0, 6], [0, 0]]; + return math.geomRotatePoints(footprint, Math.PI / 6, [0, 0]); +} + describe('math/geom', () => { describe('geomEdgeEqual', () => { it('returns false for inequal edges', () => { @@ -274,16 +307,64 @@ describe('math/geom', () => { assert.equal(math.geomGetSmallestSurroundingRectangle([]), null); }); - it('calculates a smallest surrounding rectangle', () => { - // +-- p1 ------ p3 - // | | - // p0 ------ p2 --+ - const points = [[0, -1], [5, 1], [10, -1], [15, 1]]; + it('normalizes nearly equivalent axis angles', () => { + const points = [[0, 0], [10, -1e-13], [10, 2], [0, 2], [0, 0]]; const ssr = math.geomGetSmallestSurroundingRectangle(points); assert.ok(ssr instanceof Object); - assert.deepEqual(ssr.poly, [[0, -1], [15, -1], [15, 1], [0, 1], [0, -1]]); assert.equal(ssr.angle, 0); }); + + it('calculates a smallest surrounding rectangle for an off-axis L-shaped building', () => { + // p5 --- p4 + // | | + // | p3 ------ p2 + // | | + // p0 ------------- p1 + const points = offAxisLShapedBuildingPoints(); + const ssr = math.geomGetSmallestSurroundingRectangle(points); + assert.ok(ssr instanceof Object); + assertRectangleContainsPoints(ssr, points); + + const [longSide, shortSide] = rectangleSideLengths(ssr); + assert.closeTo(longSide, 8); + assert.closeTo(shortSide, 6); + assert.closeTo(ssr.angle, Math.PI / 6); + }); + }); + + describe('geomGetLongestSurroundingRectangle', () => { + it('returns null for empty points array', () => { + assert.equal(math.geomGetLongestSurroundingRectangle([]), null); + }); + + it('calculates a longest surrounding rectangle for an off-axis L-shaped building', () => { + // p5 --- p4 + // | | + // | p3 ------ p2 + // | | + // p0 ------------- p1 + const points = offAxisLShapedBuildingPoints(); + const lsr = math.geomGetLongestSurroundingRectangle(points); + assert.ok(lsr instanceof Object); + assertRectangleContainsPoints(lsr, points); + + const [longSide, shortSide] = rectangleSideLengths(lsr); + assert.closeTo(longSide, 9.995120760870789); + assert.closeTo(shortSide, 6.559297999321455); + assert.closeTo(lsr.angle, 1.4196541601696426); + }); + + it('can choose a longer envelope than the smallest surrounding rectangle', () => { + const points = offAxisLShapedBuildingPoints(); + const smallest = math.geomGetSmallestSurroundingRectangle(points); + const longest = math.geomGetLongestSurroundingRectangle(points); + + assert.ok(smallest instanceof Object); + assert.ok(longest instanceof Object); + assert.ok(Math.abs(smallest.angle - longest.angle) > 1e-9); + assert.ok(rectangleSideLengths(longest)[0] > rectangleSideLengths(smallest)[0]); + assert.ok(rectangleArea(longest) > rectangleArea(smallest)); + }); }); describe('geomPathLength', () => {