Coverage for python / lsst / meas / base / tests.py: 13%

1# This file is part of meas_base. 

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 

22import warnings 

23 

24import numpy as np 

25 

26import lsst.geom 

27import lsst.afw.table 

28import lsst.afw.image 

29import lsst.afw.detection 

30import lsst.afw.cameraGeom.testUtils 

31import lsst.afw.geom 

32import lsst.afw.coord 

33import lsst.daf.base 

34import lsst.pex.exceptions 

35 

36from .sfm import SingleFrameMeasurementTask 

37from .forcedMeasurement import ForcedMeasurementTask 

38from ._measBaseLib import CentroidResultKey 

39 

40__all__ = ("BlendContext", "TestDataset", "AlgorithmTestCase", "TransformTestCase", 

41 "SingleFramePluginTransformSetupHelper", "ForcedPluginTransformSetupHelper", 

42 "FluxTransformTestCase", "CentroidTransformTestCase") 

43 

44 

45class BlendContext: 

46 """Context manager which adds multiple overlapping sources and a parent. 

47 

48 Notes 

49 ----- 

50 This is used as the return value for `TestDataset.addBlend`, and this is 

51 the only way it should be used. 

52 """ 

53 

54 def __init__(self, owner): 

55 self.owner = owner 

56 self.parentRecord = self.owner.catalog.addNew() 

57 self.parentImage = lsst.afw.image.ImageF(self.owner.exposure.getBBox()) 

58 self.children = [] 

59 

60 def __enter__(self): 

61 # BlendContext is its own context manager, so we just return self. 

62 return self 

63 

64 def addChild(self, instFlux, centroid, shape=None): 

65 """Add a child to the blend; return corresponding truth catalog record. 

66 

67 instFlux : `float` 

68 Total instFlux of the source to be added. 

69 centroid : `lsst.geom.Point2D` 

70 Position of the source to be added. 

71 shape : `lsst.afw.geom.Quadrupole` 

72 Second moments of the source before PSF convolution. Note that 

73 the truth catalog records post-convolution moments) 

74 """ 

75 record, image = self.owner.addSource(instFlux, centroid, shape) 

76 record.set(self.owner.keys["parent"], self.parentRecord.getId()) 

77 self.parentImage += image 

78 self.children.append((record, image)) 

79 return record 

80 

81 def __exit__(self, type_, value, tb): 

82 # We're not using the context manager for any kind of exception safety 

83 # or guarantees; we just want the nice "with" statement syntax. 

84 

85 if type_ is not None: 

86 # exception was raised; just skip all this and let it propagate 

87 return 

88 

89 # On exit, compute and set the truth values for the parent object. 

90 self.parentRecord.set(self.owner.keys["nChild"], len(self.children)) 

91 # Compute instFlux from sum of component fluxes 

92 instFlux = 0.0 

93 for record, image in self.children: 

94 instFlux += record.get(self.owner.keys["instFlux"]) 

95 self.parentRecord.set(self.owner.keys["instFlux"], instFlux) 

96 # Compute centroid from instFlux-weighted mean of component centroids 

97 x = 0.0 

98 y = 0.0 

99 for record, image in self.children: 

100 w = record.get(self.owner.keys["instFlux"])/instFlux 

101 x += record.get(self.owner.keys["centroid"].getX())*w 

102 y += record.get(self.owner.keys["centroid"].getY())*w 

103 self.parentRecord.set(self.owner.keys["centroid"], lsst.geom.Point2D(x, y)) 

104 # Compute shape from instFlux-weighted mean of offset component shapes 

105 xx = 0.0 

106 yy = 0.0 

107 xy = 0.0 

108 for record, image in self.children: 

109 w = record.get(self.owner.keys["instFlux"])/instFlux 

110 dx = record.get(self.owner.keys["centroid"].getX()) - x 

111 dy = record.get(self.owner.keys["centroid"].getY()) - y 

112 xx += (record.get(self.owner.keys["shape"].getIxx()) + dx**2)*w 

113 yy += (record.get(self.owner.keys["shape"].getIyy()) + dy**2)*w 

114 xy += (record.get(self.owner.keys["shape"].getIxy()) + dx*dy)*w 

115 self.parentRecord.set(self.owner.keys["shape"], lsst.afw.geom.Quadrupole(xx, yy, xy)) 

116 # Run detection on the parent image to get the parent Footprint. 

117 self.owner._installFootprint(self.parentRecord, self.parentImage) 

118 # Create perfect HeavyFootprints for all children; these will need to 

119 # be modified later to account for the noise we'll add to the image. 

120 deblend = lsst.afw.image.MaskedImageF(self.owner.exposure.maskedImage, True) 

121 for record, image in self.children: 

122 deblend.image.array[:, :] = image.array 

123 heavyFootprint = lsst.afw.detection.HeavyFootprintF(self.parentRecord.getFootprint(), deblend) 

