From 7332543774451725062e2b3c7738972ca0f1add8 Mon Sep 17 00:00:00 2001
From: DinoBektesevic <dinob@uw.edu>
Date: Thu, 6 Jul 2023 17:46:35 -0700
Subject: [PATCH] Implement toLayeredImage, toImageStack and run method.
 Cleanup.

Implement the methods required to run KBMOD on ImageCollection.
Cleanup ImageCollection behaviour:
*) return image collection when indexed by lists, arrays and slices
*) return Row when indexed by integer
*) return Table when sliced by columns
*) Rename the exts to processable. Alias it to a property so that
   each Standardizer can implement its own internal structure the
   way it wants (but also because I was too lazy to rename everything)
*) Fix documentation
*) Move WCS and BBOX as properties to a Standardizer - if that's
   where we need to explain why they are special that's where they
   need to live. Make them an abstractproperty and demand that the
   Standardizers return a list of None's if need be.
*) Fix forceStandardizer keyword (again).
*) Add toLayeredImage as an abstract method to the Standardizers
   Implement them for the three example Standardizers we have.
*) Add toImageStack as an abstract method to the standardizers.
   Implment them in ImageCollection
*) Add run method prototype to ImageCollection to showcase how
   we can neatly integrate with the ImageCollection to execute
   KBMOD runs.

Write an example python script showcasing most of this functionality.

TODO: tests, unittests, integrationtests all the tests.
---
 src/kbmod/__init__.py                         |   6 +-
 src/kbmod/image_collection.py                 | 247 ++++++++++++++++--
 src/kbmod/standardizer.py                     | 179 +++++++++----
 .../standardizers/butler_standardizer.py      |  75 +++++-
 .../decam_community_pipe_fits.py              |   9 +-
 src/kbmod/standardizers/fits_standardizer.py  |  44 +++-
 .../standardizers/multi_extension_fits.py     |   2 +-
 src/kbmod/standardizers/rubin_scipipe_fits.py |  13 +-
 8 files changed, 479 insertions(+), 96 deletions(-)

diff --git a/src/kbmod/__init__.py b/src/kbmod/__init__.py
index c816ba79e..b3bc7c267 100644
--- a/src/kbmod/__init__.py
+++ b/src/kbmod/__init__.py
@@ -5,7 +5,6 @@
 except ImportError:
     warnings.warn("Unable to determine the package version. " "This is likely a broken installation.")
 
-from .standardizers import *
 from . import (
     analysis,
     analysis_utils,
@@ -15,7 +14,8 @@
     jointfit_functions,
     result_list,
     run_search,
-    standardizer,
 )
 
-
+from .standardizers import *
+from .standardizer import Standardizer
+from .image_collection import ImageCollection
diff --git a/src/kbmod/image_collection.py b/src/kbmod/image_collection.py
index 140fcc1f5..94e90189d 100644
--- a/src/kbmod/image_collection.py
+++ b/src/kbmod/image_collection.py
@@ -6,18 +6,24 @@
 import os
 import glob
 import json
+import time
 
+import astropy.units as u
 from astropy.io import fits
 from astropy.time import Time
 from astropy.table import Table
 from astropy.wcs import WCS
 from astropy.utils import isiterable
+from astropy.coordinates import SkyCoord
 
 import numpy as np
 
-from kbmod.file_utils import FileUtils
-from kbmod.search import pixel_pos, layered_image
+from kbmod.search import image_stack, stack_search
 from kbmod.standardizer import Standardizer
+from kbmod.analysis_utils import PostProcess
+
+
+__all__  = ["ImageCollection", ]
 
 
 class ImageCollection:
@@ -134,14 +140,15 @@ def __init__(self, metadata, standardizers=None):
 
         if standardizers is not None:
             self.data.meta["n_entries"] = len(standardizers)
-            self._standardizers = standardizers
+            self._standardizers = np.array(standardizers)
         elif metadata.meta and "n_entries" in metadata.meta:
             n_entries = metadata.meta["n_entries"]
-            self._standardizers = [None]*n_entries
+            self._standardizers = np.full((n_entries, ), None)
         else:
             n_entries = len(np.unique(metadata["location"]))
             self.data.meta["n_entries"] = n_entries
-            self._standardizers = [None]*n_entries
+            self._standardizers = np.full((n_entries, ), None)
+            #self._standardizers = [None]*n_entries
 
         # hidden indices that track the unravelled lookup to standardizer
         # extension index. I should imagine there's a better than double-loop
@@ -157,7 +164,9 @@ def __init__(self, metadata, standardizers=None):
             no_std_map = True
             self.data["ext_idx"] = [None]*len(self.data)
 
-        if standardizers and no_std_map:
+        # standardizers are not Falsy - empty lists, Nones, empty tuples, empty
+        # arrays etc...
+        if (standardizers is not None  and any(standardizers)) and no_std_map:
             std_idxs, ext_idxs = [], []
             for i, stdFits in enumerate(standardizers):
                 for j, ext in enumerate(stdFits.exts):
