From 369ecc6fbdd0cabcd1057c27e0259d4cfb3b186e Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Wed, 29 Apr 2026 13:39:37 -0500
Subject: [PATCH 01/15] fix: bug in maxk algorithm for interpolated

---
 src/SBInterpolatedImage.cpp | 11 +++++++++--
 1 file changed, 9 insertions(+), 2 deletions(-)

diff --git a/src/SBInterpolatedImage.cpp b/src/SBInterpolatedImage.cpp
index c92672fe06..51ca6293a0 100644
--- a/src/SBInterpolatedImage.cpp
+++ b/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.

From 86d164e1d60a986f1fa9404dc88d78f056834ed1 Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Wed, 29 Apr 2026 16:38:34 -0500
Subject: [PATCH 02/15] fix: add a test

---
 .github/workflows/ci.yml        | 11 +++--
 tests/test_interpolatedimage.py | 82 +++++++++++++++++++++++++++++++++
 2 files changed, 88 insertions(+), 5 deletions(-)

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 6fda7ddaaf..6d27de48a8 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -201,11 +201,12 @@ jobs:
 
             - name: Run unit tests
               run: |
-                cd tests
-                coverage run -m pytest -v
-                pytest -v run_examples.py
-                cd ..  # N.B. This seems to happen automatically if omitted.
-                       # Less confusing to include it explicitly.
+                # cd tests
+                pytest -vvsx -k test_interpolatedimage_maxk_python
+                # coverage run -m pytest -v
+                # pytest -v run_examples.py
+                # cd ..  # N.B. This seems to happen automatically if omitted.
+                #        # Less confusing to include it explicitly.
 
             - name: Check shared library
               if: matrix.py == 3.12
diff --git a/tests/test_interpolatedimage.py b/tests/test_interpolatedimage.py
index 6d8c6dbbd4..c509bcfb3c 100644
--- a/tests/test_interpolatedimage.py
+++ b/tests/test_interpolatedimage.py
@@ -1890,5 +1890,87 @@ def test_drawreal_seg_fault():
     np.testing.assert_array_equal(image.array, 0)
 
 
+def _find_maxk(kim, thresh, count_thresh=True):
+    No2 = kim.xmax
+    dk = np.pi / No2
+
+    thresh *= thresh
+    max_ix = No2
+    n_below_thresh = 0
+    maxk_ix = 0
+
+    ix = 0
+    while ix <= max_ix:
+
+        iy = 0
+        while iy <= ix:
+            # The bottom side of the square in the lower-right quadrant.
+            norm_kval = kim(iy, -ix)
+            norm_kval = (norm_kval * norm_kval.conjugate()).real
+
+            if norm_kval <= thresh and iy != ix and ix != No2:
+                # The top side of the square in the upper-right quadrant.
+                norm_kval = kim(iy, ix)
+                norm_kval = (norm_kval * norm_kval.conjugate()).real
+
+            if norm_kval <= thresh and iy > 0:
+                # The right side of the square in the lower-right quadrant.
+                norm_kval = kim(ix, -iy)
+                norm_kval = (norm_kval * norm_kval.conjugate()).real
+
+            if norm_kval <= thresh and ix > 0 and iy != No2:
+                # The right side of the square in the upper-right quadrant.
+                norm_kval = kim(ix, iy)
+                norm_kval = (norm_kval * norm_kval.conjugate()).real
+
+            iy += 1
+
+            if norm_kval > thresh:
+                maxk_ix = ix
+                n_below_thresh = 0
+                break
+
+        if count_thresh:
+            if norm_kval <= thresh:
+                n_below_thresh += 1
+        else:
+            n_below_thresh += 1
+        if n_below_thresh == 5:
+            break
+
+        ix += 1
+
+    maxk_ix += 1
+    return maxk_ix * dk
+
+
+@timer
+def test_interpolatedimage_maxk_python():
+    # this code makes an image where there is a gap in the fourier
+    # space image of 4 pixels where pixels go above and below the
+    # maxk threshold. Four pixels is exactly the gap needed to trigger
+    # the bug in the maxk code we are testing for.
+    im = galsim.Gaussian(fwhm=0.9).drawImage(scale=0.2)
+    iim = galsim.InterpolatedImage(im)
+    orig_maxk = iim.maxk
+    kim = iim._xim.copy()
+    kim.wcs = galsim.PixelScale(1.0)
+    kim = kim.calculate_fft()
+    kx, ky = kim.get_pixel_centers()
+    kx *= kim.scale
+    ky *= kim.scale
+    thresh = iim.gsparams.maxk_threshold * kim(0,0).real
+    maxk_py = _find_maxk(kim, thresh, count_thresh=True)
+    maxk_ix = np.floor(maxk_py / kim.scale).astype(int)
+    kim[maxk_ix, maxk_ix + 4] = kim[0, 0].real * 0.1
+
+    maxk_py_false = _find_maxk(kim, thresh, count_thresh=False) / im.wcs._maxScale()
+    maxk_py_true = _find_maxk(kim, thresh, count_thresh=True) / im.wcs._maxScale()
+
+    print("galsim|pybuggy|pyfixed:", orig_maxk, maxk_py_false, maxk_py_true)
+    assert maxk_py_false != maxk_py_true
+    assert orig_maxk == maxk_py_true
+
+
 if __name__ == "__main__":
     runtests(__file__)

