add median clippign to define background

python/lsst/meas/algorithms/gp_interpolation.py

@@ -27,6 +27,15 @@ def get_name_pkl():
     file_name = f"out_test_{n_file}.pkl"
     return os.path.join(rep, file_name)
 
+@jit
+def median_with_mad_clipping(data, mad_multiplier=2.0):
+
+    median = jnp.median(data)
+    mad = jnp.median(jnp.abs(data - median))
+    clipping_range = mad_multiplier * mad
+    clipped_data = jnp.clip(data, median - clipping_range, median + clipping_range)
+    median_clipped = jnp.median(clipped_data)
+    return median_clipped
 
 @jit
 def jax_pdist_squared(X):
@@ -324,7 +333,7 @@ def interpolate_sub_masked_image(self, sub_masked_image):
         # Do GP interpolation if bad pixel found.
         else:
             # gp interpolation
-            mean = np.mean(good_pixel[:, 2:])
+            mean = median_with_mad_clipping(good_pixel[:, 2:]) # np.mean(good_pixel[:, 2:])
             sub_image_array = sub_masked_image.getVariance().array
             white_noise = np.sqrt(
                 np.mean(sub_image_array[np.isfinite(sub_image_array)])
@@ -365,7 +374,7 @@ def interpolate_sub_masked_image(self, sub_masked_image):
             updateMaskFromArray(sub_masked_image.mask, bad_pixel, self.interpBit)
             return sub_masked_image
 
-def interpolateOverDefectsGP(image, fwhm, badList, method="treegp", bin_spacing=15, threshold_subdivide=20000):
+def interpolateOverDefectsGP(image, fwhm, badList, method="treegp", bin_spacing=25, threshold_subdivide=20000):
     """
     Interpolates over defects in an image using Gaussian Process interpolation.