Coverage for python/lsst/meas/extensions/psfex/psfexPsfDeterminer.py: 15%

1# This file is part of meas_extensions_psfex.

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.

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__all__ = (

23 "PsfexPsfDeterminerConfig",

24 "PsfexPsfDeterminerTask",

25 "PsfexNoStarsError",

26 "PsfexNoGoodStarsError",

27 "PsfexTooFewGoodStarsError",

28)

30import os

31import numpy as np

33import lsst.daf.base as dafBase

34import lsst.pex.config as pexConfig

35import lsst.pex.exceptions as pexExcept

36import lsst.geom as geom

37import lsst.afw.geom.ellipses as afwEll

38import lsst.afw.display as afwDisplay

39import lsst.afw.image as afwImage

40import lsst.afw.math as afwMath

41import lsst.meas.algorithms as measAlg

42import lsst.meas.algorithms.utils as maUtils

43import lsst.meas.extensions.psfex as psfex

44from lsst.pipe.base import AlgorithmError

47class PsfexNoStarsError(AlgorithmError):

48 """Raised if no stars are available for PSF determination."""

50 def __init__(self) -> None:

51 super().__init__("No psf candidates supplied.")

53 @property

54 def metadata(self) -> dict:

55 return {}

58class PsfexNoGoodStarsError(AlgorithmError):

59 """Raised if no "good" stars are available for PSF determination.

61 Parameters

62 ----------

63 num_available_stars : `int`

64 Number of available stars for PSF determination.

65 """

67 def __init__(self, num_available_stars) -> None:

68 self._num_available_stars = num_available_stars

69 super().__init__(f"No good psf candidates to pass to psfex out of {num_available_stars} available.")

71 @property

72 def metadata(self) -> dict:

73 return {

74 "num_available_stars": self._num_available_stars,

75 }

78class PsfexTooFewGoodStarsError(AlgorithmError):

79 """Raised if too few good stars are available for PSF determination.

81 Parameters

82 ----------

83 num_available_stars : `int`

84 Number of available stars for PSF determination.

85 num_good_stars : `int`

86 Number of good stars used for PSF determination.

87 poly_ndim_initial : `int`

88 Initial number of dependency parameters (dimensions) used in

89 polynomial fitting.

90 poly_ndim_final : `int`

91 Final number of dependency parameters (dimensions) set in the PSFEx

92 model after initializing the PSF structure.

93 """

95 def __init__(

96 self,

97 num_available_stars,

98 num_good_stars,

99 poly_ndim_initial,

100 poly_ndim_final,

101 ) -> None:

102 self._num_available_stars = num_available_stars

103 self._num_good_stars = num_good_stars

104 self._poly_ndim_initial = poly_ndim_initial

105 self._poly_ndim_final = poly_ndim_final

106 super().__init__(

107 f"Failed to determine psfex psf: too few good stars ({num_good_stars}) out of "

108 f"{num_available_stars} available. To accommodate this low count, the polynomial dimension was "

109 f"lowered from {poly_ndim_initial} to {poly_ndim_final}, which is not handled by the Science "

110 "Pipelines code."

111 )

112

113 @property

114 def metadata(self) -> dict:

115 return {

116 "num_available_stars": self._num_available_stars,

117 "num_good_stars": self._num_good_stars,

118 "poly_ndim_initial": self._poly_ndim_initial,

119 "poly_ndim_final": self._poly_ndim_final,

120 }

121

122

123class PsfexPsfDeterminerConfig(measAlg.BasePsfDeterminerConfig):

124 spatialOrder = pexConfig.Field[int](

125 doc="specify spatial order for PSF kernel creation",

126 default=2,

127 check=lambda x: x >= 1,

128 )

129 sizeCellX = pexConfig.Field[int](

130 doc="size of cell used to determine PSF (pixels, column direction)",

131 default=256,

132 # minValue = 10,

133 check=lambda x: x >= 10,

134 )

