Skip to content

Commit b2468f1

Browse files
sbryngelsonsbryngelson
authored
Refactor/moller trumbore ray triangle (#1305)
1 parent 4e14d3d commit b2468f1

1 file changed

Lines changed: 90 additions & 72 deletions

File tree

‎src/common/m_model.fpp‎

Lines changed: 90 additions & 72 deletions
Original file line numberDiff line numberDiff line change
@@ -583,11 +583,17 @@ contains
583583

584584
end function f_model_is_inside
585585

586-
!> This procedure, given a cell center will determine if a point exists instide a surface
587-
!! @param ntrs Number of triangles in the model
588-
!! @param pid Patch ID od this model
586+
!> This procedure determines if a point is inside a surface using
587+
!! the generalized winding number (Jacobson et al., SIGGRAPH 2013).
588+
!! In 3D, sums the solid angle subtended by each triangle (Van
589+
!! Oosterom-Strackee formula). In 2D (p==0), sums the signed
590+
!! angle subtended by each boundary edge. Returns ~1.0 inside,
591+
!! ~0.0 outside. Unlike ray casting, this is robust to small
592+
!! triangles/edges and vertex winding order.
593+
!! @param ntrs Number of triangles in the model.
594+
!! @param pid Patch ID of this model.
589595
!! @param point Point to test.
590-
!! @return fraction The perfentage of candidate rays cast indicate that we are inside the model
596+
!! @return fraction Winding number (~1.0 inside, ~0.0 outside).
591597
function f_model_is_inside_flat(ntrs, pid, point) result(fraction)
592598

593599
$:GPU_ROUTINE(parallelism='[seq]')
@@ -597,52 +603,76 @@ contains
597603
real(wp), dimension(1:3), intent(in) :: point
598604

599605
real(wp) :: fraction
600-
type(t_ray) :: ray
601-
type(t_triangle) :: tri
602-
integer :: i, j, k, q, nInOrOut, nHits
603606

604-
! cast 26 rays from the point and count the number at leave the boundary
605-
nInOrOut = 0
606-
do i = -1, 1
607-
do j = -1, 1
608-
do k = -1, 1
609-
if (i /= 0 .or. j /= 0 .or. k /= 0) then
610-
! We cannot get inersections if the ray is exactly in line with triangle plane
611-
if (p == 0 .and. k == 0) cycle
612-
613-
! generate the ray
614-
ray%o = point
615-
ray%d(:) = [real(i, wp), real(j, wp), real(k, wp)]
616-
ray%d = ray%d/sqrt(real(abs(i) + abs(j) + abs(k), wp))
617-
618-
! count the number of intersections
619-
nHits = 0
620-
do q = 1, ntrs
621-
tri%v(:, :) = gpu_trs_v(:, :, q, pid)
622-
tri%n(1:3) = gpu_trs_n(1:3, q, pid)
623-
nHits = nHits + f_intersects_triangle(ray, tri)
624-
end do
625-
! if the ray intersected an odd number of times, we must be inside
626-
nInOrOut = nInOrOut + mod(nHits, 2)
627-
end if
628-
end do
629-
end do
630-
end do
607+
real(wp) :: r1(3), r2(3), r3(3)
608+
real(wp) :: r1_mag, r2_mag, r3_mag
609+
real(wp) :: numerator, denominator
610+
real(wp) :: d1(2), d2(2)
611+
integer :: q
612+
613+
fraction = 0.0_wp
631614

632615
if (p == 0) then
633-
! in 2D, we skipped 8 rays
634-
fraction = real(nInOrOut)/18._wp
616+
! 2D winding number: sum signed angles subtended by
617+
! each boundary edge at the query point.
618+
do q = 1, gpu_boundary_edge_count(pid)
619+
d1(1) = gpu_boundary_v(q, 1, 1, pid) - point(1)
620+
d1(2) = gpu_boundary_v(q, 1, 2, pid) - point(2)
621+
d2(1) = gpu_boundary_v(q, 2, 1, pid) - point(1)
622+
d2(2) = gpu_boundary_v(q, 2, 2, pid) - point(2)
623+
624+
! Signed angle = atan2(d1 x d2, d1 . d2)
625+
fraction = fraction + atan2( &
626+
d1(1)*d2(2) - d1(2)*d2(1), &
627+
d1(1)*d2(1) + d1(2)*d2(2))
628+
end do
629+
630+
! 2D winding number = total angle / (2*pi)
631+
fraction = fraction/(2.0_wp*acos(-1.0_wp))
635632
else
636-
fraction = real(nInOrOut)/26._wp
633+
! 3D winding number: sum solid angles via Van
634+
! Oosterom-Strackee formula.
635+
do q = 1, ntrs
636+
r1 = gpu_trs_v(1, :, q, pid) - point
637+
r2 = gpu_trs_v(2, :, q, pid) - point
638+
r3 = gpu_trs_v(3, :, q, pid) - point
639+
640+
r1_mag = sqrt(dot_product(r1, r1))
641+
r2_mag = sqrt(dot_product(r2, r2))
642+
r3_mag = sqrt(dot_product(r3, r3))
643+
644+
! Skip if query point is coincident with a vertex
645+
! (magnitudes are zero/subnormal).
646+
if (r1_mag*r2_mag*r3_mag < tiny(1.0_wp)) cycle
647+
648+
! tan(Omega/2) = numerator / denominator
649+
! numerator = scalar triple product r1 . (r2 x r3)
650+
numerator = r1(1)*(r2(2)*r3(3) - r2(3)*r3(2)) &
651+
+ r1(2)*(r2(3)*r3(1) - r2(1)*r3(3)) &
652+
+ r1(3)*(r2(1)*r3(2) - r2(2)*r3(1))
653+
654+
denominator = r1_mag*r2_mag*r3_mag &
655+
+ dot_product(r1, r2)*r3_mag &
656+
+ dot_product(r2, r3)*r1_mag &
657+
+ dot_product(r3, r1)*r2_mag
658+
659+
fraction = fraction + atan2(numerator, denominator)
660+
end do
661+
662+
! Each atan2 returns Omega/2 per triangle; divide
663+
! by 2*pi to get winding number = sum(Omega)/(4*pi).
664+
fraction = fraction/(2.0_wp*acos(-1.0_wp))
637665
end if
638666

639667
end function f_model_is_inside_flat
640668

641-
! From https://www.scratchapixel.com/lessons/3e-basic-rendering/ray-tracing-rendering-a-triangle/ray-triangle-intersection-geometric-solution.html
642-
!> This procedure checks if a ray intersects a triangle.
669+
!> This procedure checks if a ray intersects a triangle using the
670+
!! Moller-Trumbore algorithm (barycentric coordinates). Unlike the
671+
!! previous cross-product sign test, this is vertex winding-order
672+
!! independent.
643673
!! @param ray Ray.
644674
!! @param triangle Triangle.
645-
!! @return True if the ray intersects the triangle, false otherwise.
675+
!! @return 1 if the ray intersects the triangle, 0 otherwise.
646676
function f_intersects_triangle(ray, triangle) result(intersects)
647677

648678
$:GPU_ROUTINE(parallelism='[seq]')
@@ -652,46 +682,34 @@ contains
652682

653683
integer :: intersects
654684

655-
real(wp) :: N(3), P(3), C(3), edge(3), vp(3)
656-
real(wp) :: d, t, NdotRayDirection
685+
real(wp) :: edge1(3), edge2(3), h(3), s(3), q(3)
686+
real(wp) :: a, f, u, v, t
657687

658688
intersects = 0
659689

660-
N(1:3) = triangle%n(1:3)
661-
NdotRayDirection = dot_product(N(1:3), ray%d(1:3))
662-
if (abs(NdotRayDirection) < 0.0000001_wp) then
663-
return
664-
end if
690+
edge1 = triangle%v(2, :) - triangle%v(1, :)
691+
edge2 = triangle%v(3, :) - triangle%v(1, :)
692+
h = f_cross(ray%d, edge2)
693+
a = dot_product(edge1, h)
665694

666-
d = -sum(N(:)*triangle%v(1, :))
667-
t = -(sum(N(:)*ray%o(:)) + d)/NdotRayDirection
668-
if (t < 0) then
669-
return
670-
end if
695+
! Ray nearly parallel to triangle plane. In single precision
696+
! builds epsilon(1.0) ~ 1.2e-7, so use 10*epsilon as a floor.
697+
if (abs(a) < max(1e-7_wp, 10.0_wp*epsilon(1.0_wp))) return
671698

672-
P = ray%o + t*ray%d
673-
edge = triangle%v(2, :) - triangle%v(1, :)
674-
vp = P - triangle%v(1, :)
675-
C = f_cross(edge, vp)
676-
if (sum(N(:)*C(:)) < 0) then
677-
return
678-
end if
699+
f = 1.0_wp/a
700+
s = ray%o - triangle%v(1, :)
701+
u = f*dot_product(s, h)
679702

680-
edge = triangle%v(3, :) - triangle%v(2, :)
681-
vp = P - triangle%v(2, :)
682-
C = f_cross(edge, vp)
683-
if (sum(N(:)*C(:)) < 0) then
684-
return
685-
end if
703+
if (u < 0.0_wp .or. u > 1.0_wp) return
686704

687-
edge = triangle%v(1, :) - triangle%v(3, :)
688-
vp = P - triangle%v(3, :)
689-
C = f_cross(edge, vp)
690-
if (sum(N(:)*C(:)) < 0) then
691-
return
692-
end if
705+
q = f_cross(s, edge1)
706+
v = f*dot_product(ray%d, q)
707+
708+
if (v < 0.0_wp .or. u + v > 1.0_wp) return
709+
710+
t = f*dot_product(edge2, q)
693711

694-
intersects = 1
712+
if (t > 0.0_wp) intersects = 1
695713

696714
end function f_intersects_triangle
697715

0 commit comments

Comments
 (0)