@@ -195,7 +204,7 @@ def read(cls, *args, format=None, units=None, descriptions=None, **kwargs):
         return cls(metadata)
 
     @classmethod
-    def _fromStandardizers(cls, standardizers, meta=None):
+    def fromStandardizers(cls, standardizers, meta=None):
         """Create ImageCollection from a collection `Standardizers`.
 
         The `Standardizer` is "unravelled", i.e. the shared metadata is
@@ -224,7 +233,7 @@ def _fromStandardizers(cls, standardizers, meta=None):
             # ButlerStd.stdMeta. Everything that is an iterable, except for a
             # string because that could be a location key?
             unravelColumns = [key for key, val in stdMeta.items() if isiterable(val) and not isinstance(val, str)]
-            for j, ext in enumerate(stdFits.exts):
+            for j, ext in enumerate(stdFits.processable):
                 row = {}
                 for key in stdMeta.keys():
                     if key in unravelColumns:
@@ -263,7 +272,7 @@ def _fromFilepaths(cls, filepaths, forceStandardizer, **kwargs):
             Standardizer.fromFile(path=path, forceStandardizer=forceStandardizer, **kwargs)
             for path in filepaths
         ]
-        return cls._fromStandardizers(standardizers)
+        return cls.fromStandardizers(standardizers)
 
     @classmethod
     def _fromDir(cls, path, recursive, forceStandardizer, **kwargs):
@@ -357,7 +366,7 @@ def fromDatasetRefs(cls, butler, refs, **kwargs):
         standardizer_cls = Standardizer.get(standardizer="ButlerStandardizer")
         standardizer = standardizer_cls(butler, refs, **kwargs)
         meta = {"root": butler.datastore.root.geturl(), "n_entries": len(list(refs))}
-        return cls._fromStandardizers([standardizer, ], meta=meta)
+        return cls.fromStandardizers([standardizer, ], meta=meta)
 
     def fromAQueryTable(self):  # ? TBD
         pass
@@ -372,7 +381,12 @@ def __repr__(self):
         return repr(self.data).replace("Table", "ImageInfoSet")
 
     def __getitem__(self, key):
-        return self.data[key]
+        if isinstance(key, (int, str, np.integer)):
+            return self.data[self._userColumns][key]
+        elif isinstance(key, (list, np.ndarray, slice)):
+            return self.__class__(self.data[key], standardizers=self._standardizers[key])
+        else:
+            return self.data[key]
 
     def __setitem__(self, key, val):
         self.data[key] = val
@@ -417,7 +431,7 @@ def standardizers(self):
         """A list of used standardizer names."""
         return self._standardizer_names
 
-    def _get_standardizer(self, index):
+    def get_standardizer(self, index, **kwargs):
         """Get the standardizer and extension index for the selected row of the
         unravelled metadata table.
 
@@ -429,6 +443,8 @@ def _get_standardizer(self, index):
         index : `int`
             Index, as it appears in the unravelled table of metadata
             properties.
+        **kwargs : `dict`
+            Keyword arguments are passed onto the Standardizer constructor.
 
         Returns
         -------
@@ -441,13 +457,13 @@ def _get_standardizer(self, index):
         ext_idx = row["ext_idx"]
         if self._standardizers[std_idx] is None:
             std_cls = Standardizer.registry[row["std_name"]]
-            self._standardizers[std_idx] = std_cls(row["location"])
+            self._standardizers[std_idx] = std_cls(**kwargs, **row)
 
         # maybe a clever dataclass to shortcut the idx lookups on the user end?
         return {"std": self._standardizers[std_idx],
                 "ext": self.data[index]["ext_idx"]}
 
-    def get_standardizers(self, idxs):
+    def get_standardizers(self, idxs, **kwargs):
         """ Get the standardizers used to extract metadata of the selected
         rows.
 
@@ -455,6 +471,9 @@ def get_standardizers(self, idxs):
         ----------
         idx : `int` or `iterable`
             Index of the row for which to retrieve the Standardizer.
+        **kwargs : `dict`
+            Keyword arguments are passed onto the constructors of the retrieved
+            Standardizer.
 
         Returns
         -------
@@ -463,19 +482,13 @@ def get_standardizers(self, idxs):
             the extension (``ext``) that maps to the given metadata row index.
         """
         if isinstance(idxs, int):
-            return self._get_standardizer(idxs)
+            return [self.get_standardizer(idxs, **kwargs), ]
         else:
-            return [self._get_standardizer(idx) for idx in idxs]
+            return [self.get_standardizer(idx, **kwargs) for idx in idxs]
 
     ########################
     # FUNCTIONALITY (object operations, transformative functionality)
     ########################
-    def toImageStack(self):
-        pass
-
-    def plot_onsky(self):
-        pass
-
     def write(self, *args, format=None, serialize_method=None, **kwargs):
         tmpdata = self.data.copy()
 