135 sizeCellY = pexConfig.Field[int](

136 doc="size of cell used to determine PSF (pixels, row direction)",

137 default=sizeCellX.default,

138 # minValue = 10,

139 check=lambda x: x >= 10,

140 )

141 samplingSize = pexConfig.Field[float](

142 doc="Resolution of the internal PSF model relative to the pixel size; "

143 "e.g. 0.5 is equal to 2x oversampling",

144 default=0.5,

145 )

146 badMaskBits = pexConfig.ListField[str](

147 doc="List of mask bits which cause a source to be rejected as bad "

148 "N.b. INTRP is used specially in PsfCandidateSet; it means \"Contaminated by neighbour\"",

149 default=["INTRP", "SAT"],

150 )

151 psfexBasis = pexConfig.ChoiceField[str](

152 doc="BASIS value given to psfex. PIXEL_AUTO will use the requested samplingSize only if "

153 "the FWHM < 3 pixels. Otherwise, it will use samplingSize=1. PIXEL will always use the "

154 "requested samplingSize",

155 allowed={

156 "PIXEL": "Always use requested samplingSize",

157 "PIXEL_AUTO": "Only use requested samplingSize when FWHM < 3",

158 },

159 default='PIXEL_AUTO',

160 optional=False,

161 )

162 tolerance = pexConfig.Field[float](

163 doc="tolerance of spatial fitting",

164 default=1e-2,

165 )

166 lam = pexConfig.Field[float](

167 doc="floor for variance is lam*data",

168 default=0.05,

169 )

170 reducedChi2ForPsfCandidates = pexConfig.Field[float](

171 doc="for psf candidate evaluation",

172 default=2.0,

173 )

174 spatialReject = pexConfig.Field[float](

175 doc="Rejection threshold (stdev) for candidates based on spatial fit",

176 default=3.0,

177 )

178 recentroid = pexConfig.Field[bool](

179 doc="Should PSFEX be permitted to recentroid PSF candidates?",

180 default=False,

181 )

182 photometricFluxField = pexConfig.Field[str](

183 doc="Flux field to use for photometric normalization. This overrides the "

184 "``PHOTFLUX_KEY`` field for psfex. The associated flux error is "

185 "derived by appending ``Err`` to this field.",

186 default="base_CircularApertureFlux_9_0_instFlux",

187 )

188 rejectFlaggedCentroids = pexConfig.Field[bool](

189 doc="Should flagged centroids be rejected? Note that setting this to False "

190 "additionally ignores aperture flux flags since those are tied together.",

191 default=True,

192 )

193

194 def setDefaults(self):

195 super().setDefaults()

196 self.stampSize = 41

197

198 def validate(self):

199 super().validate()

200 fluxFieldEnding = self.photometricFluxField.rpartition("_")[-1]

201 if "_" not in self.photometricFluxField or not (

202 fluxFieldEnding == "flux" or fluxFieldEnding.endswith("Flux")

203 ):

204 raise ValueError(

205 "Expected `photometricFluxField` to end with '_flux' "

206 f"or '_*Flux', but got '{self.photometricFluxField}'"

207 )

208

209

210class PsfexPsfDeterminerTask(measAlg.BasePsfDeterminerTask):

211 ConfigClass = PsfexPsfDeterminerConfig

212

213 def determinePsf(self, exposure, psfCandidateList, metadata=None, flagKey=None):

214 """Determine a PSFEX PSF model for an exposure given a list of PSF

215 candidates.

216

217 Parameters

218 ----------

219 exposure: `lsst.afw.image.Exposure`

220 Exposure containing the PSF candidates.

221 psfCandidateList: iterable of `lsst.meas.algorithms.PsfCandidate`

222 Sequence of PSF candidates typically obtained by detecting sources

223 and then running them through a star selector.

224 metadata: metadata, optional

225 A home for interesting tidbits of information.

226 flagKey: `lsst.afw.table.Key`, optional

227 Schema key used to mark sources actually used in PSF determination.

228

229 Returns

230 -------

231 psf: `lsst.meas.extensions.psfex.PsfexPsf`

232 The determined PSF.

233 """

