Fix for intersection point calculation (#31)

* Fix for interaction point calc. * Better intersect expression; add new test case. * Rename intermediate vars * Add comments to intersection Co-authored-by: Júlio Hoffimann <[email protected]> * Add ref link. * Update src/polytopes/segment.jl Co-authored-by: Júlio Hoffimann <[email protected]>
JuliaGeometry · Jan 31, 2021 · 5025edf · 5025edf
1 parent 3f3fef1
commit 5025edf
Show file tree

Hide file tree

Showing 2 changed files with 13 additions and 11 deletions.
diff --git a/src/polytopes/segment.jl b/src/polytopes/segment.jl
@@ -138,18 +138,16 @@ function intersectpoint(s1::Segment{Dim,T}, s2::Segment{Dim,T}) where {Dim,T}
   x1, x2 = vertices(s1)
   y1, y2 = vertices(s2)
 
-  # solve Ax = b
-  A = (y2 - y1) - (x2 - x1)
-  b = (x1 - y1)
-  s = zero(T)
-  for i in 1:length(A)
-    if !iszero(A[i])
-      s = b[i] / A[i]
-      break
-    end
-  end
+  v1 = x1 - x2
+  v2 = y1 - y2
+  v12 = x1 - y1
+
+  # the intersection point lies in between x1 and x2 at a fraction s
+  # (https://en.wikipedia.org/wiki/Line-line_intersection#Formulas)
+  s = (v12[1] * v2[2] - v12[2] * v2[1]) /
+      (v1[1] * v2[2] - v1[2] * v2[1])
 
-  x1 + s*(x2 - x1)
+  x1 - s*v1
 end
 
 # intersection of two line segments assuming that they are colinear

diff --git a/test/polytopes.jl b/test/polytopes.jl
@@ -82,6 +82,10 @@
     s1 = Segment(P2(2,1), P2(1,2))
     s2 = Segment(P2(1,0), P2(1,1))
     @test s1 ∩ s2 === s2 ∩ s1 === nothing
+
+    s1 = Segment(P2(1.5,1.5), P2(3.0,1.5))
+    s2 = Segment(P2(3.0,1.0), P2(2.0,2.0))
+    @test s1 ∩ s2 == s2 ∩ s1 == P2(2.5,1.5)
   end
 
   @testset "Triangles" begin