intersect_RayTriangle Function

public function intersect_RayTriangle(line, triangle, intersect_p) result(isIntersect)

Function computes intersection of ray with triangle

http://geomalgorithms.com/a06-_intersect-2.html intersect_RayTriangle(): intersect a ray with a 3D triangle Input: a ray R, and a triangle T Output: *I = intersection point (when it exists) Return: -1 = triangle is degenerate (a segment or point) 0 = disjoint (no intersect) 1 = intersect in unique point I1 2 = are in the same planeint todo: when line lies in triangle, need to treat properly

Arguments

Type IntentOptional Attributes Name
type(tem_line_type), intent(in) :: line

line segment to check for interection
type(tem_triangle_type), intent(in) :: triangle

cube to check intersection of line
type(tem_point_type), intent(out), optional :: intersect_p

intersection point if there is intersection

Return Value logical


Calls

proc~~intersect_raytriangle~~CallsGraph proc~intersect_raytriangle intersect_RayTriangle proc~cross_product3d cross_product3D proc~intersect_raytriangle->proc~cross_product3d

Source Code

  function intersect_RayTriangle( line, triangle, intersect_p ) result(isIntersect)
    ! ---------------------------------------------------------------------------!
    !> line segment to check for interection
    type(tem_line_type), intent(in) :: line
    !> cube to check intersection of line
    type(tem_triangle_type), intent(in) :: triangle
    !> intersection point if there is intersection
    type( tem_point_type), optional, intent(out) :: intersect_p
    logical :: isIntersect
    ! ---------------------------------------------------------------------------!
    real(kind=rk) :: u(3), v(3), n(3)  ! triangle vectors and normal vector
    real(kind=rk) :: dir(3), w0(3), w(3)    ! ray vectors
    real(kind=rk) :: r, a, b       ! params to calc ray-plane intersect
    real(kind=rk) :: uu, uv, vv, wu, wv, D
    real(kind=rk) :: s, t
    real(kind=rk) :: temp_p(3)
    ! ---------------------------------------------------------------------------!
    isIntersect = .false. ! set not intersect as default

    ! get triangle edge vectors u & v, and plane normal n
    u(:) = triangle%nodes(:,2) - triangle%nodes(:,1)
    v(:) = triangle%nodes(:,3) - triangle%nodes(:,1)
    n = cross_product3D( u, v )

    ! triangle is degenerate
    ! do not deal with this case
    if (all(n .feq. 0._rk)) return

    dir = line%vec   ! ray direction vector
    w0 = line%origin - triangle%nodes(:,1)
    a = -dot_product( n, w0);
    b =  dot_product( n, dir);

    ! if ray parallel to triangle plane, treat it as no intersecting
    if (abs(b) < eps .and. abs(a) > tiny(a)) return
!    {   if (a == 0.0_rk)         ! ray lies in triangle plane
!          return 2;
!        else return 0;             ! ray disjoint from plane
!    }

    if (abs(b) < eps) then
      ! origin is very close to triangle plane, but parallel.
      ! Just use the origin itself as point to check.
      r = 0.0_rk
    else
      ! get intersect point of ray with triangle plane
      r = a / b
      if (r < 0._rk ) then     ! ray goes away from triangle
        return
      endif
    end if
    ! for a segment, also test if (r > 1.0) => no intersect

    ! intersect point of ray and plane
    temp_p = line%origin + r * dir

    ! is I inside T?
    uu = dot_product(u,u)
    uv = dot_product(u,v)
    vv = dot_product(v,v)
    w = temp_p - triangle%nodes(:,1)
    wu = dot_product(w,u)
    wv = dot_product(w,v)
    D = uv * uv - uu * vv

    ! get and test parametric coords
    s = (uv * wv - vv * wu) / D
    ! point is outside triangle
    if (s < 0._rk-eps .or. s > 1._rk+eps) then
      return
    endif

    t = (uv * wu - uu * wv) / D
    ! point is outside triangle
    if (t < 0._rk-eps .or. (s + t) > 1._rk+eps) then
      return
    endif

    isIntersect = .true. ! point is inside triangle
    if( present( intersect_p ) ) then
      intersect_p%coord = temp_p
    endif

  end function intersect_RayTriangle