Coverage for python/lsst/analysis/ap/imageQA.py: 78%

1# This file is part of analysis_ap. 

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 

22"""Pixel-level QA for AP difference images. 

23 

24Two utilities live here: 

25 

26- `image_diff_stats` reports basic statistics (median, MAD, stdev, kurtosis, 

27 mask-plane fractions) of a single difference image. 

28- `pixel_compare` aligns two difference images for the same data id and 

29 returns their pixel-wise difference, ratio, mask differences, and a small 

30 summary. 

31 

32Both are intended for human-driven analysis from a notebook. The first is 

33also designed to be mapped over a list of dataIds to produce a per-detector 

34quality DataFrame. 

35""" 

36 

37from __future__ import annotations 

38 

39__all__ = ["image_diff_stats", "image_diff_stats_table", "pixel_compare"] 

40 

41from dataclasses import dataclass 

42from typing import Iterable, Mapping 

43 

44import numpy as np 

45import pandas as pd 

46 

47from .imageUtils import get_visit_image 

48 

49 

50def _finite_view(array): 

51 """Return only the finite entries of ``array`` as a flat 1-D view.""" 

52 flat = np.asarray(array).ravel() 

53 return flat[np.isfinite(flat)] 

54 

55 

56def _basic_stats(array): 

57 """Median, MAD, stdev, and kurtosis of the finite entries of ``array``. 

58 

59 Kurtosis is the Pearson definition (``E[((x-mu)/sigma)**4] - 3``); zero 

60 for a Gaussian. Computed without scipy so this module has no extra deps. 

61 """ 

62 finite = _finite_view(array) 

63 if finite.size == 0: 63 ↛ 64line 63 didn't jump to line 64 because the condition on line 63 was never true

64 return {"median": np.nan, "mad": np.nan, "stddev": np.nan, "kurtosis": np.nan} 

65 median = float(np.median(finite)) 

66 mad = float(np.median(np.abs(finite - median))) 

67 stddev = float(np.std(finite)) 

68 if stddev > 0: 68 ↛ 72line 68 didn't jump to line 72 because the condition on line 68 was always true

69 z = (finite - finite.mean()) / stddev 

70 kurtosis = float(np.mean(z**4) - 3.0) 

71 else: 

72 kurtosis = np.nan 

73 return {"median": median, "mad": mad, "stddev": stddev, "kurtosis": kurtosis} 

74 

75 

76def _mask_plane_fractions(mask): 

77 """Return ``{"frac_<plane>": fraction_of_pixels_set}`` for a mask. 

78 

79 Parameters 

80 ---------- 

81 mask : `lsst.images.Mask` 

82 The mask whose planes are to be summed. 

83 

84 Returns 

85 ------- 

86 fractions : `dict` [`str`, `float`] 

87 Per-plane fraction of pixels where that plane bit is set. 

88 """ 

89 if mask.bbox.area == 0: 89 ↛ 90line 89 didn't jump to line 90 because the condition on line 89 was never true

90 return {} 

91 return {f"frac_{name}": float(np.mean(mask.get(name))) for name in mask.schema.names} 

92 

93 

94def _mask_difference(mask1, mask2): 

95 """Return where two masks differ in any plane. 

96 

97 Planes are matched by name; a plane absent from one mask counts as unset 

98 there. 

99 

100 Parameters 

101 ---------- 

102 mask1, mask2 : `lsst.images.Mask` 

103 Masks with the same shape. 

104 

105 Returns 

106 ------- 

107 differ : `numpy.ndarray` [`bool`] 

108 True where at least one plane is set in only one mask. 

109 """ 

110 differ = np.zeros(mask1.bbox.shape, dtype=bool) 

111 for name in set(mask1.schema.names) | set(mask2.schema.names): 

112 set1 = mask1.get(name) if name in mask1.schema.names else False 

113 set2 = mask2.get(name) if name in mask2.schema.names else False 

114 differ |= set1 != set2 

115 return differ 

116 

117 

118def image_diff_stats(butler, visit, detector, dataset_name="difference_image"): 

119 """Summary statistics on a single difference image. 

120 

121 Parameters 

122 ---------- 

123 butler : `lsst.daf.butler.Butler` 

124 Butler initialized with the relevant collections. 

125 visit, detector : `int` 

126 Data id selecting the exposure to load. 

127 dataset_name : `str`, optional 

128 Butler dataset type name to load. Defaults to ``"difference_image"``. 

129 

130 Returns 

131 ------- 

132 stats : `dict` 

133 ``{"visit", "detector", "median", "mad", "stddev", "kurtosis"}`` plus 

134 a ``frac_<PLANE>`` entry per mask plane on the image, named as in 

135 `lsst.images` (e.g. ``frac_SATURATED``). 

136 """ 

137 diff = get_visit_image(butler, dataset_name, {"visit": visit, "detector": detector}, 

138 image_type="difference") 

139 stats = {"visit": visit, "detector": detector} 

140 stats.update(_basic_stats(diff.image.array)) 

141 stats.update(_mask_plane_fractions(diff.mask)) 

142 return stats 

143 

144 

