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] 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 86e59b453..b8297c0a6 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 50baaeda2..d07f18008 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 d30127c86..292e31358 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 6691d49d8..70400bdc3 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 7a108a95e..94bae3386 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 6b1df1133..c27616cbf 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 42c49bd35..c135761c1 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 dfce61adc..28efbe11d 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 bf9e932bb..0603e3197 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 8647ada79..ba810af09 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 ea92205b8..8bc923623 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 900ed6384..050db096c 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 fa4ddadc2..9bd86a3c7 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 750224f52..cfd0d9a69 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 dbcf7b8ed..01cb30224 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