Coverage for python/lsst/meas/algorithms/computeRoughPsfShapelets.py: 0%

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 

22from __future__ import annotations 

23 

24__all__ = ["ComputeRoughPsfShapeletsTask", "ComputeRoughPsfShapeletsConfig"] 

25 

26import itertools 

27from typing import TYPE_CHECKING, Any, Literal 

28 

29import numpy as np 

30from logging import Logger 

31import scipy.signal 

32import scipy.stats 

33from sklearn.covariance import MinCovDet 

34from sklearn.neighbors import KernelDensity 

35 

36from lsst.afw.geom import ellipses 

37from lsst.afw.image import ImageD, ImageF, MaskedImageF, ExposureSummaryStats 

38from lsst.afw.table import Point2DKey, QuadrupoleKey, Schema, SourceCatalog 

39from lsst.geom import Box2D 

40from lsst.pex.config import Config, Field, ListField, FieldValidationError 

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

42from lsst.shapelet import LAGUERRE, ShapeletFunction, computeOffset 

43 

44from ._algorithmsLib import SpanSetMoments 

45 

46if TYPE_CHECKING: 

47 import matplotlib.axes 

48 import matplotlib.figure 

49 import matplotlib.image 

50 import matplotlib.patches 

51 

52 

53class NoStarsForShapeletsError(AlgorithmError): 

54 """Exception raised when any of the methods involving selection fail to 

55 find any usable stars (e.g. compute_raw_moments, _threshold_with_bounds, 

56 _find_first_radius_mode). 

57 """ 

58 

59 @property 

60 def metadata(self) -> dict[str, Any]: 

61 return {} 

62 

63 

64class ComputeRoughPsfShapeletsConfig(Config): 

65 bad_mask_planes = ListField( 

66 "Mask planes to identify pixels to drop from the calculation.", 

67 dtype=str, 

68 default=["SAT", "SUSPECT", "INTRP"], 

69 ) 

70 bad_pixel_max_fraction = Field( 

71 "Maximum fraction of a footprint's pixels that can be bad (according " 

72 "to bad_mask_planes) before that footprint is fully dropped.", 

73 dtype=float, 

74 default=0.25, 

75 ) 

76 bad_pixel_exclusion_radius = Field( 

77 "If a bad pixel (according to bad_mask_planes) falls within this " 

78 "radius of the unweighted centroid of a footprint, drop that footprint " 

79 "entirely.", 

80 dtype=float, 

81 default=2.0, 

82 ) 

83 max_footprint_area = Field( 

84 "Footprints with a pixel count larger than this threshold are dropped before computing moments.", 

85 dtype=int, 

86 default=10000, 

87 ) 

88 min_snr = Field( 

89 "Mininum flux S/N for inclusion in the star sample.", 

90 dtype=float, 

91 default=50.0, 

92 ) 

93 max_radius_factor = Field( 

94 "Maximium multiple of the mode of the radius distribution for inclusion in the star sample.", 

95 dtype=float, 

96 default=2.0, 

97 ) 

98 max_shape_distance = Field( 

99 "Maximum Mahalanobis distance (distance from the center of the shape " 

100 "distribution in elliptical sigma units) to select a star from the candidate " 

101 "sample, comparing the shape of that source vs. a robust estimate of the " 

102 "distribution of the shapes of all sources.", 

103 dtype=float, 

104 default=3.0, 

105 ) 

106 min_n_stars = Field( 

107 "Minimum number of stars to select. The S/N, radius, and shape distance thresholds are " 

108 "relaxed as needed to meet this target.", 

109 dtype=int, 

110 default=10, 

111 ) 

112 max_n_stars = Field( 

113 "Maximum number of stars to select. High shape-distance sources " 

114 "are dropped as needed to meet this target.", 

115 dtype=int, 

116 default=20, 

117 ) 

118 logarithmic_shapes = Field( 

119 "If True, transform the (xx, yy, xy) moments to conformal shear " 

120 "and log trace radius when selecting stars in order to map the ellipse " 

121 "parameter to a space with (-inf, inf) bounds on all quantities (but " 

122 "a less-linear relationship to the pixel data).", 

123 dtype=bool, 

124 default=False, 

125 ) 

126 radius_mode_min = Field( 

127 "Minimum radius in pixels at which to start searching for the first mode. " 

128 "This just needs be large enough to avoid unphysically small garbage (e.g. CRs).", 

129 dtype=float, 

130 default=1.0, 

131 ) 

132 radius_kde_bandwidth = Field( 

133 "Bandwidth of the Gaussian kernel density estimator used to find the " 

134 "first mode of the radius distribution.", 

135 dtype=float, 

136 default=1.0, 

137 ) 