234 psfCandidateList = self.downsampleCandidates(psfCandidateList)

235

236 import lsstDebug

237 display = lsstDebug.Info(__name__).display

238 displayExposure = display and \

239 lsstDebug.Info(__name__).displayExposure # display the Exposure + spatialCells

240 displayPsfComponents = display and \

241 lsstDebug.Info(__name__).displayPsfComponents # show the basis functions

242 showBadCandidates = display and \

243 lsstDebug.Info(__name__).showBadCandidates # Include bad candidates (meaningless, methinks)

244 displayResiduals = display and \

245 lsstDebug.Info(__name__).displayResiduals # show residuals

246 displayPsfMosaic = display and \

247 lsstDebug.Info(__name__).displayPsfMosaic # show mosaic of reconstructed PSF(x,y)

248 normalizeResiduals = lsstDebug.Info(__name__).normalizeResiduals

249 afwDisplay.setDefaultMaskTransparency(75)

250 # Normalise residuals by object amplitude

251

252 mi = exposure.getMaskedImage()

253

254 nCand = len(psfCandidateList)

255 if nCand == 0:

256 raise PsfexNoStarsError()

257 #

258 # How big should our PSF models be?

259 #

260 if display: # only needed for debug plots

261 # construct and populate a spatial cell set

262 bbox = mi.getBBox(afwImage.PARENT)

263 psfCellSet = afwMath.SpatialCellSet(bbox, self.config.sizeCellX, self.config.sizeCellY)

264 else:

265 psfCellSet = None

266

267 sizes = np.empty(nCand)

268 for i, psfCandidate in enumerate(psfCandidateList):

269 try:

270 if psfCellSet:

271 psfCellSet.insertCandidate(psfCandidate)

272 except Exception as e:

273 self.log.error("Skipping PSF candidate %d of %d: %s", i, len(psfCandidateList), e)

274 continue

275

276 source = psfCandidate.getSource()

277 quad = afwEll.Quadrupole(source.getIxx(), source.getIyy(), source.getIxy())

278 rmsSize = quad.getTraceRadius()

279 sizes[i] = rmsSize

280

281 pixKernelSize = self.config.stampSize

282 actualKernelSize = int(2*np.floor(0.5*pixKernelSize/self.config.samplingSize) + 1)

283

284 if display:

285 rms = np.median(sizes)

286 self.log.debug("Median PSF RMS size=%.2f pixels (\"FWHM\"=%.2f)",

287 rms, 2*np.sqrt(2*np.log(2))*rms)

288

289 self.log.trace("Psfex Kernel size=%.2f, Image Kernel Size=%.2f", actualKernelSize, pixKernelSize)

290

291 # -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- BEGIN PSFEX

292 #

293 # Insert the good candidates into the set

294 #

295 defaultsFile = os.path.join(os.environ["MEAS_EXTENSIONS_PSFEX_DIR"], "config", "default-lsst.psfex")

296 args_md = dafBase.PropertySet()

297 args_md.set("BASIS_TYPE", str(self.config.psfexBasis))

298 args_md.set("PSFVAR_DEGREES", str(self.config.spatialOrder))

299 args_md.set("PSF_SIZE", str(actualKernelSize))

300 args_md.set("PSF_SAMPLING", str(self.config.samplingSize))

301 args_md.set("PHOTFLUX_KEY", str(self.config.photometricFluxField))

302 args_md.set("PHOTFLUXERR_KEY", str(self.config.photometricFluxField) + "Err")

303 prefs = psfex.Prefs(defaultsFile, args_md)

304 prefs.setCommandLine([])

305 prefs.addCatalog("psfexPsfDeterminer")

306

307 prefs.use()

308 principalComponentExclusionFlag = bool(bool(psfex.Context.REMOVEHIDDEN)

309 if False else psfex.Context.KEEPHIDDEN)

310 context = psfex.Context(prefs.getContextName(), prefs.getContextGroup(),

311 prefs.getGroupDeg(), principalComponentExclusionFlag)

312 psfSet = psfex.Set(context)

