确定球体上的经纬度矩形和圆是否重叠
假设我有以下内容:
- 由最小和最大纬度和经度定义的区域(通常是“经纬度矩形”,尽管除了某些投影之外它实际上不是矩形)。
- 由中心纬度/经度和半径定义的圆
如何确定:
- 两个形状是否重叠?
- 圆是否完全包含在矩形内?
我正在寻找一个完整的公式/算法,而不是数学课本身。
Suppose I have the following:
- A region defined by minimum and maximum latitude and longitude (commonly a 'lat-long rect', though it's not actually rectangular except in certain projections).
- A circle, defined by a center lat/long and a radius
How can I determine:
- Whether the two shapes overlap?
- Whether the circle is entirely contained within the rect?
I'm looking for a complete formula/algorithm, rather than a lesson in the math, per-se.
如果你对这篇内容有疑问,欢迎到本站社区发帖提问 参与讨论,获取更多帮助,或者扫码二维码加入 Web 技术交流群。
绑定邮箱获取回复消息
由于您还没有绑定你的真实邮箱,如果其他用户或者作者回复了您的评论,将不能在第一时间通知您!
发布评论
评论(7)
警告:如果圆/“矩形”跨越球体的大部分,例如:
“矩形”:最小长 = -90deg,最大长 = +90deg,最小纬度 = +70deg ,最大纬度 = +80 度
圆:中心 = 纬度 = +85 度,长 = +160 度,半径 = 20 度(例如,如果 A 点在圆上,C 点是圆的中心,O 点是球体的中心,则角度AOC = 40 度)。
它们相交,但数学上可能有几种情况来检查相交/包含。 以下点位于上述圆上:P1=(+65 度纬度,+160 度长),P2=(+75 度纬度,-20 度长)。 P1 在“矩形”外部,P2 在“矩形”内部,因此圆/“矩形”至少有 2 个点相交。
好的,这是我对解决方案的概述:
设 C = 半径为 R 的圆心(如上所示表示为球面角)。 C 具有纬度 LATC 和经度 LONGC。 由于“矩形”这个词在这里有点误导(恒定纬度的线不是大圆的线段),我将使用术语“边界框”。
函数
InsideCircle(P)
返回 +1、0 或 -1: 如果点 P 在圆内,则返回 +1;如果点 P 在圆上,则返回 0 ,如果 P 点在圆外,则为 -1:计算 C 和任意点 P 之间的大圆距离 D(表示为球面角)将告诉您 P 是否在圆内:InsideCircle(P) =sign(RD)
(正如 Sente 中的用户 @Die 提到的,在这个论坛的其他地方已经询问过大圆距离)定义
PANG(x)
= x 的主角 = MOD(x+180deg, 360deg)-180deg。PANG(x)
始终在 -180deg 和 +180deg 之间(包含 -180deg 和 +180deg)(+180deg 应映射到 -180deg)。要定义边界框,您需要知道 4 个数字,但经度有一个小问题。 LAT1和LAT2表示边界纬度(假设LAT1<LAT2); 那里没有任何歧义。 LONG1 和 LONG2 表示经度间隔的边界经度,但这会很棘手,并且更容易将此间隔重写为中心和宽度,其中 LONGM = 该间隔的中心,LONGW = 宽度。 请注意,经度间隔始终有 2 种可能性。 您必须指定是哪种情况,无论是包括还是排除 180 度子午线,例如,从 -179 度到 +177 度的最短间隔为 LONGM = +179 度和 LONGW = 4 度,但从 -179 度到 + 的其他间隔177deg 的 LONGM = -1deg 且 LONGW = 356deg。 如果您天真地尝试与区间 [-179,177] 进行“常规”比较,您最终将使用更大的区间,而这可能不是您想要的。 顺便说一句,如果以下两个条件都成立,则具有纬度 LATP 和经度 LONGP 的点 P 位于边界框内:
如果 PTEST = union(PCORNER,PLAT,PLONG) 中的以下任何点 P 如下所述,则圆与边界框相交,但
InsideCircle()
并不都返回相同的结果:上面列出的这些点 PLAT 和 PLONG 是边界框上距离圆“最近”的点(如果角不是;我在引号中使用“最近”,从纬度/长距离的意义上来说,不是很大-圆距离),并涵盖圆的中心位于边界框边界的一侧但圆上的点“潜行穿过”边界框边界的情况。
如果 PTEST 中的所有点 P 返回
InsideCircle(P)
== +1(全部在圆内),则该圆包含整个边界框。如果 PTEST 中的所有点 P 返回
InsideCircle(P)
== -1(全部在圆之外),则圆完全包含在边界框内。否则,圆和边界框之间至少有一个交点。 请注意,这不会计算这些点的位置,尽管如果您在 PTEST 中取任意 2 个点 P1 和 P2,其中 InsideCircle(P1) = -InsideCircle(P2),那么您可以通过二等分找到交点(效率低下)。 (如果 InsideCircle(P) 返回 0,则有一个交点,尽管浮点数学中的等式通常不可信。)
可能有一种更有效的方法可以做到这一点,但上述方法应该可行。
warning: this can be tricky if the circles / "rectangles" span large portions of the sphere, e.g.:
"rectangle": min long = -90deg, max long = +90deg, min lat = +70deg, max lat = +80deg
circle: center = lat = +85deg, long = +160deg, radius = 20deg (e.g. if point A is on the circle and point C is the circle's center, and point O is the sphere's center, then angle AOC = 40deg).
These intersect but the math is likely to have several cases to check intersection/containment. The following points lie on the circle described above: P1=(+65deg lat,+160deg long), P2=(+75deg lat, -20deg long). P1 is outside the "rectangle" and P2 is inside the "rectangle" so the circle/"rectangle" intersect in at least 2 points.
OK, here's my shot at an outline of the solution:
Let C = circle center with radius R (expressed as a spherical angle as above). C has latitude LATC and longitude LONGC. Since the word "rectangle" is kind of misleading here (lines of constant latitude are not segments of great circles), I'll use the term "bounding box".
function
InsideCircle(P)
returns +1,0,or -1: +1 if point P is inside the circle, 0 if point P is on the circle, and -1 if point P is outside the circle: calculation of great-circle distance D (expressed as spherical angle) between C and any point P will tell you whether or not P is inside the circle:InsideCircle(P) = sign(R-D)
(As user @Die in Sente mentioned, great circle distances have been asked on this forum elsewhere)Define
PANG(x)
= the principal angle of x = MOD(x+180deg, 360deg)-180deg.PANG(x)
is always between -180deg and +180deg, inclusive (+180deg should map to -180deg).To define the bounding box, you need to know 4 numbers, but there's a slight issue with longitude. LAT1 and LAT2 represent the bounding latitudes (assuming LAT1 < LAT2); there's no ambiguity there. LONG1 and LONG2 represent the bounding longitudes of a longitude interval, but this gets tricky, and it's easier to rewrite this interval as a center and width, with LONGM = the center of that interval and LONGW = width. NOTE that there are always 2 possibilities for longitude intervals. You have to specify which of these cases it is, whether you are including or excluding the 180deg meridian, e.g. the shortest interval from -179deg to +177deg has LONGM = +179deg and LONGW = 4deg, but the other interval from -179deg to +177deg has LONGM = -1deg and LONGW = 356deg. If you naively try to do "regular" comparisons with the interval [-179,177] you will end up using the larger interval and that's probably not what you want. As an aside, point P, with latitude LATP and longitude LONGP, is inside the bounding box if both of the following are true:
The circle intersects the bounding box if ANY of the following points P in PTEST = union(PCORNER,PLAT,PLONG) as described below, do not all return the same result for
InsideCircle()
:These points PLAT and PLONG as listed above are the points on the bounding box that are "closest" to the circle (if the corners are not; I use "closest" in quotes, in the sense of lat/long distance and not great-circle distance), and cover the cases where the circle's center lies on one side of the bounding box's boundary but points on the circle "sneak across" the bounding box boundary.
If all points P in PTEST return
InsideCircle(P)
== +1 (all inside the circle) then the circle contains the bounding box in its entirety.If all points P in PTEST return
InsideCircle(P)
== -1 (all outside the circle) then the circle is contained entirely within the bounding box.Otherwise there is at least one point of intersection between the circle and the bounding box. Note that this does not calculate where those points are, although if you take any 2 points P1 and P2 in PTEST where InsideCircle(P1) = -InsideCircle(P2), then you could find a point of intersection (inefficiently) by bisection. (If InsideCircle(P) returns 0 then you have a point of intersection, though equality in floating-point math is generally not to be trusted.)
There's probably a more efficient way to do this but the above should work.
使用立体投影。 所有圆(特别是纬度、经度和您的圆)都映射到平面中的圆(或线)。 现在这只是一个关于平面几何中的圆和线的问题(更好的是,所有的经线都是经过 0 的线,所有的纬度都是围绕 0 的圆)
Use the Stereographic projection. All circles (specifically latitudes, longitudes and your circle) map to circles (or lines) in the plane. Now it's just a question about circles and lines in plane geometry (even better, all the longitues are lines through 0, and all the latitudes are circles around 0)
假设:
第一次检查很简单。 第二次检查只需要找到四个距离。 第三次检查只需要找到从圆心到(最近框纬度,圆心经度)的距离。
第四个检查需要找到最接近圆心的边界框的经度线。 然后找到距离圆心最远的经线所在的大圆的中心。 求从圆心到大圆心的初始方位角。 找到该轴承上距圆心的点圆半径。 如果该点位于距离圆中心最近的经度线的另一侧,则圆和边界框在该侧相交。
我觉得这应该是有缺陷的,但是我一直没能发现。
我似乎无法解决的真正问题是找到完美包含圆的边界框(对于不包含极点的圆)。 最小/最大纬度的方位似乎是圆心和圆半径/(球体周长/4)纬度的函数。 在赤道附近,它下降到 pi/2(东)或 3*pi/2(西)。 当中心接近极点且半径接近球周长/4 时,方位角接近零(北)或 pi(南)。
Assumptions:
The first check is trivial. The second check just requires finding the four distances. The third check just requires finding the distance from circle-center to (closest-box-latitude, circle-center-longitude).
The fourth check requires finding the longitude line of the bounding box that is closest to the circle-center. Then find the center of the great circle on which that longitude line rests that is furthest from circle-center. Find the initial-bearing from circle-center to the great-circle-center. Find the point circle-radius from circle-center on that bearing. If that point is on the other side of the closest-longitude-line from circle-center, then the circle and bounding box intersect on that side.
It seems to me that there should be a flaw in this, but I haven't been able to find it.
The real problem that I can't seem to solve is to find the bounding-box that perfectly contains the circle (for circles that don't contain a pole). The bearing to the latitude min/max appears to be a function of the latitude of circle-center and circle-radius/(sphere circumference/4). Near the equator, it falls to pi/2 (east) or 3*pi/2 (west). As the center approaches the pole and the radius approaches sphere-circumference/4, the bearing approach zero (north) or pi (south).
这个怎么样?
找到连接矩形中心(点
Cr
)和圆心的向量v
。 找到v
与矩形相交的点i
。 如果<代码>||i-Cr|| + r > ||v|| 然后它们相交。换句话说,矩形内部线段的长度加上圆形内部线段的长度应该大于(
v
,中心连接线段)的总长度。找到点
i
应该是棘手的部分,特别是如果它落在经度边缘上,但你应该能够比我更快地找到一些东西。编辑:此方法无法判断圆是否完全在矩形内。 为此,您需要找到从其中心到矩形所有四个边缘的距离。
编辑:上面是不正确的。正如 Federico Ramponi 所建议的那样,在某些情况下,即使在欧几里得几何中它也不起作用。 我会发布另一个答案。 请不接受这一点并随时投反对票。 我很快就会删除它。
How about this?
Find vector
v
that connects the center of the rectangle, pointCr
, to the center of the circle. Find pointi
wherev
intersects the rectangle. If||i-Cr|| + r > ||v||
then they intersect.In other words, the length of the segment inside the rectangle plus the length of the segment inside the circle should be greater than the total length (of
v
, the center-connecting line segment).Finding point
i
should be the tricky part, especially if it falls on a longitude edge, but you should be able to come up with something faster than I can.Edit: This method can't tell if the circle is completely within the rectangle. You'd need to find the distance from its center to all four of the rectangle's edges for that.
Edit: The above is incorrect. There are some cases, as Federico Ramponi has suggested, where it does not work even in Euclidean geometry. I'll post another answer. Please unaccept this and feel free to vote down. I'll delete it shortly.
这应该适用于地球上的任何点。 如果您想将其更改为不同大小的球体,只需将 kEarchRadiusKms 更改为您想要的球体半径即可。
该方法用于计算纬度点和经度点之间的距离。
我从这里得到了这个距离公式:
http://www.codeproject.com/csharp/distance Betweenlocations.asp
如果距离矩形任意顶点之间的距离小于圆的半径则圆和矩形重叠。 如果圆心和所有顶点之间的距离大于圆的半径,并且所有这些距离都小于矩形的宽度和高度,则圆应该在矩形的内部。
如果您发现问题,请随时更正我的代码,因为我确信存在一些我没有想到的情况。
另外,我不确定这是否适用于跨越半球两端的矩形,因为距离方程可能会崩溃。
This should work for any points on earth. If you want to change it to a different size sphere just change the kEarchRadiusKms to whatever radius you want for your sphere.
This method is used to calculate the distance between to lat and lon points.
I got this distance formula from here:
http://www.codeproject.com/csharp/distancebetweenlocations.asp
If the distance between any vertex of the rectangle is less than the distance of the radius of the circle then the circle and rectangle overlap. If the distance between the center of the circle and all of the vertices is greater than the radius of the circle and all of those distances are shorter than the width and height of the rectangle then the circle should be inside of the rectangle.
Feel free to correct my code if you can find a problem with it as I'm sure there some condition that I have not thought of.
Also I'm not sure if this works for a rectangle that spans the ends of the hemispheres as the distance equation might break down.
再尝试一次......
我认为解决方案是测试一组点,正如 Jason S 所建议的那样,但我不同意他选择的点,我认为这在数学上是错误的。
您需要找到纬度/经度框两侧到圆心的距离为局部最小值或最大值的点。 将这些点添加到角点集中,然后上面的算法应该是正确的。
即,让经度为 x 维度,纬度为 y 维度,让每个维度
盒子的一侧是参数曲线 P(t) = P0 + t (P1-P0) for o <= t <= 1.0,其中
P0和P1是两个相邻的角。
令 f(P) = f(Px, Py) 为距圆心的距离。
那么f(P0+t(P1-P0))就是t的距离函数:g(t)。 找到距离函数的导数为零的所有点:g'(t) == 0。(当然,丢弃解决方案超出了域 0 <= t <= 1.0)不幸的是
,这需要找到零的先验表达式,因此没有封闭形式的解决方案。 此类方程只能通过牛顿-拉夫逊迭代来求解。
好吧,我知道你想要的是代码,而不是数学。 但数学是我所拥有的一切。
One more try at this...
I think the solution is to test a set of points, just as Jason S has suggested, but I disagree with his selection of points, which I think is mathematically wrong.
You need to find the points on the sides of the lat/long box where the distance to the center of the circle is a local minimum or maximum. Add those points to the set of corners and then the algorithm above should be correct.
I.e, letting longitude be the x dimension and latitude be the y dimension, let each
side of the box be a parametric curve P(t) = P0 + t (P1-P0) for o <= t <= 1.0, where
P0 and P1 are two adjacent corners.
Let f(P) = f(P.x, P.y) be the distance from the center of the circle.
Then f (P0 + t (P1-P0)) is a distance function of t: g(t). Find all the points where the derivative of the distance function is zero: g'(t) == 0. (Discarding solutions outsize the domain 0 <= t <= 1.0, of course)
Unfortunately, this needs to find the zero of a transcendental expression, so there's no closed form solution. This type of equation can only solved by Newton-Raphson iteration.
OK, I realize that you wanted code, not the math. But the math is all I've got.
有关欧几里得几何答案,请参阅:圆-矩形碰撞检测(交集)
For the Euclidean geometry answer, see: Circle-Rectangle collision detection (intersection)