We present a one-line algorithm that tests for the intersection of two segments of a circle. We use this technique to check if two rectangular areas on the globe intersect. Although the problem is deviously bewildering, the algorithm is unbelievably simple. A stored SQL procedure implementation and a sample GIS SQL query illustrate the algorithm. The algorithm was derived on Fri, 28 Aug 1998 01:21:34 GMT.

Testing two "rectangular" areas on the globe for intersection
is markedly more complicated than it may appear. The challenge comes
from longitude's wrapping around, that is, abruptly jumping from 180
degrees (i.e., `180E`

) to -180 degrees when we cross the
International Date Line. Latitudes do not wrap around. Therefore, in
the body of this article, we concentrate on a one-dimensional problem
of checking two segments of a circle for intersection. We will
generalize to georectangles and geographical information systems in
the last section.

To identify points on a circle, we will pick one point as the
origin and will refer to the other points by their angular distance
: "longitudes". The longitudes are measured in degrees,
counter-clockwise (that is, "East-wise") from the point of origin. We
will keep the longitudes within the range `[-180,180)`

, by
adding or subtracting 360 if necessary.

A segment `X`

of a circle is specified by its left
end, `L`

and its right end, _{X}`R`

. In other words, a
segment
_{X}`[`

starts at the point `L`

, _{X}`R`

]_{X}`L`

and
stretches counter-clockwise to the point _{X}`R`

. The coordinate
-- longitude -- of the point _{X}`L`

may be less, equal or
_{X}*greater*
than that of the point `R`

. We assume that
all segments are shorter than the full circle. We call two segments intersecting if one segment contains an _{X}*internal*
point of the other.

The following examples illustrate how non-obvious the intersection testing is.

Segment A | Segment B | Intersect? |
---|---|---|

`[1, 10]` |
`[15, 20]` |
NO |

`[10, 1]` |
`[15, 20]` |
YES |

`[5, 6]` |
`[4, 7]` |
YES |

`[6, 5]` |
`[4, 7]` |
still YES |

`[178, -178]` |
`[179, -179]` |
YES |

`[-178, 178]` |
`[179, -179]` |
NO |

`[178, -178]` |
`[-170, -100]` |
NO |

`[-178, 178]` |
`[-170, -100]` |
YES |

It is quite difficult to devise a single inequality involving
`L`

, _{A}`L`

, _{B}`R`

, and _{A}`R`

that holds if and only if the segments intersect._{B}

Suppose we are given two segments: segment A `[`

and segment B `L`

, _{A}`R`

]_{A}`[`

. We assume that the longitudes, `L`

, _{B}`R`

]_{B}`L`

, _{A}`R`

, _{A}`L`

, and _{B}`R`

, are the numbers within [-180,180)._{B}

diff(where`R`

-_{A}`L`

) + diff(_{A}`L`

-_{B}`R`

) + diff(_{A}`R`

-_{B}`L`

) + diff(_{B}`L`

-_{A}`R`

) = 360 <=> the segments do not intersect_{B}

`diff(x) = x mod 360`

normalizes `x`

to be within [0,360); that isdiff(x) = x if x >= 0; 360+x otherwise

The algorithm is indeed that simple. We merely need to compute four sets of differences, adding 360 to make all of them positive. If these four numbers add to 360, the segments in question do NOT intersect.

The segments do not intersect => the sum is equal to 360.

Suppose `[`

and `L`

, _{A}`R`

]_{A}`[`

do
not intersect. If we start at the point `L`

, _{B}`R`

]_{B}`L`

and walk to the East, we will
encounter the following points_{A}

in exactly that sequence. That means the segments`L`

,_{A}`R`

,_{A}`L`

,_{B}`R`

, back to_{B}`L`

_{A}

`[``L`_{A}

, `R`_{A}

]

, `[``R`_{A}

, `L`_{B}

]

, `[``L`_{B}

, `R`_{B}

]

and `[``R`_{B}

, `L`_{A}

]

