Coverage for python / lsst / meas / algorithms / scaleVariance.py: 28%

1# 

2# LSST Data Management System 

3# Copyright 2022 AURA/LSST. 

4# 

5# This product includes software developed by the 

6# LSST Project (http://www.lsst.org/). 

7# 

8# This program is free software: you can redistribute it and/or modify 

9# it under the terms of the GNU General Public License as published by 

10# the Free Software Foundation, either version 3 of the License, or 

11# (at your option) any later version. 

12# 

13# This program is distributed in the hope that it will be useful, 

14# but WITHOUT ANY WARRANTY; without even the implied warranty of 

15# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

16# GNU General Public License for more details. 

17# 

18# You should have received a copy of the LSST License Statement and 

19# the GNU General Public License along with this program. If not, 

20# see <http://www.lsstcorp.org/LegalNotices/>. 

21# 

22from contextlib import contextmanager 

23import numpy as np 

24 

25from lsst.pex.config import Config, Field, ListField, ConfigurableField 

26from lsst.pipe.base import Task, Struct, AlgorithmError 

27from . import SubtractBackgroundTask 

28 

29__all__ = ["ScaleVarianceConfig", "ScaleVarianceTask", "ExceedsMaxVarianceScaleError"] 

30 

31 

32class ExceedsMaxVarianceScaleError(AlgorithmError): 

33 """Raised if ScaleVariance exceeds a specified threshold. 

34 

35 Parameters 

36 ---------- 

37 maxScaling: `float` 

38 Maximum variance scaling. 

39 """ 

40 def __init__(self, scaleVarianceValue, scaleVarianceLimit, **kwargs): 

41 msg = (f"Variance rescaling factor ({scaleVarianceValue}) exceeds configured limit " 

42 f"({scaleVarianceLimit})") 

43 

44 self.msg = msg 

45 self._metadata = kwargs 

46 super().__init__(msg, **kwargs) 

47 self._metadata["scaleVarianceValue"] = scaleVarianceValue 

48 self._metadata["scaleVarianceLimit"] = scaleVarianceLimit 

49 

50 def __str__(self): 

51 # Exception doesn't handle **kwargs, so we need a custom str. 

52 return f"{self.msg}: {self.metadata}" 

53 

54 @property 

55 def metadata(self): 

56 for key, value in self._metadata.items(): 

57 if not (isinstance(value, int) or isinstance(value, float) or isinstance(value, str)): 

58 raise TypeError(f"{key} is of type {type(value)}, but only (int, float, str) are allowed.") 

59 return self._metadata 

60 

61 

62class ScaleVarianceConfig(Config): 

63 background = ConfigurableField(target=SubtractBackgroundTask, doc="Background subtraction") 

64 maskPlanes = ListField[str]( 

65 default=["DETECTED", "DETECTED_NEGATIVE", "BAD", "SAT", "NO_DATA", "INTRP"], 

66 doc="Mask planes for pixels to ignore when scaling variance", 

67 ) 

68 limit = Field[float](default=10.0, doc="Maximum variance scaling value to permit") 

69 

70 def setDefaults(self): 

71 self.background.binSize = 32 

72 self.background.useApprox = False 

73 self.background.undersampleStyle = "REDUCE_INTERP_ORDER" 

74 self.background.ignoredPixelMask = ["DETECTED", "DETECTED_NEGATIVE", "BAD", "SAT", "NO_DATA", "INTRP"] 

75 

76 

77class ScaleVarianceTask(Task): 

78 """Scale the variance in a MaskedImage 

79 

80 The variance plane in a convolved or warped image (or a coadd derived 

81 from warped images) does not accurately reflect the noise properties of 

82 the image because variance has been lost to covariance. This Task 

83 attempts to correct for this by scaling the variance plane to match 

84 the observed variance in the image. This is not perfect (because we're 

85 not tracking the covariance) but it's simple and is often good enough. 

86 

87 The task implements a pixel-based and an image-based correction estimator. 

88 """ 

89 ConfigClass = ScaleVarianceConfig 

90 _DefaultName = "scaleVariance" 

91 

92 def __init__(self, *args, **kwargs): 

93 Task.__init__(self, *args, **kwargs) 

94 self.makeSubtask("background") 

95 

96 @contextmanager 

97 def subtractedBackground(self, maskedImage): 

98 """Context manager for subtracting the background 

99 

100 We need to subtract the background so that the entire image 

101 (apart from objects, which should be clipped) will have the 

102 image/sqrt(variance) distributed about zero. 

103 

104 This context manager subtracts the background, and ensures it 

105 is restored on exit. 

106 

107 Parameters 

108 ---------- 

109 maskedImage : `lsst.afw.image.MaskedImage` 

110 Image+mask+variance to have background subtracted and restored. 

111 

112 Returns 

113 ------- 

114 context : context manager 

115 Context manager that ensure the background is restored. 

116 """ 

117 bg = self.background.fitBackground(maskedImage) 

118 bgImage = bg.getImageF(self.background.config.algorithm, self.background.config.undersampleStyle) 

119 maskedImage -= bgImage 

120 try: 

121 yield 

122 finally: 

123 maskedImage += bgImage 

124 

125 def run(self, maskedImage): 

