Coverage for metadetect/lsst/metadetect.py: 18%

1import logging

2import numpy as np

3import ngmix

4from ngmix.gexceptions import BootPSFFailure

5from lsst.pex.config import (

6 Config,

7 ConfigField,

8 ConfigurableField,

9 Field,

10 FieldValidationError,

11 ListField,

12)

13from lsst.pipe.base import Task

14from lsst.meas.algorithms import SourceDetectionTask

16from .. import procflags

18from .skysub import subtract_sky_mbexp

20from .defaults import (

21 DEFAULT_STAMP_SIZE,

22 DEFAULT_SUBTRACT_SKY,

23 DEFAULT_PGAUSS_FWHM,

24)

25from . import measure

26from .metacal_exposures import get_metacal_mbexps_fixnoise

27from .util import get_integer_center, get_jacobian, override_config

29LOG = logging.getLogger('lsst_metadetect')

32def run_metadetect(

33 mbexp,

34 noise_mbexp,

35 rng,

36 mfrac_mbexp=None,

37 ormasks=None,

38 config=None,

39 show=False,

40):

41 """

42 Run metadetection on the input MultiBandObsList

44 Note that bright object masking must be applied outside of this function.

46 Parameters

47 ----------

48 mbexp: lsst.afw.image.MultibandExposure

49 The exposures to process

50 noise_mbexp: lsst.afw.image.MultibandExposure

51 The noise exposures for metacal

52 mfrac_mbexp: lsst.afw.image.MultibandExposure, optional

53 The fraction of masked exposures for the pixel; for coadds this is the

54 fraction of input images contributing to each pixel that were masked

55 ormasks: list of images, optional

56 A list of logical or masks, such as created for all images that went

57 into a coadd.

59 Note when coadding an ormask is created in the .mask attribute. But

60 this code expects the mask attribute for each exposure to be not an or

61 of all masks from the original exposures, but a mask indicating problem

62 areas such as bright objects or apodized edges.

64 In the future we expect the MultibandExposure to have an ormask

65 attribute

66 rng: np.random.RandomState

67 Random number generator

68 config: dict, optional

69 Configuration for the fitter, metacal, psf, detect, Entries

70 in this dict override defaults; see lsst_configs.py

71 show: bool, optional

72 if set to True images will be shown

74 Returns

75 -------

76 result dict

77 This is keyed by shear string 'noshear', '1p', ... or None if there was

78 a problem doing the metacal steps; this only happens if the setting

79 metacal_psf is set to 'fitgauss' and the fitting fails

80 """

82 config_override = config if config is not None else {}

83 config = MetadetectConfig()

84 config.setDefaults()

86 override_config(config, config_override)

88 config.freeze()

89 config.validate()

90 task = MetadetectTask(config=config)

91 result = task.run(

92 mbexp,

93 noise_mbexp,

94 rng,

95 mfrac_mbexp,

96 ormasks,

97 show=show,

98 )

99 return result

100

101

102class PGaussConfig(Config):

103 fwhm = Field[float](

104 doc="FWHM of the Gaussian weight function for PGauss (in arcseconds)",

105 default=DEFAULT_PGAUSS_FWHM,

106 )

107

108

109class MetacalConfig(Config):

110 types = ListField[str](

111 doc="List of artificial shears to apply.",

112 default=[

113 "noshear",

114 "1p",

115 "1m",

116 ],

117 )

118

119 def validate(self):

120 super().validate()

121 if not set(self.types).issubset({"noshear", "1p", "1m", "2p", "2m"}):

122 raise FieldValidationError(

123 self.__class__.types,

124 self,

125 "types must be a list consisting of any combinations of "

126 "{'noshear', '1p', '1m', '2p', '2m'}",

127 )

128

129

130class MetadetectConfig(Config):

131 subtract_sky = Field[bool](

132 doc="Whether to subtract the sky before running metadetect",

133 default=DEFAULT_SUBTRACT_SKY,

134 )

135

136 pgauss = ConfigField[PGaussConfig](

137 doc="PGauss config",

138 )

139

