diff --git a/src/lines.jl b/src/lines.jl index c89ae8da..60a04b56 100644 --- a/src/lines.jl +++ b/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}