From 4837e3258c6547b0f98e4cb0a64b2502285d325c Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Thu, 30 Apr 2026 06:59:55 -0500
Subject: [PATCH 03/15] test: working test!

---
 tests/test_interpolatedimage.py | 32 ++++++++++++++++++++++----------
 1 file changed, 22 insertions(+), 10 deletions(-)

diff --git a/tests/test_interpolatedimage.py b/tests/test_interpolatedimage.py
index c509bcfb3c..daf69139b4 100644
--- a/tests/test_interpolatedimage.py
+++ b/tests/test_interpolatedimage.py
@@ -1923,13 +1923,13 @@ def _find_maxk(kim, thresh, count_thresh=True):
                 norm_kval = kim(ix, iy)
                 norm_kval = (norm_kval * norm_kval.conjugate()).real
 
-            iy += 1
-
             if norm_kval > thresh:
                 maxk_ix = ix
                 n_below_thresh = 0
                 break
 
+            iy += 1
+
         if count_thresh:
             if norm_kval <= thresh:
                 n_below_thresh += 1
@@ -1950,12 +1950,10 @@ def test_interpolatedimage_maxk_python():
     # space image of 4 pixels where pixels go above and below the
     # maxk threshold. Four pixels is exactly the gap needed to trigger
     # the bug in the maxk code we are testing for.
-    im = galsim.Gaussian(fwhm=0.9).drawImage(scale=0.2)
+    im = galsim.Gaussian(fwhm=0.9 / 0.2).drawImage(scale=1)
     iim = galsim.InterpolatedImage(im)
-    orig_maxk = iim.maxk
-    kim = iim._xim.copy()
-    kim.wcs = galsim.PixelScale(1.0)
-    kim = kim.calculate_fft()
+    xim = iim._xim.copy()
+    kim = xim.calculate_fft()
     kx, ky = kim.get_pixel_centers()
     kx *= kim.scale
     ky *= kim.scale
@@ -1964,12 +1962,26 @@ def test_interpolatedimage_maxk_python():
     maxk_ix = np.floor(maxk_py / kim.scale).astype(int)
     kim[maxk_ix, maxk_ix + 4] = kim[0, 0].real * 0.1
 
+    new_im = kim.calculate_inverse_fft()
+    sbii = galsim._galsim.SBInterpolatedImage(
+        new_im._image,
+        new_im.bounds._b,
+        iim._pad_image.bounds._b,
+        iim._x_interpolant._i,
+        iim._k_interpolant._i,
+        0,
+        0,
+        iim.gsparams._gsp,
+    )
+
     maxk_py_false = _find_maxk(kim, thresh, count_thresh=False) / im.wcs._maxScale()
     maxk_py_true = _find_maxk(kim, thresh, count_thresh=True) / im.wcs._maxScale()
+    sbii.calculateMaxK(0)
+    sbii_maxk = sbii.maxK()
 
