This routine checks for triangle box overlap
this routine is conversion of c-code tribox3.c triBoxOverlap function. use separating axis theorem to test overlap between triangle and box need to test for overlap in these directions: 1) the {x,y,z}-directions (actually, since we use the AABB of the triangle we do not even need to test these) 2) normal of the triangle 3) separating axis test. crossproduct(edge from tri, {x,y,z}-directin) this gives 3x3=9 more 7tests This code is available at: http://fileadmin.cs.lth.se/cs/Personal/Tomas_Akenine-Moller/code/tribox3.txt
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real(kind=rk), | intent(in) | :: | boxCenter(3) |
box center |
||
real(kind=rk), | intent(in) | :: | boxHalfwidth(3) |
halfwidth of the box |
||
real(kind=rk), | intent(in) | :: | triNodes(3,3) |
nodes of the triangle 1st index denote x,y,z coordinates and 2nd index denote nodes |
function triBoxOverlap_loc( boxCenter, boxHalfwidth, triNodes ) &
& result(overlaps)
!--------------------------------------------------------------------------!
!> box center
real(kind=rk), intent(in) :: boxCenter(3)
!> halfwidth of the box
real(kind=rk), intent(in) :: boxHalfwidth(3)
!> nodes of the triangle
!! 1st index denote x,y,z coordinates and
!! 2nd index denote nodes
real(kind=rk), intent(in) :: triNodes(3,3)
logical :: overlaps
!--------------------------------------------------------------------------!
real(kind=rk) :: nodes(3,3)
real(kind=rk) :: dirVec1(3), dirVec2(3), dirVec3(3)
real(kind=rk) :: mmin, mmax
real(kind=rk) :: normal(3), edge1(3), edge2(3), edge3(3)
!--------------------------------------------------------------------------!
overlaps = .false.
dirVec1 = (/ 1.0_rk, 0.0_rk, 0.0_rk /)
dirVec2 = (/ 0.0_rk, 1.0_rk, 0.0_rk /)
dirVec3 = (/ 0.0_rk, 0.0_rk, 1.0_rk /)
! move everything so that the box center is in (0,0,0)
nodes(:,1) = triNodes(:,1) - boxCenter
nodes(:,2) = triNodes(:,2) - boxCenter
nodes(:,3) = triNodes(:,3) - boxCenter
! compute triangle edges
edge1 = nodes(:,2) - nodes(:,1)
edge2 = nodes(:,3) - nodes(:,2)
edge3 = nodes(:,1) - nodes(:,3)
! separating axis test
if(axistest(dirVec1, edge1, nodes, boxhalfwidth)) return
if(axistest(dirVec1, edge2, nodes, boxhalfwidth)) return
if(axistest(dirVec1, edge3, nodes, boxhalfwidth)) return
if(axistest(dirVec2, edge1, nodes, boxhalfwidth)) return
if(axistest(dirVec2, edge2, nodes, boxhalfwidth)) return
if(axistest(dirVec2, edge3, nodes, boxhalfwidth)) return
if(axistest(dirVec3, edge1, nodes, boxhalfwidth)) return
if(axistest(dirVec3, edge2, nodes, boxhalfwidth)) return
if(axistest(dirVec3, edge3, nodes, boxhalfwidth)) return
! Bullet 1:
! first test overlap in the {x,y,z}-directions
! find min, max of the triangle each direction, and test for overlap in
! that direction -- this is equivalent to testing a minimal AABB around
! the triangle against the AABB
! test in X-direction
mmin = min( nodes(1,1), nodes(1,2), nodes(1,3) )
mmax = max( nodes(1,1), nodes(1,2), nodes(1,3) )
if ( mmin > boxhalfwidth(1) .or. mmax < -boxhalfwidth(1) ) return
! test in Y-direction
mmin = min( nodes(2,1), nodes(2,2), nodes(2,3) )
mmax = max( nodes(2,1), nodes(2,2), nodes(2,3) )
if ( mmin > boxhalfwidth(2) .or. mmax < -boxhalfwidth(2) ) return
! test in Z-direction
mmin = min( nodes(3,1), nodes(3,2), nodes(3,3) )
mmax = max( nodes(3,1), nodes(3,2), nodes(3,3) )
if ( mmin > boxhalfwidth(3) .or. mmax < -boxhalfwidth(3) ) return
! Bullet 2:
! test if the box intersects the plane of the triangle
! compute plane equation of triangle: normal*x+d=0
normal = cross_product3D(edge1,edge2)
if( .not. planeBoxOverlap(normal, nodes(:,1), boxhalfwidth)) return
overlaps = .true.
end function triBoxOverlap_loc