Remove pixel interpolation functions

src/kbmod/fake_data/fake_data_creator.py

@@ -84,9 +84,9 @@ def add_fake_object(img, x, y, flux, psf=None):
     ----------
     img : `RawImage` or `LayeredImage`
         The image to modify.
-    x : `float`
+    x : `int` or `float`
         The x pixel location of the fake object.
-    y : `float`
+    y : `int` or `float`
         The y pixel location of the fake object.
     flux : `float`
         The flux value.
@@ -98,15 +98,24 @@ def add_fake_object(img, x, y, flux, psf=None):
     else:
         sci = img
 
+    # Explicitly cast to float because the indexing uses different order
+    # float integer and float.
+    x = int(x)
+    y = int(y)
+
     if psf is None:
-        sci.interpolated_add(x, y, flux)
+        if (x >= 0) and (y >= 0) and (x < sci.width) and (y < sci.height):
+            sci.set_pixel(y, x, flux + sci.get_pixel(y, x))
     else:
         dim = psf.get_dim()
-        initial_x = x - psf.get_radius()
-        initial_y = y - psf.get_radius()
+        initial_x = int(x - psf.get_radius())
+        initial_y = int(y - psf.get_radius())
         for i in range(dim):
             for j in range(dim):
-                sci.interpolated_add(float(initial_x + i), float(initial_y + j), flux * psf.get_value(i, j))
+                xp = initial_x + i
+                yp = initial_y + j
+                if (xp >= 0) and (yp >= 0) and (xp < sci.width) and (yp < sci.height):
+                    sci.set_pixel(yp, xp, flux * psf.get_value(i, j) + sci.get_pixel(yp, xp))
 
 
 def create_fake_times(num_times, t0=0.0, obs_per_day=1, intra_night_gap=0.01, inter_night_gap=1):

src/kbmod/fake_data/insert_fake_orbit.py

@@ -1,5 +1,6 @@
 from astropy.wcs import WCS
 
+from kbmod.fake_data.fake_data_creator import add_fake_object
 from kbmod.fake_data.pyoorb_helper import PyoorbOrbit
 from kbmod.search import ImageStack, LayeredImage, PSF, RawImage
 from kbmod.work_unit import WorkUnit
@@ -10,7 +11,7 @@ def safe_add_fake_detection(img, x, y, flux):
 
     Parameters
     ----------
-    img : `RawImage` or `LayeredImage`
+    img : `LayeredImage`
         The image to modify.
     x : `float`
         The x pixel location of the fake object.
@@ -35,16 +36,7 @@ def safe_add_fake_detection(img, x, y, flux):
     if img.get_mask().get_pixel(int(y), int(x)) != 0:
         return False
 
-    sci = img.get_science()
-    psf = img.get_psf()
-    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.interpolated_add(float(initial_x + i), float(initial_y + j), flux * psf.get_value(i, j))
-
+    add_fake_object(img.get_science(), x, y, flux, img.get_psf())
     return True
 
 

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
   --------
@@ -198,7 +198,7 @@ static const auto DOC_RawImage_compute_bounds = R"doc(
   bounds : `tuple`
       A ``(min, max)`` tuple.
   )doc";
-  
+
 static const auto DOC_RawImage_find_peak = R"doc(
   Returns the pixel coordinates of the maximum value.
 
@@ -264,52 +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_interpolated_add = R"doc(
-  Add the given value at a given point, to the image.
-
-  Addition is performed by determining the nearest Manhattan neighbors, weighing
-  the given value by the distance to these neighbors and then adding that value
-  at the index locations of the neighbors. Sort of like an inverse bilinear
-  interpolation.
-
-  Parameters
-  ----------
-  x : `float`
-      The x coordinate at which to add value.
-  y : `float`
-      Y coordinate.
-  value : `float`
-      Value to add.
-  )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/pydocs/stack_search_docs.h

@@ -169,7 +169,7 @@ static const auto DOC_StackSearch_get_number_total_results = R"doc(
   result : `int`
       The number of cached results.
   )doc";
-                                        
+
 static const auto DOC_StackSearch_get_results = R"doc(
   Get a batch of cached results.
 

src/kbmod/search/pydocs/trajectory_list_docs.h

@@ -143,7 +143,6 @@ static const auto DOC_TrajectoryList_sort_by_likelihood = R"doc(
   Raises a ``RuntimeError`` the data is on GPU.
   )doc";
 
-
 }  // namespace pydocs
 
 #endif /* TRAJECTORY_LIST_DOCS_ */