Coverage for python/lsst/ip/diffim/computeSpatiallySampledMetrics.py: 17%

1# This file is part of ip_diffim.

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/>.

22import numpy as np

23import scipy.signal

25import lsst.geom

27import lsst.afw.image as afwImage

28import lsst.afw.table as afwTable

29import lsst.pipe.base as pipeBase

30import lsst.pex.config as pexConfig

32from lsst.ip.diffim.utils import getPsfFwhm, angleMean, evaluateMaskFraction, getKernelCenterDisplacement

33from lsst.meas.algorithms import SkyObjectsTask

34from lsst.pex.exceptions import InvalidParameterError, RangeError

35from lsst.utils.timer import timeMethod

37import lsst.utils

39__all__ = ["SpatiallySampledMetricsConfig", "SpatiallySampledMetricsTask"]

42class SpatiallySampledMetricsConnections(pipeBase.PipelineTaskConnections,

43 dimensions=("instrument", "visit", "detector"),

44 defaultTemplates={"coaddName": "deep",

45 "warpTypeSuffix": "",

46 "fakesType": ""}):

47 science = pipeBase.connectionTypes.Input(

48 doc="Input science exposure.",

49 dimensions=("instrument", "visit", "detector"),

50 storageClass="ExposureF",

51 name="{fakesType}calexp"

52 )

53 template = pipeBase.connectionTypes.Input(

54 doc="Warped and not PSF-matched template used to create the difference image.",

55 dimensions=("instrument", "visit", "detector"),

56 storageClass="ExposureF",

57 name="{fakesType}{coaddName}Diff_templateExp",

58 )

59 difference = pipeBase.connectionTypes.Input(

60 doc="Difference image with detection mask plane filled in.",

61 dimensions=("instrument", "visit", "detector"),

62 storageClass="ExposureF",

63 name="{fakesType}{coaddName}Diff_differenceExp",

64 )

65 diaSources = pipeBase.connectionTypes.Input(

66 doc="Filtered diaSources on the difference image.",

67 dimensions=("instrument", "visit", "detector"),

68 storageClass="ArrowAstropy",

69 name="{fakesType}dia_source_detector",

70 )

71 psfMatchingKernel = pipeBase.connectionTypes.Input(

72 doc="Kernel used to PSF match the science and template images.",

73 dimensions=("instrument", "visit", "detector"),

74 storageClass="MatchingKernel",

75 name="{fakesType}{coaddName}Diff_psfMatchKernel",

76 )

77 spatiallySampledMetrics = pipeBase.connectionTypes.Output(

78 doc="Summary metrics computed at randomized locations.",

79 dimensions=("instrument", "visit", "detector"),

80 storageClass="ArrowAstropy",

81 name="{fakesType}{coaddName}Diff_spatiallySampledMetrics",

82 )

85class SpatiallySampledMetricsConfig(pipeBase.PipelineTaskConfig,

86 pipelineConnections=SpatiallySampledMetricsConnections):

87 """Config for SpatiallySampledMetricsTask

88 """

89 metricsMaskPlanes = lsst.pex.config.ListField(

90 dtype=str,

91 doc="List of mask planes to include in metrics",

92 default=('BAD', 'CLIPPED', 'CR', 'DETECTED', 'DETECTED_NEGATIVE', 'EDGE',

93 'INEXACT_PSF', 'INJECTED', 'INJECTED_TEMPLATE', 'INTRP', 'NOT_DEBLENDED',

94 'NO_DATA', 'REJECTED', 'SAT', 'SAT_TEMPLATE', 'SENSOR_EDGE', 'STREAK', 'SUSPECT',

95 'UNMASKEDNAN',

96 ),

97 )

98 metricSources = pexConfig.ConfigurableField(

99 target=SkyObjectsTask,

100 doc="Generate QA metric sources",

101 )

102

103 def setDefaults(self):

104 self.metricSources.avoidMask = ["NO_DATA", "EDGE"]

105

106

107class SpatiallySampledMetricsTask(lsst.pipe.base.PipelineTask):

108 """Detect and measure sources on a difference image.

109 """

110 ConfigClass = SpatiallySampledMetricsConfig

111 _DefaultName = "spatiallySampledMetrics"

112

113 def __init__(self, **kwargs):

114 super().__init__(**kwargs)

115

116 self.makeSubtask("metricSources")

117 self.schema = afwTable.SourceTable.makeMinimalSchema()

118 self.schema.addField(

119 "x", "F",

120 "X location of the metric evaluation.",

121 units="pixel")

122 self.schema.addField(

123 "y", "F",

124 "Y location of the metric evaluation.",

125 units="pixel")

