Repository: spacetelescope/jwst

Issue:
Outlier detection step skips execution using call() on calints
_Issue [JP-3441](https://jira.stsci.edu/browse/JP-3441) was created on JIRA by [Nikolay Nikolov](https://jira.stsci.edu/secure/ViewProfile.jspa?name=nnikolov):_

I ran the outlier detection step (version 1.12.0) from calwebb_tso3 using the call() command on a single segment calints.fits obtained in the BOTS mode on NIRSpec with the following command:

calwebb_tso3.outlier_detection_step.OutlierDetectionStep.call(svdir_st2all + froot[0] + 'calints.fits', output_dir=svdir_st3all, save_results=True)

The step executed, but produced the following warnings and skipped:

WARNING - Input is not the expected CubeModel
WARNING - Outlier detection step will be skipped

Note: When running the Tso3Pipeline(), and including the outlier detection step with Tso3Pipeline.outlier_detection.skip = False, I can see that the step executes and outputs at the terminal messages such as:

New pixels flagged as outliers: 1 (0.00%)

So the issue must be only when running the step separately. 

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/spacetelescope/jwst/pull/8473?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) Report
Attention: Patch coverage is `50.92593%` with `53 lines` in your changes are missing coverage. Please review.
> Project coverage is 56.32%. Comparing base [(`6580914`)](https://app.codecov.io/gh/spacetelescope/jwst/commit/65809147c8caca44f61cbd7b7e46bd3ed24dbe0e?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) to head [(`7d3739f`)](https://app.codecov.io/gh/spacetelescope/jwst/pull/8473?dropdown=coverage&src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope).
> Report is 13 commits behind head on master.

| [Files](https://app.codecov.io/gh/spacetelescope/jwst/pull/8473?dropdown=coverage&src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) | Patch % | Lines |
|---|---|---|
| [jwst/outlier\_detection/outlier\_detection.py](https://app.codecov.io/gh/spacetelescope/jwst/pull/8473?src=pr&el=tree&filepath=jwst%2Foutlier_detection%2Foutlier_detection.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope#diff-andzdC9vdXRsaWVyX2RldGVjdGlvbi9vdXRsaWVyX2RldGVjdGlvbi5weQ==) | 54.38% | [26 Missing :warning: ](https://app.codecov.io/gh/spacetelescope/jwst/pull/8473?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) |
| [jwst/outlier\_detection/outlier\_detection\_tso.py](https://app.codecov.io/gh/spacetelescope/jwst/pull/8473?src=pr&el=tree&filepath=jwst%2Foutlier_detection%2Foutlier_detection_tso.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope#diff-andzdC9vdXRsaWVyX2RldGVjdGlvbi9vdXRsaWVyX2RldGVjdGlvbl90c28ucHk=) | 44.44% | [20 Missing :warning: ](https://app.codecov.io/gh/spacetelescope/jwst/pull/8473?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) |
| [jwst/outlier\_detection/outlier\_detection\_ifu.py](https://app.codecov.io/gh/spacetelescope/jwst/pull/8473?src=pr&el=tree&filepath=jwst%2Foutlier_detection%2Foutlier_detection_ifu.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope#diff-andzdC9vdXRsaWVyX2RldGVjdGlvbi9vdXRsaWVyX2RldGVjdGlvbl9pZnUucHk=) | 45.45% | [6 Missing :warning: ](https://app.codecov.io/gh/spacetelescope/jwst/pull/8473?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) |
| [jwst/outlier\_detection/outlier\_detection\_spec.py](https://app.codecov.io/gh/spacetelescope/jwst/pull/8473?src=pr&el=tree&filepath=jwst%2Foutlier_detection%2Foutlier_detection_spec.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope#diff-andzdC9vdXRsaWVyX2RldGVjdGlvbi9vdXRsaWVyX2RldGVjdGlvbl9zcGVjLnB5) | 0.00% | [1 Missing :warning: ](https://app.codecov.io/gh/spacetelescope/jwst/pull/8473?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) |

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##           master    #8473      +/-   ##
==========================================
- Coverage   56.38%   56.32%   -0.06%     
==========================================
  Files         387      388       +1     
  Lines       38716    38867     +151     
==========================================
+ Hits        21830    21893      +63     
- Misses      16886    16974      +88     
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/spacetelescope/jwst/pull/8473?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope).

emolter: Regression tests [started here](https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/1436/).

Marking this ready for review even though I still need to figure out why the readthedocs build is failing.
emolter: the only regtest failure was also present [here](https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/1434/testReport/junit/jwst.regtest/test_miri_lrs_slit_spec3/_stable_deps__test_miri_lrs_slit_spec3_user_wcs_shape1_s2d_/), so I don't think this PR causes that issue.

However, I seem to have introduced a new unit test failure with the last small bugfix. Will track that down
emolter: New [regtest run started here](https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/1442/) for code changes that address comment from @braingram.  The imaging mode is still refactored, but the new version should preserve the way memory is used.  Let me know whether the refactor is desired - some of the refactoring is definitely extraneous to the TSO changes and could be reverted if not desired.
braingram: > New [regtest run started here](https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/1442/) for code changes that address comment from @braingram. The imaging mode is still refactored, but the new version should preserve the way memory is used. Let me know whether the refactor is desired - some of the refactoring is definitely extraneous to the TSO changes and could be reverted if not desired.

Thanks for putting this together. One general question I have is: what does the refactoring of the non-tso modes provide? If it's not required for the tso changes I think it makes sense to be a separate PR where the addressed issues and improvements are described.
emolter: Agree with you RE some aspects of the additional refactor, specifically the modelwise_operation decorator.  I started working on it when I misunderstood exactly what was needed for the TSO modes, and then just didn't revert it when I realized the TSO modes needed something different.  I do think the decorators could be a useful way to ensure that memory is taken care of responsibly when running certain operations, and I'd be curious to hear your thoughts on that.  However it's a solution to a problem that doesn't exist yet, and may never exist.

It does feel a bit strange to revert the use of `timewise_operation` in the imaging mode, because it is needed for the TSO data and its use in both places avoids substantial duplication.  But I am happy to do so if that is preferred.

I think splitting the imaging mode into its own class and inheriting from a base class makes sense to include in this PR because imaging, spec, and now TSO all use similar but distinct data processing workflows.  I don't know much about the coronagraphic modes but it's possible that coronagraphy might need its own, very similar, workflow as well - this is something I have been meaning to ask @hbushouse.  
braingram: Thanks for the detailed response!

> Agree with you RE some aspects of the additional refactor, specifically the modelwise_operation decorator. I started working on it when I misunderstood exactly what was needed for the TSO modes, and then just didn't revert it when I realized the TSO modes needed something different. I do think the decorators could be a useful way to ensure that memory is taken care of responsibly when running certain operations, and I'd be curious to hear your thoughts on that. However it's a solution to a problem that doesn't exist yet, and may never exist.

Thanks for the clarification. I'd be happy to discuss this.
 
> It does feel a bit strange to revert the use of `timewise_operation` in the imaging mode, because it is needed for the TSO data and its use in both places avoids substantial duplication. But I am happy to do so if that is preferred.

Perhaps we could discuss this also as I'm still struggling to see how the `get_sections` api helps with the TSO data. As an example for non-tso data if an association containing nircam images is fed into outlier detection (with the default `in_memory=False` and `resample_data=True`) the `drizzled_models` produced by the resampling will be a `ModelContainer` containing filenames. When these are passed to `create_median` the `get_sections` api is used to limit the number of in-memory models to 1 at a time while still being able to take a median across all models. If the input to `create_median` is instead a container of models (not filenames) the function becomes very inefficient where each model is cloned for each section.

Is the input tso data a 3D cube? If so, the entirety of the input will be read on the first call to `datamodels.open`. `calwebb_tso3` does [transform the input to a container of 2d images](https://github.com/spacetelescope/jwst/blob/8b254ae2b69fed42c413f4ff7163ff2f983dfa59/jwst/pipeline/calwebb_tso3.py#L100) before passing them to `outlier_detection` (do you know why it doesn't pass the cube? let's assume that can be changed as part of this improvement). Since the entire cube is in memory (and is not resampled) does it make sense to compute medians on slices from the cube (based on the size of the rolling window). This should be simpler and more efficient than trying to use the `get_sections` api.

> 
> I think splitting the imaging mode into its own class and inheriting from a base class makes sense to include in this PR because imaging, spec, and now TSO all use similar but distinct data processing workflows. I don't know much about the coronagraphic modes but it's possible that coronagraphy might need its own, very similar, workflow as well - this is something I have been meaning to ask @hbushouse.

Splitting up the modes into sub-classes sounds great to me. I'm much less familiar with the non-image paths in this step and have been mostly focused on trying to improve the performance. I certainly do not mean to derail this PR but hopefully some of my comments have been helpful.


emolter: Yet another regtest run [started here.](https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/1451/) 
braingram: Thanks! Please do consider all of my previous comments addressed. The code changes look good to me and feel free to ping me when the regtests finish if a review would be helpful.
emolter: All the regression test failures seem to be unrelated.  This is ready for your review @hbushouse.

You may wonder to yourself why there aren't any expected regtest failures for such a change, which should affect all data processed through calwebb_tso3.  The reason is that there are no regtest datasets with >25 integrations (the default rolling window), which means that as written the code will revert to using a simple median.

Another note: I did check that the unit test I wrote fails for simple median but succeeds for rolling-median.  I think the fake test data are therefore similar to the real test cases where this behavior was seen.
perrygreenfield: I think that all references to the term "drizzle" in `outlier_detection_tso.py` should be removed since there is no drizzling involved. I realize it makes it easier to reuse code from the original source code, but it is confusing to see the odd references. The methods and functions within this file shouldn't use the term. As for the use of the `build_driz_weight` function from `resample_utils.py`, that's unfortunate (I don't think it needed that name in that file either) but you can change the import to rename it to simply `build_weight`
perrygreenfield: So in looking at things, it seems that this aspect should be able to be handled by your PR entirely. Namely, that you change `calwebb_tso3.py` to simply supply a cube model to the `outlier_detection` step and make the necessary changes to the code to deal with that. It should simplify the `outlier_detection_tso.py` code. What I'm not entirely sure about is the `outlier_detection.py` code. Superficially it seems that a small change to the logic should permit as an argument a `ModelCube` instance. But I'm not sure if there are implied association-related attributes that must be dealt with in the rest of the step, but looking at the code, it doesn't seem like it is needed by anything, but I may be missing something. Brett can correct me if I'm wrong.
emolter: > So in looking at things, it seems that this aspect should be able to be handled by your PR entirely. Namely, that you change `calwebb_tso3.py` to simply supply a cube model to the `outlier_detection` step and make the necessary changes to the code to deal with that. It should simplify the `outlier_detection_tso.py` code. What I'm not entirely sure about is the `outlier_detection.py` code. Superficially it seems that a small change to the logic should permit as an argument a `ModelCube` instance. But I'm not sure if there are implied association-related attributes that must be dealt with in the rest of the step, but looking at the code, it doesn't seem like it is needed by anything, but I may be missing something. Brett can correct me if I'm wrong.

@perrygreenfield if I understand what you're suggesting, it's to remove all the logic in `calwebb_tso3` prior to the outlier detection step that builds a ModelContainer.  The following is what is currently in there:
```
            # FIXME OutlierDetectionStep can accept a cube input and will
            # unpack it as a container which looks almost identical to the below
            # code. There are some differences:
            # - the inputs here are SlitModels
            # - outlier detection will generate drizzle weight during the unpack
            #   the below code does not
            # - the code below takes a "bitwise or" whereas outlier detection
            #   assigns over the dq array (I think this is the same as I expect
            #   outlier detection to only add dq flags).

            # Perform regular outlier detection

            input_2dmodels = ModelContainer()
            for i in range(cube.data.shape[0]):
                # convert each plane of data cube into its own array
                image = datamodels.ImageModel(data=cube.data[i],
                                              err=cube.err[i], dq=cube.dq[i])
                image.update(cube)
                image.meta.wcs = cube.meta.wcs
                input_2dmodels.append(image)

            self.log.info("Performing outlier detection on input images ...")
            input_2dmodels = self.outlier_detection(input_2dmodels)

            # Transfer updated DQ values to original input observation
            for i in range(cube.data.shape[0]):
                # Update DQ arrays with those from outlier_detection step
                cube.dq[i] = np.bitwise_or(cube.dq[i], input_2dmodels[i].dq)
```

That seems like a reasonable idea, especially given the long fixme comment.  Then, you are suggesting that the `_convert_inputs()` method be bypassed inside outlier_detection for TSO data only.  Is that right?  If so, then I agree there should be no problems in outlier_detection.py, because we would only be bypassing `_convert_inputs` for TSO data and handling those data completely separately.  I may need to modify or make a new version of the `detect_outliers()` method, which at present I have re-used.

Is everyone comfortable with the differences named here, i.e., that outlier detection still generates a weight even though drizzling is not actually done, and the bitwise_or is removed in favor of handling setting dq flags entirely within the step?
melanieclarke: Just to note - there is a related ticket (JP-3441, #8009) that requests the ability to directly run the outlier detection step on a calints file, outside the pipeline.  If cube handling is moved inside the step, that should incidentally solve the other issue as well.
perrygreenfield: @emolter, yes, I think so. If it can be made to work, it would be much more straightforward, and work better with Brett's new `ModelContainer` replacement. I do worry a bit about unnoticed snags, but something like that should work. Right now it seems to go through a lot of unnecessary contortions to fit the current scheme. There isn't much code in `detect_outliers` anyway and perhaps the blot stuff can just be ignored.
emolter: > I'm to blame for the long comment. I still don't fully understand what's going on but looking at it now I'm confused by the weights :-/ For the [default](https://github.com/spacetelescope/jwst/blob/80dd865164aa8803ee939a31af0ce9a38e459551/jwst/outlier_detection/outlier_detection_step.py#L52) `ivm` weight type this should use the [var_noise](https://github.com/spacetelescope/jwst/blob/80dd865164aa8803ee939a31af0ce9a38e459551/jwst/resample/resample_utils.py#L174). Does this exist for the cubes? If so, is the `var_noise` passed on to the models during the conversion to a `ModelContainer` in `calwebb_tso3` or in `outlier_detection`? If not, it looks like the `wht` is effectively just a mask for `DO_NOT_USE` (1 for usable pixels, 0 for others). I think I must be missing something but will this interfere with the outlier detection?

I can confirm that the weights seem to just be the DQ mask, and the var_noise is not passed into the step at all, neither during `calwebb_tso3` nor within `_convert_inputs()`.  So I suppose it will indeed make a difference to bypass both, because if a CubeModel is passed into `build_driz_weight()` that has `var_noise`, then those weights will get applied.   I would have to think through whether this behavior would be acceptable or correct for outlier detection - maybe someone knows off the top of their head.
perrygreenfield: 
> I can confirm that the weights seem to just be the DQ mask, and the var_noise is not passed into the step at all, neither during `calwebb_tso3` nor within `_convert_inputs()`. So I suppose it will indeed make a difference to bypass both, because if a CubeModel is passed into `build_driz_weight()` that has `var_noise`, then those weights will get applied. I would have to think through whether this behavior would be acceptable or correct for outlier detection - maybe someone knows off the top of their head.

I want to make sure I understand the situation. In past behavior it was just the DQ mask? And with this change, what changes with regard to this, instead it has `var_noise` instead of the DQ mask? In any case, is it possible the old behavior is not desired by the instrument team?

emolter: I don't know what was desired by the instrument teams initially @perrygreenfield, all I can say is it's a bit silly to call `build_driz_weight()` just to return a copy of the DQ mask, but indeed that is what seems to be happening (for the imaging modes too) because `_convert_inputs()` doesn't save the var_noise.  My plan is to simply retain the current behavior, because this ticket has nothing to do with issues with the weights.  But it might be worth following up in the future.
emolter: [new regtests started here](https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/1477/), after adding support for passing CubeModels through outlier_detection_tso
emolter: more regression tests [started here](https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/1481/).  Hopefully accidental file name changes to crfints files have been reverted.
emolter: [Another regtest run](https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/1484/) that incorporates my attempted file name issue fix
emolter: I believe the remaining regtest failures in `test_miri_lrs_slitless` are expected: the TSO3 pipeline is taking in a .crfints file with 10 integrations, and with this code change it should also get a CubeModel with 10 integrations out of outlier detection.
hbushouse: @emolter Latest CI run had lots of failures due to the `n_ints` name needing updating in one or more unit tests.

Patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index f6895f620..0c2ba7d03 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -105,6 +105,10 @@ outlier_detection
 
 - Added tests for changes made in #8464. [#8481]
 
+- Added the option to use a rolling median instead of a simple median
+  to detect outliers in TSO data, with user-defined
+  rolling window width via the ``rolling_window_width`` parameter. [#8473]
+
 photom
 ------
 
diff --git a/docs/jwst/outlier_detection/arguments.rst b/docs/jwst/outlier_detection/arguments.rst
index a52821825..a2900bcbd 100644
--- a/docs/jwst/outlier_detection/arguments.rst
+++ b/docs/jwst/outlier_detection/arguments.rst
@@ -57,6 +57,10 @@ that control the behavior of the processing:
 ``--backg`` (float, default=0.0)
   User-specified background value to apply to the median image.
 
+``--rolling_window_width`` (int, default=25)
+  Number of integrations over which to take the median when using rolling-window
+  median for TSO observations.
+
 ``--save_intermediate_results`` (boolean, default=False)
   Specifies whether or not to save any intermediate products created
   during step processing.
diff --git a/docs/jwst/outlier_detection/index.rst b/docs/jwst/outlier_detection/index.rst
index b32398e33..b2f76757a 100644
--- a/docs/jwst/outlier_detection/index.rst
+++ b/docs/jwst/outlier_detection/index.rst
@@ -10,7 +10,8 @@ Outlier Detection
    main.rst
    arguments.rst
    outlier_detection_step.rst
-   outlier_detection.rst
+   outlier_detection_imaging.rst
+   outlier_detection_tso.rst
    outlier_detection_ifu.rst
    outlier_detection_spec.rst
 
diff --git a/docs/jwst/outlier_detection/outlier_detection.rst b/docs/jwst/outlier_detection/outlier_detection_imaging.rst
similarity index 90%
rename from docs/jwst/outlier_detection/outlier_detection.rst
rename to docs/jwst/outlier_detection/outlier_detection_imaging.rst
index 81fa2faa0..c4cf5eceb 100644
--- a/docs/jwst/outlier_detection/outlier_detection.rst
+++ b/docs/jwst/outlier_detection/outlier_detection_imaging.rst
@@ -171,37 +171,19 @@ These changes result in a minimum amount of memory usage during processing at th
 expense of reading and writing the products from disk.
 
 
+Outlier Detection for Coronagraphic Data
+----------------------------------------
+Coronagraphic data is processed in a near-identical manner to direct imaging data, but
+no resampling occurs.
+
+
 Outlier Detection for TSO data
 -------------------------------
-Time-series observations (TSO) result in input data stored as a 3D CubeModel
-where each plane in the cube represents a separate integration without changing the
-pointing.  Normal imaging data benefit from combining all integrations into a
+Normal imaging data benefit from combining all integrations into a
 single image. TSO data's value, however, comes from looking for variations from one
 integration to the next.  The outlier detection algorithm, therefore, gets run with 
-a few variations to accomodate the nature of these 3D data.
-
-#. Input data is converted from a CubeModel (3D data array) to a ModelContainer
-
-   - Each plane in the original input CubeModel gets copied to a separate model
-     in the ModelContainer
-
-#. The median image is created without resampling the input data
-
-   - All integrations are aligned already, so no resampling or shifting needs to be performed
-  
-#. A matched median gets created by combining the single median frame with the 
-   noise model for each input integration.
-
-#. Perform statistical comparison between the matched median with each input integration.  
-
-#. Update input data model DQ arrays with the mask of detected outliers.
-
-
-.. note:: 
-
-  This same set of steps also gets used to perform outlier detection on
-  coronographic data, because it too is processed as 3D (per-integration)
-  cubes.
+a few variations to accomodate the nature of these 3D data. See the 
+:ref:`TSO outlier detection <outlier-detection-tso>` documentation for details.
 
 
 Outlier Detection for IFU data
@@ -214,4 +196,9 @@ See the :ref:`IFU outlier detection <outlier-detection-ifu>` documentation for
 all the details.
 
 
+Outlier Detection for Slit data
+-------------------------------
+See the :ref:`IFU outlier detection <outlier-detection-spec>` documentation for
+details.
+
 .. automodapi:: jwst.outlier_detection.outlier_detection
diff --git a/docs/jwst/outlier_detection/outlier_detection_step.rst b/docs/jwst/outlier_detection/outlier_detection_step.rst
index 3d3338f9a..3f8b61302 100644
--- a/docs/jwst/outlier_detection/outlier_detection_step.rst
+++ b/docs/jwst/outlier_detection/outlier_detection_step.rst
@@ -35,8 +35,8 @@ input data:
      use :py:class:`~jwst.outlier_detection.outlier_detection.OutlierDetection` with
      resampling turned off as described in :ref:`outlier-detection-imaging`
    - **Time-Series Observations(TSO)**: both imaging and spectroscopic modes, use
-     :py:class:`~jwst.outlier_detection.outlier_detection.OutlierDetection` with resampling
-     turned off as described in :ref:`outlier-detection-imaging`
+     :py:class:`~jwst.outlier_detection.outlier_detection.OutlierDetectionTSO`
+     as described in :ref:`outlier-detection-tso`
    - **IFU observations**: use
      :py:class:`~jwst.outlier_detection.outlier_detection_ifu.OutlierDetectionIFU` as
      described in :ref:`outlier-detection-ifu`
diff --git a/docs/jwst/outlier_detection/outlier_detection_tso.rst b/docs/jwst/outlier_detection/outlier_detection_tso.rst
new file mode 100644
index 000000000..2f74611e3
--- /dev/null
+++ b/docs/jwst/outlier_detection/outlier_detection_tso.rst
@@ -0,0 +1,41 @@
+.. _outlier-detection-tso:
+
+Outlier Detection for TSO Data
+==============================
+
+Time-series observations (TSO) result in input data stored as a 3D CubeModel
+where each plane in the cube represents a separate integration without changing the
+pointing.  Normal imaging data benefit from combining all integrations into a
+single image. TSO data's value, however, comes from looking for variations from one
+integration to the next.  The outlier detection algorithm, therefore, gets run with 
+a few variations to accomodate the nature of these 3D data.
+
+#. Input data is converted into a CubeModel (3D data array) if a ModelContainer
+   of 2D data arrays is provided.
+
+#. The median image is created without resampling the input data
+
+   * The median image is created by combining all planes in the 
+     CubeModel pixel-by-pixel using a rolling-median algorithm, in order
+     to flag outliers integration-by-integration but preserve real time variability.
+   * The ``rolling_window_width`` parameter specifies the number of integrations over
+     which to compute the median. The default is 25. If the number of integrations
+     is less than or equal to ``rolling_window_width``, a simple median is used instead.
+   * The ``nlow`` and ``nhigh`` parameters specify how many low and high values
+     to ignore when computing the median for any given pixel.
+   * The ``maskpt`` parameter sets the percentage of the weight image values to
+     use, and any pixel with a weight below this value gets flagged as "bad" and
+     ignored when the median is taken.
+   * The rolling-median CubeModel (3D data array) is written out to disk as `_<asn_id>_median.fits`
+     if the ``save_intermediate_results`` parameter is set to True.
+   * All integrations are aligned already, so no resampling or shifting needs to be performed
+
+#. A matched median gets created by combining the single median frame with the 
+   noise model for each input integration.
+
+#. A statistical comparison is performed between the matched median with each input integration.  
+
+#. The input data model DQ arrays are updated with the mask of detected outliers.
+
+
+.. automodapi:: jwst.outlier_detection.outlier_detection_tso
\ No newline at end of file
diff --git a/jwst/outlier_detection/outlier_detection.py b/jwst/outlier_detection/outlier_detection.py
index 2a026d8e3..a44a89433 100644
--- a/jwst/outlier_detection/outlier_detection.py
+++ b/jwst/outlier_detection/outlier_detection.py
@@ -165,20 +165,13 @@ class OutlierDetection:
 
         # Perform median combination on set of drizzled mosaics
         median_model.data = self.create_median(drizzled_models)
-        if pars['save_intermediate_results']:
-            if self.outlierpars.get('asn_id', None) is None:
-                suffix_to_remove = self.resample_suffix
-            else:
-                suffix_to_remove = f"_{self.outlierpars['asn_id']}{self.resample_suffix}"
-            median_model_output_path = self.make_output_path(
-                basepath=median_model.meta.filename.replace(suffix_to_remove, '.fits'),
-                suffix='median')
-            median_model.save(median_model_output_path)
-            log.info(f"Saved model in {median_model_output_path}")
+
+        if self.outlierpars['save_intermediate_results']:
+            self.save_median(median_model) 
         else:
             # since we're not saving intermediate results if the drizzled models
             # were written to disk, remove them
-            if not pars['in_memory']:
+            if not self.outlierpars['in_memory']:
                 for fn in drizzled_models._models:
                     _remove_file(fn)
 
@@ -203,6 +196,7 @@ class OutlierDetection:
             for fn in blot_models._models:
                 _remove_file(fn)
         del median_model, blot_models
+    
 
     def create_median(self, resampled_models):
         """Create a median image from the singly resampled images.
@@ -218,35 +212,7 @@ class OutlierDetection:
 
         log.info("Computing median")
 
-        # Compute weight means without keeping datamodels for eacn input open
-        # Start by insuring that the ModelContainer does NOT open and keep each datamodel
-        ropen_orig = resampled_models._return_open
-        resampled_models._return_open = False  # turn off auto-opening of models
-        # keep track of resulting computation for each input resampled datamodel
-        weight_thresholds = []
-        # For each model, compute the bad-pixel threshold from the weight arrays
-        for resampled in resampled_models:
-            m = datamodel_open(resampled)
-            weight = m.wht
-            # necessary in order to assure that mask gets applied correctly
-            if hasattr(weight, '_mask'):
-                del weight._mask
-            mask_zero_weight = np.equal(weight, 0.)
-            mask_nans = np.isnan(weight)
-            # Combine the masks
-            weight_masked = np.ma.array(weight, mask=np.logical_or(
-                mask_zero_weight, mask_nans))
-            # Sigma-clip the unmasked data
-            weight_masked = sigma_clip(weight_masked, sigma=3, maxiters=5)
-            mean_weight = np.mean(weight_masked)
-            # Mask pixels where weight falls below maskpt percent
-            weight_threshold = mean_weight * maskpt
-            weight_thresholds.append(weight_threshold)
-            # close and delete the model, just to explicitly try to keep the memory as clean as possible
-            m.close()
-            del m
-        # Reset ModelContainer attribute to original value
-        resampled_models._return_open = ropen_orig
+        weight_thresholds = self.compute_weight_threshold_container(resampled_models, maskpt)
 
         # Now, set up buffered access to all input models
         resampled_models.set_buffer(1.0)  # Set buffer at 1Mb
@@ -282,6 +248,98 @@ class OutlierDetection:
             del resampled_sci, resampled_weight
 
         return median_image
+    
+
+    def save_median(self, median_model):
+        '''
+        Save median if requested by user
+
+        Parameters
+        ----------
+        median_model : ~jwst.datamodels.ImageModel
+            The median ImageModel or CubeModel to save
+        '''
+        
+        if self.outlierpars.get('asn_id', None) is None:
+            suffix_to_remove = self.resample_suffix
+        else:
+            suffix_to_remove = f"_{self.outlierpars['asn_id']}{self.resample_suffix}"
+        median_model_output_path = self.make_output_path(
+            basepath=median_model.meta.filename.replace(suffix_to_remove, '.fits'),
+            suffix='median')
+        median_model.save(median_model_output_path)
+        log.info(f"Saved model in {median_model_output_path}")
+
+
+    def compute_weight_threshold_container(self, resampled_models, maskpt):
+        '''
+        Compute weight means without keeping datamodels for each input open
+
+        Parameters
+        ----------
+        resampled_models : ~jwst.datamodels.ModelContainer
+            The input data models.
+
+        maskpt : float
+            The percentage of the mean weight to use as a threshold for masking.
+
+        Returns
+        -------
+        list
+            The weight thresholds for each integration.
+        '''
+
+        # Start by ensuring that the ModelContainer does NOT open and keep each datamodel
+        ropen_orig = resampled_models._return_open
+        resampled_models._return_open = True 
+        # keep track of resulting computation for each input resampled datamodel
+        weight_thresholds = []
+        # For each model, compute the bad-pixel threshold from the weight arrays
+        for resampled in resampled_models:
+            weight = resampled.wht
+            weight_threshold = self.compute_weight_threshold(weight, maskpt)
+            weight_thresholds.append(weight_threshold)
+            # close and delete the model, just to explicitly try to keep the memory as clean as possible
+            resampled.close()
+            del resampled
+        # Reset ModelContainer attribute to original value
+        resampled_models._return_open = ropen_orig
+        return weight_thresholds
+    
+
+    def compute_weight_threshold(self, weight, maskpt):
+        '''
+        Compute the weight threshold for a single image or cube.
+
+        Parameters
+        ----------
+        weight : numpy.ndarray
+            The weight array
+
+        maskpt : float
+            The percentage of the mean weight to use as a threshold for masking.
+
+        Returns
+        -------
+        float
+            The weight threshold for this integration.
+        '''
+
+        # necessary in order to assure that mask gets applied correctly
+        if hasattr(weight, '_mask'):
+            del weight._mask
+        mask_zero_weight = np.equal(weight, 0.)
+        mask_nans = np.isnan(weight)
+        # Combine the masks
+        weight_masked = np.ma.array(weight, mask=np.logical_or(
+            mask_zero_weight, mask_nans))
+        # Sigma-clip the unmasked data
+        weight_masked = sigma_clip(weight_masked, sigma=3, maxiters=5)
+        mean_weight = np.mean(weight_masked)
+        # Mask pixels where weight falls below maskpt percent
+        weight_threshold = mean_weight * maskpt
+        return weight_threshold
+
 
     def blot_median(self, median_model):
         """Blot resampled median image back to the detector images."""
@@ -364,10 +422,13 @@ def flag_cr(sci_image, blot_image, snr="5.0 4.0", scale="1.2 0.7", backg=0,
     Parameters
     ----------
     sci_image : ~jwst.datamodels.ImageModel
-        the science data
+        the science data. Can also accept a CubeModel, but only if
+        resample_data is False
 
     blot_image : ~jwst.datamodels.ImageModel
-        the blotted median image of the dithered science frames
+        the blotted median image of the dithered science frames.
+        Can also accept a CubeModel, but only if resample_data
+        is False
 
     snr : str
         Signal-to-noise ratio
@@ -381,6 +442,14 @@ def flag_cr(sci_image, blot_image, snr="5.0 4.0", scale="1.2 0.7", backg=0,
     resample_data : bool
         Boolean to indicate whether blot_image is created from resampled,
         dithered data or not
+
+    Notes
+    -----
+    Accepting a CubeModel for sci_image and blot_image with resample_data=True
+    appears to be a relatively simple extension, as the only thing that explicitly
+    relies on the dimensionality is the kernel, which could be generalized.
+    However, this is not currently needed, as CubeModels are only passed in for
+    TSO data, where resampling is always False.
     """
     snr1, snr2 = [float(val) for val in snr.split()]
     scale1, scale2 = [float(val) for val in scale.split()]
@@ -440,7 +509,7 @@ def flag_cr(sci_image, blot_image, snr="5.0 4.0", scale="1.2 0.7", backg=0,
     # Report number (and percent) of new DO_NOT_USE pixels found
     count_outlier = np.count_nonzero(sci_image.dq & DO_NOT_USE)
     count_added = count_outlier - count_existing
-    percent_cr = count_added / (sci_image.shape[0] * sci_image.shape[1]) * 100
+    percent_cr = count_added / sci_image.dq.size * 100
     log.info(f"New pixels flagged as outliers: {count_added} ({percent_cr:.2f}%)")
 
 
diff --git a/jwst/outlier_detection/outlier_detection_step.py b/jwst/outlier_detection/outlier_detection_step.py
index b7798e1aa..18140b1c6 100644
--- a/jwst/outlier_detection/outlier_detection_step.py
+++ b/jwst/outlier_detection/outlier_detection_step.py
@@ -10,11 +10,13 @@ from jwst.lib.pipe_utils import is_tso
 from jwst.outlier_detection import outlier_detection
 from jwst.outlier_detection import outlier_detection_ifu
 from jwst.outlier_detection import outlier_detection_spec
+from jwst.outlier_detection import outlier_detection_tso
 
 # Categorize all supported versions of outlier_detection
 outlier_registry = {'imaging': outlier_detection.OutlierDetection,
                     'ifu': outlier_detection_ifu.OutlierDetectionIFU,
-                    'slitspec': outlier_detection_spec.OutlierDetectionSpec
+                    'slitspec': outlier_detection_spec.OutlierDetectionSpec,
+                    'tso': outlier_detection_tso.OutlierDetectionTSO
                     }
 
 # Categorize all supported modes
@@ -61,6 +63,7 @@ class OutlierDetectionStep(Step):
         backg = float(default=0.0)
         kernel_size = string(default='7 7')
         threshold_percent = float(default=99.8)
+        rolling_window_width = integer(default=25)
         ifu_second_check = boolean(default=False)
         save_intermediate_results = boolean(default=False)
         resample_data = boolean(default=True)
@@ -112,6 +115,7 @@ class OutlierDetectionStep(Step):
                 'backg': self.backg,
                 'kernel_size': self.kernel_size,
                 'threshold_percent': self.threshold_percent,
+                'rolling_window_width': self.rolling_window_width,
                 'ifu_second_check': self.ifu_second_check,
                 'allowed_memory': self.allowed_memory,
                 'in_memory': self.in_memory,
@@ -121,7 +125,7 @@ class OutlierDetectionStep(Step):
                 'make_output_path': self.make_output_path,
             }
 
-            # Add logic here to select which version of OutlierDetection
+            # Select which version of OutlierDetection
             # needs to be used depending on the input data
             if self.input_container:
                 single_model = self.input_models[0]
@@ -130,10 +134,13 @@ class OutlierDetectionStep(Step):
             exptype = single_model.meta.exposure.type
             self.check_input()
 
-            # check for TSO models first
-            if is_tso(single_model) or exptype in CORON_IMAGE_MODES:
-                # algorithm selected for TSO data (no resampling)
-                pars['resample_data'] = False  # force resampling off...
+            if is_tso(single_model):
+                # force resampling off and use rolling median
+                pars['resample_data'] = False
+                detection_step = outlier_registry['tso']
+            elif exptype in CORON_IMAGE_MODES:
+                # force resampling off but use same workflow as imaging
+                pars['resample_data'] = False
                 detection_step = outlier_registry['imaging']
             elif exptype in IMAGE_MODES:
                 # imaging with resampling
@@ -173,6 +180,7 @@ class OutlierDetectionStep(Step):
                 self.input_models.meta.cal_step.outlier_detection = state
             return self.input_models
 
+
     def check_input(self):
         """Use this method to determine whether input is valid or not."""
         if self.input_container:
@@ -198,7 +206,7 @@ class OutlierDetectionStep(Step):
     def _check_input_cube(self):
         """Check to see whether input is the expected CubeModel object."""
         ninputs = self.input_models.shape[0]
-        if not isinstance(self.input_models, datamodels.CubeModel):
+        if type(self.input_models) not in [datamodels.CubeModel, datamodels.SlitModel]:
             self.log.warning("Input is not the expected CubeModel")
             self.log.warning("Outlier detection step will be skipped")
             self.valid_input = False
diff --git a/jwst/outlier_detection/outlier_detection_tso.py b/jwst/outlier_detection/outlier_detection_tso.py
new file mode 100644
index 000000000..4333adbbb
--- /dev/null
+++ b/jwst/outlier_detection/outlier_detection_tso.py
@@ -0,0 +1,172 @@
+import numpy as np
+from .outlier_detection import OutlierDetection, flag_cr
+from jwst.resample.resample_utils import build_mask
+
+from jwst import datamodels as dm
+import logging
+log = logging.getLogger(__name__)
+log.setLevel(logging.DEBUG)
+
+
+__all__ = ["OutlierDetectionTSO"]
+
+
+class OutlierDetectionTSO(OutlierDetection):
+    """Class to flag outlier pixels in DQ of TSO data. Works similarly to
+    imaging outlier detection, but does not resample and uses a rolling median."""
+
+    def __init__(self, input_model: dm.SlitModel | dm.CubeModel | dm.ModelContainer, **pars):
+        """Initialize class for TSO data processing.
+
+        Parameters
+        ----------
+        input_models : ~jwst.datamodels.CubeModel or ~jwst.datamodels.ModelContainer
+            The input TSO data cube. if a ModelContainer is passed in, it is assumed
+            to be a list of ImageModels, and gets converted to a CubeModel.
+
+        """
+        super().__init__(input_model, **pars)
+        if isinstance(self.inputs, dm.ModelContainer):
+            raise TypeError("OutlierDetectionTSO does not support ModelContainer input.")
+        
+        self.input_models = self.inputs
+        self.converted = False
+
+
+    def do_detection(self):
+        """Flag outlier pixels in DQ of input images."""
+        weighted_cube = self.weight_no_resample()
+
+        maskpt = self.outlierpars.get('maskpt', 0.7)
+        weight_threshold = self.compute_weight_threshold(weighted_cube.wht, maskpt)
+
+        rolling_width = self.outlierpars.get("rolling_window_width", 25)
+        if (rolling_width > 1) and (rolling_width < weighted_cube.shape[0]):
+            medians = self.compute_rolling_median(weighted_cube, weight_threshold, w=rolling_width)
+
+        else:
+            medians = np.nanmedian(weighted_cube.data, axis=0)
+            # this is a 2-D array, need to repeat it into the time axis
+            # for consistent shape with rolling median case
+            medians = np.broadcast_to(medians, weighted_cube.shape)
+
+        # turn this into a cube model for saving and passing to flag_cr
+        median_model = dm.CubeModel(data=medians)
+
+        # Save median model if pars['save_intermediate_results'] is True
+        # this will be a CubeModel with rolling median values.
+        if self.outlierpars['save_intermediate_results']:
+            with dm.open(weighted_cube[0]) as dm0:
+                median_model.update(dm0)
+            self.save_median(median_model)
+
+        # no need for blotting, resample is turned off for TSO
+        # go straight to outlier detection
+        log.info("Flagging outliers")
+        flag_cr(self.input_models, median_model, **self.outlierpars)
+
+
+    def weight_no_resample(self):
+        """
+        give weights to model without resampling
+
+        Notes
+        -----
+        Prior to PR #8473, the `build_driz_weight` function was used to
+        create the weights for the input models for TSO data. However, that
+        function was simply returning a copy of the DQ array because the
+        var_noise was not being passed in by calwebb_tso3. As of PR #8473,
+        a cube model that includes the var_noise is passed into TSO
+        outlier detection, so `build_driz_weight` would weight the cube model
+        by the variance. Therefore `build_driz_weight` was removed in order to
+        preserve the original behavior. If it is determined later that exposure
+        time or inverse variance weighting should be used here, build_driz_weight
+        should be re-implemented.
+        """
+        weighted_cube = self.input_models.copy()
+        dqmask = build_mask(self.input_models.dq, self.outlierpars['good_bits'])
+        weighted_cube.wht = dqmask.astype(np.float32)
+        return weighted_cube
+
+
+    def compute_rolling_median(self, model: dm.CubeModel, weight_threshold: np.ndarray, w: int=25) -> np.ndarray:
+        '''
+        Set bad and low-weight data to NaN, then compute the rolling median over the time axis.
+
+        Parameters
+        ----------
+        model : ~jwst.datamodels.CubeModel
+            The input cube model
+
+        weight_threshold : np.ndarray
+            The weight thresholds for each integration.
+
+        w : int
+            The window size for the rolling median.
+
+        Returns
+        -------
+        np.ndarray
+            The rolling median of the input data. Same dimensions as input.
+        '''
+
+        sci = model.data
+        weight = model.wht
+        badmask = np.less(weight, weight_threshold)
+        log.debug("Percentage of pixels with low weight: {}".format(
+            np.sum(badmask) / weight.size * 100))
+
+        # Fill resampled_sci array with nan's where mask values are True
+        sci[badmask] = np.nan
+        del badmask
+
+        if w > sci.shape[0]:
+            raise ValueError("Window size must be less than the number of integrations.")
+        meds = moving_median_over_zeroth_axis(sci, w)
+
+        del sci
+        return meds
+
+
+def moving_median_over_zeroth_axis(x: np.ndarray, w: int) -> np.ndarray:
+    """
+    Calculate the median of a moving window over the zeroth axis of an N-d array.
+    Algorithm works by expanding the array into an additional dimension
+    where the new axis has the same length as the window size. Each entry in that
+    axis is a copy of the original array shifted by 1 with respect to the previous
+    entry, such that the rolling median is simply the median over the new axis.
+    modified from https://stackoverflow.com/a/71154394, see link for more details.
+
+    Parameters
+    ----------
+    x : np.ndarray
+        The input array.
+
+    w : int
+        The window size.
+
+    Returns
+    -------
+    np.ndarray
+        The rolling median of the input array. Same dimensions as input.
+
+    Notes
+    -----
+    Returns NaN for the first and last w//2 elements.
+    """
+    if w <= 1:
+        raise ValueError("Rolling median window size must be greater than 1.")
+    shifted = np.zeros((x.shape[0] + w - 1, w, *x.shape[1:])) * np.nan
+    for idx in range(w - 1):
+        shifted[idx:-w + idx + 1, idx] = x
+    shifted[idx + 1:, idx + 1] = x
+    medians = np.median(shifted, axis=1)
+    for idx in range(w - 1):
+        medians[idx] = np.median(shifted[idx, :idx + 1])
+        medians[-idx - 1] = np.median(shifted[-idx - 1, -idx - 1:])
+    medians = medians[(w - 1) // 2:-(w - 1) // 2]
+
+    # Fill in the edges with the nearest valid value
+    medians[:w // 2] = medians[w // 2]
+    medians[-w // 2:] = medians[-w // 2]
+    return medians
diff --git a/jwst/pipeline/calwebb_tso3.py b/jwst/pipeline/calwebb_tso3.py
index 3c38265b9..91ad8b505 100644
--- a/jwst/pipeline/calwebb_tso3.py
+++ b/jwst/pipeline/calwebb_tso3.py
@@ -5,7 +5,6 @@ from astropy.table import vstack
 
 from stdatamodels.jwst import datamodels
 
-from jwst.datamodels import ModelContainer
 
 from ..stpipe import Pipeline
 
@@ -78,7 +77,7 @@ class Tso3Pipeline(Pipeline):
 
         # Input may consist of multiple exposures, so loop over each of them
         input_exptype = None
-        for cube in input_models:
+        for i, cube in enumerate(input_models):
             if input_exptype is None:
                 input_exptype = cube.meta.exposure.type
 
@@ -87,51 +86,26 @@ class Tso3Pipeline(Pipeline):
                 self.log.warning('Input data are 2D; skipping outlier_detection')
                 break
 
-            # FIXME OutlierDetectionStep can accept a cube input and will
-            # unpack it as a container which looks almost identical to the below
-            # code. There are some differences:
-            # - the inputs here are SlitModels
-            # - outlier detection will generate drizzle weight during the unpack
-            #   the below code does not
-            # - the code below takes a "bitwise or" whereas outlier detection
-            #   assigns over the dq array (I think this is the same as I expect
-            #   outlier detection to only add dq flags).
-
-            # Perform regular outlier detection
-
-            input_2dmodels = ModelContainer()
-            for i in range(cube.data.shape[0]):
-                # convert each plane of data cube into its own array
-                image = datamodels.ImageModel(data=cube.data[i],
-                                              err=cube.err[i], dq=cube.dq[i])
-                image.update(cube)
-                image.meta.wcs = cube.meta.wcs
-                input_2dmodels.append(image)
-
             self.log.info("Performing outlier detection on input images ...")
-            input_2dmodels = self.outlier_detection(input_2dmodels)
-
-            # Transfer updated DQ values to original input observation
-            for i in range(cube.data.shape[0]):
-                # Update DQ arrays with those from outlier_detection step
-                cube.dq[i] = np.bitwise_or(cube.dq[i], input_2dmodels[i].dq)
-
-            cube.meta.cal_step.outlier_detection = \
-                input_2dmodels[0].meta.cal_step.outlier_detection
+            cube = self.outlier_detection(cube)
 
-            del input_2dmodels
-
-        # Save crfints products
-        if input_models[0].meta.cal_step.outlier_detection == 'COMPLETE':
-            self.log.info("Saving crfints products with updated DQ arrays ...")
-            for cube in input_models:
+            # Save crfints products
+            if cube.meta.cal_step.outlier_detection == 'COMPLETE':
+                self.log.info("Saving crfints products with updated DQ arrays ...")
                 # preserve output filename
                 original_filename = cube.meta.filename
+
+                # ensure output filename will not have duplicate asn_id
+                if "_"+input_models.meta.asn_table.asn_id in original_filename:
+                    original_filename = original_filename.replace(
+                        "_"+input_models.meta.asn_table.asn_id, ''
+                    )
                 self.save_model(
                     cube, output_file=original_filename, suffix='crfints',
                     asn_id=input_models.meta.asn_table.asn_id
                 )
                 cube.meta.filename = original_filename
+            input_models[i] = cube
 
         # Create final photometry results as a single output
         # regardless of how many input members there may be


Test patch:
diff --git a/jwst/outlier_detection/tests/test_medfilt.py b/jwst/outlier_detection/tests/test_algorithms.py
similarity index 68%
rename from jwst/outlier_detection/tests/test_medfilt.py
rename to jwst/outlier_detection/tests/test_algorithms.py
index 13cbd7822..932b66026 100644
--- a/jwst/outlier_detection/tests/test_medfilt.py
+++ b/jwst/outlier_detection/tests/test_algorithms.py
@@ -5,7 +5,7 @@ import numpy as np
 import scipy.signal
 
 from jwst.outlier_detection.outlier_detection_ifu import medfilt
-
+from jwst.outlier_detection.outlier_detection_tso import moving_median_over_zeroth_axis
 
 
 @pytest.mark.parametrize("shape,kern_size", [
@@ -40,3 +40,16 @@ def test_medfilt_nan(arr, kern_size, expected):
         )
         result = medfilt(arr, kern_size)
     np.testing.assert_allclose(result, expected)
+
+
+def test_rolling_median():
+
+    time_axis = np.arange(10)
+    expected_time_axis = np.array([1, 1, 2, 3, 4, 5, 6, 7, 8, 8])
+    spatial_axis = np.ones((5, 5))
+    arr = time_axis[:, np.newaxis, np.newaxis] * spatial_axis[np.newaxis, :, :]
+
+    w = 3
+    result = moving_median_over_zeroth_axis(arr, w)
+    expected = expected_time_axis[:, np.newaxis, np.newaxis] * spatial_axis[np.newaxis, :, :]
+    assert np.allclose(result, expected)
diff --git a/jwst/outlier_detection/tests/test_outlier_detection.py b/jwst/outlier_detection/tests/test_outlier_detection.py
index 00ea0502a..7757ce61d 100644
--- a/jwst/outlier_detection/tests/test_outlier_detection.py
+++ b/jwst/outlier_detection/tests/test_outlier_detection.py
@@ -148,7 +148,7 @@ def we_many_sci(
         tsci.data = np.random.normal(loc=background, size=shape, scale=sigma)
         # Add a source in the center
         tsci.data[7, 7] += signal
-        tsci.meta.filename = f"foo{i+2}_cal.fits"
+        tsci.meta.filename = f"foo{i + 2}_cal.fits"
         all_sci.append(tsci)
 
     return all_sci
@@ -300,9 +300,9 @@ def test_outlier_step_image_weak_CR_dither(exptype, tmp_cwd):
     assert result[0].dq[12, 12] == OUTLIER_DO_NOT_USE
 
 
-@pytest.mark.parametrize("exptype, tsovisit", exptypes_tso + exptypes_coron)
-def test_outlier_step_image_weak_CR_nodither(exptype, tsovisit, tmp_cwd):
-    """Test whole step with an outlier for TSO & coronagraphic modes"""
+@pytest.mark.parametrize("exptype, tsovisit", exptypes_coron)
+def test_outlier_step_image_weak_CR_coron(exptype, tsovisit, tmp_cwd):
+    """Test whole step with an outlier for coronagraphic modes"""
     bkg = 1.5
     sig = 0.02
     container = ModelContainer(
@@ -327,3 +327,48 @@ def test_outlier_step_image_weak_CR_nodither(exptype, tsovisit, tmp_cwd):
 
     # Verify CR is flagged
     assert result[0].dq[12, 12] == OUTLIER_DO_NOT_USE
+
+
+@pytest.mark.parametrize("exptype, tsovisit", exptypes_tso)
+def test_outlier_step_weak_cr_tso(exptype, tsovisit):
+    '''Test outlier detection with rolling median on time-varying source
+    This test fails if rolling_window_width is set to 100, i.e., take simple median
+    '''
+    bkg = 1.5
+    sig = 0.02
+    rolling_window_width = 7
+    numsci = 50
+    signal = 7.0
+    im = we_many_sci(
+        numsci=numsci, background=bkg, sigma=sig, signal=signal, exptype=exptype, tsovisit=tsovisit
+    )
+    cube_data = np.array([i.data for i in im])
+    cube_err = np.array([i.err for i in im])
+    cube_dq = np.array([i.dq for i in im])
+    cube_var_noise = np.array([i.var_rnoise for i in im])
+    cube = datamodels.CubeModel(data=cube_data, err=cube_err, dq=cube_dq, var_noise=cube_var_noise)
+
+    # update metadata of cube to match the first image
+    cube.meta = im[0].meta
+
+    # Drop a weak CR on the science array
+    cr_timestep = 5
+    cube.data[cr_timestep, 12, 12] = bkg + sig * 10
+
+    # make time variability that has larger total amplitude than
+    # the CR signal but deviations frame-by-frame are smaller
+    real_time_variability = signal * np.cos(np.linspace(0, np.pi, numsci))
+    for i, model in enumerate(im):
+        model.data[7, 7] += real_time_variability[i]
+        model.err[7, 7] = np.sqrt(sig ** 2 + model.data[7, 7])
+
+    result = OutlierDetectionStep.call(cube, rolling_window_width=rolling_window_width)
+
+    # Make sure nothing changed in SCI array
+    np.testing.assert_allclose(cube.data, result.data)
+
+    # Verify source is not flagged
+    assert np.all(result.dq[:, 7, 7] == datamodels.dqflags.pixel["GOOD"])
+
+    # Verify CR is flagged
+    assert result.dq[cr_timestep, 12, 12] == OUTLIER_DO_NOT_USE