图形学之最小矩形覆盖(JS实现)

213 阅读3分钟

最小矩形覆盖

题面

给定一些点的坐标,要求能够覆盖所有点的最小面积的矩形,输出所求矩形的面积和四个顶点的坐标

分析

首先可以先求凸包,因为覆盖了凸包上的顶点,凸包内的顶点也一定能被覆盖

结论:这个矩形的一条边一定与凸包的一条边重合。

然后对于凸包的每一条边,我们通过旋转卡壳找到最左侧的点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),
    };
  }
}

参考文献

旋转卡壳

P3187 [HNOI2007] 最小矩形覆盖