Merge branch 'main' into remove_legacy_visualization

src/kbmod/image_collection.py

@@ -270,6 +270,10 @@ def fromTargets(cls, tgts, force=None, config=None, **kwargs):
             of the registered standardizers can be provided. Optionally,
             provide the `Standardizer` class itself in which case it will be
             called for each target in the iterable.
+        config : `~StandardizerConfig`, `dict` or `None`, optional
+            Standardizer configuration or dictionary containing the config
+            parameters for standardization. When `None` default values for the
+            appropriate `Standardizer` will be used.
         **kwargs : `dict`
             Remaining keyword arguments are passed to the `Standardizer`.
 
@@ -296,10 +300,14 @@ def fromDir(cls, dirpath, recursive=False, force=None, config=None, **kwargs):
         recursive : `bool`
             If the location is a local filesystem directory, scan it
             recursively including all sub-directories.
-        forceStandardizer : `Standardizer` or `None`
+        force : `Standardizer` or `None`
             If `None`, when applicable, determine the correct `Standardizer` to
             use automatically. Otherwise force the use of the given
             `Standardizer`.
+        config : `~StandardizerConfig`, `dict` or `None`, optional
+            Standardizer configuration or dictionary containing the config
+            parameters for standardization. When `None` default values for the
+            appropriate `Standardizer` will be used.
         **kwargs : `dict`
             Remaining kwargs, not listed here, are passed onwards to
             the underlying `Standardizer`.

src/kbmod/search/layered_image.cpp