138 shapelet_order = Field( 

139 "Order of the shapelet expansion fit to the stars.", 

140 dtype=int, 

141 default=6, 

142 ) 

143 shapelet_scale_factor = Field( 

144 "Scale factor to apply to the moments ellipse when computing the ellipse for the shapelet basis.", 

145 dtype=float, 

146 default=1.0, 

147 ) 

148 shapelet_circular_basis = Field( 

149 "Whether to use a circular shapelet basis with the same moments trace instead of an elliptical one.", 

150 dtype=bool, 

151 default=True, 

152 ) 

153 non_gaussian_non_atmosphere_index_list = ListField( 

154 doc="Index list of non-Gaussian, non-atmospheric contributors in the list of shapelet " 

155 "coefficients associated with shapelet_order (default accommodates an order of 6)", 

156 dtype=int, 

157 default=list(range(6, 10)) + list(range(15, 24)), 

158 ) 

159 

160 def validate(self) -> None: 

161 if self.min_n_stars > self.max_n_stars: 

162 raise ValueError( 

163 f"min_n_stars={self.min_n_stars} is greater than max_n_stars={self.max_n_stars}." 

164 ) 

165 if self.shapelet_order < 0: 

166 raise ValueError(f"shapelet order {self.shapelet_order} must be nonnegative.") 

167 if self.shapelet_scale_factor <= 0.0: 

168 raise ValueError(f"shapelet scale factor {self.shapelet_scale_factor} must be positive.") 

169 # Check that the shapelet decomposition order results in at least as 

170 # many decomposition coefficients as the maximum index identified in 

171 # the non_gaussian_non_atmosphere_index_list config. 

172 n_coeff = 0 

173 for i_offset in range(self.shapelet_order + 1): 

174 n_coeff += (i_offset + 1) 

175 if max(self.non_gaussian_non_atmosphere_index_list) > n_coeff: 

176 raise FieldValidationError( 

177 ComputeRoughPsfShapeletsConfig, 

178 self, 

179 "The largest coefficient index defined in config.non_gaussian_non_atmosphere_index_list " 

180 f"({len(self.non_gaussian_non_atmosphere_index_list)}) is larger than the number " 

181 f"of coefficients expected ({n_coeff}) from the shapelet decomposition order " 

182 f"({self.shapelet_order})." 

183 ) 

184 

185 

186class ComputeRoughPsfShapeletsTask(Task): 

187 """A task that computes a rough shapelet expansion of the PSF from a set 

188 of high S/N detections. 

189 

190 Notes 

191 ----- 

192 This task is expected to be run early in single-epoch processing - just 

193 after background subtraction and an initial high S/N detection phase, and 

194 before any deblending or measurement - in order to identify out-of-focus 

195 or otherwise bad PSFs. 

196 

197 Given a background-subtracted `lsst.afw.image.MaskedImage`, an 

198 `lsst.afw.table.SourceCatalog` with footprints attached, and a random 

199 number generator seed, the `run` method will: 

200 

201 - Compute the *unweighted* 0th-2nd moments of every non-child source over 

202 the footprint (except certain configurable masked pixels). This is 

203 delegated to the `compute_raw_moments` method (which uses the C++ 

204 `SpanSetMoments` class for the pixel-level processing). Unweighted 

205 moments are used to avoid "latching onto" a small piece of PSF 

206 substructure, but can be much noiser than the Gaussian-weighed moments 

207 we usually use. 

208 

209 - Select a "candidate" sample of sources with successfully measured 

210 moments that satisfy a S/N cut and a radius cut (determined from the 

211 first mode of the radius distribution, via kernel density estimation), 

212 and then use a robust covariance estimator (`scikit_learn.MinCovDet`) 

213 to select presumed isolated stars that are close to the center of that 

214 distribution, in 3-parameter shape space. This is delegated to the 

215 `select_stars` method. 

216 

217 - Fit a single shapelet expansion to the selected stars. This is 

218 mostly delegated to the `SpanSetMoments.fit_shapelets` method. 

219 

220 The radial shapelet terms at 0th, 2nd, and 4th order are expected to form 

221 a space in which donut-shaped PSFs are well-separated from those with 

222 monotonic profiles. Other terms *may* be useful in identifying other 

223 kinds of undesirable PSF structure. 

224 """ 

225 

226 ConfigClass = ComputeRoughPsfShapeletsConfig 

227 _DefaultName = "computeRoughPsfShapelets" 

228 config: ComputeRoughPsfShapeletsConfig 

229 

230 def __init__( 

231 self, 

232 config: ComputeRoughPsfShapeletsConfig | None = None, 

233 *, 

234 schema: Schema, 

235 **kwargs: Any, 

236 ): 

237 super().__init__(config=config, **kwargs) 

238 self.schema = schema 

