polyhedron V1.0

src/PointInPoly.jl

@@ -1,5 +1,7 @@
 module PointInPoly
 export pinpoly 
 
-
+include("polygon.jl")
+include("polyhedron.jl")
+include("utils.jl")
 end

src/polyhedron.jl

@@ -2,23 +2,53 @@ const INF = 1.0e9
 const EPS = 1.0e-9
 const BIAS = EPS * [ 0.09816024370339371,  0.6196549438024468,  0.9755418207471769]
 
+@inline function cross_product(a::NTuple{3,Float64}, b::NTuple{3,Float64})
+    return Float64[a[2]*b[3]-b[2]*a[3], a[3]*b[1]-a[1]*b[3], a[1]*b[2]-b[1]*a[2]] 
+end
+
+function between(p::NTuple{3,Float64}, point1::NTuple{3,Float64}, point2::NTuple{3,Float64})
+    ans = true
+    for i = 1:3
+        if !(point1[i] < p[i] < point2[i])
+            ans = false
+            break
+        end
+    end
+    return ans
+end
+
 """
-Determine whether a point lies inside a 3D polyhedron. Return 1 if inside, 0 if outside, -1 if exactly on the edge. 
+Determine whether a point lies inside a 3D polyhedron. Return 1 if inside, 0 if outside or exactly on the boundary. 
 
-`pinpoly(nodes::Vector{NTuple{3, Real}}, faces::Vector{NTuple{3,Int}}, point::NTuple{3, Real})`
+`pinpoly(node_x::NTuple{N,Float64}, node_y::NTuple{N,Float64}, node_z::NTuple{N,Float64}, faces::NTuple{N,NTuple{3,Int}}, point::NTuple{3, Float64}) where N`
 
-`nodes`: Vector of the node positions of the polyhedron.
+`nodes_x/y/z`: Vectors of the node positions of the polyhedron.
 
-`faces`: Vector of the triangular faces on the boundary.
+`faces`: Vector of the node numbers of the triangular faces on the boundary.
  
 `point`: Position of the point.
 
 This algorithm was proposed by W. Randolph Franklin: https://wrf.ecse.rpi.edu//Research/Short_Notes/pnpoly.html.
 """  
-function pinpoly(nodes::Vector{Vector{Float64}}, faces::Vector{NTuple{3,Int}}, point::NTuple{3, Float64})
+function pinpoly(node_x::NTuple{N,Float64}, node_y::NTuple{N,Float64}, node_z::NTuple{N,Float64}, faces::NTuple{N,NTuple{3,Int}}, point::NTuple{3, Float64}) where N
+    # Check if the point is out of the bounding box.
+    if !between(point, (minimum(node_x), minimum(node_y), minimum(node_z)), (maximum(node_x), maximum(node_y), maximum(node_z)))
+        return 0
+    end
+
     c = 0
-    for face in faces
-        cross = pcrossface(nodes[face[1]], nodes[face[2]], nodes[face[3]], point)
+    for i in eachindex(faces)
+        face = faces[i]
+        cross = pcrossface(
+            (node_x[face[1]], node_y[face[1]], node_z[face[1]]), 
+            (node_x[face[2]], node_y[face[2]], node_z[face[2]]), 
+            (node_x[face[3]], node_y[face[3]], node_z[face[3]]), 
+            point
+            )
+        # println(cross)
+        # println(((node_x[face[1]], node_y[face[1]], node_z[face[1]]), 
+        # (node_x[face[2]], node_y[face[2]], node_z[face[2]]), 
+        # (node_x[face[3]], node_y[face[3]], node_z[face[3]])))
         if cross == -1
             return -1
         end
