跳转至

计算几何——固定半径圆能覆盖平面上的最大点数

给定一个半径为r的圆,给定一系列平面上的点,并给出横纵坐标,问圆最多能覆盖多少个点。点在圆上或者在圆内都被视为覆盖。

  • POJ 1981 Circle and Points
  • LeetCode 1453.圆形靶内的最大飞镖数量

暴力枚举

一种思路是从平面点里面任意选出两个点,以这两个点的距离为弦,然后判断其他点是否在圆内或圆上。时间复杂度O(n^3)

平面扫描

以每一个点做为圆心画一个圆,枚举与其相交的其他圆,保存交点和角度,按角度排序,再次扫描一遍。时间复杂度O(n^2 \log n)

假设我们选取圆A作为分析的起点,以半径r做圆,圆A和圆E相交于点K和点J,如果圆心放在弧JK上,那么必然保证点A和点E可以被圆覆盖。所以问题转化成,求与圆A相交的圆中,圆A的某一段弧被最多的圆包含,上图中,弧OJ就被圆H和圆E包含,那么圆心放在弧OJ上,必然保证点A,点E和点H被圆覆盖。

接下来是如何去描述这个弧,过点J做AE的中垂线,交点为P,显然PE = \frac{AE}{2},其中AS所在的直线为水平线,我们设\angle AEJ = \phi, \angle SAE = \theta,我们把极坐标的原点建立在点A,因为弧在圆A上,半径是固定的,那么只需要弧对应的圆心角的起始角度和终止角度即可唯一确定一段弧。

比如描述弧JK,我们只需要求出\angle KAS, \angle SAJ即可,显然\angle KAS = \theta - \phi, \angle SAJ = \theta + \phi,那么根据几何关系求出两个角度就不难了。

现在考虑当E在过A的水平线下方的情况:

发现起始角度是\theta- \phi,终止角度是\phi + \theta

另外两种情况只需要转换参考系,起始角度和终止角度的表达形式都不变。在计算\theta的时候用到了函数

double atan2(double y, double x);
\operatorname{atan} 2(y, x)=\left\{\begin{array}{ll} \arctan \left(\frac{y}{x}\right) & x>0 \\ \arctan \left(\frac{y}{x}\right)+\pi & y \geq 0, x<0 \\ \arctan \left(\frac{y}{x}\right)-\pi & y<0, x<0 \\ +\frac{\pi}{2} & y>0, x=0 \\ -\frac{\pi}{2} & y<0, x=0 \\ \text { undefined } & y=0, x=0 \end{array}\right.

现在问题是如何找出哪一段弧被最多的圆包含,这就变成了一个扫描线问题。

一个小细节就是在排序的时候需要考虑当角度相同的情况,第一个样例就是一个很好的例子,假设考虑第一个点(-2,0),那么与其相交的三段弧的范围是:[\frac{\pi}{2}, 0], [0,0], [0, \frac{\pi}{2}],很容易看出原点(0,0)就是结果。在排序的过程中,会有4个0存在,但是从实际角度来看,在角度相同的情况,最先应该扫描起点,所以在排序的过程中增加了一个判断。

class Solution {
    struct Node
    {
        double angle;
        bool isStartPoint;

        bool operator<(const Node & obj) const
        {
            return (angle < obj.angle) || (angle == obj.angle && isStartPoint);
        }
    };

public:
    int numPoints(vector<vector<int>>& points, int r) {
        std::ios_base::sync_with_stdio(false);
        cin.tie(NULL);
        cout.tie(NULL);

        int n = points.size();
        int res = 1;
        vector<Node> seq(n * n);

        for (int i = 0; i < n; ++i) {
            int m = 0;
            for (int j = 0; j < n; ++j) {
                double d = dis(points[i], points[j]);
                if (i != j && d <= 2.0 * r) {
                    double phi = acos(d / 2 / r);
                    double theta = atan2(points[j][1] * 1.0 - points[i][1], points[j][0] * 1.0 - points[i][0]);
                    seq[m].angle = theta - phi; seq[m++].isStartPoint = true;
                    seq[m].angle = theta + phi; seq[m++].isStartPoint = false;
                }
            }

            sort(seq.begin(), seq.begin() + m);

            int cnt = 1;
            for (int j = 0; j < m; ++j) {
                if (seq[j].isStartPoint) ++cnt;
                else --cnt;
                res = max(res, cnt);
            }
        }

        return res;
    }

    inline double dis(vector<int> & p1, vector<int> & p2)
    {
        return sqrt((p1[0] - p2[0]) * (p1[0] - p2[0]) * 1.0 
            + (p1[1] - p2[1]) * (p1[1] - p2[1]) * 1.0);
    }
};