126 """Rescale the variance in a maskedImage in place. 

127 

128 Parameters 

129 ---------- 

130 maskedImage : `lsst.afw.image.MaskedImage` 

131 Image for which to determine the variance rescaling factor. The image 

132 is modified in place. 

133 

134 Returns 

135 ------- 

136 factor : `float` 

137 Variance rescaling factor. 

138 

139 Raises 

140 ------ 

141 ExceedsMaxVarianceScaleError 

142 If the estimated variance rescaling factor by both methods exceed the 

143 configured limit. 

144 

145 Notes 

146 ----- 

147 The task calculates and applies the pixel-based correction unless 

148 it is over the ``config.limit`` threshold. In this case, the image-based 

149 method is applied. 

150 """ 

151 with self.subtractedBackground(maskedImage): 

152 factor = self.pixelBased(maskedImage) 

153 if factor > self.config.limit: 

154 self.log.warning("Pixel-based variance rescaling factor (%f) exceeds configured limit (%f); " 

155 "trying image-based method", factor, self.config.limit) 

156 factor = self.imageBased(maskedImage) 

157 if factor > self.config.limit: 

158 raise ExceedsMaxVarianceScaleError(factor, self.config.limit) 

159 self.log.info("Renormalizing variance by %f", factor) 

160 maskedImage.variance *= factor 

161 return factor 

162 

163 def computeScaleFactors(self, maskedImage): 

164 """Calculate and return both variance scaling factors without modifying the image. 

165 

166 Parameters 

167 ---------- 

168 maskedImage : `lsst.afw.image.MaskedImage` 

169 Image for which to determine the variance rescaling factor. 

170 

171 Returns 

172 ------- 

173 R : `lsst.pipe.base.Struct` 

174 - ``pixelFactor`` : `float` The pixel based variance rescaling factor 

175 or 1 if all pixels are masked or invalid. 

176 - ``imageFactor`` : `float` The image based variance rescaling factor 

177 or 1 if all pixels are masked or invalid. 

178 """ 

179 with self.subtractedBackground(maskedImage): 

180 pixelFactor = self.pixelBased(maskedImage) 

181 imageFactor = self.imageBased(maskedImage) 

182 return Struct(pixelFactor=pixelFactor, imageFactor=imageFactor) 

183 

184 def pixelBased(self, maskedImage): 

185 """Determine the variance rescaling factor from pixel statistics 

186 

187 We calculate SNR = image/sqrt(variance), and the distribution 

188 for most of the background-subtracted image should have a standard 

189 deviation of unity. We use the interquartile range as a robust estimator 

190 of the SNR standard deviation. The variance rescaling factor is the 

191 factor that brings that distribution to have unit standard deviation. 

192 

193 This may not work well if the image has a lot of structure in it, as 

194 the assumptions are violated. In that case, use an alternate 

195 method. 

196 

197 Parameters 

198 ---------- 

199 maskedImage : `lsst.afw.image.MaskedImage` 

200 Image for which to determine the variance rescaling factor. 

201 

202 Returns 

203 ------- 

204 factor : `float` 

205 Variance rescaling factor or 1 if all pixels are masked or non-finite. 

206 

207 """ 

208 maskVal = maskedImage.mask.getPlaneBitMask(self.config.maskPlanes) 

209 isGood = (((maskedImage.mask.array & maskVal) == 0) 

210 & np.isfinite(maskedImage.image.array) 

211 & np.isfinite(maskedImage.variance.array) 

212 & (maskedImage.variance.array > 0)) 

213 

214 nGood = np.sum(isGood) 

215 self.log.debug("Number of selected background pixels: %d of %d.", nGood, isGood.size) 

216 if nGood < 2: 

217 # Not enough good data, np.percentile needs at least 2 points 

218 # to estimate a range 

219 return 1.0 

220 # Robust measurement of stdev using inter-quartile range 

221 snr = maskedImage.image.array[isGood]/np.sqrt(maskedImage.variance.array[isGood]) 

222 q1, q3 = np.percentile(snr, (25, 75)) 

223 stdev = 0.74*(q3 - q1) 

224 return stdev**2 

225 

226 def imageBased(self, maskedImage): 

227 """Determine the variance rescaling factor from image statistics 

228 

229 We calculate average(SNR) = stdev(image)/median(variance), and 

230 the value should be unity. We use the interquartile range as a robust 

231 estimator of the stdev. The variance rescaling factor is the 

232 factor that brings this value to unity. 

233 

234 This may not work well if the pixels from which we measure the 

235 standard deviation of the image are not effectively the same pixels 

236 from which we measure the median of the variance. In that case, use 

237 an alternate method. 

238 

239 Parameters 

240 ---------- 

241 maskedImage : `lsst.afw.image.MaskedImage` 

242 Image for which to determine the variance rescaling factor. 

243 

244 Returns 

245 ------- 

246 factor : `float` 

247 Variance rescaling factor or 1 if all pixels are masked or non-finite. 

248 """ 

249 maskVal = maskedImage.mask.getPlaneBitMask(self.config.maskPlanes) 

250 isGood = (((maskedImage.mask.array & maskVal) == 0) 

251 & np.isfinite(maskedImage.image.array) 

252 & np.isfinite(maskedImage.variance.array) 

253 & (maskedImage.variance.array > 0)) 

254 nGood = np.sum(isGood) 

255 self.log.debug("Number of selected background pixels: %d of %d.", nGood, isGood.size) 

256 if nGood < 2: 

257 # Not enough good data, np.percentile needs at least 2 points 

258 # to estimate a range 

259 return 1.0 

260 # Robust measurement of stdev 

261 q1, q3 = np.percentile(maskedImage.image.array[isGood], (25, 75)) 

262 ratio = 0.74*(q3 - q1)/np.sqrt(np.median(maskedImage.variance.array[isGood])) 

263 return ratio**2