239 self.shapelets_base_name = "RoughPsfShapelets" 

240 self._flux_key = schema.addField( 

241 self.shapelets_base_name + "_flux", 

242 type=float, 

243 doc="Unweighted zeroth moment.", 

244 ) 

245 self._flux_err_key = schema.addField( 

246 self.shapelets_base_name + "_fluxErr", 

247 type=float, 

248 doc="Uncertainty on the unweighted zeroth moment.", 

249 ) 

250 self._center_key = Point2DKey.addFields( 

251 schema, self.shapelets_base_name, "Center from unweighted first moments.", "pixels" 

252 ) 

253 self._shape_key = QuadrupoleKey.addFields( 

254 schema, self.shapelets_base_name, "Shape from unweighted second moments." 

255 ) 

256 self._flag_key = schema.addField( 

257 self.shapelets_base_name + "_flag", 

258 type="Flag", 

259 doc="Flag set if the raw PSF moments were not computed.", 

260 ) 

261 self._candidate_key = schema.addField( 

262 self.shapelets_base_name + "_candidate", 

263 type="Flag", 

264 doc="Flag set if this source passed the radius_fraction cut (see configuration).", 

265 ) 

266 self._used_key = schema.addField( 

267 self.shapelets_base_name + "_used", 

268 type="Flag", 

269 doc=( 

270 "Flag set if this source passed the radius_fraction and shape_distance cuts " 

271 "(see configuration) and was used to fit the shapelet expansion." 

272 ), 

273 ) 

274 

275 def run(self, *, masked_image: MaskedImageF, catalog: SourceCatalog, seed: int) -> Struct: 

276 """Compute raw moments, select stars, and fit a shapelet expansion to 

277 them. 

278 

279 Parameters 

280 ---------- 

281 masked_image 

282 Masked image to measure on. Must be background-subtracted. 

283 catalog 

284 Catalog of detections to extract footprints from and fill output 

285 columns of. Its schema must be a superset of ``self.schema``. 

286 seed 

287 A random-number generator seed, used for the robust covariance 

288 estimator. 

289 

290 Returns 

291 ------- 

292 `lsst.pipe.base.Struct` 

293 A struct of results containing: 

294 

295 - ``shapelet`` (`lsst.shapelet.ShapeletFunction`): A 

296 Gauss-Laguerre (polar shaplet) expansion of the PSF. 

297 - ``radial`` (`list` [`float`]): the purely radial coefficients 

298 of the shapelet expansion. 

299 - all attributes returned by the `select_stars` method. 

300 """ 

301 moments = self.compute_raw_moments(masked_image=masked_image, catalog=catalog) 

302 result = self.select_stars(catalog, seed=seed) 

303 star_moments = [moments[star_id] for star_id in result.star_ids] 

304 result.shapelet = SpanSetMoments.fit_shapelets( 

305 masked_image, 

306 star_moments, 

307 self.config.shapelet_order, 

308 self.config.shapelet_scale_factor, 

309 self.config.shapelet_circular_basis, 

310 ) 

311 result.shapelet.getEllipse().setCore(result.mean_shape) 

312 result.shapelet.changeBasisType(LAGUERRE) 

313 result.radial = result.shapelet.getCoefficients()[ 

314 [computeOffset(i) for i in range(0, self.config.shapelet_order + 1, 2)] 

315 ] 

316 self.log.info("Rough PSF shapelet radial terms: %s.", result.radial) 

317 result = self.compute_shapelets_only_iq_score(result) 

318 return result 

319 

320 def compute_raw_moments( 

321 self, *, masked_image: MaskedImageF, catalog: SourceCatalog 

322 ) -> dict[int, SpanSetMoments]: 

323 """Compute the unweighted moments of the footprints in a catalog. 

324 

325 Parameters 

326 ---------- 

327 masked_image 

328 Masked image to measure on. Must be background-subtracted. 

329 catalog 

330 Catalog of detections to extract footprints from and fill output 

331 columns of. Its schema must be a superset of ``self.schema``. 

332 

333 Returns 

334 ------- 

335 `dict` [`int`, `SpanSetMoments`] 

336 Objects used to construct and hold the unweighted moments and the 

337 pixel region used to computed them, keyed by source ID. 

338 """ 

339 bitmask = masked_image.mask.getPlaneBitMask(self.config.bad_mask_planes) 

340 all_moments: dict[int, SpanSetMoments] = {} 

341 for record in catalog: 

342 if record.getParent() != 0: 

343 record.set(self._flag_key, True) 

344 self.log.debug("Skipping child source %s", record.getId()) 

345 continue 

346 footprint_spans = record.getFootprint().getSpans() 

347 if footprint_spans.getArea() > self.config.max_footprint_area: 