124 record.setFootprint(heavyFootprint) 

125 

126 

127class TestDataset: 

128 """A simulated dataset consisuting of test image and truth catalog. 

129 

130 TestDataset creates an idealized image made of pure Gaussians (including a 

131 Gaussian PSF), with simple noise and idealized Footprints/HeavyFootprints 

132 that simulated the outputs of detection and deblending. Multiple noise 

133 realizations can be created from the same underlying sources, allowing 

134 uncertainty estimates to be verified via Monte Carlo. 

135 

136 Parameters 

137 ---------- 

138 bbox : `lsst.geom.Box2I` or `lsst.geom.Box2D` 

139 Bounding box of the test image. 

140 threshold : `float` 

141 Threshold absolute value used to determine footprints for 

142 simulated sources. This thresholding will be applied before noise is 

143 actually added to images (or before the noise level is even known), so 

144 this will necessarily produce somewhat artificial footprints. 

145 exposure : `lsst.afw.image.ExposureF` 

146 The image to which test sources should be added. Ownership should 

147 be considered transferred from the caller to the TestDataset. 

148 Must have a Gaussian PSF for truth catalog shapes to be exact. 

149 **kwds 

150 Keyword arguments forwarded to makeEmptyExposure if exposure is `None`. 

151 

152 Notes 

153 ----- 

154 Typical usage: 

155 

156 .. code-block: py 

157 

158 bbox = lsst.geom.Box2I(lsst.geom.Point2I(0,0), lsst.geom.Point2I(100, 

159 100)) 

160 dataset = TestDataset(bbox) 

161 dataset.addSource(instFlux=1E5, centroid=lsst.geom.Point2D(25, 26)) 

162 dataset.addSource(instFlux=2E5, centroid=lsst.geom.Point2D(75, 24), 

163 shape=lsst.afw.geom.Quadrupole(8, 7, 2)) 

164 with dataset.addBlend() as family: 

165 family.addChild(instFlux=2E5, centroid=lsst.geom.Point2D(50, 72)) 

166 family.addChild(instFlux=1.5E5, centroid=lsst.geom.Point2D(51, 74)) 

167 exposure, catalog = dataset.realize(noise=100.0, 

168 schema=TestDataset.makeMinimalSchema()) 

169 """ 

170 

171 def __init__(self, bbox, threshold=10.0, exposure=None, **kwds): 

172 if exposure is None: 

173 exposure = self.makeEmptyExposure(bbox, **kwds) 

174 self.threshold = lsst.afw.detection.Threshold(threshold, lsst.afw.detection.Threshold.VALUE) 

175 self.exposure = exposure 

176 self.psfShape = self.exposure.getPsf().computeShape(bbox.getCenter()) 

177 self.schema = self.makeMinimalSchema() 

178 self.catalog = lsst.afw.table.SourceCatalog(self.schema) 

179 

180 @classmethod 

181 def makeMinimalSchema(cls): 

182 """Return the minimal schema needed to hold truth catalog fields. 

183 

184 Notes 

185 ----- 

186 When `TestDataset.realize` is called, the schema must include at least 

187 these fields. Usually it will include additional fields for 

188 measurement algorithm outputs, allowing the same catalog to be used 

189 for both truth values (the fields from the minimal schema) and the 

190 measurements. 

191 """ 

192 if not hasattr(cls, "_schema"): 

193 schema = lsst.afw.table.SourceTable.makeMinimalSchema() 

194 cls.keys = {} 

195 cls.keys["coordErr"] = lsst.afw.table.CoordKey.addErrorFields(schema) 

196 cls.keys["parent"] = schema.find("parent").key 

197 cls.keys["nChild"] = schema.addField("deblend_nChild", type=np.int32) 

198 cls.keys["instFlux"] = schema.addField("truth_instFlux", type=np.float64, 

199 doc="true instFlux", units="count") 

200 cls.keys["instFluxErr"] = schema.addField("truth_instFluxErr", type=np.float64, 

201 doc="true instFluxErr", units="count") 

202 cls.keys["centroid"] = lsst.afw.table.Point2DKey.addFields( 

203 schema, "truth", "true simulated centroid", "pixel" 

204 ) 

205 cls.keys["centroid_sigma"] = lsst.afw.table.CovarianceMatrix2fKey.addFields( 

206 schema, "truth", ['x', 'y'], "pixel" 

207 ) 

208 cls.keys["centroid_flag"] = schema.addField("truth_flag", type="Flag", 

209 doc="set if the object is a star") 

210 cls.keys["shape"] = lsst.afw.table.QuadrupoleKey.addFields( 

211 schema, "truth", "true shape after PSF convolution", lsst.afw.table.CoordinateType.PIXEL 

212 ) 

213 cls.keys["isStar"] = schema.addField("truth_isStar", type="Flag", 

214 doc="set if the object is a star") 

