JP-3997: Enable interpolation in pixel map creation for imaging data in resample

changes/10137.outlier_detection.rst

@@ -0,0 +1 @@
+Add ``pixmap_stepsize`` and ``pixmap_order`` parameters to allow interpolated calculation of pixel map for imaging mosaics when resampling

changes/10137.resample.rst

@@ -0,0 +1 @@
+Add ``pixmap_stepsize`` and ``pixmap_order`` parameters to allow interpolated calculation of pixel map for imaging mosaics; improves runtime of calwebb_image3 by 20 percent for a 30-image mosaic tested.

jwst/outlier_detection/imaging.py

@@ -35,6 +35,8 @@ def detect_outliers(
     fillval,
     in_memory,
     make_output_path,
+    pixmap_stepsize=1,
+    pixmap_order=1,
 ):
     """
     Flag outliers in imaging data.
@@ -78,6 +80,12 @@ def detect_outliers(
     make_output_path : function
         The functools.partial instance to pass to save_blot. Must be
         specified if save_blot is True.
+    pixmap_stepsize : float, optional
+        Indicates the spacing at which WCS is evaluated when computing pixel map.
+        Larger step sizes result in faster performance at the cost of accuracy.
+        Interpolation is only performed if ``pixmap_stepsize > 1``. Default is 1.
+    pixmap_order : int, optional
+        Interpolating spline order for pixel map computation. Must be 1 or 3. Default is 1.
 
     Returns
     -------
@@ -105,6 +113,8 @@ def detect_outliers(
             enable_ctx=False,
             enable_var=False,
             compute_err=None,
+            pixmap_order=pixmap_order,
+            pixmap_stepsize=pixmap_stepsize,
         )
         median_data, median_wcs = median_with_resampling(
             input_models,
@@ -139,6 +149,8 @@ def detect_outliers(
                     backg,
                     save_blot=save_intermediate_results,
                     make_output_path=make_output_path,
+                    pixmap_stepsize=pixmap_stepsize,
+                    pixmap_order=pixmap_order,
                 )
             else:
                 flag_model_crs(image, median_data, snr1)

jwst/outlier_detection/outlier_detection_step.py

@@ -56,6 +56,8 @@ class OutlierDetectionStep(Step):
         good_bits = string(default="~DO_NOT_USE")  # DQ flags to allow
         search_output_file = boolean(default=False)
         in_memory = boolean(default=True) # in_memory flag ignored if run within the pipeline; set at pipeline level instead
+        pixmap_stepsize = float(default=1.0)  # Interpolation step size for pixel map; interpolation is used for stepsize > 1
+        pixmap_order = integer(default=1)  # Spline order for pixel mapping, must be 1 or 3
     """  # noqa: E501
 
     def process(self, input_data):
@@ -136,6 +138,8 @@ def process(self, input_data):
                 self.fillval,
                 self.in_memory,
                 self.make_output_path,
