Fix incorrect implementation of bilinear interpolation (#240)

* Fix incorrect implementation of bilinear interpolation. Both the pixel index calculation and the weights calculation were wrong

* Added a unit test to ensure bilinear mapmaking stays correct

Author: amaurea

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<NonTiled, Interpol> {
                     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<NonTiled, NearestNeighbor>::GetPixels(int i_det, int
 template<>
 inline int Pixelizor2_Flat<NonTiled, Bilinear>::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<Tiled, NearestNeighbor>::GetPixels(int i_det, int i_t
 template<>
 inline int Pixelizor2_Flat<Tiled, Bilinear>::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++) {

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