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

Calculating the fraction of the area of multiple squares overlapped by a circle

这是一道基于我遇到的编程题的几何题。基本上,我有一个充满纬度和经度点的 MySQL 数据库,彼此间隔 1 公里,对应于居住在每个点周围平方公里内的人口。然后我想知道每个网格被重叠的任意大小的圆所占的相对比例,这样我就可以算出有多少人大致生活在给定的圆内。

这是问题的一种形式的实际示例(距离未按比例):

我想知道居住在 X 点半径范围内的人口数量。我的数据库发现它的 A 点和 B 点的条目与 X 点足够接近,因此具有相关性。此示例中的点 A 类似于 40.7458、-74.0375,点 B 类似于 40.7458、-74.0292。从 A 和 B 到其网格边缘的每条绿线代表 0.5 km,因此 A 和 B 周围的灰色圆圈分别代表 1 km^2。

X 点大约在 40.744,-74.032,半径(紫色)为 0.05 公里。

现在我可以使用地理三角函数轻松计算显示的红线。所以我知道 AX 线大约是 0.504 公里,距离线 BX 大约是 0.309 公里,不管我怎么想。

所以我的问题是:计算 X 周围的紫色圆圈所占网格 A 的分数和网格 B 的分数的可靠方法是什么?

最终我将计算人口总数并将它们乘以这个分数。所以在这种情况下,周围1km^2的格子对应的是9561人,B周围的格子是10763人。因此,如果我知道(只是假设地)X 周围的半径覆盖了 A 面积的 1% 和 B 面积的 3%,我就可以对该覆盖的总人口做出合理的粗略估计通过将 A 和 B 人口乘以各自的分数然后将它们相加来圈出。

我只用上面的两个正方形完成了,但是根据半径的大小(可以是任意的),可能会有很多可能的正方形,就像这样,使它成为一个更普遍的问题:

在某些情况下,很容易找出所讨论的正方形网格 100% 被半径包围,原则上这很容易(例如,如果 AX 之间的距离小于 X 周围的半径,我知道我不需要做任何进一步的数学运算)。

现在,很容易找出哪些点在圆的范围内。但我有点想弄清楚它们相应区域的分数是多少。

感谢您的帮助。

有了足够的分类(如下所示),所有计算都可以简化为 原始 计算,即提供 angular 区域[=图像中描绘的橙色区域的 41=]

y0 > 0时,如上图,不管x0是正还是负,都可以准确计算出橙色区域为x0到[=13的积分=] 的 sqrt(r^2 - y^2) 减去 rectangular 面积 (x1 - x0) * (y1 - y0). 积分有一个众所周知的封闭表达式,因此不需要使用任何数值算法来计算它。

圆形和正方形之间的其他交点可以简化为矩形和右 angular 形状的组合,如上面涂成橙色的那个。例如,下图中由水平和垂直橙色光线分隔的交点可以通过将红色矩形的面积加上两个 angular 形状:蓝色和绿色来表示。

蓝色区域是直接应用上面确定的原始情况(劣质矩形坍塌成空)的结果。绿色区域也可以用同样的方式测量,一旦负 y 坐标被它的绝对值取代(另一个 y0)。

应用这些想法可以列举所有情况。基本上,应该考虑这样一种情况,即正方形的一个、两个、三个或四个角位于圆内,而其余(如果有的话)位于圆外。枚举本身就是一个问题,但至少在理论上可以通过考虑相对较少的"typical"个配置来解决。

对于所描述的每种情况,必须计算一些矩形和 angular 区域的分解,并将部分加起来(或减去),如上面的三色示例所示。每个部分的面积将减少到 rectangular 或 primitive angular areas.

要将这条攻击路线转化为可行的算法,还需要做大量工作。更深入的分析可以阐明如何最大限度地减少要考虑的 "typical" 配置的数量。如果不是,我认为要考虑的组合数量,无论多大,都应该是可控的。

如果您的问题有一个大概的答案,您可以使用另一种编程更简单的技术。这道题的思路全部简化为计算正方形和圆形的交点面积。我在另一个答案中没有解释这一点,但是找到可能截取圆的正方形应该不是问题,否则请告诉我们。

计算路口大概面积的思路很简单。在正方形中随机生成足够多的点,并检查其中有多少点属于圆圈。圆中的点数与正方形中随机点总数之间的比率将给出交点相对于正方形面积的比例。

现在,鉴于您必须对圆周围的所有正方形重复相同的例程(即中心到圆心的距离与圆的半径相差不大的正方形),您可以重新使用随机通过将它们从一个正方形转换到另一个正方形来点。

如果此方法不适用于您的问题,我不想详细说明,所以我只想说明生成均匀分布在正方形中的随机点是相当容易的。您只需要为 x 坐标生成随机数,并独立地为 y 生成随机数。然后只考虑所有对 (x, y)。然后,对于每个(x, y)验证是否(x - a)^2 + (y - b)^2 <= r^2,其中(a, b)代表圆心,r代表半径。

我最终想出了一个非常好的近似解决方案,我认为。这是它在 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 度角。所以它大致计算出它落在这些极端之间的位置,然后使用它来计算大致覆盖的网格正方形的百分比。它在我的测试中给出了合理的结果。

对于我使用的正方形和圆形的大小,这似乎足够了?

我在 Processing.js 中编写了一个小应用程序来尝试帮助我解决这个问题。无需全部解释,您可以通过查看此屏幕截图了解算法 "thinking":

基本上,如果圆圈是黄色的,则意味着它已经确定它是 100% 进入的,如果是红色的,则表示它已经被快速筛选为 100% 退出。其他情况是边缘情况。点下方的数字(范围从 0 到 1)是使用上述方法计算的(四舍五入的)覆盖率百分比,而其下方的数字是上述代码中使用的计算出的 theta 值。

就我的目的而言,我认为这种近似是可行的。