@@ -522,3 +535,193 @@ def get_duration(self):
         """
         # maybe timespan?
         return self.data["mjd"][-1] - self.data["mjd"][0]
+
+    def toImageStack(self):
+        """Return an `~kbmod.search.image_stack` object for processing with
+        KBMOD.
+
+        Returns
+        -------
+        imageStack : `~kbmod.search.image_stack`
+            Image stack for processing with KBMOD.
+        """
+        # unpack the layred image list to flatten the array
+        # this is so stupidly costly because we have an internal array
+        # representation that doesn't interface with numpy via ndarray it makes
+        # a copy every time
+        layeredImages = [img for std in self._standardizers for img in std.toLayeredImage()]
+        return image_stack(layeredImages)
+
+    def _calc_suggested_angle(self, wcs, center_pixel=(1000, 2000), step=12):
+        """Projects an unit-vector parallel with the ecliptic onto the image
+        and calculates the angle of the projected unit-vector in the pixel
+        space.
+
+        Parameters
+        ----------
+        wcs : ``astropy.wcs.WCS``
+            World Coordinate System object.
+        center_pixel : tuple, array-like
+            Pixel coordinates of image center.
+        step : ``float`` or ``int``
+            Size of step, in arcseconds, used to find the pixel coordinates of
+                the second pixel in the image parallel to the ecliptic.
+
+        Returns
+        -------
+        suggested_angle : ``float``
+            Angle the projected unit-vector parallel to the ecliptic
+            closes with the image axes. Used to transform the specified
+            search angles, with respect to the ecliptic, to search angles
+            within the image.
+
+        Note
+        ----
+        It is not neccessary to calculate this angle for each image in an
+        image set if they have all been warped to a common WCS.
+
+        See Also
+        --------
+        run_search.do_gpu_search
+        """
+        # pick a starting pixel approximately near the center of the image
+        # convert it to ecliptic coordinates
+        start_pixel = np.array(center_pixel)
+        start_pixel_coord = SkyCoord.from_pixel(start_pixel[0], start_pixel[1], wcs)
+        start_ecliptic_coord = start_pixel_coord.geocentrictrueecliptic
+
+        # pick a guess pixel by moving parallel to the ecliptic
+        # convert it to pixel coordinates for the given WCS
+        guess_ecliptic_coord = SkyCoord(
+            start_ecliptic_coord.lon + step * u.arcsec,
+            start_ecliptic_coord.lat,
+            frame="geocentrictrueecliptic",
+        )
+        guess_pixel_coord = guess_ecliptic_coord.to_pixel(wcs)
+
+        # calculate the distance, in pixel coordinates, between the guess and
+        # the start pixel. Calculate the angle that represents in the image.
+        x_dist, y_dist = np.array(guess_pixel_coord) - start_pixel
+        return np.arctan2(y_dist, x_dist)
+
+    def run(self, config):
+        """Run KBMOD on the images in collection.
+
+        Parameters
+        ----------
+        config : `~kbmod.configuration.KBMODConfig`
+            Processing configuration
+
+        Returns
+        -------
+        results : `kbmod.results.ResultList`
+            KBMOD search results.
+
+        Notes
+        -----
+        Requires WCS.
+        """
+        imageStack = self.toImageStack()
+
+        # Compute the ecliptic angle for the images. Assume they are all the
+        # same size? Technically that is currently a requirement, although it's
+        # not explicit (can this be in C++ code?)
+        center_pixel = (imageStack.get_width()/2, imageStack.get_height()/2)
+        suggested_angle = self._calc_suggested_angle(self.wcs[0], center_pixel)
+
+        # Set up the post processing data structure.
+        kb_post_process = PostProcess(config, self.data["mjd"].data)
+
+        # Perform the actual search.
+        search = stack_search(imageStack)
+        # search, search_params = self.do_gpu_search(search, img_info,
+        #                                            suggested_angle, kb_post_process)
+        # not sure why these were separated, I guess it made it look neater?
+        # definitely doesn't feel like everything is in place if there are so
+        # many ifs for a config - feels like that should be a config job?
+        # Anyhow, I'll be lazy and just unravel this here.
+        search_params = {}
+
+        # Run the grid search
+        # Set min and max values for angle and velocity
+        if config["average_angle"] == None:
+            average_angle = suggested_angle
+        else:
+            average_angle = config["average_angle"]
+        ang_min = average_angle - config["ang_arr"][0]
+        ang_max = average_angle + config["ang_arr"][1]
+        vel_min = config["v_arr"][0]
+        vel_max = config["v_arr"][1]
+        search_params["ang_lims"] = [ang_min, ang_max]
+        search_params["vel_lims"] = [vel_min, vel_max]
+
+        # Set the search bounds.
+        if config["x_pixel_bounds"] and len(config["x_pixel_bounds"]) == 2:
+            search.set_start_bounds_x(config["x_pixel_bounds"][0], config["x_pixel_bounds"][1])
+        elif config["x_pixel_buffer"] and config["x_pixel_buffer"] > 0:
+            width = search.get_image_stack().get_width()
+            search.set_start_bounds_x(-config["x_pixel_buffer"], width + config["x_pixel_buffer"])
+
+        if config["y_pixel_bounds"] and len(config["y_pixel_bounds"]) == 2:
+            search.set_start_bounds_y(config["y_pixel_bounds"][0], config["y_pixel_bounds"][1])
+        elif config["y_pixel_buffer"] and config["y_pixel_buffer"] > 0:
+            height = search.get_image_stack().get_height()
+            search.set_start_bounds_y(-config["y_pixel_buffer"], height + config["y_pixel_buffer"])
+
+        # If we are using barycentric corrections, compute the parameters and
+        # enable it in the search function. This can't be not-none atm because
+        # I hadn't copied bary_corr over....
+        if config["bary_dist"] is not None:
+            bary_corr = self._calc_barycentric_corr(img_info, config["bary_dist"])
+            # print average barycentric velocity for debugging
+
+            mjd_range = img_info.get_duration()
+            bary_vx = bary_corr[-1, 0] / mjd_range
+            bary_vy = bary_corr[-1, 3] / mjd_range
+            bary_v = np.sqrt(bary_vx * bary_vx + bary_vy * bary_vy)
+            bary_ang = np.arctan2(bary_vy, bary_vx)
+            print("Average Velocity from Barycentric Correction", bary_v, "pix/day", bary_ang, "angle")
+            search.enable_corr(bary_corr.flatten())
+
+        search_start = time.time()
+        print("Starting Search")
+        print("---------------------------------------")
+        param_headers = (
+            "Ecliptic Angle",
+            "Min. Search Angle",
+            "Max Search Angle",
+            "Min Velocity",
+            "Max Velocity",
+        )
+        param_values = (suggested_angle, *search_params["ang_lims"], *search_params["vel_lims"])
+        for header, val in zip(param_headers, param_values):
+            print("%s = %.4f" % (header, val))
+
+        # If we are using gpu_filtering, enable it and set the parameters.
+        if config["gpu_filter"]:
+            print("Using in-line GPU sigmaG filtering methods", flush=True)
+            coeff = post_process._find_sigmaG_coeff(config["sigmaG_lims"])
+            search.enable_gpu_sigmag_filter(
+                np.array(config["sigmaG_lims"]) / 100.0,
+                coeff,
+                config["lh_level"],
+            )
+
+        # If we are using an encoded image representation on GPU, enable it and
+        # set the parameters.
+        if config["encode_psi_bytes"] > 0 or config["encode_phi_bytes"] > 0:
+            search.enable_gpu_encoding(config["encode_psi_bytes"], config["encode_phi_bytes"])
+
+        # Enable debugging.
+        if config["debug"]:
+            search.set_debug(config["debug"])
+
+        search.search(
+            int(config["ang_arr"][2]),
+            int(config["v_arr"][2]),
+            *search_params["ang_lims"],
+            *search_params["vel_lims"],
+            int(config["num_obs"]),
+        )
+        print("Search finished in {0:.3f}s".format(time.time() - search_start), flush=True)
+        return search, search_params
diff --git a/src/kbmod/standardizer.py b/src/kbmod/standardizer.py
index 9a7a06662..c24a0e0a6 100644
--- a/src/kbmod/standardizer.py
+++ b/src/kbmod/standardizer.py
@@ -1,5 +1,7 @@
 import abc
 
+from kbmod.search import layered_image, raw_image, psf
+
 
 __all__ = ["Standardizer",]
 
@@ -41,6 +43,9 @@ class Standardizer(abc.ABC):
     multiple science exposures is expected to unravel that metadata and
     associate it with each of the science exposures individualy.
 
+    The data for which this standardization and unravelling will occur needs to
+    be appended to the items in ``processable`` attribute of theclass.
+
     This class is an abstract base class to serve as a recipe for implementing
     a Standardizer specialized for processing an particular dataset. Do not
     attempt to instantiate this class directly.
@@ -52,8 +57,15 @@ class Standardizer(abc.ABC):
 
     Attributes
     ----------
-    loc : `str`
+    location : `str`
         Location of the file being standardized.
+    processable : `list`
+        List of standardizers internal units that can be processed by KBMOD.
+        For example, for FITS files on local filesystem these could be the
+        AstroPy header units (`~astropy.io.fits.HDUList` elements), for Vera C.
+        Rubin their internal `Exposure*` objects etc. The number of processable
+        units will match the number of returned  standardized images, variances,
+        masks and PSFs.
     wcs : `list`
         WCSs associated with the processable image data.
     bbox : `list`
@@ -71,26 +83,6 @@ class Standardizer(abc.ABC):
     standardizers with low priority when processing an upload.
     """
 
