Remove remaining interpolate functions

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<std::optional<Index>, 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<std::optional<Index>, 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<float, float> coords, const std::pair<unsigned, unsigned> 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

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

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``.

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<float, 4> 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<std::optional<Index>, 4> neighbors;
-        std::array<float, 4> 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

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.

src/kbmod/search/stamp_creator.cpp

@@ -40,22 +40,20 @@ std::vector<RawImage> 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<bool>& 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<bool>& 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<bool>& use_index) {
     return create_summed_image(create_stamps(stack, trj, radius, false /*=keep_no_data*/, use_index));

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

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()

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

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)