Merge pull request #390 from dirac-institute/cleanups

Remove interpolation from stamp generation

src/kbmod/search/pydocs/raw_image_docs.h

@@ -187,8 +187,6 @@ static const auto DOC_RawImage_create_stamp = R"doc(
       The y value of the center of the stamp.
   radius : `int`
       The stamp radius. Width = 2*radius+1.
-  interpolate : `bool`
-      A Boolean indicating whether to interpolate pixel values.
   keep_no_data : `bool`
       A Boolean indicating whether to preserve NO_DATA tags or to
       replace them with 0.0.

src/kbmod/search/pydocs/stamp_creator_docs.h

@@ -19,8 +19,6 @@ static const auto DOC_StampCreator_create_stamps = R"doc(
       The trajectory to project to each time.
   radius : `int`
       The stamp radius. Width = 2*radius+1.
-  interpolate : `bool`
-      A Boolean indicating whether to interpolate pixel values.
   keep_no_data : `bool`
       A Boolean indicating whether to preserve NO_DATA tags or to
       replace them with 0.0.

src/kbmod/search/raw_image.cpp

@@ -110,7 +110,7 @@ float RawImage::interpolate(const Point& p) const {
     return total / sumWeights;
 }
 
-RawImage RawImage::create_stamp(const Point& p, const int radius, const bool interpolate,
+RawImage RawImage::create_stamp(const Point& p, const int radius,
                                 const bool keep_no_data) const {
     if (radius < 0) throw std::runtime_error("stamp radius must be at least 0");
 
@@ -120,23 +120,12 @@ RawImage RawImage::create_stamp(const Point& p, const int radius, const bool int
     // the pixel grid to coordinate system transformation.
     auto [corner, anchor, w, h] = indexing::anchored_block({(int)p.y, (int)p.x}, radius, width, height);
 
-    Image stamp;
-    if (keep_no_data)
-        stamp = Image::Constant(dim, dim, NO_DATA);
-    else
-        stamp = Image::Zero(dim, dim);
+    Image stamp = Image::Constant(dim, dim, NO_DATA);
     stamp.block(anchor.i, anchor.j, h, w) = image.block(corner.i, corner.j, h, w);
 
-    if (interpolate) {
-        for (int yoff = 0; yoff < dim; ++yoff) {
-            for (int xoff = 0; xoff < dim; ++xoff) {
-                // 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.x + static_cast<float>(xoff - radius), p.y + static_cast<float>(yoff - radius)});
-            }
-        }
-    }
+    if (!keep_no_data)
+        stamp = (stamp.array() == NO_DATA).select(0.0, stamp);
+
     return RawImage(stamp);
 }
 
@@ -229,7 +218,7 @@ void RawImage::convolve(PSF psf) {
 
 void RawImage::apply_mask(int flags, const RawImage& mask) {
     for (unsigned int j = 0; j < height; ++j) {
-	    for (unsigned int i = 0; i < width; ++i) {
+        for (unsigned int i = 0; i < width; ++i) {
             int pix_flags = static_cast<int>(mask.image(j, i));
             if ((flags & pix_flags) != 0) {
                 image(j, i) = NO_DATA;
@@ -635,8 +624,8 @@ static void raw_image_bindings(py::module& m) {
             // python interface adapters
             .def(
                     "create_stamp",
-                    [](rie& cls, float x, float y, int radius, bool interp, bool keep_no_data) {
-                        return cls.create_stamp({x, y}, radius, interp, keep_no_data);
+                    [](rie& cls, float x, float y, int radius, bool keep_no_data) {
+                        return cls.create_stamp({x, y}, radius, keep_no_data);
                     })
             .def(
                     "interpolate",

src/kbmod/search/raw_image.h

@@ -79,8 +79,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(const Point& p, const int radius, const bool interpolate,
-                          const bool keep_no_data) const;
+    RawImage create_stamp(const Point& p, const int radius, const bool keep_no_data) const;
 
     // pixel modifiers
     void add(const Index& idx, const float value);

src/kbmod/search/stamp_creator.cpp

@@ -11,8 +11,7 @@ void deviceGetCoadds(const unsigned int num_images, const unsigned int width, co
 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()");
     }
@@ -26,42 +25,42 @@ std::vector<RawImage> StampCreator::create_stamps(ImageStack& stack, const Traje
             float time = stack.get_zeroed_time(i);
             Point pos{trj.x + time * trj.vx, trj.y + time * trj.vy};
             RawImage& img = stack.get_single_image(i).get_science();
-            stamps.push_back(img.create_stamp(pos, radius, interpolate, keep_no_data));
+            stamps.push_back(img.create_stamp(pos, radius, keep_no_data));
         }
     }
     return stamps;
 }
 
-// For stamps used for visualization we interpolate the pixel values, replace
-// NO_DATA tages with zeros, and return all the stamps (regardless of whether
-// individual timesteps have been filtered).
+// For stamps used for visualization we replace NO_DATA tages with zeros
+// and return all the stamps (regardless of whether 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));
+            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));
+            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));
+            create_stamps(stack, trj, radius, false /*=keep_no_data*/, use_index));
 }
 
 std::vector<RawImage> StampCreator::get_coadded_stamps(ImageStack& stack, std::vector<Trajectory>& t_array,
@@ -86,7 +85,7 @@ std::vector<RawImage> StampCreator::get_coadded_stamps_cpu(ImageStack& stack,
 
     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) {

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

tests/test_raw_image.py

@@ -371,10 +371,34 @@ def test_grow_mask(self):
     def test_make_stamp(self):
         """Test stamp creation."""
         img = RawImage(self.array)
-        stamp = img.create_stamp(2.5, 2.5, 2, True, False)
+        stamp = img.create_stamp(2.5, 2.5, 2, False)
         self.assertEqual(stamp.image.shape, (5, 5))
         self.assertTrue((stamp.image == self.array[0:5, 0:5]).all())
 
+        # Test a stamp that is not at the corner.
+        stamp = img.create_stamp(8.5, 5.5, 1, False)
+        self.assertEqual(stamp.image.shape, (3, 3))
+        self.assertTrue((stamp.image == self.array[4:7, 7:10]).all())
+
+        # Test a stamp with NO_DATA.
+        img2 = RawImage(self.masked_array)
+        stamp = img2.create_stamp(7.5, 5.5, 1, True)
+        self.assertEqual(stamp.image.shape, (3, 3))
+        self.assertTrue((stamp.image == self.masked_array[4:7, 6:9]).all())
+
+        # Test a stamp with NO_DATA and replacement.
+        img2 = RawImage(self.masked_array)
+        stamp = img2.create_stamp(7.5, 5.5, 2, False)
+        self.assertEqual(stamp.image.shape, (5, 5))
+        expected = np.copy(self.masked_array[3:8, 5:10])
+        expected[expected == KB_NO_DATA] = 0.0
+        self.assertTrue((stamp.image == expected).all())
+
+        # Test a stamp that goes out of bounds.
+        stamp = img.create_stamp(0.5, 11.5, 1, False)
+        expected = np.array([[0.0, 100.0, 101.0], [0.0, 110.0, 111.0], [0.0, 0.0, 0.0]])
+        self.assertTrue((stamp.image == expected).all())
+
     def test_read_write_file(self):
         """Test file writes and reads correctly."""
         img = RawImage(self.array, 10.0)

tests/test_search.py

@@ -263,9 +263,9 @@ def test_sci_viz_stamps(self):
 
             # Compute the interpolated pixel value at the projected location.
             t = times[i]
-            x = float(self.trj.x) + self.trj.vx * t
-            y = float(self.trj.y) + self.trj.vy * t
-            pixVal = self.imlist[i].get_science().interpolate(x, y)
+            x = int(self.trj.x + self.trj.vx * t)
+            y = int(self.trj.y + self.trj.vy * t)
+            pixVal = self.imlist[i].get_science().get_pixel(y, x)
             if pixVal == KB_NO_DATA:
                 pivVal = 0.0