348 record.set(self._flag_key, True) 

349 self.log.debug( 

350 "Skipping source %s with footprint area %d > %d.", 

351 record.getId(), 

352 footprint_spans.getArea(), 

353 self.config.max_footprint_area, 

354 ) 

355 continue 

356 moments = SpanSetMoments.compute( 

357 record.getFootprint().getSpans(), 

358 masked_image=masked_image, 

359 bad_bitmask=bitmask, 

360 bad_pixel_max_fraction=self.config.bad_pixel_max_fraction, 

361 bad_pixel_exclusion_radius=self.config.bad_pixel_exclusion_radius, 

362 ) 

363 record.set(self._flux_key, moments.flux) 

364 record.set(self._flux_err_key, moments.variance**0.5) 

365 record.set(self._center_key, moments.center) 

366 record.set(self._shape_key, moments.shape) 

367 record.set(self._flag_key, moments.any_flags_set) 

368 if not moments.any_flags_set: 

369 all_moments[record.getId()] = moments 

370 if all_moments: 

371 self.log.verbose( 

372 "Successfully measured raw moments for %d of %d sources.", 

373 len(all_moments), 

374 len(catalog), 

375 ) 

376 else: 

377 raise NoStarsForShapeletsError("No raw moments could be measured.") 

378 return all_moments 

379 

380 def select_stars(self, catalog: SourceCatalog, seed: int) -> Struct: 

381 """Select probable stars from the distribution of second moments. 

382 

383 Parameters 

384 ---------- 

385 catalog 

386 Catalog of detections to extract footprints from and fill output 

387 columns of. Its schema must be a superset of ``self.schema``. 

388 seed 

389 A random-number generator seed, used for the robust covariance 

390 estimator. 

391 

392 Returns 

393 ------- 

394 `lsst.pipe.base.Struct` 

395 A struct of results containing: 

396 

397 - ``star_ids`` (`numpy.ndarray`): the source IDs that are expected 

398 to be stars. 

399 - ``mean_shape`` (`lsst.afw.geom.ellipses.BaseCore`): the mean of 

400 the shape distribution. 

401 - ``shape_covariance`` (`numpy.ndarray`): the covariance of the 

402 distribution of shapes; a 3x3 matrix. This uses the same 

403 parameterization of the shapes as ``mean_shape``. 

404 - ``radius_cut`` (`float`): the indended radius cut (i.e. the 

405 mode of the radius distribution multipled by the 

406 ``radius_factor`` configuration option). 

407 - ``radius_kde`` (`sklearn.neighbors.KernelDensity`): kernel 

408 density estimator on the radius distribution, used to determine 

409 the radius cut. 

410 """ 

411 # Cut on flags and SNR first. 

412 indices = np.arange(len(catalog), dtype=int)[np.logical_not(catalog[self._flag_key])] 

413 signalToNoise = [] 

414 for index in indices: 

415 if (np.isfinite(catalog[self._flux_err_key][index]) 

416 and np.isfinite(catalog[self._flux_err_key][index]) 

417 and catalog[self._flux_err_key][index] != 0.0): 

418 signalToNoise.append(catalog[self._flux_key][index] / catalog[self._flux_err_key][index]) 

419 else: 

420 signalToNoise.append(np.nan) 

421 signalToNoise = np.array(signalToNoise) 

422 indices = indices[ 

423 self._threshold_with_bounds( 

424 signalToNoise, 

425 threshold=self.config.min_snr, 

426 min_count=self.config.min_n_stars, 

427 max_count=len(catalog), 

428 name="S/N", 

429 kind=">", 

430 ) 

431 ] 

432 # Cut on radius next. 

433 radii = np.zeros(indices.size, dtype=np.float64) 

434 for n, index in enumerate(indices): 

435 record = catalog[index] 

436 shape = record.get(self._shape_key) 

437 radii[n] = shape.getTraceRadius() 

438 radius_mode, radius_kde = self._find_first_radius_mode(radii) 

439 radius_cut = self.config.max_radius_factor * radius_mode 

440 indices = indices[ 

441 self._threshold_with_bounds( 

442 radii, 

443 threshold=radius_cut, 

444 min_count=self.config.min_n_stars, 

445 max_count=len(catalog), 

446 name="radius", 

447 kind="<", 

448 ) 

449 ] 

450 shape_data = np.zeros((len(indices), 3), dtype=np.float64) 

451 for n, index in enumerate(indices): 

452 record = catalog[index] 

453 record.set(self._candidate_key, True) 

454 shape = record.get(self._shape_key) 

455 if self.config.logarithmic_shapes: 

456 shape = ellipses.SeparableConformalShearLogTraceRadius(shape) 

457 shape_data[n, :] = shape.getParameterVector() 

