Merge pull request #5 from GalSim-developers/merge-main-fix-tests-1

test: merge latest changes from main into JAX tests

Author: beckermr

src/SBInterpolatedImage.cpp

@@ -1174,9 +1174,11 @@ namespace galsim {
         for(int ix=0; ix<=max_ix; ++ix) {
             xdbg<<"Start search for ix = "<<ix<<std::endl;
             // Search along the two sides with either kx = ix or ky = ix.
+            double norm_kval = INFINITY; // something we know is above threshold
+
             for(int iy=0; iy<=ix; ++iy) {
                 // The bottom side of the square in the lower-right quadrant.
-                double norm_kval = fast_norm((*_kimage)(iy,-ix));
+                norm_kval = fast_norm((*_kimage)(iy,-ix));
                 xdbg<<"norm_kval at "<<iy<<','<<-ix<<" = "<<norm_kval<<std::endl;
                 if (norm_kval <= thresh && iy != ix && ix != No2) {
                     // The top side of the square in the upper-right quadrant.
@@ -1205,8 +1207,13 @@ namespace galsim {
                 }
             }
             xdbg<<"Done ix = "<<ix<<".  Current count = "<<n_below_thresh<<std::endl;
+            // we can get here either if the loop above finishes normally or if it hits
+            // the break statement. So we have to test the last value of norm_kval to
+            // see if the break statement was hit. If it was not, then we can increment
+            // the counter for the # below thresh.
+            if (norm_kval <= thresh) ++n_below_thresh;
             // If we get through 5 rows with nothing above the threshold, stop looking.
-            if (++n_below_thresh == 5) break;
+            if (n_below_thresh == 5) break;
         }
         xdbg<<"Finished.  maxk_ix = "<<maxk_ix<<std::endl;
         // Add 1 to get the first row that is below the threshold.

tests/test_interpolatedimage.py

@@ -1930,5 +1930,48 @@ def test_drawreal_seg_fault():
         np.testing.assert_array_equal(image.array, 0)
 
 
+@timer
+def test_interpolatedimage_maxk_kspace_pixel_gap():
+    # this code makes an image where there is a gap in the fourier
+    # space image of a certain number pixels where pixels go above
+    # and below the maxk threshold. At >five pixels, galsim should
+    # ignore the gap, but less than that it should increase maxk.
+
+    print("\n| offset | orig       | new                |")
+    print("|--------|------------|--------------------|")
+    for offset in [0, 3, 4, 5, 6, 7]:
+        im = galsim.Gaussian(fwhm=4.5).drawImage(scale=1)
+        iim = galsim.InterpolatedImage(im, scale=1)
+        orig_maxk = iim.maxk
+
+        kim = iim._xim.copy().calculate_fft()
+        kx, ky = kim.get_pixel_centers()
+        kx *= kim.scale
+        ky *= kim.scale
+        # this is the last pixel above threshold. galsim adds 1
+        # to the last pixel it finds above threshold to compute orig_maxk
+        # and so we subtract 1
+        maxk_ix = np.floor(orig_maxk / kim.scale).astype(int) - 1
+        if offset > 0:
+            kim[maxk_ix + offset, maxk_ix] = kim[0, 0].real
+        new_im = kim.calculate_inverse_fft()
+        new_maxk = galsim.InterpolatedImage(new_im, scale=1, pad_factor=1).maxk
+
+        print("| % 6d | %10.6f | %18.6f |" % (
+            offset,
+            orig_maxk,
+            new_maxk),
+        )
+
+        if offset <= 5:
+            # for offsets <=5, we should get an offset of offset pixels
+            # in the location of maxk
+            np.testing.assert_allclose(new_maxk - orig_maxk, offset * kim.scale, atol=1e-12, rtol=0)
+        else:
+            np.testing.assert_allclose(new_maxk, orig_maxk, atol=1e-12, rtol=0)
+
+    print("|--------|------------|--------------------|")
+
+
 if __name__ == "__main__":
     runtests(__file__)