code: # -----------------------------------------------------------
  # Code adapted from: https://github.com/akanazawa/cmr/blob/master/utils/geometry.py
  # 
  # MIT License
  # 
  # Copyright (c) 2018 akanazawa
  # 
  # Permission is hereby granted, free of charge, to any person obtaining a copy
  # of this software and associated documentation files (the "Software"), to deal
  # in the Software without restriction, including without limitation the rights
  # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  # copies of the Software, and to permit persons to whom the Software is
  # furnished to do so, subject to the following conditions:
  # 
  # The above copyright notice and this permission notice shall be included in all
  # copies or substantial portions of the Software.
  # 
  # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  # SOFTWARE.
  # -----------------------------------------------------------
  
  # Geometry stuff.
  from __future__ import absolute_import
  from __future__ import division
  from __future__ import print_function
  
  import numpy as np
  
  def triangle_direction_intersection(tri, trg):
      '''
      Finds where an origin-centered ray going in direction trg intersects a triangle.
      Args:
          tri: 3 X 3 vertex locations. tri[0, :] is 0th vertex.
      Returns:
          alpha, beta, gamma
      '''
      p0 = np.copy(tri[0, :])
      # Don't normalize
      d1 = np.copy(tri[1, :]) - p0;
      d2 = np.copy(tri[2, :]) - p0;
      d = trg / np.linalg.norm(trg)
  
      mat = np.stack([d1, d2, d], axis=1)
  
      try:
        inv_mat = np.linalg.inv(mat)
      except np.linalg.LinAlgError:
        return False, 0
  
      a_b_mg = -1*np.matmul(inv_mat, p0)
      is_valid = (a_b_mg[0] >= 0) and (a_b_mg[1] >= 0) and ((a_b_mg[0] + a_b_mg[1]) <= 1) and (a_b_mg[2] < 0)
      if is_valid:
          return True, -a_b_mg[2]*d
      else:
          return False, 0
  
  
  def project_verts_on_mesh(verts, mesh_verts, mesh_faces):
      verts_out = np.copy(verts)
      for nv in range(verts.shape[0]):
          max_norm = 0
          vert = np.copy(verts_out[nv, :])
          for f in range(mesh_faces.shape[0]):
              face = mesh_faces[f]
              tri = mesh_verts[face, :]
              # is_v=True if it does intersect and returns the point
              is_v, pt = triangle_direction_intersection(tri, vert)
              # Take the furthest away intersection point
              if is_v and np.linalg.norm(pt) > max_norm:
                  max_norm = np.linalg.norm(pt)
                  verts_out[nv, :] = pt
  
      return verts_out
  
apis: numpy.matmul numpy.linalg.inv numpy.copy numpy.stack numpy.linalg.norm