diff --git a/src/kbmod/search/bindings.cpp b/src/kbmod/search/bindings.cpp index f0d2d1b04..4175c9639 100644 --- a/src/kbmod/search/bindings.cpp +++ b/src/kbmod/search/bindings.cpp @@ -18,7 +18,6 @@ namespace py = pybind11; #include "kernel_testing_helpers.cpp" #include "psi_phi_array.cpp" - PYBIND11_MODULE(search, m) { m.attr("KB_NO_DATA") = pybind11::float_(search::NO_DATA); m.attr("HAS_GPU") = pybind11::bool_(search::HAVE_GPU); diff --git a/src/kbmod/search/layered_image.cpp b/src/kbmod/search/layered_image.cpp index efe335fdd..24f961e9d 100644 --- a/src/kbmod/search/layered_image.cpp +++ b/src/kbmod/search/layered_image.cpp @@ -229,7 +229,7 @@ RawImage LayeredImage::generate_psi_image() { const int num_pixels = get_npixels(); for (int p = 0; p < num_pixels; ++p) { float var_pix = var_array[p]; - if (var_pix != NO_DATA) { + if (var_pix != NO_DATA && var_pix != 0.0 && sci_array[p] != NO_DATA) { result_arr[p] = sci_array[p] / var_pix; } else { result_arr[p] = NO_DATA; @@ -251,7 +251,7 @@ RawImage LayeredImage::generate_phi_image() { const int num_pixels = get_npixels(); for (int p = 0; p < num_pixels; ++p) { float var_pix = var_array[p]; - if (var_pix != NO_DATA) { + if (var_pix != NO_DATA && var_pix != 0.0) { result_arr[p] = 1.0 / var_pix; } else { result_arr[p] = NO_DATA; @@ -280,28 +280,26 @@ static void layered_image_bindings(py::module& m) { .def("contains", &li::contains, pydocs::DOC_LayeredImage_cointains) .def("get_science_pixel", &li::get_science_pixel, pydocs::DOC_LayeredImage_get_science_pixel) .def("get_variance_pixel", &li::get_variance_pixel, pydocs::DOC_LayeredImage_get_variance_pixel) - .def( - "contains", - [](li& cls, int i, int j) { - return cls.contains({i, j}); - }) - .def( - "get_science_pixel", - [](li& cls, int i, int j) { - return cls.get_science_pixel({i, j}); - }) - .def( - "get_variance_pixel", - [](li& cls, int i, int j) { - return cls.get_variance_pixel({i, j}); - }) + .def("contains", + [](li& cls, int i, int j) { + return cls.contains({i, j}); + }) + .def("get_science_pixel", + [](li& cls, int i, int j) { + return cls.get_science_pixel({i, j}); + }) + .def("get_variance_pixel", + [](li& cls, int i, int j) { + return cls.get_variance_pixel({i, j}); + }) .def("set_psf", &li::set_psf, pydocs::DOC_LayeredImage_set_psf) .def("get_psf", &li::get_psf, py::return_value_policy::reference_internal, pydocs::DOC_LayeredImage_get_psf) .def("binarize_mask", &li::binarize_mask, pydocs::DOC_LayeredImage_binarize_mask) .def("apply_mask", &li::apply_mask, pydocs::DOC_LayeredImage_apply_mask) .def("union_masks", &li::union_masks, pydocs::DOC_LayeredImage_union_masks) - .def("union_threshold_masking", &li::union_threshold_masking, pydocs::DOC_LayeredImage_union_threshold_masking) + .def("union_threshold_masking", &li::union_threshold_masking, + pydocs::DOC_LayeredImage_union_threshold_masking) .def("sub_template", &li::subtract_template, pydocs::DOC_LayeredImage_sub_template) .def("save_layers", &li::save_layers, pydocs::DOC_LayeredImage_save_layers) .def("get_science", &li::get_science, py::return_value_policy::reference_internal, diff --git a/src/kbmod/search/layered_image.h b/src/kbmod/search/layered_image.h index 11fe159ee..b2e11e421 100644 --- a/src/kbmod/search/layered_image.h +++ b/src/kbmod/search/layered_image.h @@ -42,13 +42,13 @@ class LayeredImage { // Getter functions for the pixels of the science and variance layers that check // the mask layer for any set bits. inline float get_science_pixel(const Index& idx) const { - return contains(idx) ? (mask.get_pixel(idx) == 0 ? science.get_pixel(idx) : NO_DATA) : NO_DATA; + // The get_pixel() functions perform the bounds checking and will return NO_DATA for out of bounds. + return mask.get_pixel(idx) == 0 ? science.get_pixel(idx) : NO_DATA; } inline float get_variance_pixel(const Index& idx) const { - return contains(idx) ? - (mask.get_pixel(idx) == 0 ? variance.get_pixel(idx) : NO_DATA) : - NO_DATA; + // The get_pixel() functions perform the bounds checking and will return NO_DATA for out of bounds. + return mask.get_pixel(idx) == 0 ? variance.get_pixel(idx) : NO_DATA; } inline bool contains(const Index& idx) const { diff --git a/src/kbmod/search/psi_phi_array.cpp b/src/kbmod/search/psi_phi_array.cpp index 866fe258f..3f2656680 100644 --- a/src/kbmod/search/psi_phi_array.cpp +++ b/src/kbmod/search/psi_phi_array.cpp @@ -177,11 +177,12 @@ void set_encode_cpu_psi_phi_array(PsiPhiArray& data, const std::vector throw std::runtime_error("CPU PsiPhi already allocated."); } if (debug) { - printf("Allocating CPU memory for encoded PsiPhi array using %lu bytes.\n", data.get_total_array_size()); + printf("Allocating CPU memory for encoded PsiPhi array using %lu bytes.\n", + data.get_total_array_size()); } T* encoded = (T*)malloc(data.get_total_array_size()); if (encoded == nullptr) { - throw std::runtime_error("Unable to allocate space for CPU PsiPhi array."); + throw std::runtime_error("Unable to allocate space for CPU PsiPhi array."); } // Create a safe maximum that is slightly less than the true max to avoid @@ -224,7 +225,7 @@ void set_float_cpu_psi_phi_array(PsiPhiArray& data, const std::vector& } float* encoded = (float*)malloc(data.get_total_array_size()); if (encoded == nullptr) { - throw std::runtime_error("Unable to allocate space for CPU PsiPhi array."); + throw std::runtime_error("Unable to allocate space for CPU PsiPhi array."); } int current_index = 0; @@ -266,12 +267,10 @@ void fill_psi_phi_array(PsiPhiArray& result_data, int num_bytes, const std::vect result_data.set_phi_scaling(phi_params[0], phi_params[1], phi_params[2]); if (debug) { - printf("Encoding psi to %i bytes min=%f, max=%f, scale=%f\n", - result_data.get_num_bytes(), psi_params[0], psi_params[1], - psi_params[2]); - printf("Encoding phi to %i bytes min=%f, max=%f, scale=%f\n", - result_data.get_num_bytes(), phi_params[0], phi_params[1], - phi_params[2]); + printf("Encoding psi to %i bytes min=%f, max=%f, scale=%f\n", result_data.get_num_bytes(), + psi_params[0], psi_params[1], psi_params[2]); + printf("Encoding phi to %i bytes min=%f, max=%f, scale=%f\n", result_data.get_num_bytes(), + phi_params[0], phi_params[1], phi_params[2]); } // Do the local encoding. @@ -281,7 +280,9 @@ void fill_psi_phi_array(PsiPhiArray& result_data, int num_bytes, const std::vect set_encode_cpu_psi_phi_array(result_data, psi_imgs, phi_imgs, debug); } } else { - if (debug) { printf("Encoding psi and phi as floats.\n"); } + if (debug) { + printf("Encoding psi and phi as floats.\n"); + } // Just interleave psi and phi images. set_float_cpu_psi_phi_array(result_data, psi_imgs, phi_imgs, debug); } @@ -289,7 +290,8 @@ void fill_psi_phi_array(PsiPhiArray& result_data, int num_bytes, const std::vect #ifdef HAVE_CUDA // Create a copy of the encoded data in GPU memory. if (debug) { - printf("Allocating GPU memory for PsiPhi array using %lu bytes.\n", result_data.get_total_array_size()); + printf("Allocating GPU memory for PsiPhi array using %lu bytes.\n", + result_data.get_total_array_size()); } device_allocate_psi_phi_array(&result_data); if (result_data.get_gpu_array_ptr() == nullptr) { diff --git a/src/kbmod/search/psi_phi_array_utils.h b/src/kbmod/search/psi_phi_array_utils.h index d06a4483c..c363f65d4 100644 --- a/src/kbmod/search/psi_phi_array_utils.h +++ b/src/kbmod/search/psi_phi_array_utils.h @@ -26,7 +26,7 @@ namespace search { std::array compute_scale_params_from_image_vect(const std::vector& imgs, int num_bytes); void fill_psi_phi_array(PsiPhiArray& result_data, int num_bytes, const std::vector& psi_imgs, - const std::vector& phi_imgs, bool debug=false); + const std::vector& phi_imgs, bool debug = false); } /* namespace search */ diff --git a/src/kbmod/search/pydocs/layered_image_docs.h b/src/kbmod/search/pydocs/layered_image_docs.h index cd54decd4..72b924dd7 100644 --- a/src/kbmod/search/pydocs/layered_image_docs.h +++ b/src/kbmod/search/pydocs/layered_image_docs.h @@ -238,11 +238,29 @@ static const auto DOC_LayeredImage_get_variance_pixel = R"doc( )doc"; static const auto DOC_LayeredImage_generate_psi_image = R"doc( - todo + Generates the full psi image where the value of each pixel p in the + resulting image is science[p] / variance[p]. To handle masked bits + apply_mask() must be called before the psi image is generated. Otherwise, + all pixels are used. + Convolves the resulting image with the PSF. + + Returns + ------- + result : `kbmod.RawImage` + A ``RawImage`` of the same dimensions as the ``LayeredImage``. )doc"; static const auto DOC_LayeredImage_generate_phi_image = R"doc( - todo + Generates the full phi image where the value of each pixel p in the + resulting image is 1.0 / variance[p]. To handle masked bits + apply_mask() must be called before the phi image is generated. Otherwise, + all pixels are used. + Convolves the resulting image with the PSF. + + Returns + ------- + result : `kbmod.RawImage` + A ``RawImage`` of the same dimensions as the ``LayeredImage``. )doc"; } // namespace pydocs diff --git a/src/kbmod/search/raw_image.cpp b/src/kbmod/search/raw_image.cpp index 1dcb13214..a2e541355 100644 --- a/src/kbmod/search/raw_image.cpp +++ b/src/kbmod/search/raw_image.cpp @@ -110,8 +110,7 @@ float RawImage::interpolate(const Point& p) const { return total / sumWeights; } -RawImage RawImage::create_stamp(const Point& p, const int radius, - const bool keep_no_data) const { +RawImage RawImage::create_stamp(const Point& p, const int radius, const bool keep_no_data) const { if (radius < 0) throw std::runtime_error("stamp radius must be at least 0"); const int dim = radius * 2 + 1; @@ -123,8 +122,7 @@ RawImage RawImage::create_stamp(const Point& p, const int radius, Image stamp = Image::Constant(dim, dim, NO_DATA); stamp.block(anchor.i, anchor.j, h, w) = image.block(corner.i, corner.j, h, w); - if (!keep_no_data) - stamp = (stamp.array() == NO_DATA).select(0.0, stamp); + if (!keep_no_data) stamp = (stamp.array() == NO_DATA).select(0.0, stamp); return RawImage(stamp); } @@ -580,29 +578,25 @@ static void raw_image_bindings(py::module& m) { .def("set_pixel", &rie::set_pixel, pydocs::DOC_RawImage_set_pixel) .def("set_all", &rie::set_all, pydocs::DOC_RawImage_set_all) // python interface adapters (avoids having to construct Index & Point) - .def( - "get_pixel", - [](rie& cls, int i, int j) { - return cls.get_pixel({i, j}); - }) - .def( - "pixel_has_data", - [](rie& cls, int i, int j) { - return cls.pixel_has_data({i, j}); - }) - .def( - "set_pixel", - [](rie& cls, int i, int j, double val) { - cls.set_pixel({i, j}, val); - }) + .def("get_pixel", + [](rie& cls, int i, int j) { + return cls.get_pixel({i, j}); + }) + .def("pixel_has_data", + [](rie& cls, int i, int j) { + return cls.pixel_has_data({i, j}); + }) + .def("set_pixel", + [](rie& cls, int i, int j, double val) { + cls.set_pixel({i, j}, val); + }) // methods .def("l2_allclose", &rie::l2_allclose, pydocs::DOC_RawImage_l2_allclose) .def("compute_bounds", &rie::compute_bounds, pydocs::DOC_RawImage_compute_bounds) .def("find_peak", &rie::find_peak, pydocs::DOC_RawImage_find_peak) .def("find_central_moments", &rie::find_central_moments, - pydocs::DOC_RawImage_find_central_moments) - .def("center_is_local_max", &rie::center_is_local_max, - pydocs::DOC_RawImage_center_is_local_max) + pydocs::DOC_RawImage_find_central_moments) + .def("center_is_local_max", &rie::center_is_local_max, pydocs::DOC_RawImage_center_is_local_max) .def("create_stamp", &rie::create_stamp, pydocs::DOC_RawImage_create_stamp) .def("interpolate", &rie::interpolate, pydocs::DOC_RawImage_interpolate) .def("interpolated_add", &rie::interpolated_add, pydocs::DOC_RawImage_interpolated_add) @@ -620,21 +614,17 @@ static void raw_image_bindings(py::module& m) { .def("append_fits_extension", &rie::append_to_fits, pydocs::DOC_RawImage_append_to_fits) .def("load_fits", &rie::from_fits, pydocs::DOC_RawImage_load_fits) // python interface adapters - .def( - "create_stamp", - [](rie& cls, float x, float y, int radius, bool keep_no_data) { - return cls.create_stamp({x, y}, radius, keep_no_data); - }) - .def( - "interpolate", - [](rie& cls, float x, float y) { - return cls.interpolate({x, y}); - }) - .def( - "interpolated_add", - [](rie& cls, float x, float y, float val) { - cls.interpolated_add({x, y}, val); - }); + .def("create_stamp", + [](rie& cls, float x, float y, int radius, bool keep_no_data) { + return cls.create_stamp({x, y}, radius, keep_no_data); + }) + .def("interpolate", + [](rie& cls, float x, float y) { + return cls.interpolate({x, y}); + }) + .def("interpolated_add", [](rie& cls, float x, float y, float val) { + cls.interpolated_add({x, y}, val); + }); } #endif diff --git a/tests/test_layered_image.py b/tests/test_layered_image.py index bfc45c368..da8c8fe9a 100644 --- a/tests/test_layered_image.py +++ b/tests/test_layered_image.py @@ -44,6 +44,29 @@ def test_create(self): self.assertGreaterEqual(science.get_pixel(y, x), -100.0) self.assertLessEqual(science.get_pixel(y, x), 100.0) + # Check direct lookup of pixel values matches the RawImage lookup. + self.assertEqual(science.get_pixel(y, x), self.image.get_science_pixel(y, x)) + self.assertEqual(variance.get_pixel(y, x), self.image.get_variance_pixel(y, x)) + + # Check that the LayeredImage pixel lookups work with a masked pixel. + # But the the mask was not applied yet to the images themselves. + mask.set_pixel(5, 6, 1) + self.assertGreater(science.get_pixel(5, 6), KB_NO_DATA) + self.assertGreater(variance.get_pixel(5, 6), KB_NO_DATA) + self.assertEqual(self.image.get_science_pixel(5, 6), KB_NO_DATA) + self.assertEqual(self.image.get_variance_pixel(5, 6), KB_NO_DATA) + + # Test that out of bounds pixel lookups are handled correctly. + self.assertEqual(self.image.get_science_pixel(-1, 1), KB_NO_DATA) + self.assertEqual(self.image.get_science_pixel(1, -1), KB_NO_DATA) + self.assertEqual(self.image.get_science_pixel(self.image.get_height() + 1, 1), KB_NO_DATA) + self.assertEqual(self.image.get_science_pixel(1, self.image.get_width() + 1), KB_NO_DATA) + + self.assertEqual(self.image.get_variance_pixel(-1, 1), KB_NO_DATA) + self.assertEqual(self.image.get_variance_pixel(1, -1), KB_NO_DATA) + self.assertEqual(self.image.get_variance_pixel(self.image.get_height() + 1, 1), KB_NO_DATA) + self.assertEqual(self.image.get_variance_pixel(1, self.image.get_width() + 1), KB_NO_DATA) + def test_create_from_layers(self): sci = RawImage(30, 40) for y in range(40): @@ -285,7 +308,11 @@ def test_psi_and_phi_image(self): for x in range(6): sci.set_pixel(y, x, float(x)) var.set_pixel(y, x, float(y + 1)) + + # Mask a single pixel and set another to variance of zero. + sci.set_pixel(3, 1, KB_NO_DATA) var.set_pixel(3, 1, KB_NO_DATA) + var.set_pixel(3, 2, 0.0) # Generate and check psi and phi images. psi = img.generate_psi_image() @@ -298,12 +325,15 @@ def test_psi_and_phi_image(self): for y in range(5): for x in range(6): - has_data = not (x == 1 and y == 3) + has_data = y != 3 or x == 0 or x > 2 self.assertEqual(psi.pixel_has_data(y, x), has_data) self.assertEqual(phi.pixel_has_data(y, x), has_data) - if x != 1 or y != 3: + if has_data: self.assertAlmostEqual(psi.get_pixel(y, x), x / (y + 1)) self.assertAlmostEqual(phi.get_pixel(y, x), 1.0 / (y + 1)) + else: + self.assertEqual(psi.get_pixel(y, x), KB_NO_DATA) + self.assertEqual(phi.get_pixel(y, x), KB_NO_DATA) def test_subtract_template(self): sci = self.image.get_science()