126 self.metricSources.skySourceKey = self.schema.addField("sky_source", type="Flag",

127 doc="Metric evaluation objects.")

128 self.schema.addField(

129 "source_density", "F",

130 "Density of diaSources at location.",

131 units="count/degree^2")

132 self.schema.addField(

133 "dipole_density", "F",

134 "Density of dipoles at location.",

135 units="count/degree^2")

136 self.schema.addField(

137 "dipole_direction", "F",

138 "Mean dipole orientation.",

139 units="radian")

140 self.schema.addField(

141 "dipole_separation", "F",

142 "Mean dipole separation.",

143 units="pixel")

144 self.schema.addField(

145 "template_value", "F",

146 "Median of template at location.",

147 units="nJy")

148 self.schema.addField(

149 "template_variance", "F",

150 "Median of template variance at location.",

151 units="nJy^2")

152 self.schema.addField(

153 "science_value", "F",

154 "Median of science at location.",

155 units="nJy")

156 self.schema.addField(

157 "science_variance", "F",

158 "Median of science variance at location.",

159 units="nJy^2")

160 self.schema.addField(

161 "diffim_value", "F",

162 "Median of diffim at location.",

163 units="nJy")

164 self.schema.addField(

165 "diffim_variance", "F",

166 "Median of diffim variance at location.",

167 units="nJy^2")

168 self.schema.addField(

169 "diffim_chi2PerPix", "F",

170 "Robust normalized noise of diffim at location:"

171 " (1.4826*MAD(image))^2 / median(variance), evaluated on background"

172 " pixels (DETECTED, DETECTED_NEGATIVE, BAD, SAT, EDGE, NO_DATA excluded)."

173 " Expected ~1.0 for a well-decorrelated diffim; values >>1 indicate"

174 " residual structure, <<1 indicates over-estimated variance.")

175 self.schema.addField(

176 "science_psfSize", "F",

177 "Width of the science image PSF at location.",

178 units="pixel")

179 self.schema.addField(

180 "template_psfSize", "F",

181 "Width of the template image PSF at location.",

182 units="pixel")

183 for maskPlane in self.config.metricsMaskPlanes:

184 self.schema.addField(

185 "%s_mask_fraction"%maskPlane.lower(), "F",

186 "Fraction of pixels with %s mask"%maskPlane

187 )

188 self.schema.addField(

189 "psfMatchingKernel_sum", "F",

190 "PSF matching kernel sum at location.")

191 self.schema.addField(

192 "psfMatchingKernel_dx", "F",

193 "PSF matching kernel centroid offset in x at location.",

194 units="pixel")

195 self.schema.addField(

196 "psfMatchingKernel_dy", "F",

197 "PSF matching kernel centroid offset in y at location.",

198 units="pixel")

199 self.schema.addField(

200 "psfMatchingKernel_length", "F",

201 "PSF matching kernel centroid offset module.",

202 units="arcsecond")

203 self.schema.addField(

204 "psfMatchingKernel_position_angle", "F",

205 "PSF matching kernel centroid offset position angle.",

206 units="radian")

207 self.schema.addField(

208 "psfMatchingKernel_direction", "F",

209 "PSF matching kernel centroid offset direction in detector plane.",

210 units="radian")

211 self.schema.addField(

212 "psfMatchingKernel_residualNorm", "F",

213 "Shape-only PSF match residual at location:"

214 "L2 norm of (K-convolved template PSF - science PSF),"

215 " relative to the science PSF L2 norm. Larger values indicate worse"

216 " PSF matching. Assumes the kernel was solved to convolve the template.")

217

218 @timeMethod

219 def run(self, science, template, difference, diaSources, psfMatchingKernel):

220 """Calculate difference image metrics on specific locations across the images

221

222 Parameters

223 ----------

224 science : `lsst.afw.image.ExposureF`

225 Science exposure that the template was subtracted from.

226 template : `lsst.afw.image.ExposureF`

227 Warped and non PSF-matched template that was used produce

228 the difference image.

229 difference : `lsst.afw.image.ExposureF`

230 Result of subtracting template from the science image.

231 diaSources : `lsst.afw.table.SourceCatalog`

232 The catalog of detected sources.

233 psfMatchingKernel : `~lsst.afw.math.LinearCombinationKernel`

234 The PSF matching kernel of the subtraction to evaluate.

235

236 Returns

237 -------

238 results : `lsst.pipe.base.Struct`

239 ``spatiallySampledMetrics`` : `astropy.table.Table`

240 Image quality metrics spatially sampled locations.

241 """

242

243 idFactory = lsst.meas.base.IdGenerator().make_table_id_factory()

