Repository: spacetelescope/jwst

Issue:
Background over-subtracted 
I was running the pipeline on a set of simulated MIRI images (3 dithers). After finishing stage 1 & 2 procedures, I got three images with background ~ 8. Then I run stage 3, and merged the three images together. The merged image has background ~ -8, instead of ~0. 

I investigated each step in stage 3, i.e. tweakreg -> skymatch -> outlier_detection -> resample. After skymatch, the three images still have background ~ 8. But after outlier_detection, the background became ~ 0. After resample (which merges the three images), the background became ~ -8. So, outlier_detection and resample both subtracted the background, leading to an over-subtraction. 

Any idea what's wrong?

Thanks,
Guang



Discussion:
codecov[bot]: # [Codecov](https://codecov.io/gh/spacetelescope/jwst/pull/4858?src=pr&el=h1) Report
> Merging [#4858](https://codecov.io/gh/spacetelescope/jwst/pull/4858?src=pr&el=desc) into [master](https://codecov.io/gh/spacetelescope/jwst/commit/03f66f6572c947768ab590d49c80ff315dec9c31&el=desc) will **increase** coverage by `0.57%`.
> The diff coverage is `71.42%`.

[![Impacted file tree graph](https://codecov.io/gh/spacetelescope/jwst/pull/4858/graphs/tree.svg?width=650&height=150&src=pr&token=Utf5Zs9g7z)](https://codecov.io/gh/spacetelescope/jwst/pull/4858?src=pr&el=tree)

```diff
@@            Coverage Diff             @@
##           master    #4858      +/-   ##
==========================================
+ Coverage   52.37%   52.94%   +0.57%     
==========================================
  Files         401      401              
  Lines       35428    35454      +26     
  Branches     5485     5489       +4     
==========================================
+ Hits        18554    18772     +218     
+ Misses      15747    15523     -224     
- Partials     1127     1159      +32     
```

| Flag | Coverage Δ | |
|---|---|---|
| #unit | `52.94% <71.42%> (+0.57%)` | :arrow_up: |

| [Impacted Files](https://codecov.io/gh/spacetelescope/jwst/pull/4858?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [jwst/assign\_wcs/pointing.py](https://codecov.io/gh/spacetelescope/jwst/pull/4858/diff?src=pr&el=tree#diff-andzdC9hc3NpZ25fd2NzL3BvaW50aW5nLnB5) | `86.95% <71.42%> (+1.24%)` | :arrow_up: |
| [jwst/lib/set\_telescope\_pointing.py](https://codecov.io/gh/spacetelescope/jwst/pull/4858/diff?src=pr&el=tree#diff-andzdC9saWIvc2V0X3RlbGVzY29wZV9wb2ludGluZy5weQ==) | `85.25% <0.00%> (-3.88%)` | :arrow_down: |
| [jwst/resample/gwcs\_drizzle.py](https://codecov.io/gh/spacetelescope/jwst/pull/4858/diff?src=pr&el=tree#diff-andzdC9yZXNhbXBsZS9nd2NzX2RyaXp6bGUucHk=) | `57.00% <0.00%> (+2.00%)` | :arrow_up: |
| [jwst/resample/resample\_utils.py](https://codecov.io/gh/spacetelescope/jwst/pull/4858/diff?src=pr&el=tree#diff-andzdC9yZXNhbXBsZS9yZXNhbXBsZV91dGlscy5weQ==) | `56.48% <0.00%> (+12.03%)` | :arrow_up: |
| [jwst/assign\_wcs/util.py](https://codecov.io/gh/spacetelescope/jwst/pull/4858/diff?src=pr&el=tree#diff-andzdC9hc3NpZ25fd2NzL3V0aWwucHk=) | `61.93% <0.00%> (+13.92%)` | :arrow_up: |
| [jwst/outlier\_detection/outlier\_detection\_step.py](https://codecov.io/gh/spacetelescope/jwst/pull/4858/diff?src=pr&el=tree#diff-andzdC9vdXRsaWVyX2RldGVjdGlvbi9vdXRsaWVyX2RldGVjdGlvbl9zdGVwLnB5) | `40.17% <0.00%> (+21.36%)` | :arrow_up: |
| [jwst/outlier\_detection/outlier\_detection.py](https://codecov.io/gh/spacetelescope/jwst/pull/4858/diff?src=pr&el=tree#diff-andzdC9vdXRsaWVyX2RldGVjdGlvbi9vdXRsaWVyX2RldGVjdGlvbi5weQ==) | `72.09% <0.00%> (+31.62%)` | :arrow_up: |
| [jwst/resample/resample.py](https://codecov.io/gh/spacetelescope/jwst/pull/4858/diff?src=pr&el=tree#diff-andzdC9yZXNhbXBsZS9yZXNhbXBsZS5weQ==) | `78.16% <0.00%> (+59.77%)` | :arrow_up: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/spacetelescope/jwst/pull/4858?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/spacetelescope/jwst/pull/4858?src=pr&el=footer). Last update [03f66f6...372c92a](https://codecov.io/gh/spacetelescope/jwst/pull/4858?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

hbushouse: Needs change log entries.
jdavies-st: Before merging, I'm running regression tests on this branch.  I expect it will affect all `crf`, `s2d` and `i2d` products.
jdavies-st: And it looks like this fix causes a bunch of errors for the spectra modes as generally these get passed a `SourceModelContainer` that is a `MultiExposureModel` masquerading as a `ModelContainer`, and it does not have a `copy()` method.  So I'l have to implement a more subtle fix further down in the code.
hbushouse: Well of course! Nothing is ever that easy. I guess we need to rescind your regal status.
jdavies-st: It looks like it is now a 2-line fix.  The background level gets subtracted once in `resample`, and once in `flag_cr`.  The issue here is that both operate on the inputs in-place, and do not make copies.  *tsk* *tsk*

Patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index 157cab225..2cf457872 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -20,6 +20,11 @@ flat_field
 
 - Implement applying the inverse operation. [#5302]
 
+outlier_detection
+-----------------
+
+- Fix bug where background was being subtracted on the input data [#4858]
+
 pathloss
 --------
 
@@ -41,6 +46,7 @@ pipeline
 
 - Implement master background subtraction in Spec2Pipeline for NIRSpec MOS data [#5302]
 
+
 0.17.0 (2020-08-28)
 ===================
 
diff --git a/jwst/assign_wcs/pointing.py b/jwst/assign_wcs/pointing.py
index 115f3b7c3..bf697aa96 100644
--- a/jwst/assign_wcs/pointing.py
+++ b/jwst/assign_wcs/pointing.py
@@ -86,17 +86,6 @@ def _roll_angle_from_matrix(matrix, v2, v3):
     return new_roll
 
 
-#def create_fitswcs_transform(input_model):
-    #"""
-
-    #"""
-    #ff = fits_support.to_fits(input_model._instance, input_model._schema)
-    #hdu = fits_support.get_hdu(ff._hdulist, "PRIMARY")
-    #header = hdu.header
-    #transform = gwutils.make_fitswcs_transform(header)
-    #return transform
-
-
 def wcsinfo_from_model(input_model):
     """
     Create a dict {wcs_keyword: array_of_values} pairs from a data model.
@@ -223,19 +212,26 @@ def create_fitswcs(inp, input_frame=None):
         wcsinfo = wcsinfo_from_model(inp)
         wavetable = None
         spatial_axes, spectral_axes, unknown = gwutils.get_axes(wcsinfo)
-        sp_axis = spectral_axes[0]
-        if wcsinfo['CTYPE'][sp_axis] == 'WAVE-TAB':
-            wavetable = inp.wavetable
-        transform = fitswcs_transform_from_model(wcsinfo, wavetable)
+        if spectral_axes:
+            sp_axis = spectral_axes[0]
+            if wcsinfo['CTYPE'][sp_axis] == 'WAVE-TAB':
+                wavetable = inp.wavetable
+        transform = fitswcs_transform_from_model(wcsinfo, wavetable=wavetable)
         output_frame = frame_from_model(wcsinfo)
-    #elif isinstance(inp, str):
-        #transform = create_fitswcs_transform(inp)
-        #output_frame = frame_from_fits(inp)
     else:
         raise TypeError("Input is expected to be a DataModel instance or a FITS file.")
 
     if input_frame is None:
-        input_frame = "detector"
+        wcsaxes = wcsinfo['WCSAXES']
+        if wcsaxes == 2:
+            input_frame = cf.Frame2D(name="detector")
+        elif wcsaxes == 3:
+            input_frame = cf.CoordinateFrame(name="detector", naxes=3,
+                axes_order=(0, 1, 2), unit=(u.pix, u.pix, u.pix),
+                axes_type=["SPATIAL", "SPATIAL", "SPECTRAL"],
+                axes_names=('x', 'y', 'z'), axis_physical_types=None)
+        else:
+            raise TypeError(f"WCSAXES is expected to be 2 or 3, instead it is {wcsaxes}")
     pipeline = [(input_frame, transform),
                (output_frame, None)]
 
diff --git a/jwst/assign_wcs/util.py b/jwst/assign_wcs/util.py
index 5db3e77da..8c239d7e2 100644
--- a/jwst/assign_wcs/util.py
+++ b/jwst/assign_wcs/util.py
@@ -259,13 +259,8 @@ def wcs_from_footprints(dmodels, refmodel=None, transform=None, bounding_box=Non
 
     out_frame = refmodel.meta.wcs.output_frame
     input_frame = dmodels[0].meta.wcs.input_frame
-    wnew = wcs_from_fiducial(fiducial, coordinate_frame=out_frame, projection=prj, transform=transform)
-    # temporary fix before gwcs 0.14 is released
-    # wnew = wcs_from_fiducial(fiducial, coordinate_frame=out_frame, projection=prj,
-    #                          transform=transform, input_frame=input_frame)
-    pipe = wnew.pipeline[:]
-    pipe[0] = (input_frame, pipe[0][1])
-    wnew = WCS(pipe)
+    wnew = wcs_from_fiducial(fiducial, coordinate_frame=out_frame, projection=prj,
+                             transform=transform, input_frame=input_frame)
 
     footprints = [w.footprint().T for w in wcslist]
     domain_bounds = np.hstack([wnew.backward_transform(*f) for f in footprints])
diff --git a/jwst/outlier_detection/outlier_detection.py b/jwst/outlier_detection/outlier_detection.py
index 0915e14cb..fd5b6a6da 100644
--- a/jwst/outlier_detection/outlier_detection.py
+++ b/jwst/outlier_detection/outlier_detection.py
@@ -350,7 +350,9 @@ class OutlierDetection:
         """
 
         for image, blot in zip(self.input_models, blot_models):
-            flag_cr(image, blot, **self.outlierpars)
+            image_copy = image.copy()
+            flag_cr(image_copy, blot, **self.outlierpars)
+            image.dq = image_copy.dq.copy()
 
         if self.converted:
             # Make sure actual input gets updated with new results
@@ -390,12 +392,13 @@ def flag_cr(sci_image, blot_image, **pars):
     snr1, snr2 = [float(val) for val in pars.get('snr', '5.0 4.0').split()]
     scl1, scl2 = [float(val) for val in pars.get('scale', '1.2 0.7').split()]
 
-    # Get background level if it has been subtracted
-    if (sci_image.meta.background.subtracted is True and
+    # Get background level of science data if it has not been subtracted, so it
+    # can be added into the level of the blotted data, which has been
+    # background-subtracted
+    if (sci_image.meta.background.subtracted is False and
         sci_image.meta.background.level is not None):
         subtracted_background = sci_image.meta.background.level
-        log.debug(f"Including subtracted background ({subtracted_background}) "
-            "back into blotted image")
+        log.debug(f"Adding background level {subtracted_background} to blotted image")
     else:
         # No subtracted background.  Allow user-set value, which defaults to 0
         subtracted_background = backg
@@ -418,7 +421,6 @@ def flag_cr(sci_image, blot_image, **pars):
     #
     # Model the noise and create a CR mask
     diff_noise = np.abs(sci_data - blot_data)
-    # ta = np.sqrt(np.abs(blot_data + subtracted_background) + rn ** 2)
     ta = np.sqrt(np.abs(blot_data + subtracted_background) + err_data ** 2)
     t2 = scl1 * blot_deriv + snr1 * ta
     tmp1 = np.logical_not(np.greater(diff_noise, t2))
diff --git a/jwst/outlier_detection/outlier_detection_step.py b/jwst/outlier_detection/outlier_detection_step.py
index 6fb3f1144..dd7a9bf45 100644
--- a/jwst/outlier_detection/outlier_detection_step.py
+++ b/jwst/outlier_detection/outlier_detection_step.py
@@ -37,7 +37,7 @@ class OutlierDetectionStep(Step):
 
     Parameters
     -----------
-    user_input : asn file or ModelContainer
+    input_data : asn file or ModelContainer
         Single filename association table, or a datamodels.ModelContainer.
 
     """
@@ -64,9 +64,9 @@ class OutlierDetectionStep(Step):
         search_output_file = boolean(default=False)
     """
 
-    def process(self, user_input):
+    def process(self, input_data):
         """Perform outlier detection processing on input data."""
-        with datamodels.open(user_input) as input_models:
+        with datamodels.open(input_data) as input_models:
             self.input_models = input_models
             if not isinstance(self.input_models, datamodels.ModelContainer):
                 self.input_container = False
diff --git a/jwst/resample/resample.py b/jwst/resample/resample.py
index a1ceb0ac5..41c30da34 100644
--- a/jwst/resample/resample.py
+++ b/jwst/resample/resample.py
@@ -126,7 +126,9 @@ class ResampleData:
                                             kernel=self.drizpars['kernel'],
                                             fillval=self.drizpars['fillval'])
 
-            for n, img in enumerate(exposure):
+            for n, img_exp in enumerate(exposure):
+                # Make a copy as this is operated on in-place below
+                img = img_exp.copy()
                 exposure_times['start'].append(img.meta.exposure.start_time)
                 exposure_times['end'].append(img.meta.exposure.end_time)
 


Test patch:
diff --git a/jwst/assign_wcs/tests/test_wcs.py b/jwst/assign_wcs/tests/test_wcs.py
index 0589e21a3..ba0579549 100644
--- a/jwst/assign_wcs/tests/test_wcs.py
+++ b/jwst/assign_wcs/tests/test_wcs.py
@@ -4,13 +4,13 @@ from astropy import units as u
 from astropy import wcs
 from asdf.tests import helpers
 
-from .. import pointing
-from ...transforms import models
-from ...datamodels import ImageModel, fits_support
+from jwst.assign_wcs import pointing
+from jwst.transforms import models
+from jwst.datamodels import ImageModel, CubeModel, fits_support
 
 
-
-def _create_model_2d():
+@pytest.fixture
+def create_model_2d():
     im = ImageModel()
     im.meta.wcsinfo.crpix1 = 2.5
     im.meta.wcsinfo.crpix2 = 3
@@ -28,8 +28,23 @@ def _create_model_2d():
     return im
 
 
-def _create_model_3d():
-    im = _create_model_2d()
+@pytest.fixture
+def create_model_3d():
+    im = CubeModel()
+    im.meta.wcsinfo.crpix1 = 2.5
+    im.meta.wcsinfo.crpix2 = 3
+    im.meta.wcsinfo.crval1 = 5.6
+    im.meta.wcsinfo.crval2 = -72.3
+    im.meta.wcsinfo.wcsaxes = 2
+    im.meta.wcsinfo.cunit1 = 'deg'
+    im.meta.wcsinfo.cunit2 = 'deg'
+    im.meta.wcsinfo.ctype1 = 'RA---TAN'
+    im.meta.wcsinfo.ctype2 = 'DEC--TAN'
+    im.meta.wcsinfo.pc1_1 = 1.
+    im.meta.wcsinfo.pc1_2 = 0
+    im.meta.wcsinfo.pc2_1 = 0.
+    im.meta.wcsinfo.pc2_2 = 1.
+
     im.meta.wcsinfo.crpix3 = 1
     im.meta.wcsinfo.crval3 = 100
     im.meta.wcsinfo.wcsaxes = 3
@@ -83,12 +98,13 @@ def test_v23_to_sky():
     assert_allclose(radec, expected_ra_dec, atol=1e-10)
 
 
-def test_frame_from_model(tmpdir):
+def test_frame_from_model_3d(tmpdir, create_model_3d):
     """ Tests creating a frame from a data model. """
     # Test CompositeFrame initialization (celestial and spectral)
-    im = _create_model_3d()
+    im = create_model_3d
     frame = pointing.frame_from_model(im)
     radec, lam = frame.coordinates(1, 2, 3)
+
     assert_allclose(radec.spherical.lon.value, 1)
     assert_allclose(radec.spherical.lat.value, 2)
     u.allclose(lam, 3 * u.um)
@@ -97,25 +113,30 @@ def test_frame_from_model(tmpdir):
     im.meta.wcsinfo.ctype1 = 'ALPHA1A'
     im.meta.wcsinfo.ctype2 = 'BETA1A'
     frame = pointing.frame_from_model(im)
+
     assert frame.frames[1].name == 'ALPHA1A_BETA1A'
     assert frame.frames[1].axes_names == ('ALPHA1A', 'BETA1A')
 
     tree = {'frame': frame}
     helpers.assert_roundtrip_tree(tree, tmpdir)
 
+
+def test_frame_from_model_2d(tmpdir, create_model_2d):
+    """ Tests creating a frame from a data model. """
     # Test 2D spatial custom frame
-    im = _create_model_2d()
+    im = create_model_2d
     frame = pointing.frame_from_model(im)
+
     assert frame.name == "sky"
     assert frame.axes_names == ("RA", "DEC")
+
     tree = {'frame': frame}
     helpers.assert_roundtrip_tree(tree, tmpdir)
 
 
-@pytest.mark.filterwarnings("ignore: The WCS transformation has more axes")
-def test_create_fits_wcs(tmpdir):
-    """ Test GWCS vs FITS WCS results. """
-    im = _create_model_3d()
+def test_create_fitswcs(tmpdir, create_model_3d):
+    """GWCS from create_fitswcs function and astropy.wcs give same result"""
+    im = create_model_3d
     w3d = pointing.create_fitswcs(im)
     gra, gdec, glam = w3d(1, 1, 1)
 
@@ -124,6 +145,8 @@ def test_create_fits_wcs(tmpdir):
     w = wcs.WCS(hdu.header)
     wcel = w.sub(['celestial'])
     ra, dec = wcel.all_pix2world(1, 1, 0)
+
+    # Check that astropy.wcs.WCS and gwcs.WCS give same result
     assert_allclose((ra, dec), (gra, gdec))
 
     # test serialization
diff --git a/jwst/outlier_detection/tests/test_outlier_detection.py b/jwst/outlier_detection/tests/test_outlier_detection.py
index 568c2c3ff..60f3c9fc0 100644
--- a/jwst/outlier_detection/tests/test_outlier_detection.py
+++ b/jwst/outlier_detection/tests/test_outlier_detection.py
@@ -2,26 +2,41 @@ import pytest
 import numpy as np
 from scipy.ndimage.filters import gaussian_filter
 
+from jwst.outlier_detection import OutlierDetectionStep
 from jwst.outlier_detection.outlier_detection import flag_cr
 from jwst import datamodels
+from jwst.assign_wcs.pointing import create_fitswcs
+
+
+OUTLIER_DO_NOT_USE = np.bitwise_or(datamodels.dqflags.pixel["DO_NOT_USE"],
+                                    datamodels.dqflags.pixel["OUTLIER"])
 
 
 @pytest.fixture
 def sci_blot_image_pair():
     """Provide a science and blotted ImageModel pair."""
-    shape = (10, 10)
+    shape = (20, 20)
+    sigma = 0.02
+    background = 3
+
     sci = datamodels.ImageModel(shape)
 
     # Populate keywords
     sci.meta.exposure.exposure_time = 1
+    sci.meta.background.subtracted = False
+    sci.meta.background.level = background
 
-    # Add poisson noise to image data
-    p = np.random.poisson(size=shape, lam=1e3)
-    sci.data = p / p.mean() - 1
+    sci.data = np.random.normal(loc=background, size=shape, scale=sigma)
+    sci.err = np.zeros(shape) + sigma
 
-    # The blot image is just a smoothed version of the science image
+    # Add a source in the center
+    sci.data[10, 10] += 20 * sigma
+
+    # The blot image is just a smoothed version of the science image that has
+    # its background subtracted
     blot = sci.copy()
     blot.data = gaussian_filter(blot.data, sigma=3)
+    blot.data -= background
 
     return sci, blot
 
@@ -31,26 +46,129 @@ def test_flag_cr(sci_blot_image_pair):
     sci, blot = sci_blot_image_pair
     assert (sci.dq == 0).all()
 
-    # Add some background
-    sci.data += 3
-    blot.data += 3
-
-    # Drop a CR on the science array
-    sci.data[5, 5] += 10
+    # Drop some CRs on the science array
+    sci.data[3, 3] += 100
+    sci.data[3, 7] += 1e3
+    sci.data[7, 3] += 1e4
+    sci.data[7, 7] += 1e5
 
+    # run flag_cr() which updates in-place.  Copy sci first.
+    data_copy = sci.data.copy()
     flag_cr(sci, blot)
-    assert sci.dq[5, 5] > 0
 
+    # Make sure science data array is unchanged after flag_cr()
+    np.testing.assert_allclose(sci.data, data_copy)
 
-def test_flag_cr_with_subtracted_background(sci_blot_image_pair):
-    """Test the flag_cr function on background-subtracted data"""
-    sci, blot = sci_blot_image_pair
+    # Verify that both DQ flags are set in the DQ array for all outliers
+    assert sci.dq[3, 3] == OUTLIER_DO_NOT_USE
+    assert sci.dq[3, 7] == OUTLIER_DO_NOT_USE
+    assert sci.dq[7, 3] == OUTLIER_DO_NOT_USE
+    assert sci.dq[7, 7] == OUTLIER_DO_NOT_USE
 
-    sci.meta.background.subtracted = True
-    sci.meta.background.level = 3
+    # Verify the source wasn't flagged
+    assert sci.dq[10, 10] == datamodels.dqflags.pixel["GOOD"]
+
+
+@pytest.fixture
+def we_three_sci():
+    """Provide 3 science images with different noise but identical source
+    and same background level"""
+    shape = (20, 20)
+    sigma = 0.02
+    background = 1.5
+    signal = 200 * sigma
+
+    sci1 = datamodels.ImageModel(shape)
+
+    # Populate keywords
+    sci1.meta.instrument.name = "MIRI"
+    sci1.meta.instrument.detector = "MIRIMAGE"
+    sci1.meta.exposure.type = "MIR_IMAGE"
+    sci1.meta.observation.date = "2020-01-01"
+    sci1.meta.observation.time = "00:00:00"
+    sci1.meta.telescope = "JWST"
+    sci1.meta.exposure.exposure_time = 1
+    sci1.meta.wcsinfo.wcsaxes = 2
+    sci1.meta.wcsinfo.ctype1 = "RA---TAN"
+    sci1.meta.wcsinfo.ctype2 = "DEC--TAN"
+    sci1.meta.wcsinfo.cdelt1 = 3e-6
+    sci1.meta.wcsinfo.cdelt2 = 3e-6
+    sci1.meta.wcsinfo.roll_ref = 0
+    sci1.meta.wcsinfo.v3yangle = 0
+    sci1.meta.wcsinfo.vparity = -1
+    sci1.meta.wcsinfo.pc1_1 = 1
+    sci1.meta.wcsinfo.pc1_2 = 0
+    sci1.meta.wcsinfo.pc2_1 = 0
+    sci1.meta.wcsinfo.pc2_2 = 1
+    sci1.meta.wcsinfo.crpix1 = 5
+    sci1.meta.wcsinfo.crpix2 = 5
+    sci1.meta.wcsinfo.crval1 = 0
+    sci1.meta.wcsinfo.crval2 = 0
+    sci1.meta.wcsinfo.cunit1 = "deg"
+    sci1.meta.wcsinfo.cunit2 = "deg"
+    sci1.meta.background.subtracted = False
+    sci1.meta.background.level = background
+
+    # Replace the FITS-type WCS with an Identity WCS
+    sci1.meta.wcs = create_fitswcs(sci1)
+
+
+    sci1.err = np.zeros(shape) + sigma
+
+    # Make copies with different noise
+    sci2 = sci1.copy()
+    sci3 = sci1.copy()
+
+    # Populate data background with random noise
+    sci1.data = np.random.normal(loc=background, size=shape, scale=sigma)
+    sci2.data = np.random.normal(loc=background, size=shape, scale=sigma)
+    sci3.data = np.random.normal(loc=background, size=shape, scale=sigma)
+
+    sci1.meta.filename = "foo1_cal.fits"
+    sci2.meta.filename = "foo2_cal.fits"
+    sci3.meta.filename = "foo3_cal.fits"
+
+    # Add a source in the centers
+    sci1.data[7, 7] += signal
+    sci2.data[7, 7] += signal
+    sci3.data[7, 7] += signal
+
+    return sci1, sci2, sci3
+
+
+def test_outlier_step_no_outliers(we_three_sci):
+    """Test whole step, no outliers"""
+    container = datamodels.ModelContainer(list(we_three_sci))
+    pristine = container.copy()
+    result = OutlierDetectionStep.call(container)
+
+    # Make sure nothing changed in SCI and DQ arrays
+    for image, uncorrected in zip(pristine, container):
+        np.testing.assert_allclose(image.data, uncorrected.data)
+        np.testing.assert_allclose(image.dq, uncorrected.dq)
+
+    # Make sure nothing changed in SCI and DQ arrays
+    for image, corrected in zip(container, result):
+        np.testing.assert_allclose(image.data, corrected.data)
+        np.testing.assert_allclose(image.dq, corrected.dq)
+
+
+def test_outlier_step(we_three_sci):
+    """Test whole step with an outlier"""
+    container = datamodels.ModelContainer(list(we_three_sci))
 
     # Drop a CR on the science array
-    sci.data[5, 5] += 10
+    container[0].data[12, 12] += 1e3
 
-    flag_cr(sci, blot)
-    assert sci.dq[5, 5] > 0
+    result = OutlierDetectionStep.call(container)
+
+    # Make sure nothing changed in SCI array
+    for image, corrected in zip(container, result):
+        np.testing.assert_allclose(image.data, corrected.data)
+
+    # Verify source is not flagged
+    for r in result:
+        assert r.dq[7, 7] == datamodels.dqflags.pixel["GOOD"]
+
+    # Verify CR is flagged
+    assert result[0].dq[12, 12] == OUTLIER_DO_NOT_USE