Merge pull request #3630 from JuliaReach/schillic/3629

#3629 -  Fix `intersection` of `LineSegments`

src/Sets/LineSegment/intersection.jl

@@ -34,24 +34,34 @@ function intersection(LS1::LineSegment, LS2::LineSegment)
     m = intersection(L1, L2)
     N = promote_type(eltype(LS1), eltype(LS2))
     if m == L1
-        # determine which segment is in both
-        p1 = max(min(LS1.p[1], LS1.q[1]), min(LS2.p[1], LS2.q[1]))
-        p2 = max(min(LS1.p[2], LS1.q[2]), min(LS2.p[2], LS2.q[2]))
-        q1 = min(max(LS1.p[1], LS1.q[1]), max(LS2.p[1], LS2.q[1]))
-        q2 = min(max(LS1.p[2], LS1.q[2]), max(LS2.p[2], LS2.q[2]))
-        if _isapprox(p1, q1) && _isapprox(p2, q2)
-            return Singleton([p1, p2])  # edges have a point in common
+        # line segments are on the same line
 
-        elseif _leq(p1, q1) && _leq(p2, q2)
-            return LineSegment([p1, p2], [q1, q2])
+        # check that the line segments intersect
+        @inbounds begin
+            # box approximation
+            l, r = extrema((LS1.p[1], LS1.q[1]))
+            b, t = extrema((LS1.p[2], LS1.q[2]))
+            # check that the other line segment has at least one point inside the box
+            if !(l <= LS2.p[1] <= r && b <= LS2.p[2] <= t) &&
+               !(l <= LS2.q[1] <= r && b <= LS2.q[2] <= t)
+                return EmptySet{N}(2)
+            end
+        end
 
+        # find the middle two points of the four end points
+        points = [LS1.p, LS1.q, LS2.p, LS2.q]
+        sorted_points = sort(points; by=p -> (p[1], p[2]))
+        @inbounds begin
+            mid_point1 = sorted_points[2]
+            mid_point2 = sorted_points[3]
+        end
+        if _isapprox(mid_point1, mid_point2)
+            return Singleton(mid_point1)  # only one point in common
         else
-            return EmptySet{N}(2)  # no intersection
+            return LineSegment(mid_point1, mid_point2)
         end
-
     elseif m isa Singleton && m.element ∈ LS1 && m.element ∈ LS2
         return m  # the intersection point between the lines is in the segments
-
     else
         return EmptySet{N}(2)  # no intersection
     end

test/Sets/LineSegment.jl

@@ -123,34 +123,38 @@ for N in [Float64, Rational{Int}, Float32]
     @test !is_intersection_empty(l1, l1_copy) && !intersection_empty && point ∈ l1
 
     # intersection
-    s1 = LineSegment(N[-5.0, -5.0], N[5.0, 5.0])
+    s1 = LineSegment(N[-5, -5], N[5, 5])
     # collinear shifted down
-    s2 = LineSegment(N[-6.0, -6.0], N[4.0, 4.0])
-    @test intersection(s1, s2) isa LineSegment
-    @test isapprox(intersection(s1, s2).p, N[-5, -5])
-    @test isapprox(intersection(s1, s2).q, N[4, 4])
+    s2 = LineSegment(N[-6, -6], N[4, 4])
+    cap = intersection(s1, s2)
+    @test cap isa LineSegment
+    @test ispermutation([cap.p, cap.q], [N[-5, -5], N[4, 4]])
     # parallel, not intersect
-    s3 = LineSegment(N[-5.0, -4.0], N[4.0, 5.0])
-    @test intersection(s1, s3) isa EmptySet
+    s3 = LineSegment(N[-5, -4], N[4, 5])
+    @test intersection(s1, s3) == EmptySet{N}(2)
     # intersect outside of segment
-    s4 = LineSegment(N[0.0, 10.0], N[6.0, 5.0])
-    @test intersection(s1, s4) isa EmptySet
+    s4 = LineSegment(N[0, 10], N[6, 5])
+    @test intersection(s1, s4) == EmptySet{N}(2)
     # intersect in segment
-    s5 = LineSegment(N[5.0, -5.0], N[-5.0, 5.0])
+    s5 = LineSegment(N[5, -5], N[-5, 5])
     @test intersection(s1, s5) isa Singleton
     @test isapprox(intersection(s1, s5).element, N[0, 0])
     # parallel, no points in common
-    s6 = LineSegment(N[10.0, 10.0], N[11.0, 11.0])
-    @test intersection(s1, s6) isa EmptySet
+    s6 = LineSegment(N[10, 10], N[11, 11])
+    @test intersection(s1, s6) == EmptySet{N}(2)
     # parallel one point in common
-    s7 = LineSegment(N[5.0, 5.0], N[6.0, 6.0])
+    s7 = LineSegment(N[5, 5], N[6, 6])
     @test intersection(s1, s7) isa Singleton
     @test isapprox(intersection(s1, s7).element, N[5, 5])
-    s8 = LineSegment(N[0.0, 0.0], N[1.0, 0.0])
-    s9 = LineSegment(N[0.0, 0.0], N[2.0, 0.0])
-    @test intersection(s9, s8) isa LineSegment
-    @test isapprox(intersection(s9, s8).p, N[0, 0])
-    @test isapprox(intersection(s9, s8).q, N[1, 0])
+    s8 = LineSegment(N[0, 0], N[1, 0])
+    s9 = LineSegment(N[0, 0], N[2, 0])
+    cap = intersection(s8, s9)
+    @test cap isa LineSegment
+    @test ispermutation([cap.p, cap.q], [N[0, 0], N[1, 0]])
+    # intersect in segment, different orientation
+    s10 = LineSegment(N[-1, 2], N[2, -1])
+    s11 = LineSegment(N[0, 1], N[1, 0])
+    @test intersection(s10, s11) == s11
 
     # subset
     l = LineSegment(N[1, 1], N[2, 2])