计算多个正方形与圆重叠的面积分数

Question

这是一个基于编程问题的几何问题。我有一个MySQL数据库，里面存储了一些经纬度点，这些点之间距离为1公里，对应着每个点周围1平方公里范围内居住的人口。现在我想知道任意大小的圆形覆盖这些网格时所占比例，以便大致计算出圆形范围内居住的人口数量。

以下是一个实际问题的实例（距离不按比例）：

我想知道在X点半径范围内的人口数量。我的数据库确定点A和B与点X足够接近。在这个例子中，点A的位置大约是40.7458，-74.0375，点B的位置大约是40.7458，-74.0292。从A和B到其网格边缘的每条绿线代表0.5公里，因此A和B周围的灰色圆分别代表1平方公里。

点X位于40.744，-74.032左右，具有0.05公里的半径。现在我可以使用地理三角函数轻松计算出所示的红线。因此，我知道线AX大约为0.504公里，线BX的距离大约为0.309公里。

那么我的问题是：如何计算紫色圆圈内的网格A和网格B所占比例呢？

最终，我将通过这个比例乘以总人口数量。在这种情况下，1平方公里的网格对应着9561人口，而B周围的网格对应着10763人口。因此，如果我知道（只是假设）X周围的半径覆盖了A面积的1％和B面积的3％，我可以通过将A和B的人口分别乘以它们各自的比例并相加来合理地估算圆形范围内的总人口数量。

对于上面的两个正方形，我已经计算出了比例，但是根据半径的大小可能会有许多可能的正方形，如下图所示，这使得它成为一个更一般的问题：

在某些情况下，如果很容易确定正方形网格被半径完全包含，那么这个问题原则上就很简单（例如，如果AX之间的距离小于X周围的半径，那么我就知道不需要进行进一步的数学计算）。

现在，确定哪些点在圆圈范围内很容易，但是我无法确定它们各自所占比例。谢谢您的帮助。

- nucleon

---

在计算机图形学中，当人们想要绘制抗锯齿的东西时，他们使用一个称为 alpha 的值来描述形状覆盖了多少像素。并非所有人都认为像素是简单的正方形，但许多人确实这样认为。因此，寻找计算围绕边界覆盖范围的圆形扫描转换算法可能会有所帮助。无论如何，与最近中心点的水平和垂直距离的分离将比直线距离更有用，因为前者将该点与分隔每个网格点关联区域的边缘相关联。 - MvG

谢谢 - 这是一个我之前没有想到的非常有趣的方法。谷歌搜索“圆形扫描转换算法”会出现很多有希望的网站。很高兴知道这不是一个简单的问题，而不仅仅是我太蠢了。 - nucleon

Answer 1

我最终得出了一个相当不错的近似解决方案。以下是 PHP 代码：

//$p is an array of latitude, longitude, value, and distance from the centerpoint 
//$cx,$cy are the lat/lon of the center point, $cr is the radius of the circle
//$pdist is the distance from each node to its edge (in this case, .5 km, since it is a 1km x 1km grid)

function sum_circle($p, $cx, $cy, $cr, $pdist) {
    $total = 0; //initialize the total
    $hyp = sqrt(($pdist*$pdist)+($pdist*$pdist)); //hypotenuse of distance

    for($i=0; $i<count($p); $i++) { //cycle over all points
        $px = $p[$i][0];    //x value of point
        $py = $p[$i][1];    //y value of point
        $pv = $p[$i][2];    //associated value of point (e.g. population)
        $dist = $p[$i][3];  //calculated distance of point coordinate to centerpoint

        //first, the easy case — items that are well outside the maximum distance
        if($dist>$cr+$hyp) { //if the distance is greater than circle radius plus the hypoteneuse
            $per = 0; //then use 0% of its associated value
        } else if($dist+$hyp<=$cr) { //other easy case - completely inside circle (distance + hypotenuse <= radius)
            $per = 1; //then use 100% of its associated value
        } else { //the edge cases
            $mx = ($cx-$px); $my = ($cy-$py); //calculate the angle of the difference
            $theta = abs(rad2deg(atan2($my,$mx)));
            $theta = abs((($theta + 89) % 90 + 90) % 90 - 89); //reduce it to a positive degree between 0 and 90
            $tf = abs(1-($theta/45)); //this basically makes it so that if the angle is close to 45, it returns 0, 
                                      //if it is close to 0 or 90, it returns 1
            $hyp_adjust = ($hyp*(1-$tf)+($pdist*$tf)); //now we create a mixed value that is weighted by whether the
                                                        //hypotenuse or the distance between cells should be used                                                   

            $per = ($cr-$dist+$hyp_adjust)/100; //lastly, we use the above numbers to estimate what percentage of 
                                                //the square associated with the centerpoint is covered
            if($per>1) $per = 1; //normalize for over 100% or under 0%
            if($per<0) $per = 0; 
        }
        $total+=$per*$pv;   //add the value multiplied by the percentage to the total
    }  
    return $total;
}