244

245 spatiallySampledMetrics = afwTable.SourceCatalog(self.schema)

246 spatiallySampledMetrics.getTable().setIdFactory(idFactory)

247

248 self.metricSources.run(mask=science.mask, seed=difference.info.id, catalog=spatiallySampledMetrics)

249

250 metricsMaskPlanes = []

251 for maskPlane in self.config.metricsMaskPlanes:

252 try:

253 metricsMaskPlanes.append(maskPlane)

254 except InvalidParameterError:

255 self.log.info("Unable to calculate metrics for mask plane %s: not in image"%maskPlane)

256

257 for src in spatiallySampledMetrics:

258 self._evaluateLocalMetric(src, science, template, difference, diaSources,

259 metricsMaskPlanes=metricsMaskPlanes,

260 psfMatchingKernel=psfMatchingKernel)

261 spatiallySampledMetrics = spatiallySampledMetrics.copy(deep=True).asAstropy()

262 return pipeBase.Struct(spatiallySampledMetrics=spatiallySampledMetrics)

263

264 def _evaluateLocalMetric(self, src, science, template, difference, diaSources,

265 metricsMaskPlanes, psfMatchingKernel):

266 """Calculate image quality metrics at spatially sampled locations.

267

268 Parameters

269 ----------

270 src : `lsst.afw.table.SourceRecord`

271 The source record to be updated with metric calculations.

272 diaSources : `lsst.afw.table.SourceCatalog`

273 The catalog of detected sources.

274 science : `lsst.afw.image.Exposure`

275 The science image.

276 difference : `lsst.afw.image.Exposure`

277 Result of subtracting template from the science image.

278 metricsMaskPlanes : `list` of `str`

279 Mask planes to calculate metrics from.

280 psfMatchingKernel : `~lsst.afw.math.LinearCombinationKernel`

281 The PSF matching kernel of the subtraction to evaluate.

282 """

283 bbox = src.getFootprint().getBBox()

284 pix = bbox.getCenter()

285 src.set('science_psfSize', getPsfFwhm(science.psf, position=pix))

286 try:

287 src.set('template_psfSize', getPsfFwhm(template.psf, position=pix))

288 except (InvalidParameterError, RangeError):

289 src.set('template_psfSize', np.nan)

290

291 metricRegionSize = 100

292 bbox.grow(metricRegionSize)

293 bbox = bbox.clippedTo(science.getBBox())

294 nPix = bbox.getArea()

295 pixScale = science.wcs.getPixelScale(bbox.getCenter())

296 area = nPix*pixScale.asDegrees()**2

297 peak = src.getFootprint().getPeaks()[0]

298 src.set('x', peak['i_x'])

299 src.set('y', peak['i_y'])

300 src.setCoord(science.wcs.pixelToSky(peak['i_x'], peak['i_y']))

301 selectSources = diaSources[bbox.contains(diaSources['x'], diaSources['y'])]

302 sourceDensity = len(selectSources)/area

303 dipoleSources = selectSources[selectSources["isDipole"]]

304 dipoleDensity = len(dipoleSources)/area

305

306 if dipoleSources:

307 meanDipoleOrientation = angleMean(dipoleSources["dipoleAngle"])

308 src.set('dipole_direction', meanDipoleOrientation)

309 meanDipoleSeparation = np.mean(dipoleSources["dipoleLength"])

310 src.set('dipole_separation', meanDipoleSeparation)

311

312 templateVal = np.median(template[bbox].image.array)

313 templateVar = np.median(template[bbox].variance.array)

314 scienceVal = np.median(science[bbox].image.array)

315 scienceVar = np.median(science[bbox].variance.array)

316 diffimVal = np.median(difference[bbox].image.array)

317 diffimVar = np.median(difference[bbox].variance.array)

318 src.set('source_density', sourceDensity)

319 src.set('dipole_density', dipoleDensity)

320 src.set('template_value', templateVal)

321 src.set('template_variance', templateVar)

322 src.set('science_value', scienceVal)

323 src.set('science_variance', scienceVar)

324 src.set('diffim_value', diffimVal)

325 src.set('diffim_variance', diffimVar)

326 src.set('diffim_chi2PerPix', self._diffimChi2PerPix(difference[bbox]))

327 for maskPlane in metricsMaskPlanes:

328 src.set("%s_mask_fraction"%maskPlane.lower(),

329 evaluateMaskFraction(difference.mask[bbox], maskPlane)

330 )

331

332 krnlSum, dx, dy, direction, length = getKernelCenterDisplacement(

333 psfMatchingKernel, src.get('x'), src.get('y'))

