From b6ffa344e0c4f4c280d48ff78fce764e36df68d7 Mon Sep 17 00:00:00 2001 From: Jeremy Kubica <104161096+jeremykubica@users.noreply.github.com> Date: Fri, 20 Sep 2024 09:02:54 -0400 Subject: [PATCH 1/3] Move obstime from RawImage to LayeredImage --- src/kbmod/file_utils.py | 4 -- src/kbmod/search/layered_image.cpp | 15 +++++-- src/kbmod/search/layered_image.h | 8 ++-- src/kbmod/search/pydocs/layered_image_docs.h | 2 + src/kbmod/search/pydocs/raw_image_docs.h | 14 ++---- src/kbmod/search/raw_image.cpp | 21 +++------ src/kbmod/search/raw_image.h | 9 +--- .../fits_standardizers/fits_standardizer.py | 2 +- src/kbmod/work_unit.py | 44 ++++++------------- tests/test_butlerstd.py | 3 -- tests/test_fake_data_creator.py | 2 +- tests/test_psi_phi_array.py | 8 ++-- tests/test_raw_image.py | 26 ++++------- tests/test_standardizer.py | 3 -- tests/test_work_unit.py | 2 +- 15 files changed, 60 insertions(+), 103 deletions(-) diff --git a/src/kbmod/file_utils.py b/src/kbmod/file_utils.py index 86e59b45..b8297c0a 100644 --- a/src/kbmod/file_utils.py +++ b/src/kbmod/file_utils.py @@ -75,10 +75,6 @@ class FileUtils: Examples -------- - * Load an external file of visit ID to timestamp mappings. - - ``time_dict = FileUtils.load_time_dictionary("kbmod/data/demo_times.dat")`` - * Load the results of a KBMOD run as trajectory objects. ``FileUtils.load_results_file_as_trajectories("results_DEMO.txt")`` diff --git a/src/kbmod/search/layered_image.cpp b/src/kbmod/search/layered_image.cpp index 50baaeda..d07f1800 100644 --- a/src/kbmod/search/layered_image.cpp +++ b/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_
  • (m, "LayeredImage", pydocs::DOC_LayeredImage) - .def(py::init()) + .def(py::init(), py::arg("sci"), + py::arg("var"), py::arg("msk"), py::arg("psf"), py::arg("obs_time") = -1.0) .def(py::init(), py::arg("sci").noconvert(true), py::arg("var").noconvert(true), py::arg("msk").noconvert(true), py::arg("psf"), py::arg("obs_time")) diff --git a/src/kbmod/search/layered_image.h b/src/kbmod/search/layered_image.h index d30127c8..292e3135 100644 --- a/src/kbmod/search/layered_image.h +++ b/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; diff --git a/src/kbmod/search/pydocs/layered_image_docs.h b/src/kbmod/search/pydocs/layered_image_docs.h index 6691d49d..70400bdc 100644 --- a/src/kbmod/search/pydocs/layered_image_docs.h +++ b/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 ------ diff --git a/src/kbmod/search/pydocs/raw_image_docs.h b/src/kbmod/search/pydocs/raw_image_docs.h index 7a108a95..94bae338 100644 --- a/src/kbmod/search/pydocs/raw_image_docs.h +++ b/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.], diff --git a/src/kbmod/search/raw_image.cpp b/src/kbmod/search/raw_image.cpp index 6b1df113..c27616cb 100644 --- a/src/kbmod/search/raw_image.cpp +++ b/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,14 +25,12 @@ 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 @@ -41,7 +38,6 @@ 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_(m, "RawImage", pydocs::DOC_RawImage) .def(py::init<>()) .def(py::init()) - .def(py::init(), py::arg("img").noconvert(true), - py::arg("obs_time") = -1.0d) - .def(py::init(), py::arg("w"), py::arg("h"), - py::arg("value") = 0.0f, py::arg("obs_time") = -1.0d) + .def(py::init(), py::arg("img").noconvert(true)) + .def(py::init(), 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 diff --git a/src/kbmod/search/raw_image.h b/src/kbmod/search/raw_image.h index 42c49bd3..c135761c 100644 --- a/src/kbmod/search/raw_image.h +++ b/src/kbmod/search/raw_image.h @@ -26,8 +26,8 @@ using ImageIRef = Eigen::Ref; 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; }; diff --git a/src/kbmod/standardizers/fits_standardizers/fits_standardizer.py b/src/kbmod/standardizers/fits_standardizers/fits_standardizer.py index dfce61ad..28efbe11 100644 --- a/src/kbmod/standardizers/fits_standardizers/fits_standardizer.py +++ b/src/kbmod/standardizers/fits_standardizers/fits_standardizer.py @@ -412,5 +412,5 @@ def toLayeredImage(self): for sci, var, mask, psf, t in zip(sciences, variances, masks, psfs, mjds): # 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 diff --git a/src/kbmod/work_unit.py b/src/kbmod/work_unit.py index bf9e932b..0603e319 100644 --- a/src/kbmod/work_unit.py +++ b/src/kbmod/work_unit.py @@ -389,22 +389,23 @@ def to_fits(self, filename, overwrite=False): for i in range(self.im_stack.img_count()): layered = self.im_stack.get_single_image(i) + obstime = layered.get_obstime() c_indices = self._per_image_indices[i] n_indices = len(c_indices) img_wcs = self.get_wcs(i) - sci_hdu = raw_image_to_hdu(layered.get_science(), img_wcs) + sci_hdu = raw_image_to_hdu(layered.get_science(), obstime, img_wcs) sci_hdu.name = f"SCI_{i}" sci_hdu.header["NIND"] = n_indices for j in range(n_indices): sci_hdu.header[f"IND_{j}"] = c_indices[j] hdul.append(sci_hdu) - var_hdu = raw_image_to_hdu(layered.get_variance()) + var_hdu = raw_image_to_hdu(layered.get_variance(), obstime) var_hdu.name = f"VAR_{i}" hdul.append(var_hdu) - msk_hdu = raw_image_to_hdu(layered.get_mask()) + msk_hdu = raw_image_to_hdu(layered.get_mask(), obstime) msk_hdu.name = f"MSK_{i}" hdul.append(msk_hdu) @@ -463,23 +464,24 @@ def to_sharded_fits(self, filename, directory, overwrite=False): for i in range(self.im_stack.img_count()): layered = self.im_stack.get_single_image(i) + obstime = layered.get_obstime() c_indices = self._per_image_indices[i] n_indices = len(c_indices) sub_hdul = fits.HDUList() img_wcs = self.get_wcs(i) - sci_hdu = raw_image_to_hdu(layered.get_science(), img_wcs) + sci_hdu = raw_image_to_hdu(layered.get_science(), obstime, img_wcs) sci_hdu.name = f"SCI_{i}" sci_hdu.header["NIND"] = n_indices for j in range(n_indices): sci_hdu.header[f"IND_{j}"] = c_indices[j] sub_hdul.append(sci_hdu) - var_hdu = raw_image_to_hdu(layered.get_variance()) + var_hdu = raw_image_to_hdu(layered.get_variance(), obstime) var_hdu.name = f"VAR_{i}" sub_hdul.append(var_hdu) - msk_hdu = raw_image_to_hdu(layered.get_mask()) + msk_hdu = raw_image_to_hdu(layered.get_mask(), obstime) msk_hdu.name = f"MSK_{i}" sub_hdul.append(msk_hdu) @@ -835,13 +837,15 @@ def load_layered_image_from_shard(file_path): return img -def raw_image_to_hdu(img, wcs=None): +def raw_image_to_hdu(img, obstime, wcs=None): """Helper function that creates a HDU out of RawImage. Parameters ---------- img : `RawImage` The RawImage to convert. + obstime : `float` + The time of the observation. wcs : `astropy.wcs.WCS` An optional WCS to include in the header. @@ -857,28 +861,6 @@ def raw_image_to_hdu(img, wcs=None): append_wcs_to_hdu_header(wcs, hdu.header) # Set the time stamp. - hdu.header["MJD"] = img.obstime - return hdu - - -def hdu_to_raw_image(hdu): - """Helper function that creates a RawImage from a HDU. - - Parameters - ---------- - hdu : `astropy.io.fits.hdu.image.ImageHDU` - The image extension. + hdu.header["MJD"] = obstime - Returns - ------- - img : `RawImage` or None - The RawImage if there is valid data and None otherwise. - """ - img = None - if isinstance(hdu, fits.hdu.image.ImageHDU): - # This will be a copy whenever dtype != np.single including when - # endianness doesn't match the native float. - img = RawImage(hdu.data.astype(np.single)) - if "MJD" in hdu.header: - img.obstime = hdu.header["MJD"] - return img + return hdu diff --git a/tests/test_butlerstd.py b/tests/test_butlerstd.py index 8647ada7..ba810af0 100644 --- a/tests/test_butlerstd.py +++ b/tests/test_butlerstd.py @@ -223,9 +223,6 @@ def test_to_layered_image(self): # times can only be compred approximately, because sometimes we # calculate the time in the middle of the exposure self.assertAlmostEqual(expected_mjd, img.get_obstime(), 2) - self.assertAlmostEqual(expected_mjd, img.get_science().obstime, 2) - self.assertAlmostEqual(expected_mjd, img.get_variance().obstime, 2) - self.assertAlmostEqual(expected_mjd, img.get_mask().obstime, 2) if __name__ == "__main__": diff --git a/tests/test_fake_data_creator.py b/tests/test_fake_data_creator.py index ea92205b..8bc92362 100644 --- a/tests/test_fake_data_creator.py +++ b/tests/test_fake_data_creator.py @@ -25,7 +25,7 @@ def test_create_fake_times(self): self.assertAlmostEqual(times2[i], expected[i]) def test_add_fake_object(self): - img = RawImage(20, 10, 0.0, 1.0) # All zero image. + img = RawImage(20, 10, 0.0) # All zero image. p = PSF(np.full((3, 3), 1.0 / 9.0)) # Equal PSF. add_fake_object(img, 5.5, 3.5, 100.0, p) diff --git a/tests/test_psi_phi_array.py b/tests/test_psi_phi_array.py index 900ed638..050db096 100644 --- a/tests/test_psi_phi_array.py +++ b/tests/test_psi_phi_array.py @@ -30,14 +30,14 @@ def setUp(self): psi_1_vals = np.arange(0, self.width * self.height, dtype=np.single) psi_1_arr = psi_1_vals.reshape(self.height, self.width) - self.psi_1 = RawImage(img=psi_1_arr, obs_time=1.0) + self.psi_1 = RawImage(img=psi_1_arr) psi_2_vals = np.arange(self.width * self.height, 2 * self.width * self.height, dtype=np.single) psi_2_arr = psi_2_vals.reshape(self.height, self.width) - self.psi_2 = RawImage(img=psi_2_arr, obs_time=2.0) + self.psi_2 = RawImage(img=psi_2_arr) - self.phi_1 = RawImage(np.full((self.height, self.width), 0.1, dtype=np.single), obs_time=1.0) - self.phi_2 = RawImage(np.full((self.height, self.width), 0.2, dtype=np.single), obs_time=2.0) + self.phi_1 = RawImage(np.full((self.height, self.width), 0.1, dtype=np.single)) + self.phi_2 = RawImage(np.full((self.height, self.width), 0.2, dtype=np.single)) self.zeroed_times = [0.0, 1.0] diff --git a/tests/test_raw_image.py b/tests/test_raw_image.py index fa4ddadc..9bd86a3c 100644 --- a/tests/test_raw_image.py +++ b/tests/test_raw_image.py @@ -37,55 +37,47 @@ def test_create(self): img = RawImage() self.assertEqual(img.width, 0) self.assertEqual(img.height, 0) - self.assertEqual(img.obstime, -1.0) # from NumPy arrays - img = RawImage(img=self.array, obs_time=10.0) + img = RawImage(img=self.array) self.assertEqual(img.image.shape, (self.height, self.width)) - self.assertEqual(img.obstime, 10.0) self.assertEqual(img.npixels, self.width * self.height) self.assertTrue((img.image == self.array).all()) img2 = RawImage(img=self.array) self.assertTrue((img2.image == img.image).all()) - self.assertEqual(img2.obstime, -1.0) # from dimensions img = RawImage(self.width, self.height) self.assertEqual(img.image.shape, (self.height, self.width)) - self.assertEqual(img.obstime, -1.0) self.assertTrue((img.image == 0).all()) # dimensions and optional values img = RawImage(self.height, self.width, 10) self.assertTrue((img.image == 10).all()) - img = RawImage(self.height, self.width, 10, 12.5) + img = RawImage(self.height, self.width, 10) self.assertTrue((img.image == 10).all()) - self.assertEqual(img.obstime, 12.5) - img = RawImage(self.height, self.width, value=7.5, obs_time=12.5) + img = RawImage(self.height, self.width, value=7.5) self.assertTrue((img.image == 7.5).all()) - self.assertEqual(img.obstime, 12.5) # copy constructor, set the old image to all zeros and change the time. - img = RawImage(img=self.array, obs_time=10.0) + img = RawImage(img=self.array) img2 = RawImage(img) img.set_all(0.0) - img.obstime = 1.0 self.assertTrue((img2.image == self.array).all()) - self.assertEqual(img2.obstime, 10.0) def test_pixel_getters(self): """Test RawImage masked pixel value getters""" - img = RawImage(img=self.array, obs_time=10.0) + img = RawImage(img=self.array) self.assertFalse(pixel_value_valid(img.get_pixel(-1, 5))) self.assertFalse(pixel_value_valid(img.get_pixel(5, self.width))) self.assertFalse(pixel_value_valid(img.get_pixel(5, -1))) self.assertFalse(pixel_value_valid(img.get_pixel(self.height, 5))) def test_contains(self): - img = RawImage(img=self.array, obs_time=10.0) + img = RawImage(img=self.array) self.assertTrue(img.contains_index(0, 0)) self.assertTrue(img.contains_index(1, 2)) self.assertFalse(img.contains_index(1, -1)) @@ -102,7 +94,7 @@ def test_contains(self): self.assertFalse(img.contains_point(1.0, self.height + 1e-4)) def test_validity_checker(self): - img = RawImage(img=np.array([[0, 0], [0, 0]]).astype(np.float32), obs_time=10.0) + img = RawImage(img=np.array([[0, 0], [0, 0]]).astype(np.float32)) self.assertTrue(img.pixel_has_data(0, 0)) img.set_pixel(0, 0, np.nan) @@ -119,7 +111,7 @@ def test_validity_checker(self): def test_interpolated_add(self): """Test that we can add values to the pixel.""" - img = RawImage(img=self.array, obs_time=10.0) + img = RawImage(img=self.array) # Get the original value using (r, c) lookup. org_val17 = img.get_pixel(1, 7) @@ -153,7 +145,7 @@ def test_interpolated_add(self): def test_approx_equal(self): """Test RawImage pixel value setters.""" - img = RawImage(img=self.array, obs_time=10.0) + img = RawImage(img=self.array) # This test is testing L^\infy norm closeness. Eigen isApprox uses L2 # norm closeness. diff --git a/tests/test_standardizer.py b/tests/test_standardizer.py index 750224f5..cfd0d9a6 100644 --- a/tests/test_standardizer.py +++ b/tests/test_standardizer.py @@ -364,9 +364,6 @@ def test_to_layered_image(self): # Test that we correctly set metadata self.assertEqual(expected_mjd, img.get_obstime()) - self.assertEqual(expected_mjd, img.get_science().obstime) - self.assertEqual(expected_mjd, img.get_variance().obstime) - self.assertEqual(expected_mjd, img.get_mask().obstime) if __name__ == "__main__": diff --git a/tests/test_work_unit.py b/tests/test_work_unit.py index dbcf7b8e..01cb3022 100644 --- a/tests/test_work_unit.py +++ b/tests/test_work_unit.py @@ -15,7 +15,7 @@ import kbmod.search as kb from kbmod.reprojection_utils import fit_barycentric_wcs from kbmod.wcs_utils import make_fake_wcs, wcs_fits_equal -from kbmod.work_unit import hdu_to_raw_image, raw_image_to_hdu, WorkUnit +from kbmod.work_unit import raw_image_to_hdu, WorkUnit import numpy.testing as npt From 88eccdd4a9053efa66441d6c811adce903652107 Mon Sep 17 00:00:00 2001 From: Jeremy Kubica <104161096+jeremykubica@users.noreply.github.com> Date: Fri, 27 Sep 2024 14:44:42 -0400 Subject: [PATCH 2/3] Fix the logging message for coadd generation --- src/kbmod/search/stamp_creator.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/kbmod/search/stamp_creator.cpp b/src/kbmod/search/stamp_creator.cpp index c64ff1e2..08c10a2c 100644 --- a/src/kbmod/search/stamp_creator.cpp +++ b/src/kbmod/search/stamp_creator.cpp @@ -65,8 +65,8 @@ std::vector StampCreator::get_coadded_stamps(ImageStack& stack, std::v std::vector>& 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)) { From b70bf8c4c0ed4fa04b85d031518c48e615602f0b Mon Sep 17 00:00:00 2001 From: Jeremy Kubica <104161096+jeremykubica@users.noreply.github.com> Date: Sun, 29 Sep 2024 15:36:16 -0400 Subject: [PATCH 3/3] Fix segfault in weighted variance stamps --- src/kbmod/image_collection.py | 10 +++++++++- src/kbmod/search/stamp_creator.cpp | 16 ++++++++++------ .../fits_standardizers/fits_standardizer.py | 4 ++++ 3 files changed, 23 insertions(+), 7 deletions(-) diff --git a/src/kbmod/image_collection.py b/src/kbmod/image_collection.py index 8dc900a3..dce44354 100644 --- a/src/kbmod/image_collection.py +++ b/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`. diff --git a/src/kbmod/search/stamp_creator.cpp b/src/kbmod/search/stamp_creator.cpp index 08c10a2c..c5d1341a 100644 --- a/src/kbmod/search/stamp_creator.cpp +++ b/src/kbmod/search/stamp_creator.cpp @@ -14,7 +14,7 @@ std::vector StampCreator::create_stamps(ImageStack& stack, const Traje bool keep_no_data, const std::vector& 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 stamps; unsigned int num_times = stack.img_count(); @@ -277,7 +277,7 @@ std::vector StampCreator::create_variance_stamps(ImageStack& stack, co int radius, const std::vector& 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 stamps; unsigned int num_times = stack.img_count(); @@ -297,22 +297,26 @@ std::vector StampCreator::create_variance_stamps(ImageStack& stack, co RawImage StampCreator::get_variance_weighted_stamp(ImageStack& stack, const Trajectory& trj, int radius, const std::vector& 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 empty_vect; std::vector sci_stamps = create_stamps(stack, trj, radius, true /*=keep_no_data*/, use_index); std::vector 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)) { diff --git a/src/kbmod/standardizers/fits_standardizers/fits_standardizer.py b/src/kbmod/standardizers/fits_standardizers/fits_standardizer.py index dfce61ad..ae29b21d 100644 --- a/src/kbmod/standardizers/fits_standardizers/fits_standardizer.py +++ b/src/kbmod/standardizers/fits_standardizers/fits_standardizer.py @@ -410,6 +410,10 @@ 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))