最小矩形覆盖
题面
给定一些点的坐标,要求能够覆盖所有点的最小面积的矩形,输出所求矩形的面积和四个顶点的坐标
分析
首先可以先求凸包,因为覆盖了凸包上的顶点,凸包内的顶点也一定能被覆盖
结论:这个矩形的一条边一定与凸包的一条边重合。
然后对于凸包的每一条边,我们通过旋转卡壳找到最左侧的点l,最右侧的点r,最高点p,通过这在矩形四条边上的点就可以求得矩形的四个顶点,就得到了我们要求的矩形。
什么是旋转卡壳
旋转卡壳(Rotating Calipers,也称「旋转卡尺」)算法,在凸包算法的基础上,通过枚举凸包上某一条边的同时维护其他需要的点(这里就是需要得到上面提到的最左侧的点l,最右侧的点r和最上面的点u),能够在线性时间内求解如凸包直径、最小矩形覆盖等和凸包性质相关的问题。
代码实现
const length = points.length;
if (length < 3) return [];
function next(i: number) {
return (i + 1) % length;
}
const rects: Point[][] = [];
let smallestRect: Point[] = [];
for (let i = 0; i < length; i++) {
let l = next(i),
r = next(i),
up = next(i);
const p1 = points[i];
const p2 = points[next(i)];
while (
Vector.crossProduct2D(
GeoPoint.offset(p1, p2),
GeoPoint.offset(p1, points[up])
) -
Vector.crossProduct2D(
GeoPoint.offset(p1, p2),
GeoPoint.offset(p1, points[next(up)])
) <
0
)
up = next(up);
while (
Vector.dotProduct2D(
GeoPoint.offset(p1, p2),
GeoPoint.offset(p1, points[r])
) -
Vector.dotProduct2D(
GeoPoint.offset(p1, p2),
GeoPoint.offset(p1, points[next(r)])
) <
0
)
r = next(r);
if (!i) l = r;
while (
Vector.dotProduct2D(
GeoPoint.offset(p1, p2),
GeoPoint.offset(p1, points[1])
) -
Vector.dotProduct2D(
GeoPoint.offset(p1, p2),
GeoPoint.offset(p1, points[next(1)])
) >
0
)
l = next(1);
const r1 = new GeoLine(p1, p2).ptProjection(points[r]);
const r2 = new GeoLine(p1, p2).ptProjection(points[l]);
const k = GeoLine.getK(p1, p2);
const line = {
a: GeoLine.getK(p1, p2),
b: -1,
c: points[up].y - k * points[up].x,
};
let r3: Point;
let r4: Point;
if (k == Infinity || k == -Infinity) {
const distance = points[up].x - p1.x;
r3 = {
x: r1.x + distance,
y: r1.y,
};
r4 = {
x: r2.x + distance,
y: r2.y,
};
} else {
r3 = GeoLine.normalLineFeet(line, r1);
r4 = GeoLine.normalLineFeet(line, r2);
}
rects.push([r1, r2, r3, r4]);
}
let area = Infinity;
for (const rect of rects) {
const currentArea =
GeoPoint.distance(rect[0], rect[1]) * GeoPoint.distance(rect[1], rect[2]);
if (currentArea < area) {
smallestRect = rect;
area = currentArea;
}
}
return smallestRect;
}
工具类
type Line2 = {
p1: Point;
p2: Point;
};
type Point = {
x: number;
y: number;
};
type Line3 = {
a: number;
b: number;
c: number;
};
const DEFAULT_TOL = 0.0000001;
class GeoPoint {
static offset(p1: Point, p2: Point): Point {
return { x: p2.x - p1.x, y: p2.y - p1.y };
}
static distance(p1: Point, p2: Point): number {
if (!p1 || !p2) return 0;
const x = p1.x - p2.x;
const y = p1.y - p2.y;
return Math.sqrt(x * x + y * y);
}
static getLeftBottomPoint(points: Point[]): Point {
let pole = points[0];
for (const point of points) {
if (point.y < pole.y || (point.y == pole.y && point.x < pole.x)) {
pole = point;
continue;
}
}
return pole;
}
static defaultPoint(): Point {
return {
x: 0,
y: 0,
};
}
static pointEqual(p1: Point, p2: Point) {
return p1.x == p2.x && p1.y == p2.y;
}
}
class Vector {
static crossProduct2D(p1: Point, p2: Point) {
return p1.x * p2.y - p1.y * p2.x;
}
static dotProduct2D(p1: Point, p2: Point) {
return p1.x * p2.x + p1.y * p2.y;
}
}
class GeoLine {
public p1: Point;
public p2: Point;
constructor(p1: Point, p2: Point) {
this.p1 = p1;
this.p2 = p2;
}
public ptProjection(pt: Point, tol: number = DEFAULT_TOL): Point {
const A: Point = this.p1,
B: Point = this.p2;
const v: Point = {
x: B.x - A.x,
y: B.y - A.y,
};
const t =
(v.x * (pt.x - A.x) + v.y * (pt.y - A.y)) / (v.x * v.x + v.y * v.y);
return {
x: A.x + v.x * t,
y: A.y + v.y * t,
};
}
static getK(P1: Point, p2: Point) {
return (p2.y - P1.y) / (p2.x - P1.x);
}
static normalLineFeet(norLine: Line3, pos: Point) {
const { a, b, c } = norLine;
const { x, y } = pos;
return {
x: (b * b * x - a * b * y - a * c) / (a * a + b * b),
y: (-a * b * x + a * a * y - b * c) / (a * +b * b),
};
}
}
参考文献