313 psfSet.setVigSize(pixKernelSize, pixKernelSize)

314 psfSet.setFwhm(2*np.sqrt(2*np.log(2))*np.median(sizes))

315 psfSet.setRecentroid(self.config.recentroid)

316

317 catindex, ext = 0, 0

318 backnoise2 = afwMath.makeStatistics(mi.getImage(), afwMath.VARIANCECLIP).getValue()

319 ccd = exposure.getDetector()

320 if ccd:

321 gain = np.mean(np.array([a.getGain() for a in ccd]))

322 else:

323 gain = 1.0

324 self.log.warning("Setting gain to %g", gain)

325

326 contextvalp = []

327 for i, key in enumerate(context.getName()):

328 if key[0] == ':':

329 try:

330 contextvalp.append(exposure.getMetadata().getScalar(key[1:]))

331 except KeyError as e:

332 raise RuntimeError("%s parameter not found in the header of %s" %

333 (key[1:], prefs.getContextName())) from e

334 else:

335 try:

336 contextvalp.append(np.array([psfCandidateList[_].getSource().get(key)

337 for _ in range(nCand)]))

338 except KeyError as e:

339 raise RuntimeError("%s parameter not found" % (key,)) from e

340 psfSet.setContextname(i, key)

341

342 if display:

343 frame = 0

344 if displayExposure:

345 disp = afwDisplay.Display(frame=frame)

346 disp.mtv(exposure, title="psf determination")

347

348 badBits = mi.getMask().getPlaneBitMask(self.config.badMaskBits)

349 fluxName = prefs.getPhotfluxRkey()

350 fluxFlagName = fluxName.rpartition("_")[0] + "_flag"

351

352 xpos, ypos = [], []

353 for i, psfCandidate in enumerate(psfCandidateList):

354 source = psfCandidate.getSource()

355

356 # skip sources with bad centroids

357 xc, yc = source.getX(), source.getY()

358 if not np.isfinite(xc) or not np.isfinite(yc):

359 continue

360 if self.config.rejectFlaggedCentroids:

361 # centroids might be finite but still failing

362 if source.get("slot_Centroid_flag"):

363 continue

364 # skip flagged sources

365 if source.get(fluxFlagName):

366 continue

367 # skip nonfinite and negative sources

368 flux = source.get(fluxName)

369 if flux < 0 or not np.isfinite(flux):

370 continue

371

372 try:

373 pstamp = psfCandidate.getMaskedImage(pixKernelSize, pixKernelSize).clone()

374 except pexExcept.LengthError:

375 self.log.warning("Could not get stamp image for psfCandidate: %s with kernel size: %s",

376 psfCandidate, pixKernelSize)

377 continue

378

379 # From this point, we're configuring the "sample" (PSFEx's version

380 # of a PSF candidate).

381 # Having created the sample, we must proceed to configure it, and

382 # then fini (finalize), or it will be malformed.

383 try:

384 sample = psfSet.newSample()

385 sample.setCatindex(catindex)

386 sample.setExtindex(ext)

387 sample.setObjindex(i)

388

389 imArray = pstamp.getImage().getArray()

390 imArray[np.where(np.bitwise_and(pstamp.getMask().getArray(), badBits))] = \

391 -2*psfex.BIG

392 sample.setVig(imArray)

393

394 sample.setNorm(flux)

395 sample.setBacknoise2(backnoise2)

396 sample.setGain(gain)

397 sample.setX(xc)

398 sample.setY(yc)

399 sample.setFluxrad(sizes[i])

400

401 for j in range(psfSet.getNcontext()):

402 sample.setContext(j, float(contextvalp[j][i]))

403 except Exception as e:

404 self.log.error("Exception when processing sample at (%f,%f): %s", xc, yc, e)

405 continue

406 else:

407 psfSet.finiSample(sample)

408

409 xpos.append(xc) # for QA

410 ypos.append(yc)

411

412 if displayExposure:

413 with disp.Buffering():

414 disp.dot("o", xc, yc, ctype=afwDisplay.CYAN, size=4)