-    print("galsim|pybuggy|pyfixed:", orig_maxk, maxk_py_false, maxk_py_true)
-    assert maxk_py_false != maxk_py_true
-    assert orig_maxk == maxk_py_true
+    print("galsim|pybuggy|pyfixed:", sbii_maxk, maxk_py_false, maxk_py_true)
+    assert np.abs(maxk_py_false - maxk_py_true) >= kim.scale
+    np.testing.assert_allclose(sbii_maxk, maxk_py_true, atol=1e-12, rtol=0)
 
 
 if __name__ == "__main__":

From f3d3c878e5a5e35719edd72997f96877611b574d Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Thu, 30 Apr 2026 07:58:49 -0500
Subject: [PATCH 04/15] test: better test of code

---
 tests/test_interpolatedimage.py | 133 +++++++++++---------------------
 1 file changed, 45 insertions(+), 88 deletions(-)

diff --git a/tests/test_interpolatedimage.py b/tests/test_interpolatedimage.py
index daf69139b4..ef0a1d8a3f 100644
--- a/tests/test_interpolatedimage.py
+++ b/tests/test_interpolatedimage.py
@@ -1890,98 +1890,55 @@ def test_drawreal_seg_fault():
     np.testing.assert_array_equal(image.array, 0)
 
 
-def _find_maxk(kim, thresh, count_thresh=True):
-    No2 = kim.xmax
-    dk = np.pi / No2
-
-    thresh *= thresh
-    max_ix = No2
-    n_below_thresh = 0
-    maxk_ix = 0
-
-    ix = 0
-    while ix <= max_ix:
-
-        iy = 0
-        while iy <= ix:
-            # The bottom side of the square in the lower-right quadrant.
-            norm_kval = kim(iy, -ix)
-            norm_kval = (norm_kval * norm_kval.conjugate()).real
-
-            if norm_kval <= thresh and iy != ix and ix != No2:
-                # The top side of the square in the upper-right quadrant.
-                norm_kval = kim(iy, ix)
-                norm_kval = (norm_kval * norm_kval.conjugate()).real
-
-            if norm_kval <= thresh and iy > 0:
-                # The right side of the square in the lower-right quadrant.
-                norm_kval = kim(ix, -iy)
-                norm_kval = (norm_kval * norm_kval.conjugate()).real
-
-            if norm_kval <= thresh and ix > 0 and iy != No2:
-                # The right side of the square in the upper-right quadrant.
-                norm_kval = kim(ix, iy)
-                norm_kval = (norm_kval * norm_kval.conjugate()).real
-
-            if norm_kval > thresh:
-                maxk_ix = ix
-                n_below_thresh = 0
-                break
-
-            iy += 1
-
-        if count_thresh:
-            if norm_kval <= thresh:
-                n_below_thresh += 1
-        else:
-            n_below_thresh += 1
-        if n_below_thresh == 5:
-            break
-
-        ix += 1
-
-    maxk_ix += 1
-    return maxk_ix * dk
-
-
 @timer
 def test_interpolatedimage_maxk_python():
     # this code makes an image where there is a gap in the fourier
