code: def bilinear_ray_update(
      image, image_update, theta, ray_position, projected_value
  ):
      """Compute the update along a ray using bilinear interpolation.
  
      Parameters
      ----------
      image : 2D array, dtype=float
          Current reconstruction estimate.
      image_update : 2D array, dtype=float
          Array of same shape as ``image``. Updates will be added to this array.
      theta : float
          Angle of the projection.
      ray_position : float
          Position of the ray within the projection.
      projected_value : float
          Projected value (from the sinogram).
  
      Returns
      -------
      deviation :
          Deviation before updating the image.
      """
      ray_sum, weight_norm = bilinear_ray_sum(
          image, theta, ray_position
      )
      if weight_norm > 0.0:
          deviation = -(ray_sum - projected_value) / weight_norm
      else:
          deviation = 0.0
      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))
      # beta for equiripple Hamming window
      hamming_beta = 0.46164
  
      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)
              hamming_window = (1 - hamming_beta) - hamming_beta * cos(
                  2 * np.pi * k / (Ns - 1)
              )
              if i > 0 and j > 0:
                  image_update[i, j] += (
                      deviation
                      * (1.0 - di)
                      * (1.0 - dj)
                      * ds
                      * hamming_window
                  )
              if i > 0 and j < image.shape[1] - 1:
                  image_update[i, j + 1] += (
                      deviation * (1.0 - di) * dj * ds * hamming_window
                  )
              if i < image.shape[0] - 1 and j > 0:
                  image_update[i + 1, j] += (
                      deviation * di * (1 - dj) * ds * hamming_window
                  )
              if i < image.shape[0] - 1 and j < image.shape[1] - 1:
                  image_update[i + 1, j + 1] += (
                      deviation * di * dj * ds * hamming_window
                  )
  
      return deviation
  
imports: import numpy
  from numpy import cos, sin, floor, ceil, sqrt
label_desc: Use image node_id 2 for calling a global function with example usage: bilinear_ray_update(image, image_update, theta, ray_position, projected_value) and returns: deviation