-    @abc.abstractmethod
-    def __init__(self, location, *args, **kwargs):
-        self.location = location
-
-    def __init_subclass__(cls, **kwargs):
-        # Registers subclasses of this class with set `name` class
-        # parameters as availible standardizers. Note the specific
-        # standardizer has to be imported before this class since the
-        # mechanism is triggered at definition time.
-        name = getattr(cls, "name", False)
-        if name and name is not None:
-            super().__init_subclass__(**kwargs)
-            Standardizer.registry[cls.name] = cls
-
-    def __str__(self):
-        return f"{self.name}({self.location}, {self.exts})"
-
-    def __repr__(self):
-        return f"{self.__class__.__name__}({self.location})"
-
     @classmethod
     def get(cls, tgt=None, standardizer=None):
         """Get the standardizer class that can handle given file.
@@ -229,6 +221,70 @@ def canStandardize(self, tgt):
         """
         raise NotImplementedError()
 
+    def __init_subclass__(cls, **kwargs):
+        # Registers subclasses of this class with set `name` class
+        # parameters as availible standardizers. Note the specific
+        # standardizer has to be imported before this class since the
+        # mechanism is triggered at definition time.
+        name = getattr(cls, "name", False)
+        if name and name is not None:
+            super().__init_subclass__(**kwargs)
+            Standardizer.registry[cls.name] = cls
+
+    @abc.abstractmethod
+    def __init__(self, location, *args, **kwargs):
+        self.location = location
+
+    def __str__(self):
+        return f"{self.name}({self.location}, {self.exts})"
+
+    def __repr__(self):
+        return f"{self.__class__.__name__}({self.location})"
+
+    # I was a bit lazy and didn't feel like rekajiggering the whole code just
+    # to replace .exts with .processable - although it would save us some
+    # overhead.
+    @abc.abstractproperty
+    def processable(self):
+        """A list of processable units (i.e. images that can be processed by
+        KBMOD) extracted from the given target data source."""
+        raise NotImplementedError()
+
+    # all these should probably be named in plural - but lord allmighty that's
+    # a lot of references to update and it feels so unnatural - help?
+    @abc.abstractproperty
+    def wcs(self):
+        """A list of WCS's or `None` for each entry marked as processable.
+
+        Expected to be an WCS object, but when it is not possible to construct
+        one `None` will be returned in its place.
+        """
+        raise NotImplementedError()
+
+    @abc.abstractproperty
+    def bbox(self):
+        """A list of bounding boxes or `None` for each entry marked as
+        processable.
+
+        ========== ============================================================
+        Key        Description
+        ========== ============================================================
+        center_ra  On-sky (ICRS) right ascension coordinate, expressed in
+                   decimal degrees, of the center of an rectangular bounding
+                   box.
+        center_dec On-sky (ICRS) declination coordinate, expressed in decimal
+                   degrees, of the center of an rectangular bounding box.
+        corner_ra  On-sky (ICRS) right ascension coordinate, expressed in
+                   decimal degrees, of an corner of a rectangular bounding box.
+        corner_dec On-sky (ICRS) declination coordinate, expressed in decimal
+                   degrees, of an corner of a rectangular bounding box.
+        ========== ============================================================
+
+        Often, when WCS can not be constructed, the BBOX can not be
+        constructed. Then `None` is expected to be returned in it's place.
+        """
+        raise NotImplementedError()
+
     @abc.abstractmethod
     def standardizeWCS(self):
         """Creates an WCS for every processable science image.
@@ -263,7 +319,7 @@ def standardizeMetadata(self):
         """Standardizes required and optional metadata from the given data.
 
         The standardized metadata is expected to be a dict containing the