-    # space image of 4 pixels where pixels go above and below the
-    # maxk threshold. Four pixels is exactly the gap needed to trigger
-    # the bug in the maxk code we are testing for.
-    im = galsim.Gaussian(fwhm=0.9 / 0.2).drawImage(scale=1)
-    iim = galsim.InterpolatedImage(im)
-    xim = iim._xim.copy()
-    kim = xim.calculate_fft()
-    kx, ky = kim.get_pixel_centers()
-    kx *= kim.scale
-    ky *= kim.scale
-    thresh = iim.gsparams.maxk_threshold * kim(0,0).real
-    maxk_py = _find_maxk(kim, thresh, count_thresh=True)
-    maxk_ix = np.floor(maxk_py / kim.scale).astype(int)
-    kim[maxk_ix, maxk_ix + 4] = kim[0, 0].real * 0.1
-
-    new_im = kim.calculate_inverse_fft()
-    sbii = galsim._galsim.SBInterpolatedImage(
-        new_im._image,
-        new_im.bounds._b,
-        iim._pad_image.bounds._b,
-        iim._x_interpolant._i,
-        iim._k_interpolant._i,
-        0,
-        0,
-        iim.gsparams._gsp,
-    )
-
-    maxk_py_false = _find_maxk(kim, thresh, count_thresh=False) / im.wcs._maxScale()
-    maxk_py_true = _find_maxk(kim, thresh, count_thresh=True) / im.wcs._maxScale()
-    sbii.calculateMaxK(0)
-    sbii_maxk = sbii.maxK()
+    # 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.
+
+    def _compute_maxk_cpp(xim, iim):
+        # this little function exists only to invoke the C++
+        # maxk code...
+        # we use copies to avoid side effects
+        ikim = xim.copy()
+        sbii = galsim._galsim.SBInterpolatedImage(
+            ikim._image,
+            ikim.bounds._b,
+            iim._pad_image.copy().bounds._b,
+            iim._x_interpolant._i,
+            iim._k_interpolant._i,
+            0,
+            0,
+            iim.gsparams._gsp,
+        )
+        sbii.calculateMaxK(0)  # this call is needed to invoke the C++ code
+        return sbii.maxK()
 
-    print("galsim|pybuggy|pyfixed:", sbii_maxk, maxk_py_false, maxk_py_true)
-    assert np.abs(maxk_py_false - maxk_py_true) >= kim.scale
-    np.testing.assert_allclose(sbii_maxk, maxk_py_true, atol=1e-12, rtol=0)
+    im = galsim.Gaussian(fwhm=0.9 / 0.2).drawImage(scale=1)
+    iim = galsim.InterpolatedImage(im, scale=1)
+    orig_maxk = _compute_maxk_cpp(iim._xim, iim)
+
+    for offset in [3, 4, 5, 6, 7]:
+        kim = iim._xim.copy().calculate_fft()
+        kx, ky = kim.get_pixel_centers()
+        kx *= kim.scale
+        ky *= kim.scale
+        # maxk_ix is the last pixel above threshold
+        # we subtract 1 since galsim adds 1 in the C++ code
+        maxk_ix = np.floor(orig_maxk / kim.scale).astype(int) - 1
+        kim[maxk_ix, maxk_ix + offset] = kim[0, 0].real * 0.1
+        new_im = kim.calculate_inverse_fft()
+        new_maxk = _compute_maxk_cpp(new_im, iim)
+
+        print("offset|orig|new:", offset, orig_maxk, new_maxk)
+
+        if offset < 5:
+            # for offsets smaller than 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)
 
 
 if __name__ == "__main__":

From e52a9bcc7d233d62ecf81307d4b424d79328cf0a Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Thu, 30 Apr 2026 08:20:15 -0500
Subject: [PATCH 05/15] test: better test for five pixel gaps

---
 tests/test_interpolatedimage.py | 20 +++++++++++---------
 1 file changed, 11 insertions(+), 9 deletions(-)

diff --git a/tests/test_interpolatedimage.py b/tests/test_interpolatedimage.py
index ef0a1d8a3f..e65e24f14b 100644
--- a/tests/test_interpolatedimage.py
+++ b/tests/test_interpolatedimage.py
@@ -1915,26 +1915,28 @@ def _compute_maxk_cpp(xim, iim):
         sbii.calculateMaxK(0)  # this call is needed to invoke the C++ code
         return sbii.maxK()
 
-    im = galsim.Gaussian(fwhm=0.9 / 0.2).drawImage(scale=1)
-    iim = galsim.InterpolatedImage(im, scale=1)
-    orig_maxk = _compute_maxk_cpp(iim._xim, iim)
-
+    print(" ")
     for offset in [3, 4, 5, 6, 7]:
+        im = galsim.Gaussian(fwhm=0.9 / 0.2).drawImage(scale=1)
+        iim = galsim.InterpolatedImage(im, scale=1)
+        orig_maxk = _compute_maxk_cpp(iim._xim, iim)
+
         kim = iim._xim.copy().calculate_fft()
         kx, ky = kim.get_pixel_centers()
         kx *= kim.scale
         ky *= kim.scale
