Skip to content

Commit 2050861

Browse files
committed
add geometry intersects algorithm
1 parent 56968fd commit 2050861

File tree

3 files changed

+366
-90
lines changed

3 files changed

+366
-90
lines changed

src/commonMain/kotlin/com/jillesvangurp/geo/GeoGeometry.kt

Lines changed: 39 additions & 88 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,8 @@
2424
package com.jillesvangurp.geo
2525

2626
import com.jillesvangurp.geojson.*
27+
import com.jillesvangurp.geojson.latitude
28+
import com.jillesvangurp.geojson.longitude
2729
import kotlin.math.*
2830

2931

@@ -220,6 +222,7 @@ class GeoGeometry {
220222
return polygonContains(latitude, longitude, *polygonCoordinatesPoints[0])
221223
}
222224

225+
223226
/**
224227
* Determine whether a point is contained in a polygon. Note, technically
225228
* the points that make up the polygon are not contained by it.
@@ -231,102 +234,50 @@ class GeoGeometry {
231234
* [longitude,latitude]
232235
* @return true if the polygon contains the coordinate
233236
*/
234-
fun polygonContains(latitude: Double, longitude: Double, vararg polygonPoints: PointCoordinates): Boolean {
237+
fun polygonContains(
238+
latitude: Double,
239+
longitude: Double,
240+
vararg polygonPoints: PointCoordinates
241+
): Boolean {
242+
fun wrapLongitude(diff: Double): Double = when {
243+
diff > 180 -> diff - 360
244+
diff < -180 -> diff + 360
245+
else -> diff
246+
}
235247
validate(latitude, longitude, false)
248+
require(polygonPoints.size >= 3) { "a polygon must have at least three points" }
236249

237-
if (polygonPoints.size <= 2) {
238-
throw IllegalArgumentException("a polygon must have at least three points")
239-
}
240-
@Suppress("UNCHECKED_CAST") val bbox = boundingBox(polygonPoints as Array<DoubleArray>)
241-
if (!bboxContains(bbox, latitude, longitude)) {
242-
// outside the containing bbox
250+
251+
// 1) Shift longitudes so test-point’s lon → 0, deals with antimeridian
252+
val norm = polygonPoints.map {
253+
doubleArrayOf(
254+
wrapLongitude(it[0] - longitude),
255+
it[1]
256+
)
257+
}.toTypedArray()
258+
259+
// 2) Quick bbox check around (0, latitude)
260+
val bbox = boundingBox(norm)
261+
if (!bboxContains(bbox, latitude, 0.0)) {
243262
return false
244263
}
264+
if(polygonPoints.firstOrNull() { it.latitude ==latitude && it.longitude==longitude } != null) {
265+
// edge case for vertices
266+
return true
267+
}
245268

269+
// 3) Ray-cast from (0,latitude) eastward
246270
var hits = 0
247-
248-
var lastLatitude = polygonPoints[polygonPoints.size - 1][1]
249-
var lastLongitude = polygonPoints[polygonPoints.size - 1][0]
250-
var currentLatitude: Double
251-
var currentLongitude: Double
252-
253-
// Walk the edges of the polygon
254-
var i = 0
255-
while (i < polygonPoints.size) {
256-
currentLatitude = polygonPoints[i][1]
257-
currentLongitude = polygonPoints[i][0]
258-
259-
if (currentLongitude == lastLongitude) {
260-
lastLatitude = currentLatitude
261-
lastLongitude = currentLongitude
262-
i++
263-
continue
271+
val y = latitude
272+
for (i in norm.indices) {
273+
val (x1, y1) = norm[i]
274+
val (x2, y2) = norm[(i + 1) % norm.size]
275+
if ((y1 > y) != (y2 > y)) {
276+
val xInt = x1 + (y - y1) * (x2 - x1) / (y2 - y1)
277+
if (xInt > 0) hits++
264278
}
265-
266-
val leftLatitude: Double
267-
if (currentLatitude < lastLatitude) {
268-
if (latitude >= lastLatitude) {
269-
lastLatitude = currentLatitude
270-
lastLongitude = currentLongitude
271-
i++
272-
continue
273-
}
274-
leftLatitude = currentLatitude
275-
} else {
276-
if (latitude >= currentLatitude) {
277-
lastLatitude = currentLatitude
278-
lastLongitude = currentLongitude
279-
i++
280-
continue
281-
}
282-
leftLatitude = lastLatitude
283-
}
284-
285-
val test1: Double
286-
val test2: Double
287-
if (currentLongitude < lastLongitude) {
288-
if (longitude < currentLongitude || longitude >= lastLongitude) {
289-
lastLatitude = currentLatitude
290-
lastLongitude = currentLongitude
291-
i++
292-
continue
293-
}
294-
if (latitude < leftLatitude) {
295-
hits++
296-
lastLatitude = currentLatitude
297-
lastLongitude = currentLongitude
298-
i++
299-
continue
300-
}
301-
test1 = latitude - currentLatitude
302-
test2 = longitude - currentLongitude
303-
} else {
304-
if (longitude < lastLongitude || longitude >= currentLongitude) {
305-
lastLatitude = currentLatitude
306-
lastLongitude = currentLongitude
307-
i++
308-
continue
309-
}
310-
if (latitude < leftLatitude) {
311-
hits++
312-
lastLatitude = currentLatitude
313-
lastLongitude = currentLongitude
314-
i++
315-
continue
316-
}
317-
test1 = latitude - lastLatitude
318-
test2 = longitude - lastLongitude
319-
}
320-
321-
if (test1 < test2 / (lastLongitude - currentLongitude) * (lastLatitude - currentLatitude)) {
322-
hits++
323-
}
324-
lastLatitude = currentLatitude
325-
lastLongitude = currentLongitude
326-
i++
327279
}
328-
329-
return hits % 2 == 1
280+
return hits and 1 == 1
330281
}
331282

332283
/**
Lines changed: 100 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,100 @@
1+
package com.jillesvangurp.geojson
2+
3+
import com.jillesvangurp.geo.GeoGeometry
4+
5+
fun Geometry.intersects(other: Geometry): Boolean {
6+
val bbox1 = this.boundingBox()
7+
val bbox2 = other.boundingBox()
8+
if (!bboxesIntersect(bbox1, bbox2)) return false
9+
10+
return when (this) {
11+
is Geometry.Point -> other.contains(this.coordinates ?: return false)
12+
is Geometry.MultiPoint -> this.coordinates?.any { other.contains(it) } == true
13+
is Geometry.LineString -> intersectsLine(this.coordinates ?: return false, other)
14+
is Geometry.MultiLineString -> this.coordinates?.any { intersectsLine(it, other) } == true
15+
is Geometry.Polygon -> intersectsPolygon(this.coordinates ?: return false, other)
16+
is Geometry.MultiPolygon -> this.coordinates?.any { intersectsPolygon(it, other) } == true
17+
is Geometry.GeometryCollection -> this.geometries.any { it.intersects(other) }
18+
}
19+
}
20+
21+
private fun bboxesIntersect(b1: BoundingBox, b2: BoundingBox): Boolean {
22+
// Pick the shorter east-going arc for each raw interval
23+
fun normalize(w: Double, e: Double): Pair<Double, Double> {
24+
// span measured modulo 360
25+
val span = ( (e - w + 360) % 360 )
26+
return if (span > 180) e to w else w to e
27+
}
28+
29+
fun lonIntersects(w1: Double, e1: Double, w2: Double, e2: Double): Boolean {
30+
val (nW1, nE1) = normalize(w1, e1)
31+
val (nW2, nE2) = normalize(w2, e2)
32+
// now split only if it still wraps
33+
val arcs1 = if (nW1 <= nE1) listOf(nW1 to nE1)
34+
else listOf(nW1 to 180.0, -180.0 to nE1)
35+
val arcs2 = if (nW2 <= nE2) listOf(nW2 to nE2)
36+
else listOf(nW2 to 180.0, -180.0 to nE2)
37+
38+
return arcs1.any { (a0, a1) ->
39+
arcs2.any { (b0, b1) ->
40+
!(a0 > b1 || a1 < b0)
41+
}
42+
}
43+
}
44+
45+
val lonOverlap = lonIntersects(
46+
b1.westLongitude, b1.eastLongitude,
47+
b2.westLongitude, b2.eastLongitude
48+
)
49+
val latOverlap = !(b1.northLatitude < b2.southLatitude ||
50+
b1.southLatitude > b2.northLatitude)
51+
52+
return lonOverlap && latOverlap
53+
}
54+
55+
private fun intersectsLine(line: LineStringCoordinates, other: Geometry): Boolean {
56+
val segments = line.zipWithNextCompat()
57+
return segments.any { (start, end) ->
58+
when (other) {
59+
is Geometry.Point -> other.coordinates?.onLineSegment(start, end) ?: false
60+
is Geometry.LineString -> other.coordinates?.zipWithNextCompat()?.any { (oStart, oEnd) ->
61+
GeoGeometry.linesCross(start, end, oStart, oEnd)
62+
} == true
63+
is Geometry.MultiLineString -> other.coordinates?.any { oLine ->
64+
oLine.zipWithNextCompat().any { (oStart, oEnd) ->
65+
GeoGeometry.linesCross(start, end, oStart, oEnd)
66+
}
67+
} == true
68+
is Geometry.Polygon -> other.outerCoordinates.zipWithNextCompat().any { (oStart, oEnd) ->
69+
GeoGeometry.linesCross(start, end, oStart, oEnd)
70+
} || other.contains(start)
71+
is Geometry.MultiPolygon -> other.coordinates?.any { poly ->
72+
poly.first().zipWithNextCompat().any { (oStart, oEnd) ->
73+
GeoGeometry.linesCross(start, end, oStart, oEnd)
74+
} || GeoGeometry.polygonContains(start.latitude, start.longitude, poly)
75+
} == true
76+
is Geometry.GeometryCollection -> other.geometries.any { intersectsLine(line, it) }
77+
else -> false
78+
}
79+
}
80+
}
81+
82+
private fun intersectsPolygon(polygon: PolygonCoordinates, other: Geometry): Boolean {
83+
val outer = polygon.first()
84+
return outer.zipWithNextCompat().any { (start, end) ->
85+
intersectsLine(arrayOf(start, end), other)
86+
} || when (other) {
87+
is Geometry.Point -> other.coordinates?.let { GeoGeometry.polygonContains(it.latitude, it.longitude, polygon) } ?: false
88+
is Geometry.MultiPoint -> other.coordinates?.any { GeoGeometry.polygonContains(it.latitude, it.longitude, polygon) } == true
89+
else -> false
90+
}
91+
}
92+
93+
private fun <T> Array<T>.zipWithNextCompat(): List<Pair<T, T>> {
94+
// multiplatform doesn't have zipWithNext
95+
val result = mutableListOf<Pair<T, T>>()
96+
for (i in 0 until this.size - 1) {
97+
result.add(this[i] to this[i + 1])
98+
}
99+
return result
100+
}

0 commit comments

Comments
 (0)