@@ -2,8 +2,9 @@
 
 namespace search {
 
-LayeredImage::LayeredImage(const RawImage& sci, const RawImage& var, const RawImage& msk, const PSF& psf)
-        : psf(psf) {
+LayeredImage::LayeredImage(const RawImage& sci, const RawImage& var, const RawImage& msk, const PSF& psf,
+                           double obs_time)
+        : obstime(obs_time), psf(psf) {
     // Get the dimensions of the science layer and check for consistency with
     // the other two layers.
     width = sci.get_width();
@@ -20,9 +21,10 @@ LayeredImage::LayeredImage(const RawImage& sci, const RawImage& var, const RawIm
 }
 
 LayeredImage::LayeredImage(Image& sci, Image& var, Image& msk, PSF& psf, double obs_time)
-        : science(sci, obs_time), variance(var, obs_time), mask(msk, obs_time), psf(psf) {
+        : science(sci), variance(var), mask(msk), psf(psf) {
     width = science.get_width();
     height = science.get_height();
+    obstime = obs_time;
     assert_sizes_equal(variance.get_width(), width, "variance layer width");
     assert_sizes_equal(variance.get_height(), height, "variance layer height");
     assert_sizes_equal(mask.get_width(), width, "mask layer width");
@@ -33,6 +35,7 @@ LayeredImage::LayeredImage(Image& sci, Image& var, Image& msk, PSF& psf, double
 LayeredImage::LayeredImage(const LayeredImage& source) noexcept {
     width = source.width;
     height = source.height;
+    obstime = source.obstime;
     science = source.science;
     mask = source.mask;
     variance = source.variance;
@@ -43,6 +46,7 @@ LayeredImage::LayeredImage(const LayeredImage& source) noexcept {
 LayeredImage::LayeredImage(LayeredImage&& source) noexcept
         : width(source.width),
           height(source.height),
+          obstime(source.obstime),
           science(std::move(source.science)),
           mask(std::move(source.mask)),
           variance(std::move(source.variance)),
@@ -52,6 +56,7 @@ LayeredImage::LayeredImage(LayeredImage&& source) noexcept
 LayeredImage& LayeredImage::operator=(const LayeredImage& source) noexcept {
     width = source.width;
     height = source.height;
+    obstime = source.obstime;
     science = source.science;
     mask = source.mask;
     variance = source.variance;
@@ -64,6 +69,7 @@ LayeredImage& LayeredImage::operator=(LayeredImage&& source) noexcept {
     if (this != &source) {
         width = source.width;
         height = source.height;
+        obstime = source.obstime;
         science = std::move(source.science);
         mask = std::move(source.mask);
         variance = std::move(source.variance);
@@ -252,7 +258,8 @@ static void layered_image_bindings(py::module& m) {
     using pf = search::PSF;
 
     py::class_<li>(m, "LayeredImage", pydocs::DOC_LayeredImage)
-            .def(py::init<const ri&, const ri&, const ri&, pf&>())
+            .def(py::init<const ri&, const ri&, const ri&, const pf&, double>(), py::arg("sci"),
+                 py::arg("var"), py::arg("msk"), py::arg("psf"), py::arg("obs_time") = -1.0)
             .def(py::init<search::Image&, search::Image&, search::Image&, pf&, double>(),
                  py::arg("sci").noconvert(true), py::arg("var").noconvert(true),
                  py::arg("msk").noconvert(true), py::arg("psf"), py::arg("obs_time"))

src/kbmod/search/layered_image.h

@@ -14,7 +14,8 @@
 namespace search {
 class LayeredImage {
 public:
-    explicit LayeredImage(const RawImage& sci, const RawImage& var, const RawImage& msk, const PSF& psf);
+    explicit LayeredImage(const RawImage& sci, const RawImage& var, const RawImage& msk, const PSF& psf,
+                          double obs_time = -1.0);
 
     // Build a layered image from the underlying matrices, taking ownership of the image data.
     explicit LayeredImage(Image& sci, Image& var, Image& msk, PSF& psf, double obs_time);
@@ -33,8 +34,8 @@ class LayeredImage {
     unsigned get_width() const { return width; }
     unsigned get_height() const { return height; }
     uint64_t get_npixels() const { return width * height; }
-    double get_obstime() const { return science.get_obstime(); }
-    void set_obstime(double obstime) { science.set_obstime(obstime); }
+    double get_obstime() const { return obstime; }
+    void set_obstime(double new_obstime) { obstime = new_obstime; }
 
     // Getter functions for the data in the individual layers.
     RawImage& get_science() { return science; }
@@ -92,6 +93,7 @@ class LayeredImage {
 private:
     unsigned width;
     unsigned height;
+    double obstime;
 
     PSF psf;
     RawImage science;

src/kbmod/search/pydocs/layered_image_docs.h

@@ -15,6 +15,8 @@ static const auto DOC_LayeredImage = R"doc(
       The `RawImage` for the mask layer.
   psf : `PSF`
       The PSF for the image.
+  obstime : `float`
+      The time of the image.
 
   Raises
   ------

src/kbmod/search/pydocs/raw_image_docs.h

@@ -5,23 +5,19 @@ namespace pydocs {
 static const auto DOC_RawImage = R"doc(
   Image and the time it was observed at.
 
-  Instantiated from a pair of obstime and image, or from image dimensions and
-  obstime. When instantiating from image dimensions and obstime, it's possible
-  to provide a default value used to fill the array with. Otherwise the array is
-  filled with zeros.
+  Instantiated from an image or from image dimensions and and a value (which
+  defaults to zero).
 
   Parameters
   ----------
   image : `numpy.array`, optional
       Image, row-major a 2D array. The array *must* be of dtype `numpy.single`.
-  obstime : `float`, optional
-      MJD stamp, time the image was observed at.
   width : `int`, optional
       Width of the image.
   height : `int`, optional
       Height of the image.
   value : `float`, optional
-      When instantiated from dimensions and obstime, value that fills the array.
+      When instantiated from dimensions, value that fills the array.
       Default is 0.
 
   Attributes
@@ -32,8 +28,6 @@ static const auto DOC_RawImage = R"doc(
       Image width, in pixels.
   npixels : `int`
       Number of pixels in the image, equivalent to ``width*height``.
-  obstime : `float`
-      MJD time of observation.
   image : `np.array[np,single]`
       Image array.
 
@@ -61,7 +55,7 @@ static const auto DOC_RawImage = R"doc(
   >>> from kbmod.search import RawImage
   >>> import numpy as np
   >>> ri = RawImage()
-  >>> ri = RawImage(w=2, h=3, value=1, obstime=10)
+  >>> ri = RawImage(w=2, h=3, value=1)
   >>> ri.image
   array([[1., 1.],
        [1., 1.],

src/kbmod/search/raw_image.cpp

@@ -4,17 +4,16 @@ namespace search {
 using Index = indexing::Index;
 using Point = indexing::Point;
 
-RawImage::RawImage() : width(0), height(0), obstime(-1.0), image() {}
+RawImage::RawImage() : width(0), height(0), image() {}
 
-RawImage::RawImage(Image& img, double obs_time) {
+RawImage::RawImage(Image& img) {
     image = std::move(img);
     height = image.rows();
     width = image.cols();
-    obstime = obs_time;
 }
 
-RawImage::RawImage(unsigned w, unsigned h, float value, double obs_time)
-        : width(w), height(h), obstime(obs_time) {
+RawImage::RawImage(unsigned w, unsigned h, float value)
+        : width(w), height(h) {
     if (value != 0.0f)
         image = Image::Constant(height, width, value);
     else
@@ -26,22 +25,19 @@ RawImage::RawImage(const RawImage& old) noexcept {
     width = old.get_width();
     height = old.get_height();
     image = old.get_image();
-    obstime = old.get_obstime();
 }
 
 // Move constructor
 RawImage::RawImage(RawImage&& source) noexcept
         : width(source.width),
           height(source.height),
-          obstime(source.obstime),
           image(std::move(source.image)) {}
 
 // Copy assignment
 RawImage& RawImage::operator=(const RawImage& source) noexcept {
     width = source.width;
     height = source.height;
     image = source.image;
-    obstime = source.obstime;
     return *this;
 }
 
@@ -51,7 +47,6 @@ RawImage& RawImage::operator=(RawImage&& source) noexcept {
         width = source.width;
         height = source.height;
         image = std::move(source.image);
-        obstime = source.obstime;
     }
     return *this;
 }
@@ -484,15 +479,13 @@ static void raw_image_bindings(py::module& m) {
     py::class_<rie>(m, "RawImage", pydocs::DOC_RawImage)
             .def(py::init<>())
             .def(py::init<search::RawImage&>())
-            .def(py::init<search::Image&, double>(), py::arg("img").noconvert(true),
-                 py::arg("obs_time") = -1.0d)
-            .def(py::init<unsigned, unsigned, float, double>(), py::arg("w"), py::arg("h"),
-                 py::arg("value") = 0.0f, py::arg("obs_time") = -1.0d)
+            .def(py::init<search::Image&>(), py::arg("img").noconvert(true))
+            .def(py::init<unsigned, unsigned, float>(), py::arg("w"), py::arg("h"),
+                 py::arg("value") = 0.0f)
             // attributes and properties
             .def_property_readonly("height", &rie::get_height)
             .def_property_readonly("width", &rie::get_width)
             .def_property_readonly("npixels", &rie::get_npixels)
-            .def_property("obstime", &rie::get_obstime, &rie::set_obstime)
             .def_property("image", py::overload_cast<>(&rie::get_image, py::const_), &rie::set_image)
             .def_property("imref", py::overload_cast<>(&rie::get_image), &rie::set_image)
             // pixel accessors and setters

src/kbmod/search/raw_image.h

@@ -26,8 +26,8 @@ using ImageIRef = Eigen::Ref<Image>;
 class RawImage {
 public:
     explicit RawImage();
-    explicit RawImage(Image& img, double obs_time = -1.0);
-    explicit RawImage(unsigned w, unsigned h, float value = 0.0, double obs_time = -1.0);
+    explicit RawImage(Image& img);
+    explicit RawImage(unsigned w, unsigned h, float value = 0.0);
 
     RawImage(const RawImage& old) noexcept;  // Copy constructor
     RawImage(RawImage&& source) noexcept;    // Move constructor
@@ -70,10 +70,6 @@ class RawImage {
 
     void replace_masked_values(float value = 0.0);
 
-    // Functions for locally storing the image time.
-    double get_obstime() const { return obstime; }
-    void set_obstime(double new_time) { obstime = new_time; }
-
     // this will be a raw pointer to the underlying array
     // we use this to copy to GPU and nowhere else!
     float* data() { return image.data(); }
@@ -129,7 +125,6 @@ class RawImage {
 private:
     unsigned width;
     unsigned height;
-    double obstime;
     Image image;
 };
 

src/kbmod/search/stamp_creator.cpp

@@ -14,7 +14,7 @@ std::vector<RawImage> StampCreator::create_stamps(ImageStack& stack, const Traje
                                                   bool keep_no_data, const std::vector<bool>& use_index) {
     if (use_index.size() > 0)
         assert_sizes_equal(use_index.size(), stack.img_count(), "create_stamps() use_index");
-    bool use_all_stamps = use_index.size() == 0;
+    bool use_all_stamps = (use_index.size() == 0);
 
     std::vector<RawImage> stamps;
     unsigned int num_times = stack.img_count();
@@ -65,8 +65,8 @@ std::vector<RawImage> StampCreator::get_coadded_stamps(ImageStack& stack, std::v
                                                        std::vector<std::vector<bool>>& use_index_vect,
                                                        const StampParameters& params, bool use_gpu) {
     logging::Logger* rs_logger = logging::getLogger("kbmod.search.stamp_creator");
-    rs_logger->info("Performing stamp filtering on " + std::to_string(t_array.size()) + " trajectories.");
-    DebugTimer timer = DebugTimer("stamp filtering", rs_logger);
+    rs_logger->info("Generating co_added stamps on " + std::to_string(t_array.size()) + " trajectories.");
+    DebugTimer timer = DebugTimer("coadd generating", rs_logger);
 
     // We use the GPU if we have it for everything except STAMP_VAR_WEIGHTED which is CPU only.
     if (use_gpu && (params.stamp_type != STAMP_VAR_WEIGHTED)) {
@@ -277,7 +277,7 @@ std::vector<RawImage> StampCreator::create_variance_stamps(ImageStack& stack, co
                                                            int radius, const std::vector<bool>& use_index) {
     if (use_index.size() > 0)
         assert_sizes_equal(use_index.size(), stack.img_count(), "create_stamps() use_index");
-    bool use_all_stamps = use_index.size() == 0;
+    bool use_all_stamps = (use_index.size() == 0);
 
     std::vector<RawImage> stamps;
     unsigned int num_times = stack.img_count();
@@ -297,22 +297,26 @@ std::vector<RawImage> StampCreator::create_variance_stamps(ImageStack& stack, co
 
 RawImage StampCreator::get_variance_weighted_stamp(ImageStack& stack, const Trajectory& trj, int radius,
                                                    const std::vector<bool>& use_index) {
+    if (radius < 0) throw std::runtime_error("Invalid stamp radius. Must be >= 0.");
     unsigned int num_images = stack.img_count();
     if (num_images == 0) throw std::runtime_error("Unable to create mean image given 0 images.");
     unsigned int stamp_width = 2 * radius + 1;
 
     // Make the stamps for each time step.
-    std::vector<bool> empty_vect;
     std::vector<RawImage> sci_stamps = create_stamps(stack, trj, radius, true /*=keep_no_data*/, use_index);
     std::vector<RawImage> var_stamps = create_variance_stamps(stack, trj, radius, use_index);
+    if (sci_stamps.size() != var_stamps.size()) {
+        throw std::runtime_error("Mismatched number of stamps returned.");
+    }
+    num_images = sci_stamps.size();
 
     // Do the weighted mean.
     Image result = Image::Zero(stamp_width, stamp_width);
-    for (int y = 0; y < stamp_width; ++y) {
-        for (int x = 0; x < stamp_width; ++x) {
+    for (unsigned int y = 0; y < stamp_width; ++y) {
+        for (unsigned int x = 0; x < stamp_width; ++x) {
             float sum = 0.0;
             float scale = 0.0;
-            for (int i = 0; i < num_images; ++i) {
+            for (unsigned int i = 0; i < num_images; ++i) {
                 float sci_val = sci_stamps[i].get_pixel({y, x});
                 float var_val = var_stamps[i].get_pixel({y, x});
                 if (pixel_value_valid(sci_val) && pixel_value_valid(var_val) && (var_val != 0.0)) {

src/kbmod/standardizers/fits_standardizers/fits_standardizer.py

@@ -410,7 +410,11 @@ def toLayeredImage(self):
         # copy. TODO: fix when/if CPP stuff is fixed.
         imgs = []
         for sci, var, mask, psf, t in zip(sciences, variances, masks, psfs, mjds):
+            # Make sure the science and variance layers are float32.
+            sci = sci.astype(np.float32)
+            var = var.astype(np.float32)
+
             # Converts nd array mask from bool to np.float32
             mask = mask.astype(np.float32)
-            imgs.append(LayeredImage(RawImage(sci, t), RawImage(var, t), RawImage(mask, t), psf))
+            imgs.append(LayeredImage(RawImage(sci), RawImage(var), RawImage(mask), psf, obs_time=t))
         return imgs