140 detect = ConfigurableField(

141 doc="Detection config",

142 target=SourceDetectionTask,

143 )

144

145 metacal = ConfigField[MetacalConfig](

146 doc="Metacal config",

147 )

148

149 stamp_size = Field[int](

150 doc="The stamp size to use for measurements",

151 default=DEFAULT_STAMP_SIZE,

152 )

153

154 shear_bands = ListField[str](

155 doc="List of bands to use for shear measurements. Default is to use all bands.",

156 default=None,

157 optional=True,

158 )

159

160 def setDefaults(self):

161 super().setDefaults()

162

163 def validate(self):

164 super().validate()

165

166

167class MetadetectTask(Task):

168 ConfigClass = MetadetectConfig

169 _DefaultName = "metadetect"

170

171 def __init__(self, *args, **kwargs):

172 super().__init__(*args, **kwargs)

173 self.makeSubtask("detect")

174

175 def run(

176 self,

177 mbexp,

178 noise_mbexp,

179 rng,

180 mfrac_mbexp=None,

181 ormasks=None,

182 show=False,

183 ):

184 # This is to support methods that are not yet refactored.

185 config = self.config.toDict()

186 config['detect']['thresh'] = self.detect.config.thresholdValue

187

188 ormask = combine_ormasks(mbexp, ormasks)

189 mfrac, wgts = get_mfrac_mbexp(mbexp=mbexp, mfrac_mbexp=mfrac_mbexp)

190

191 if self.config.subtract_sky:

192 subtract_sky_mbexp(

193 mbexp=mbexp, thresh=self.config.detect.thresholdValue

194 )

195

196 psf_stats = fit_original_psfs_mbexp(

197 mbexp=mbexp,

198 wgts=wgts,

199 rng=rng,

200 )

201

202 metacal_types = config['metacal'].get('types', None)

203

204 mdict, _ = get_metacal_mbexps_fixnoise(

205 mbexp=mbexp,

206 noise_mbexp=noise_mbexp,

207 types=metacal_types,

208 )

209

210 result = {}

211 for shear_str, mcal_mbexp in mdict.items():

212 res = detect_deblend_and_measure(

213 mbexp=mcal_mbexp,

214 config=config,

215 rng=rng,

216 show=show,

217 )

218

219 if res is not None:

220 band = mcal_mbexp.bands[0]

221 exp = mcal_mbexp[band]

222

223 add_mfrac(config=config, mfrac=mfrac, res=res, exp=exp)

224 add_ormask(ormask, res)

225 add_original_psf(psf_stats, res)

226

227 result[shear_str] = res

228

229 return result

230

231

232def detect_deblend_and_measure(

233 mbexp,

234 config,

235 rng,

236 show=False,

237):

238 """

239 run detection, deblending and measurements.

240

241 Note deblending is always run in a hierarchical detection process, but the

242 deblending is only used for getting centers, and because there is currently

243 no other way to split footprints

244

245 Parameters

246 ----------

247 mbexp: lsst.afw.image.MultibandExposure

248 The metacal'ed exposures to process

249 config: dict, optional

250 Configuration for the fitter, metacal, psf, detect, Entries

251 in this dict override defaults; see lsst_configs.py

252 rng: np.random.RandomState

253 Random number generator

254 show: bool, optional

255 If set to True, show images during processing

256 """

257

258 LOG.info('measuring with blended stamps')

259

260 sources, detexp = measure.detect_and_deblend(

261 mbexp=mbexp,

262 rng=rng,

263 thresh=config['detect']['thresh'],

264 show=show,

265 )

266

267 results = measure.measure(

268 mbexp=mbexp,

269 detexp=detexp,

270 sources=sources,

271 config=config,

272 rng=rng,

273 )

274

275 return results

276

277

278def add_mfrac(config, mfrac, res, exp):

279 """

280 calculate and add mfrac to the input result array

281 """

282 if np.any(mfrac > 0):

283 # we are using the positions with the metacal shear removed for

284 # this.

285

286 cen, _ = get_integer_center(

287 wcs=exp.getWcs(),

288 bbox=exp.getBBox(),

289 as_double=True,

290 )

