diff --git a/include/lsst/meas/algorithms/CorrelatedNoiseDetection.h b/include/lsst/meas/algorithms/CorrelatedNoiseDetection.h
new file mode 100644
index 000000000..20f0cd263
--- /dev/null
+++ b/include/lsst/meas/algorithms/CorrelatedNoiseDetection.h
@@ -0,0 +1,121 @@
+// -*- lsst-c++ -*-
+/*
+ * LSST Data Management System
+ * Copyright 2017 LSST/AURA.
+ *
+ * This product includes software developed by the
+ * LSST Project (http://www.lsst.org/).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the LSST License Statement and
+ * the GNU General Public License along with this program. If not,
+ * see .
+ */
+
+#ifndef LSST_MEAS_ALGORITHMS_CorrelatedNoiseDetection_h_INCLUDED
+#define LSST_MEAS_ALGORITHMS_CorrelatedNoiseDetection_h_INCLUDED
+
+#include
+#include
+
+#include "lsst/afw/image/MaskedImage.h"
+
+namespace lsst { namespace meas { namespace algorithms {
+
+/**
+ * Estimate the noise correlation kernel of an image, assuming it is stationary.
+ *
+ * @param[in] image The image to be measured. The empirical covariance
+ * will be divided by the values in image's variance plane,
+ * and its mask will be used to reject pixels containing
+ * objects or artifacts.
+ * @param[in] radius Distance in pixels to which the correlation is measured
+ * on either side; the returned image will have dimensions
+ * (2*radius + 1, 2*radius + 1).
+ * @param[in] badBitMask A bit mask indicating pixels that should not be included
+ * in the measurement. Should generally include at least
+ * DETECTED.
+ *
+ * @return the noise correlation kernel: an image in which the central pixel
+ * represents the fraction of the total variance/covariance in the variance,
+ * and neighboring pixels contain the correlation at different offsets.
+ */
+afw::image::Image measureCorrelationKernel(
+ afw::image::MaskedImage const & image,
+ int radius,
+ afw::image::MaskPixel badBitMask
+);
+
+/**
+ * Estimate the noise correlation kernel of an image, assuming it is stationary.
+ *
+ * @param[in] image The image to be measured. The empirical covariance
+ * will be divided by the values in image's variance plane,
+ * and its mask will be used to reject pixels containing
+ * objects or artifacts.
+ * @param[in] radius Distance in pixels to which the correlation is measured
+ * on either side; the returned image will have dimensions
+ * (2*radius + 1, 2*radius + 1).
+ * @param[in] badMaskPlanes A list of mask planenames indicating pixels that
+ * should not be included in the measurement.
+ * Should generally include at least DETECTED.
+ *
+ * @return the noise correlation kernel: an image in which the central pixel
+ * represents the fraction of the total variance/covariance in the variance,
+ * and neighboring pixels contain the correlation at different offsets.
+ */
+afw::image::Image measureCorrelationKernel(
+ afw::image::MaskedImage const & image,
+ int radius,
+ std::vector const & badMaskPlanes
+);
+
+
+/**
+ * Fit an optimal detection kernel for a PSF that corrects for correlated noise.
+ *
+ * This is the same as the kernel that yields the PSF when convolved with the
+ * correlation kernel.
+ *
+ * The returned kernel cannot be used in detection in quite the same way as the
+ * PSF is used in detection on images with noise correlated noise.
+ * In the uncorrelated noise case with image @f$z@f$, PSF @f$\phi@f$, and
+ * per-pixel variance @f$\sigma^2@f$, the significance image is
+ * @f[
+ * \nu = \frac{\phi \ast z}{\sigma \sqrt{\phi \cdot \phi}}
+ * @f]
+ * In the correlated noise case, with @f$\psi@f$ the kernel computed by this
+ * function, the significance image is
+ * @f[
+ * \nu = \frac{\psi \ast z}{\sigma \sqrt{\phi \cdot \psi}}
+ * @f]
+ * (note the difference in the denominator).
+ *
+ * @param[in] psf Image of the PSF model.
+ * @param[in] correlation Noise correlation kernel, of the sort estimated by
+ * measureCorrelationKernel.
+ * @param[in] radius Radius of the detection kernel image in pixels.
+ * The returned image will have dimensions
+ * (2*radius + 1, 2*radius + 1).
+ *
+ * @return an image containing the optimal detection kernel.
+ */
+afw::image::Image fitGeneralDetectionKernel(
+ afw::image::Image const & psf,
+ afw::image::Image const & correlation,
+ int radius
+);
+
+
+}}} // namespace lsst::meas::algorithms
+
+#endif // !LSST_MEAS_ALGORITHMS_CorrelatedNoiseDetection_h_INCLUDED
diff --git a/python/lsst/meas/algorithms/SConscript b/python/lsst/meas/algorithms/SConscript
index 7796502fd..022805ab1 100644
--- a/python/lsst/meas/algorithms/SConscript
+++ b/python/lsst/meas/algorithms/SConscript
@@ -4,6 +4,7 @@ scripts.BasicSConscript.pybind11(["binnedWcs",
"cr",
"coaddBoundedField",
"coaddPsf/coaddPsf",
+ "correlatedNoiseDetection",
"doubleGaussianPsf",
"imagePsf",
"interp",
diff --git a/python/lsst/meas/algorithms/__init__.py b/python/lsst/meas/algorithms/__init__.py
index cbf9a59a1..1342db2b8 100644
--- a/python/lsst/meas/algorithms/__init__.py
+++ b/python/lsst/meas/algorithms/__init__.py
@@ -63,6 +63,7 @@
from .loadIndexedReferenceObjects import *
from .indexerRegistry import *
from .reserveSourcesTask import *
+from .correlatedNoiseDetection import *
from .version import *
diff --git a/python/lsst/meas/algorithms/correlatedNoiseDetection.cc b/python/lsst/meas/algorithms/correlatedNoiseDetection.cc
new file mode 100644
index 000000000..4ab1af358
--- /dev/null
+++ b/python/lsst/meas/algorithms/correlatedNoiseDetection.cc
@@ -0,0 +1,66 @@
+/*
+ * LSST Data Management System
+ * Copyright 2017 LSST/AURA.
+ *
+ * This product includes software developed by the
+ * LSST Project (http://www.lsst.org/).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the LSST License Statement and
+ * the GNU General Public License along with this program. If not,
+ * see .
+ */
+#include "pybind11/pybind11.h"
+#include "pybind11/stl.h"
+
+#include "lsst/meas/algorithms/CorrelatedNoiseDetection.h"
+
+namespace py = pybind11;
+using namespace pybind11::literals;
+
+namespace lsst {
+namespace meas {
+namespace algorithms {
+
+PYBIND11_PLUGIN(correlatedNoiseDetection) {
+ py::module::import("lsst.afw.image");
+
+ py::module mod("correlatedNoiseDetection");
+
+ mod.def(
+ "measureCorrelationKernel",
+ (afw::image::Image (*)(
+ afw::image::MaskedImage const &,
+ int,
+ afw::image::MaskPixel
+ ))&measureCorrelationKernel,
+ "image"_a, "radius"_a, "badBitMask"_a
+ );
+
+ mod.def(
+ "measureCorrelationKernel",
+ (afw::image::Image (*)(
+ afw::image::MaskedImage const &,
+ int,
+ std::vector const &
+ ))&measureCorrelationKernel,
+ "image"_a, "radius"_a, "badMaskPlanes"_a
+ );
+
+ mod.def("fitGeneralDetectionKernel", &fitGeneralDetectionKernel, "psf"_a, "correlations"_a, "radius"_a);
+
+ return mod.ptr();
+}
+
+} // algorithms
+} // meas
+} // lsst
diff --git a/src/CorrelatedNoiseDetection.cc b/src/CorrelatedNoiseDetection.cc
new file mode 100644
index 000000000..f0fca6927
--- /dev/null
+++ b/src/CorrelatedNoiseDetection.cc
@@ -0,0 +1,192 @@
+// -*- LSST-C++ -*-
+/*
+ * LSST Data Management System
+ * Copyright 2017 LSST/AURA.
+ *
+ * This product includes software developed by the
+ * LSST Project (http://www.lsst.org/).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the LSST License Statement and
+ * the GNU General Public License along with this program. If not,
+ * see .
+ */
+
+#include "lsst/pex/exceptions.h"
+#include "lsst/afw/math/LeastSquares.h"
+#include "lsst/meas/algorithms/CorrelatedNoiseDetection.h"
+
+namespace lsst { namespace meas { namespace algorithms {
+
+
+afw::image::Image measureCorrelationKernel(
+ afw::image::MaskedImage const & mi,
+ int radius,
+ afw::image::MaskPixel badBitMask
+) {
+ afw::image::Image result(2*radius + 1, 2*radius + 1);
+ afw::image::Image count(result.getDimensions());
+ int const width = mi.getWidth();
+ int const height = mi.getHeight();
+ // iterate over pixels in the MaskedImage, skipping any that meet our bad mask criteria
+ for (int y1 = 0; y1 < height; ++y1) {
+ int const y2a = std::max(0, y1 - radius);
+ int const y2b = std::min(height, y1 + radius + 1);
+ for (int x1 = 0; x1 < width; ++x1) {
+ if ((*mi.getMask())(x1, y1) & badBitMask) {
+ continue;
+ }
+ float z = (*mi.getImage())(x1, y1) / (*mi.getVariance())(x1, y1);
+ // iterate over neighboring pixels, with bounds set to avoid image boundaries
+ int const x2a = std::max(0, x1 - radius);
+ int const x2b = std::min(height, x1 + radius + 1);
+ for (int y2 = y2a; y2 < y2b; ++y2) {
+ auto miIter = mi.row_begin(y2) + x2a;
+ auto outIter = result.row_begin(radius + y2 - y1) + radius + x2a - x1;
+ auto countIter = count.row_begin(radius + y2 - y1) + radius + x2a - x1;
+ for (int x2 = x2a; x2 < x2b; ++x2, ++miIter, ++outIter, ++countIter) {
+ if (miIter.mask() & badBitMask) {
+ continue;
+ }
+ *outIter += z*miIter.image();
+ *countIter += 1;
+ }
+ }
+ }
+ }
+ result.getArray().deep() /= count.getArray();
+ result.setXY0(-radius, -radius);
+ return result;
+}
+
+afw::image::Image measureCorrelationKernel(
+ afw::image::MaskedImage const & image,
+ int radius,
+ std::vector const & badMaskPlanes
+) {
+ afw::image::MaskPixel mask = 0x0;
+ for (auto plane : badMaskPlanes) {
+ mask |= image.getMask()->getPlaneBitMask(plane);
+ }
+ return measureCorrelationKernel(image, radius, mask);
+}
+
+
+namespace {
+
+afw::geom::Extent2I getOddBoxHalfWidth(afw::geom::Box2I const & box, std::string const & name) {
+ afw::geom::Extent2I r((box.getWidth() - 1) / 2, (box.getHeight() - 1) / 2);
+ if (r.getX()*2 + 1 != box.getWidth()) {
+ throw LSST_EXCEPT(
+ pex::exceptions::InvalidParameterError,
+ name + " image width must be an odd integer"
+ );
+ }
+ if (r.getY()*2 + 1 != box.getHeight()) {
+ throw LSST_EXCEPT(
+ pex::exceptions::InvalidParameterError,
+ name + " image height must be an odd integer"
+ );
+ }
+ return r;
+}
+
+} // anonymous
+
+
+afw::image::Image fitGeneralDetectionKernel(
+ afw::image::Image const & psf,
+ afw::image::Image const & correlation,
+ int radius
+) {
+ auto const psfR = getOddBoxHalfWidth(psf.getBBox(), "PSF");
+ auto const corrR = getOddBoxHalfWidth(correlation.getBBox(), "Correlation kernel");
+ afw::geom::Extent2I const outR(radius, radius);
+ if (psfR.getX() <= corrR.getX()) {
+ throw LSST_EXCEPT(
+ pex::exceptions::InvalidParameterError,
+ "PSF image width must be greater than correlation kernel width"
+ );
+ }
+ if (psfR.getY() <= corrR.getY()) {
+ throw LSST_EXCEPT(
+ pex::exceptions::InvalidParameterError,
+ "PSF image height must be greater than correlation kernel height"
+ );
+ }
+ if (psfR.getX() <= outR.getX()) {
+ throw LSST_EXCEPT(
+ pex::exceptions::InvalidParameterError,
+ "PSF image width must be greater than output kernel width"
+ );
+ }
+ if (psfR.getY() <= outR.getY()) {
+ throw LSST_EXCEPT(
+ pex::exceptions::InvalidParameterError,
+ "PSF image height must be greater than output kernel height"
+ );
+ }
+
+ int const psfN = (psfR.getX()*2 + 1)*(psfR.getY()*2 + 1);
+ int const outN = (outR.getX()*2 + 1)*(outR.getY()*2 + 1);
+
+ // Get locators at the center of each image, so we can use indices with
+ // the origin at the center and make the code easier to read.
+ auto psfL = psf.xy_at(psfR.getX(), psfR.getY());
+ auto corrL = correlation.xy_at(corrR.getX(), corrR.getY());
+
+ // rhs is the PSF image, with rows and columns flattened into one dimension
+ Eigen::VectorXd rhs = Eigen::VectorXd::Zero(psfN);
+
+ // matrix represents convolution with the correlated noise kernel image
+ Eigen::MatrixXd matrix = Eigen::MatrixXd::Zero(psfN, outN);
+ int xyN = 0;
+ for (int y = -psfR.getY(); y <= psfR.getY(); ++y) {
+ for (int x = -psfR.getX(); x <= psfR.getX(); ++x) {
+ int ijN = 0;
+ for (int i = -outR.getY(); i <= outR.getY(); ++i) {
+ for (int j = -outR.getX(); j <= outR.getX(); ++j) {
+ // Could move these checks into the inner loop bounds for performance,
+ // but this is easier to read and probably fast enough.
+ if (std::abs(y - i) <= corrR.getY() && std::abs(x - j) <= corrR.getX()) {
+ matrix(xyN, ijN) = corrL(x - j, y - i);
+ }
+ ++ijN;
+ }
+ }
+ rhs[xyN] = psfL(x, y);
+ ++xyN;
+ }
+ }
+
+ // solve for the kernel that produces the PSF when convolved with the
+ // noise correlation kernel
+ auto lstsq = afw::math::LeastSquares::fromDesignMatrix(matrix, rhs);
+ auto solution = lstsq.getSolution();
+
+ // copy the result from the flattened solution vector into an image
+ afw::image::Image result(outR.getX()*2 + 1, outR.getY()*2 + 1);
+ auto outL = result.xy_at(outR.getX(), outR.getY());
+ xyN = 0;
+ for (int y = -outR.getY(); y <= outR.getY(); ++y) {
+ for (int x = -outR.getX(); x <= outR.getX(); ++x) {
+ outL(x, y) = solution[xyN];
+ ++xyN;
+ }
+ }
+
+ result.setXY0(-outR.getX(), -outR.getY());
+ return result;
+}
+
+
+}}} // namespace lsst::meas::algorithms
diff --git a/tests/test_correlatedNoiseDetection.py b/tests/test_correlatedNoiseDetection.py
new file mode 100755
index 000000000..4f0cf8ccf
--- /dev/null
+++ b/tests/test_correlatedNoiseDetection.py
@@ -0,0 +1,103 @@
+#
+# LSST Data Management System
+#
+# Copyright 2017 AURA/LSST.
+#
+# This product includes software developed by the
+# LSST Project (http://www.lsst.org/).
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the LSST License Statement and
+# the GNU General Public License along with this program. If not,
+# see .
+#
+from __future__ import absolute_import, division, print_function
+
+import unittest
+import numpy as np
+import scipy.ndimage
+
+import lsst.afw.image
+import lsst.meas.algorithms
+import lsst.utils.tests
+
+try:
+ type(display)
+except NameError:
+ display = False
+
+
+class CorrelatedNoiseDetectionTestCase(lsst.utils.tests.TestCase):
+ """Test utilities for dealing with correlated noise when running detection."""
+
+ def testMeasureCorrelatationKernel(self):
+ random = np.random.RandomState(50)
+ n = 300
+ mi = lsst.afw.image.MaskedImage(n, n, dtype=np.float32)
+ rms = 5.0
+ mi.image.array[:, :] = rms*random.randn(n, n)
+ mi.variance.array[:, :] = rms**2
+ plane = "DETECTED"
+ bitmask = mi.mask.getPlaneBitMask(plane)
+ # Randomly mask some pixels and give them higher values so we'll know if
+ # masking them out doesn't work.
+ m = (random.randn(n, n) > 1)
+ mi.mask.array[:, :] = bitmask*m
+ mi.image.array[m] += 10.0
+ radius = 5
+ kernel1 = lsst.meas.algorithms.measureCorrelationKernel(mi, radius=radius, badBitMask=bitmask)
+ kernel2 = lsst.meas.algorithms.measureCorrelationKernel(mi, radius=radius, badMaskPlanes=[plane])
+ if display:
+ from lsst.afw.display import Display
+ d = Display(frame=0)
+ d.mtv(kernel1)
+ self.assertEqual(kernel1.array.shape, (2*radius + 1, 2*radius + 1))
+ self.assertImagesEqual(kernel1, kernel2)
+ self.assertFloatsAlmostEqual(kernel1.array[5, 5], 1.0, atol=1E-2)
+ self.assertFloatsAlmostEqual(kernel1.array[:5, :], 0.0, atol=1E-2)
+ self.assertFloatsAlmostEqual(kernel1.array[6:, :], 0.0, atol=1E-2)
+ self.assertFloatsAlmostEqual(kernel1.array[:, :5], 0.0, atol=1E-2)
+ self.assertFloatsAlmostEqual(kernel1.array[:, 6:], 0.0, atol=1E-2)
+
+ def testFitGeneralDetectionKernel(self):
+ psfImage = lsst.afw.detection.GaussianPsf(41, 41, sigma=5.0).computeKernelImage()
+ corrImage = lsst.afw.detection.GaussianPsf(15, 15, sigma=2.0).computeKernelImage().convertF()
+ detImage = lsst.meas.algorithms.fitGeneralDetectionKernel(psfImage, corrImage, radius=19)
+ paddedArray = np.zeros((41, 41), dtype=np.float64)
+ paddedArray[1:-1, 1:-1] = detImage.array
+ checkArray = scipy.ndimage.convolve(paddedArray, corrImage.array, mode='constant')
+ checkImage = lsst.afw.image.ImageD(checkArray.shape[1], checkArray.shape[0])
+ checkImage.array[:, :] = checkArray
+ if display:
+ from lsst.afw.display import Display
+ d1 = Display(frame=1)
+ d1.mtv(psfImage)
+ d2 = Display(frame=2)
+ d2.mtv(corrImage)
+ d3 = Display(frame=3)
+ d3.mtv(detImage)
+ d4 = Display(frame=4)
+ d4.mtv(checkImage)
+ self.assertImagesAlmostEqual(psfImage, checkImage, atol=1E-5)
+
+
+class TestMemory(lsst.utils.tests.MemoryTestCase):
+ pass
+
+
+def setup_module(module):
+ lsst.utils.tests.init()
+
+
+if __name__ == "__main__":
+ lsst.utils.tests.init()
+ unittest.main()