Coverage for python/lsst/analysis/ap/imageQA.py: 23%
72 statements
« prev ^ index » next coverage.py v7.16.0, created at 2026-09-16 10:41 +0000
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/>.
22"""Pixel-level QA for AP difference images.
24Two utilities live here:
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 XOR, and a small summary.
31Both are intended for human-driven analysis from a notebook. The first is
32also designed to be mapped over a list of dataIds to produce a per-detector
33quality DataFrame.
34"""
36from __future__ import annotations
38__all__ = ["image_diff_stats", "image_diff_stats_table", "pixel_compare"]
40from dataclasses import dataclass
41from typing import Iterable, Mapping
43import numpy as np
44import pandas as pd
47def _finite_view(array):
48 """Return only the finite entries of ``array`` as a flat 1-D view."""
49 flat = np.asarray(array).ravel()
50 return flat[np.isfinite(flat)]
53def _basic_stats(array):
54 """Median, MAD, stdev, and kurtosis of the finite entries of ``array``.
56 Kurtosis is the Pearson definition (``E[((x-mu)/sigma)**4] - 3``); zero
57 for a Gaussian. Computed without scipy so this module has no extra deps.
58 """
59 finite = _finite_view(array)
60 if finite.size == 0:
61 return {"median": np.nan, "mad": np.nan, "stddev": np.nan, "kurtosis": np.nan}
62 median = float(np.median(finite))
63 mad = float(np.median(np.abs(finite - median)))
64 stddev = float(np.std(finite))
65 if stddev > 0:
66 z = (finite - finite.mean()) / stddev
67 kurtosis = float(np.mean(z**4) - 3.0)
68 else:
69 kurtosis = np.nan
70 return {"median": median, "mad": mad, "stddev": stddev, "kurtosis": kurtosis}
73def _mask_plane_fractions(mask):
74 """Return ``{plane_name: fraction_of_pixels_set}`` for an afw mask.
76 Parameters
77 ----------
78 mask : `lsst.afw.image.Mask`
79 The mask whose planes are to be summed.
81 Returns
82 -------
83 fractions : `dict` [`str`, `float`]
84 Per-plane fraction of pixels where that plane bit is set.
85 """
86 arr = mask.array
87 npix = arr.size
88 if npix == 0:
89 return {}
90 out = {}
91 for plane_name, plane_bit in mask.getMaskPlaneDict().items():
92 bitmask = np.uint64(1) << np.uint64(plane_bit)
93 out[f"frac_{plane_name}"] = float(np.sum((arr & bitmask) != 0) / npix)
94 return out
97def image_diff_stats(butler, visit, detector, dataset_name="difference_image"):
98 """Summary statistics on a single difference image.
100 Parameters
101 ----------
102 butler : `lsst.daf.butler.Butler`
103 Butler initialized with the relevant collections.
104 visit, detector : `int`
105 Data id selecting the exposure to load.
106 dataset_name : `str`, optional
107 Butler dataset type name to load. Defaults to ``"difference_image"``.
109 Returns
110 -------
111 stats : `dict`
112 ``{"visit", "detector", "median", "mad", "stddev", "kurtosis"}`` plus
113 a ``frac_<PLANE>`` entry per mask plane on the exposure.
114 """
115 diff = butler.get(dataset_name, {"visit": visit, "detector": detector})
116 stats = {"visit": visit, "detector": detector}
117 stats.update(_basic_stats(diff.image.array))
118 stats.update(_mask_plane_fractions(diff.mask))
119 return stats
122def image_diff_stats_table(butler,
123 data_ids: Iterable[Mapping[str, int]],
124 dataset_name="difference_image"):
125 """Run `image_diff_stats` over many data ids and return a DataFrame.
127 Parameters
128 ----------
129 butler : `lsst.daf.butler.Butler`
130 data_ids : iterable of mapping
131 Each mapping must contain at least ``visit`` and ``detector`` keys.
132 dataset_name : `str`, optional
133 See `image_diff_stats`.
135 Returns
136 -------
137 table : `pandas.DataFrame`
138 One row per dataId, indexed by ``(visit, detector)``. Mask-plane
139 columns may be missing on some rows if the corresponding plane is not
140 defined on every exposure; pandas fills those with NaN.
141 """
142 rows = []
143 for data_id in data_ids:
144 rows.append(image_diff_stats(butler,
145 data_id["visit"],
146 data_id["detector"],
147 dataset_name=dataset_name))
148 if not rows:
149 return pd.DataFrame()
150 return pd.DataFrame(rows).set_index(["visit", "detector"]).sort_index()
153@dataclass
154class PixelCompareResult:
155 """Container for `pixel_compare` outputs.
157 Attributes
158 ----------
159 img1, img2 : `lsst.afw.image.Exposure`
160 The two loaded exposures, in case the caller wants their masks/WCS.
161 diff : `numpy.ndarray`
162 Pixel array of ``img1 - img2``.
163 ratio : `numpy.ndarray`
164 Pixel array of ``img1 / img2``; NaN where ``img2 == 0``.
165 mask_diff : `numpy.ndarray`
166 XOR of the two mask arrays: nonzero pixels are where the mask planes
167 differ between the two images.
168 summary : `dict`
169 ``{"visit", "detector", "diff_median", "diff_mad", "diff_stddev",
170 "n_mask_pixels_changed", "frac_mask_pixels_changed"}``.
171 """
172 img1: object
173 img2: object
174 diff: np.ndarray
175 ratio: np.ndarray
176 mask_diff: np.ndarray
177 summary: dict
180def pixel_compare(butler1, butler2, visit, detector,
181 dataset_name="difference_image"):
182 """Compare two difference images for the same data id, pixel by pixel.
184 Pull the same dataId from two butlers (or two collections of one butler)
185 and inspect the residuals.
187 Parameters
188 ----------
189 butler1, butler2 : `lsst.daf.butler.Butler`
190 Two butlers, possibly the same instance with different collections.
191 visit, detector : `int`
192 Data id to load from both butlers.
193 dataset_name : `str`, optional
194 Butler dataset type name to load from each butler. Defaults to
195 ``"difference_image"``.
197 Returns
198 -------
199 result : `PixelCompareResult`
201 Raises
202 ------
203 ValueError
204 If the two loaded images differ in shape (cannot be aligned by simple
205 subtraction).
206 """
207 img1 = butler1.get(dataset_name, {"visit": visit, "detector": detector})
208 img2 = butler2.get(dataset_name, {"visit": visit, "detector": detector})
210 a1 = img1.image.array
211 a2 = img2.image.array
212 if a1.shape != a2.shape:
213 raise ValueError(f"Image shapes differ for visit={visit} detector={detector}: "
214 f"{a1.shape} vs {a2.shape}")
216 diff = a1 - a2
217 with np.errstate(divide="ignore", invalid="ignore"):
218 ratio = np.where(a2 != 0, a1 / a2, np.nan)
220 mask_diff = img1.mask.array ^ img2.mask.array
222 finite = _finite_view(diff)
223 if finite.size:
224 diff_median = float(np.median(finite))
225 diff_mad = float(np.median(np.abs(finite - diff_median)))
226 diff_stddev = float(np.std(finite))
227 else:
228 diff_median = diff_mad = diff_stddev = np.nan
230 n_changed = int(np.count_nonzero(mask_diff))
231 summary = {
232 "visit": visit,
233 "detector": detector,
234 "diff_median": diff_median,
235 "diff_mad": diff_mad,
236 "diff_stddev": diff_stddev,
237 "n_mask_pixels_changed": n_changed,
238 "frac_mask_pixels_changed": (float(n_changed / mask_diff.size)
239 if mask_diff.size else 0.0),
240 }
241 return PixelCompareResult(img1=img1, img2=img2, diff=diff, ratio=ratio,
242 mask_diff=mask_diff, summary=summary)