text: # copied from https://gist.github.com/dendenxu/ee5008acb5607195582e7983a384e644#file-moller_trumbore-py-L8
  
  
  import torch
  
  
  def normalize(x: torch.Tensor, eps: float = 1e-8) -> torch.Tensor:
      return x / (x.norm(dim=-1, keepdim=True) + eps)
  
  
  def inside_mesh(ray_o: torch.Tensor, verts: torch.Tensor, faces: torch.Tensor) -> torch.Tensor:
      """
      Check whether a bunch of points is inside the mesh defined by verts and faces
      effectively calculating their occupancy values
  
      Parameters
      ----------
      ray_o : torch.Tensor(float), (n_rays, 3)
      verts : torch.Tensor(float), (n_verts, 3)
      faces : torch.Tensor(long), (n_faces, 3)
      """
      ray_d = torch.rand_like(ray_o)  # (n_rays, 3)
      ray_d = normalize(ray_d)  # (n_rays, 3)
      u, v, t = moller_trumbore(ray_o, ray_d, verts[faces])  # (n_rays, n_faces, 3)
      intersection = ((t >= 0.0) * (u >= 0.0) * (v >= 0.0) * ((u + v) <= 1.0)).bool()  # (n_rays, n_faces)
      inside = (intersection.count_nonzero(dim=-1) % 2).bool()  # if mod 2 is 0, outside, inside is false
      return inside
  
  
  def moller_trumbore(ray_o: torch.Tensor, ray_d: torch.Tensor, tris: torch.Tensor, eps=1e-8):
      """
      The Moller Trumbore algorithm for fast ray triangle intersection
      Naive batch implementation (m rays and n triangles at the same time)
  
      Parameters
      ----------
      ray_o : torch.Tensor, (n_rays, 3)
      ray_d : torch.Tensor, (n_rays, 3)
      faces : torch.Tensor, (n_faces, 3, 3)
      """
      E1 = tris[:, 1] - tris[:, 0]  # vector of edge 1 on triangle (n_faces, 3)
      E2 = tris[:, 2] - tris[:, 0]  # vector of edge 2 on triangle (n_faces, 3)
  
      # batch cross product
      N = torch.cross(E1, E2)  # normal to E1 and E2, automatically batched to (n_faces, 3)
  
      invdet = 1. / -(torch.einsum('md,nd->mn', ray_d, N) + eps)  # inverse determinant (n_faces, 3)
  
      # __import__('ipdb').set_trace()
      A0 = ray_o[:, None] - tris[None, :, 0]  # (n_rays, 3) - (n_faces, 3) -> (n_rays, n_faces, 3) automatic broadcast
      DA0 = torch.cross(A0, ray_d[:, None].expand(
          *A0.shape))  # (n_rays, n_faces, 3) x (n_rays, 3) -> (n_rays, n_faces, 3) no automatic broadcast
  
      u = torch.einsum('mnd,nd->mn', DA0, E2) * invdet
      v = -torch.einsum('mnd,nd->mn', DA0, E1) * invdet
      t = torch.einsum('mnd,nd->mn', A0, N) * invdet
  
      return u, v, t
  