Fix #168 (#185)

* Add broken tests for #168

* Fix #168

---------

Co-authored-by: Simon <[email protected]>

Author: cmcaine

src/lines.jl

@@ -2,66 +2,47 @@
 """
     intersects(a::Line, b::Line) -> Bool, Point
 
-Intersection of 2 line segmens `a` and `b`.
-Returns intersection_found::Bool, intersection_point::Point
+Intersection of 2 line segments `a` and `b`.
+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}
     T = promote_type(T1, T2)
-    v1, v2 = a
-    v3, v4 = b
-    MT = Mat{2,2,T,4}
     p0 = zero(Point2{T})
 
-    verticalA = v1[1] == v2[1]
-    verticalB = v3[1] == v4[1]
-
-    # if a segment is vertical the linear algebra might have trouble
-    # so we will rotate the segments such that neither is vertical
-    dorotation = verticalA || verticalB
-
-    if dorotation
-        θ = T(0.0)
-        if verticalA && verticalB
-            θ = T(π / 4)
-        elseif verticalA || verticalB # obviously true, but make it clear
-            θ34 = -atan(v4[2] - v3[2], v4[1] - v3[1])
-            θ12 = -atan(v2[2] - v1[2], v2[1] - v1[1])
-            θ = verticalA ? θ34 : θ12
-            θ = abs(θ) == T(0) ? (θ12 + θ34) / 2 : θ
-            θ = abs(θ) == T(pi) ? (θ12 + θ34) / 2 : θ
-        end
-        rotation = MT(cos(θ), sin(θ), -sin(θ), cos(θ))
-        v1 = rotation * v1
-        v2 = rotation * v2
-        v3 = rotation * v3
-        v4 = rotation * v4
-    end
+    x1, y1 = a[1]
+    x2, y2 = a[2]
+    x3, y3 = b[1]
+    x4, y4 = b[2]
 
-    a = det(MT(v1[1] - v2[1], v1[2] - v2[2], v3[1] - v4[1], v3[2] - v4[2]))
+    denominator = ((y4 - y3) * (x2 - x1)) - ((x4 - x3) * (y2 - y1))
+    numerator_a = ((x4 - x3) * (y1 - y3)) - ((y4 - y3) * (x1 - x3))
+    numerator_b = ((x2 - x1) * (y1 - y3)) - ((y2 - y1) * (x1 - x3))
 
-    (abs(a) < eps(T)) && return false, p0 # Lines are parallel
+    if denominator == 0
+        # no intersection: lines are parallel
+        return false, p0
+    end
 
-    d1 = det(MT(v1[1], v1[2], v2[1], v2[2]))
-    d2 = det(MT(v3[1], v3[2], v4[1], v4[2]))
-    x = det(MT(d1, v1[1] - v2[1], d2, v3[1] - v4[1])) / a
-    y = det(MT(d1, v1[2] - v2[2], d2, v3[2] - v4[2])) / a
+    # If we ever need to know if the lines are coincident, we can get that too:
+    # denominator == numerator_a == numerator_b == 0 && return :coincident_lines
 
-    (x < prevfloat(min(v1[1], v2[1])) || x > nextfloat(max(v1[1], v2[1]))) &&
-        return false, p0
-    (y < prevfloat(min(v1[2], v2[2])) || y > nextfloat(max(v1[2], v2[2]))) &&
-        return false, p0
-    (x < prevfloat(min(v3[1], v4[1])) || x > nextfloat(max(v3[1], v4[1]))) &&
-        return false, p0
-    (y < prevfloat(min(v3[2], v4[2])) || y > nextfloat(max(v3[2], v4[2]))) &&
-        return false, p0
+    # unknown_a and b tell us how far along the line segment the intersection is.
+    unknown_a = numerator_a / denominator
+    unknown_b = numerator_b / denominator
 
-    point = Point2{T}(x, y)
-    # don't forget to rotate the answer back
-    if dorotation
-        point = transpose(rotation) * point
+    # 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
+        # Substituting an unknown back lets us find the intersection point.
+        x = x1 + (unknown_a * (x2 - x1))
+        y = y1 + (unknown_a * (y2 - y1))
+        return true, Point2{T}(x, y)
     end
 
-    return true, point
+    # lines intersect, but outside of at least one of these line segments.
+    return false, p0
 end
 
 function simple_concat(vec::AbstractVector, range, endpoint::P) where {P}

test/runtests.jl

@@ -629,6 +629,20 @@ end
     found, point = intersects(d, e)
     @test found && point ≈ Point(0.0, 4.0)
     @test intersects(a, f) === (false, Point(0.0, 0.0))
+
+    # issue #168
+    # If these tests fail then you can increase the tolerance on the checks so
+    # long as you know what you're doing :)
+    line_helper(a, b, c, d) = Line(Point(a, b), Point(c, d))
+    b, loc = intersects(line_helper(-3.1, 15.588457268119894, 3.1, 15.588457268119894),
+                        line_helper(2.0866025403784354, 17.37050807568877, -4.0866025403784505, 13.806406460551015))
+    @test b
+    @test loc ≈ Point(-1.0000000000000058, 15.588457268119894)
+
+    b, loc = intersects(line_helper(5743.933982822018, 150.0, 5885.355339059327, -50.0),
+                        line_helper(5760.0, 100.0, 5760.0, 140.0))
+    @test b
+    @test loc ≈ Point(5760.0, 127.27922061357884)
 end
 
 @testset "Offsetintegers" begin