diff --git a/src/stcal/outlier_detection/utils.py b/src/stcal/outlier_detection/utils.py
index 2e55fc128..abebd0136 100644
--- a/src/stcal/outlier_detection/utils.py
+++ b/src/stcal/outlier_detection/utils.py
@@ -73,32 +73,7 @@ def compute_weight_threshold(weight, maskpt):
     return weight_threshold
 
 
-def _valid_abs_deriv(array):
-    """Take the absolute derivate of a numpy array."""
-    out = np.zeros_like(array)  # use same dtype as input
-
-    # compute row-wise absolute diffference
-    d = np.abs(np.diff(array, axis=0))
-    out[1:] = d  # no need to do max yet
-    # since these are absolute differences |r0-r1| = |r1-r0|
-    # make a view of the target portion of the array
-    v = out[:-1]
-    # compute an in-place maximum
-    np.putmask(v, d > v, d)
-
-    # compute col-wise absolute difference
-    d = np.abs(np.diff(array, axis=1))
-    v = out[:, 1:]
-    np.putmask(v, d > v, d)
-    v = out[:, :-1]
-    np.putmask(v, d > v, d)
-    return out
-
-
 def _abs_deriv(array):
-    # FIXME this assumes off-edge pixels are 0
-    # FIXME this upcasts to float64
-    # FIXME _valid_abs_deriv fixes the above issues and is more efficient
     # but fixing the bugs will likely change the output
     tmp = np.zeros(array.shape, dtype=np.float64)
     out = np.zeros(array.shape, dtype=np.float64)
@@ -123,7 +98,6 @@ def _absolute_subtract(array, tmp, out):
     return tmp, out
 
 
