code: def natural_neighbor_to_points(points, values, xi):
      r
      tri = Delaunay(points)
  
      members, triangle_info = geometry.find_natural_neighbors(tri, xi)
  
      img = np.empty(shape=(xi.shape[0]), dtype=values.dtype)
      img.fill(np.nan)
  
      for ind, (grid, neighbors) in enumerate(members.items()):
  
          if len(neighbors) > 0:
  
              points_transposed = np.array(points).transpose()
              img[ind] = natural_neighbor_point(points_transposed[0], points_transposed[1],
                                                values, xi[grid], tri, neighbors, triangle_info)
  
      return img
docs: r"""Generate a natural neighbor interpolation to the given points.
  
      This assigns values to the given interpolation points using the Liang and Hale
      [Liang2010]_. approach.
  
      Parameters
      ----------
      points: array_like, shape (n, 2)
          Coordinates of the data points.
      values: array_like, shape (n,)
          Values of the data points.
      xi: array_like, shape (M, 2)
          Points to interpolate the data onto.
  
      Returns
      -------
      img: (M,) ndarray
          Array representing the interpolated values for each input point in `xi`
  
      See Also
      --------
      natural_neighbor_to_grid