Added test for ray_triangle_intersect #20

COMODO-research · Apr 1, 2024 · 87303d6 · 87303d6
1 parent ba8d2a0
commit 87303d6
Show file tree

Hide file tree

Showing 3 changed files with 89 additions and 9 deletions.
diff --git a/examples/demo_ray_triangle_intersect.jl b/examples/demo_ray_triangle_intersect.jl
@@ -25,7 +25,7 @@ hp1 = poly!(ax1,M,color=:white, shading = FastShading, transparency=true,strokec
 for x = range(-1.25,1.25,np)
     ray_origin = GeometryBasics.Point3{Float64}(x,0.25*sin(x*pi),1.25)
     ray_vector = Vec3{Float64}(0.0,0.0,-2)
-    P,indIntersect = ray_triangle_intersect(F,V,ray_origin,ray_vector;rayType = :ray, triSide = 1)
+    P,indIntersect = ray_triangle_intersect(F,V,ray_origin,ray_vector; rayType = :ray, triSide = 1)
     scatter!(ax1,ray_origin,markersize = markerSize,color=:blue)
     scatter!(ax1,ray_origin.+ray_vector,markersize = markerSize,color=:red)
     scatter!(ax1,P,markersize = markerSize,color=:green)
@@ -40,7 +40,7 @@ hp1 = poly!(ax1,M,color=:white, shading = FastShading, transparency=true,strokec
 for x = range(-1.25,1.25,np)
     ray_origin = GeometryBasics.Point3{Float64}(x,0.25*sin(x*pi),1.25)
     ray_vector = Vec3{Float64}(0.0,0.0,-2)
-    P,indIntersect = ray_triangle_intersect(F,V,ray_origin,ray_vector;rayType = :ray, triSide = 0)
+    P,indIntersect = ray_triangle_intersect(F,V,ray_origin,ray_vector; rayType = :ray, triSide = 0)
     scatter!(ax1,ray_origin,markersize = markerSize,color=:blue)
     scatter!(ax1,ray_origin.+ray_vector,markersize = markerSize,color=:red)
     scatter!(ax1,P,markersize = markerSize,color=:green)
@@ -55,7 +55,7 @@ hp1 = poly!(ax1,M,color=:white, shading = FastShading, transparency=true,strokec
 for x = range(-1.25,1.25,np)
     ray_origin = GeometryBasics.Point3{Float64}(x,0.25*sin(x*pi),1.25)
     ray_vector = Vec3{Float64}(0.0,0.0,-2)
-    P,indIntersect = ray_triangle_intersect(F,V,ray_origin,ray_vector;rayType = :ray, triSide = -1)
+    P,indIntersect = ray_triangle_intersect(F,V,ray_origin,ray_vector; rayType = :ray, triSide = -1)
     scatter!(ax1,ray_origin,markersize = markerSize,color=:blue)
     scatter!(ax1,ray_origin.+ray_vector,markersize = markerSize,color=:red)
     scatter!(ax1,P,markersize = markerSize,color=:green)
@@ -72,7 +72,7 @@ hp1 = poly!(ax1,M,color=:white, shading = FastShading, transparency=true,strokec
 for x = range(-1.25,1.25,np)
     ray_origin = GeometryBasics.Point3{Float64}(x,0.25*sin(x*pi),1.25)
     ray_vector = Vec3{Float64}(0.0,0.0,-2)
-    P,indIntersect = ray_triangle_intersect(F,V,ray_origin,ray_vector;rayType = :line, triSide = 1)
+    P,indIntersect = ray_triangle_intersect(F,V,ray_origin,ray_vector; rayType = :line, triSide = 1)
     scatter!(ax1,ray_origin,markersize = markerSize,color=:blue)
     scatter!(ax1,ray_origin.+ray_vector,markersize = markerSize,color=:red)
     scatter!(ax1,P,markersize = markerSize,color=:green)
@@ -87,7 +87,7 @@ hp1 = poly!(ax1,M,color=:white, shading = FastShading, transparency=true,strokec
 for x = range(-1.25,1.25,np)
     ray_origin = GeometryBasics.Point3{Float64}(x,0.25*sin(x*pi),1.25)
     ray_vector = Vec3{Float64}(0.0,0.0,-2)
-    P,indIntersect = ray_triangle_intersect(F,V,ray_origin,ray_vector;rayType = :line, triSide = 0)
+    P,indIntersect = ray_triangle_intersect(F,V,ray_origin,ray_vector; rayType = :line, triSide = 0)
     scatter!(ax1,ray_origin,markersize = markerSize,color=:blue)
     scatter!(ax1,ray_origin.+ray_vector,markersize = markerSize,color=:red)
     scatter!(ax1,P,markersize = markerSize,color=:green)
@@ -102,7 +102,7 @@ hp1 = poly!(ax1,M,color=:white, shading = FastShading, transparency=true,strokec
 for x = range(-1.25,1.25,np)
     ray_origin = GeometryBasics.Point3{Float64}(x,0.25*sin(x*pi),1.25)
     ray_vector = Vec3{Float64}(0.0,0.0,-2)
-    P,indIntersect = ray_triangle_intersect(F,V,ray_origin,ray_vector;rayType = :line, triSide = -1)
+    P,indIntersect = ray_triangle_intersect(F,V,ray_origin,ray_vector; rayType = :line, triSide = -1)
     scatter!(ax1,ray_origin,markersize = markerSize,color=:blue)
     scatter!(ax1,ray_origin.+ray_vector,markersize = markerSize,color=:red)
     scatter!(ax1,P,markersize = markerSize,color=:green)

diff --git a/src/functions.jl b/src/functions.jl
@@ -1145,8 +1145,6 @@ method both refines and smoothes the geometry through spline approximation.
 # References
 [Charles Loop, Smooth Subdivision Surfaces Based on Triangles M.S. Mathematics Thesis, University of Utah. 1987.](https://www.microsoft.com/en-us/research/wp-content/uploads/2016/02/thesis-10.pdf)
 [Jos Stam, Charles Loop, Quad/Triangle Subdivision, doi: 10.1111/1467-8659.t01-2-00647](https://doi.org/10.1111/1467-8659.t01-2-00647)
-
-
 """
 function subtri(F,V,n; method = :linear)
 
@@ -1226,6 +1224,7 @@ function subtri(F,V,n; method = :linear)
     end
 end
 
+
 function subquad(F,V,n; method=:linear)
     if iszero(n)
         return F,V
@@ -1289,6 +1288,7 @@ function subquad(F,V,n; method=:linear)
     end
 end
 
+
 # Create geodesic dome
 function geosphere(n,r)
     M = platonicsolid(4,r)
@@ -1305,6 +1305,7 @@ function geosphere(n,r)
     return F,V
 end
 
+
 function hexbox(boxDim::Vector{T},boxEl::Vector{Int64}) where T <: Real
     boxNod = boxEl.+1 # Number of nodes in each direction
     numElements = prod(boxEl) # Total number of elements
@@ -1378,6 +1379,7 @@ function hexbox(boxDim::Vector{T},boxEl::Vector{Int64}) where T <: Real
     return E,V,F,Fb,CFb_type
 end
 
+
 function con_face_edge(F,E_uni=nothing,indReverse=nothing)
     if isnothing(E_uni) | isnothing(indReverse)
         E = meshedges(F)
@@ -1386,6 +1388,7 @@ function con_face_edge(F,E_uni=nothing,indReverse=nothing)
     return [Vector{Int64}(a) for a ∈ eachrow(reshape(indReverse,length(F),length(F[1])))] # [indReverse[[1,2,3].+ (i-1)*3] for i ∈ eachindex(F)]
 end
 
+
 function con_edge_face(F,E_uni=nothing,indReverse=nothing)
     if isnothing(E_uni) || isnothing(indReverse)
         E = meshedges(F)
@@ -1402,6 +1405,7 @@ function con_edge_face(F,E_uni=nothing,indReverse=nothing)
     return con_E2F
 end
 
+
 function con_face_face(F,E_uni=nothing,indReverse=nothing,con_E2F=nothing,con_F2E=nothing)
     if length(F)>1 # More than one face so compute connectivity
         if isnothing(E_uni)| isnothing(indReverse)
@@ -1428,6 +1432,7 @@ function con_face_face(F,E_uni=nothing,indReverse=nothing,con_E2F=nothing,con_F2
     end
 end
 
+
 function con_face_face_v(F,V=nothing,con_V2F=nothing)
     if length(F)>1 # More than one face so compute connectivity
         if isnothing(con_V2F) 
@@ -1447,6 +1452,7 @@ function con_face_face_v(F,V=nothing,con_V2F=nothing)
     end
 end
 
+
 function con_vertex_simplex(F,V=nothing)
     if isnothing(V)
         n = maximum(reduce(vcat,F))
@@ -1462,14 +1468,17 @@ function con_vertex_simplex(F,V=nothing)
     return con_V2F
 end
 
+
 function con_vertex_face(F,V=nothing)
     return con_vertex_simplex(F,V)
 end
 
+
 function con_vertex_edge(E,V=nothing)
     return con_vertex_simplex(E,V)
 end
 
+
 function con_edge_edge(E_uni,con_V2E=nothing)
     if isnothing(con_V2E)
         con_V2E = con_vertex_edge(E_uni) 
@@ -1485,6 +1494,7 @@ function con_edge_edge(E_uni,con_V2E=nothing)
     return con_E2E
 end
 
+
 function con_vertex_vertex_f(F,V=nothing,con_V2F=nothing)
     if isnothing(V)
         n = maximum(reduce(vcat,F))
@@ -1510,6 +1520,7 @@ function con_vertex_vertex_f(F,V=nothing,con_V2F=nothing)
 
 end
 
+
 function con_vertex_vertex(E,V=nothing,con_V2E=nothing)
     if isnothing(V)
         n = maximum(reduce(vcat,E))
@@ -2378,6 +2389,7 @@ function ray_triangle_intersect(F::Vector{TriangleFace{Int64}},V,ray_origin,ray_
     end
     return P,indIntersect
 end
+
 function ray_triangle_intersect(f::TriangleFace{Int64},V,ray_origin,ray_vector; rayType = :ray, triSide = 1, tolEps = eps(Float64))
 
     # Edge vectors

diff --git a/test/runtests.jl b/test/runtests.jl
@@ -2061,7 +2061,75 @@ end
 end
 
 
-@testset "separate vertices" begin
+# @testset "distseedpoints" verbose = true begin
+
+# end
+
+@testset "ray_triangle_intersect" verbose = true begin
+    eps_level = 0.001
+
+    # Single cube
+    r = 2 * sqrt(3) / 2
+    M = cube(r)
+    F = faces(M)
+    V = coordinates(M)
+    F = quad2tri(F,V,convert_method = :forward)
+    @testset "ray" begin 
+        ray_origin = GeometryBasics.Point3{Float64}(0.25,0.0,1.5) # Slight off so we hit one triangle, not two at the edge
+        ray_vector = Vec3{Float64}(0.0,0.0,-1)
+
+        P,indIntersect = ray_triangle_intersect(F,V,ray_origin,ray_vector; rayType = :ray, triSide = 1)
+        @test isapprox(P,Point3{Float64}[[0.25, 0.0, 1.0]],atol=eps_level)
+        @test isa(indIntersect,Vector{Int64}) # indIntersect == 3
+
+        P,indIntersect = ray_triangle_intersect(F,V,ray_origin,ray_vector; rayType = :ray, triSide = 0)
+        @test isapprox(P,Point3{Float64}[[0.25, 0.0, -1.0], [0.25, 0.0, 1.0]],atol=eps_level)
+        @test isa(indIntersect,Vector{Int64})
+
+        P,indIntersect = ray_triangle_intersect(F,V,ray_origin,ray_vector; rayType = :ray, triSide = -1)
+        @test isapprox(P,Point3{Float64}[[0.25, 0.0, -1.0]],atol=eps_level)
+        @test isa(indIntersect,Vector{Int64})
+
+        ray_origin = GeometryBasics.Point3{Float64}(0.0,0.0,1.5) # At centre so hits an edge between two triangles
+        ray_vector = Vec3{Float64}(0.0,0.0,-1)
+        P,indIntersect = ray_triangle_intersect(F,V,ray_origin,ray_vector; rayType = :ray, triSide = 0, tolEps = 1e-3)
+        @test isapprox(P,Point3{Float64}[[0.0, 0.0, -1.0], [0.0, 0.0, -1.0], 
+                                         [0.0, 0.0,  1.0], [0.0, 0.0,  1.0]],atol=eps_level)
+        @test isa(indIntersect,Vector{Int64})
+    end
+
+    @testset "line type" begin 
+        ray_origin = GeometryBasics.Point3{Float64}(0.25,0.0,1.5) # Slight off so we hit one triangle, not two at the edge
+        ray_vector = Vec3{Float64}(0.0,0.0,-1) # Shorst so only one hit
+
+        P,indIntersect = ray_triangle_intersect(F,V,ray_origin,ray_vector; rayType = :line, triSide = 1)
+        @test isapprox(P,Point3{Float64}[[0.25, 0.0, 1.0]],atol=eps_level)
+        @test isa(indIntersect,Vector{Int64}) # indIntersect == 3
+
+        ray_vector = Vec3{Float64}(0.0,0.0,-3) # Long so two hits potentially
+
+        P,indIntersect = ray_triangle_intersect(F,V,ray_origin,ray_vector; rayType = :line, triSide = 1)
+        @test isapprox(P,Point3{Float64}[[0.25, 0.0, 1.0]],atol=eps_level)
+        @test isa(indIntersect,Vector{Int64}) # indIntersect == 3
+
+        P,indIntersect = ray_triangle_intersect(F,V,ray_origin,ray_vector; rayType = :line, triSide = 0)
+        @test isapprox(P,Point3{Float64}[[0.25, 0.0, -1.0], [0.25, 0.0, 1.0]],atol=eps_level)
+        @test isa(indIntersect,Vector{Int64})
+
+        P,indIntersect = ray_triangle_intersect(F,V,ray_origin,ray_vector; rayType = :line, triSide = -1)
+        @test isapprox(P,Point3{Float64}[[0.25, 0.0, -1.0]],atol=eps_level)
+        @test isa(indIntersect,Vector{Int64})
+
+        ray_origin = GeometryBasics.Point3{Float64}(0.0,0.0,1.5) # At centre so hits an edge between two triangles
+        P,indIntersect = ray_triangle_intersect(F,V,ray_origin,ray_vector; rayType = :line, triSide = 0, tolEps = 1e-3)
+        @test isapprox(P,Point3{Float64}[[0.0, 0.0, -1.0], [0.0, 0.0, -1.0], 
+                                         [0.0, 0.0,  1.0], [0.0, 0.0,  1.0]],atol=eps_level)
+        @test isa(indIntersect,Vector{Int64})
+    end
+end
+
+
+@testset "separate_vertices" begin
 
     eps_level = 0.001
     r = 2 * sqrt(3) / 2