291 jac = get_jacobian(exp=exp, cen=cen)

292

293 res['mfrac'] = measure_weighted_mfrac(

294 mfrac=mfrac,

295 x=res['col'] - res['col0'],

296 y=res['row'] - res['row0'],

297 jac=jac,

298 fwhm=config.get('mfrac_fwhm', None),

299 )

300 else:

301 res['mfrac'] = 0

302

303

304def measure_weighted_mfrac(

305 *,

306 mfrac,

307 x,

308 y,

309 jac,

310 fwhm,

311):

312 """

313 Measure a Gaussian-weighted average of an image.

314

315 This function is meant to be used with images that represent the fraction

316 of single-epoch images that are masked in each pixel of a coadd. It

317 computes a Gaussian-weighted average of the image at a list of locations.

318

319 Parameters

320 ----------

321 mfrac : np.ndarray

322 The input image with which to compute the weighted averages.

323 x : np.ndarray

324 The input x/col values for the positions at which to compute the

325 weighted average.

326 y : np.ndarray

327 The input y/row values for the positions at which to compute the

328 weighted average.

329 box_sizes : np.ndarray

330 The size of the stamp to use to measure the weighted average. Should be

331 big enough to hold 2 * `fwhm`.

332 jac: ngmix.Jacobian

333 The jacobian for the data

334 fwhm : float or None

335 The FWHM of the Gaussian aperture in arcseconds. If None, a default

336 of 1.2 is used.

337

338 Returns

339 -------

340 mfracs : np.ndarray

341 The weighted averages at each input location.

342 """

343

344 if fwhm is None:

345 fwhm = 1.2

346

347 ny, nx = mfrac.shape

348

349 gauss_wgt = ngmix.GMixModel(

350 [0, 0, 0, 0, ngmix.moments.fwhm_to_T(fwhm), 1],

351 'gauss',

352 )

353 sigma = ngmix.moments.fwhm_to_sigma(fwhm)

354 box_rad = int(round(sigma * 5))

355

356 mfracs = []

357 for i in range(x.shape[0]):

358 ix = int(np.floor(x[i] + 0.5))

359 iy = int(np.floor(y[i] + 0.5))

360

361 xstart = ix - box_rad

362 xend = ix + box_rad + 1

363 ystart = iy - box_rad

364 yend = iy + box_rad + 1

365

366 if xstart < 0:

367 xstart = 0

368 if ystart < 0:

369 ystart = 0

370 if xend > nx:

371 xend = nx

372 if yend > ny:

373 yend = ny

374

375 sub_mfrac = mfrac[xstart:xend, ystart:yend]

376 if sub_mfrac.size == 0:

377 mfracs.append(1.0)

378 else:

379 cy, cx = (y[i] - ystart, x[i] - xstart)

380 this_jac = jac.copy()

381 this_jac.set_cen(row=cy, col=cx)

382

383 obs = ngmix.Observation(

384 image=sub_mfrac,

385 jacobian=this_jac,

386 )

387

388 stats = gauss_wgt.get_weighted_sums(obs, maxrad=box_rad)

389

390 # this is the weighted average in the image using the

391 # Gaussian as the weight.

392 mfracs.append(stats["sums"][5] / stats["wsum"])

393

394 return np.array(mfracs)

395

396

397def add_ormask(ormask, res):

398 """

399 copy in ormask values using the row, col positions

400 """

401 for i in range(res.size):

402 row_diff = res['row'][i] - res['row0'][i]

403 col_diff = res['col'][i] - res['col0'][i]

404 local_row = int(np.floor(row_diff + 0.5))

405 local_col = int(np.floor(col_diff + 0.5))

406

407 res['ormask'][i] = ormask[local_row, local_col]

408

409

410def add_original_psf(psf_stats, res):

411 """

412 copy in psf results

413 """

414 res['psfrec_flags'][:] = psf_stats['flags']

415 res['psfrec_g'][:, 0] = psf_stats['g1']

416 res['psfrec_g'][:, 1] = psf_stats['g2']

