From 4de6e337a61f4980d30aa0ccc4e45b30804c3294 Mon Sep 17 00:00:00 2001 From: Jeremy Kubica <104161096+jeremykubica@users.noreply.github.com> Date: Fri, 22 Sep 2023 10:17:03 -0400 Subject: [PATCH 1/5] Move fake object creation out of LayeredImage --- README.md | 10 ++++--- src/kbmod/fake_data_creator.py | 44 ++++++++++++++++++++++++++++--- src/kbmod/search/ImageStack.h | 1 + src/kbmod/search/LayeredImage.cpp | 34 +++++------------------- src/kbmod/search/LayeredImage.h | 5 ---- src/kbmod/search/bindings.cpp | 7 +++-- tests/regression_test.py | 3 ++- tests/test_analysis_utils.py | 11 +++++--- tests/test_bilinear_interp.py | 3 ++- tests/test_image_stack.py | 3 ++- tests/test_layered_image.py | 25 +++--------------- tests/test_readme_example.py | 8 +++++- tests/test_search.py | 5 +++- tests/test_search_encode.py | 5 +++- tests/test_search_filter.py | 5 +++- tests/test_stamp_parity.py | 5 +++- 16 files changed, 100 insertions(+), 74 deletions(-) diff --git a/README.md b/README.md index 33031a747..31983e377 100644 --- a/README.md +++ b/README.md @@ -111,9 +111,13 @@ velocity = (2, 0) # Inject object into images for im in imgs: dt = im.get_obstime() - t0 - im.add_object(position[0] + dt * velocity[0], - position[1] + dt * velocity[1], - flux) + add_fake_object( + im, + position[0] + dt * velocity[0], + position[1] + dt * velocity[1], + flux, + psf, + ) # Create a new image stack with the inserted object. stack = kb.image_stack(imgs) diff --git a/src/kbmod/fake_data_creator.py b/src/kbmod/fake_data_creator.py index 3f1e8b50f..0754ec3ae 100644 --- a/src/kbmod/fake_data_creator.py +++ b/src/kbmod/fake_data_creator.py @@ -15,6 +15,44 @@ from kbmod.search import * +def add_fake_object(img, x, y, flux, psf=None): + """Add a fake object to a LayeredImage or RawImage + + Parameters + ---------- + img : RawImage or LayeredImage + The image to modify. + x : float + The x pixel location of the fake object. + y : float + The y pixel location of the fake object. + flux : float + The flux value. + psf : PointSpreadFunc + The PSF for the image. + """ + if type(img) is LayeredImage: + sci = img.get_science() + else: + sci = img + + if psf is None: + sci.add_pixel_interp(x, y, flux) + else: + dim = psf.get_dim() + initial_x = x - psf.get_radius() + initial_y = y - psf.get_radius() + + for i in range(dim): + for j in range(dim): + sci.add_pixel_interp(initial_x + i, initial_y + j, flux * psf.get_value(i, j)) + + # The python/C++ interface requires us to explicitly re-set the science + # image in a LayeredImage. + if sci is not img: + img.set_science(sci) + + class FakeDataSet: """This class creates fake data sets for testing and demo notebooks.""" @@ -107,9 +145,9 @@ def insert_object(self, trj): # Get the image for the timestep, add the object, and # re-set the image. This last step needs to be done # explicitly because of how pybind handles references. - current_layered_image = self.stack.get_single_image(i) - current_layered_image.add_object(px, py, trj.flux) - self.stack.set_single_image(i, current_layered_image) + current = self.stack.get_single_image(i) + add_fake_object(current, px, py, trj.flux, current.get_psf()) + self.stack.set_single_image(i, current) # Save the trajectory into the internal list. self.trajectories.append(trj) diff --git a/src/kbmod/search/ImageStack.h b/src/kbmod/search/ImageStack.h index c9e013c9e..f7a84d667 100644 --- a/src/kbmod/search/ImageStack.h +++ b/src/kbmod/search/ImageStack.h @@ -23,6 +23,7 @@ class ImageStack { public: ImageStack(const std::vector& filenames, const std::vector& psfs); ImageStack(const std::vector& imgs); + ImageStack(const std::vector &science, const std::vector &variance); // Simple getters. unsigned imgCount() const { return images.size(); } diff --git a/src/kbmod/search/LayeredImage.cpp b/src/kbmod/search/LayeredImage.cpp index 6f6f453af..2d7d5019a 100644 --- a/src/kbmod/search/LayeredImage.cpp +++ b/src/kbmod/search/LayeredImage.cpp @@ -9,9 +9,7 @@ namespace search { -LayeredImage::LayeredImage(std::string path, const PointSpreadFunc& psf) : psf(psf), psf_sq(psf) { - psf_sq.squarePSF(); - +LayeredImage::LayeredImage(std::string path, const PointSpreadFunc& psf) : psf(psf) { int f_begin = path.find_last_of("/"); int f_end = path.find_last_of(".fits") - 4; filename = path.substr(f_begin, f_end - f_begin); @@ -35,7 +33,7 @@ LayeredImage::LayeredImage(std::string path, const PointSpreadFunc& psf) : psf(p LayeredImage::LayeredImage(const RawImage& sci, const RawImage& var, const RawImage& msk, const PointSpreadFunc& psf) - : psf(psf), psf_sq(psf) { + : psf(psf) { // Get the dimensions of the science layer and check for consistency with // the other two layers. width = sci.getWidth(); @@ -45,9 +43,6 @@ LayeredImage::LayeredImage(const RawImage& sci, const RawImage& var, const RawIm if (width != msk.getWidth() or height != msk.getHeight()) throw std::runtime_error("Science and Mask layers are not the same size."); - // Set the remaining variables. - psf_sq.squarePSF(); - // Copy the image layers. science = sci; mask = msk; @@ -60,11 +55,10 @@ LayeredImage::LayeredImage(std::string name, int w, int h, float noise_stdev, fl LayeredImage::LayeredImage(std::string name, int w, int h, float noise_stdev, float pixel_variance, double time, const PointSpreadFunc& psf, int seed) - : psf(psf), psf_sq(psf) { + : psf(psf) { filename = name; width = w; height = h; - psf_sq.squarePSF(); std::vector raw_sci(width * height); std::random_device r; @@ -84,24 +78,6 @@ LayeredImage::LayeredImage(std::string name, int w, int h, float noise_stdev, fl void LayeredImage::setPSF(const PointSpreadFunc& new_psf) { psf = new_psf; - psf_sq = new_psf; - psf_sq.squarePSF(); -} - -void LayeredImage::addObject(float x, float y, float flux) { - const std::vector& k = psf.getKernel(); - int dim = psf.getDim(); - float initial_x = x - static_cast(psf.getRadius()); - float initial_y = y - static_cast(psf.getRadius()); - - int count = 0; - for (int i = 0; i < dim; ++i) { - for (int j = 0; j < dim; ++j) { - science.addPixelInterp(initial_x + static_cast(i), initial_y + static_cast(j), - flux * k[count]); - count++; - } - } } void LayeredImage::growMask(int steps) { @@ -111,6 +87,10 @@ void LayeredImage::growMask(int steps) { void LayeredImage::convolvePSF() { science.convolve(psf); + + // Square the PSF use that on the variance image. + psf_sq = new_psf; + psf_sq.squarePSF(); variance.convolve(psf_sq); } diff --git a/src/kbmod/search/LayeredImage.h b/src/kbmod/search/LayeredImage.h index b258944e4..f1635559e 100644 --- a/src/kbmod/search/LayeredImage.h +++ b/src/kbmod/search/LayeredImage.h @@ -36,7 +36,6 @@ class LayeredImage { // Set an image specific point spread function. void setPSF(const PointSpreadFunc& psf); const PointSpreadFunc& getPSF() const { return psf; } - const PointSpreadFunc& getPSFSQ() const { return psf_sq; } // Basic getter functions for image data. std::string getName() const { return filename; } @@ -65,9 +64,6 @@ class LayeredImage { // Subtracts a template image from the science layer. void subtractTemplate(const RawImage& sub_template); - // Adds an (artificial) object to the image (science) data. - void addObject(float x, float y, float flux); - // Saves the data in each later to a file. void saveLayers(const std::string& path); @@ -91,7 +87,6 @@ class LayeredImage { unsigned height; PointSpreadFunc psf; - PointSpreadFunc psf_sq; RawImage science; RawImage mask; RawImage variance; diff --git a/src/kbmod/search/bindings.cpp b/src/kbmod/search/bindings.cpp index 4fe03ebe1..c3e8100df 100644 --- a/src/kbmod/search/bindings.cpp +++ b/src/kbmod/search/bindings.cpp @@ -100,6 +100,7 @@ PYBIND11_MODULE(search, m) { .def("create_stamp", &ri::createStamp) .def("set_pixel", &ri::setPixel, "Set the value of a given pixel.") .def("add_pixel", &ri::addToPixel, "Add to the value of a given pixel.") + .def("add_pixel_interp", &ri::addPixelInterp, "Add to the interpolated value of a given pixel." .def("apply_mask", &ri::applyMask) .def("grow_mask", &ri::growMask) .def("pixel_has_data", &ri::pixelHasData, @@ -137,9 +138,6 @@ PYBIND11_MODULE(search, m) { )pbdoc") .def(py::init()) .def(py::init()) - .def("set_psf", &li::setPSF, "Sets the PSF object.") - .def("get_psf", &li::getPSF, "Returns the PSF object.") - .def("get_psfsq", &li::getPSFSQ) .def("apply_mask_flags", &li::applyMaskFlags) .def("apply_mask_threshold", &li::applyMaskThreshold) .def("sub_template", &li::subtractTemplate) @@ -151,7 +149,6 @@ PYBIND11_MODULE(search, m) { .def("set_mask", &li::setMask) .def("set_variance", &li::setVariance) .def("convolve_psf", &li::convolvePSF) - .def("add_object", &li::addObject) .def("grow_mask", &li::growMask) .def("get_name", &li::getName, "Returns the name of the layered image.") .def("get_width", &li::getWidth, "Returns the image's width in pixels.") @@ -159,6 +156,8 @@ PYBIND11_MODULE(search, m) { .def("get_npixels", &li::getNPixels, "Returns the image's total number of pixels.") .def("get_obstime", &li::getObstime, "Get the image's observation time.") .def("set_obstime", &li::setObstime, "Set the image's observation time.") + .def("get_psf", &li::getPSF, "Get the image's PSF.") + .def("set_psf", &li::getPSF, "Set the image's PSF.") .def("generate_psi_image", &li::generatePsiImage) .def("generate_phi_image", &li::generatePhiImage); py::class_(m, "image_stack") diff --git a/tests/regression_test.py b/tests/regression_test.py index 39a13c14f..1bc8c6479 100644 --- a/tests/regression_test.py +++ b/tests/regression_test.py @@ -12,6 +12,7 @@ import numpy as np from astropy.io import fits +from kbmod.fake_data_creator import add_fake_object from kbmod.file_utils import * from kbmod.run_search import run_search from kbmod.search import * @@ -213,7 +214,7 @@ def make_fake_image_stack(times, trjs, psf_vals): for trj in trjs: px = trj.x + time * trj.x_v + 0.5 py = trj.y + time * trj.y_v + 0.5 - img.add_object(px, py, trj.flux) + add_fake_object(img, px, py, trj.flux, p) imlist.append(img) stack = image_stack(imlist) diff --git a/tests/test_analysis_utils.py b/tests/test_analysis_utils.py index fe1d8a374..4f4bfade4 100644 --- a/tests/test_analysis_utils.py +++ b/tests/test_analysis_utils.py @@ -1,6 +1,7 @@ import unittest from kbmod.analysis_utils import * +from kbmod.fake_data_creator import add_fake_object from kbmod.result_list import * from kbmod.search import * @@ -179,10 +180,12 @@ def test_apply_stamp_filter(self): for i in range(self.img_count): time = i / self.img_count - self.imlist[i].add_object( + add_fake_object( + self.imlist[i], self.start_x + time * self.x_vel + 0.5, self.start_y + time * self.y_vel + 0.5, self.object_flux, + self.p, ) stack = image_stack(self.imlist) @@ -228,10 +231,12 @@ def test_apply_stamp_filter_2(self): for i in range(self.img_count): time = i / self.img_count - self.imlist[i].add_object( + add_fake_object( + self.imlist[i], self.start_x + time * self.x_vel, self.start_y + time * self.y_vel, self.object_flux, + self.p, ) stack = image_stack(self.imlist) @@ -346,7 +351,7 @@ def test_load_and_filter_results_lh(self): # Add the objects. for j, trj in enumerate(trjs): - im.add_object(trj.x, trj.y, fluxes[j]) + add_fake_object(im, trj.x, trj.y, fluxes[j], self.p) # Append the image. imlist.append(im) diff --git a/tests/test_bilinear_interp.py b/tests/test_bilinear_interp.py index ed0b071b9..c4dd8dae0 100644 --- a/tests/test_bilinear_interp.py +++ b/tests/test_bilinear_interp.py @@ -2,6 +2,7 @@ import numpy +from kbmod.fake_data_creator import add_fake_object import kbmod.search as kb @@ -12,7 +13,7 @@ def setUp(self): self.images = [] for c in range(self.im_count): im = kb.layered_image(str(c), 10, 10, 0.0, 1.0, c, p) - im.add_object(2 + c * 0.5 + 0.5, 2 + c * 0.5 + 0.5, 1) + add_fake_object(im, 2 + c * 0.5 + 0.5, 2 + c * 0.5 + 0.5, 1, p) self.images.append(im) def test_pixels(self): diff --git a/tests/test_image_stack.py b/tests/test_image_stack.py index e2d707b9e..9965ee370 100644 --- a/tests/test_image_stack.py +++ b/tests/test_image_stack.py @@ -1,6 +1,7 @@ import tempfile import unittest +from kbmod.fake_data_creator import add_fake_object from kbmod.search import * @@ -150,7 +151,7 @@ def test_different_psfs(self): last_val = -100.0 for i in range(self.num_images): img = self.im_stack.get_single_image(i) - img.add_object(10, 20, 500.0) + add_fake_object(im, 10, 20, 500.0, self.p[i]) sci = img.get_science() pix_val = sci.get_pixel(10, 20) diff --git a/tests/test_layered_image.py b/tests/test_layered_image.py index 539c8a895..e48fc8fdb 100644 --- a/tests/test_layered_image.py +++ b/tests/test_layered_image.py @@ -3,6 +3,7 @@ from astropy.io import fits +from kbmod.fake_data_creator import add_fake_object from kbmod.search import * @@ -72,28 +73,15 @@ def test_create_from_layers(self): self.assertEqual(variance.get_pixel(x, y), 1.0) self.assertAlmostEqual(science.get_pixel(x, y), x + 30.0 * y) - def test_add_object(self): - science = self.image.get_science() - science_50_50 = science.get_pixel(50, 50) - self.image.add_object(50, 50, 500.0) - - science = self.image.get_science() - self.assertLess(science_50_50, science.get_pixel(50, 50)) - def test_overwrite_psf(self): p1 = self.image.get_psf() self.assertEqual(p1.get_size(), 25) self.assertEqual(p1.get_dim(), 5) self.assertEqual(p1.get_radius(), 2) - psq1 = self.image.get_psfsq() - self.assertEqual(psq1.get_size(), 25) - self.assertEqual(psq1.get_dim(), 5) - self.assertEqual(psq1.get_radius(), 2) - # Get the science pixel with the original PSF blurring. science_org = self.image.get_science() - self.image.add_object(50, 50, 500.0) + add_fake_object(self.image, 50, 50, 500.0, p1) science_pixel_psf1 = self.image.get_science().get_pixel(50, 50) # Change the PSF. @@ -105,15 +93,10 @@ def test_overwrite_psf(self): self.assertEqual(p2.get_dim(), 1) self.assertEqual(p2.get_radius(), 0) - psq2 = self.image.get_psfsq() - self.assertEqual(psq2.get_size(), 1) - self.assertEqual(psq2.get_dim(), 1) - self.assertEqual(psq2.get_radius(), 0) - # Check that the science pixel with the new PSF blurring is # larger (because the PSF is tighter). self.image.set_science(science_org) - self.image.add_object(50, 50, 500.0) + add_fake_object(self.image, 50, 50, 500.0, p2) science_pixel_psf2 = self.image.get_science().get_pixel(50, 50) self.assertLess(science_pixel_psf1, science_pixel_psf2) @@ -122,7 +105,7 @@ def test_mask_threshold(self): threshold = 20.0 # Add an object brighter than the threshold. - self.image.add_object(50, 50, 500.0) + add_fake_object(self.image, 50, 50, 500.0, self.p) # Find all the pixels that should be masked. science = self.image.get_science() diff --git a/tests/test_readme_example.py b/tests/test_readme_example.py index c4da5849f..dc428545e 100644 --- a/tests/test_readme_example.py +++ b/tests/test_readme_example.py @@ -26,7 +26,13 @@ def test_make_and_copy(self): # Inject object into images for im in imgs: dt = im.get_obstime() - t0 - im.add_object(position[0] + dt * velocity[0], position[1] + dt * velocity[1], flux) + add_fake_object( + im, + position[0] + dt * velocity[0], + position[1] + dt * velocity[1], + flux, + psf, + ) # Create a new image stack with the inserted object. stack = kb.image_stack(imgs) diff --git a/tests/test_search.py b/tests/test_search.py index 2ad095db7..06cffb3f9 100644 --- a/tests/test_search.py +++ b/tests/test_search.py @@ -2,6 +2,7 @@ import numpy as np +from kbmod.fake_data_creator import add_fake_object from kbmod.search import * @@ -56,10 +57,12 @@ def setUp(self): im = layered_image( str(i), self.dim_x, self.dim_y, self.noise_level, self.variance, time, self.p, i ) - im.add_object( + add_fake_object( + im, self.start_x + time * self.x_vel + 0.5, self.start_y + time * self.y_vel + 0.5, self.object_flux, + self.p, ) # Mask a pixel in half the images. diff --git a/tests/test_search_encode.py b/tests/test_search_encode.py index 51d48a267..4e1f0941c 100644 --- a/tests/test_search_encode.py +++ b/tests/test_search_encode.py @@ -2,6 +2,7 @@ import numpy as np +from kbmod.fake_data_creator import add_fake_object from kbmod.search import * @@ -54,10 +55,12 @@ def setUp(self): im = layered_image( str(i), self.dim_x, self.dim_y, self.noise_level, self.variance, time, self.p, i ) - im.add_object( + add_fake_object( + im, self.start_x + time * self.x_vel + 0.5, self.start_y + time * self.y_vel + 0.5, self.object_flux, + self.p, ) self.imlist.append(im) self.stack = image_stack(self.imlist) diff --git a/tests/test_search_filter.py b/tests/test_search_filter.py index 62ea3f686..1c8eb9aca 100644 --- a/tests/test_search_filter.py +++ b/tests/test_search_filter.py @@ -2,6 +2,7 @@ import numpy as np +from kbmod.fake_data_creator import add_fake_object from kbmod.search import * @@ -54,10 +55,12 @@ def setUp(self): im = layered_image( str(i), self.dim_x, self.dim_y, self.noise_level, self.variance, time, self.p, i ) - im.add_object( + add_fake_object( + im, self.start_x + time * self.x_vel + 0.5, self.start_y + time * self.y_vel + 0.5, self.object_flux, + self.p, ) self.imlist.append(im) self.stack = image_stack(self.imlist) diff --git a/tests/test_stamp_parity.py b/tests/test_stamp_parity.py index c65915d60..a3d245880 100644 --- a/tests/test_stamp_parity.py +++ b/tests/test_stamp_parity.py @@ -8,6 +8,7 @@ import numpy as np +from kbmod.fake_data_creator import add_fake_object from kbmod.search import * @@ -59,10 +60,12 @@ def setUp(self): im = layered_image( str(i), self.dim_x, self.dim_y, self.noise_level, self.variance, time, self.p, i ) - im.add_object( + add_fake_object( + im, self.start_x + time * self.x_vel + 0.5, self.start_y + time * self.y_vel + 0.5, self.object_flux, + self.p, ) # Mask a pixel in half the images. From 2b34f6e26b57557d1777bf0b8d0b50cce3ec3ddd Mon Sep 17 00:00:00 2001 From: Jeremy Kubica <104161096+jeremykubica@users.noreply.github.com> Date: Fri, 22 Sep 2023 10:36:25 -0400 Subject: [PATCH 2/5] bug fixes --- src/kbmod/fake_data_creator.py | 2 +- src/kbmod/search/ImageStack.h | 1 - src/kbmod/search/LayeredImage.cpp | 6 ++++-- src/kbmod/search/bindings.cpp | 6 +++--- tests/test_image_stack.py | 2 +- tests/test_search.py | 8 +++++--- 6 files changed, 14 insertions(+), 11 deletions(-) diff --git a/src/kbmod/fake_data_creator.py b/src/kbmod/fake_data_creator.py index 0754ec3ae..3d50c7902 100644 --- a/src/kbmod/fake_data_creator.py +++ b/src/kbmod/fake_data_creator.py @@ -31,7 +31,7 @@ def add_fake_object(img, x, y, flux, psf=None): psf : PointSpreadFunc The PSF for the image. """ - if type(img) is LayeredImage: + if type(img) is layered_image: sci = img.get_science() else: sci = img diff --git a/src/kbmod/search/ImageStack.h b/src/kbmod/search/ImageStack.h index f7a84d667..c9e013c9e 100644 --- a/src/kbmod/search/ImageStack.h +++ b/src/kbmod/search/ImageStack.h @@ -23,7 +23,6 @@ class ImageStack { public: ImageStack(const std::vector& filenames, const std::vector& psfs); ImageStack(const std::vector& imgs); - ImageStack(const std::vector &science, const std::vector &variance); // Simple getters. unsigned imgCount() const { return images.size(); } diff --git a/src/kbmod/search/LayeredImage.cpp b/src/kbmod/search/LayeredImage.cpp index 2d7d5019a..8a678d881 100644 --- a/src/kbmod/search/LayeredImage.cpp +++ b/src/kbmod/search/LayeredImage.cpp @@ -89,7 +89,7 @@ void LayeredImage::convolvePSF() { science.convolve(psf); // Square the PSF use that on the variance image. - psf_sq = new_psf; + PointSpreadFunc psf_sq = PointSpreadFunc(psf); // Copy psf_sq.squarePSF(); variance.convolve(psf_sq); } @@ -219,7 +219,9 @@ RawImage LayeredImage::generatePhiImage() { } // Convolve with the PSF squared. - result.convolve(getPSFSQ()); + PointSpreadFunc psf_sq = PointSpreadFunc(psf); // Copy + psf_sq.squarePSF(); + result.convolve(psf_sq); return result; } diff --git a/src/kbmod/search/bindings.cpp b/src/kbmod/search/bindings.cpp index c3e8100df..385d93ff5 100644 --- a/src/kbmod/search/bindings.cpp +++ b/src/kbmod/search/bindings.cpp @@ -100,7 +100,7 @@ PYBIND11_MODULE(search, m) { .def("create_stamp", &ri::createStamp) .def("set_pixel", &ri::setPixel, "Set the value of a given pixel.") .def("add_pixel", &ri::addToPixel, "Add to the value of a given pixel.") - .def("add_pixel_interp", &ri::addPixelInterp, "Add to the interpolated value of a given pixel." + .def("add_pixel_interp", &ri::addPixelInterp, "Add to the interpolated value of a given pixel.") .def("apply_mask", &ri::applyMask) .def("grow_mask", &ri::growMask) .def("pixel_has_data", &ri::pixelHasData, @@ -138,6 +138,8 @@ PYBIND11_MODULE(search, m) { )pbdoc") .def(py::init()) .def(py::init()) + .def("set_psf", &li::setPSF, "Sets the PSF object.") + .def("get_psf", &li::getPSF, "Returns the PSF object.") .def("apply_mask_flags", &li::applyMaskFlags) .def("apply_mask_threshold", &li::applyMaskThreshold) .def("sub_template", &li::subtractTemplate) @@ -156,8 +158,6 @@ PYBIND11_MODULE(search, m) { .def("get_npixels", &li::getNPixels, "Returns the image's total number of pixels.") .def("get_obstime", &li::getObstime, "Get the image's observation time.") .def("set_obstime", &li::setObstime, "Set the image's observation time.") - .def("get_psf", &li::getPSF, "Get the image's PSF.") - .def("set_psf", &li::getPSF, "Set the image's PSF.") .def("generate_psi_image", &li::generatePsiImage) .def("generate_phi_image", &li::generatePhiImage); py::class_(m, "image_stack") diff --git a/tests/test_image_stack.py b/tests/test_image_stack.py index 9965ee370..2183c8dc9 100644 --- a/tests/test_image_stack.py +++ b/tests/test_image_stack.py @@ -151,7 +151,7 @@ def test_different_psfs(self): last_val = -100.0 for i in range(self.num_images): img = self.im_stack.get_single_image(i) - add_fake_object(im, 10, 20, 500.0, self.p[i]) + add_fake_object(img, 10, 20, 500.0, self.p[i]) sci = img.get_science() pix_val = sci.get_pixel(10, 20) diff --git a/tests/test_search.py b/tests/test_search.py index 06cffb3f9..411d86a65 100644 --- a/tests/test_search.py +++ b/tests/test_search.py @@ -95,11 +95,11 @@ def test_psiphi(self): # Image1 has a single object. image1 = layered_image("test1", 5, 10, 2.0, 4.0, 1.0, p) - image1.add_object(3.5, 2.5, 400.0) + add_fake_object(image1, 3.5, 2.5, 400.0, p) # Image2 has a single object and a masked pixel. image2 = layered_image("test2", 5, 10, 2.0, 4.0, 2.0, p) - image2.add_object(2.5, 4.5, 400.0) + add_fake_object(image2, 2.5, 4.5, 400.0, p) mask = image2.get_mask() mask.set_pixel(4, 9, 1) image2.set_mask(mask) @@ -215,10 +215,12 @@ def test_results_off_chip(self): im = layered_image( str(i), self.dim_x, self.dim_y, self.noise_level, self.variance, time, self.p, i ) - im.add_object( + add_fake_object( + im, trj.x + time * trj.x_v + 0.5, trj.y + time * trj.y_v + 0.5, self.object_flux, + self.p, ) imlist.append(im) stack = image_stack(imlist) From b01f0857cb8079230528a90482caa346ea709bed Mon Sep 17 00:00:00 2001 From: Jeremy Kubica <104161096+jeremykubica@users.noreply.github.com> Date: Fri, 22 Sep 2023 11:09:30 -0400 Subject: [PATCH 3/5] Add a no-op PSF --- src/kbmod/search/PointSpreadFunc.cpp | 7 +++++++ src/kbmod/search/PointSpreadFunc.h | 1 + src/kbmod/search/bindings.cpp | 1 + tests/test_layered_image.py | 4 ++-- tests/test_psf.py | 10 ++++++++++ 5 files changed, 21 insertions(+), 2 deletions(-) diff --git a/src/kbmod/search/PointSpreadFunc.cpp b/src/kbmod/search/PointSpreadFunc.cpp index 3038460a1..dd786aa94 100644 --- a/src/kbmod/search/PointSpreadFunc.cpp +++ b/src/kbmod/search/PointSpreadFunc.cpp @@ -9,6 +9,13 @@ namespace search { +PointSpreadFunc::PointSpreadFunc() : kernel(1, 1.0) { + dim = 1; + radius = 0; + width = 1e-12; + sum = 1.0; +} + PointSpreadFunc::PointSpreadFunc(float stdev) { width = stdev; float simple_gauss[MAX_KERNEL_RADIUS]; diff --git a/src/kbmod/search/PointSpreadFunc.h b/src/kbmod/search/PointSpreadFunc.h index c60498431..d270ae6ca 100644 --- a/src/kbmod/search/PointSpreadFunc.h +++ b/src/kbmod/search/PointSpreadFunc.h @@ -26,6 +26,7 @@ namespace search { class PointSpreadFunc { public: + PointSpreadFunc(); // Create a no-op PSF. PointSpreadFunc(float stdev); PointSpreadFunc(const PointSpreadFunc& other); // Copy constructor PointSpreadFunc(PointSpreadFunc&& other); // Move constructor diff --git a/src/kbmod/search/bindings.cpp b/src/kbmod/search/bindings.cpp index 385d93ff5..3cdc7985c 100644 --- a/src/kbmod/search/bindings.cpp +++ b/src/kbmod/search/bindings.cpp @@ -55,6 +55,7 @@ PYBIND11_MODULE(search, m) { 2, {m.getDim(), m.getDim()}, {sizeof(float) * m.getDim(), sizeof(float)}); }) + .def(py::init<>()) .def(py::init()) .def(py::init>()) .def(py::init()) diff --git a/tests/test_layered_image.py b/tests/test_layered_image.py index e48fc8fdb..2c4dad479 100644 --- a/tests/test_layered_image.py +++ b/tests/test_layered_image.py @@ -84,8 +84,8 @@ def test_overwrite_psf(self): add_fake_object(self.image, 50, 50, 500.0, p1) science_pixel_psf1 = self.image.get_science().get_pixel(50, 50) - # Change the PSF. - self.image.set_psf(psf(0.0001)) + # Change the PSF to a no-op. + self.image.set_psf(psf()) # Check that we retrieve the correct PSF. p2 = self.image.get_psf() diff --git a/tests/test_psf.py b/tests/test_psf.py index 0cae25267..01e2aee0b 100644 --- a/tests/test_psf.py +++ b/tests/test_psf.py @@ -9,6 +9,16 @@ def setUp(self): sigma_list = range(self.psf_count) self.psf_list = [psf(x / 5 + 0.2) for x in sigma_list] + def test_make_noop(self): + psf0 = psf() + self.assertEqual(psf1.get_size(), 1) + self.assertEqual(psf1.get_dim(), 1) + self.assertEqual(psf1.get_radius(), 0) + + kernel0 = psf0.get_kernel() + self.assertEqual(len(kernel0), 1) + self.assertEqual(kernel0[0], 1.0) + def test_make_and_copy(self): psf1 = psf(1.0) self.assertEqual(psf1.get_size(), 25) From b62c5ad8cd0e0ef19349a8839998fe76f20e656a Mon Sep 17 00:00:00 2001 From: Jeremy Kubica <104161096+jeremykubica@users.noreply.github.com> Date: Fri, 22 Sep 2023 11:31:44 -0400 Subject: [PATCH 4/5] Add ability to convolve with a given PSF --- src/kbmod/search/LayeredImage.cpp | 22 +++++++++++++--------- src/kbmod/search/LayeredImage.h | 3 +++ src/kbmod/search/bindings.cpp | 3 ++- tests/test_layered_image.py | 26 ++++++++++++++++++++++++++ 4 files changed, 44 insertions(+), 10 deletions(-) diff --git a/src/kbmod/search/LayeredImage.cpp b/src/kbmod/search/LayeredImage.cpp index 8a678d881..d1cfba0ce 100644 --- a/src/kbmod/search/LayeredImage.cpp +++ b/src/kbmod/search/LayeredImage.cpp @@ -85,13 +85,17 @@ void LayeredImage::growMask(int steps) { variance.growMask(steps); } -void LayeredImage::convolvePSF() { - science.convolve(psf); +void LayeredImage::convolveGivenPSF(const PointSpreadFunc& given_psf) { + science.convolve(given_psf); // Square the PSF use that on the variance image. - PointSpreadFunc psf_sq = PointSpreadFunc(psf); // Copy - psf_sq.squarePSF(); - variance.convolve(psf_sq); + PointSpreadFunc psfsq = PointSpreadFunc(given_psf); // Copy + psfsq.squarePSF(); + variance.convolve(psfsq); +} + +void LayeredImage::convolvePSF() { + convolveGivenPSF(psf); } void LayeredImage::applyMaskFlags(int flags, const std::vector& exceptions) { @@ -197,7 +201,7 @@ RawImage LayeredImage::generatePsiImage() { } // Convolve with the PSF. - result.convolve(getPSF()); + result.convolve(psf); return result; } @@ -219,9 +223,9 @@ RawImage LayeredImage::generatePhiImage() { } // Convolve with the PSF squared. - PointSpreadFunc psf_sq = PointSpreadFunc(psf); // Copy - psf_sq.squarePSF(); - result.convolve(psf_sq); + PointSpreadFunc psfsq = PointSpreadFunc(psf); // Copy + psfsq.squarePSF(); + result.convolve(psfsq); return result; } diff --git a/src/kbmod/search/LayeredImage.h b/src/kbmod/search/LayeredImage.h index f1635559e..4a4cb2781 100644 --- a/src/kbmod/search/LayeredImage.h +++ b/src/kbmod/search/LayeredImage.h @@ -72,7 +72,10 @@ class LayeredImage { void setMask(RawImage& im); void setVariance(RawImage& im); + // Convolve with a given PSF or the default one. void convolvePSF(); + void convolveGivenPSF(const PointSpreadFunc& psf); + virtual ~LayeredImage(){}; // Generate psi and phi images from the science and variance layers. diff --git a/src/kbmod/search/bindings.cpp b/src/kbmod/search/bindings.cpp index 3cdc7985c..ebf6c49b1 100644 --- a/src/kbmod/search/bindings.cpp +++ b/src/kbmod/search/bindings.cpp @@ -151,7 +151,8 @@ PYBIND11_MODULE(search, m) { .def("set_science", &li::setScience) .def("set_mask", &li::setMask) .def("set_variance", &li::setVariance) - .def("convolve_psf", &li::convolvePSF) + .def("convolve_psf", &li::convolvePSF, "Convolve each layer with the object's PSF.") + .def("convolve_given_psf", &li::convolveGivenPSF, "Convolve each layer with a given PSF.") .def("grow_mask", &li::growMask) .def("get_name", &li::getName, "Returns the name of the layered image.") .def("get_width", &li::getWidth, "Returns the image's width in pixels.") diff --git a/tests/test_layered_image.py b/tests/test_layered_image.py index 2c4dad479..5d7cf43cd 100644 --- a/tests/test_layered_image.py +++ b/tests/test_layered_image.py @@ -73,6 +73,32 @@ def test_create_from_layers(self): self.assertEqual(variance.get_pixel(x, y), 1.0) self.assertAlmostEqual(science.get_pixel(x, y), x + 30.0 * y) + def test_convolve_psf(self): + sci0 = self.image.get_science() + var0 = self.image.get_variance() + msk0 = self.image.get_mask() + + # Create a copy of the image. + img_b = layered_image(sci, var, mask, self.p) + + # A no-op PSF does not change the image. + img_b.convolve_given_psf(psf()) + sci1 = img_b.get_science() + var1 = img_b.get_variance() + for y in range(img2.get_height()): + for x in range(img2.get_width()): + self.assertAlmostEqual(sci0.get_pixel(x, y), sci1.get_pixel(x, y)) + self.assertAlmostEqual(var0.get_pixel(x, y), var1.get_pixel(x, y)) + + # The default PSF (stdev=1.0) DOES have the image. + img_b.convolve_psf(psf()) + sci1 = img_b.get_science() + var1 = img_b.get_variance() + for y in range(img2.get_height()): + for x in range(img2.get_width()): + self.assertNotAlmostEqual(sci0.get_pixel(x, y), sci1.get_pixel(x, y)) + self.assertNotAlmostEqual(var0.get_pixel(x, y), var1.get_pixel(x, y)) + def test_overwrite_psf(self): p1 = self.image.get_psf() self.assertEqual(p1.get_size(), 25) From 1401864a7ffd0dff7083e6005aed161641e8ef92 Mon Sep 17 00:00:00 2001 From: Jeremy Kubica <104161096+jeremykubica@users.noreply.github.com> Date: Fri, 22 Sep 2023 11:39:48 -0400 Subject: [PATCH 5/5] Bug fixes --- tests/test_layered_image.py | 12 ++++++------ tests/test_psf.py | 6 +++--- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/tests/test_layered_image.py b/tests/test_layered_image.py index 5d7cf43cd..334233e4f 100644 --- a/tests/test_layered_image.py +++ b/tests/test_layered_image.py @@ -79,23 +79,23 @@ def test_convolve_psf(self): msk0 = self.image.get_mask() # Create a copy of the image. - img_b = layered_image(sci, var, mask, self.p) + img_b = layered_image(sci0, var0, msk0, self.p) # A no-op PSF does not change the image. img_b.convolve_given_psf(psf()) sci1 = img_b.get_science() var1 = img_b.get_variance() - for y in range(img2.get_height()): - for x in range(img2.get_width()): + for y in range(img_b.get_height()): + for x in range(img_b.get_width()): self.assertAlmostEqual(sci0.get_pixel(x, y), sci1.get_pixel(x, y)) self.assertAlmostEqual(var0.get_pixel(x, y), var1.get_pixel(x, y)) # The default PSF (stdev=1.0) DOES have the image. - img_b.convolve_psf(psf()) + img_b.convolve_psf() sci1 = img_b.get_science() var1 = img_b.get_variance() - for y in range(img2.get_height()): - for x in range(img2.get_width()): + for y in range(img_b.get_height()): + for x in range(img_b.get_width()): self.assertNotAlmostEqual(sci0.get_pixel(x, y), sci1.get_pixel(x, y)) self.assertNotAlmostEqual(var0.get_pixel(x, y), var1.get_pixel(x, y)) diff --git a/tests/test_psf.py b/tests/test_psf.py index 01e2aee0b..4d088c90d 100644 --- a/tests/test_psf.py +++ b/tests/test_psf.py @@ -11,9 +11,9 @@ def setUp(self): def test_make_noop(self): psf0 = psf() - self.assertEqual(psf1.get_size(), 1) - self.assertEqual(psf1.get_dim(), 1) - self.assertEqual(psf1.get_radius(), 0) + self.assertEqual(psf0.get_size(), 1) + self.assertEqual(psf0.get_dim(), 1) + self.assertEqual(psf0.get_radius(), 0) kernel0 = psf0.get_kernel() self.assertEqual(len(kernel0), 1)