@@ -29,23 +59,29 @@ function pinpoly(nodes::Vector{Vector{Float64}}, faces::Vector{NTuple{3,Int}}, p
     return c
 end
 
-function pcrossface(node1::Vector{Float64}, node2::Vector{Float64}, node3::Vector{Float64}, point::NTuple{3, Float64})
+function pcrossface(nodeA::NTuple{3, Float64}, nodeB::NTuple{3, Float64}, nodeC::NTuple{3, Float64}, point::NTuple{3, Float64})
+
+    node1 = collect(nodeA)
+    node2 = collect(nodeB)
+    node3 = collect(nodeC)
+    O = collect(point)
+
     e1 = node2 - node1
     e2 = node3 - node1
 
-    O = collect(point)
     D = Float64[INF, 0., 0.]
 
-    P = cross_product(D, e2)
+    P = cross_product(Tuple(D), Tuple(e2))
 
     det = e1' * P
-    # If the determinant is near zero, the ray lies in the plane of the triangle.
+    # If the determinant is near zero, the ray lies parallel to the plane of the triangle.
     if abs(det) == 0.0
-        # 已知射线在三角形平面内。
+        # 已知射线与三角形共面，即参数t无穷大。
         # Determine if the point is in the triangle.
-        if pintriangle(node1, node2, node3, point) != 0
-            return -1
-        end
+        # if pintriangle(node1, node2, node3, point) != 0
+        #     return -1
+        # end
+
         return 0
     end
     inv_det = 1.0 / det
@@ -56,40 +92,30 @@ function pcrossface(node1::Vector{Float64}, node2::Vector{Float64}, node3::Vecto
         return 0
     end
 
-    Q = cross_product(T, e1)
+    Q = cross_product(Tuple(T), Tuple(e1))
     v = D' * Q * inv_det
     if v < 0.0 || u + v > 1.0
         return 0
     end
 
     # 已知射线与三角形平面的交点不在三角形外。
     t = e2' * Q * inv_det
-    if abs(t) == 0.0
+    #注意射线是单向的，排除反向相交情况。
+    if t < 0.
+        return 0
+    end
+    if t == 0.
         return -1
     end
 
     # 已知射线与三角形平面的交点不在三角形外，且射线起点不在三角形平面内。
     if u == 0.0 || v == 0.0 || u+v == 1.0
-        return pcrossface(node1, node2, node3, point + BIAS)
+        return pcrossface(nodeA, nodeB, nodeC, (point[1] + BIAS[1], point[2] + BIAS[2], point[3] + BIAS[3]))
     end
 
     # 已知射线与三角形平面的交点在三角形内（不包括边上），且射线起点不在三角形平面内。
     return 1     
 end
 
-function pintriangle(A::Vector{Float64}, B::Vector{Float64}, C::Vector{Float64}, point::NTuple{3, Float64})
-    P = collect(point)
-    PA = P - A
-    PB = P - B
-    PC = P - C
 
-    t1 = sign(cross_product(PA, PB))
-    t2 = sign(cross_product(PB, PC))
-    t3 = sign(cross_product(PC, PA))
 
-    if t1==t2 && t2==t3
-        return 1
-    else
-        return 0
-    end
-end

src/utils.jl

@@ -7,3 +7,14 @@ function cross_product(a::Vector{T} where T <: Number, b::Vector{T} where T <: N
         error("undef dim")
     end
 end
+
+function betweeneq(p::NTuple{3,Float64}, point1::NTuple{3,Float64}, point2::NTuple{3,Float64})
+    ans = true
+    for i = 1:3
+        if !(point1[i] <= p[i] <= point2[i])
+            ans = false
+            break
+        end
+    end
+    return ans
+end

test/example.jl

@@ -0,0 +1,10 @@
+include("../src/polyhedron.jl")
+
+node_x=Tuple(Float64[0,1,0,0])
+node_y=Tuple(Float64[0,0,1,0])
+node_z=Tuple(Float64[0,0.001,0,1])
+faces = ((1,2,3), (1,2,4), (2,3,4), (1,3,4))
+points = ((0.01, 0.01, 0.01), (0.01,0.01,0.), (0.01, 0.01, -0.01))
+for point in points
+    println(point," -> ", pinpoly(node_x, node_y, node_z, faces, point))
+end