-        following keys in data:
+        following keys and values:
 
         ======== ==============================================================
         Key      Description
@@ -286,6 +342,19 @@ def standardizeMetadata(self):
         -------
         metadata : `dict`
             Standardized metadata.
+
+
+        Notes
+        -----
+        If a target being standardized contains only 1 image that can be
+        standardized and processed by KBMOD, the values in the dictionary can
+        be strings, integers, floats, and other such single-values non-iterable
+        types, or an iterable type containing only 1 element.
+        When the target being standardized contains multiple images/extensions
+        that are processable by KBMOD the single-valued dict entries will be
+        considered as shared by all the processable images/extension/units of
+        the target. The iterable values in the dictionary must match the number
+        and order of the units marked as processable.
         """
         raise NotImplementedError()
 
@@ -299,8 +368,8 @@ def standardizeScienceImage(self):
 
         Returns
         -------
-        image : `np.array`
-            Science image.
+        image : `list[~np.array]`
+            Science images.
         """
         raise NotImplemented()
 
@@ -320,8 +389,8 @@ def standardizeVarianceImage(self):
 
         Returns
         -------
-        variance : `np.array`
-            Variance image.
+        variance : `list[`np.array]`
+            Variance images.
         """
         raise NotImplementedError()
 
@@ -337,50 +406,50 @@ def standardizeMaskImage(self):
 
         Returns
         -------
-        mask : `np.array`
-            Mask image.
+        mask : `list[~np.array]`
+            Mask images.
         """
         raise NotImplementedError()
 
     # no idea really what to do bout this one? AFAIK almost no images come
     # with PSFs in them
-    #@abc.abstractmethod
-    #def standardizePSF(self):
-    #    raise NotImplementedError()
-
-    def standardize(self):
-        """Invokes standardize metadata, image, variance, mask and PSF and
-        returns the results as a dictionary.
+    @abc.abstractmethod
+    def standardizePSF(self):
+        """Returns PSF for each extension marked as processable.
 
         Returns
         -------
-        standardizedData : `dict`
-            Dictionary with standardized ``meta``, ``science``, ``variance``,
-            ``mask`` and ``psf`` values.
+        psfs : `list[~kbmod.search.psf]`
+            List of `~kbmod.search.psf` objects.
         """
-        std = {"meta": self.standardizeStddata()}
-        std.update({"science": self.standardizeScienceImage()})
-        std.update({"variance": self.standardizeVarianceImage()})
-        std.update({"mask": self.standardizeMaskImage()})
-        std.update({"psf": self.standardizePSF()})
+        raise NotImplementedError()
 
-        return std
 
+    @abc.abstractmethod
     def toLayeredImage(self):
         """Run metadata standardization methods. These include header
         and bounding box standardization.
 
         Notes
         -----
-        Implementation is processor-specific.
+        Implementation is standardizer-specific.
+        """
+        raise NotImplementedError()
+
+    def standardize(self):
+        """Invokes standardize metadata, image, variance, mask and PSF and
+        returns the results as a dictionary.
+
+        Returns
+        -------
+        standardizedData : `dict`
+            Dictionary with standardized ``meta``, ``science``, ``variance``,
+            ``mask`` and ``psf`` values.
         """
-        # note that these will be lists so we need a new constructor
-        # They will be equal length(maybe?) and each is a different
-        # detector - i.e. not an image stack. Because we don't use
-        # ndarray this will be a copy operation and not very fast.
-        header = self.standardizeHeader()
-        masks = self.standardizeMask()
-        variances = self.standardizeVariance()
-        #psfs = self.standardizePSF()
-
-        return LayeredImage(self.loc, self.exts, masks, variances, psfs, header["mjd"])
+        std = {"meta": self.standardizeMetadata()}
+        std["science"] = self.standardizeScienceImage()
+        std["variance"] = self.standardizeVarianceImage()
+        std["mask"] = self.standardizeMaskImage()
+        std["psf"] = self.standardizePSF()
+
+        return std
diff --git a/src/kbmod/standardizers/butler_standardizer.py b/src/kbmod/standardizers/butler_standardizer.py
index d81e608c4..17b2ce8eb 100644
--- a/src/kbmod/standardizers/butler_standardizer.py
+++ b/src/kbmod/standardizers/butler_standardizer.py
@@ -1,6 +1,12 @@
+import astropy.nddata as bitmask
 from astropy.wcs import WCS
 
+import numpy as np
+
+from scipy.signal import convolve2d
+
 from kbmod.standardizer import Standardizer
+from kbmod.search import layered_image, raw_image, psf
 
 
 __all__ = ["ButlerStandardizer"]
@@ -28,16 +34,58 @@ def canStandardize(self, tgt):
             return True, tgt
         return False, []
 
-    def __init__(self, butler, datasetRefs):
+    def __init__(self, butler, datasetRefs=None, id=None, **kwargs):
+        # this requires an invovled logic sequence where
+        # 1) we need to check if butler is Butler, datasetRefs a list of
+        #    DatasetRef - proceed as is
+        # 2) if butler is a string - instantiate butler
+        # 3) if the datRefs are a list of DatasetRefs - great
+        # 4) else construct DatasetRef's
+        #    4.1) They are DataId objects - butler.registry.getDataset(id)
+        #    4.2) if the datRefs are strings and UUIDs - b.r.getDataset(id)
+        #    4.3) mappings must specify datasetType, instrument, detector,
+        #         visit, and collection, then:
+        #         butler.registry.queryDatasets(datasetType='goodSeeingDiff_differenceTempExp',
+        #                                       dataId= {"instrument":"DECam", "detector":35, "visit":974895},
+        #                                       collections=["collection"/datref.run, ])
+        #        I guess majority of these we can put in our table by default,
+        #        certainly the number of columns we need is growing compared to
+        #        metadata columns for the user - maybe separate into two tables
+        #        and then flatten before writing and serialize the metadata
+        #        about which columns go where into `meta`, but this is getting
+        #        a bit too far atm, let's just raise an error in _get_std for now
+        # Another question to answer here is the construction and lazy loading
+        # of standardizers. FitsStandardizers have 1 input init param - location
+        # so in _get_std in image collection we just pass that. A more generic
+        # way would be to pass the whole row. The issue here is that __init__
+        # is now having to deal with handling a mapping and that makes it less
+        # obvious for the end user. Optionally, we can give a mandatory
+        # **kwarg to every __init__ and rely on column naming and unpacking.
+        # this is not the best, explicit better than implicit, but also now
+        # the __init__ method kwargs need to be named the same as columns and
+        # be optional too. For example does this init make sense to you?
+        # I don't think I would be able to write this without knowing how the
+        # whole thing works...
+        #breakpoint()
         self.butler = butler
-        self.refs = datasetRefs
+        self.location = butler.datastore.root
+        if datasetRefs is None:
+            # import this here for now, for optimization
+            from lsst.daf.butler.core import DatasetId
+            self.refs = [butler.registry.getDataset(DatasetId(id)),]
+        else:
+            self.refs = datasetRefs
         # wth? why is it a data reference if it doesn't reference data without
         # being explicit about it? Why `collection` and `collections` don't
         # behave the same way? Why is this a thing! AAARRRGGGHH!
-        self.exts = [butler.get(ref, collections=[ref.run, ]) for ref in datasetRefs]
+        self.exts = [butler.get(ref, collections=[ref.run, ]) for ref in self.refs]
         self._wcs = []
         self._bbox = []
 
+    @property
+    def processable(self):
+        return self.exts
+
     @property
     def wcs(self):
         if not self._wcs:
@@ -113,7 +161,7 @@ def _maskSingleExp(self, exp):
         )
 
         brigthness_threshold = exp.image.array.mean() - exp.image.array.std()
-        threshold_mask = exp.image.array > brightness_threshold
+        threshold_mask = exp.image.array > brigthness_threshold
 
         net_mask = bit_mask & threshold_mask
 
@@ -160,6 +208,7 @@ def standardizeMetadata(self):
 
         metadata["mjd"] = [e.visitInfo.date.toAstropy().mjd for e in self.exts]
         metadata["filter"] = [e.info.getFilter().physicalLabel for e in self.exts]
+        metadata["id"] = [str(r.id) for r in self.refs]
         metadata["exp_id"] = [e.info.id for e in self.exts]
 
         # it's also like super dificult to extract any information out of the
@@ -196,3 +245,21 @@ def standardizeVarianceImage(self):
 
     def standardizeMaskImage(self):
         return (self._maskSingleExp(exp) for exp in self.exts)
+
+    def standardizePSF(self):
+        return (psf(1) for e in self.exts)
+
+    def toLayeredImage(self):
+        meta = self.standardizeMetadata()
+        sciences = self.standardizeScienceImage()
+        variances = self.standardizeVarianceImage()
+        masks = self.standardizeMaskImage()
+
+        psfs = self.standardizePSF()
+
+        # guaranteed to exist, i.e. safe to access
+        mjds = meta["mjd"]
+        imgs = []
+        for sci, var, mask, psf, t in zip(sciences, variances, masks, psfs, mjds):
+            imgs.append(layered_image(raw_image(sci), raw_image(var), raw_image(mask), t, psf))
+        return imgs
diff --git a/src/kbmod/standardizers/decam_community_pipe_fits.py b/src/kbmod/standardizers/decam_community_pipe_fits.py
index a6d23631f..69778d15e 100644
--- a/src/kbmod/standardizers/decam_community_pipe_fits.py
+++ b/src/kbmod/standardizers/decam_community_pipe_fits.py
@@ -3,6 +3,11 @@
 
 from astropy.stats import sigma_clipped_stats
 from astropy.io import fits
+from astropy.nddata import bitmask
+
+import numpy as np
+
+from scipy.signal import convolve2d
 
 from .multi_extension_fits import MultiExtensionFits
 from astro_metadata_translator import ObservationInfo
@@ -84,8 +89,8 @@ def _getScienceImages(cls, hdulist):
 
         return exts
 
-    def __init__(self, location):
-        super().__init__(location)
+    def __init__(self, location, **kwargs):
+        super().__init__(location, **kwargs)
 
         # self.exts is filled in super() but that could include tracking and
         # focus chips too (as they fall under image-like conditions). We want
diff --git a/src/kbmod/standardizers/fits_standardizer.py b/src/kbmod/standardizers/fits_standardizer.py
index b7c04ad56..61261f763 100644
--- a/src/kbmod/standardizers/fits_standardizer.py
+++ b/src/kbmod/standardizers/fits_standardizer.py
@@ -6,6 +6,7 @@
 from astropy.wcs import WCS
 
 from kbmod.standardizer import Standardizer
+from kbmod.search import layered_image, raw_image, psf
 
 
 __all__ = ["FitsStandardizer",]
@@ -50,7 +51,6 @@ class FitsStandardizer(Standardizer):
     """File extensions this processor can handle."""
 
     @classmethod
-    @abstractmethod
     def canStandardize(cls, tgt):
         # docstring inherited from Standardizer; TODO: check it's True
         canProcess, hdulist = False, None
@@ -85,6 +85,10 @@ def __init__(self, location):
         self._wcs = []
         self._bbox = []
 
+    @property
+    def processable(self):
+        return self.exts
+
     @property
     def wcs(self):
         if not self._wcs:
@@ -237,7 +241,7 @@ def translateHeader(self, *args, **kwargs):
 
     def standardizeMetadata(self):
         metadata = self.translateHeader()
-        metadata.update({"location": self.location})
+        metadata["location"] = self.location
         metadata.update({"wcs": self.wcs, "bbox": self.bbox})
 
         # calculate the pointing from the bbox or wcs if they exist?
@@ -262,4 +266,40 @@ def standardizeMetadata(self):
         return metadata
 
     def standardizeScienceImage(self):
+        # the assumption here is that all Exts are AstroPy HDU objects
         return (ext.data for ext in self.exts)
+
+    def standardizePSF(self):
+        return (psf(1) for e in self.exts)
+
+    def toLayeredImage(self):
+        """Returns a list of `~kbmod.search.layered_image` objects for each
+        entry marked as processable.
+
+        Returns
+        -------
+        layeredImage : `list`
+            Layered image objects.
+        """
+        meta = self.standardizeMetadata()
+        sciences = self.standardizeScienceImage()
+        variances = self.standardizeVarianceImage()
+        masks = self.standardizeMaskImage()
+
+        psfs = self.standardizePSF()
+
+        # guaranteed to exist, i.e. safe to access, but isn't unravelled here
+        # or potentially it is - we can't tell?
+        if isinstance(meta["mjd"], (list, tuple)):
+            mjds = meta["mjd"]
+        else:
+            mjds = (meta["mjd"] for e in self.exts)
+
+        imgs = []
+        for sci, var, mask, psf, t in zip(sciences, variances, masks, psfs, mjds):
+            imgs.append(layered_image(raw_image(sci), raw_image(var), raw_image(mask), t, psf))
+
+        return imgs
+
+
+
diff --git a/src/kbmod/standardizers/multi_extension_fits.py b/src/kbmod/standardizers/multi_extension_fits.py
index 50927199e..9d47867bb 100644
--- a/src/kbmod/standardizers/multi_extension_fits.py
+++ b/src/kbmod/standardizers/multi_extension_fits.py
@@ -95,7 +95,7 @@ def canStandardize(cls, location):
         canStandardize = parentCanStandardize and len(hdulist) > 1
         return canStandardize, hdulist
 
-    def __init__(self, location, set_exts=False):
+    def __init__(self, location, set_exts=False, **kwargs):
         super().__init__(location)
 
         if set_exts:
diff --git a/src/kbmod/standardizers/rubin_scipipe_fits.py b/src/kbmod/standardizers/rubin_scipipe_fits.py
index 8b7014cc5..7ca674401 100644
--- a/src/kbmod/standardizers/rubin_scipipe_fits.py
+++ b/src/kbmod/standardizers/rubin_scipipe_fits.py
@@ -7,6 +7,7 @@
 from astropy.time import Time
 from astropy.nddata import bitmask
 
+import numpy as np
 from scipy.signal import convolve2d
 
 from kbmod.standardizers import MultiExtensionFits
@@ -64,8 +65,8 @@ def canStandardize(cls, location):
         canStandardize = parentCanStandardize and isRubin
         return canStandardize, hdulist
 
-    def __init__(self, location):
-        super().__init__(location)
+    def __init__(self, location, **kwargs):
+        super().__init__(location, **kwargs)
         # this is the only science-image header for Rubin
         self.exts = [self.hdulist[1], ]
 
@@ -86,7 +87,6 @@ def translateHeader(self):
         return standardizedHeader
 
     def standardizeMaskImage(self):
-
         idx = self.hdulist.index_of("MASK")
         bitfield = self.hdulist[idx].data
         bit_mask = bitmask.bitfield_to_boolean_mask(
@@ -97,7 +97,7 @@ def standardizeMaskImage(self):
 
         idx = self.hdulist.index_of("IMAGE")
         image = self.hdulist[idx].data
-        brigthness_threshold = image.mean() - image.std()
+        brightness_threshold = image.mean() - image.std()
         threshold_mask = image > brightness_threshold
 
         net_mask = threshold_mask & bit_mask
@@ -106,8 +106,7 @@ def standardizeMaskImage(self):
         grow_kernel = np.ones((11, 11))
         grown_mask = convolve2d(net_mask, grow_kernel, mode="same")
 
-        return grown_mask
+        return [grown_mask, ]
 
     def standardizeVarianceImage(self):
-        idx = self.hdulist.index_of("VARIANCE")
-        return self.hdulist[idx].data
+        return [self.hdulist["VARIANCE"].data, ]