215 cls.keys["negative"] = schema.addField("is_negative", type="Flag", 

216 doc="set if source was detected as significantly negative" 

217 ) 

218 schema.getAliasMap().set("slot_Shape", "truth") 

219 schema.getAliasMap().set("slot_Centroid", "truth") 

220 schema.getAliasMap().set("slot_ModelFlux", "truth") 

221 cls._schema = schema 

222 schema = lsst.afw.table.Schema(cls._schema) 

223 schema.disconnectAliases() 

224 return schema 

225 

226 @staticmethod 

227 def makePerturbedWcs(oldWcs, minScaleFactor=1.2, maxScaleFactor=1.5, 

228 minRotation=None, maxRotation=None, 

229 minRefShift=None, maxRefShift=None, 

230 minPixShift=2.0, maxPixShift=4.0, randomSeed=1): 

231 """Return a perturbed version of the input WCS. 

232 

233 Create a new undistorted TAN WCS that is similar but not identical to 

234 another, with random scaling, rotation, and offset (in both pixel 

235 position and reference position). 

236 

237 Parameters 

238 ---------- 

239 oldWcs : `lsst.afw.geom.SkyWcs` 

240 The input WCS. 

241 minScaleFactor : `float` 

242 Minimum scale factor to apply to the input WCS. 

243 maxScaleFactor : `float` 

244 Maximum scale factor to apply to the input WCS. 

245 minRotation : `lsst.geom.Angle` or `None` 

246 Minimum rotation to apply to the input WCS. If `None`, defaults to 

247 30 degrees. 

248 maxRotation : `lsst.geom.Angle` or `None` 

249 Minimum rotation to apply to the input WCS. If `None`, defaults to 

250 60 degrees. 

251 minRefShift : `lsst.geom.Angle` or `None` 

252 Miniumum shift to apply to the input WCS reference value. If 

253 `None`, defaults to 0.5 arcsec. 

254 maxRefShift : `lsst.geom.Angle` or `None` 

255 Miniumum shift to apply to the input WCS reference value. If 

256 `None`, defaults to 1.0 arcsec. 

257 minPixShift : `float` 

258 Minimum shift to apply to the input WCS reference pixel. 

259 maxPixShift : `float` 

260 Maximum shift to apply to the input WCS reference pixel. 

261 randomSeed : `int` 

262 Random seed. 

263 

264 Returns 

265 ------- 

266 newWcs : `lsst.afw.geom.SkyWcs` 

267 A perturbed version of the input WCS. 

268 

269 Notes 

270 ----- 

271 The maximum and minimum arguments are interpreted as absolute values 

272 for a split range that covers both positive and negative values (as 

273 this method is used in testing, it is typically most important to 

274 avoid perturbations near zero). Scale factors are treated somewhat 

275 differently: the actual scale factor is chosen between 

276 ``minScaleFactor`` and ``maxScaleFactor`` OR (``1/maxScaleFactor``) 

277 and (``1/minScaleFactor``). 

278 

279 The default range for rotation is 30-60 degrees, and the default range 

280 for reference shift is 0.5-1.0 arcseconds (these cannot be safely 

281 included directly as default values because Angle objects are 

282 mutable). 

283 

284 The random number generator is primed with the seed given. If 

285 `None`, a seed is automatically chosen. 

286 """ 

287 random_state = np.random.RandomState(randomSeed) 

288 if minRotation is None: 

289 minRotation = 30.0*lsst.geom.degrees 

290 if maxRotation is None: 

291 maxRotation = 60.0*lsst.geom.degrees 

292 if minRefShift is None: 

293 minRefShift = 0.5*lsst.geom.arcseconds 

294 if maxRefShift is None: 

295 maxRefShift = 1.0*lsst.geom.arcseconds 

296 

297 def splitRandom(min1, max1, min2=None, max2=None): 

298 if min2 is None: 

299 min2 = -max1 

300 if max2 is None: 

301 max2 = -min1 

302 if random_state.uniform() > 0.5: 

303 return float(random_state.uniform(min1, max1)) 

304 else: 

305 return float(random_state.uniform(min2, max2)) 

306 # Generate random perturbations 

307 scaleFactor = splitRandom(minScaleFactor, maxScaleFactor, 1.0/maxScaleFactor, 1.0/minScaleFactor) 

308 rotation = splitRandom(minRotation.asRadians(), maxRotation.asRadians())*lsst.geom.radians 

309 refShiftRa = splitRandom(minRefShift.asRadians(), maxRefShift.asRadians())*lsst.geom.radians 

310 refShiftDec = splitRandom(minRefShift.asRadians(), maxRefShift.asRadians())*lsst.geom.radians 

311 pixShiftX = splitRandom(minPixShift, maxPixShift) 

312 pixShiftY = splitRandom(minPixShift, maxPixShift) 

313 # Compute new CD matrix 

314 oldTransform = lsst.geom.LinearTransform(oldWcs.getCdMatrix()) 

