code: def ray_triangle_candidates(ray_origins,
                              ray_directions,
                              tree):
      
      ray_bounding = ray_bounds(ray_origins=ray_origins,
                                ray_directions=ray_directions,
                                bounds=tree.bounds)
      ray_candidates = [[]] * len(ray_origins)
      ray_id = [[]] * len(ray_origins)
  
      for i, bounds in enumerate(ray_bounding):
          ray_candidates[i] = np.array(list(tree.intersection(bounds)),
                                       dtype=np.int)
          ray_id[i] = np.ones(len(ray_candidates[i]), dtype=np.int) * i
  
      ray_id = np.hstack(ray_id)
      ray_candidates = np.hstack(ray_candidates)
  
      return ray_candidates, ray_id
docs: Do broad- phase search for triangles that the rays
      may intersect.
  
      Does this by creating a bounding box for the ray as it
      passes through the volume occupied by the tree
  
      Parameters
      ------------
      ray_origins:      (m,3) float, ray origin points
      ray_directions:   (m,3) float, ray direction vectors
      tree:             rtree object, contains AABB of each triangle
  
      Returns
      ----------
      ray_candidates: (n,) int, triangle indexes
      ray_id:         (n,) int, corresponding ray index for a triangle candidate