334

335 point1 = lsst.geom.SpherePoint(

336 src.get('coord_ra'), src.get('coord_dec'),

337 lsst.geom.radians)

338 point2 = science.wcs.pixelToSky(src.get('x') + dx, src.get('y') + dy)

339 bearing = point1.bearingTo(point2)

340 pa_ref_angle = lsst.geom.Angle(np.pi/2, lsst.geom.radians)

341 pa = pa_ref_angle - bearing

342 # Wrap around to get Delta_RA from -pi to +pi

343 pa = pa.wrapCtr()

344 position_angle = pa.asRadians()

345

346 src.set('psfMatchingKernel_sum', krnlSum)

347 src.set('psfMatchingKernel_dx', dx)

348 src.set('psfMatchingKernel_dy', dy)

349 src.set('psfMatchingKernel_length', length*pixScale.asArcseconds())

350 src.set('psfMatchingKernel_position_angle', position_angle) # in E of N position angle

351 src.set('psfMatchingKernel_direction', direction) # direction offset in detector

352

353 src.set('psfMatchingKernel_residualNorm',

354 self._psfMatchResidualNorm(psfMatchingKernel, science.psf, template.psf,

355 src.get('x'), src.get('y')))

356

357 def _diffimChi2PerPix(self, difference):

358 """Robust normalized noise of the difference image.

359

360 Computes ``(1.4826 * MAD(image))^2 / median(variance)`` on

361 background-only pixels. The 1.4826 factor rescales the

362 Median Absolute Deviation to a Gaussian-sigma estimate, so the

363 ratio has expectation 1.0 on pure noise.

364

365 Returns NaN if no usable pixels remain or the median variance

366 is non-positive.

367 """

368 image = difference.image.array

369 variance = difference.variance.array

370 mask = difference.mask

371 excludePlanes = [p for p in ("DETECTED", "DETECTED_NEGATIVE", "BAD",

372 "SAT", "EDGE", "NO_DATA")

373 if p in mask.getMaskPlaneDict()]

374 if excludePlanes:

375 badBits = mask.getPlaneBitMask(excludePlanes)

376 good = (mask.array & badBits) == 0

377 else:

378 good = np.ones(image.shape, dtype=bool)

379 good &= np.isfinite(image) & np.isfinite(variance) & (variance > 0)

380 if not np.any(good):

381 return np.nan

382 varMed = np.median(variance[good])

383 if not np.isfinite(varMed) or varMed <= 0:

384 return np.nan

385 imgGood = image[good]

386 mad = np.median(np.abs(imgGood - np.median(imgGood)))

387 return float((1.4826*mad)**2/varMed)

388

389 def _psfMatchResidualNorm(self, kernel, sciencePsf, templatePsf, x, y):

390 """Relative L2 norm of the PSF-match residual at (x, y).

391

392 Convolves the template PSF with the matching kernel evaluated at

393 (x, y) and compares to the science PSF at the same position, both

394 renormalized to unit sum so the result captures shape mismatch only.

395 The kernel is assumed to convolve the template.

396

397 Returns NaN if either PSF cannot be evaluated at the position or

398 the resulting images cannot be normalized.

399 """

400 point = lsst.geom.Point2D(x, y)

401 try:

402 psfSci = sciencePsf.computeKernelImage(point).array

403 psfTmp = templatePsf.computeKernelImage(point).array

404 except (InvalidParameterError, RangeError):

405 return np.nan

406

407 kImage = afwImage.ImageD(kernel.getDimensions())

408 kernel.computeImage(kImage, doNormalize=True, x=x, y=y)

409 matched = scipy.signal.fftconvolve(psfTmp, kImage.array, mode='same')

410

411 matchedSum = matched.sum()

412 psfSciSum = psfSci.sum()

413 if matchedSum <= 0 or psfSciSum <= 0:

414 return np.nan

415 matched = matched/matchedSum

416 psfSci = psfSci/psfSciSum

417

418 # PSF stamps may differ in size between the two exposures; crop

419 # both to a common centered region before differencing.

420 h = min(matched.shape[0], psfSci.shape[0])

421 w = min(matched.shape[1], psfSci.shape[1])

422

423 def _crop(a):

424 sy = (a.shape[0] - h)//2

425 sx = (a.shape[1] - w)//2

426 return a[sy:sy + h, sx:sx + w]

427 matched = _crop(matched)

428 psfSci = _crop(psfSci)

429

430 sciNorm = np.sqrt(np.sum(psfSci**2))

431 if sciNorm == 0:

432 return np.nan

433 return float(np.sqrt(np.sum((matched - psfSci)**2))/sciNorm)