are
all disjoint and together make up the full circle. The function `diff(x) `

is designed so that `diff(``R`_{X}

-`L`_{X}

)

is the length of a segment `[``L`_{X}

, `R`_{X}

]

. The (angular) length is a
number within [0,360). The sum `diff(``R`_{A}

-`L`_{A}

) + diff(`L`_{B}

-`R`_{A}

) + diff(`R`_{B}

-`L`_{B}

) + diff(`L`_{A}

-`R`_{B}

)

is then the length of the full
circle.The sum is equal to 360 => the segments do not intersect.

We will prove the second part of the correctness theorem by contradiction. Let us assume that the sum is equal to 360 but the segments still intersect. There are three major cases to consider: (i) segment B is entirely within segment A (or the other way around); (ii) Segment A contains `L`

but not _{B}`R`

; (iii) Segment A contains both endpoints _{B}`L`

and _{B}`R`

, but not the whole _{B}`[`

.`L`

, _{B}`R`

]_{B}

Case 1: Segment B is entirely within Segment A. If we walk from `L`

to the East we will pass the points _{A}`L`

, _{B}`R`

, _{B}`R`

, and come back to _{A}`L`

. Noting that _{A}`diff(-x) + diff(x) = 360`

for `x > 0`

, we write:

diff(Thus we have`L`

-_{B}`R`

) = 360 - diff(_{A}`R`

-_{A}`L`

) = 360 - diff(_{B}`R`

-_{B}`L`

) - diff(_{B}`R`

-_{A}`R`

) diff(_{B}`L`

-_{A}`R`

) = 360 - diff(_{B}`R`

-_{B}`L`

) = 360 - diff(_{A}`R`

-_{B}`L`

) - diff(_{B}`L`

-_{B}`L`

) diff(_{A}`R`

-_{A}`L`

) = diff(_{A}`L`

-_{B}`L`

) + diff(_{A}`R`

-_{B}`L`

) + diff(_{B}`R`

-_{A}`R`

)_{B}

diff(which is contrary to our assumption that this sum is 360.`R`

-_{A}`L`

) + diff(_{A}`L`

-_{B}`R`

) + diff(_{A}`R`

-_{B}`L`

) + diff(_{B}`L`

-_{A}`R`

) = 720_{B}

Case 2: Segment A contains `L`

but not _{B}`R`

. This implies that _{B}`R`

must be distinct from _{B}`L`

and _{A}`L`

must be distinct from _{B}`R`

. The same walk as in Case 1 gives us the following sequence of points: _{A}`L`

, _{A}`L`

, _{B}`R`

, _{A}`R`

, _{B}`L`

. Therefore _{A}

diff(which means that the sum`L`

-_{A}`R`

) = 360 - diff(_{B}`R`

-_{B}`L`

) = 360 - diff(_{A}`L`

-_{B}`L`

) - diff(_{A}`R`

-_{A}`L`

) - diff(_{B}`R`

-_{B}`R`

)_{A}

diff(and not 360 as we assumed initially.`R`

-_{A}`L`

) + diff(_{A}`L`

-_{B}`R`

) + diff(_{A}`R`

-_{B}`L`

) + diff(_{B}`L`

-_{A}`R`

) = 720_{B}

Case 3: Segment A contains both points `L`

and _{B}`R`

, but not all the points in between; for example: _{B}`[6, 5]`

and `[4, 7]`

. The walk to the East starting from `L`

goes through the points _{A}`L`

, _{A}`R`

, _{B}`L`

, _{B}`R`

and back to _{A}`L`

. Let us consider the equality _{A}

diff(The right-hand side of the equality looks almost the same as the sum we want to evaluate:`R`

-_{B}`L`

) + diff(_{A}`L`

-_{B}`R`

) + diff(_{B}`R`

-_{A}`L`

) + diff(_{B}`L`

-_{A}`R`

) = 360_{A}