-        # maxk_ix is the last pixel above threshold
-        # we subtract 1 since galsim adds 1 in the C++ code
+        # 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
-        kim[maxk_ix, maxk_ix + offset] = kim[0, 0].real * 0.1
+        kim[maxk_ix, maxk_ix + offset] = kim[0, 0].real
         new_im = kim.calculate_inverse_fft()
         new_maxk = _compute_maxk_cpp(new_im, iim)
 
         print("offset|orig|new:", offset, orig_maxk, new_maxk)
 
-        if offset < 5:
-            # for offsets smaller than 5, we should get an offset of offset pixels
+        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:

From b254bbf8b7047c5e15230e18150f7ee7817f33e1 Mon Sep 17 00:00:00 2001
From: "Matthew R. Becker" <beckermr@users.noreply.github.com>
Date: Thu, 30 Apr 2026 07:03:13 -0500
Subject: [PATCH 06/15] Apply suggestion from @beckermr

---
 .github/workflows/ci.yml | 11 +++++------
 1 file changed, 5 insertions(+), 6 deletions(-)

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 6d27de48a8..6fda7ddaaf 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -201,12 +201,11 @@ jobs:
 
             - name: Run unit tests
               run: |
-                # cd tests
-                pytest -vvsx -k test_interpolatedimage_maxk_python
-                # coverage run -m pytest -v
-                # pytest -v run_examples.py
-                # cd ..  # N.B. This seems to happen automatically if omitted.
-                #        # Less confusing to include it explicitly.
+                cd tests
+                coverage run -m pytest -v
+                pytest -v run_examples.py
+                cd ..  # N.B. This seems to happen automatically if omitted.
+                       # Less confusing to include it explicitly.
 
             - name: Check shared library
               if: matrix.py == 3.12

From eb07d6d4a22efe8717e64fbd353b32b6a0fbe267 Mon Sep 17 00:00:00 2001
From: "Matthew R. Becker" <beckermr@users.noreply.github.com>
Date: Thu, 30 Apr 2026 08:23:03 -0500
Subject: [PATCH 07/15] Apply suggestion from @beckermr

---
 tests/test_interpolatedimage.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_interpolatedimage.py b/tests/test_interpolatedimage.py
index e65e24f14b..9442fbf407 100644
--- a/tests/test_interpolatedimage.py
+++ b/tests/test_interpolatedimage.py
@@ -1891,7 +1891,7 @@ def test_drawreal_seg_fault():
 
 
 @timer
-def test_interpolatedimage_maxk_python():
+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

From 47f8eea20036e42f456e7975e3e7c42c3af0dfac Mon Sep 17 00:00:00 2001
From: "Matthew R. Becker" <beckermr@users.noreply.github.com>
Date: Thu, 30 Apr 2026 08:23:38 -0500
Subject: [PATCH 08/15] Apply suggestion from @beckermr

---
 tests/test_interpolatedimage.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_interpolatedimage.py b/tests/test_interpolatedimage.py
index 9442fbf407..53cbfe6a51 100644
--- a/tests/test_interpolatedimage.py
+++ b/tests/test_interpolatedimage.py
@@ -1894,7 +1894,7 @@ def test_drawreal_seg_fault():
 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
+    # and below the maxk threshold. At >five pixels, galsim should
     # ignore the gap, but less than that it should increase maxk.
 
     def _compute_maxk_cpp(xim, iim):

From 27ac5df9f27505608ea36fc164f856ce7774b339 Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Thu, 30 Apr 2026 13:41:47 -0500
Subject: [PATCH 09/15] test: add code to print table

---
 tests/test_interpolatedimage.py | 14 +++++++++++---
 1 file changed, 11 insertions(+), 3 deletions(-)

diff --git a/tests/test_interpolatedimage.py b/tests/test_interpolatedimage.py
index 53cbfe6a51..8ec80f5bef 100644
--- a/tests/test_interpolatedimage.py
+++ b/tests/test_interpolatedimage.py
@@ -1915,11 +1915,12 @@ def _compute_maxk_cpp(xim, iim):
         sbii.calculateMaxK(0)  # this call is needed to invoke the C++ code
         return sbii.maxK()
 