315 rTransform = lsst.geom.LinearTransform.makeRotation(rotation) 

316 sTransform = lsst.geom.LinearTransform.makeScaling(scaleFactor) 

317 newTransform = oldTransform*rTransform*sTransform 

318 matrix = newTransform.getMatrix() 

319 # Compute new coordinate reference pixel (CRVAL) 

320 oldSkyOrigin = oldWcs.getSkyOrigin() 

321 newSkyOrigin = lsst.geom.SpherePoint(oldSkyOrigin.getRa() + refShiftRa, 

322 oldSkyOrigin.getDec() + refShiftDec) 

323 # Compute new pixel reference pixel (CRPIX) 

324 oldPixOrigin = oldWcs.getPixelOrigin() 

325 newPixOrigin = lsst.geom.Point2D(oldPixOrigin.getX() + pixShiftX, 

326 oldPixOrigin.getY() + pixShiftY) 

327 return lsst.afw.geom.makeSkyWcs(crpix=newPixOrigin, crval=newSkyOrigin, cdMatrix=matrix) 

328 

329 @staticmethod 

330 def makeEmptyExposure(bbox, wcs=None, crval=None, cdelt=None, psfSigma=2.0, psfDim=17, calibration=4, 

331 visitId=1234, mjd=60000.0, detector=1): 

332 """Create an Exposure, with a PhotoCalib, Wcs, and Psf, but no pixel values. 

333 

334 Parameters 

335 ---------- 

336 bbox : `lsst.geom.Box2I` or `lsst.geom.Box2D` 

337 Bounding box of the image in image coordinates. 

338 wcs : `lsst.afw.geom.SkyWcs`, optional 

339 New WCS for the exposure (created from CRVAL and CDELT if `None`). 

340 crval : `lsst.afw.geom.SpherePoint`, optional 

341 ICRS center of the TAN WCS attached to the image. If `None`, (45 

342 degrees, 45 degrees) is assumed. 

343 cdelt : `lsst.geom.Angle`, optional 

344 Pixel scale of the image. If `None`, 0.2 arcsec is assumed. 

345 psfSigma : `float`, optional 

346 Radius (sigma) of the Gaussian PSF attached to the image 

347 psfDim : `int`, optional 

348 Width and height of the image's Gaussian PSF attached to the image 

349 calibration : `float`, optional 

350 The spatially-constant calibration (in nJy/count) to set the 

351 PhotoCalib of the exposure. 

352 visitId : `int`, optional 

353 Visit id to store in VisitInfo. 

354 mjd : `float`, optional 

355 Modified Julian Date of this exposure to store in VisitInfo. 

356 detector: `int`, optional 

357 Detector id to assign to the attached Detector object. 

358 

359 Returns 

360 ------- 

361 exposure : `lsst.age.image.ExposureF` 

362 An empty image. 

363 """ 

364 if wcs is None: 

365 if crval is None: 

366 crval = lsst.geom.SpherePoint(45.0, 45.0, lsst.geom.degrees) 

367 if cdelt is None: 

368 cdelt = 0.2*lsst.geom.arcseconds 

369 crpix = lsst.geom.Box2D(bbox).getCenter() 

370 wcs = lsst.afw.geom.makeSkyWcs(crpix=crpix, crval=crval, 

371 cdMatrix=lsst.afw.geom.makeCdMatrix(scale=cdelt)) 

372 exposure = lsst.afw.image.ExposureF(bbox) 

373 psf = lsst.afw.detection.GaussianPsf(psfDim, psfDim, psfSigma) 

374 photoCalib = lsst.afw.image.PhotoCalib(calibration) 

375 visitInfo = lsst.afw.image.VisitInfo(id=visitId, 

376 exposureTime=30.0, 

377 date=lsst.daf.base.DateTime(mjd), 

378 observatory=lsst.afw.coord.Observatory(11.1*lsst.geom.degrees, 

379 22.2*lsst.geom.degrees, 

380 0.333), 

381 hasSimulatedContent=True) 

382 exposure.setWcs(wcs) 

383 exposure.setPsf(psf) 

384 exposure.setPhotoCalib(photoCalib) 

385 exposure.info.setVisitInfo(visitInfo) 

386 exposure.info.setDetector(lsst.afw.cameraGeom.testUtils.DetectorWrapper(id=detector).detector) 

387 return exposure 

388 

389 @staticmethod 

390 def drawGaussian(bbox, instFlux, ellipse): 

391 """Create an image of an elliptical Gaussian. 

392 

393 Parameters 

394 ---------- 

395 bbox : `lsst.geom.Box2I` or `lsst.geom.Box2D` 

396 Bounding box of image to create. 

397 instFlux : `float` 

398 Total instrumental flux of the Gaussian (normalized analytically, 

399 not using pixel values). 

400 ellipse : `lsst.afw.geom.Ellipse` 

401 Defines the centroid and shape. 

402 

403 Returns 

404 ------- 

405 image : `lsst.afw.image.ImageF` 

406 An image of the Gaussian. 

407 """ 

