From 4cfb5ed9fbf3cc6ee3555d03afdb7ab8a51919d1 Mon Sep 17 00:00:00 2001 From: Jeremy Kubica <104161096+jeremykubica@users.noreply.github.com> Date: Mon, 16 Dec 2024 09:28:57 -0500 Subject: [PATCH] Remove remaining interpolate functions --- src/kbmod/search/geom.h | 67 ------------------------ src/kbmod/search/pydocs/geom_docs.h | 43 --------------- src/kbmod/search/pydocs/raw_image_docs.h | 34 ++---------- src/kbmod/search/raw_image.cpp | 65 +---------------------- src/kbmod/search/raw_image.h | 5 -- src/kbmod/search/stamp_creator.cpp | 10 ++-- tests/test_end_to_end.py | 2 +- tests/test_geom.py | 33 ------------ tests/test_image_stack.py | 4 -- tests/test_stamp_creator.py | 2 +- 10 files changed, 11 insertions(+), 254 deletions(-) diff --git a/src/kbmod/search/geom.h b/src/kbmod/search/geom.h index dfbfc874e..679e64bc5 100644 --- a/src/kbmod/search/geom.h +++ b/src/kbmod/search/geom.h @@ -157,56 +157,6 @@ inline Rectangle anchored_block(const Index& idx, const int r, const unsigned wi return {{top, left}, {anchor_top, anchor_left}, (unsigned)rangej, (unsigned)rangei}; } -// returns top-right-bot-left (clockwise) corners around an Index. -inline auto manhattan_neighbors(const Index& idx, const unsigned width, const unsigned height) { - std::array, 4> idxs; - - // if conditions are not met, the element is not set. Because of optional - // type, casting a not-set element to bool returns false. Pybind11 - // evaluates it as a None. See `interpolate` in RawImage for an example. - // top bot - if (idx.j >= 0 && idx.j < width) { - if (idx.i - 1 >= 0 && idx.i - 1 < height) idxs[0] = {idx.i - 1, idx.j}; - if (idx.i + 1 >= 0 && idx.i + 1 < height) idxs[2] = {idx.i + 1, idx.j}; - } - - // right left - if (idx.i >= 0 && idx.i < height) { - if (idx.j + 1 >= 0 && idx.j + 1 < width) idxs[1] = {idx.i, idx.j + 1}; - if (idx.j - 1 >= 0 && idx.j - 1 < width) idxs[3] = {idx.i, idx.j - 1}; - } - return idxs; -} - -// Note the distinct contextual distinction between manhattan neighborhood of -// Index and Point. Point returns closes **pixel indices** that are neighbors. -// This includes the pixel the Point is residing within. This is not the case -// for Index, which will never return itself as a neighbor. -inline auto manhattan_neighbors(const Point& p, const unsigned width, const unsigned height) { - std::array, 4> idxs; - - // The index in which the point resides. - // Almost always top-left corner, except when - // point is negative or (0, 0) - auto idx = p.to_index(); - - // top-left bot-left - if (idx.j >= 0 && idx.j < width) { - if (idx.i >= 0 && idx.i < height) idxs[0] = {idx.i, idx.j}; - if (idx.i + 1 >= 0 && idx.i + 1 < height) idxs[3] = {idx.i + 1, idx.j}; - } - - // top-right - if (idx.i >= 0 && idx.i < height) - if (idx.j + 1 >= 0 && idx.j + 1 < width) idxs[1] = {idx.i, idx.j + 1}; - - // bot-right - if ((idx.i + 1 >= 0 && idx.i + 1 < width) && (idx.j + 1 >= 0 && idx.j + 1 < width)) - idxs[2] = {idx.i + 1, idx.j + 1}; - - return idxs; -} - #ifdef Py_PYTHON_H static void index_bindings(py::module& m) { PYBIND11_NUMPY_DTYPE(Index, i, j); @@ -312,23 +262,6 @@ static void geom_functions(py::module& m) { return anchored_block({idx.first, idx.second}, r, shape.second, shape.first); }, pydocs::DOC_anchored_block); - - // Safe to cast to int as "ij" implies user is using indices, i.e. ints. - // CPP can be so odd, why not return a 1 or, you know... a bool? Mostly for - // testing purposes. - m.def( - "manhattan_neighbors", - [](const std::pair coords, const std::pair shape, - const std::string indexing) { - if (indexing.compare("ij") == 0) - return manhattan_neighbors(Index{(int)coords.first, (int)coords.second}, shape.second, - shape.first); - else if (indexing.compare("xy") == 0) - return manhattan_neighbors(Point{coords.first, coords.second}, shape.second, shape.first); - else - throw std::domain_error("Expected 'ij' or 'xy' got " + indexing + " instead."); - }, - pydocs::DOC_manhattan_neighbors); } #endif // Py_PYTHON_H } // namespace indexing diff --git a/src/kbmod/search/pydocs/geom_docs.h b/src/kbmod/search/pydocs/geom_docs.h index a28460250..1759d6723 100644 --- a/src/kbmod/search/pydocs/geom_docs.h +++ b/src/kbmod/search/pydocs/geom_docs.h @@ -154,49 +154,6 @@ static const auto DOC_anchored_block = R"doc( [-1., 10., 11.]]) )doc"; -static const auto DOC_manhattan_neighbors = R"doc( - Returns a list of nearest neighbors as determined by Manhattan distance of 1. - - Indexing scheme ``ij`` handles the input as `Index`, i.e. the closest - neighbors are top, right, bot, and left indices when not on the edge of an - array. When ``xy``, the input is treated as a `Point`, i.e. real Cartesian - coordinates. In this case the closest neighbors are the closest pixels. Pixels - are array elements understood to span the whole integer range and which center - coordinates lie on a half-integer grid. - For example, origin of an image is the pixel with index ``(0, 0)``, it spans - the range ``[0, 1]`` and its center is the point ``(0.5, 0.5)``. So the - closest neighbor to a `Point(0.6, 0.6)` is `Index(0, 0)` and `Point(1, 1)` is - equidistant from indices ``[(0, 0), (0, 1), (1, 0), (1, 1)]``. - - Parameters - ---------- - coords : `tuple` - Center, around which the neighbors will be returned. - shape : `tuple` - Dimensions of the image/array. - indexing : `str` - Indexing scheme, ``ij`` or ``xy``. - - Returns - ------- - neighbors : `list[Index | None]` - List of indices in clockwise order starting from top or top-left (as - appropriate), or `None` when the returned `Index` would lie outside of the - array/ - - Raises - ------ - ValueError - when indexing is not ``ij`` or ``xy`` - - Examples - -------- - >>> shape = (10, 10) - >>> manhattan_neighbors((0, 0), shape, "ij") - [None, (0, 1), (1, 0), None] - >>> manhattan_neighbors((1, 1), shape, "xy") - [(0, 0), (0, 1), (1, 1), (1, 0)] - )doc"; - } // namespace pydocs #endif // GEOM_DOCS diff --git a/src/kbmod/search/pydocs/raw_image_docs.h b/src/kbmod/search/pydocs/raw_image_docs.h index c1ddf197e..2253bbf55 100644 --- a/src/kbmod/search/pydocs/raw_image_docs.h +++ b/src/kbmod/search/pydocs/raw_image_docs.h @@ -46,9 +46,9 @@ static const auto DOC_RawImage = R"doc( representing indices to a 2D matrix, usually expressed with the ``(i, j)`` convention. Pixel accessing or setting methods of `RawImage`, such as `get_pixel`, use the indexing convention. This matches NumPy convention. Other - methods, such as `interpolate` or `add_fake_object`, however, use the `(x, y)` - convention which is the reversed NumPy convention. Refer to individual methods - signature and docstring to see which one they use. + methods, such as `add_fake_object`, however, use the `(x, y)` convention which + is the reversed NumPy convention. Refer to individual methods signature and docstring + to see which one they use. Examples -------- @@ -264,34 +264,6 @@ static const auto DOC_RawImage_create_stamp = R"doc( The stamp. )doc"; -static const auto DOC_RawImage_interpolate = R"doc( - Get the interoplated value of a point. - - Parameters - ---------- - x : `float` - The x-coordinate, the abscissa. - y : `float` - The y-coordinate, the ordinate. - - Returns - ------- - value : `float` - Bilinearly interpolated value at that point. - - )doc"; - -static const auto DOC_RawImage_get_interp_neighbors_and_weights = R"doc( - Returns a tuple of Manhattan neighbors and interpolation weights. - - Parameters - ---------- - x : `float` - The x coordinate at which to add value. - y : `float` - Y coordinate. - )doc"; - static const auto DOC_RawImage_apply_mask = R"doc( Applies a mask to the RawImage by comparing the given bit vector with the values in the mask layer and marking pixels ``NO_DATA``. diff --git a/src/kbmod/search/raw_image.cpp b/src/kbmod/search/raw_image.cpp index 57e808b91..72360f81a 100644 --- a/src/kbmod/search/raw_image.cpp +++ b/src/kbmod/search/raw_image.cpp @@ -63,55 +63,6 @@ bool RawImage::l2_allclose(const RawImage& img_b, float atol) const { return true; } -inline auto RawImage::get_interp_neighbors_and_weights(const Point& p) const { - // top-left, top right, bot-right, bot-left - auto neighbors = indexing::manhattan_neighbors(p, width, height); - - float tmp_integral; - float x_frac = std::modf(p.x - 0.5f, &tmp_integral); - float y_frac = std::modf(p.y - 0.5f, &tmp_integral); - - // weights for top-left, top right, bot-right, bot-left - std::array weights{ - (1 - x_frac) * (1 - y_frac), - x_frac * (1 - y_frac), - x_frac * y_frac, - (1 - x_frac) * y_frac, - }; - - struct NeighborsAndWeights { - std::array, 4> neighbors; - std::array weights; - }; - - return NeighborsAndWeights{neighbors, weights}; -} - -float RawImage::interpolate(const Point& p) const { - if (!contains(p)) return NO_DATA; - - auto [neighbors, weights] = get_interp_neighbors_and_weights(p); - - // this is the way apparently the way it's supposed to be done - // too bad the loop couldn't have been like - // for (auto& [neighbor, weight] : std::views::zip(neigbors, weights)) - // https://en.cppreference.com/w/cpp/ranges/zip_view - // https://en.cppreference.com/w/cpp/utility/optional - float sumWeights = 0.0; - float total = 0.0; - for (int i = 0; i < neighbors.size(); i++) { - if (auto idx = neighbors[i]) { - if (pixel_has_data(*idx)) { - sumWeights += weights[i]; - total += weights[i] * image(idx->i, idx->j); - } - } - } - - if (sumWeights == 0.0) return NO_DATA; - return total / sumWeights; -} - void RawImage::replace_masked_values(float value) { for (unsigned int y = 0; y < height; ++y) { for (unsigned int x = 0; x < width; ++x) { @@ -494,24 +445,12 @@ static void raw_image_bindings(py::module& m) { 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( - "get_interp_neighbors_and_weights", - [](rie& cls, float x, float y) { - auto tmp = cls.get_interp_neighbors_and_weights({x, y}); - return py::make_tuple(tmp.neighbors, tmp.weights); - }, - pydocs::DOC_RawImage_get_interp_neighbors_and_weights) .def("apply_mask", &rie::apply_mask, pydocs::DOC_RawImage_apply_mask) .def("convolve_gpu", &rie::convolve, pydocs::DOC_RawImage_convolve_gpu) .def("convolve_cpu", &rie::convolve_cpu, pydocs::DOC_RawImage_convolve_cpu) // 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("create_stamp", [](rie& cls, float x, float y, int radius, bool keep_no_data) { + return cls.create_stamp({x, y}, radius, keep_no_data); }); } #endif diff --git a/src/kbmod/search/raw_image.h b/src/kbmod/search/raw_image.h index c98bfd08e..69262bfdb 100644 --- a/src/kbmod/search/raw_image.h +++ b/src/kbmod/search/raw_image.h @@ -78,11 +78,6 @@ class RawImage { // (NaNs) as equal if they appear in both images. bool l2_allclose(const RawImage& imgB, float atol) const; - // Get the interpolated brightness of a real values point - // using the four neighboring array. - inline auto get_interp_neighbors_and_weights(const Point& p) const; - float interpolate(const Point& p) const; - // Create a "stamp" image of a give radius (width=2*radius+1) about the // given point. // keep_no_data indicates whether to use the NO_DATA flag or replace with 0.0. diff --git a/src/kbmod/search/stamp_creator.cpp b/src/kbmod/search/stamp_creator.cpp index c5d1341af..6ff9df36e 100644 --- a/src/kbmod/search/stamp_creator.cpp +++ b/src/kbmod/search/stamp_creator.cpp @@ -40,22 +40,20 @@ std::vector StampCreator::get_stamps(ImageStack& stack, const Trajecto return create_stamps(stack, t, radius, false /*=keep_no_data*/, empty_vect); } -// For creating coadded stamps, we do not interpolate the pixel values and keep -// invalid pixels tagged (so we can filter it out of mean/median). +// For creating coadded stamps, we keep invalid pixels tagged (so we can filter it out of mean/median). RawImage StampCreator::get_median_stamp(ImageStack& stack, const Trajectory& trj, int radius, const std::vector& use_index) { return create_median_image(create_stamps(stack, trj, radius, true /*=keep_no_data*/, use_index)); } -// For creating coadded stamps, we do not interpolate the pixel values and keep -// invalid pixels tagged (so we can filter it out of mean/median). +// For creating coadded stamps, we keep invalid pixels tagged (so we can filter it out of mean/median). RawImage StampCreator::get_mean_stamp(ImageStack& stack, const Trajectory& trj, int radius, const std::vector& use_index) { return create_mean_image(create_stamps(stack, trj, radius, true /*=keep_no_data*/, use_index)); } -// For creating summed stamps, we do not interpolate the pixel values and replace -// invalid pixels with zero (which is the same as filtering it out for the sum). +// For creating summed stamps, we replace invalid pixels with zero (which is the same as +// filtering it out for the sum). RawImage StampCreator::get_summed_stamp(ImageStack& stack, const Trajectory& trj, int radius, const std::vector& use_index) { return create_summed_image(create_stamps(stack, trj, radius, false /*=keep_no_data*/, use_index)); diff --git a/tests/test_end_to_end.py b/tests/test_end_to_end.py index ec6aa20ba..1c955f281 100644 --- a/tests/test_end_to_end.py +++ b/tests/test_end_to_end.py @@ -15,7 +15,7 @@ # this is the first test to actually test things like get_all_stamps from # analysis utils. For now stamps have to be RawImages (because methods like -# interpolate and convolve are defined to work on RawImage and not as funciton) +# convolve are defined to work on RawImage and not as funciton) # so it makes sense to duplicate all this functionality to return np arrays # (instead of RawImages), but hopefully we can deduplicate all this by making # these operations into functions and calling on the .image attribute diff --git a/tests/test_geom.py b/tests/test_geom.py index 6b3efe00a..f48fb1b6e 100644 --- a/tests/test_geom.py +++ b/tests/test_geom.py @@ -8,7 +8,6 @@ Rectangle, centered_range, anchored_block, - manhattan_neighbors, ) @@ -203,38 +202,6 @@ def test_anchored_block(self): ) self.assertEqual(rect.anchor, Index(2, 0)) - def test_manhattan_neighbors(self): - """Test returned values and ordering of manhattan neighbors.""" - shape = (5, 5) - - # fmt: off - # First for Index - # top right bot left - topleft = [None, (0, 1), (1, 0), None] - topright = [None, None, (1, 4), (0, 3)] - botleft = [(3, 0), (4, 1), None, None] - botright = [(3, 4), None, None, (4, 3)] - first = [(0, 1), (1, 2), (2, 1), (1, 0)] - self.assertEqual(topleft, manhattan_neighbors((0, 0), shape, "ij")) - self.assertEqual(topright, manhattan_neighbors((0, 4), shape, "ij")) - self.assertEqual(botright, manhattan_neighbors((4, 4), shape, "ij")) - self.assertEqual(botleft, manhattan_neighbors((4, 0), shape, "ij")) - self.assertEqual(first, manhattan_neighbors((1, 1), shape, "ij")) - - # then for Point - note the semantic difference of "neighboring" - # top-l top-r bot-r bot-l - topleft = [(0, 0), (0, 1), (1, 1), (1, 0)] - topright = [(3, 0), (3, 1), (4, 1), (4, 0)] - botleft = [(0, 3), (0, 4), (1, 4), (1, 3)] - botright = [(3, 3), (3, 4), (4, 4), (4, 3)] - first = [(0, 0), (0, 1), (1, 1), (1, 0)] - self.assertEqual(topleft, manhattan_neighbors((0, 0), shape, "xy")) - self.assertEqual(topright, manhattan_neighbors((0, 4), shape, "xy")) - self.assertEqual(botleft, manhattan_neighbors((4, 0), shape, "xy")) - self.assertEqual(botright, manhattan_neighbors((4, 4), shape, "xy")) - self.assertEqual(first, manhattan_neighbors((1, 1), shape, "xy")) - # fmt: on - if __name__ == "__main__": unittest.main() diff --git a/tests/test_image_stack.py b/tests/test_image_stack.py index a3eb346da..99405a04a 100644 --- a/tests/test_image_stack.py +++ b/tests/test_image_stack.py @@ -115,10 +115,6 @@ def test_make_global_mask(self): else: self.assertEqual(mask.get_pixel(y, x), 0) - # WOW, this is the first test that caught the fact that interpolated_add - # called add, and that add had flipped i and j by accident. The first one. - # TODO: more clean understandable tests for basic functionality, these big - # are super hard to debug.... def test_different_psfs(self): # Add a stationary fake object to each image. Then test that # the flux at each time is monotonically increasing (because diff --git a/tests/test_stamp_creator.py b/tests/test_stamp_creator.py index 58a9567ab..7a5099f2f 100644 --- a/tests/test_stamp_creator.py +++ b/tests/test_stamp_creator.py @@ -218,7 +218,7 @@ def test_sci_viz_stamps(self): self.assertEqual(sci_stamps[i].width, 5) self.assertEqual(sci_stamps[i].height, 5) - # Compute the interpolated pixel value at the projected location. + # Compute the pixel value at the projected location. x = self.trj.get_x_index(times[i]) y = self.trj.get_y_index(times[i]) pixVal = self.imlist[i].get_science().get_pixel(y, x)