-    print(" ")
+    print("\n| offset | orig       | new via direct C++ | new via galsim II |")
+    print("|--------|------------|--------------------|-------------------|")
     for offset in [3, 4, 5, 6, 7]:
         im = galsim.Gaussian(fwhm=0.9 / 0.2).drawImage(scale=1)
         iim = galsim.InterpolatedImage(im, scale=1)
-        orig_maxk = _compute_maxk_cpp(iim._xim, iim)
+        orig_maxk = iim.maxk
 
         kim = iim._xim.copy().calculate_fft()
         kx, ky = kim.get_pixel_centers()
@@ -1933,7 +1934,12 @@ def _compute_maxk_cpp(xim, iim):
         new_im = kim.calculate_inverse_fft()
         new_maxk = _compute_maxk_cpp(new_im, iim)
 
-        print("offset|orig|new:", offset, orig_maxk, new_maxk)
+        print("| % 6d | %10.6f | %18.6f | %17.6f |" % (
+            offset,
+            orig_maxk,
+            new_maxk,
+            galsim.InterpolatedImage(new_im, scale=1).maxk),
+        )
 
         if offset <= 5:
             # for offsets <=5, we should get an offset of offset pixels
@@ -1942,6 +1948,8 @@ def _compute_maxk_cpp(xim, iim):
         else:
             np.testing.assert_allclose(new_maxk, orig_maxk, atol=1e-12, rtol=0)
 
+    print("|--------|------------|--------------------|-------------------|")
+
 
 if __name__ == "__main__":
     runtests(__file__)

From 60cd2dc9dc99c016680fc7e0c093514c6500da34 Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Thu, 30 Apr 2026 13:57:24 -0500
Subject: [PATCH 10/15] test: add some code showing how the FFTs relate

---
 tests/test_interpolatedimage.py | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/tests/test_interpolatedimage.py b/tests/test_interpolatedimage.py
index 8ec80f5bef..1ab41dea48 100644
--- a/tests/test_interpolatedimage.py
+++ b/tests/test_interpolatedimage.py
@@ -1917,8 +1917,10 @@ def _compute_maxk_cpp(xim, iim):
 
     print("\n| offset | orig       | new via direct C++ | new via galsim II |")
     print("|--------|------------|--------------------|-------------------|")
-    for offset in [3, 4, 5, 6, 7]:
+    for offset in [0, 3, 4, 5, 6, 7]:
         im = galsim.Gaussian(fwhm=0.9 / 0.2).drawImage(scale=1)
+        new_im = im.copy().calculate_fft().calculate_inverse_fft()
+        np.testing.assert_allclose(im.array, new_im[im.bounds].array, atol=1e-6, rtol=1e-6)
         iim = galsim.InterpolatedImage(im, scale=1)
         orig_maxk = iim.maxk
 
@@ -1930,15 +1932,19 @@ def _compute_maxk_cpp(xim, iim):
         # 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
-        kim[maxk_ix, maxk_ix + offset] = kim[0, 0].real
+        if offset > 0:
+            kim[maxk_ix, maxk_ix + offset] = kim[0, 0].real * iim.gsparams.maxk_threshold * 2.0
         new_im = kim.calculate_inverse_fft()
         new_maxk = _compute_maxk_cpp(new_im, iim)
 
+        if offset == 0:
+            np.testing.assert_allclose(im.array, new_im[iim._image.bounds].array, atol=1e-6, rtol=1e-6)
+
         print("| % 6d | %10.6f | %18.6f | %17.6f |" % (
             offset,
             orig_maxk,
             new_maxk,
-            galsim.InterpolatedImage(new_im, scale=1).maxk),
+            galsim.InterpolatedImage(new_im[iim._image.bounds], scale=1).maxk),
         )
 
         if offset <= 5:

From 1f9ad23b649fe500fb55617c9126ae7150435dc8 Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Thu, 30 Apr 2026 14:06:03 -0500
Subject: [PATCH 11/15] test: put back what generates weird results

