Merge raw_image_eigen into codebase, fix tests.

src/kbmod/analysis_utils.py

@@ -325,7 +325,7 @@ def apply_stamp_filter(
                 kb.HAS_GPU and len(trj_slice) > 100,
             )
             for ind, stamp in enumerate(stamps_slice):
-                if stamp.get_width() > 1:
+                if stamp.width > 1:
                     result_list.results[ind + start_idx].stamp = np.array(stamp)
                     all_valid_inds.append(ind + start_idx)
 

src/kbmod/fake_data_creator.py

@@ -17,7 +17,7 @@
 from kbmod.work_unit import WorkUnit
 
 
-def add_fake_object(img, x, y, flux, psf=None):
+def add_fake_object(img, y, x, flux, psf=None):
     """Add a fake object to a LayeredImage or RawImage
 
     Parameters
@@ -39,21 +39,20 @@ def add_fake_object(img, x, y, flux, psf=None):
         sci = img
 
     if psf is None:
-        sci.add_pixel_interp(x, y, flux)
+        sci.interpolated_add(y, x, 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))
+                sci.interpolated_add(float(initial_y + j), float(initial_x + i), flux * psf.get_value(i, j))
 
 
 class FakeDataSet:
     """This class creates fake data sets for testing and demo notebooks."""
 
-    def __init__(self, width, height, num_times, noise_level=2.0, psf_val=0.5, obs_per_day=3, use_seed=False):
+    def __init__(self, height, width, num_times, noise_level=2.0, psf_val=0.5, obs_per_day=3, use_seed=False):
         """The constructor.
 
         Parameters

src/kbmod/filters/stamp_filters.py

@@ -106,8 +106,8 @@ def keep_row(self, row: ResultRow):
         stamp = row.stamp.reshape([self.width, self.width])
         peak_pos = RawImage(stamp).find_peak(True)
         return (
-            abs(peak_pos.x - self.stamp_radius) < self.x_thresh
-            and abs(peak_pos.y - self.stamp_radius) < self.y_thresh
+            abs(peak_pos.i - self.stamp_radius) < self.x_thresh
+            and abs(peak_pos.j - self.stamp_radius) < self.y_thresh
         )
 
 

src/kbmod/search/bindings.cpp

@@ -10,7 +10,6 @@ namespace py = pybind11;
 
 #include "psf.cpp"
 #include "raw_image.cpp"
-#include "raw_image_eigen.cpp"
 #include "layered_image.cpp"
 #include "image_stack.cpp"
 #include "stack_search.cpp"
