DRAFT PR for discussion

src/kbmod/fake_data_creator.py

@@ -39,15 +39,15 @@ def add_fake_object(img, x, y, flux, psf=None):
         sci = img
 
     if psf is None:
-        sci.add_pixel_interp(x, y, flux)
+        sci.add_pixel(x, y, 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.add_pixel(initial_x + i, initial_y + j, flux * psf.get_value(i, j))
 
 
 class FakeDataSet:
@@ -136,8 +136,8 @@ def insert_object(self, trj):
 
         for i in range(self.num_times):
             dt = self.times[i] - t0
-            px = trj.x + dt * trj.vx + 0.5
-            py = trj.y + dt * trj.vy + 0.5
+            px = trj.get_x_pos(dt)
+            py = trj.get_y_pos(dt)
 
             # Get the image for the timestep, add the object, and
             # re-set the image. This last step needs to be done

src/kbmod/search/common.h

@@ -23,8 +23,8 @@ enum StampType { STAMP_SUM = 0, STAMP_MEAN, STAMP_MEDIAN };
 
 // The position (in pixels) of a trajectory.
 struct PixelPos {
-    float x;
-    float y;
+    int x;
+    int y;
 
     const std::string to_string() const { return "x: " + std::to_string(x) + " y: " + std::to_string(y); }
 
@@ -52,9 +52,9 @@ struct Trajectory {
     short obs_count;
 
     // Get pixel positions from a zero-shifted time.
-    float get_x_pos(float time) const { return x + time * vx; }
-    float get_y_pos(float time) const { return y + time * vy; }
-    PixelPos get_pos(float time) const { return {x + time * vx, y + time * vy}; }
+    int get_x_pos(float time) const { return (int)floor(x + time * vx); }
+    int get_y_pos(float time) const { return (int)floor(y + time * vy); }
+    PixelPos get_pos(float time) const { return { get_x_pos(time), get_y_pos(time)}; }
 
     // I can't believe string::format is not a thing until C++ 20
     const std::string to_string() const {

src/kbmod/search/kernels.cu

@@ -147,8 +147,8 @@ __global__ void searchFilterImages(int num_images, int width, int height, void *
         for (int i = 0; i < num_images; ++i) {
             // Predict the trajectory's position.
             float curr_time = image_data.image_times[i];
-            int current_x = x + int(curr_trj.vx * curr_time + 0.5);
-            int current_y = y + int(curr_trj.vy * curr_time + 0.5);
+            int current_x = x + (int)(floor(curr_trj.vx * curr_time));
+            int current_y = y + (int)(floor(curr_trj.vy * curr_time));
 
             // Test if trajectory goes out of the image, in which case we do not
             // look up a pixel value for this time step (allowing trajectories to
@@ -424,8 +424,8 @@ __global__ void deviceGetCoaddStamp(int num_images, int width, int height, float
 
         // Predict the trajectory's position.
         float curr_time = image_data.image_times[t];
-        int current_x = int(trj.x + trj.vx * curr_time);
-        int current_y = int(trj.y + trj.vy * curr_time);
+        int current_x = trj.x + (int)(floor(trj.vx * curr_time));
+        int current_y = trj.y + (int)(floor(trj.vy * curr_time));
 
         // Get the stamp and add it to the list of values.
         int img_x = current_x - params.radius + stamp_x;

src/kbmod/search/pydocs/common_docs.h

@@ -35,7 +35,7 @@ namespace pydocs {
 
   Returns
   -------
-  `float`
+  `int`
      The predicted x position (in pixels).
   )doc";
 
@@ -49,7 +49,7 @@ namespace pydocs {
 
   Returns
   -------
-  `float`
+  `int`
      The predicted y position (in pixels).
   )doc";
 
@@ -72,10 +72,10 @@ namespace pydocs {
 
   Attributes
   ----------
-    x : `float`
-        An x position on an image (in fractional pixels).
-    y : `float`
-        An x position on an image (in fractional pixels).
+    x : `int`
+        An x position on an image (in pixels).
+    y : `int`
+        An x position on an image (in pixels).
   )doc";
 
   static const auto DOC_ImageMoments = R"doc(

src/kbmod/search/pydocs/raw_image_docs.h

@@ -81,11 +81,6 @@ namespace pydocs{
   Add to the raw value of a given pixel.
   )doc";
 
-  static const auto DOC_RawImage_add_pixel_interp = R"doc(
-  Add to the value calculated by bilinear interpolation
-  of the neighborhood of the given pixel position.
-  )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. Modifies the image in-place.
@@ -121,10 +116,6 @@ namespace pydocs{
   Returns the value of a pixel.
   )doc";
 
-  static const auto DOC_RawImage_get_pixel_interp =  R"doc(
-  Get the interoplated value of a pixel.
-  )doc";
-
   static const auto DOC_RawImage_convolve = R"doc(
   Convolve the image with a PSF.
   )doc";

src/kbmod/search/raw_image.cpp

@@ -198,19 +198,14 @@ void RawImage::append_layer_to_file(const std::string& filename) {
     fits_report_error(stderr, status);
 }
 
-RawImage RawImage::create_stamp(float x, float y, int radius, bool interpolate, bool keep_no_data) const {
+RawImage RawImage::create_stamp(float x, float y, int radius, bool keep_no_data) const {
     if (radius < 0) throw std::runtime_error("stamp radius must be at least 0");
 
     int dim = radius * 2 + 1;
     RawImage stamp(dim, dim);
     for (int yoff = 0; yoff < dim; ++yoff) {
         for (int xoff = 0; xoff < dim; ++xoff) {
-            float pix_val;
-            if (interpolate)
-                pix_val = get_pixel_interp(x + static_cast<float>(xoff - radius),
-                                           y + static_cast<float>(yoff - radius));
-            else
-                pix_val = get_pixel(static_cast<int>(x) + xoff - radius, static_cast<int>(y) + yoff - radius);
+            float pix_val = get_pixel(static_cast<int>(x) + xoff - radius, static_cast<int>(y) + yoff - radius);
             if ((pix_val == NO_DATA) && !keep_no_data) pix_val = 0.0;
             stamp.set_pixel(xoff, yoff, pix_val);
         }
@@ -320,99 +315,13 @@ void RawImage::grow_mask(int steps) {
     }
 }
 
-std::vector<float> RawImage::bilinear_interp(float x, float y) 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};
-}
-
-void RawImage::add_pixel_interp(float x, float y, float value) {
-    // Interpolation values
-    std::vector<float> iv = bilinear_interp(x, y);
-
-    add_to_pixel(iv[0], iv[1], value * iv[2]);
-    add_to_pixel(iv[3], iv[4], value * iv[5]);
-    add_to_pixel(iv[6], iv[7], value * iv[8]);
-    add_to_pixel(iv[9], iv[10], value * iv[11]);
-}
-
 void RawImage::add_to_pixel(float fx, float fy, float value) {
     assert(fx - floor(fx) == 0.0 && fy - floor(fy) == 0.0);
     int x = static_cast<int>(fx);
     int y = static_cast<int>(fy);
     if (x >= 0 && x < width && y >= 0 && y < height) pixels[y * width + x] += value;
 }
 
-float RawImage::get_pixel_interp(float x, float y) const {
-    if ((x < 0.0 || y < 0.0) || (x > static_cast<float>(width) || y > static_cast<float>(height)))
-        return NO_DATA;
-    std::vector<float> iv = bilinear_interp(x, y);
-    float a = get_pixel(iv[0], iv[1]);
-    float b = get_pixel(iv[3], iv[4]);
-    float c = get_pixel(iv[6], iv[7]);
-    float d = get_pixel(iv[9], iv[10]);
-    float interpSum = 0.0;
-    float total = 0.0;
-    if (a != NO_DATA) {
-        interpSum += iv[2];
-        total += a * iv[2];
-    }
-    if (b != NO_DATA) {
-        interpSum += iv[5];
-        total += b * iv[5];
-    }
-    if (c != NO_DATA) {
-        interpSum += iv[8];
-        total += c * iv[8];
-    }
-    if (d != NO_DATA) {
-        interpSum += iv[11];
-        total += d * iv[11];
-    }
-    if (interpSum == 0.0) {
-        return NO_DATA;
-    } else {
-        return total / interpSum;
-    }
-}
-
 void RawImage::set_all_pix(float value) {
     for (auto& p : pixels) p = value;
 }
@@ -664,13 +573,11 @@ static void raw_image_bindings(py::module& m) {
             .def("create_stamp", &ri::create_stamp, pydocs::DOC_RawImage_create_stamp)
             .def("set_pixel", &ri::set_pixel, pydocs::DOC_RawImage_set_pixel)
             .def("add_pixel", &ri::add_to_pixel, pydocs::DOC_RawImage_add_pixel)
-            .def("add_pixel_interp", &ri::add_pixel_interp, pydocs::DOC_RawImage_add_pixel_interp)
             .def("apply_mask", &ri::apply_mask, pydocs::DOC_RawImage_apply_mask)
             .def("grow_mask", &ri::grow_mask, pydocs::DOC_RawImage_grow_mask)
             .def("pixel_has_data", &ri::pixel_has_data, pydocs::DOC_RawImage_pixel_has_data)
             .def("set_all", &ri::set_all_pix, pydocs::DOC_RawImage_set_all)
             .def("get_pixel", &ri::get_pixel, pydocs::DOC_RawImage_get_pixel)
-            .def("get_pixel_interp", &ri::get_pixel_interp, pydocs::DOC_RawImage_get_pixel_interp)
             .def("convolve", &ri::convolve, pydocs::DOC_RawImage_convolve)
             .def("convolve_cpu", &ri::convolve_cpu, pydocs::DOC_RawImage_convolve_cpu)
             .def("load_fits", &ri::load_from_file, pydocs::DOC_RawImage_load_fits)

src/kbmod/search/raw_image.h

@@ -51,10 +51,6 @@ class RawImage {
     const std::vector<float>& get_pixels() const { return pixels; }
     float* data() { return pixels.data(); }  // Get pointer to pixels
 
-    // Get the interpolated brightness of a real values point
-    // using the four neighboring pixels.
-    float get_pixel_interp(float x, float y) const;
-
     // Check if two raw images are approximately equal.
     bool approx_equal(const RawImage& imgB, float atol) const;
 
@@ -74,8 +70,6 @@ class RawImage {
 
     void set_all_pix(float value);
     void add_to_pixel(float fx, float fy, float value);
-    void add_pixel_interp(float x, float y, float value);
-    std::vector<float> bilinear_interp(float x, float y) const;
 
     // Grow the area of masked pixels.
     void grow_mask(int steps);
@@ -95,7 +89,7 @@ class RawImage {
     // 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.
-    RawImage create_stamp(float x, float y, int radius, bool interpolate, bool keep_no_data) const;
+    RawImage create_stamp(float x, float y, int radius, bool keep_no_data) const;
 
     // The maximum value of the image and return the coordinates. The parameter
     // furthest_from_center indicates whether to break ties using the peak further

src/kbmod/search/stack_search.cpp

@@ -325,7 +325,7 @@ std::vector<RawImage> StackSearch::get_coadded_stamps_cpu(std::vector<Trajectory
 
     for (int i = 0; i < num_trajectories; ++i) {
         std::vector<RawImage> stamps =
-                StampCreator::create_stamps(stack, t_array[i], params.radius, false, true, use_index_vect[i]);
+                StampCreator::create_stamps(stack, t_array[i], params.radius, true, use_index_vect[i]);
 
         RawImage coadd(1, 1);
         switch (params.stamp_type) {
@@ -425,7 +425,8 @@ std::vector<float> StackSearch::create_curves(Trajectory t, const std::vector<Ra
         /* Do not use get_pixel_interp(), because results from create_curves must
          * be able to recover the same likelihoods as the ones reported by the
          * gpu search.*/
-        float pix_val = imgs[i].get_pixel(t.x + int(times[i] * t.vx + 0.5), t.y + int(times[i] * t.vy + 0.5));
+        PixelPos p = t.get_pos(times[i]);
+        float pix_val = imgs[i].get_pixel(p.x, p.y);
         if (pix_val == NO_DATA) pix_val = 0.0;
         lightcurve.push_back(pix_val);
     }

src/kbmod/search/stamp_creator.cpp

@@ -4,8 +4,7 @@ namespace search {
 StampCreator::StampCreator() {}
 
 std::vector<RawImage> StampCreator::create_stamps(ImageStack& stack, const Trajectory& trj, int radius,
-                                                  bool interpolate, bool keep_no_data,
-                                                  const std::vector<bool>& use_index) {
+                                                   bool keep_no_data, const std::vector<bool>& use_index) {
     if (use_index.size() > 0 && use_index.size() != stack.img_count()) {
         throw std::runtime_error("Wrong size use_index passed into create_stamps()");
     }
@@ -17,9 +16,9 @@ std::vector<RawImage> StampCreator::create_stamps(ImageStack& stack, const Traje
         if (use_all_stamps || use_index[i]) {
             // Calculate the trajectory position.
             float time = stack.get_zeroed_time(i);
-            PixelPos pos = {trj.x + time * trj.vx, trj.y + time * trj.vy};
+            PixelPos pos = trj.get_pos(time);
             RawImage& img = stack.get_single_image(i).get_science();
-            stamps.push_back(img.create_stamp(pos.x, pos.y, radius, interpolate, keep_no_data));
+            stamps.push_back(img.create_stamp(pos.x, pos.y, radius, keep_no_data));
         }
     }
     return stamps;
@@ -30,31 +29,27 @@ std::vector<RawImage> StampCreator::create_stamps(ImageStack& stack, const Traje
 // individual timesteps have been filtered).
 std::vector<RawImage> StampCreator::get_stamps(ImageStack& stack, const Trajectory& t, int radius) {
     std::vector<bool> empty_vect;
-    return create_stamps(stack, t, radius, true /*=interpolate*/, false /*=keep_no_data*/, empty_vect);
+    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
 // NO_DATA 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, false /*=interpolate*/, true /*=keep_no_data*/, use_index));
+                                       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
 // NO_DATA 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, false /*=interpolate*/, true /*=keep_no_data*/, use_index));
+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 NO_DATA
 // 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 /*=interpolate*/, false /*=keep_no_data*/, use_index));
+                                       const std::vector<bool>& use_index) {
+    return create_summed_image(create_stamps(stack, trj, radius, false /*=keep_no_data*/, use_index));
 }
 
 #ifdef Py_PYTHON_H

src/kbmod/search/stamp_creator.h

@@ -15,13 +15,12 @@ class StampCreator {
     StampCreator();
 
     // Functions science stamps for filtering, visualization, etc. User can specify
-    // the radius of the stamp, whether to interpolate among pixels, whether to keep NO_DATA values
-    // or replace them with zero, and what indices to use.
+    // the radius of the stamp, whether to keep NO_DATA values or replace them with zero,
+    // and what indices to use.
     // The indices to use are indicated by use_index: a vector<bool> indicating whether to use
     // each time step. An empty (size=0) vector will use all time steps.
     static std::vector<RawImage> create_stamps(ImageStack& stack, const Trajectory& trj, int radius,
-                                               bool interpolate, bool keep_no_data,
-                                               const std::vector<bool>& use_index);
+                                        bool keep_no_data, const std::vector<bool>& use_index);
 
     static std::vector<RawImage> get_stamps(ImageStack& stack, const Trajectory& t, int radius);