code: def bilinear_ray_sum(image, theta, ray_position):
      """
      Compute the projection of an image along a ray.
  
      Parameters
      ----------
      image : 2D array, dtype=float
          Image to project.
      theta : float
          Angle of the projection
      ray_position : float
          Position of the ray within the projection
  
      Returns
      -------
      projected_value : float
          Ray sum along the projection
      norm_of_weights :
          A measure of how long the ray's path through the reconstruction
          circle was
      """
      theta = theta / 180.0 * np.pi
      radius = image.shape[0] // 2 - 1
      projection_center = image.shape[0] // 2
      rotation_center = image.shape[0] // 2
      # (s, t) is the (x, y) system rotated by theta
      t = ray_position - projection_center
      # s0 is the half-length of the ray's path in the reconstruction circle
      s0 = (
          sqrt(radius * radius - t * t)
          if radius * radius >= t * t
          else 0.0
      )
      Ns = 2 * int(ceil(2 * s0))  # number of steps along the ray
      ray_sum = 0.0
      weight_norm = 0.0
  
      if Ns > 0:
          # step length between samples
          ds = 2 * s0 / Ns
          dx = -ds * cos(theta)
          dy = -ds * sin(theta)
          # point of entry of the ray into the reconstruction circle
          x0 = s0 * cos(theta) - t * sin(theta)
          y0 = s0 * sin(theta) + t * cos(theta)
          for k in range(Ns + 1):
              x = x0 + k * dx
              y = y0 + k * dy
              index_i = x + rotation_center
              index_j = y + rotation_center
              i = int(floor(index_i))
              j = int(floor(index_j))
              di = index_i - floor(index_i)
              dj = index_j - floor(index_j)
              # Use linear interpolation between values
              # Where values fall outside the array, assume zero
              if i > 0 and j > 0:
                  weight = (1.0 - di) * (1.0 - dj) * ds
                  ray_sum += weight * image[i, j]
                  weight_norm += weight * weight
              if i > 0 and j < image.shape[1] - 1:
                  weight = (1.0 - di) * dj * ds
                  ray_sum += weight * image[i, j + 1]
                  weight_norm += weight * weight
              if i < image.shape[0] - 1 and j > 0:
                  weight = di * (1 - dj) * ds
                  ray_sum += weight * image[i + 1, j]
                  weight_norm += weight * weight
              if i < image.shape[0] - 1 and j < image.shape[1] - 1:
                  weight = di * dj * ds
                  ray_sum += weight * image[i + 1, j + 1]
                  weight_norm += weight * weight
  
      return ray_sum, weight_norm
  
imports: import numpy
  from numpy import cos, sin, floor, ceil, sqrt
label_desc: Use image node_id 1 for calling a global function with example usage: bilinear_ray_sum(image, theta, ray_position) and returns: ray_sum,  weight_norm