408 x, y = np.meshgrid(np.arange(bbox.getBeginX(), bbox.getEndX()), 

409 np.arange(bbox.getBeginY(), bbox.getEndY())) 

410 t = ellipse.getGridTransform() 

411 xt = t[t.XX] * x + t[t.XY] * y + t[t.X] 

412 yt = t[t.YX] * x + t[t.YY] * y + t[t.Y] 

413 image = lsst.afw.image.ImageF(bbox) 

414 image.array[:, :] = np.exp(-0.5*(xt**2 + yt**2))*instFlux/(2.0*ellipse.getCore().getArea()) 

415 return image 

416 

417 def _installFootprint(self, record, image, setPeakSignificance=True, negative=False): 

418 """Create simulated Footprint and add it to a truth catalog record. 

419 """ 

420 schema = lsst.afw.detection.PeakTable.makeMinimalSchema() 

421 if setPeakSignificance: 

422 schema.addField("significance", type=float, 

423 doc="Ratio of peak value to configured standard deviation.") 

424 if negative: 

425 self.threshold.setPolarity(False) 

426 # Run detection on the single-source image 

427 fpSet = lsst.afw.detection.FootprintSet(image, self.threshold, peakSchema=schema) 

428 # the call below to the FootprintSet ctor is actually a grow operation 

429 fpSet = lsst.afw.detection.FootprintSet(fpSet, int(self.psfShape.getDeterminantRadius() + 1.0), True) 

430 if setPeakSignificance: 

431 # This isn't a traditional significance, since we're using the VALUE 

432 # threshold type, but it's the best we can do in that case. 

433 for footprint in fpSet.getFootprints(): 

434 footprint.updatePeakSignificance(self.threshold.getValue()) 

435 if negative: 

436 self.threshold.setPolarity(True) 

437 # Update the full exposure's mask plane to indicate the detection 

438 fpSet.setMask(self.exposure.mask, "DETECTED") 

439 # Attach the new footprint to the exposure 

440 if len(fpSet.getFootprints()) > 1: 

441 raise RuntimeError("Threshold value results in multiple Footprints for a single object") 

442 if len(fpSet.getFootprints()) == 0: 

443 raise RuntimeError("Threshold value results in zero Footprints for object; " 

444 "did you forget to set negative=True for a negative source?") 

445 record.setFootprint(fpSet.getFootprints()[0]) 

446 

447 def addSource(self, instFlux, centroid, shape=None, setPeakSignificance=True, 

448 negative=False): 

449 """Add a source to the simulation. 

450 

451 To insert a point source with a given signal-to-noise (sn), the total 

452 ``instFlux`` should be: ``sn*noise*psf_scale``, where ``noise`` is the 

453 noise you will pass to ``realize()``, and 

454 ``psf_scale=sqrt(4*pi*r^2)``, where ``r`` is the width of the PSF. 

455 

456 Parameters 

457 ---------- 

458 instFlux : `float` 

459 Total instFlux of the source to be added. 

460 centroid : `lsst.geom.Point2D` 

461 Position of the source to be added. 

462 shape : `lsst.afw.geom.Quadrupole` 

463 Second moments of the source before PSF convolution. Note that the 

464 truth catalog records post-convolution moments. If `None`, a point 

465 source will be added. 

466 setPeakSignificance : `bool` 

467 Set the ``significance`` field for peaks in the footprints? 

468 See ``lsst.meas.algorithms.SourceDetectionTask.setPeakSignificance`` 

469 for how this field is computed for real datasets. 

470 negative : `bool` 

471 Treat this as a negative source, using a negative polarity 

472 Threshold to create the footprint. `instFlux` should be negative 

473 if this is True. 

474 

475 Returns 

476 ------- 

477 record : `lsst.afw.table.SourceRecord` 

478 A truth catalog record. 

479 image : `lsst.afw.image.ImageF` 

480 Single-source image corresponding to the new source. 

481 """ 

482 # Create and set the truth catalog fields 

483 record = self.catalog.addNew() 

484 record.set(self.keys["instFlux"], instFlux) 

485 record.set(self.keys["instFluxErr"], 0) 

486 record.set(self.keys["centroid"], centroid) 

487 record.set(self.keys["negative"], negative) 

488 covariance = np.random.normal(0, 0.1, 4).reshape(2, 2) 

489 covariance[0, 1] = covariance[1, 0] # CovarianceMatrixKey assumes symmetric x_y_Cov 

490 record.set(self.keys["centroid_sigma"], covariance.astype(np.float32)) 

491 if shape is None: 

492 record.set(self.keys["isStar"], True) 

493 fullShape = self.psfShape 

494 else: 

495 record.set(self.keys["isStar"], False) 