@@ -27,10 +26,9 @@ PYBIND11_MODULE(search, m) {
     .export_values();
   indexing::index_bindings(m);
   indexing::point_bindings(m);
-  indexing::rectangle_bindings(m);
+  indexing::anchored_rectangle_bindings(m);
   search::psf_bindings(m);
   search::raw_image_bindings(m);
-  search::raw_image_eigen_bindings(m);
   search::layered_image_bindings(m);
   search::image_stack_bindings(m);
   search::stack_search_bindings(m);
@@ -42,10 +40,6 @@ PYBIND11_MODULE(search, m) {
   m.def("create_median_image", &search::create_median_image);
   m.def("create_summed_image", &search::create_summed_image);
   m.def("create_mean_image", &search::create_mean_image);
-  // Functions from raw_image_eigen.cpp
-  m.def("create_median_image_eigen", &search::create_median_image_eigen);
-  m.def("create_summed_image_eigen", &search::create_summed_image_eigen);
-  m.def("create_mean_image_eigen", &search::create_mean_image_eigen);
   // Functions from filtering.cpp
   m.def("sigmag_filtered_indices", &search::sigmaGFilteredIndices);
   m.def("calculate_likelihood_psi_phi", &search::calculateLikelihoodFromPsiPhi);

src/kbmod/search/common.h

@@ -1,11 +1,15 @@
 #ifndef COMMON_H_
 #define COMMON_H_
 
+#include <assert.h>
+#include <string>
+
+#include "pydocs/common_docs.h"
+
+
 // assert(condition, message if !condition)
 #define assertm(exp, msg) assert(((void)msg, exp))
 
-#include <string>
-#include "pydocs/common_docs.h"
 
 namespace search {
 #ifdef HAVE_CUDA

src/kbmod/search/geom.h

@@ -3,6 +3,7 @@
 
 #include <iostream>
 #include <optional>
+#include <utility>  // pair
 #include <array>
 #include <vector>
 #include <assert.h>
@@ -20,11 +21,11 @@ namespace indexing {
     int i;
 
     const std::string to_string() const {
-      return "Index(" + std::to_string(i) + ", " + std::to_string(j) +")";
+      return "Index(" + std::to_string(j) + ", " + std::to_string(i) +")";
     }
 
     const std::string to_yaml() const {
-      return "{x: " + std::to_string(i) + " y: " + std::to_string(j) + "}";
+      return "{j: " + std::to_string(j) + " i: " + std::to_string(i) + "}";
     }
 
     friend std::ostream& operator<<(std::ostream& os, const Index& rc);
@@ -41,15 +42,15 @@ namespace indexing {
     float x;
 
     const Index to_index() const {
-      return {(int)floor(y), (int)floor(x)};
+      return {(int)(floor(y-0.5f)+0.5f), (int)(floor(x-0.5f)+0.5f)};
     }
 
     const std::string to_string() const {
-      return "Point(" + std::to_string(x) + ", " + std::to_string(y) +")";
+      return "Point(" + std::to_string(y) + ", " + std::to_string(x) +")";
     }
 
     const std::string to_yaml() const {
-      return "{x: " + std::to_string(x) + " y: " + std::to_string(y) + "}";
+      return "{x: " + std::to_string(y) + " y: " + std::to_string(x) + "}";
     }
 
     friend std::ostream& operator<<(std::ostream& os, const Point& rc);
@@ -61,44 +62,76 @@ namespace indexing {
   }
 
 
-  struct Rectangle{
-    Index idx;
-    unsigned width;
+  // A rectangle that also contains it's corner's origin Index with respect to
+  // another, larger rectangle.
+  struct AnchoredRectangle{
+    Index corner;
+    Index anchor;
+
     unsigned height;
+    unsigned width;
 
     const std::string to_string() const {
-      return "Rectangle(" + idx.to_yaml() + ", dx: " + std::to_string(width) + \
-        ", dy: " + std::to_string(height) + ")";
+      return "AnchoredRectangle(corner:" + corner.to_string() +
+        ", anchor:" + anchor.to_string() +
+        ", height: " + std::to_string(height) +
+        ", width: " + std::to_string(width) + ")";
     }
 
     const std::string to_yaml() const {
-      return "{idx: " + idx.to_yaml() +      \
-        " width: " + std::to_string(width) + \
-        " height: " + std::to_string(height);
+      return "{corner: " + corner.to_yaml() +
+        " anchor: " + anchor.to_yaml() +
+        " height: " + std::to_string(height) +
+        " width: " + std::to_string(width) + "}";
     }
 
-    friend std::ostream& operator<<(std::ostream& os, const Rectangle& rc);
+    friend std::ostream& operator<<(std::ostream& os, const AnchoredRectangle& rc);
   };
 
-  std::ostream& operator<<(std::ostream& os, const Rectangle& rc){
+  std::ostream& operator<<(std::ostream& os, const AnchoredRectangle& rc){
     os << rc.to_string();
     return os;
   }
 
 
-  inline Rectangle centered_block(const Index& idx,
-                                  const int r,
-                                  const unsigned width,
-                                  const unsigned height) {
-    int left_x = ((idx.i-r >= 0) && (idx.i-r < width)) ? idx.i-r : idx.i;
-    int right_x = ((idx.i+r >= 0) && (idx.i+r < width)) ? idx.i+r : width - idx.i;
-    int top_y = ((idx.j-r >= 0) && (idx.j-r < height)) ? idx.j-r : idx.j;
-    int bot_y = ((idx.j+r >= 0) && (idx.j+r < height)) ? idx.j+r : height - idx.i;
-    assertm(bot_y - top_y > 0, "Rectangle has negative height!");
-    assertm(right_x - left_x > 0, "Rectangle has negative width!");
-    unsigned dx = right_x - left_x + 1;
-    unsigned dy = bot_y - top_y + 1;
-    return {{top_y, left_x}, dy, dx};
+  // Given an index component `idx_val`, and a radius `r` around it, returns the
+  // start and end index components, and length of the range, that fit in the
+  // [0, max_range] limits.
+  inline std::tuple<int, int, int> centered_range(int idx_val, const int r,
+                                                  const int max_range) {
+    // pin start to the []0, max_range] range
+    int start = std::max(0, idx_val-r);
+    start = std::min(start, max_range);
+
+    // pin end to the [0, max_range]
+    int end = std::max(0, idx_val+r);
+    end = std::min(max_range, idx_val+r);
+
+    // range is inclusive of the first element, and can not be longer than start
+    // minus max range
+    int length = end-start+1;
+    length = std::min(length, max_range-start);
+
+    return std::make_tuple(start, end, length);
+  }
+
+
+  // get an Eigen block coordinates (top-left, height, width) anchored inside a
+  // square matrix of dimensions 2*r+1 (anchor top-left).
+  inline AnchoredRectangle anchored_block(const Index& idx,
+                                          const int r,
+                                          const unsigned height,
+                                          const unsigned width) {
+    auto [top, bot, rangey] = centered_range(idx.j, r, height);
+    auto [left, right, rangex] = centered_range(idx.i, r, width);
+    assertm(rangey > 0, "Selected block lies outside of the image limits.");
+    assertm(rangex > 0, "Selected block lies outside of the image limits.");
+
+    int anchor_top = std::max(r-idx.j, 0);
+    int anchor_left = std::max(r-idx.i, 0);
+
+    // now it's safe to cast ranges to unsigned
+    return {{top, left}, {anchor_top, anchor_left}, (unsigned)rangey, (unsigned)rangex};
   }
 
 
@@ -142,15 +175,15 @@ namespace indexing {
     // point is negative or (0, 0)
     auto idx = p.to_index();
 
-    // top-left bot-left
+    // top-left bot-right
     if (idx.i >= 0 && idx.i<width){
       if (idx.j >= 0 && idx.j<height)
         idxs[0] = {idx.j, idx.i};
       if (idx.j+1 >= 0 && idx.j+1<height)
         idxs[1] = {idx.j+1, idx.i};
     }
 
-    // top-right
+    // bot-right
     if (idx.j >= 0 && idx.j<height)
       if (idx.i+1 >= 0 && idx.i+1<width)
         idxs[3] = {idx.j, idx.i+1};
@@ -164,38 +197,6 @@ namespace indexing {
   }
 
 
-  /*#ifndef NDEBUG
-  // these are helper functions not used in the code, but help with debugging
-  inline std::vector<Index> all_neighbors(const Index& idx,
-  const unsigned width,
-  const unsigned height) {
-  auto res = manhattan_neighbors(idx, width, height);
-
-  // top-left
-  if ((idx.i-1 >= 0 && idx.i-1<width ? true : false) &&
-  (idx.j-1 >= 0 && idx.j-1<height ? true : false))
-  idxs.push_back(Index(idx.i-1, idx.j-1));
-
-  // top-right
-  if ((idx.i+1 >= 0 && idx.i+1<width ? true : false) &&
-  (idx.j-1 >= 0 && idx.j-1<height ? true : false))
-  idxs.push_back(Index(idx.i+1, idx.j-1));
-
-  // bot left
-  if ((idx.i-1 >= 0 && idx.i-1<width ? true : false) &&
-  (idx.j+1 >= 0 && idx.j+1<height ? true : false))
-  idxs.push_back(Index(idx.i-1, idx.j+1));
-
-  // bot right
-  if ((idx.i+1 >= 0 && idx.i+1<width ? true : false) &&
-  (idx.j+1 >= 0 && idx.j+1<height ? true : false))
-  idxs.push_back(Index(idx.i+1, idx.j+11));
-
-  return idxs;
-  }
-  #endif // NDEBUG
-  */
-
 #ifdef Py_PYTHON_H
   static void index_bindings(py::module &m) {
     py::class_<Index>(m, "Index")
@@ -213,28 +214,35 @@ namespace indexing {
       .def(py::init<float, float>())
       .def_readwrite("x", &Point::x)
       .def_readwrite("y", &Point::y)
+      .def("to_index", &Point::to_index)
       .def("to_yaml", &Point::to_yaml)
       .def("__repr__", &Point::to_string)
       .def("__str__", &Point::to_string);
   }
 
-  static void rectangle_bindings(py::module &m) {
-    py::class_<Rectangle>(m, "Rectangle")
-      .def(py::init<Index, unsigned, unsigned>())
-      .def(py::init( [](int i, int j, unsigned width, unsigned height)
-      { return Rectangle{{i, j}, width, height }; } ))
+  static void anchored_rectangle_bindings(py::module &m) {
+    py::class_<AnchoredRectangle>(m, "AnchoredRectangle")
+      .def(py::init<Index, Index, unsigned, unsigned>())
+      .def(py::init( [](std::pair<int, int> corner, std::pair<int, int> anchor, unsigned height, unsigned width) {
+        return AnchoredRectangle{{corner.first, corner.second}, {anchor.first, anchor.second}, height, width};
+      }))
+      .def(py::init( [](std::pair<int, int> corner, unsigned height, unsigned width) {
+        return AnchoredRectangle{{corner.first, corner.second}, {0, 0}, height, width};
+      }))
       .def(py::init<int, int, float, float>())
-      .def_readwrite("width", &Rectangle::width)
-      .def_readwrite("height", &Rectangle::height)
+      .def_readwrite("width", &AnchoredRectangle::width)
+      .def_readwrite("height", &AnchoredRectangle::height)
+      .def_readwrite("corner", &AnchoredRectangle::corner)
+      .def_readwrite("anchor", &AnchoredRectangle::anchor)
       .def_property("i",
-                    /*get*/ [](Rectangle& rect) { return rect.idx.i;},
-                    /*set*/ [](Rectangle& rect, int value) { rect.idx.i = value; })
+                    /*get*/ [](AnchoredRectangle& rect) { return rect.corner.i;},
+                    /*set*/ [](AnchoredRectangle& rect, int value) { rect.corner.i = value; })
       .def_property("j",
-                    /*get*/ [](Rectangle& rect) { return rect.idx.j;},
-                    /*set*/ [](Rectangle& rect, int value) { rect.idx.j = value; })
-      .def("to_yaml", &Rectangle::to_yaml)
-      .def("__repr__", &Rectangle::to_string)
-      .def("__str__", &Rectangle::to_string);
+                    /*get*/ [](AnchoredRectangle& rect) { return rect.corner.j;},
+                    /*set*/ [](AnchoredRectangle& rect, int value) { rect.corner.j = value; })
+      .def("to_yaml", &AnchoredRectangle::to_yaml)
+      .def("__repr__", &AnchoredRectangle::to_string)
+      .def("__str__", &AnchoredRectangle::to_string);
   }
 #endif // Py_PYTHON_H