From 49c36da159b1b1a7ba29d0f75876ba016403ba44 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sigurd=20Kirkevold=20N=C3=A6ss?= Date: Fri, 17 Apr 2026 11:34:51 -0700 Subject: [PATCH 1/2] Fix incorrect implementation of bilinear interpolation. Both the pixel index calculation and the weights calculation were wrong --- src/Projection.cxx | 31 +++++++++++++++++++------------ 1 file changed, 19 insertions(+), 12 deletions(-) diff --git a/src/Projection.cxx b/src/Projection.cxx index aabf5599..fd1b2298 100644 --- a/src/Projection.cxx +++ b/src/Projection.cxx @@ -51,6 +51,8 @@ inline bool isNone(const bp::object &pyo) return (pyo.ptr() == Py_None); } +int ifloor(double x) { return int(x)-int(x<0); } +int iround(double x) { return ifloor(x+0.5); } // ProjEng template system // @@ -812,6 +814,7 @@ class Pixelizor2_Flat { double iy0=0., double ix0=0.) { naxis[0] = ny; naxis[1] = nx; + // Note, this y,x order is the opposite of what wcslib uses cdelt[0] = dy; cdelt[1] = dx; crpix[0] = iy0; @@ -930,12 +933,16 @@ inline int Pixelizor2_Flat::GetPixels(int i_det, int template<> inline int Pixelizor2_Flat::GetPixels(int i_det, int i_time, const double *coords, int pixinds[interp_count][index_count], FSIGNAL pixweights[interp_count]) { // For bilinear mapmaking we need to visit the four bounding pixels - double x = coords[0] / cdelt[1] + crpix[1] - 1 + 0.5; - double y = coords[1] / cdelt[0] + crpix[0] - 1 + 0.5; - int x1 = int(x)-int(x<0); - int y1 = int(y)-int(y<0); - double wx[2] = {x-x1, 1-(x-x1)}; - double wy[2] = {y-y1, 1-(y-y1)}; + // 0-based pixel coordinate + double x = coords[0] / cdelt[1] + crpix[1] - 1; + double y = coords[1] / cdelt[0] + crpix[0] - 1; + // index of pixel to the left of this point. The pixel to the right of it + // will be that number +1 + int x1 = ifloor(x); + int y1 = ifloor(y); + // Weight of before and after pixels. Sum to 1. + double wx[2] = {1-(x-x1), x-x1}; + double wy[2] = {1-(y-y1), y-y1}; // Loop through the our cases int iout = 0; for(int iy = y1; iy < y1+2; iy++) { @@ -1141,12 +1148,12 @@ inline int Pixelizor2_Flat::GetPixels(int i_det, int i_t template<> inline int Pixelizor2_Flat::GetPixels(int i_det, int i_time, const double *coords, int pixinds[interp_count][index_count], FSIGNAL pixweights[interp_count]) { // For bilinear mapmaking we need to visit the four bounding pixels - double x = coords[0] / parent_pix.cdelt[1] + parent_pix.crpix[1] - 1 + 0.5; - double y = coords[1] / parent_pix.cdelt[0] + parent_pix.crpix[0] - 1 + 0.5; - int x1 = int(x); - int y1 = int(y); - double wx[2] = {x-x1, 1-(x-x1)}; - double wy[2] = {y-y1, 1-(y-y1)}; + double x = coords[0] / parent_pix.cdelt[1] + parent_pix.crpix[1] - 1; + double y = coords[1] / parent_pix.cdelt[0] + parent_pix.crpix[0] - 1; + int x1 = ifloor(x); + int y1 = ifloor(y); + double wx[2] = {1-(x-x1), x-x1}; + double wy[2] = {1-(y-y1), y-y1}; // Loop through the our cases int iout = 0; for(int iy = y1; iy < y1+2; iy++) { From 377901bd593349982bc4fa5275b916f8b71417b6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sigurd=20Kirkevold=20N=C3=A6ss?= Date: Mon, 20 Apr 2026 11:29:19 -0700 Subject: [PATCH 2/2] Added a unit test to ensure bilinear mapmaking stays correct --- test/test_proj_eng.py | 46 ++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 45 insertions(+), 1 deletion(-) diff --git a/test/test_proj_eng.py b/test/test_proj_eng.py index a3eb1cf7..5a2a7b2e 100644 --- a/test/test_proj_eng.py +++ b/test/test_proj_eng.py @@ -7,7 +7,7 @@ # Don't require pixell for testing try: - from pixell import enmap + from pixell import enmap, wcsutils pixell_found = True except ModuleNotFoundError: pixell_found = False @@ -150,6 +150,50 @@ def test_20_threads(self): if interpol == 'nearest': self.assertEqual(counts1.sum(), 0) + @requires_pixell + def test_30_bilin(self): + """1-sample single boresight-pointed detector with no coordinate transformation""" + # Trivial geometry where pixel coordinate and sky coordinate are the same thing + shape = (2,2) + wcs = wcsutils.explicit(crval=[0,0], crpix=[1,1], cdelt=[1,1], ctype=["RA---CAR","DEC--CAR"]) + # Cases we will consider + cases = [ + # A single sample hitting (0.1, 0.7) + np.array([[0.1],[0.7]]), + # A single sample with integer coordinates + np.array([[1.0],[0.0]]), + ] + for x, y in cases: + whole = np.all(x==np.round(x)) and np.all(y==np.round(y)) + csl = proj.CelestialSightLine.for_lonlat(x*DEG, y*DEG) + fp = proj.FocalPlane.from_xieta([0.0],[0.0]) + asm = proj.Assembly.attach(csl, fp) + dtype = np.float32 + tod = np.ones((1,1), dtype) + # Expected result for NN and bilin + wy = np.sum(np.round([1-y,y]),1) + wx = np.sum(np.round([1-x,x]),1) + targ_nn = wy[:,None]*wx[None,:] + wy = np.sum([1-y,y],1) + wx = np.sum([1-x,x],1) + targ_li = wy[:,None]*wx[None,:] + # Nearest neighbor untiled + p = proj.Projectionist.for_geom(shape, wcs, interpol="nearest") + m_nn = p.to_map(tod, asm, comps="T")[0] + assert np.allclose(m_nn, targ_nn) + # Bilinear untiled + p = proj.Projectionist.for_geom(shape, wcs, interpol="bilinear") + m_li = p.to_map(tod, asm, comps="T")[0] + assert np.allclose(m_li, targ_li) + if whole: assert np.allclose(m_nn, m_li) + # Nearest neighbor tiled + p = proj.Projectionist.for_tiled(shape, wcs, (10, 10), interpol="nearest") + m_nn = p.to_map(tod, asm, comps="T")[0][0] + assert np.allclose(m_nn, targ_nn) + p = proj.Projectionist.for_tiled(shape, wcs, (10, 10), interpol="bilinear") + m_li = p.to_map(tod, asm, comps="T")[0][0] + assert np.allclose(m_li, targ_li) + if whole: assert np.allclose(m_nn, m_li) if __name__ == '__main__': unittest.main()