LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
noise_covariance.py
Go to the documentation of this file.
1# This file is part of meas_algorithms.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
21#
22import itertools
23import warnings
24from contextlib import contextmanager
25from dataclasses import dataclass
26
27import lsst.afw.image
28import numpy as np
29from lsst.meas.algorithms import SubtractBackgroundTask
30from lsst.pex.config import Config, ConfigurableField, Field, ListField
31from lsst.pipe.base import Task
32
33__all__ = (
34 "ComputeNoiseCorrelationConfig",
35 "ComputeNoiseCorrelationTask",
36 "CorrelationMatrix",
37)
38
39
40@dataclass(frozen=True)
42 """A class holding correlation coefficients for a set of background pixels.
43
44 CorrelationMatrix is a dataclass that is initialized with a numpy ndarray
45 and provides some convenience methods for accessing the matrix elements.
46 A CorrelationMatrix instance is callable wth two integer values x and y,
47 which returns the <I(m,n) I(m+x, n+y) / sqrt( V(m,n) V(m+x,n+y) )>, where
48 I is the image, V is the variance plane and < > denotes the expectation
49 operator.
50
51 Parameters
52 ----------
53 array : `numpy.ndarray`
54 The matrix of correlation coefficients.
55 """
56
57 array: np.ndarray
58
59 @property
60 def shape(self) -> tuple[int, int]:
61 """The shape of the correlation matrix."""
62 return self.array.shape
63
64 def __call__(self, x: int, y: int) -> float:
65 return self.array[x, y]
66
67
69 background = ConfigurableField(
70 target=SubtractBackgroundTask, doc="Background subtraction"
71 )
72 maskPlanes = ListField[str](
73 default=["DETECTED", "DETECTED_NEGATIVE", "BAD", "SAT", "NO_DATA", "INTRP"],
74 doc="Mask planes for pixels to ignore when calculating correlations",
75 )
76 size = Field[int](default=5, doc="Size of the correlation matrix to produce")
77 scaleEmpiricalVariance = Field[bool](
78 default=False,
79 doc=(
80 "Scale down the correlation coefficients x by the empirical variance of the background "
81 "in addition to the variance plane?"
82 ),
83 )
84 subtractEmpiricalMean = Field[bool](
85 default=False, doc="Subtract the empirical mean in addition to the background?"
86 )
87
88 def setDefaults(self):
89 self.background.binSize = 32
90 self.background.useApprox = False
91 self.background.undersampleStyle = "REDUCE_INTERP_ORDER"
92 self.background.ignoredPixelMask = [
93 "DETECTED",
94 "DETECTED_NEGATIVE",
95 "BAD",
96 "SAT",
97 "NO_DATA",
98 "INTRP",
99 ]
100
101
103 """Compute the noise correlation coefficients in a MaskedImage
104
105 The variance plane in a convolved or warped image (or a coadd derived
106 from warped images) does not accurately reflect the noise properties of
107 the image because variance has been lost to covariance. This Task computes
108 a matrix of correlation coefficients of a desired size. It assumes that the
109 noise is (at least the correlation coefficients are) stationary and uses
110 spatial averaging to compute the correlation coefficients.
111 """
112
113 ConfigClass = ComputeNoiseCorrelationConfig
114 _DefaultName = "computeNoiseCorrelation"
115
116 def __init__(self, *args, **kwargs):
117 super().__init__(*args, **kwargs)
118 self.makeSubtask("background")
119 self.background: SubtractBackgroundTask
120
121 @contextmanager
122 def subtractedBackground(self, maskedImage: lsst.afw.image.MaskedImage):
123 """Context manager for subtracting the background
124
125 We need to subtract the background so that the entire image
126 (apart from objects, which should be clipped) will have the
127 image/sqrt(variance) distributed about zero with unit variance.
128 This context manager subtracts the background, and ensures it
129 is restored on exit.
130
131 Parameters
132 ----------
133 maskedImage : `lsst.afw.image.MaskedImage`
134 Image+mask+variance to have background subtracted and restored.
135
136 Returns
137 -------
138 context : context manager
139 Context manager that ensure the background is restored.
140 """
141 bg = self.background.fitBackground(maskedImage)
142 bgImage = bg.getImageF(
143 self.background.config.algorithm, self.background.config.undersampleStyle
144 )
145 maskedImage -= bgImage
146 try:
147 yield
148 finally:
149 maskedImage += bgImage
150
151 def run(
152 self,
153 maskedImage: lsst.afw.image.MaskedImage,
154 refMaskedImage: lsst.afw.image.MaskedImage | None = None,
155 ) -> CorrelationMatrix:
156 """Compute the correlation matrix from a maskedImage.
157
158 Parameters
159 ----------
160 maskedImage : `~lsst.afw.image.MaskedImage`
161 Image for which to determine the correlation matrix.
162 refMaskedImage : `~lsst.afw.image.MaskedImage`, optional
163 Image from which to determine which pixels to mask.
164 If None, it defaults to ``maskedImage``.
165
166 Returns
167 -------
168 corr_matrix : `CorrelationMatrix`
169 Correlation matrix of the maskedImage.
170
171 Raises
172 ------
173 RuntimeError
174 Raised if ``refMaskedImage`` is provided and does not have the same
175 dimensions as ``maskedImage``.
176 """
177 with self.subtractedBackground(maskedImage):
178 if refMaskedImage is None:
179 refMaskedImage = maskedImage
180 elif refMaskedImage.getDimensions() != maskedImage.getDimensions():
181 raise RuntimeError(
182 "Reference image has different dimensions than input image"
183 )
184
185 corr_matrix = self._pixelBased(maskedImage, refMaskedImage)
186
187 return corr_matrix
188
190 self,
191 maskedImage: lsst.afw.image.MaskedImage,
192 refMaskedImage: lsst.afw.image.MaskedImage,
193 ) -> CorrelationMatrix:
194 """Determine correlation coefficients between pixels
195
196 This is the concrete routine that does the computation.
197
198 Parameters
199 ----------
200 maskedImage : `~lsst.afw.image.MaskedImage`
201 Image for which to determine the variance rescaling factor.
202 refMaskedImage : `~lsst.afw.image.MaskedImage`
203 Image from which to determine which pixels to mask.
204
205 Returns
206 -------
207 corr_matrix : `CorrelationMatrix`
208 Correlation matrix of the maskedImage.
209 """
210 maskVal = refMaskedImage.mask.getPlaneBitMask(self.config.maskPlanes)
211 isGood = (
212 ((refMaskedImage.mask.array & maskVal) == 0)
213 & np.isfinite(refMaskedImage.image.array)
214 & np.isfinite(refMaskedImage.variance.array)
215 & (refMaskedImage.variance.array > 0)
216 )
217
218 nGood = np.sum(isGood)
219 self.log.debug(
220 "Number of selected background pixels: %d of %d.", nGood, isGood.size
221 )
222
223 # Catch any RuntimeWarning that may arise by dividing by zero.
224 # This is okay since we handle the zero variance case immediately.
225 with warnings.catch_warnings():
226 warnings.simplefilter("ignore", category=RuntimeWarning)
227 normalized_arr = maskedImage.image.array / np.sqrt(
228 maskedImage.variance.array
229 )
230 normalized_arr[~isGood] = np.nan
231
232 corr_matrix = np.empty(
233 (self.config.size + 1, self.config.size + 1), dtype=np.float32
234 )
235
236 for dx, dy in itertools.product(
237 range(self.config.size + 1), range(self.config.size + 1)
238 ):
239 sliceX = slice(None, -dx) if dx != 0 else slice(None, None)
240 sliceY = slice(None, -dy) if dy != 0 else slice(None, None)
241 arr1 = normalized_arr[sliceX, sliceY]
242
243 sliceX = slice(dx, None) if dx != 0 else slice(None, None)
244 sliceY = slice(dy, None) if dy != 0 else slice(None, None)
245 arr2 = normalized_arr[sliceX, sliceY]
246
247 if self.config.subtractEmpiricalMean:
248 arr1 -= np.nanmean(arr1)
249 arr2 -= np.nanmean(arr2)
250 if self.config.scaleEmpiricalVariance:
251 # Do not use nanstd direct, as it will subtract the
252 # empirical mean regardless of config set.
253 arr1 /= np.nanmean(arr1**2) ** 0.5
254 arr2 /= np.nanmean(arr2**2) ** 0.5
255
256 cov = np.nanmean(arr1 * arr2)
257 # Adjust for the bias in the estimator. Temporary reference:
258 # https://en.wikipedia.org/wiki/Pearson_correlation_coefficient#Practical_issues
259 # TODO: Explain this in the DMTN-215 (DM-33418).
260 cov *= 1.0 + 0.5 * (1 - cov**2) / (~np.isnan(arr1 * arr2)).sum()
261
262 corr_matrix[dx, dy] = cov
263
264 return CorrelationMatrix(corr_matrix)
A class to manipulate images, masks, and variance as a single object.
Definition MaskedImage.h:74
CorrelationMatrix run(self, lsst.afw.image.MaskedImage maskedImage, lsst.afw.image.MaskedImage|None refMaskedImage=None)
subtractedBackground(self, lsst.afw.image.MaskedImage maskedImage)
CorrelationMatrix _pixelBased(self, lsst.afw.image.MaskedImage maskedImage, lsst.afw.image.MaskedImage refMaskedImage)