496 fullShape = shape.convolve(self.psfShape) 

497 record.set(self.keys["shape"], fullShape) 

498 # Create an image containing just this source 

499 image = self.drawGaussian(self.exposure.getBBox(), instFlux, 

500 lsst.afw.geom.Ellipse(fullShape, centroid)) 

501 # Generate a footprint for this source 

502 self._installFootprint(record, image, setPeakSignificance, negative=negative) 

503 # Actually add the source to the full exposure 

504 self.exposure.image.array[:, :] += image.array 

505 return record, image 

506 

507 def addBlend(self): 

508 """Return a context manager which can add a blend of multiple sources. 

509 

510 Notes 

511 ----- 

512 Note that nothing stops you from creating overlapping sources just using the addSource() method, 

513 but addBlend() is necesssary to create a parent object and deblended HeavyFootprints of the type 

514 produced by the detection and deblending pipelines. 

515 

516 Examples 

517 -------- 

518 .. code-block: py 

519 d = TestDataset(...) 

520 with d.addBlend() as b: 

521 b.addChild(flux1, centroid1) 

522 b.addChild(flux2, centroid2, shape2) 

523 """ 

524 return BlendContext(self) 

525 

526 def transform(self, wcs, **kwds): 

527 """Copy this dataset transformed to a new WCS, with new Psf and PhotoCalib. 

528 

529 Parameters 

530 ---------- 

531 wcs : `lsst.afw.geom.SkyWcs` 

532 WCS for the new dataset. 

533 **kwds 

534 Additional keyword arguments passed on to 

535 `TestDataset.makeEmptyExposure`. If not specified, these revert 

536 to the defaults for `~TestDataset.makeEmptyExposure`, not the 

537 values in the current dataset. 

538 

539 Returns 

540 ------- 

541 newDataset : `TestDataset` 

542 Transformed copy of this dataset. 

543 """ 

544 bboxD = lsst.geom.Box2D() 

545 xyt = lsst.afw.geom.makeWcsPairTransform(self.exposure.getWcs(), wcs) 

546 for corner in lsst.geom.Box2D(self.exposure.getBBox()).getCorners(): 

547 bboxD.include(xyt.applyForward(lsst.geom.Point2D(corner))) 

548 bboxI = lsst.geom.Box2I(bboxD) 

549 result = TestDataset(bbox=bboxI, wcs=wcs, **kwds) 

550 oldPhotoCalib = self.exposure.getPhotoCalib() 

551 newPhotoCalib = result.exposure.getPhotoCalib() 

552 oldPsfShape = self.exposure.getPsf().computeShape(bboxD.getCenter()) 

553 for record in self.catalog: 

554 if record.get(self.keys["nChild"]): 

555 raise NotImplementedError("Transforming blended sources in TestDatasets is not supported") 

556 magnitude = oldPhotoCalib.instFluxToMagnitude(record.get(self.keys["instFlux"])) 

557 newFlux = newPhotoCalib.magnitudeToInstFlux(magnitude) 

558 oldCentroid = record.get(self.keys["centroid"]) 

559 newCentroid = xyt.applyForward(oldCentroid) 

560 if record.get(self.keys["isStar"]): 

561 newDeconvolvedShape = None 

562 else: 

563 affine = lsst.afw.geom.linearizeTransform(xyt, oldCentroid) 

564 oldFullShape = record.get(self.keys["shape"]) 

565 oldDeconvolvedShape = lsst.afw.geom.Quadrupole( 

566 oldFullShape.getIxx() - oldPsfShape.getIxx(), 

567 oldFullShape.getIyy() - oldPsfShape.getIyy(), 

568 oldFullShape.getIxy() - oldPsfShape.getIxy(), 

569 False 

570 ) 

571 newDeconvolvedShape = oldDeconvolvedShape.transform(affine.getLinear()) 

572 result.addSource(newFlux, newCentroid, newDeconvolvedShape) 

573 return result 

574 

575 def realize(self, noise, schema, randomSeed=1): 

576 r"""Simulate an exposure and detection catalog for this dataset. 

577 

578 The simulation includes noise, and the detection catalog includes 

579 `~lsst.afw.detection.heavyFootprint.HeavyFootprint`\ s. 

580 

581 Parameters 

582 ---------- 

583 noise : `float` 

584 Standard deviation of noise to be added to the exposure. The 

585 noise will be Gaussian and constant, appropriate for the 

586 sky-limited regime. 

587 schema : `lsst.afw.table.Schema` 

588 Schema of the new catalog to be created. Must start with 

589 ``self.schema`` (i.e. ``schema.contains(self.schema)`` must be 

590 `True`), but typically contains fields for already-configured 

591 measurement algorithms as well. 

592 randomSeed : `int`, optional 

593 Seed for the random number generator. 

594 If `None`, a seed is chosen automatically. 

595 

596 Returns 

597 ------- 

598 `exposure` : `lsst.afw.image.ExposureF` 

599 Simulated image. 

600 `catalog` : `lsst.afw.table.SourceCatalog` 

601 Simulated detection catalog. 

602 """ 

