Add direct correlation correction unit test.

 * Add sourceless test image generation modification to `makeFakeImages`.
 * Add direct correlation estimation function `estimatePixelCorrelation`.
 * Add new test case `testNoiseDiffimCorrection`.

tests/test_imageDecorrelation.py

@@ -122,17 +122,18 @@ def makeFakeImages(size=(256, 256), svar=0.04, tvar=0.04, psf1=3.3, psf2=2.2, of
     im1 = np.random.normal(scale=np.sqrt(svar), size=x0im.shape)  # variance of science image
     im2 = np.random.normal(scale=np.sqrt(tvar), size=x0im.shape)  # variance of template
 
-    fluxes = np.random.uniform(50, 30000, n_sources)
-    xposns = np.random.uniform(xim.min()+16, xim.max()-5, n_sources)
-    yposns = np.random.uniform(yim.min()+16, yim.max()-5, n_sources)
+    if n_sources > 0:
+        fluxes = np.random.uniform(50, 30000, n_sources)
+        xposns = np.random.uniform(xim.min()+16, xim.max()-5, n_sources)
+        yposns = np.random.uniform(yim.min()+16, yim.max()-5, n_sources)
 
-    # Make the source closest to the center of the image the one that increases in flux
-    ind = np.argmin(xposns**2. + yposns**2.)
+        # Make the source closest to the center of the image the one that increases in flux
+        ind = np.argmin(xposns**2. + yposns**2.)
 
-    # vary the y-width of psf across x-axis of science image (zero means no variation):
-    psf1_yvary = psf_yvary_factor * (yim.mean() - yposns) / yim.max()
-    if verbose:
-        print('PSF y spatial-variation:', psf1_yvary.min(), psf1_yvary.max())
+        # vary the y-width of psf across x-axis of science image (zero means no variation):
+        psf1_yvary = psf_yvary_factor * (yim.mean() - yposns) / yim.max()
+        if verbose:
+            print('PSF y spatial-variation:', psf1_yvary.min(), psf1_yvary.max())
 
     for i in range(n_sources):
         flux = fluxes[i]
@@ -209,8 +210,48 @@ def makeExposure(imgArray, psfArray, imgVariance):
     return im1ex, im2ex
 
 
+def estimatePixelCorrelation(B, nDist=40, convEdge=17):
+    """Estimate correlation as a function of pixel distance in the image
+    by sampling pixel pairs.
+
+    Parameters
+    ----------
+    B : `numpy.ndarray` of N x N `float` elements
+        Noise only image with zero pixel expectation value and identical variance
+        in all pixels. Must have equal dimensions.
+    nDist : `int`, optional
+        Estimated distances goes from 0 to nDist-1.
+        nDist must be smaller than the half dimensions of B.
+    convEdge : `int`, optional
+        Edge width where convolution did not happen.
+
+    Returns
+    -------
+    S : `numpy.ndarray` of nDist `float` elements
+        Correlation from 0 to nDist-1 pix distance. Pixels are normed by their
+        variance estimation. S[0], the autocorrelation, should be close to 1.
+    """
+    S = np.zeros(nDist, dtype=float)
+    nSample = 10000
+    # Cannot use nDist wide edge, otherwise 2nd pixel can go off the image.
+    # Don't bother with it.
+    A = B/np.sqrt(np.mean(B[convEdge:-convEdge, convEdge:-convEdge] *
+                          B[convEdge:-convEdge, convEdge:-convEdge]))
+    lEdge = nDist + convEdge
+    rEdge = B.shape[0] - lEdge
+    for r in range(nDist):
+        ind1 = np.random.randint(lEdge, rEdge, (2, nSample))
+        ind2 = np.copy(ind1)
+        # generate delta x,y in random directions uniformly
+        c_dxy = np.exp(2.j*np.pi*np.random.random(nSample))
+        ind2[0] += np.around(np.real(c_dxy)*r).astype(int)
+        ind2[1] += np.around(np.imag(c_dxy)*r).astype(int)
+        S[r] = np.sum(A[ind1[0], ind1[1]] * A[ind2[0], ind2[1]])/nSample
+    return S
+
+
 class DiffimCorrectionTest(lsst.utils.tests.TestCase):
-    """!A test case for the diffim image decorrelation algorithm.
+    """A test case for the diffim image decorrelation algorithm.
     """
 
     def setUp(self):
@@ -226,7 +267,7 @@ def setUp(self):
                                                        "NO_DATA", "DETECTED_NEGATIVE"]))
 
     def _setUpImages(self, svar=0.04, tvar=0.04, varyPsf=0.):
-        """!Generate a fake aligned template and science image.
+        """Generate a fake aligned template and science image.
         """
 
         self.svar = svar  # variance of noise in science image
@@ -236,6 +277,17 @@ def _setUpImages(self, svar=0.04, tvar=0.04, varyPsf=0.):
             = makeFakeImages(svar=self.svar, tvar=self.tvar, psf1=self.psf1_sigma, psf2=self.psf2_sigma,
                              n_sources=50, psf_yvary_factor=varyPsf, verbose=False)
 
+    def _setUpSourcelessImages(self, svar, tvar):
+        """Generate noise only template and science images.
+        """
+
+        self.svar = svar  # variance of noise in science image
+        self.tvar = tvar  # variance of noise in template image
+
+        self.im1ex, self.im2ex = makeFakeImages(
+            svar=self.svar, tvar=self.tvar, psf1=self.psf1_sigma, psf2=self.psf2_sigma,
+            n_sources=0, seed=22, varSourceChange=0, psf_yvary_factor=0)
+
     def _computeVarianceMean(self, maskedIm):
         statObj = afwMath.makeStatistics(maskedIm.getVariance(),
                                          maskedIm.getMask(), afwMath.MEANCLIP,
@@ -347,6 +399,38 @@ def testDiffimCorrection(self):
         # Template variance is higher than that of the science img.
         self._testDiffimCorrection(svar=0.04, tvar=0.08)
 
+    def testNoiseDiffimCorrection(self):
+        """Test correction by estimating correlation directly on a noise difference image.
+
+        Notes
+        ------
+
+        See `lsst-dm/diffimTests` notebook `DM-24371_correlation_estimate.ipynb`
+        for further details of how the correlation looks like in the uncorrected
+        and corrected cases and where the tolerance numbers come from.
+        """
+        svar = 1.
+        tvar = 100.
+        self._setUpSourcelessImages(svar=svar, tvar=tvar)
+        diffExp, mKernel, expected_var = self._makeAndTestUncorrectedDiffim()
+        corrected_diffExp = self._runDecorrelationTask(diffExp, mKernel)
+
+        rho_sci = estimatePixelCorrelation(self.im1ex.getImage().getArray())
+        rho_rawdiff = estimatePixelCorrelation(diffExp.getImage().getArray())
+        rho_corrdiff = estimatePixelCorrelation(corrected_diffExp.getImage().getArray())
+
+        self.assertFloatsAlmostEqual(rho_sci[0], 1., atol=0.1, rtol=None)
+        self.assertFloatsAlmostEqual(rho_rawdiff[0], 1., atol=0.1, rtol=None)
+        self.assertFloatsAlmostEqual(rho_corrdiff[0], 1., atol=0.1, rtol=None)
+
+        self.assertFloatsAlmostEqual(rho_sci[1:], 0., atol=0.1, rtol=None)
+
+        self.assertGreater(rho_rawdiff[1], 0.2)
+        self.assertGreater(rho_rawdiff[2], 0.2)
+        self.assertGreater(rho_rawdiff[3], 0.2)
+
+        self.assertFloatsAlmostEqual(rho_corrdiff[1:], 0., atol=0.1, rtol=None)
+
     def _runDecorrelationTaskMapReduced(self, diffExp, mKernel):
         """ Run decorrelation using the imageMapReducer.
         """