458 shape_dist = MinCovDet(random_state=seed).fit(shape_data) 

459 m_distances = shape_dist.mahalanobis(shape_data) 

460 indices = indices[ 

461 self._threshold_with_bounds( 

462 m_distances, 

463 threshold=self.config.max_shape_distance, 

464 min_count=self.config.min_n_stars, 

465 max_count=self.config.max_n_stars, 

466 name="Mahalanobis distance", 

467 kind="<", 

468 ) 

469 ] 

470 for index in indices: 

471 catalog[index].set(self._used_key, True) 

472 star_ids = catalog["id"][indices] 

473 if self.config.logarithmic_shapes: 

474 mean_shape = ellipses.SeparableConformalShearLogTraceRadius(shape_dist.location_) 

475 else: 

476 mean_shape = ellipses.Quadrupole(shape_dist.location_) 

477 return Struct( 

478 star_ids=star_ids, 

479 mean_shape=mean_shape, 

480 shape_covariance=shape_dist.covariance_, 

481 radius_cut=radius_cut, 

482 radius_kde=radius_kde, 

483 catalog=catalog, 

484 ) 

485 

486 def compute_shapelets_only_iq_score(self, results: Struct) -> Struct: 

487 """Compute the image quality "score" that only includes the 

488 non-atmospheric coefficents. 

489 

490 The "shapelet_only_iq_score" is a dimensionless image quality score 

491 as determined from the shapelets decomposition that only includes 

492 power from the non-atmospheric coefficents from the shapelet 

493 decomposition. The score spans the range [0.0, 1.0] with lower 

494 values indicating better image quality. 

495 

496 Parameters 

497 ---------- 

498 results 

499 Result struct returned by `run`. To be updated in-place. 

500 

501 Returns 

502 ------- 

503 `lsst.pipe.base.Struct` 

504 Additional results struct components include: 

505 - ``shapelets_summary`` (`lsst.afw.image.ExposureSummaryStats`): 

506 An ExposureSummaryStats object with the following parameters 

507 initialized: 

508 - ``shapeletsCoeffs`` (`list`[`float`]): The coefficients 

509 from the shapelet decomposition. 

510 - ``nShapeletsStar`` (`int`): The number of sources used in the 

511 shapelet decomposition. 

512 - ``shapeletsOnlyIqScore`` (`float`): The image quality score 

513 which only includes power from the non-atmospheric 

514 decomposition coefficients. 

515 - ``shapeletsStarEMedian`` (`float`): The median ellipticity 

516 (sqrt(starE1**2.0 + starE2**2.0)) of the stars used in the 

517 shapelet decomposition. 

518 - ``shapeletsStarUnNormalizedEMedian`` (`float`): The median 

519 un-normalized ellipticity 

520 (sqrt((starXX - starYY)**2.0 + (2.0*starXY)**2.0)) of the 

521 stars used in the shapelet decomposition. 

522 """ 

523 shapelets_summary = ExposureSummaryStats() 

524 results.non_gaussian_non_atmosphere_index_list = \ 

525 self.config.non_gaussian_non_atmosphere_index_list 

526 shapelets_coeffs = results.shapelet.getCoefficients() 

527 shapelets_coeffs_list = [float(coeff) for coeff in shapelets_coeffs] 

528 shapelets_summary.shapeletsCoeffs = shapelets_coeffs_list 

529 shapelets_used_cat = results.catalog[ 

530 (results.catalog[self.shapelets_base_name + "_used"]) 

531 & (~results.catalog[self.shapelets_base_name + "_flag"])].copy(deep=True) 

532 shapelets_summary.nShapeletsStar = len(shapelets_used_cat) 

533 

534 shapelets_star_xx = shapelets_used_cat[self.shapelets_base_name + "_xx"] 

535 shapelets_star_yy = shapelets_used_cat[self.shapelets_base_name + "_yy"] 

536 shapelets_star_xy = shapelets_used_cat[self.shapelets_base_name + "_xy"] 

537 

538 # Compute a shapelets score based on the power in the 

539 # non-atmospheric decomposition coefficients. 

540 total_power = np.sum(shapelets_coeffs**2.0) 

541 non_atm_power = np.sum( 

542 shapelets_coeffs[self.config.non_gaussian_non_atmosphere_index_list]**2.0 

543 ) 

544 shapelets_only_iq_score = np.where(total_power > 0, non_atm_power/total_power, 0.0) 

545 shapelets_summary.shapeletsOnlyIqScore = float(shapelets_only_iq_score) 

546 

547 self.log.info("n_shapelets_star = %d; shapelets_only_iq_score = %.5f", 

548 len(shapelets_used_cat), shapelets_only_iq_score) 