145def image_diff_stats_table(butler, 

146 data_ids: Iterable[Mapping[str, int]], 

147 dataset_name="difference_image"): 

148 """Run `image_diff_stats` over many data ids and return a DataFrame. 

149 

150 Parameters 

151 ---------- 

152 butler : `lsst.daf.butler.Butler` 

153 data_ids : iterable of mapping 

154 Each mapping must contain at least ``visit`` and ``detector`` keys. 

155 dataset_name : `str`, optional 

156 See `image_diff_stats`. 

157 

158 Returns 

159 ------- 

160 table : `pandas.DataFrame` 

161 One row per dataId, indexed by ``(visit, detector)``. Mask-plane 

162 columns may be missing on some rows if the corresponding plane is not 

163 defined on every exposure; pandas fills those with NaN. 

164 """ 

165 rows = [] 

166 for data_id in data_ids: 

167 rows.append(image_diff_stats(butler, 

168 data_id["visit"], 

169 data_id["detector"], 

170 dataset_name=dataset_name)) 

171 if not rows: 

172 return pd.DataFrame() 

173 return pd.DataFrame(rows).set_index(["visit", "detector"]).sort_index() 

174 

175 

176@dataclass 

177class PixelCompareResult: 

178 """Container for `pixel_compare` outputs. 

179 

180 Attributes 

181 ---------- 

182 img1, img2 : `lsst.images.DifferenceImage` 

183 The two loaded images, in case the caller wants their masks/WCS. 

184 diff : `numpy.ndarray` 

185 Pixel array of ``img1 - img2``. 

186 ratio : `numpy.ndarray` 

187 Pixel array of ``img1 / img2``; NaN where ``img2 == 0``. 

188 mask_diff : `numpy.ndarray` [`bool`] 

189 True where at least one mask plane, matched by name, is set in only 

190 one of the two images. 

191 summary : `dict` 

192 ``{"visit", "detector", "diff_median", "diff_mad", "diff_stddev", 

193 "n_mask_pixels_changed", "frac_mask_pixels_changed"}``. 

194 """ 

195 img1: object 

196 img2: object 

197 diff: np.ndarray 

198 ratio: np.ndarray 

199 mask_diff: np.ndarray 

200 summary: dict 

201 

202 

203def pixel_compare(butler1, butler2, visit, detector, 

204 dataset_name="difference_image"): 

205 """Compare two difference images for the same data id, pixel by pixel. 

206 

207 Pull the same dataId from two butlers (or two collections of one butler) 

208 and inspect the residuals. A legacy image is converted to `lsst.images` 

209 form when read, so it can be compared with an `lsst.images` one. 

210 

211 Parameters 

212 ---------- 

213 butler1, butler2 : `lsst.daf.butler.Butler` 

214 Two butlers, possibly the same instance with different collections. 

215 visit, detector : `int` 

216 Data id to load from both butlers. 

217 dataset_name : `str`, optional 

218 Butler dataset type name to load from each butler. Defaults to 

219 ``"difference_image"``. 

220 

221 Returns 

222 ------- 

223 result : `PixelCompareResult` 

224 

225 Raises 

226 ------ 

227 ValueError 

228 If the two loaded images differ in shape (cannot be aligned by simple 

229 subtraction). 

230 """ 

231 data_id = {"visit": visit, "detector": detector} 

232 img1 = get_visit_image(butler1, dataset_name, data_id, image_type="difference") 

233 img2 = get_visit_image(butler2, dataset_name, data_id, image_type="difference") 

234 

235 a1 = img1.image.array 

236 a2 = img2.image.array 

237 if a1.shape != a2.shape: 237 ↛ 238line 237 didn't jump to line 238 because the condition on line 237 was never true

238 raise ValueError(f"Image shapes differ for visit={visit} detector={detector}: " 

239 f"{a1.shape} vs {a2.shape}") 

240 

241 diff = a1 - a2 

242 with np.errstate(divide="ignore", invalid="ignore"): 

243 ratio = np.where(a2 != 0, a1 / a2, np.nan) 

244 

245 mask_diff = _mask_difference(img1.mask, img2.mask) 

246 

247 finite = _finite_view(diff) 

248 if finite.size: 248 ↛ 253line 248 didn't jump to line 253 because the condition on line 248 was always true

249 diff_median = float(np.median(finite)) 

250 diff_mad = float(np.median(np.abs(finite - diff_median))) 

251 diff_stddev = float(np.std(finite)) 

252 else: 

253 diff_median = diff_mad = diff_stddev = np.nan 

254 

255 n_changed = int(np.count_nonzero(mask_diff)) 

256 summary = { 

257 "visit": visit, 

258 "detector": detector, 

259 "diff_median": diff_median, 

260 "diff_mad": diff_mad, 

261 "diff_stddev": diff_stddev, 

262 "n_mask_pixels_changed": n_changed, 

263 "frac_mask_pixels_changed": (float(n_changed / mask_diff.size) 

264 if mask_diff.size else 0.0), 

265 } 

266 return PixelCompareResult(img1=img1, img2=img2, diff=diff, ratio=ratio, 

267 mask_diff=mask_diff, summary=summary)