Use padded singular values for building J matrix

The matrix J was being built using the un-padded values, whereas the
padded values ought to be used (as called for by DeForest (2004)). This
caused incorrect averaging weights to be computed for some sample
values, in transformations where the padding of singular values is
required.

For the implemented test cases, the effect of this fix is quite small.

reproject/adaptive/deforest.pyx

@@ -281,8 +281,6 @@ def map_coordinates(double[:,:] source, double[:,:] target, Ci, int max_samples_
     cdef double[:,:] J = np.zeros((2, 2))
     cdef double[:,:] U = np.zeros((2, 2))
     cdef double[:] s = np.zeros((2,))
-    cdef double[:] s_padded = np.zeros((2,))
-    cdef double[:] si = np.zeros((2,))
     cdef double[:,:] V = np.zeros((2, 2))
     cdef int samples_width
     cdef double[:] transformed = np.zeros((2,))
@@ -322,12 +320,14 @@ def map_coordinates(double[:,:] source, double[:,:] target, Ci, int max_samples_
 
                 # Find and pad the singular values of the Jacobian.
                 svd2x2_decompose(Ji, U, s, V)
-                s_padded[0] = max(1.0, s[0])
-                s_padded[1] = max(1.0, s[1])
-                si[0] = 1.0/s[0]
-                si[1] = 1.0/s[1]
-                svd2x2_compose(V, si, U, J)
-                svd2x2_compose(U, s_padded, V, Ji_padded)
+                s[0] = max(1.0, s[0])
+                s[1] = max(1.0, s[1])
+                svd2x2_compose(U, s, V, Ji_padded)
+                # Build J, the inverse of Ji, by using 1/s and swapping the
+                # order of U and V.
+                s[0] = 1.0/s[0]
+                s[1] = 1.0/s[1]
+                svd2x2_compose(V, s, U, J)
 
                 # We'll need to sample some number of input images to set this
                 # output pixel. Later on, we'll compute weights to assign to

reproject/adaptive/tests/test_core.py

@@ -103,7 +103,7 @@ def test_reproject_adaptive_2d_rotated(center_jacobian, roundtrip_coords):
 
 
 @pytest.mark.parametrize('roundtrip_coords', (False, True))
-def test_reproject_adaptive_high_aliasing_potential(roundtrip_coords):
+def test_reproject_adaptive_high_aliasing_potential_rotation(roundtrip_coords):
     # Generate sample data with vertical stripes alternating with every column
     data_in = np.arange(40*40).reshape((40, 40))
     data_in = (data_in) % 2
@@ -177,6 +177,36 @@ def test_reproject_adaptive_high_aliasing_potential(roundtrip_coords):
     np.testing.assert_allclose(array_out, data_ref)
 
 
+@pytest.mark.parametrize('roundtrip_coords', (False, True))
+def test_reproject_adaptive_high_aliasing_potential_shearing(roundtrip_coords):
+    # Generate sample data with vertical stripes alternating with every column
+    data_in = np.arange(40*40).reshape((40, 40))
+    data_in = (data_in) % 2
+
+    # Set up the input image coordinates, defining pixel coordinates as world
+    # coordinates (with an offset)
+    wcs_in = WCS(naxis=2)
+    wcs_in.wcs.crpix = 21, 21
+    wcs_in.wcs.crval = 0, 0
+    wcs_in.wcs.cdelt = 1, 1
+
+    # Set up the output image coordinates, with shearing in both x and y
+    wcs_out = WCS(naxis=2)
+    wcs_out.wcs.crpix = 3, 5
+    wcs_out.wcs.crval = 0, 0
+    wcs_out.wcs.cdelt = 1, 1
+    wcs_out.wcs.pc = np.array([[1, 1], [0, 1]]) @ np.array([[1, 0], [0.5, 1]])
+
+    array_out = reproject_adaptive((data_in, wcs_in),
+                                   wcs_out, shape_out=(11, 6),
+                                   return_footprint=False,
+                                   roundtrip_coords=False)
+
+    # We should get values close to 0.5 (and average of the 1s and 0s in the
+    # input image). This is as opposed to values near 0 or 1, which would
+    # indicate incorrect averaging of sampled points.
+    np.testing.assert_allclose(array_out, 0.5, atol=0.1, rtol=0)
+
 def prepare_test_data(file_format):
     pytest.importorskip('sunpy', minversion='2.1.0')
     from sunpy.map import Map