这似乎很有效并且速度很快（即使在边缘情况下它确实使用了一些三角函数）。基本逻辑是，在计算边缘情况时，圆形半径的两个极端可能性是，圆形半径要么正好垂直于网格，要么正好与其呈45度角。因此，它大致确定了它在这些极端情况之间的位置，然后使用它来大致计算覆盖网格正方形的百分比。在我的测试中，它给出了合理的结果。

对于我使用的正方形和圆形的大小，这似乎是足够的？

我写了一个小应用程序来尝试帮助我解决这个问题。不解释全部内容，您可以通过查看此屏幕截图来了解算法的“思考”方式： 基本上，如果圆形是黄色的，那就意味着它已经确定它是100％的，如果是红色的，则已经快速地被筛选为100％的。其他情况是边缘情况。点下面的数字（从0到1）是使用上述方法计算得出的（四舍五入的）覆盖百分比，而下面的数字是上述代码中使用的计算θ值。

对于我的目的，我认为这个近似值是可行的。

- nucleon

Answer 2

通过足够的分类（如下所示），所有计算都可以简化为一个“原始”计算，即提供图像中橙色区域的“角面积”的计算。

当y0>0，如上图所示，并且无论x0是正数还是负数，橙色区域可以精确计算为从x0到x1的积分sqrt(r^2-y^2)减去矩形面积(x1-x0)*(y1-y0)。该积分具有众所周知的闭合表达式，因此不需要使用任何数值算法进行计算。

圆与正方形之间的其他相交部分可以简化为组合矩形和直角形，如上图中涂成橙色的部分。例如，由以下图片中的水平和垂直橙色射线限定的相交部分可以通过求红色矩形加上两个角形（蓝色和绿色）的面积来表示。

蓝色区域是通过直接应用上述原始情况得出的（其中下部矩形折叠为零）。同样地，绿色区域也可以用相同的方法进行测量，只需将负的y坐标替换为它的绝对值（另一个y为0）。

应用这些思想，可以枚举所有情况。基本上，应该考虑正方形中有一个、两个、三个或四个角落在圆形内，而其余的（如果有的话）则落在外面的情况。枚举本身就是一个问题，但至少在理论上，可以通过考虑相对较少的“典型”配置来解决。

对于每个按描述枚举的情况，都必须计算一些矩形和角度区域的分解，并将部分相加（或相减），如上述三色示例所示。每个部分的面积将缩小到矩形或原始角度区域。

要将这种攻击方式转化为工作算法需要做相当多的工作。更深入的分析可以揭示如何最小化要考虑的“典型”配置数量。如果不能，我认为要考虑的组合数量，无论多大，都应该是可管理的。

- Leandro Caniglia

Answer 3

如果您的问题可以得到近似答案，那么还有另一种更简单的编程技巧可供使用。这个问题的整个思路就是计算正方形和圆形的交集面积。我在我的其他回答中没有解释这一点，但是找到可能截取圆形的正方形不应该是一个问题，否则，请告诉我们。

计算交集的近似面积的想法非常简单。在正方形中随机生成足够多的点，并检查其中有多少属于圆形。圆形中的点数与正方形中随机点的总数之比将给出相对于正方形面积的交集比例。

现在，考虑到您必须为所有围绕圆形的正方形（即，中心距离圆心不太不同于圆半径的正方形）重复相同的例程，因此您可以通过将它们从一个正方形转换到另一个正方形来重复使用随机点。

如果这种方法不适合您的问题，我就不想深入讨论了，但是让我简单指出，在正方形中生成均匀分布的随机点相当容易。您只需要为x坐标和独立的y坐标生成随机数即可。然后考虑所有的(x, y)对。然后，对于每个(x, y)，验证是否(x - a)^2 + (y - b)^2 <= r^2，其中(a, b)代表圆心，r代表半径。

- Leandro Caniglia