603 random_state = np.random.RandomState(randomSeed) 

604 assert schema.contains(self.schema) 

605 mapper = lsst.afw.table.SchemaMapper(self.schema) 

606 mapper.addMinimalSchema(self.schema, True) 

607 exposure = self.exposure.clone() 

608 exposure.variance.array[:, :] = noise**2 

609 exposure.image.array[:, :] += random_state.randn(exposure.height, exposure.width)*noise 

610 catalog = lsst.afw.table.SourceCatalog(schema) 

611 catalog.extend(self.catalog, mapper=mapper) 

612 # Loop over sources and generate new HeavyFootprints that divide up 

613 # the noisy pixels, not the ideal no-noise pixels. 

614 for record in catalog: 

615 # parent objects have non-Heavy Footprints, which don't need to be 

616 # updated after adding noise. 

617 if record.getParent() == 0: 

618 continue 

619 # get flattened arrays that correspond to the no-noise and noisy 

620 # parent images 

621 parent = catalog.find(record.getParent()) 

622 footprint = parent.getFootprint() 

623 parentFluxArrayNoNoise = np.zeros(footprint.getArea(), dtype=np.float32) 

624 footprint.spans.flatten(parentFluxArrayNoNoise, self.exposure.image.array, self.exposure.getXY0()) 

625 parentFluxArrayNoisy = np.zeros(footprint.getArea(), dtype=np.float32) 

626 footprint.spans.flatten(parentFluxArrayNoisy, exposure.image.array, exposure.getXY0()) 

627 oldHeavy = record.getFootprint() 

628 fraction = (oldHeavy.getImageArray() / parentFluxArrayNoNoise) 

629 # N.B. this isn't a copy ctor - it's a copy from a vanilla 

630 # Footprint, so it doesn't copy the arrays we don't want to 

631 # change, and hence we have to do that ourselves below. 

632 newHeavy = lsst.afw.detection.HeavyFootprintF(oldHeavy) 

633 newHeavy.getImageArray()[:] = parentFluxArrayNoisy*fraction 

634 newHeavy.getMaskArray()[:] = oldHeavy.getMaskArray() 

635 newHeavy.getVarianceArray()[:] = oldHeavy.getVarianceArray() 

636 record.setFootprint(newHeavy) 

637 lsst.afw.table.updateSourceCoords(exposure.wcs, catalog) 

638 return exposure, catalog 

639 

640 

641class AlgorithmTestCase: 

642 """Base class for tests of measurement tasks. 

643 """ 

644 def makeSingleFrameMeasurementConfig(self, plugin=None, dependencies=()): 

645 """Create an instance of `SingleFrameMeasurementTask.ConfigClass`. 

646 

647 Only the specified plugin and its dependencies will be run; the 

648 Centroid, Shape, and ModelFlux slots will be set to the truth fields 

649 generated by the `TestDataset` class. 

650 

651 Parameters 

652 ---------- 

653 plugin : `str` 

654 Name of measurement plugin to enable. 

655 dependencies : iterable of `str`, optional 

656 Names of dependencies of the measurement plugin. 

657 

658 Returns 

659 ------- 

660 config : `SingleFrameMeasurementTask.ConfigClass` 

661 The resulting task configuration. 

662 """ 

663 config = SingleFrameMeasurementTask.ConfigClass() 

664 with warnings.catch_warnings(): 

665 warnings.filterwarnings("ignore", message="ignoreSlotPluginChecks", category=FutureWarning) 

666 config = SingleFrameMeasurementTask.ConfigClass(ignoreSlotPluginChecks=True) 

667 config.slots.centroid = "truth" 

668 config.slots.shape = "truth" 

669 config.slots.modelFlux = None 

670 config.slots.apFlux = None 

671 config.slots.psfFlux = None 

672 config.slots.gaussianFlux = None 

673 config.slots.calibFlux = None 

674 config.plugins.names = (plugin,) + tuple(dependencies) 

675 return config 

676 

677 def makeSingleFrameMeasurementTask(self, plugin=None, dependencies=(), config=None, schema=None, 

678 algMetadata=None): 

679 """Create a configured instance of `SingleFrameMeasurementTask`. 

680 

681 Parameters 

682 ---------- 

683 plugin : `str`, optional 

684 Name of measurement plugin to enable. If `None`, a configuration 

685 must be supplied as the ``config`` parameter. If both are 

686 specified, ``config`` takes precedence. 

687 dependencies : iterable of `str`, optional 

688 Names of dependencies of the specified measurement plugin. 

689 config : `SingleFrameMeasurementTask.ConfigClass`, optional 

690 Configuration for the task. If `None`, a measurement plugin must 

691 be supplied as the ``plugin`` paramter. If both are specified, 

692 ``config`` takes precedence. 

693 schema : `lsst.afw.table.Schema`, optional 

694 Measurement table schema. If `None`, a default schema is 

695 generated. 

696 algMetadata : `lsst.daf.base.PropertyList`, optional 

697 Measurement algorithm metadata. If `None`, a default container 

698 will be generated. 

699 

700 Returns 

701 ------- 

702 task : `SingleFrameMeasurementTask` 

703 A configured instance of the measurement task. 

704 """ 