417 res['psfrec_T'][:] = psf_stats['T']

418

419

420def combine_ormasks(mbexp, ormasks):

421 """

422 logical or together all the ormasks, or if ormasks is None create zeroed

423 versions for each band

424 """

425 if ormasks is None:

426 bands = mbexp.bands

427 dims = mbexp[bands[0]].image.array.shape

428 ormask = np.zeros(dims, dtype='i4')

429 else:

430 for imask, tormask in enumerate(ormasks):

431 if imask == 0:

432 ormask = tormask.copy()

433 else:

434 ormask |= tormask

435

436 return ormask

437

438

439def get_mfrac_mbexp(mbexp, mfrac_mbexp):

440 """

441 set the masked fraction image, averaged over all bands

442

443 Parameters

444 ----------

445 mbexp: lsst.afw.image.MultibandExposure

446 The exposures to process

447

448 Returns

449 -------

450 mfrac: array

451 """

452 wgts = []

453

454 wsum = 0.0

455

456 mfrac = None

457 for exp, mfrac_exp in zip(mbexp, mfrac_mbexp):

458 varray = exp.variance.array

459 w = np.where(np.isfinite(varray) & (varray > 0))

460

461 if w[0].size == 0:

462 raise ValueError('no variance are finite')

463

464 var = np.median(exp.variance.array[w])

465 wgt = 1 / var

466 wgts.append(wgt)

467 wsum += wgt

468

469 if mfrac is None:

470 mfrac = wgt * mfrac_exp.image.array

471 else:

472 mfrac += wgt * mfrac_exp.image.array

473

474 mfrac *= 1.0 / wsum

475

476 return mfrac, wgts

477

478

479def fit_original_psfs_mbexp(mbexp, rng, wgts):

480 """

481 fit the original psfs at the center of the image and get the mean g1,g2,T

482 across all bands

483

484 This can fail and flags will be set, but we proceed

485 """

486

487 assert len(wgts) == len(mbexp)

488 wsum = sum(wgts)

489 if wsum <= 0:

490 raise ValueError(f'got sum(wgts) = {wsum}')

491

492 fitter = ngmix.admom.AdmomFitter(rng=rng)

493 guesser = ngmix.guessers.GMixPSFGuesser(

494 rng=rng,

495 ngauss=1,

496 guess_from_moms=True,

497 )

498 runner = ngmix.runners.PSFRunner(fitter=fitter, guesser=guesser, ntry=4)

499

500 try:

501 g1sum = 0.0

502 g2sum = 0.0

503 Tsum = 0.0

504

505 for exp, wgt in zip(mbexp, wgts):

506 cen, _ = get_integer_center(

507 wcs=exp.getWcs(),

508 bbox=exp.getBBox(),

509 as_double=True,

510 )

511 jac = get_jacobian(exp=exp, cen=cen)

512

513 psf_im = measure.extract_psf_image(exp, cen)

514

515 psf_cen = (np.array(psf_im.shape) - 1.0) / 2.0

516 psf_jacob = jac.copy()

517 psf_jacob.set_cen(row=psf_cen[0], col=psf_cen[1])

518

519 psf_obs = ngmix.Observation(

520 psf_im,

521 jacobian=psf_jacob,

522 )

523 res = runner.go(obs=psf_obs)

524 if res['flags'] != 0:

525 raise BootPSFFailure('failed to fit psf')

526

527 g1, g2 = res['e']

528 T = res['T']

529

530 g1sum += g1 * wgt

531 g2sum += g2 * wgt

532 Tsum += T * wgt

533

534 g1 = g1sum / wsum

535 g2 = g2sum / wsum

536 T = Tsum / wsum

537

538 flags = 0

539

540 except BootPSFFailure:

541 flags = procflags.PSF_FAILURE

542 g1 = -9999.0

543 g2 = -9999.0

544 T = -9999.0

545

546 return {

547 'flags': flags,

548 'g1': g1,

549 'g2': g2,

550 'T': T,

551 }

Coverage for metadetect / lsst / metadetect.py: 18%

192 statements