-# TODO add tests
 def flag_crs(
     sci_data,
     sci_err,
@@ -135,7 +109,6 @@ def flag_crs(
     return np.greater(np.abs(sci_data - blot_data), snr * np.nan_to_num(sci_err))
 
 
-# TODO add tests
 def flag_resampled_crs(
     sci_data,
     sci_err,
@@ -212,10 +185,6 @@ def flag_resampled_crs(
     return mask1_smoothed & mask2
 
 
-# FIXME (or fixed) interp and sinscl were "options" only when provided
-# as part of the step spec (which becomes outlierpars). As neither was
-# in the spec (and providing unknown arguments causes an error), these
-# were never configurable and always defaulted to linear and 1.0
 def gwcs_blot(median_data, median_wcs, blot_data, blot_wcs, pix_ratio):
     """
     Resample the output/resampled image to recreate an input image based on
@@ -229,7 +198,6 @@ def gwcs_blot(median_data, median_wcs, blot_data, blot_wcs, pix_ratio):
         Datamodel containing header and WCS to define the 'blotted' image
     """
     # Compute the mapping between the input and output pixel coordinates
-    # TODO stcal.alignment.resample_utils.calc_pixmap does not work here
     pixmap = calc_gwcs_pixmap(blot_wcs, median_wcs, blot_data.shape)
     log.debug("Pixmap shape: {}".format(pixmap[:, :, 0].shape))
     log.debug("Sci shape: {}".format(blot_data.shape))
@@ -248,20 +216,16 @@ def gwcs_blot(median_data, median_wcs, blot_data, blot_wcs, pix_ratio):
     return outsci
 
 
-# TODO tests, duplicate in resample, resample_utils
-def calc_gwcs_pixmap(in_wcs, out_wcs, shape):
+def calc_gwcs_pixmap(in_wcs, out_wcs, in_shape):
     """ Return a pixel grid map from input frame to output frame.
     """
-    bb = wcs_bbox_from_shape(shape)
+    bb = wcs_bbox_from_shape(in_shape)
     log.debug("Bounding box from data shape: {}".format(bb))
 
     grid = gwcs.wcstools.grid_from_bounding_box(bb)
-    pixmap = np.dstack(reproject(in_wcs, out_wcs)(grid[0], grid[1]))
-
-    return pixmap
+    return np.dstack(reproject(in_wcs, out_wcs)(grid[0], grid[1]))
 
 
-# TODO tests, duplicate in resample, assign_wcs, resample_utils
 def reproject(wcs1, wcs2):
     """
     Given two WCSs or transforms return a function which takes pixel
diff --git a/tests/outlier_detection/test_utils.py b/tests/outlier_detection/test_utils.py
index 6e0302119..bcfdaf724 100644
--- a/tests/outlier_detection/test_utils.py
+++ b/tests/outlier_detection/test_utils.py
@@ -1,12 +1,19 @@
 import warnings
 
+import gwcs
 import pytest
 import numpy as np
 import scipy.signal
+from astropy.modeling import models
 
 from stcal.outlier_detection.utils import (
     _abs_deriv,
     compute_weight_threshold,
+    flag_crs,
+    flag_resampled_crs,
+    gwcs_blot,
+    calc_gwcs_pixmap,
+    reproject,
     medfilt,
 )
 
@@ -26,7 +33,7 @@ def test_abs_deriv_single_value(shape, diff):
     np.testing.assert_allclose(result, expected)
 
 
-@pytest.mark.skip(reason="_abs_deriv has edge effects due to treating off-edge pixels as 0")
+@pytest.mark.skip(reason="_abs_deriv has edge effects due to treating off-edge pixels as 0: see JP-3683")
 @pytest.mark.parametrize("nrows,ncols", [(5, 5), (7, 11), (17, 13)])
 def test_abs_deriv_range(nrows, ncols):
     arr = np.arange(nrows * ncols).reshape(nrows, ncols)
@@ -65,6 +72,79 @@ def test_compute_weight_threshold_zeros():
     np.testing.assert_allclose(result, 21)
 
 
+def test_flag_crs():
+    sci = np.zeros((10, 10), dtype=np.float32)
+    err = np.ones_like(sci)
+    blot = np.zeros_like(sci)
+    # add a cr
+    sci[2, 3] = 10
+    crs = flag_crs(sci, err, blot, 1)
+    ys, xs = np.where(crs)
+    np.testing.assert_equal(ys, 2)
+    np.testing.assert_equal(xs, 3)
+
+
+def test_flag_resampled_crs():
+    sci = np.zeros((10, 10), dtype=np.float32)
+    err = np.ones_like(sci)
+    blot = np.zeros_like(sci)
+    # add a cr
+    sci[2, 3] = 10
+
+    snr1, snr2 = 5, 4
+    scale1, scale2 = 1.2, 0.7
+    backg = 0.0
+    resample = True
+    crs = flag_resampled_crs(sci, err, blot, snr1, snr2, scale1, scale2, backg, resample)
+    np.testing.assert_equal(ys, 2)
+    np.testing.assert_equal(xs, 3)
+
+
+def test_gwcs_blot():
+    # set up a very simple wcs that scales by 1x
+    output_frame = gwcs.Frame2D(name="world")
+    forward_transform = models.Scale(1) & models.Scale(1)
+
+    median_data = np.arange(100, dtype=np.float32).reshape((10, 10))
+    median_wcs = gwcs.WCS(forward_transform, output_frame=output_frame)
+    blot_data = np.zeros((5, 5), dtype=np.float32)
+    blot_wcs = gwcs.WCS(forward_transform, output_frame=output_frame)
+    pix_ratio = 1.0
+
+    blotted = gwcs_blot(median_data, median_wcs, blot_data, blot_wcs, pix_ratio)
+    # since the median data is larger and the wcs are equivalent the blot
+    # will window the data to the shape of the blot data
+    assert blotted.shape == blot_data.shape
+    np.testing.assert_equal(blotted, median_data[:blot_data.shape[0], :blot_data.shape[1]])
+
+
+def test_calc_gwcs_pixmap():
+    # generate 2 wcses with different scales
+    output_frame = gwcs.Frame2D(name="world")
+    in_transform = models.Scale(1) & models.Scale(1)
+    out_transform = models.Scale(2) & models.Scale(2)
+    in_wcs = gwcs.WCS(in_transform, output_frame=output_frame)
+    out_wcs = gwcs.WCS(out_transform, output_frame=output_frame)
+    in_shape = (3, 4)
+    pixmap = calc_gwcs_pixmap(in_wcs, out_wcs, in_shape)
+    # we expect given the 2x scale difference to have a pixmap
+    # with pixel coordinates / 2
+    # use mgrid to generate these coordinates (and reshuffle to match the pixmap)
+    expected = np.swapaxes(np.mgrid[:4, :3] / 2., 0, 2)
+    np.testing.assert_equal(pixmap, expected)
+
+
+def test_reproject():
+    # generate 2 wcses with different scales
+    output_frame = gwcs.Frame2D(name="world")
+    wcs1 = gwcs.WCS(models.Scale(1) & models.Scale(1), output_frame=output_frame)
+    wcs2 = gwcs.WCS(models.Scale(2) & models.Scale(2), output_frame=output_frame)
+    project = reproject(wcs1, wcs2)
+    pys, pxs = project(np.array([3]), np.array([1]))
+    np.testing.assert_equal(pys, 1.5)
+    np.testing.assert_equal(pxs, 0.5)
+
+
 @pytest.mark.parametrize("shape,kern_size", [
     ([7, 7], [3, 3]),
     ([7, 7], [3, 1]),