705 if config is None: 

706 if plugin is None: 

707 raise ValueError("Either plugin or config argument must not be None") 

708 config = self.makeSingleFrameMeasurementConfig(plugin=plugin, dependencies=dependencies) 

709 if schema is None: 

710 schema = TestDataset.makeMinimalSchema() 

711 # Clear all aliases so only those defined by config are set. 

712 schema.setAliasMap(None) 

713 if algMetadata is None: 

714 algMetadata = lsst.daf.base.PropertyList() 

715 return SingleFrameMeasurementTask(schema=schema, algMetadata=algMetadata, config=config) 

716 

717 def makeForcedMeasurementConfig(self, plugin=None, dependencies=()): 

718 """Create an instance of `ForcedMeasurementTask.ConfigClass`. 

719 

720 In addition to the plugins specified in the plugin and dependencies 

721 arguments, the `TransformedCentroid` and `TransformedShape` plugins 

722 will be run and used as the centroid and shape slots; these simply 

723 transform the reference catalog centroid and shape to the measurement 

724 coordinate system. 

725 

726 Parameters 

727 ---------- 

728 plugin : `str` 

729 Name of measurement plugin to enable. 

730 dependencies : iterable of `str`, optional 

731 Names of dependencies of the measurement plugin. 

732 

733 Returns 

734 ------- 

735 config : `ForcedMeasurementTask.ConfigClass` 

736 The resulting task configuration. 

737 """ 

738 

739 config = ForcedMeasurementTask.ConfigClass() 

740 config.slots.centroid = "base_TransformedCentroid" 

741 config.slots.shape = "base_TransformedShape" 

742 config.slots.modelFlux = None 

743 config.slots.apFlux = None 

744 config.slots.psfFlux = None 

745 config.slots.gaussianFlux = None 

746 config.plugins.names = (plugin,) + tuple(dependencies) + ("base_TransformedCentroid", 

747 "base_TransformedShape") 

748 return config 

749 

750 def makeForcedMeasurementTask(self, plugin=None, dependencies=(), config=None, refSchema=None, 

751 algMetadata=None): 

752 """Create a configured instance of `ForcedMeasurementTask`. 

753 

754 Parameters 

755 ---------- 

756 plugin : `str`, optional 

757 Name of measurement plugin to enable. If `None`, a configuration 

758 must be supplied as the ``config`` parameter. If both are 

759 specified, ``config`` takes precedence. 

760 dependencies : iterable of `str`, optional 

761 Names of dependencies of the specified measurement plugin. 

762 config : `SingleFrameMeasurementTask.ConfigClass`, optional 

763 Configuration for the task. If `None`, a measurement plugin must 

764 be supplied as the ``plugin`` paramter. If both are specified, 

765 ``config`` takes precedence. 

766 refSchema : `lsst.afw.table.Schema`, optional 

767 Reference table schema. If `None`, a default schema is 

768 generated. 

769 algMetadata : `lsst.daf.base.PropertyList`, optional 

770 Measurement algorithm metadata. If `None`, a default container 

771 will be generated. 

772 

773 Returns 

774 ------- 

775 task : `ForcedMeasurementTask` 

776 A configured instance of the measurement task. 

777 """ 

778 if config is None: 

779 if plugin is None: 

780 raise ValueError("Either plugin or config argument must not be None") 

781 config = self.makeForcedMeasurementConfig(plugin=plugin, dependencies=dependencies) 

782 if refSchema is None: 

783 refSchema = TestDataset.makeMinimalSchema() 

784 if algMetadata is None: 

785 algMetadata = lsst.daf.base.PropertyList() 

786 return ForcedMeasurementTask(refSchema=refSchema, algMetadata=algMetadata, config=config) 

787 

788 

789class TransformTestCase: 

790 """Base class for testing measurement transformations. 

791 

792 Notes 

793 ----- 

794 We test both that the transform itself operates successfully (fluxes are 

795 converted to magnitudes, flags are propagated properly) and that the 

796 transform is registered as the default for the appropriate measurement 

797 algorithms. 

798 

799 In the simple case of one-measurement-per-transformation, the developer 

800 need not directly write any tests themselves: simply customizing the class 

801 variables is all that is required. More complex measurements (e.g. 

802 multiple aperture fluxes) require extra effort. 

803 """ 

804 name = "MeasurementTransformTest" 

805 """The name used for the measurement algorithm (str). 

806