549 shapelets_star_e1 = (shapelets_star_xx - shapelets_star_yy)/(shapelets_star_xx + shapelets_star_yy) 

550 shapelets_star_e2 = 2*shapelets_star_xy/(shapelets_star_xx + shapelets_star_yy) 

551 

552 shapelets_star_e = np.sqrt(shapelets_star_e1**2.0 + shapelets_star_e2**2.0) 

553 shapelets_star_unnormalized_e = np.sqrt( 

554 (shapelets_star_xx - shapelets_star_yy)**2.0 + (2.0*shapelets_star_xy)**2.0) 

555 shapelets_star_e_median = np.median(shapelets_star_e) 

556 shapelets_star_unnormalized_e_median = np.median(shapelets_star_unnormalized_e) 

557 shapelets_summary.shapeletsStarEMedian = float(shapelets_star_e_median) 

558 shapelets_summary.shapeletsStarUnNormalizedEMedian = float(shapelets_star_unnormalized_e_median) 

559 results.shapelets_summary = shapelets_summary 

560 

561 return results 

562 

563 @staticmethod 

564 def compute_shapelets_iq_score(shapelets_results: Struct, catalog: SourceCatalog, 

565 centroid_offset_scale: float, shapelets_base_name: str, 

566 log: Logger | None = None) -> Struct: 

567 """Compute the image quality "score" that includes power from the 

568 median centroid offset in addition to that of the non-atmospheric 

569 coefficents. 

570 

571 The "shapelets_iq_score" is a dimensionless image quality score that 

572 includes power from the median centroid offset between those used in 

573 shapelet decomposition and those in the centroid slot of the provided 

574 source catalog, in addition to the non-atmospheric coefficents from 

575 the shapelet decomposition. The score spans the range [0.0, 1.0] with 

576 lower values indicating better image quality. When present (non-nan) 

577 the shapelets_iq_score should be considered authoritative over the 

578 shapelets_only_iq_score. 

579 

580 Parameters 

581 ---------- 

582 shapelets_results 

583 Results struct as returned by `run` method of 

584 ComputeRoughPsfShapeletsTask. To be updated in-place. 

585 catalog 

586 Catalog of detections to extract footprints from and fill output 

587 columns of. Its schema must be a superset of ``self.schema``. 

588 Must contain the columns added by ComputeRoughPsfShapeletsTask. 

589 centroid_offset_scale 

590 Value (in pixels) by which to scale the centroid offset to set 

591 and appropriate amount by which the centroid offsets can 

592 contribute to the computed power. 

593 shapelets_base_name 

594 Base name for the catalog columns added by 

595 ComputeRoughPsfShapeletsTask. 

596 log 

597 Logger to write to. 

598 

599 Returns 

600 ------- 

601 `lsst.pipe.base.Struct` 

602 Additional results struct components include: 

603 - ``shapelets_summary`` (`lsst.afw.image.ExposureSummaryStats`): 

604 An ExposureSummaryStats object with the following parameters 

605 initialized: 

606 - ``shapeletsIqScore`` (`float`): the image quality score which 

607 includes power from the centroid shift in addition to the 

608 non-atmospheric decomposition coefficients. 

609 - ``centroidDiffShapeletsVsSlotMedian`` (`float`): The median 

610 centroid difference 

611 (sqrt((slot_x - shapelet_x)**2 + (slot_y - shapelet_y)**2)) 

612 for sources used in the shapelet decomposition (pixels). 

613 """ 

614 if shapelets_base_name + "_used" not in catalog.schema.getNames(): 

615 raise RuntimeError("Catalog must contain columns added by ComputeRoughPsfShapeletsTask " 

616 f"(e.g.{shapelets_base_name}_used).") 

617 shapelets_summary = shapelets_results.shapelets_summary 

618 shapelets_coeffs = np.array(shapelets_summary.shapeletsCoeffs) 

619 shapelets_used_cat = catalog[(catalog[shapelets_base_name + "_used"]) 

620 & (~catalog[shapelets_base_name + "_flag"])].copy(deep=True) 

621 shapelets_used_finite_centroid_cat = shapelets_used_cat[( 

622 (np.isfinite(shapelets_used_cat["slot_Centroid_x"])) 

623 & (np.isfinite(shapelets_used_cat["slot_Centroid_y"])))] 

624 min_finite_centroid_sources = 3 

625 if len(shapelets_used_finite_centroid_cat) < min_finite_centroid_sources: 

626 if log is not None: 

627 log.warning("Not enough sources used in the shapelet decomposition have finite slot " 

628 "centroid measurements (%d < %d). Returning without setting shapelets_iq_score.", 

629 len(shapelets_used_finite_centroid_cat), min_finite_centroid_sources) 

630 shapelets_results.shapelets_summary = shapelets_summary 

