Calculates intersection length of a ray and an axis aligned bounding box (AABB)
rc: Center of AABB
dr: [length, width, height] of AABB
r0: starting point of ray
d0: direction of ray
intersect: Intersection length of ray and AABB
ri: Optional, ray enterence point
rf: Optional, ray exit point
IDL> aabb_intersect, [0,0,0], [1,1,1], [-1,0,0], [1,0,0], intersect, ri, rf IDL> print, intersect 1.0 IDL> print, ri -0.5 0.0 0.0 IDL> print, rf 0.5 0.0 0.0
PRO aabb_intersect, rc, dr, r0, d0, intersect, r_enter, r_exit ;+#aabb_intersect ;+Calculates intersection length of a ray and an axis aligned bounding box (AABB) ;+*** ;+##Input Arguments ;+ **rc**: Center of AABB ;+ ;+ **dr**: [length, width, height] of AABB ;+ ;+ **r0**: starting point of ray ;+ ;+ **d0**: direction of ray ;+ ;+##Output Arguments ;+ **intersect**: Intersection length of ray and AABB ;+ ;+ **ri**: Optional, ray enterence point ;+ ;+ **rf**: Optional, ray exit point ;+ ;+##Example Usage ;+```idl ;+IDL> aabb_intersect, [0,0,0], [1,1,1], [-1,0,0], [1,0,0], intersect, ri, rf ;+IDL> print, intersect ;+ 1.0 ;+IDL> print, ri ;+ -0.5 0.0 0.0 ;+IDL> print, rf ;+ 0.5 0.0 0.0 ;+``` v0 = d0/sqrt(total(d0*d0)) ;; There are 6 sides to a cube/grid side_inter = dblarr(6) ;; Intersection points of ray with planes defined by grid ipnts = dblarr(3,6) ;; Find whether ray intersects each side for i=0L,5 do begin j = fix(floor(i/2)) ind = where([0,1,2] ne j) if abs(v0[j]) gt 0 then begin ;; Intersection point with plane ipnts[*,i] = r0 + v0*( ( (rc[j] + ( (i mod 2)-0.5)*dr[j] ) - r0[j])/v0[j] ) ;; Check if point on plane is within grid side if abs(ipnts[ind[0],i] - rc[ind[0]]) le 0.5*dr[ind[0]] and $ abs(ipnts[ind[1],i] - rc[ind[1]]) le 0.5*dr[ind[1]] then side_inter[i]=1 endif endfor intersect = 0.0 r_enter = r0 r_exit = r0 w = where(side_inter ne 0,nw) if nw ge 2 then begin ;;Find two unique intersection points nunique = 0 for i=0,nw-2 do begin if total(ipnts[*,w[0]] eq ipnts[*,w[i+1]]) ne 3 then begin w = [w[0],w[i+1]] nunique = 2 break end end if nunique eq 2 then begin vi = ipnts[*,w[1]]-ipnts[*,w[0]] vi = vi/sqrt(total(vi*vi)) dot_prod = total(v0*vi) if dot_prod gt 0.0 then begin r_enter = ipnts[*,w[0]] r_exit = ipnts[*,w[1]] endif else begin r_enter = ipnts[*,w[1]] r_exit = ipnts[*,w[0]] endelse ;; Calculate intersection length intersect = sqrt(total((r_exit-r_enter)^2.0)) endif endif END