Skip to content

Add geomGetDominantSurroundingRectangle and rename SSR type #308

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
bhousel merged 1 commit into from
May 5, 2026
Merged

Add geomGetDominantSurroundingRectangle and rename SSR type #308

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
145 changes: 104 additions & 41 deletions packages/math/src/geom.ts
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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';


Expand Down Expand Up @@ -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<number>();
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
Expand Down
93 changes: 87 additions & 6 deletions packages/math/test/geom.test.js
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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', () => {
Expand Down Expand Up @@ -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', () => {
Expand Down