631 return shapelets_results 

632 

633 centroid_diff_shapelets_vs_slot = np.sqrt( 

634 (shapelets_used_finite_centroid_cat["slot_Centroid_x"] 

635 - shapelets_used_finite_centroid_cat[shapelets_base_name + "_x"])**2.0 

636 + (shapelets_used_finite_centroid_cat["slot_Centroid_y"] 

637 - shapelets_used_finite_centroid_cat[shapelets_base_name + "_y"])**2.0 

638 ) 

639 centroid_diff_shapelets_vs_slot_median = np.nanmedian(centroid_diff_shapelets_vs_slot) 

640 shapelets_summary.centroidDiffShapeletsVsSlotMedian = float(centroid_diff_shapelets_vs_slot_median) 

641 

642 # Use a fixed number of pixels (one roughly estimated to be the 

643 # largest shift possible considering the maximum footprint size 

644 # imposed in ComputeRoughPsfShapeletsTask has been determined to 

645 # work well), to scale the centroid offset between the slot value 

646 # and that determined by ComputeRoughPsfShapeletsTask. 

647 centroid_diff_scaled = centroid_diff_shapelets_vs_slot/centroid_offset_scale 

648 centroid_diff_scaled_median = np.median(centroid_diff_scaled) 

649 

650 # Compute a shapelets score that includes power from the centroid shift 

651 # in addition to the non-atmospheric decomposition coefficients. 

652 total_power_plus_scaled = (np.sum(shapelets_coeffs**2.0) 

653 + centroid_diff_scaled_median**2.0) 

654 non_atm_power_plus_scaled = ( 

655 np.sum(shapelets_coeffs[shapelets_results.non_gaussian_non_atmosphere_index_list]**2.0) 

656 + centroid_diff_scaled_median**2.0) 

657 

658 shapelets_iq_score = np.where( 

659 total_power_plus_scaled > 0, non_atm_power_plus_scaled/total_power_plus_scaled, 0.0 

660 ) 

661 shapelets_summary.shapeletsIqScore = float(shapelets_iq_score) 

662 

663 if log is not None: 

664 log.info("shapelets_iq_score = %.5f; centroid_diff_shapelets_vs_slot_median = %.3f (pixels), " 

665 "centroid_diff_scaled_median (by %.1f pixels) = %.3f", 

666 shapelets_iq_score, centroid_diff_shapelets_vs_slot_median, 

667 centroid_offset_scale, centroid_diff_scaled_median) 

668 

669 shapelets_results.shapelets_summary = shapelets_summary 

670 return shapelets_results 

671 

672 def plot_selection( 

673 self, figure: matplotlib.figure.Figure, *, catalog: SourceCatalog, results: Struct 

674 ) -> None: 

675 """Create plots of the shape distribution space used to select stars. 

676 

677 Parameters 

678 ---------- 

679 figure 

680 Matplotlib figure to plot to. 

681 catalog 

682 Catalog of sources with columns populated by the `run` method (at 

683 least through the `select_stars` step). 

684 results 

685 Result struct returned by `run` or `select_stars`. 

686 """ 

687 from matplotlib.lines import Line2D 

688 

689 shape_data = np.zeros((len(catalog), 3), dtype=np.float64) 

690 radii = np.zeros(len(catalog), dtype=np.float64) 

691 for n, record in enumerate(catalog): 

692 if record[self._flag_key]: 

693 continue 

694 shape = record[self._shape_key] 

695 if self.config.logarithmic_shapes: 

696 shape = ellipses.SeparableConformalShearLogTraceRadius(shape) 

697 shape_data[n, :] = shape.getParameterVector() 

698 radii[n] = shape.getTraceRadius() 

699 used_mask = catalog[self._used_key] 

700 candidate_mask = np.logical_and(catalog[self._candidate_key], np.logical_not(used_mask)) 

701 measured_mask = np.logical_and( 

702 np.logical_not(catalog[self._flag_key]), np.logical_not(catalog[self._candidate_key]) 

703 ) 

704 # Set up the axes. 

705 axes = figure.subplot_mosaic( 

706 [ 

707 ["radius", "radius", "radius"], 

708 ["hist0", ".", "."], 

709 ["scatter01", "hist1", "."], 

710 ["scatter02", "scatter12", "hist2"], 

711 ], 

712 gridspec_kw=dict(bottom=0.2, top=1.0), 

713 ) 

714 axes["scatter01"].sharex(axes["hist0"]) 

715 axes["scatter02"].sharex(axes["hist0"]) 

716 axes["scatter12"].sharex(axes["hist1"]) 

717 axes["scatter12"].sharey(axes["scatter02"]) 

