test and clean up line intersection code

Author: ffreyer

src/lines.jl

@@ -7,7 +7,7 @@ Returns `(intersection_found::Bool, intersection_point::Point)`
 """
 # 2D Line-segment intersection algorithm by Paul Bourke and many others.
 # http://paulbourke.net/geometry/pointlineplane/
-function intersects(a::Line{2,T1}, b::Line{2,T2}) where {T1,T2}
+function intersects(a::Line{2,T1}, b::Line{2,T2}; eps = 0) where {T1,T2}
     T = promote_type(T1, T2)
     p0 = zero(Point2{T})
 
@@ -34,7 +34,7 @@ function intersects(a::Line{2,T1}, b::Line{2,T2}) where {T1,T2}
 
     # Values between [0, 1] mean the intersection point of the lines rests on
     # both of the line segments.
-    if 0 <= unknown_a <= 1 && 0 <= unknown_b <= 1
+    if eps <= unknown_a <= 1-eps && eps <= unknown_b <= 1-eps
         # Substituting an unknown back lets us find the intersection point.
         x = x1 + (unknown_a * (x2 - x1))
         y = y1 + (unknown_a * (y2 - y1))
@@ -62,27 +62,29 @@ end
 """
     self_intersections(points::AbstractVector{<:Point})
 
-Finds all self intersections of polygon `points`
+Finds all self intersections of in a continuous line described by `points`.
+
+Note that if two points are the same, they will generate a self intersection
+unless they are the first and last point or part of consecutive segments.
 """
-function self_intersections(points::AbstractVector{<:Point})
+function self_intersections(points::AbstractVector{<:VecTypes{D, T}}) where {D, T}
     sections = similar(points, 0)
-    intersections = Int[]
-
-    wraparound(i) = mod1(i, length(points) - 1)
-
-    for (i, (a, b)) in enumerate(consecutive_pairs(points))
-        for (j, (a2, b2)) in enumerate(consecutive_pairs(points))
-            is1, is2 = wraparound(i + 1), wraparound(i - 1)
-            if i != j &&
-               is1 != j &&
-               is2 != j &&
-               !(i in intersections) &&
-               !(j in intersections)
-                intersected, p = intersects(Line(a, b), Line(a2, b2))
-                if intersected
-                    push!(intersections, i, j)
-                    push!(sections, p)
-                end
+    intersections = Tuple{Int, Int}[]
+
+    N = length(points)
+
+    for i in 1:N-3
+        a = points[i]; b = points[i+1]
+        # i+1 == j describes consecutive segments which are always "intersecting"
+        # at point i+1/j. Skip those (start at i+2)
+        # Special case: We assume points[1] == points[end] so 1 -> 2 and N-1 -> N
+        # always "intersect" at 1/N. Skip this too (end at N-2 in this case)
+        for j in i+2 : N-1 - (i == 1)
+            a2 = points[j]; b2 = points[j+1]
+            intersected, p = intersects(Line(a, b), Line(a2, b2))
+            if intersected
+                push!(intersections, (i, j))
+                push!(sections, p)
             end
         end
     end
@@ -95,15 +97,14 @@ end
 Splits polygon `points` into it's self intersecting parts. Only 1 intersection
 is handled right now.
 """
-function split_intersections(points::AbstractVector{<:Point})
+function split_intersections(points::AbstractVector{<:VecTypes{N, T}}) where {N, T}
     intersections, sections = self_intersections(points)
     return if isempty(intersections)
         return [points]
-    elseif length(intersections) == 2 && length(sections) == 1
-        a, b = intersections
+    elseif length(intersections) == 1 && length(sections) == 1
+        a, b = intersections[1]
         p = sections[1]
-        a, b = min(a, b), max(a, b)
-        poly1 = simple_concat(points, (a + 1):(b - 1), p)
+        poly1 = simple_concat(points, (a + 1):b, p)
         poly2 = simple_concat(points, (b + 1):(length(points) + a), p)
         return [poly1, poly2]
     else

test/runtests.jl

@@ -258,6 +258,8 @@ end
         found, point = intersects(d, e)
         @test found && point ≈ Point(0.0, 4.0)
         @test intersects(a, f) === (false, Point(0.0, 0.0))
+        @test intersects(a, d) === (true, Point(0.0, 0.0))
+        @test intersects(a, d, eps = 1e-6) === (false, Point(0.0, 0.0))
 
         # issue #168
         # If these tests fail then you can increase the tolerance on the checks so
@@ -273,6 +275,28 @@ end
         @test b
         @test loc ≈ Point(5760.0, 127.27922061357884)
     end
+
+    ps = [Point2f(1), Point2f(2)]
+    @test GeometryBasics.simple_concat(ps, 2:2, Point2f(3)) == [Point2f(2), Point2f(3)]
+    ps = [Point2f(i) for i in 1:4]
+    @test collect(GeometryBasics.consecutive_pairs(ps)) == collect(zip(ps[1:end-1], ps[2:end]))
+
+    ps = Point2f[(0,0), (1,0), (0,1), (1,2), (0,2), (1,1), (0,0)]
+    idxs, ips = self_intersections(ps)
+    @test idxs == [(2, 6), (3, 5)]
+    @test ips == [Point2f(0.5), Point2f(0.5, 1.5)]
+
+    ps = [Point2f(cos(x), sin(x)) for x in 0:4pi/5:4pi+0.1]
+    idxs, ips = self_intersections(ps)
+    @test idxs == [(1, 3), (1, 4), (2, 4), (2, 5), (3, 5)]
+    @test all(ips .≈ Point2f[(0.30901694, 0.2245140), (-0.118034005, 0.36327127), (-0.38196602, 0), (-0.118033946, -0.3632713), (0.309017, -0.22451389)])
+
+    @test_throws ErrorException split_intersections(ps)
+    ps = Point2f[(0,0), (1,0), (0,1), (1,1), (0, 0)]
+    idxs, ips = self_intersections(ps)
+    sps = split_intersections(ps)
+    @test sps[1] == [ps[3], ps[4], ips[1]]
+    @test sps[2] == [ps[5], ps[1], ps[2], ips[1]]
 end