+                self.pixmap_stepsize,
+                self.pixmap_order,
             )
         elif mode == "spec":
             result_models = spec.detect_outliers(

jwst/outlier_detection/tests/test_imaging.py

@@ -279,3 +279,23 @@ def test_skip_one_exposure_imaging():
     # Input is not modified
     assert result is not model
     assert model.meta.cal_step.outlier_detection is None
+
+
+def test_outlier_step_pixmap_interpolation(mirimage_three_sci):
+    """Test that outlier detection runs with non-default pixmap interpolation settings."""
+    container = ModelLibrary(mirimage_three_sci)
+
+    # Drop a CR on the science array
+    with container:
+        zeroth = container.borrow(0)
+        zeroth.data[12, 12] += 1
+        container.shelve(zeroth)
+
+    # Run outlier detection with non-default pixmap interpolation settings
+    OutlierDetectionStep.call(container, in_memory=True, pixmap_stepsize=5, pixmap_order=3)
+
+    # Verify CR is flagged
+    with container:
+        zeroth = container.borrow(0)
+        assert zeroth.dq[12, 12] == helpers.OUTLIER_DO_NOT_USE
+        container.shelve(zeroth, modify=False)

jwst/outlier_detection/utils.py

@@ -332,6 +332,8 @@ def flag_resampled_model_crs(
     median_err=None,
     save_blot=False,
     make_output_path=None,
+    pixmap_stepsize=1,
+    pixmap_order=1,
 ):
     """
     Flag outliers in a resampled model, updating DQ array in place.
@@ -364,13 +366,22 @@ def flag_resampled_model_crs(
     make_output_path : function
         The functools.partial instance to pass to save_blot. Must be
         specified if save_blot is True.
+    pixmap_stepsize : float, optional
+        Indicates the spacing at which WCS is evaluated when computing pixel map.
+        WCS coordinates of the full pixel map is computed by interpolating over
+        this sparse pixel map when ``pixmap_stepsize > 1``. Larger step sizes result in
+        faster performance at the cost of accuracy. Default is 1.
+    pixmap_order : int, optional
+        Interpolating spline order for pixel map computation. Must be 1 or 3. Default is 1.
     """
     blot = gwcs_blot(
         median_data=median_data,
         median_wcs=median_wcs,
         blot_shape=input_model.data.shape,
         blot_wcs=input_model.meta.wcs,
         fillval=np.nan,
+        pixmap_stepsize=pixmap_stepsize,
+        pixmap_order=pixmap_order,
     )
     if median_err is not None:
         blot_err = gwcs_blot(
@@ -379,6 +390,8 @@ def flag_resampled_model_crs(
             blot_shape=input_model.data.shape,
             blot_wcs=input_model.meta.wcs,
             fillval=np.nan,
+            pixmap_stepsize=pixmap_stepsize,
+            pixmap_order=pixmap_order,
         )
     else:
         blot_err = None
@@ -505,6 +518,8 @@ def flag_crs_in_models_with_resampling(
             median_err=median_err,
             save_blot=save_blot,
             make_output_path=make_output_path,
+            pixmap_stepsize=1,
+            pixmap_order=1,
         )
 
 

jwst/resample/resample.py

@@ -56,6 +56,8 @@ def __init__(
         report_var=True,
         compute_err=None,
         asn_id=None,
+        pixmap_stepsize=1,
+        pixmap_order=1,
     ):
         """
         Initialize the ResampleImage object.
@@ -267,6 +269,15 @@ def __init__(
         asn_id : str, None, optional
             The association id. The id is what appears in
             the :ref:`asn-jwst-naming`.
+
+        pixmap_stepsize : float, optional
+            Indicates the spacing at which WCS is evaluated when computing pixel map.
+            WCS coordinates of the full pixel map is computed by interpolating over
+            this sparse pixel map when ``pixmap_stepsize > 1``. Larger step sizes result in
+            faster performance at the cost of accuracy. Default is 1.
+
+        pixmap_order : int, optional
+            Interpolating spline order for pixel map computation. Must be 1 or 3. Default is 1.
         """
         self.input_models = input_models
         self.output_jwst_model = None
@@ -340,6 +351,8 @@ def __init__(
             enable_ctx=enable_ctx,
             enable_var=enable_var,
             compute_err=compute_err,
+            pixmap_stepsize=pixmap_stepsize,
+            pixmap_order=pixmap_order,
         )
 
     def input_model_to_dict(self, model, weight_type, enable_var, compute_err):

jwst/resample/resample_spec.py

@@ -207,7 +207,13 @@ def __init__(self, input_models, good_bits=0, output_wcs=None, wcs_pars=None, **
         library = ModelLibrary(input_models, on_disk=False)
 
         super().__init__(
-            library, good_bits=good_bits, output_wcs=output_wcs_dict, wcs_pars=None, **kwargs
+            library,
+            good_bits=good_bits,
+            output_wcs=output_wcs_dict,
+            wcs_pars=None,
+            pixmap_stepsize=1,
+            pixmap_order=1,
+            **kwargs,
         )
         self.intermediate_suffix = "outlier_s2d"
 

jwst/resample/resample_step.py

@@ -40,6 +40,8 @@ class ResampleStep(Step):
         enable_ctx = boolean(default=True)  # Compute and report the context array
         enable_err = boolean(default=True)  # Compute and report the err array
         report_var = boolean(default=True)  # Report the variance array
+        pixmap_stepsize = float(default=1.0)  # Interpolation step size for pixel map; interpolation is used for stepsize > 1
+        pixmap_order = integer(default=1)  # Spline order for pixel mapping, must be 1 or 3
     """  # noqa: E501
 
     reference_file_types: list = []
@@ -211,6 +213,8 @@ def get_drizpars(self):
             "weight_type": self.weight_type,
             "good_bits": GOOD_BITS,
             "blendheaders": self.blendheaders,
+            "pixmap_stepsize": self.pixmap_stepsize,
+            "pixmap_order": self.pixmap_order,
         }
 
         # Custom output WCS parameters.

jwst/resample/tests/test_resample_step.py

@@ -1704,3 +1704,23 @@ def test_spec_skip_cube():
     # Input is not modified
     assert result is not model
     assert model.meta.cal_step.resample is None
+
+
+def test_resample_imaging_pixmap_interpolation(nircam_rate):
+    """Test that resample gives similar results with non-default pixmap interpolation settings."""
+    img = AssignWcsStep.call(nircam_rate, sip_approx=False)
+    # give the data some structure
+    img.data = np.random.default_rng(seed=77).random(img.data.shape)
+    img.var_rnoise = np.ones_like(img.data) * 1e-3
+
+    # resampling rotates the image a little bit, i.e. this is indeed nontrivial
+    ref = ResampleStep.call(img, pixmap_order=1, pixmap_stepsize=1)
+    res = ResampleStep.call(img, pixmap_order=3, pixmap_stepsize=10)
+
+    # catch issue where bad inputs can cause all-NaN output when variance is zero everywhere
+    assert not np.all(np.isnan(res.data))
+    # ensure results are very similar
+    assert_allclose(res.data, ref.data, rtol=1.0e-6, atol=1.0e-9)
+    # ensure results are not identical (i.e. pixmap settings actually did something)
+    with pytest.raises(AssertionError):
+        assert_allclose(res.data, ref.data)