---
 tests/test_interpolatedimage.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_interpolatedimage.py b/tests/test_interpolatedimage.py
index 1ab41dea48..10ef93eed3 100644
--- a/tests/test_interpolatedimage.py
+++ b/tests/test_interpolatedimage.py
@@ -1933,7 +1933,7 @@ def _compute_maxk_cpp(xim, iim):
         # and so we subtract 1
         maxk_ix = np.floor(orig_maxk / kim.scale).astype(int) - 1
         if offset > 0:
-            kim[maxk_ix, maxk_ix + offset] = kim[0, 0].real * iim.gsparams.maxk_threshold * 2.0
+            kim[maxk_ix, maxk_ix + offset] = kim[0, 0].real
         new_im = kim.calculate_inverse_fft()
         new_maxk = _compute_maxk_cpp(new_im, iim)
 

From 3be432260ccce42ff6cecedfff52ecccbf835b41 Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Thu, 30 Apr 2026 14:12:02 -0500
Subject: [PATCH 12/15] test: cleanup bounds handling, weirdness remains

---
 tests/test_interpolatedimage.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/tests/test_interpolatedimage.py b/tests/test_interpolatedimage.py
index 10ef93eed3..e54b48f331 100644
--- a/tests/test_interpolatedimage.py
+++ b/tests/test_interpolatedimage.py
@@ -1936,15 +1936,16 @@ def _compute_maxk_cpp(xim, iim):
             kim[maxk_ix, maxk_ix + offset] = kim[0, 0].real
         new_im = kim.calculate_inverse_fft()
         new_maxk = _compute_maxk_cpp(new_im, iim)
+        new_im = new_im[iim._image.bounds]
 
         if offset == 0:
-            np.testing.assert_allclose(im.array, new_im[iim._image.bounds].array, atol=1e-6, rtol=1e-6)
+            np.testing.assert_allclose(im.array, new_im.array, atol=1e-6, rtol=1e-6)
 
         print("| % 6d | %10.6f | %18.6f | %17.6f |" % (
             offset,
             orig_maxk,
             new_maxk,
-            galsim.InterpolatedImage(new_im[iim._image.bounds], scale=1).maxk),
+            galsim.InterpolatedImage(new_im, scale=1).maxk),
         )
 
         if offset <= 5:

From 69afde08b9338850ea4a4b9eebb0264a008e0159 Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Fri, 1 May 2026 22:32:51 -0500
Subject: [PATCH 13/15] fix: use simpler test code

---
 tests/test_interpolatedimage.py | 40 +++++++--------------------------
 1 file changed, 8 insertions(+), 32 deletions(-)

diff --git a/tests/test_interpolatedimage.py b/tests/test_interpolatedimage.py
index e54b48f331..34f24ceba1 100644
--- a/tests/test_interpolatedimage.py
+++ b/tests/test_interpolatedimage.py
@@ -1897,30 +1897,10 @@ def test_interpolatedimage_maxk_kspace_pixel_gap():
     # and below the maxk threshold. At >five pixels, galsim should
     # ignore the gap, but less than that it should increase maxk.
 
-    def _compute_maxk_cpp(xim, iim):
-        # this little function exists only to invoke the C++
-        # maxk code...
-        # we use copies to avoid side effects
-        ikim = xim.copy()
-        sbii = galsim._galsim.SBInterpolatedImage(
-            ikim._image,
-            ikim.bounds._b,
-            iim._pad_image.copy().bounds._b,
-            iim._x_interpolant._i,
-            iim._k_interpolant._i,
-            0,
-            0,
-            iim.gsparams._gsp,
-        )
-        sbii.calculateMaxK(0)  # this call is needed to invoke the C++ code
-        return sbii.maxK()
-
-    print("\n| offset | orig       | new via direct C++ | new via galsim II |")
-    print("|--------|------------|--------------------|-------------------|")
+    print("\n| offset | orig       | new                |")
+    print("|--------|------------|--------------------|")
     for offset in [0, 3, 4, 5, 6, 7]:
         im = galsim.Gaussian(fwhm=0.9 / 0.2).drawImage(scale=1)