415

416 if psfSet.getNsample() == 0:

417 raise PsfexNoGoodStarsError(num_available_stars=nCand)

418

419 # ---- Update min and max and then the scaling

420 for i in range(psfSet.getNcontext()):

421 cmin = contextvalp[i].min()

422 cmax = contextvalp[i].max()

423 psfSet.setContextScale(i, cmax - cmin)

424 psfSet.setContextOffset(i, (cmin + cmax)/2.0)

425

426 # Don't waste memory!

427 psfSet.trimMemory()

428

429 # =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- END PSFEX

430 #

431 # Do a PSFEX decomposition of those PSF candidates

432 #

433 fields = []

434 field = psfex.Field("Unknown")

435 field.addExt(exposure.getWcs(), exposure.getWidth(), exposure.getHeight(), psfSet.getNsample())

436 field.finalize()

437

438 fields.append(field)

439

440 sets = []

441 sets.append(psfSet)

442

443 psfex.makeit(fields, sets)

444 psfs = field.getPsfs()

445

446 # Flag which objects were actually used in psfex by

447 good_indices = []

448 for i in range(sets[0].getNsample()):

449 index = sets[0].getSample(i).getObjindex()

450 if index > -1:

451 good_indices.append(index)

452

453 if flagKey is not None:

454 for i, psfCandidate in enumerate(psfCandidateList):

455 source = psfCandidate.getSource()

456 if i in good_indices:

457 source.set(flagKey, True)

458

459 xpos = np.array(xpos)

460 ypos = np.array(ypos)

461 numGoodStars = len(good_indices)

462 avgX, avgY = np.mean(xpos), np.mean(ypos)

463

464 psf = psfex.PsfexPsf(psfs[0], geom.Point2D(avgX, avgY))

465

466 # If there are too few stars, the PSFEx psf model will reduce the order

467 # to 0, which the Science Pipelines code cannot handle (see

468 # https://github.com/lsst/meas_extensions_psfex/blob/f0d5218b5446faf5e39edc30e31d2e6f673ef294/src/PsfexPsf.cc#L118

469 # ).

470 if (ndim := psf.getNdim()) == 0:

471 raise PsfexTooFewGoodStarsError(

472 num_available_stars=nCand,

473 num_good_stars=numGoodStars,

474 poly_ndim_initial=psfSet.getNcontext(),

475 poly_ndim_final=ndim,

476 )

477

478 #

479 # Display code for debugging

480 #

481 if display:

482 assert psfCellSet is not None

483

484 if displayExposure:

485 maUtils.showPsfSpatialCells(exposure, psfCellSet, showChi2=True,

486 symb="o", ctype=afwDisplay.YELLOW, ctypeBad=afwDisplay.RED,

487 size=8, display=disp)

488 if displayResiduals:

489 disp4 = afwDisplay.Display(frame=4)

490 maUtils.showPsfCandidates(exposure, psfCellSet, psf=psf, display=disp4,

491 normalize=normalizeResiduals,

492 showBadCandidates=showBadCandidates)

493 if displayPsfComponents:

494 disp6 = afwDisplay.Display(frame=6)

495 maUtils.showPsf(psf, display=disp6)

496 if displayPsfMosaic:

497 disp7 = afwDisplay.Display(frame=7)

498 maUtils.showPsfMosaic(exposure, psf, display=disp7, showFwhm=True)

499 disp.scale('linear', 0, 1)

500 #

501 # Generate some QA information

502 #

503 # Count PSF stars

504 #

505 if metadata is not None:

506 metadata["spatialFitChi2"] = np.nan

507 metadata["numAvailStars"] = nCand

508 metadata["numGoodStars"] = numGoodStars

509 metadata["avgX"] = avgX

510 metadata["avgY"] = avgY

511

512 return psf, psfCellSet

513

514

515measAlg.psfDeterminerRegistry.register("psfex", PsfexPsfDeterminerTask)

Coverage for python / lsst / meas / extensions / psfex / psfexPsfDeterminer.py: 15%

244 statements