diff --git a/include/eigen b/include/eigen
index 6d829e766..c01ff4531 160000
--- a/include/eigen
+++ b/include/eigen
@@ -1 +1 @@
-Subproject commit 6d829e766ff1b1ab867d93631163cbc63ed5798f
+Subproject commit c01ff45312582b2ea896ee307a49165ca4790332
diff --git a/lib/libcfitsio.a b/lib/libcfitsio.a
deleted file mode 100644
index 6d79f323b..000000000
Binary files a/lib/libcfitsio.a and /dev/null differ
diff --git a/lib/pybind11 b/lib/pybind11
index 3efe9d4cb..741d86f2e 160000
--- a/lib/pybind11
+++ b/lib/pybind11
@@ -1 +1 @@
-Subproject commit 3efe9d4cb5d7314faee722205e560b8b932aed9e
+Subproject commit 741d86f2e3527b667ba85d273a5eea19a0978ef5
diff --git a/pyproject.toml b/pyproject.toml
index 97311fe5b..55077fa27 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -32,15 +32,17 @@ classifiers = [
 ]
 dynamic = ["version"]
 dependencies = [
-    "astropy>=5.1",
-    "astroquery>=0.4.6",
-    "joblib>=1.4",
-    "matplotlib>=3.5",
-    "numpy<2.0",
-    "pandas>=1.5.1",
+    "astropy>=6.1",
+    "astroquery>=0.4.7",
+    "dask[array]>=2024.4.1",  # Needed for astropy library dependency
+    "jax",
+    "jaxlib",
+    "matplotlib>=3.9",
+    "numpy>=2.0",
+    "pandas>=2.2",  # Needed for scikit_learn
     "reproject",
-    "scipy>=1.9.2",
-    "scikit_learn>=1.0.0",
+    "scipy>=1.13",
+    "scikit_learn>=1.5.0",
     "koffi>=0.1.1",
     "jplephem",
     "PyYAML>=6.0"
@@ -54,15 +56,12 @@ Changelog = "https://epyc.astro.washington.edu/~kbmod/project_details/release_no
 
 [project.optional-dependencies]
 analysis = [
-    "tensorflow>=2.9",
-    "matplotlib>=3.6.1",
+    "ephem>=4.1",
     "ipywidgets>=8.0",
-    "ephem>=4.1"
 ]
 docs = [
-    "sphinx",
-    # https://github.com/spatialaudio/nbsphinx/issues/655
-    "sphinx-gallery<0.11",
+    "sphinx>=8.0",
+    "sphinx-gallery",
     "sphinx-codeautolink",
     "sphinx-design",
     "numpydoc",
diff --git a/src/kbmod/file_utils.py b/src/kbmod/file_utils.py
deleted file mode 100644
index fcc1fdde9..000000000
--- a/src/kbmod/file_utils.py
+++ /dev/null
@@ -1,70 +0,0 @@
-"""A collection of utility functions for working with files in KBMOD."""
-
-import csv
-import re
-from collections import OrderedDict
-from itertools import product
-from math import copysign
-from pathlib import Path
-
-import astropy.units as u
-import numpy as np
-from astropy.coordinates import *
-from astropy.io import fits
-from astropy.time import Time
-
-import kbmod.search as kb
-from kbmod.search import LayeredImage
-from kbmod.trajectory_utils import trajectory_from_np_object
-
-
-def load_deccam_layered_image(filename, psf):
-    """Load a layered image from the legacy deccam format.
-
-    Parameters
-    ----------
-    filename : `str`
-        The name of the file to load.
-    psf : `PSF`
-        The PSF to use for the image.
-
-    Returns
-    -------
-    img : `LayeredImage`
-        The loaded image.
-
-    Raises
-    ------
-    Raises a ``FileNotFoundError`` if the file does not exist.
-    Raises a ``ValueError`` if any of the validation checks fail.
-    """
-    if not Path(filename).is_file():
-        raise FileNotFoundError(f"{filename} not found")
-
-    img = None
-    with fits.open(filename) as hdul:
-        if len(hdul) < 4:
-            raise ValueError("Not enough extensions for legacy deccam format")
-
-        # Extract the obstime trying from a few keys and a few extensions.
-        obstime = -1.0
-        for key, ext in product(["MJD", "DATE-AVG", "MJD-OBS"], [0, 1]):
-            if key in hdul[ext].header:
-                value = hdul[ext].header[key]
-                if type(value) is float:
-                    obstime = value
-                    break
-                if type(value) is str:
-                    timesys = hdul[ext].header.get("TIMESYS", "UTC").lower()
-                    obstime = Time(value, scale=timesys).mjd
-                    break
-
-        img = LayeredImage(
-            hdul[1].data.astype(np.float32),  # Science
-            hdul[3].data.astype(np.float32),  # Variance
-            hdul[2].data.astype(np.float32),  # Mask
-            psf,
-            obstime,
-        )
-
-    return img
diff --git a/src/kbmod/filters/sigma_g_filter.py b/src/kbmod/filters/sigma_g_filter.py
index c862cddfe..87cdd67ea 100644
--- a/src/kbmod/filters/sigma_g_filter.py
+++ b/src/kbmod/filters/sigma_g_filter.py
@@ -5,6 +5,9 @@
 by Smotherman et. al. 2021
 """
 
+from functools import partial
+from jax import jit, vmap
+import jax.numpy as jnp
 import logging
 import numpy as np
 from scipy.special import erfinv
@@ -12,9 +15,50 @@
 from kbmod.results import Results
 from kbmod.search import DebugTimer
 
+
 logger = logging.getLogger(__name__)
 
 
+def sigma_g_jax(data, low_bnd, high_bnd, n_sigma, coeff, clip_negative):
+    """The core function for performing a sigma filtering on a series of data points
+    with clipped_negative. These are typically likelihoods for KBMOD.
+
+    Parameters
+    ----------
+    data : `numpy.ndarray`
+        A length T matrix of data points for filtering.
+    low_bnd : `float`
+        The lower bound of the interval to use to estimate the standard deviation.
+    high_bnd : `float`
+        The upper bound of the interval to use to estimate the standard deviation.
+    n_sigma : `float`
+        The number of standard deviations to use for the bound.
+    coeff : `float`
+        The precomputed coefficient based on the given bounds.
+    clip_negative : `bool`
+        A Boolean indicating whether to use negative values when computing
+        standard deviation.
+
+    Returns
+    -------
+    index_valid : `numpy.ndarray`
+        A length T array of Booleans indicating if each point is valid (True)
+        or has been filtered (False).
+    """
+    # Compute the percentiles for this array of values. If we are clipping the negatives then only
+    # use the positive points.
+    masked_data = jnp.where((not clip_negative) | (data > 0.0), data, jnp.nan)
+    lower_per, median, upper_per = jnp.nanpercentile(masked_data, jnp.array([low_bnd, 50, high_bnd]))
+
+    # Compute the bounds for each row, enforcing a minimum gap in case all the
+    # points are identical (upper_per == lower_per).
+    delta = upper_per - lower_per
+    nSigmaG = n_sigma * coeff * jnp.where(delta > 1e-8, delta, 1e-8)
+
+    index_valid = jnp.isfinite(data) & (data <= median + nSigmaG) & (data >= median - nSigmaG)
+    return index_valid
+
+
 class SigmaGClipping:
     """This class contains the basic information for performing SigmaG clipping.
 
@@ -41,9 +85,21 @@ def __init__(self, low_bnd=25, high_bnd=75, n_sigma=2, clip_negative=False):
 
         self.low_bnd = low_bnd
         self.high_bnd = high_bnd
-        self.clip_negative = clip_negative
         self.n_sigma = n_sigma
         self.coeff = SigmaGClipping.find_sigma_g_coeff(low_bnd, high_bnd)
+        self.clip_negative = clip_negative
+
+        # Create compiled vmapped functions that applies the Sigma G filtering
+        # with the given parameters.
+        base_fn = partial(
+            sigma_g_jax,
+            low_bnd=self.low_bnd,
+            high_bnd=self.high_bnd,
+            n_sigma=self.n_sigma,
+            coeff=self.coeff,
+            clip_negative=self.clip_negative,
+        )
+        self.sigma_g_jax_fn = vmap(jit(base_fn))
 
     @staticmethod
     def find_sigma_g_coeff(low_bnd, high_bnd):
@@ -107,15 +163,7 @@ def compute_clipped_sigma_g(self, lh):
         sigmaG = self.coeff * delta
         nSigmaG = self.n_sigma * sigmaG
 
-        # Its unclear why we only filter zeros for one of the two cases, but leaving the logic in
-        # to stay consistent with the original code.
-        if self.clip_negative:
-            good_index = np.where(
-                np.logical_and(lh != 0, np.logical_and(lh > median - nSigmaG, lh < median + nSigmaG))
-            )[0]
-        else:
-            good_index = np.where(np.logical_and(lh > median - nSigmaG, lh < median + nSigmaG))[0]
-
+        good_index = np.where(np.logical_and(lh > median - nSigmaG, lh < median + nSigmaG))[0]
         return good_index
 
     def compute_clipped_sigma_g_matrix(self, lh):
@@ -134,33 +182,8 @@ def compute_clipped_sigma_g_matrix(self, lh):
             A N x T matrix of Booleans indicating if each point is valid (True)
             or has been filtered (False).
         """
-        if self.clip_negative:
-            # We mask out the values less than zero so they are not used in the median computation.
-            masked_lh = np.copy(lh)
-            masked_lh[lh <= 0] = np.nan
-            lower_per, median, upper_per = np.nanpercentile(
-                masked_lh, [self.low_bnd, 50, self.high_bnd], axis=1
-            )
-        else:
-            lower_per, median, upper_per = np.nanpercentile(lh, [self.low_bnd, 50, self.high_bnd], axis=1)
-
-        # Compute the bounds for each row, enforcing a minimum gap in case all the
-        # points are identical (upper_per == lower_per).
-        delta = upper_per - lower_per
-        delta[delta < 1e-8] = 1e-8
-        nSigmaG = self.n_sigma * self.coeff * delta
-
-        num_cols = lh.shape[1]
-        lower_bnd = np.repeat(np.array([median - nSigmaG]).T, num_cols, axis=1)
-        upper_bnd = np.repeat(np.array([median + nSigmaG]).T, num_cols, axis=1)
-
-        # Its unclear why we only filter zeros for one of the two cases, but leaving the logic in
-        # to stay consistent with the original code.
-        if self.clip_negative:
-            index_valid = np.isfinite(lh) & (lh != 0) & (lh < upper_bnd) & (lh > lower_bnd)
-        else:
-            index_valid = np.isfinite(lh) & (lh < upper_bnd) & (lh > lower_bnd)
-        return index_valid
+        inds_valid = self.sigma_g_jax_fn(jnp.asarray(lh)).block_until_ready()
+        return inds_valid
 
 
 def apply_clipped_sigma_g(clipper, result_data):
@@ -183,4 +206,3 @@ def apply_clipped_sigma_g(clipper, result_data):
     obs_valid = clipper.compute_clipped_sigma_g_matrix(lh)
     result_data.update_obs_valid(obs_valid)
     filter_timer.stop()
-    return
diff --git a/src/kbmod/results.py b/src/kbmod/results.py
index 15b6ebac6..f45460d4f 100644
--- a/src/kbmod/results.py
+++ b/src/kbmod/results.py
@@ -331,12 +331,15 @@ def _update_likelihood(self):
         ------
         Raises an IndexError if the necessary columns are missing.
         """
+        num_rows = len(self.table)
+        if num_rows == 0:
+            return  # Nothing to do for an empty table
+
         if "psi_curve" not in self.table.colnames:
             raise IndexError("Missing column 'phi_curve'. Use add_psi_phi_data()")
         if "phi_curve" not in self.table.colnames:
             raise IndexError("Missing column 'phi_curve'. Use add_psi_phi_data()")
 
-        num_rows = len(self.table)
         num_times = len(self.table["phi_curve"][0])
         if "obs_valid" in self.table.colnames:
             phi_sum = (self.table["phi_curve"] * self.table["obs_valid"]).sum(axis=1)
diff --git a/src/kbmod/search/image_stack.cpp b/src/kbmod/search/image_stack.cpp
index ea64ce36e..c8e13b5b0 100644
--- a/src/kbmod/search/image_stack.cpp
+++ b/src/kbmod/search/image_stack.cpp
@@ -144,33 +144,6 @@ void ImageStack::clear_from_gpu() {
     data_on_gpu = false;
 }
 
-RawImage ImageStack::make_global_mask(int flags, int threshold) {
-    uint64_t npixels = get_npixels();
-
-    // Start with an empty global mask.
-    RawImage global_mask = RawImage(width, height);
-    global_mask.set_all(0.0);
-
-    // For each pixel count the number of images where it is masked.
-    std::vector<int> counts(npixels, 0);
-    for (unsigned int img = 0; img < images.size(); ++img) {
-        auto imgMask = images[img].get_mask().get_image().reshaped();
-
-        // Count the number of times a pixel has any of the given flags
-        for (uint64_t pixel = 0; pixel < npixels; ++pixel) {
-            if ((flags & static_cast<int>(imgMask[pixel])) != 0) counts[pixel]++;
-        }
-    }
-
-    // Set all pixels below threshold to 0 and all above to 1
-    auto global_m = global_mask.get_image().reshaped();
-    for (uint64_t p = 0; p < npixels; ++p) {
-        global_m[p] = counts[p] < threshold ? 0.0 : 1.0;
-    }
-
-    return global_mask;
-}
-
 #ifdef Py_PYTHON_H
 static void image_stack_bindings(py::module& m) {
     using is = search::ImageStack;
@@ -194,7 +167,6 @@ static void image_stack_bindings(py::module& m) {
             .def("build_zeroed_times", &is::build_zeroed_times, pydocs::DOC_ImageStack_build_zeroed_times)
             .def("sort_by_time", &is::sort_by_time, pydocs::DOC_ImageStack_sort_by_time)
             .def("img_count", &is::img_count, pydocs::DOC_ImageStack_img_count)
-            .def("make_global_mask", &is::make_global_mask, pydocs::DOC_ImageStack_make_global_mask)
             .def("convolve_psf", &is::convolve_psf, pydocs::DOC_ImageStack_convolve_psf)
             .def("get_width", &is::get_width, pydocs::DOC_ImageStack_get_width)
             .def("get_height", &is::get_height, pydocs::DOC_ImageStack_get_height)
diff --git a/src/kbmod/search/image_stack.h b/src/kbmod/search/image_stack.h
index 20c59abc4..0eb5310d2 100644
--- a/src/kbmod/search/image_stack.h
+++ b/src/kbmod/search/image_stack.h
@@ -47,9 +47,6 @@ class ImageStack {
 
     void convolve_psf();
 
-    // Make and return a global mask.
-    RawImage make_global_mask(int flags, int threshold);
-
     virtual ~ImageStack();
 
     // Functions to handle transfering data to/from GPU.
diff --git a/src/kbmod/search/layered_image.cpp b/src/kbmod/search/layered_image.cpp
index e6f7eab51..19affd537 100644
--- a/src/kbmod/search/layered_image.cpp
+++ b/src/kbmod/search/layered_image.cpp
@@ -119,22 +119,6 @@ void LayeredImage::apply_mask(int flags) {
     variance.apply_mask(flags, mask);
 }
 
-void LayeredImage::subtract_template(RawImage& sub_template) {
-    assert_sizes_equal(sub_template.get_width(), width, "template width");
-    assert_sizes_equal(sub_template.get_height(), height, "template height");
-    const uint64_t num_pixels = get_npixels();
-
-    logging::getLogger("kbmod.search.layered_image")->debug("Subtracting template image.");
-
-    float* sci_pixels = science.data();
-    float* tem_pixels = sub_template.data();
-    for (uint64_t i = 0; i < num_pixels; ++i) {
-        if (pixel_value_valid(sci_pixels[i]) && pixel_value_valid(tem_pixels[i])) {
-            sci_pixels[i] -= tem_pixels[i];
-        }
-    }
-}
-
 void LayeredImage::set_science(RawImage& im) {
     assert_sizes_equal(im.get_width(), width, "science layer width");
     assert_sizes_equal(im.get_height(), height, "science layer height");
@@ -212,19 +196,6 @@ RawImage LayeredImage::generate_phi_image() {
     return result;
 }
 
-double LayeredImage::compute_fraction_masked() const {
-    double masked_count = 0.0;
-    double total_count = 0.0;
-
-    for (int j = 0; j < height; ++j) {
-        for (int i = 0; i < width; ++i) {
-            if (!science_pixel_has_data({j, i})) masked_count += 1.0;
-            total_count++;
-        }
-    }
-    return masked_count / total_count;
-}
-
 #ifdef Py_PYTHON_H
 static void layered_image_bindings(py::module& m) {
     using li = search::LayeredImage;
@@ -268,7 +239,6 @@ static void layered_image_bindings(py::module& m) {
                  })
             .def("binarize_mask", &li::binarize_mask, pydocs::DOC_LayeredImage_binarize_mask)
             .def("apply_mask", &li::apply_mask, pydocs::DOC_LayeredImage_apply_mask)
-            .def("sub_template", &li::subtract_template, pydocs::DOC_LayeredImage_sub_template)
             .def("get_science", &li::get_science, py::return_value_policy::reference_internal,
                  pydocs::DOC_LayeredImage_get_science)
             .def("get_mask", &li::get_mask, py::return_value_policy::reference_internal,
@@ -285,8 +255,6 @@ static void layered_image_bindings(py::module& m) {
             .def("get_npixels", &li::get_npixels, pydocs::DOC_LayeredImage_get_npixels)
             .def("get_obstime", &li::get_obstime, pydocs::DOC_LayeredImage_get_obstime)
             .def("set_obstime", &li::set_obstime, pydocs::DOC_LayeredImage_set_obstime)
-            .def("compute_fraction_masked", &li::compute_fraction_masked,
-                 pydocs::DOC_LayeredImage_compute_fraction_masked)
             .def("generate_psi_image", &li::generate_psi_image, pydocs::DOC_LayeredImage_generate_psi_image)
             .def("generate_phi_image", &li::generate_phi_image, pydocs::DOC_LayeredImage_generate_phi_image);
 }
diff --git a/src/kbmod/search/layered_image.h b/src/kbmod/search/layered_image.h
index 5b5edf156..570448541 100644
--- a/src/kbmod/search/layered_image.h
+++ b/src/kbmod/search/layered_image.h
@@ -68,9 +68,6 @@ class LayeredImage {
     void binarize_mask(int flags_to_keep);
     void apply_mask(int flags);
 
-    // Subtracts a template image from the science layer.
-    void subtract_template(RawImage& sub_template);
-
     // Setter functions for the individual layers.
     void set_science(RawImage& im);
     void set_mask(RawImage& im);
@@ -86,9 +83,6 @@ class LayeredImage {
     RawImage generate_psi_image();
     RawImage generate_phi_image();
 
-    // Debugging and statistics functions.
-    double compute_fraction_masked() const;
-
 private:
     unsigned width;
     unsigned height;
diff --git a/src/kbmod/search/pydocs/image_stack_docs.h b/src/kbmod/search/pydocs/image_stack_docs.h
index b6f8344e5..a500175c3 100644
--- a/src/kbmod/search/pydocs/image_stack_docs.h
+++ b/src/kbmod/search/pydocs/image_stack_docs.h
@@ -130,27 +130,6 @@ static const auto DOC_ImageStack_sort_by_time = R"doc(
   Raises a ``RuntimeError`` if the input image is the data is currently on GPU.    
   )doc";
 
-static const auto DOC_ImageStack_make_global_mask = R"doc(
-  Create a new global mask from a set of flags and a threshold.
-  The global mask marks a pixel as masked if and only if it is masked
-  by one of the given flags in at least ``threshold`` individual images.
-  The returned mask is binary.
-
-  Parameters
-  ----------
-  flags : `int`
-      A bit mask of mask flags to use when counting.
-  threshold : `int`
-      The minimum number of images in which a pixel must be masked to be
-      part of the global mask.
-
-  Returns
-  -------
-  global_mask : `RawImage`
-      A RawImage containing the global mask with 1 for masked pixels
-      and 0 for unmasked pixels.
-  )doc";
-
 static const auto DOC_ImageStack_convolve_psf = R"doc(
   Convolves each image (science and variance layers) with the PSF
   stored in the LayeredImage object.
diff --git a/src/kbmod/search/pydocs/layered_image_docs.h b/src/kbmod/search/pydocs/layered_image_docs.h
index 02ce08ba0..0936157be 100644
--- a/src/kbmod/search/pydocs/layered_image_docs.h
+++ b/src/kbmod/search/pydocs/layered_image_docs.h
@@ -66,10 +66,6 @@ static const auto DOC_LayeredImage_apply_mask = R"doc(
   call binarize_mask() first.
   )doc";
 
-static const auto DOC_LayeredImage_sub_template = R"doc(
-  Subtract given image template
-  )doc";
-
 static const auto DOC_LayeredImage_get_science = R"doc(
   Returns the science layer raw_image.
   )doc";
@@ -222,16 +218,6 @@ static const auto DOC_LayeredImage_generate_phi_image = R"doc(
       A ``RawImage`` of the same dimensions as the ``LayeredImage``.
   )doc";
 
-static const auto DOC_LayeredImage_compute_fraction_masked = R"doc(
-  Computes the fraction of pixels in the layered image that are masked.
-  This can be used for debugging purposes.
-
-  Returns
-  -------
-  value : `double`
-      The fraction of pixels in the science image that are masked.
-  )doc";
-
 }  // namespace pydocs
 
 #endif /* LAYEREDIMAGE_DOCS  */
diff --git a/src/kbmod/trajectory_utils.py b/src/kbmod/trajectory_utils.py
index 9163964e7..423febad7 100644
--- a/src/kbmod/trajectory_utils.py
+++ b/src/kbmod/trajectory_utils.py
@@ -19,6 +19,55 @@
 from kbmod.search import Trajectory
 
 
+def predict_pixel_locations(times, x0, vx, centered=True, as_int=True):
+    """A vectorized Python implementation of the logic to predict the
+    pixel locations from a starting pixel and a velocity.
+
+    Parameters
+    ----------
+    times : list-like
+        The length T list of zero-shifted times.
+    x0 : list-like
+        The length R list of starting pixel locations.
+    vx : list-like
+        The length R list of pixel velocities (in pixels per day) for each
+        trajectory.
+    centered : `bool`
+        Shift the prediction to be at the center of the pixel
+        (e.g. xp = x0 + vx * time + 0.5f).
+        Default = True.
+    as_int : `bool`
+        Return the predictions as integers.
+        Default = True.
+
+    Returns
+    -------
+    pos : `numpy.ndarray`
+        A R x T matrix where R is the number of trajectories (length of x0 and vx)
+        and T is the number of times.
+    """
+    # Make sure everything is a numpy array.
+    times = np.asarray(times)
+    x0 = np.asarray(x0)
+    vx = np.asarray(vx)
+
+    # Check the arrays are the same length.
+    if len(x0) != len(vx):
+        raise ValueError(f"x0 and vx must be same size. Found {len(x0)} vs {len(vx)}")
+
+    # Compute the (floating point) predicted pixel position.
+    pos = vx[:, np.newaxis] * times[np.newaxis, :] + x0[:, np.newaxis]
+    if centered:
+        pos = pos + 0.5
+
+    # If returned as int, we do not use an explicit floor in order to stay
+    # consistent with the existing implementation.
+    if as_int:
+        pos = pos.astype(int)
+
+    return pos
+
+
 def make_trajectory_from_ra_dec(ra, dec, v_ra, v_dec, wcs):
     """Create a trajectory object from (RA, dec) information.
 
diff --git a/tests/test_image_stack.py b/tests/test_image_stack.py
index a3eb346da..6ccade166 100644
--- a/tests/test_image_stack.py
+++ b/tests/test_image_stack.py
@@ -98,23 +98,6 @@ def test_times(self):
         for i in range(self.num_images):
             self.assertEqual(times[i], 2.0 * i)
 
-    def test_make_global_mask(self):
-        # Mask a single point in each image with flag=2.
-        for i in range(self.num_images):
-            self.images[i].get_mask().set_pixel(5, 5, 2)
-
-        # Apply the global mask for flag=1 and a threshold of the bit set
-        # in at least one mask. Note this ignores flag=2 for the count.
-        mask = self.im_stack.make_global_mask(1, 1)
-
-        # Check that the correct pixels are masked
-        for y in range(self.im_stack.get_height()):
-            for x in range(self.im_stack.get_width()):
-                if y == 10 and x >= 10 and x <= 10 + (self.num_images - 1):
-                    self.assertEqual(mask.get_pixel(y, x), 1)
-                else:
-                    self.assertEqual(mask.get_pixel(y, x), 0)
-
     # WOW, this is the first test that caught the fact that interpolated_add
     # called add, and that add had flipped i and j by accident. The first one.
     # TODO: more clean understandable tests for basic functionality, these big
diff --git a/tests/test_layered_image.py b/tests/test_layered_image.py
index e2b44bcb6..e17a43f18 100644
--- a/tests/test_layered_image.py
+++ b/tests/test_layered_image.py
@@ -181,20 +181,6 @@ def test_mask_pixel(self):
                 self.assertEqual(pixel_value_valid(pix_val), expected)
                 self.assertEqual(self.image.science_pixel_has_data(y, x), expected)
 
-    def test_compute_fraction_masked(self):
-        total_pixels = self.width * self.height
-
-        # Mask 50 pixels
-        for y in range(0, 10):
-            for x in range(0, 5):
-                self.image.mask_pixel(y, x)
-        self.assertAlmostEqual(self.image.compute_fraction_masked(), 50.0 / total_pixels)
-
-        # Mask another 25 pixels.
-        for x in range(3, 28):
-            self.image.mask_pixel(12, x)
-        self.assertAlmostEqual(self.image.compute_fraction_masked(), 75.0 / total_pixels)
-
     def test_binarize_mask(self):
         # Mask out a range of pixels.
         mask = self.image.get_mask()
@@ -284,39 +270,6 @@ def test_psi_and_phi_image(self):
                 else:
                     self.assertFalse(pixel_value_valid(phi.get_pixel(y, x)))
 
-    def test_subtract_template(self):
-        sci = self.image.get_science()
-        sci.mask_pixel(7, 10)
-        sci.mask_pixel(7, 11)
-        sci.set_pixel(7, 12, math.nan)
-        sci.set_pixel(7, 13, np.nan)
-        old_sci = RawImage(sci.image.copy())  # Make a copy.
-
-        template = RawImage(self.image.get_width(), self.image.get_height())
-        template.set_all(0.0)
-        for h in range(sci.height):
-            for w in range(4, sci.width):
-                template.set_pixel(h, w, 0.01 * h)
-        self.image.sub_template(template)
-
-        for y in range(sci.height):
-            for x in range(sci.width):
-                if y == 7 and (x >= 10 and x <= 13):
-                    self.assertFalse(sci.pixel_has_data(y, x))
-                elif x < 4:
-                    val1 = old_sci.get_pixel(y, x)
-                    val2 = sci.get_pixel(y, x)
-                    self.assertEqual(val1, val2)
-                else:
-                    val1 = old_sci.get_pixel(y, x) - 0.01 * y
-                    val2 = sci.get_pixel(y, x)
-                    self.assertAlmostEqual(val1, val2, delta=1e-5)
-
-        # Test that we fail when the template size does not match.
-        template2 = RawImage(self.image.get_width(), self.image.get_height() + 1)
-        template2.set_all(0.0)
-        self.assertRaises(RuntimeError, self.image.sub_template, template2)
-
 
 if __name__ == "__main__":
     unittest.main()
diff --git a/tests/test_results.py b/tests/test_results.py
index 56bc73d73..21a5babfd 100644
--- a/tests/test_results.py
+++ b/tests/test_results.py
@@ -53,6 +53,9 @@ def test_empty(self):
         self.assertEqual(len(table.colnames), 7)
         self.assertEqual(table.get_num_times(), 0)
 
+        # Check that we don't crash on updating the likelihoods.
+        table._update_likelihood()
+
     def test_from_trajectories(self):
         table = Results.from_trajectories(self.trj_list)
         self.assertEqual(len(table), self.num_entries)
diff --git a/tests/test_sigma_g_filter.py b/tests/test_sigma_g_filter.py
index edddcffe0..031e4e147 100644
--- a/tests/test_sigma_g_filter.py
+++ b/tests/test_sigma_g_filter.py
@@ -84,7 +84,7 @@ def test_sigma_g_negative_clipping(self):
         params = SigmaGClipping(clip_negative=True)
         result = params.compute_clipped_sigma_g(lh)
         for i in range(num_points):
-            self.assertEqual(i in result, i > 2 and i != 5 and i != 14)
+            self.assertEqual(i in result, i > 2 and i != 14)
 
     def test_sigma_g_all_negative_clipping(self):
         num_points = 10
@@ -106,10 +106,11 @@ def test_sigma_g_clipping_matrix_negative_clipping(self):
         expected = np.array(
             [
                 [True] * num_points,
-                [False, False, False, True, True, False] + [True] * (num_points - 6),
+                [False, False, False] + [True] * (num_points - 3),
                 [False] * num_points,
             ]
         )
+
         sigma_g = SigmaGClipping(clip_negative=True)
 
         # Surpress the warning we get from encountering a row of all NaNs.
@@ -144,6 +145,37 @@ def test_apply_clipped_sigma_g_results(self):
             for j in range(i, num_times):
                 self.assertTrue(valid[j])
 
+    def test_sigmag_parity(self):
+        """Test that we get the same results when using the batch and the non-batch methods."""
+        num_tests = 20
+
+        # Run the test with differing numbers of points and with/without clipping.
+        for num_obs in [10, 20, 50]:
+            for clipped in [True, False]:
+                for num_extreme in [0, 1, 2, 3]:
+                    with self.subTest(
+                        num_obs_used=num_obs, use_clipped=clipped, num_extreme_used=num_extreme
+                    ):
+                        # Generate the data from a fixed random seed (same for every subtest).
+                        rng = np.random.default_rng(100)
+                        data = 10.0 * rng.random((num_tests, num_obs)) - 0.5
+
+                        # Add extreme values for each row.
+                        for row in range(num_tests):
+                            for ext_num in range(num_extreme):
+                                idx = int(num_obs * rng.random())
+                                data[row, idx] = 100.0 * rng.random() - 50.0
+
+                            clipper = SigmaGClipping(25, 75, clip_negative=clipped)
+
+                        batch_res = clipper.compute_clipped_sigma_g_matrix(data)
+                        for row in range(num_tests):
+                            # Compute the individual results (as indices) and convert
+                            # those into a vector of bools for comparison.
+                            ind_res = clipper.compute_clipped_sigma_g(data[row])
+                            ind_bools = [(idx in ind_res) for idx in range(num_obs)]
+                            self.assertTrue(np.array_equal(batch_res[row], ind_bools))
+
     def test_sigmag_computation(self):
         self.assertAlmostEqual(SigmaGClipping.find_sigma_g_coeff(25.0, 75.0), 0.7413, delta=0.001)
         self.assertRaises(ValueError, SigmaGClipping.find_sigma_g_coeff, -1.0, 75.0)
diff --git a/tests/test_trajectory_utils.py b/tests/test_trajectory_utils.py
index 53a5c84d0..896e35806 100644
--- a/tests/test_trajectory_utils.py
+++ b/tests/test_trajectory_utils.py
@@ -8,6 +8,30 @@
 
 
 class test_trajectory_utils(unittest.TestCase):
+    def test_predict_pixel_locations(self):
+        x0 = [0.0, 1.0, 2.0]
+        vx = [1.0, 0.0, -1.0]
+        times = [0, 1, 2, 3]
+
+        # When centering and using an integer, the last point for x0=2.0 vx=-1.0
+        # is at -0.5 which is rounded to 0.0 because we are using int truncation
+        # instead of an explicit floor for consistency.
+        pos = predict_pixel_locations(times, x0, vx, centered=True, as_int=True)
+        expected = np.array([[0, 1, 2, 3], [1, 1, 1, 1], [2, 1, 0, 0]])
+        self.assertTrue(np.array_equal(pos, expected))
+
+        pos = predict_pixel_locations(times, x0, vx, centered=False, as_int=True)
+        expected = np.array([[0, 1, 2, 3], [1, 1, 1, 1], [2, 1, 0, -1]])
+        self.assertTrue(np.array_equal(pos, expected))
+
+        pos = predict_pixel_locations(times, x0, vx, centered=True, as_int=False)
+        expected = np.array([[0.5, 1.5, 2.5, 3.5], [1.5, 1.5, 1.5, 1.5], [2.5, 1.5, 0.5, -0.5]])
+        self.assertTrue(np.array_equal(pos, expected))
+
+        pos = predict_pixel_locations(times, x0, vx, centered=False, as_int=False)
+        expected = np.array([[0.0, 1.0, 2.0, 3.0], [1.0, 1.0, 1.0, 1.0], [2.0, 1.0, 0.0, -1.0]])
+        self.assertTrue(np.array_equal(pos, expected))
+
     def test_predict_skypos(self):
         # Create a fake WCS with a known pointing.
         my_wcs = WCS(naxis=2)