-        new_im = im.copy().calculate_fft().calculate_inverse_fft()
-        np.testing.assert_allclose(im.array, new_im[im.bounds].array, atol=1e-6, rtol=1e-6)
         iim = galsim.InterpolatedImage(im, scale=1)
         orig_maxk = iim.maxk
 
@@ -1933,19 +1913,15 @@ def _compute_maxk_cpp(xim, iim):
         # and so we subtract 1
         maxk_ix = np.floor(orig_maxk / kim.scale).astype(int) - 1
         if offset > 0:
-            kim[maxk_ix, maxk_ix + offset] = kim[0, 0].real
+            val = kim[maxk_ix + offset, maxk_ix]
+            kim[maxk_ix + offset, maxk_ix] = val / np.abs(val) * kim[0, 0].real
         new_im = kim.calculate_inverse_fft()
-        new_maxk = _compute_maxk_cpp(new_im, iim)
-        new_im = new_im[iim._image.bounds]
-
-        if offset == 0:
-            np.testing.assert_allclose(im.array, new_im.array, atol=1e-6, rtol=1e-6)
+        new_maxk = galsim.InterpolatedImage(new_im, scale=1, pad_factor=1).maxk
 
-        print("| % 6d | %10.6f | %18.6f | %17.6f |" % (
+        print("| % 6d | %10.6f | %18.6f |" % (
             offset,
             orig_maxk,
-            new_maxk,
-            galsim.InterpolatedImage(new_im, scale=1).maxk),
+            new_maxk),
         )
 
         if offset <= 5:
@@ -1955,7 +1931,7 @@ def _compute_maxk_cpp(xim, iim):
         else:
             np.testing.assert_allclose(new_maxk, orig_maxk, atol=1e-12, rtol=0)
 
-    print("|--------|------------|--------------------|-------------------|")
+    print("|--------|------------|--------------------|")
 
 
 if __name__ == "__main__":

From f36ce82b5bf4b53d711b9ddde172c8768a406fcf Mon Sep 17 00:00:00 2001
From: "Matthew R. Becker" <beckermr@users.noreply.github.com>
Date: Mon, 4 May 2026 08:33:08 -0500
Subject: [PATCH 14/15] Update test_interpolatedimage.py

Co-authored-by: Mike Jarvis <michael@jarvis.net>
---
 tests/test_interpolatedimage.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_interpolatedimage.py b/tests/test_interpolatedimage.py
index 34f24ceba1..8c0a845685 100644
--- a/tests/test_interpolatedimage.py
+++ b/tests/test_interpolatedimage.py
@@ -1900,7 +1900,7 @@ def test_interpolatedimage_maxk_kspace_pixel_gap():
     print("\n| offset | orig       | new                |")
     print("|--------|------------|--------------------|")
     for offset in [0, 3, 4, 5, 6, 7]:
-        im = galsim.Gaussian(fwhm=0.9 / 0.2).drawImage(scale=1)
+        im = galsim.Gaussian(fwhm=4.5).drawImage(scale=1)
         iim = galsim.InterpolatedImage(im, scale=1)
         orig_maxk = iim.maxk
 

From 3ebca0ed2af83c10a0212b922b403b8d4446fc88 Mon Sep 17 00:00:00 2001
From: "Matthew R. Becker" <beckermr@users.noreply.github.com>
Date: Mon, 4 May 2026 08:34:57 -0500
Subject: [PATCH 15/15] Update test_interpolatedimage.py

---
 tests/test_interpolatedimage.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/tests/test_interpolatedimage.py b/tests/test_interpolatedimage.py
index 8c0a845685..61dd925e04 100644
--- a/tests/test_interpolatedimage.py
+++ b/tests/test_interpolatedimage.py
@@ -1913,8 +1913,7 @@ def test_interpolatedimage_maxk_kspace_pixel_gap():
         # and so we subtract 1
         maxk_ix = np.floor(orig_maxk / kim.scale).astype(int) - 1
         if offset > 0:
-            val = kim[maxk_ix + offset, maxk_ix]
-            kim[maxk_ix + offset, maxk_ix] = val / np.abs(val) * kim[0, 0].real
+            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