docs: def line_integration(image, angle, decay, sigma):
      '''Integrate the image along the given angle
  
      DIC images are the directional derivative of the underlying
      image. This filter reconstructs the original image by integrating
      along that direction.
  
      image - a 2-dimensional array
  
      angle - shear angle in radians. We integrate perpendicular to this angle
  
      decay - an exponential decay applied to the integration
  
      sigma - the standard deviation of a Gaussian which is used to
              smooth the image in the direction parallel to the shear angle.
      '''
      #
      # Normalize the image so that the mean is zero
      #
      normalized = image - np.mean(image)
      #
      # Rotate the image so the J direction is perpendicular to the shear angle.
      #
      rotated = scind.rotate(normalized, -angle)
      #
      # Smooth in only the i direction
      #
      smoothed = scind.gaussian_filter1d(rotated, sigma) if sigma > 0 else rotated
      #
      # We want img_out[:,j+1] to be img_out[:,j] * decay + img[j+1]
      # Could be done by convolution with a ramp, maybe in FFT domain,
      # but we just do a bunch of steps here.
      #
      result_fwd = smoothed.copy()
      for i in range(1,result_fwd.shape[0]):
          result_fwd[i] += result_fwd[i-1] * decay
      result_rev = smoothed.copy()
      for i in reversed(range(result_rev.shape[0]-1)):
          result_rev[i] += result_rev[i+1] * decay
      result = (result_fwd - result_rev) / 2
      #
      # Rotate and chop result
      #
      result = scind.rotate(result, angle)
      ipad = int((result.shape[0] - image.shape[0]) / 2)
      jpad = int((result.shape[1] - image.shape[1]) / 2)
      result = result[ipad:(ipad + image.shape[0]),
                      jpad:(jpad + image.shape[1])]
      #
      # Scale the resultant image similarly to the output.
      #
      img_min, img_max = np.min(image), np.max(image)
      result_min, result_max = np.min(result), np.max(result)
      if (img_min == img_max) or (result_min == result_max):
          return np.zeros(result.shape)
      result = (result - result_min) / (result_max - result_min)
      result = img_min + result * (img_max - img_min)
      return result