code: def nearby_faces(mesh, points):
      
      points = np.asanyarray(points, dtype=np.float64)
      if not util.is_shape(points, (-1, 3)):
          raise ValueError('points must be (n,3)!')
  
      # an r-tree containing the axis aligned bounding box for every triangle
      rtree = mesh.triangles_tree
      # a kd-tree containing every vertex of the mesh
      kdtree = mesh.kdtree
  
      # query the distance to the nearest vertex to get AABB of a sphere
      distance_vertex = kdtree.query(points)[0].reshape((-1, 1))
      distance_vertex += tol.merge
  
      # axis aligned bounds
      bounds = np.column_stack((points - distance_vertex,
                                points + distance_vertex))
  
      # faces that intersect axis aligned bounding box
      candidates = [list(rtree.intersection(b)) for b in bounds]
  
      return candidates
docs: For each point find nearby faces relatively quickly.
  
      The closest point on the mesh to the queried point is guaranteed to be
      on one of the faces listed.
  
      Does this by finding the nearest vertex on the mesh to each point, and
      then returns all the faces that intersect the axis aligned bounding box
      centered at the queried point and extending to the nearest vertex.
  
      Parameters
      ----------
      mesh : Trimesh object
      points : (n,3) float , points in space
  
      Returns
      -----------
      candidates : (points,) int, sequence of indexes for mesh.faces