718 for tk in itertools.chain( 

719 axes["hist0"].get_xticklabels(), 

720 axes["scatter01"].get_xticklabels(), 

721 axes["hist1"].get_xticklabels(), 

722 axes["scatter12"].get_yticklabels(), 

723 ): 

724 tk.set_visible(False) 

725 # Move hist y axes to the outside. 

726 for ax in [axes["hist0"], axes["hist1"], axes["hist2"]]: 

727 ax.yaxis.set_label_position("right") 

728 ax.yaxis.tick_right() 

729 # Add labels to outer axes. 

730 names = ["Ixx", "Iyy", "Ixy"] if not self.config.logarithmic_shapes else ["𝜂1", "𝜂2", "ln(r)"] 

731 axes["scatter02"].set_xlabel(names[0]) 

732 axes["scatter12"].set_xlabel(names[1]) 

733 axes["hist2"].set_xlabel(names[2]) 

734 axes["scatter01"].set_ylabel(names[1]) 

735 axes["scatter02"].set_ylabel(names[2]) 

736 # Make the plots. 

737 mu = results.mean_shape.getParameterVector() 

738 sigma = np.diagonal(results.shape_covariance) ** 0.5 

739 lower_bounds = [max(mu[i] - 3 * sigma[i], min(shape_data[:, i])) for i in range(3)] 

740 upper_bounds = [min(mu[i] + 3 * sigma[i], max(shape_data[:, i])) for i in range(3)] 

741 grids = [np.linspace(lower_bounds[i], upper_bounds[i], 50) for i in range(3)] 

742 axes["radius"].axvline(results.radius_cut, color="k") 

743 for color, mask, alpha in [ 

744 ("grey", measured_mask, 0.5), 

745 ("blue", candidate_mask, 0.75), 

746 ("green", used_mask, 1.0), 

747 ]: 

748 axes["radius"].hist( 

749 radii[mask], 

750 color=color, 

751 alpha=alpha, 

752 bins=16, 

753 histtype="step", 

754 range=(radii.min(), 15.0), 

755 density=True, 

756 ) 

757 for i in range(3): 

758 axes[f"hist{i}"].hist( 

759 shape_data[mask, i], 

760 bins=16, 

761 range=(lower_bounds[i], upper_bounds[i]), 

762 density=True, 

763 color=color, 

764 histtype="step", 

765 alpha=alpha, 

766 ) 

767 for j in range(3): 

768 if (ax := axes.get(f"scatter{i}{j}")) is not None: 

769 ax.scatter( 

770 shape_data[mask, i], 

771 shape_data[mask, j], 

772 c=color, 

773 s=4, 

774 alpha=alpha, 

775 edgecolors=None, 

776 ) 

777 for i in range(3): 

778 axes[f"hist{i}"].plot( 

779 grids[i], scipy.stats.norm.pdf(grids[i], loc=mu[i], scale=sigma[i]), "k", alpha=0.5 

780 ) 

781 axes[f"hist{i}"].set_xlim(lower_bounds[i], upper_bounds[i]) 

782 for j in range(3): 

783 if (ax := axes.get(f"scatter{i}{j}")) is not None: 

784 sigma_ellipse = ellipses.Quadrupole( 

785 results.shape_covariance[i, i], 

786 results.shape_covariance[j, j], 

787 results.shape_covariance[i, j], 

788 ) 

789 for factor in [1, 2, 3]: 

790 self._draw_ellipse( 

791 ax, 

792 sigma_ellipse, 

793 x=mu[i], 

794 y=mu[j], 

795 scale=factor, 

796 fill=False, 

797 edgecolor="k", 

798 alpha=0.5, 

799 ) 

800 ax.set_xlim(lower_bounds[i], upper_bounds[i]) 

801 ax.set_ylim(lower_bounds[j], upper_bounds[j]) 

802 figure.legend( 

803 [ 

804 Line2D([], [], color="green", alpha=1.0), 

805 Line2D([], [], color="blue", alpha=0.75), 

806 Line2D([], [], color="gray", alpha=0.5), 

807 ], 

808 [ 

809 f"{self.shapelets_base_name}_used ({used_mask.sum()})", 

810 f"{self.shapelets_base_name}_candidate & ~{self.shapelets_base_name}_used " 

811 "({candidate_mask.sum()})", 

812 f"~{self.shapelets_base_name}_flag & ~{self.shapelets_base_name}_candidate " 

813 f"({measured_mask.sum()})", 

814 ], 

815 loc="lower center", 

816 ) 

817 return figure 

818 

819 def plot_shapelets( 

820 self, 

821 figure: matplotlib.figure.Figure, 

822 *, 

823 image: ImageF, 

824 catalog: SourceCatalog, 

825 results: Struct, 

826 n_stars: int = 3, 

827 stamp_size: float = 2.0,