diff --git a/src/kbmod/search/bindings.cpp b/src/kbmod/search/bindings.cpp
index 4ffcee70..cdc12932 100644
--- a/src/kbmod/search/bindings.cpp
+++ b/src/kbmod/search/bindings.cpp
@@ -3,6 +3,9 @@
 #include <pybind11/eigen.h>
 #include <pybind11/stl.h>
 
+#include "common.h"
+#include "geom.h"
+
 #include "psf.cpp"
 #include "raw_image.cpp"
 #include "raw_image_eigen.cpp"
@@ -10,10 +13,6 @@
 #include "image_stack.cpp"
 #include "stack_search.cpp"
 #include "filtering.cpp"
-#include "common.h"
-
-using pp = search::PixelPos;
-using std::to_string;
 
 
 PYBIND11_MODULE(search, m) {
@@ -24,6 +23,9 @@ PYBIND11_MODULE(search, m) {
     .value("STAMP_MEAN", search::StampType::STAMP_MEAN)
     .value("STAMP_MEDIAN", search::StampType::STAMP_MEDIAN)
     .export_values();
+  indexing::index_bindings(m);
+  indexing::point_bindings(m);
+  indexing::rectangle_bindings(m);
   search::psf_bindings(m);
   search::raw_image_bindings(m);
   search::raw_image_eigen_bindings(m);
@@ -31,7 +33,6 @@ PYBIND11_MODULE(search, m) {
   search::image_stack_bindings(m);
   search::stack_search_bindings(m);
   search::trajectory_bindings(m);
-  search::index_bindings(m);
   search::pixel_pos_bindings(m);
   search::image_moments_bindings(m);
   search::stamp_parameters_bindings(m);
diff --git a/src/kbmod/search/common.h b/src/kbmod/search/common.h
index 453edf86..5006603c 100644
--- a/src/kbmod/search/common.h
+++ b/src/kbmod/search/common.h
@@ -1,6 +1,8 @@
 #ifndef COMMON_H_
 #define COMMON_H_
 
+// assert(condition, message if !condition)
+#define assertm(exp, msg) assert(((void)msg, exp))
 
 #include <string>
 #include "pydocs/common_docs.h"
diff --git a/src/kbmod/search/geom.h b/src/kbmod/search/geom.h
new file mode 100644
index 00000000..b56eaf16
--- /dev/null
+++ b/src/kbmod/search/geom.h
@@ -0,0 +1,243 @@
+#ifndef GEOM_H_
+#define GEOM_H_
+
+#include <iostream>
+#include <optional>
+#include <array>
+#include <vector>
+#include <assert.h>
+
+#include <Eigen/Core>
+
+#include "common.h"
+//#include "pydocs/geom_docs.h"
+
+
+namespace indexing {
+
+  struct Index{
+    int j;
+    int i;
+
+    const std::string to_string() const {
+      return "Index(" + std::to_string(i) + ", " + std::to_string(j) +")";
+    }
+
+    const std::string to_yaml() const {
+      return "{x: " + std::to_string(i) + " y: " + std::to_string(j) + "}";
+    }
+
+    friend std::ostream& operator<<(std::ostream& os, const Index& rc);
+  };
+
+  std::ostream& operator<<(std::ostream& os, const Index& rc){
+    os << rc.to_string();
+    return os;
+  }
+
+
+  struct Point{
+    float y;
+    float x;
+
+    const Index to_index() const {
+      return {(int)floor(y), (int)floor(x)};
+    }
+
+    const std::string to_string() const {
+      return "Point(" + std::to_string(x) + ", " + std::to_string(y) +")";
+    }
+
+    const std::string to_yaml() const {
+      return "{x: " + std::to_string(x) + " y: " + std::to_string(y) + "}";
+    }
+
+    friend std::ostream& operator<<(std::ostream& os, const Point& rc);
+  };
+
+  std::ostream& operator<<(std::ostream& os, const Point& rc){
+    os << rc.to_string();
+    return os;
+  }
+
+
+  struct Rectangle{
+    Index idx;
+    unsigned width;
+    unsigned height;
+
+    const std::string to_string() const {
+      return "Rectangle(" + idx.to_yaml() + ", dx: " + std::to_string(width) + \
+        ", dy: " + std::to_string(height) + ")";
+    }
+
+    const std::string to_yaml() const {
+      return "{idx: " + idx.to_yaml() +      \
+        " width: " + std::to_string(width) + \
+        " height: " + std::to_string(height);
+    }
+
+    friend std::ostream& operator<<(std::ostream& os, const Rectangle& rc);
+  };
+
+  std::ostream& operator<<(std::ostream& os, const Rectangle& 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};
+  }
+
+
+  inline auto manhattan_neighbors(const Index& idx,
+                                  const unsigned width,
+                                  const unsigned height) {
+    std::array<std::optional<Index>, 4> idxs;
+
+    // top bot
+    if (idx.i >= 0 && idx.i<width){
+      if (idx.j-1 >= 0 && idx.j-1<height)
+        idxs[0] = {idx.j-1, idx.i};
+      if (idx.j+1 >= 0 && idx.j+1<height)
+        idxs[1] = {idx.j+1, idx.i};
+    }
+
+    // left right
+    if (idx.j >= 0 && idx.j<height){
+      if (idx.i-1 >= 0 && idx.i-1<width)
+        idxs[2] = {idx.j, idx.i-1};
+      if (idx.i+1 >= 0 && idx.i+1<width)
+        idxs[3] = {idx.j, idx.i+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.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
+    if (idx.j >= 0 && idx.j<height)
+      if (idx.i+1 >= 0 && idx.i+1<width)
+        idxs[3] = {idx.j, idx.i+1};
+
+    // bot-right
+    if ((idx.j+1 >= 0 && idx.j+1<width) &&
+        (idx.i+1 >= 0 && idx.i+1<width))
+      idxs[2] = {idx.j+1, idx.i+1};
+
+    return idxs;
+  }
+
+
+  /*#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")
+      .def(py::init<int, int>())
+      .def(py::init<float, float>())
+      .def_readwrite("i", &Index::i)
+      .def_readwrite("j", &Index::j)
+      .def("to_yaml", &Index::to_yaml)
+      .def("__repr__", &Index::to_string)
+      .def("__str__", &Index::to_string);
+  }
+
+  static void point_bindings(py::module &m) {
+    py::class_<Point>(m, "Point")
+      .def(py::init<float, float>())
+      .def_readwrite("x", &Point::x)
+      .def_readwrite("y", &Point::y)
+      .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 }; } ))
+      .def(py::init<int, int, float, float>())
+      .def_readwrite("width", &Rectangle::width)
+      .def_readwrite("height", &Rectangle::height)
+      .def_property("i",
+                    /*get*/ [](Rectangle& rect) { return rect.idx.i;},
+                    /*set*/ [](Rectangle& rect, int value) { rect.idx.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);
+  }
+#endif // Py_PYTHON_H
+
+} // namespace indexing
+
+#endif // GEOM_H
diff --git a/src/kbmod/search/raw_image_eigen.cpp b/src/kbmod/search/raw_image_eigen.cpp
index 5bbadf36..4b8b8e66 100644
--- a/src/kbmod/search/raw_image_eigen.cpp
+++ b/src/kbmod/search/raw_image_eigen.cpp
@@ -2,6 +2,10 @@
 
 
 namespace search {
+  using Index = indexing::Index;
+  using Point = indexing::Point;
+  using Rect = indexing::Rectangle;
+
   RawImageEigen::RawImageEigen()
     : width(0),
       height(0),
@@ -73,139 +77,103 @@ namespace search {
   }
 
 
-  std::array<float, 12> RawImageEigen::get_interp_neighbors_and_weights(float y, float x) const {
-    // Linear interpolation
-    // Find the 4 pixels (aPix, bPix, cPix, dPix)
-    // that the corners (a, b, c, d) of the
-    // new pixel land in, and blend into those
-
-    // Returns a vector with 4 pixel locations
-    // and their interpolation value
-
-    // Top right
-    float ax = x + 0.5;
-    float ay = y + 0.5;
-    float a_px = floor(ax);
-    float a_py = floor(ay);
-    float a_amt = (ax - a_px) * (ay - a_py);
-
-    // Bottom right
-    float bx = x + 0.5;
-    float by = y - 0.5;
-    float b_px = floor(bx);
-    float b_py = floor(by);
-    float b_amt = (bx - b_px) * (b_py + 1.0 - by);
-
-    // Bottom left
-    float cx = x - 0.5;
-    float cy = y - 0.5;
-    float c_px = floor(cx);
-    float c_py = floor(cy);
-    float c_amt = (c_px + 1.0 - cx) * (c_py + 1.0 - cy);
-
-    // Top left
-    float dx = x - 0.5;
-    float dy = y + 0.5;
-    float d_px = floor(dx);
-    float d_py = floor(dy);
-    float d_amt = (d_px + 1.0 - dx) * (dy - d_py);
-
-    // make sure the right amount has been distributed
-    float diff = std::abs(a_amt + b_amt + c_amt + d_amt - 1.0);
-    if (diff > 0.01) std::cout << "warning: bilinear_interpSum == " << diff << "\n";
-    return {a_px, a_py, a_amt, b_px, b_py, b_amt, c_px, c_py, c_amt, d_px, d_py, d_amt};
+  inline auto RawImageEigen::get_interp_neighbors_and_weights(const Point& p) const {
+    // top/bot left; top/bot right
+    auto neighbors = indexing::manhattan_neighbors(p, width, height);
+
+    float tmp_integral;
+    float x_frac = std::modf(p.x, &tmp_integral);
+    float y_frac = std::modf(p.y, &tmp_integral);
+
+    // weights for top/bot left; top/bot right
+    std::array<float, 4> weights{
+      (1-x_frac) * (1-y_frac) +
+         x_frac  * (1-y_frac) +
+      (1-x_frac) *    y_frac  +
+         x_frac  *    y_frac
+    };
+
+    struct NeighborsAndWeights{
+      std::array<std::optional<Index>, 4> neighbors;
+      std::array<float, 4> weights;
+    };
+
+    return NeighborsAndWeights{neighbors, weights};
   }
 
 
-  float RawImageEigen::interpolate(const float y, const float x) const {
-    if ((x < 0.0 || y < 0.0) || (x > static_cast<float>(width) || y > static_cast<float>(height)))
+  float RawImageEigen::interpolate(const Point& p) const {
+    if (!contains(p))
       return NO_DATA;
 
-    auto [a_x, a_y, w_a, b_x, b_y, w_b, c_x, c_y, w_c, d_x, d_y, w_d] = get_interp_neighbors_and_weights(y, x);
+    auto [neighbors, weights] = get_interp_neighbors_and_weights(p);
 
-    float a = get_pixel(a_y, a_x);
-    float b = get_pixel(b_y, b_x);
-    float c = get_pixel(c_y, c_x);
-    float d = get_pixel(d_y, d_x);
-    float interpSum = 0.0;
+    // 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;
-    if (a != NO_DATA) {
-      interpSum += w_a;
-      total += a * w_a;
-    }
-    if (b != NO_DATA) {
-      interpSum += w_b;
-      total += b * w_b;
-    }
-    if (c != NO_DATA) {
-      interpSum += w_c;
-      total += c * w_c;
-    }
-    if (d != NO_DATA) {
-      interpSum += w_d;
-      total += d * w_d;
+    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->j, idx->i);
+        }
+      }
     }
-    if (interpSum == 0.0) {
+
+    if (sumWeights == 0.0)
       return NO_DATA;
-    } else {
-      return total / interpSum;
-    }
+    return total / sumWeights;
   }
 
 
-  Image RawImageEigen::create_stamp(const float y, const float x, const int radius,
-                                    const bool interpolate, const bool keep_no_data) const {
+  Image RawImageEigen::create_stamp(const Point& p,
+                                    const int radius,
+                                    const bool interpolate,
+                                    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;
-    Index idx(x, y, false);
-    auto [i, dimx, j, dimy] = idx.centered_block(radius, width, height);
-    Image stamp = image.block(j, i, dimy, dimx);
+    auto [idx, dimx, dimy] = indexing::centered_block(p.to_index(), radius,
+                                                      width, height);
+    Image stamp = image.block(idx.j, idx.i, dimy, dimx);
 
     if (interpolate) {
       for (int yoff = 0; yoff < dim; ++yoff) {
         for (int xoff = 0; xoff < dim; ++xoff) {
-          stamp(yoff, xoff) = this->interpolate(y + static_cast<float>(yoff - radius), x + static_cast<float>(xoff - radius));
+          // I think the {} create a temporary, but I don't know how bad that is
+          // would it be the same if we had interpolate just take 2 floats?
+          stamp(yoff, xoff) = this->interpolate({
+              p.y + static_cast<float>(yoff - radius),
+              p.x + static_cast<float>(xoff - radius)
+            });
         }
       }
     }
-
     return stamp;
   }
 
 
-  void RawImageEigen::add(const float x, const float y, const float value) {
-    // In the original implementation this just does array[x, y] += value
-    // even though that doesn't really make any sense? The cast is probably
-    // implicit. This is likely slower than before, because the floor is used
-    // not a cast. Floor is more correct, but I added the constructor nonetheless
-    Index p(x, y, false);
-    if (contains(p))
-      image(p.i, p.j) += value;
-  }
-
-
-  void RawImageEigen::add(const int x, const int y, const float value) {
-    Index p(x, y);
-    if (contains(p))
-      image(p.i, p.j) += value;
+  inline void RawImageEigen::add(const Index& idx, const float value) {
+    if (contains(idx))
+      image(idx.i, idx.j) += value;
   }
 
 
-  void RawImageEigen::add(const Index p, const float value) {
-    if (contains(p))
-      image(p.i, p.j) += value;
+  inline void RawImageEigen::add(const Point& p, const float value) {
+    add(p.to_index(), value);
   }
 
 
-  void RawImageEigen::interpolated_add(const float x, const float y, const float value) {
-    // interpolation weights and neighbors
-    auto [a_x, a_y, w_a, b_x, b_y, w_b, c_x, c_y, w_c, d_x, d_y, w_d] = get_interp_neighbors_and_weights(x, y);
-    add(a_x, a_y, value * w_a);
-    add(b_x, b_y, value * w_b);
-    add(c_x, c_y, value * w_c);
-    add(d_x, d_y, value * w_d);
+  void RawImageEigen::interpolated_add(const Point& p, const float value) {
+    auto [neighbors, weights] = get_interp_neighbors_and_weights(p);
+    for (int i=0; i<neighbors.size(); i++)
+      if (auto& idx = neighbors[i])
+        add(*idx, weights[i]);
   }
 
 
@@ -227,16 +195,6 @@ namespace search {
   }
 
 
-  // multiple times now I've encountered the issue where I can't reduce
-  // simple array operations to method calls in Eigen because masking is
-  // baked into the image in a way that it doesn't represent the identity
-  // of the operation in question. This forces us to unravel the arrays
-  // in for loops and check if NO_DATA. I wonder if resorting to Eigens
-  // expressions would be faster than this. I also wonder how "dangerous"
-  // it would be. How often, let's say, do we expect to see exactly 0 in real data?
-  // Then convolution could be as easy as:
-  // Image zeroed = (array.array() == NO_DATA).select(0, array.array());
-  // result(i, j) = zeroed.block<krows, kcols>(i, j).cwiseProduct(kernel).sum()
   void RawImageEigen::convolve_cpu(PSF& psf) {
     Image result = Image::Zero(height, width);
 
@@ -257,6 +215,7 @@ namespace search {
           for (int i = -psf_rad; i <= psf_rad; i++) {
             if ((x + i >= 0) && (x + i < width) && (y + j >= 0) && (y + j < height)) {
               float current_pixel = image(y+j, x+i);
+              // note that convention for index access is flipped for PSF
               if (current_pixel != NO_DATA) {
                 float current_psf = psf.get_value(i + psf_rad, j + psf_rad);
                 psf_portion += current_psf;
@@ -624,7 +583,8 @@ namespace search {
         float sum = 0.0;
         float count = 0.0;
         for (int i = 0; i < num_images; ++i) {
-          float pix_val = images[i].get_pixel(y, x);
+          // again, temporary is created? How bad is that?
+          float pix_val = images[i].get_pixel({y, x});
           // of course...
           if ((pix_val != NO_DATA) && (!std::isnan(pix_val))) {
             count += 1.0;
@@ -632,13 +592,11 @@ namespace search {
           }
         }
 
-        if (count > 0.0) {
-          result.set_pixel(y, x, sum / count);
-        } else {
-          // We use a 0.0 value if there is no data to allow for visualization
-          // and value based filtering.
-          result.set_pixel(y, x, 0.0);
-        }
+        // again, temporaries in the set_pixel, how bad is this
+        if (count > 0.0)
+          result.set_pixel({y, x}, sum / count);
+        else
+          result.set_pixel({y, x}, 0.0); // use 0 for visualization purposes
       } // for x
     } // for y
     return result;
@@ -646,21 +604,6 @@ namespace search {
 
 
 #ifdef Py_PYTHON_H
-  static void index_bindings(py::module &m) {
-    py::class_<Index>(m, "Index")
-      .def(py::init<int, int>())
-      .def(py::init<float, float>())
-      .def_readwrite("i", &Index::i)
-      .def_readwrite("j", &Index::j)
-      .def("__iter__", [](const Index &idx){
-        std::vector<unsigned> tmp{idx.i, idx.j};
-        return py::make_iterator(tmp.begin(), tmp.end());
-      }, py::keep_alive<0, 1>())
-      .def("__repr__", [] (const Index &idx) { return idx.to_string(); })
-      .def("__str__", &Index::to_string);
-  }
-
-
   static void raw_image_eigen_bindings(py::module &m) {
     using rie = search::RawImageEigen;
 
@@ -671,21 +614,34 @@ namespace search {
            py::arg("img").noconvert(true), py::arg("obs_time")=-1.0d)
       .def(py::init<unsigned, unsigned, float, double>(),
            py::arg("w"), py::arg("h"), py::arg("value")=0.0f, py::arg("obs_time")=-1.0d)
+      // attributes and properties
       .def_property_readonly("height", &rie::get_height)
       .def_property_readonly("width", &rie::get_width)
       .def_property_readonly("npixels", &rie::get_npixels)
       .def_property("obstime", &rie::get_obstime, &rie::set_obstime)
       .def_property("image", py::overload_cast<>(&rie::get_image, py::const_), &rie::set_image)
       .def_property("imref", py::overload_cast<>(&rie::get_image), &rie::set_image)
+      // pixel accessors and setters
       .def("get_pixel", &rie::get_pixel)
       .def("pixel_has_data", &rie::pixel_has_data)
       .def("set_pixel", &rie::set_pixel)
       .def("set_all", &rie::set_all)
+      // python interface adapters (avoids having to construct Index & Point)
+      .def("get_pixel", [](rie& cls, int j, int i){
+        return cls.get_pixel({j, i});
+      })
+      .def("pixel_has_data", [](rie& cls, int j, int i){
+        return cls.pixel_has_data({j, i});
+      })
+      .def("set_pixel", [](rie& cls, int j, int i, float val){
+        cls.set_pixel({j, i}, val);
+      })
+      // methods
       .def("l2_allclose", &rie::l2_allclose)
       .def("compute_bounds", &rie::compute_bounds)
       .def("find_peak", &rie::find_peak)
       .def("find_central_moments", &rie::find_central_moments)
-      .def("create_stamp", &rie::create_stamp, py::return_value_policy::copy)
+      .def("create_stamp", &rie::create_stamp)
       .def("interpolate", &rie::interpolate)
       .def("interpolated_add", &rie::interpolated_add)
       .def("apply_mask", &rie::apply_mask)
@@ -694,7 +650,12 @@ namespace search {
       .def("convolve_cpu", &rie::convolve_cpu)
       .def("save_fits", &rie::save_fits)
       .def("append_fits_extension", &rie::append_fits_extension)
-      .def("load_fits", &rie::load_fits);
+      .def("load_fits", &rie::load_fits)
+      // python interface adapters
+      .def("create_stamp", [](rie& cls, float y, float x, int radius,
+                              bool interp, bool keep_no_data){
+        return cls.create_stamp({y, x}, radius, interp, keep_no_data );
+      });
   }
 #endif
 
diff --git a/src/kbmod/search/raw_image_eigen.h b/src/kbmod/search/raw_image_eigen.h
index c104d519..cfe4ba32 100644
--- a/src/kbmod/search/raw_image_eigen.h
+++ b/src/kbmod/search/raw_image_eigen.h
@@ -1,103 +1,25 @@
 #ifndef RAWIMAGEEIGEN_H_
 #define RAWIMAGEEIGEN_H_
 
-
 #include <vector>
 #include <fitsio.h>
 #include <float.h>
 #include <iostream>
 #include <stdexcept>
+#include <assert.h>
+
 #include <Eigen/Core>
 
 #include "common.h"
+#include "geom.h"
 #include "psf.h"
 #include "pydocs/raw_image_docs.h"
 
 
 namespace search {
-
-  struct Index {
-    unsigned i;
-    unsigned j;
-
-    Index(unsigned x, unsigned y)
-      : i(x), j(y) {}
-
-    Index(int x, int y)
-      : i(x), j(y) {}
-
-    Index(float x, float y)
-      : i(floor(x)), j(floor(y)) {}
-
-    Index(float x, float y, bool floor)
-      : i(static_cast<unsigned>(x)), j(static_cast<int>(y)) {}
-
-    const std::string to_string() const {
-      return "Index(" + std::to_string(i) + ", " + std::to_string(j) +")";
-    }
-
-    const std::string to_yaml() const {
-      return "{x: " + std::to_string(i) + " y: " + std::to_string(j) + "}";
-    }
-
-    std::array<unsigned, 4> centered_block(const unsigned size, const unsigned width, const unsigned height) const {
-      unsigned d = size;
-      unsigned left_x = ((i-d >= 0) && (i-d < width)) ? i-d : i;
-      unsigned right_x = ((i+d >= 0) && (i+d < width)) ? i+d : width - i;
-      unsigned top_y = ((j-d >= 0) && (j-d < height)) ? j-d : j;
-      unsigned bot_y = ((j+d >= 0) && (j+d < height)) ? j+d : height - i;
-
-      unsigned dx = right_x - left_x + 1;
-      unsigned dy = bot_y - top_y + 1;
-      return {left_x, dx, top_y, dy};
-    }
-
-    friend std::ostream& operator<<(std::ostream& os, const Index& rc);
-  };
-
-  std::ostream& operator<<(std::ostream& os, const Index& rc){
-    os << rc.to_string();
-    return os;
-  }
-
-
-  struct Point{
-    float x;
-    float y;
-
-    Point(float xd, float yd)
-      : x(xd), y(yd) {}
-
-    Index to_index(){
-      return Index(x, y);
-    }
-
-    std::array<Point, 4> nearest_pixel_coords()  {
-      return {{
-        {(float)floor(x-0.5f)+0.5f, (float)floor(y+0.5f)+0.5f},
-        {(float)floor(x+0.5f)+0.5f, (float)floor(y+0.5f)+0.5f},
-        {(float)floor(x-0.5f)+0.5f, (float)floor(y-0.5f)+0.5f},
-        {(float)floor(x+0.5f)+0.5f, (float)floor(y-0.5f)+0.5f}
-        }};
-    }
-
-    std::array<Index, 4> nearest_pixel_idxs(){
-      return {{
-        {x-0.5f, y+0.5f},
-        {x+0.5f, y+0.5f},
-        {x-0.5f, y-0.5f},
-        {x+0.5f, y-0.5f}
-        }};
-    };
-
-    friend std::ostream& operator<<(std::ostream& os, const Point& rc);
-  };
-
-  std::ostream& operator<<(std::ostream& os, const Point& rc){
-    os << "x: " << rc.x << " y: " << rc.y;
-    return os;
-  }
-
+  using Index = indexing::Index;
+  using Point = indexing::Point;
+  using Rect = indexing::Rectangle;
 
   using Image = Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
   using ImageI = Eigen::Matrix<int, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
@@ -125,28 +47,29 @@ namespace search {
     Image& get_image() { return image; }
     void set_image(Image& other) { image = other; }
 
-    // Inline pixel functions.
-    float get_pixel(int y, int x) const {
-      return (x >= 0 && x < width && y >= 0 && y < height) ? image(y, x) : NO_DATA;
-    }
 
-    bool pixel_has_data(int y, int x) const {
-      return (x >= 0 && x < width && y >= 0 && y < height) ? image(y, x) != NO_DATA : false;
+    inline bool contains(const Index& idx) const {
+      return idx.i>=0 && idx.i<width && idx.j>=0 && idx.j<height;
     }
 
-    bool contains(const float x, const float y) const{
-      return x>=0 && x<=width && y>=0 && y<=height;
+    inline bool contains(const Point& p) const {
+      return p.x>=0 && p.y<width && p.y>=0 && p.y<height;
     }
 
-    bool contains(const int i, const int j) const{
-      return i>=0 && i<=width && j>=0 && j<=height;
+    inline float get_pixel(const Index& idx) const {
+      return contains(idx) ? image(idx.j, idx.i) : NO_DATA;
     }
 
-    bool contains(const Index idx) const{
-      return idx.i>=0 && idx.i<=width && idx.j>=0 && idx.j<=height;
+    bool pixel_has_data(const Index& idx) const {
+      return get_pixel(idx) != NO_DATA ? true : false;
     }
 
-    void set_pixel(unsigned j, unsigned i, float value) { image(j, i) = value; }
+    void set_pixel(const Index& idx, float value) {
+      // we should probably be letting Eigen freak out about setting an impossible
+      // index instead of silently just nod doing it; but this is how it is
+      if (contains(idx))
+        image(idx.j, idx.i) = value;
+    }
 
     // Functions for locally storing the image time.
     double get_obstime() const { return obstime; }
@@ -162,22 +85,19 @@ namespace search {
 
     // Get the interpolated brightness of a real values point
     // using the four neighboring array.
-    std::array<float, 12> get_interp_neighbors_and_weights(const float x, const float y) const;
-    float interpolate(const float x, const float y) const;
-    float interpolate2(const float x, const float y) const;
+    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.
-    Image create_stamp(const float x, float y,
-                       const int radius, const bool interpolate,
-                       const bool keep_no_data) const;
+    Image create_stamp(const Point& p,  const int radius,
+                       const bool interpolate, const bool keep_no_data) const;
 
     // pixel modifiers
-    void add(const float x, const float y, const float value);
-    void add(const int x, const int y, const float value);
-    void add(const Index p, const float value);
-    void interpolated_add(const float x, const float y, const float value);
+    void add(const Index& idx, const float value);
+    void add(const Point& p, const float value);
+    void interpolated_add(const Point& p, const float value);
 
     // Compute the min and max bounds of values in the image.
     std::array<float, 2> compute_bounds() const;