diff(only all the differences are in "reverse":`R`

-_{A}`L`

) + diff(_{A}`L`

-_{B}`R`

) + diff(_{A}`R`

-_{B}`L`

) + diff(_{B}`L`

-_{A}`R`

) ,_{B}

`diff(``R`_{A}

-`L`_{A}

)

rather than `diff(``L`_{A}

-`R`_{A}

)

, etc. Keeping in mind that `diff(x) = 360 - diff(-x) `

for `x > 0`

, we obtaindiff(if all four end points are distinct. If`R`

-_{A}`L`

) + diff(_{A}`L`

-_{B}`R`

) + diff(_{A}`R`

-_{B}`L`

) + diff(_{B}`L`

-_{A}`R`

) = 1080_{B}

`L`_{A}

and `R`_{B}

are the same, the sum will be equal to 720. Ditto if `L`_{B}

and `R`_{A}

are the same. In any event, the result contradicts our initial assumption. Note: since we assumed that segments A and B intersect, one segment must contain an internal point of the other. This fact implies that `L`_{A}

and `R`_{B}

cannot be the same when `L`_{B}

and `R`_{A}

are the same.QED.

The intersection algorithm has been used in a production geographical information system. One of the databases stores satellite imagery. A satellite image is characterized, among other attributes, by its georectangular bounding box. A user can select an area of interest (again, as a georectangle) and find out relevant satellite images, i.e., the images whose bounding box intersects the area of interest. If the area of interest is a circle or a complex polygon, we can inscribe it into a georectangle and use the algorithm above. The test may in this case yield a false positive, but it will never return a false negative. Therefore, the georectangle intersection approach can be used as the first, winnowing stage in checking for intersection of complex areas on the globe.

The following is a query template excerpted from the production
code. The query returns the bounding boxes and the descriptions of all
satellite images that intersect a user-defined georectangle. The latter is
specified by its two corner-points: `(NLAT, WLON)`

and `(SLAT,ELON)`

.

SELECT lat_n, lat_s, lon_w, lon_e, descr FROM ImageDescr WHERE (SELECT image_descr FROM Timestamps) >= '2002-06-07 16:47' AND lat_s <= 'NLAT' AND lat_n >= 'SLAT' AND q_lon_intersect(lon_w,lon_e,WLON,ELON) <> 360 ORDER BY lat_n, lat_s, lon_w, lon_e FOR READ ONLY;

The query relies on the following stored procedure `q_lon_intersect()`

, which implements the segment intersection algorithm verbatim:

-- This q_lon_intersect(lon_w,lon_e,Query_WLON,Query_ELON) -- procedure determines if two circle segments [lon_w, lon_e] and -- [Query_WLON, Query_ELON] intersect. The procedure works correctly -- for lon_w < lon_e and lon_w >= lon_e (and similarly for Query_WLON -- and Query_ELON). -- The two circle segments intersect when the result of this procedure -- is NOT equal to 360. CREATE PROCEDURE q_lon_intersect (La SMALLFLOAT, Ra SMALLFLOAT, Lb SMALLFLOAT, Rb SMALLFLOAT) returning SMALLFLOAT; DEFINE d1, d2, d3, d4 SMALLFLOAT; LET d1, d2, d3, d4 = Ra-La, Lb-Ra, Rb-Lb, La-Rb; RETURN CASE WHEN d1 < 0 THEN 360 + d1 ELSE d1 END + CASE WHEN d2 < 0 THEN 360 + d2 ELSE d2 END + CASE WHEN d3 < 0 THEN 360 + d3 ELSE d3 END + CASE WHEN d4 < 0 THEN 360 + d4 ELSE d4 END; END PROCEDURE DOCUMENT 'Compute a SMALLFLOAT that can tell if two circle segments', 'intersect. If they do, the result of this procedure is', 'NOT equal to 360.' WITH LISTING IN '/tmp/sql.complile.log';

This site's top page is **http://okmij.org/ftp/**

Your comments, problem reports, questions are very welcome!

Converted from SXML by SXML->HTML