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

1import logging

2import numpy as np

3import ngmix

4from ngmix.gexceptions import BootPSFFailure

5from lsst.pex.config import (

6 Config,

7 ChoiceField,

8 ConfigField,

9 ConfigurableField,

10 Field,

11 FieldValidationError,

12 ListField,

13)

14from lsst.pipe.base import Task

15from lsst.meas.algorithms import SourceDetectionTask

17from .. import procflags

19from .skysub import subtract_sky_mbexp

21from .defaults import (

22 DEFAULT_FWHM_REG,

23 DEFAULT_FWHM_SMOOTH,

24 DEFAULT_STAMP_SIZES,

25 DEFAULT_SUBTRACT_SKY,

26 DEFAULT_WEIGHT_FWHMS,

27)

28from . import measure

29from .metacal_exposures import get_metacal_mbexps_fixnoise

30from .util import get_integer_center, get_jacobian, override_config

32LOG = logging.getLogger('lsst_metadetect')

35def run_metadetect(

36 mbexp, noise_mbexp, rng, mfrac_mbexp=None, ormasks=None, config=None, show=False,

37):

38 """

39 Run metadetection on the input MultiBandObsList

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

43 Parameters

44 ----------

45 mbexp: lsst.afw.image.MultibandExposure

46 The exposures to process

47 noise_mbexp: lsst.afw.image.MultibandExposure

48 The noise exposures for metacal

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

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

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

52 ormasks: list of images, optional

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

54 into a coadd.

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

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

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

59 areas such as bright objects or apodized edges.

61 In the future we expect the MultibandExposure to have an ormask attribute

62 rng: np.random.RandomState

63 Random number generator

64 config: dict, optional

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

66 in this dict override defaults; see lsst_configs.py

67 show: bool, optional

68 if set to True images will be shown

70 Returns

71 -------

72 result dict

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

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

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

76 """

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

79 config = MetadetectConfig()

80 config.setDefaults()

82 override_config(config, config_override)

84 config.freeze()

85 config.validate()

86 task = MetadetectTask(config=config)

87 result = task.run(mbexp, noise_mbexp, rng, mfrac_mbexp, ormasks, show=show,)

88 return result

91class WeightConfig(Config):

92 fwhm = Field[float](

93 doc="FWHM of the Gaussian weight function (in arcseconds)",

94 default=None,

95 )

96 fwhm_smooth = Field[float](

97 doc="FWHM of the Gaussian smoothing function (in arcseconds)",

98 default=DEFAULT_FWHM_SMOOTH,

99 )

100 fwhm_reg = Field[float](

101 doc="FWHM of the Gaussian regularization function (in arcseconds)",

102 default=DEFAULT_FWHM_REG,

103 )

104

105

106class LmParseConfig(Config):

107 maxfev = Field[int](

108 doc="Maximum number of function evaluations",

109 default=4000,

110 )

111 ftol = Field[float](

112 doc="Relative error desired in the sum of squares",

113 default=1.0e-5,

114 )

115 xtol = Field[float](

116 doc="Relative error desired in the approximate solution",

117 default=1.0e-5,

118 )

119

120

121class PsfConfig(Config):

122 """Config class for fitting original PSFs."""

123 model = ChoiceField[str](

124 doc="Model for fitting original PSFs",

125 default="am",

126 allowed={

127 "am": "Adaptive moments", # TODO: What are the other possible values?

128 "gauss": "Gaussian",

129 },

130 optional=False,

131 )

132 ntry = Field[int](

133 doc="Number of tries for fitting original PSFs",

134 default=4,

135 )

136 lm_pars = ConfigField[LmParseConfig](

137 doc="Config for Levenberg-Marquardt fitting",

138 )

139

140

141class MetacalConfig(Config):

142 use_noise_image = Field[bool](

143 doc="Whether to use the noise image",

144 default=True,

145 )

146 psf = ChoiceField[str](

147 doc="PSF model",

148 default="fitgauss",

149 allowed={

150 "fitgauss": "Fit a Gaussian to the PSF",

151 }

152 )

153 types = ListField[str](

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

155 default=["noshear", "1p", "1m",],

156 )

157

158 def validate(self):

159 super().validate()

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

161 raise FieldValidationError(

162 self.__class__.types,

163 self,

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

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

166 )

167

168

169class MetadetectConfig(Config):

170 subtract_sky = Field[bool](

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

172 default=DEFAULT_SUBTRACT_SKY,

173 )

174 meas_type = ChoiceField[str](

175 doc="Measurement type",

176 default="wmom",

177 allowed={

178 "wmom": "Weighted moments",

179 "ksigma": "Fourier moments", # TODO: improve docstring

180 "pgauss": "Pre-seeing moments",

181 "am": "Adaptive moments",

182 "gauss": "Gaussian Fit",

183 }

184 )

185

186 weight = ConfigField[WeightConfig](

187 doc="FWHM of the Gaussian weight function (in arcseconds)",

188 )

189

190 psf = ConfigField[PsfConfig](

191 doc="Config for fitting the original PSFs",

192 )

193

194 detect = ConfigurableField(

195 doc="Detection config",

196 target=SourceDetectionTask,

197 )

198 metacal = ConfigField[MetacalConfig](

199 doc="Config for metacal",

200 )

201

202 @property

203 def stamp_size(self):

204 # Caution: This provides a backdoor entry to change the behavior

205 # by changing the DEFAULT_STAMP_SIZES.

206 return DEFAULT_STAMP_SIZES[self.meas_type]

207

208 def setDefaults(self):

209 super().setDefaults()

210 self.weight.fwhm = DEFAULT_WEIGHT_FWHMS.get(self.meas_type, None)

211

212 def validate(self):

213 super().validate()

214

215 if self.stamp_size != DEFAULT_STAMP_SIZES[self.meas_type]:

216 raise FieldValidationError(

217 self.__class__.stamp_size,

218 self,

219 f"stamp_size must be {DEFAULT_STAMP_SIZES[self.meas_type]} "

220 f"for meas_type {self.meas_type}",

221 )

222

223

224class MetadetectTask(Task):

225 ConfigClass = MetadetectConfig

226 _DefaultName = "metadetect"

227

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

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

230 self.makeSubtask("detect")

231

232 def run(

233 self,

234 mbexp,

235 noise_mbexp,

236 rng,

237 mfrac_mbexp=None,

238 ormasks=None,

239 show=False,

240 ):

241

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

243 config = self.config.toDict()

244 # Because this is a property and not a Field, we set this explicitly.

245 config['stamp_size'] = self.config.stamp_size

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

247

248 ormask = combine_ormasks(mbexp, ormasks)

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

250

251 if self.config.subtract_sky:

252 subtract_sky_mbexp(mbexp=mbexp, thresh=self.config.detect.thresholdValue)

253

254 psf_stats = fit_original_psfs_mbexp(

255 mbexp=mbexp,

256 wgts=wgts,

257 rng=rng,

258 )

259

260 fitter = get_fitter(config, rng=rng)

261

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

263

264 mdict, noise_mdict = get_metacal_mbexps_fixnoise(

265 mbexp=mbexp, noise_mbexp=noise_mbexp, types=metacal_types,

266 )

267

268 result = {}

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

270

271 res = detect_deblend_and_measure(

272 mbexp=mcal_mbexp,

273 fitter=fitter,

274 config=config,

275 rng=rng,

276 show=show,

277 )

278

279 if res is not None:

280 band = mcal_mbexp.bands[0]

281 exp = mcal_mbexp[band]

282

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

284 add_ormask(ormask, res)

285 add_original_psf(psf_stats, res)

286

287 result[shear_str] = res

288

289 return result

290

291

292def detect_deblend_and_measure(

293 mbexp,

294 fitter,

295 config,

296 rng,

297 show=False,

298):

299 """

300 run detection, deblending and measurements.

301

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

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

304 no other way to split footprints

305

306 Parameters

307 ----------

308 mbexp: lsst.afw.image.MultibandExposure

309 The metacal'ed exposures to process

310 fitter: e.g. ngmix.gaussmom.GaussMom or ngmix.ksigmamom.KSigmaMom

311 For calculating moments

312 config: dict, optional

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

314 in this dict override defaults; see lsst_configs.py

315 thresh: float

316 The detection threshold in units of the sky noise

317 stamp_size: int

318 Size for postage stamps.

319 show: bool, optional

320 If set to True, show images during processing

321 """

322

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

324

325 sources, detexp = measure.detect_and_deblend(

326 mbexp=mbexp,

327 rng=rng,

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

329 show=show,

330 )

331

332 if config['weight'] is not None:

333 fwhm_reg = config['weight'].get('fwhm_reg', 0)

334 else:

335 fwhm_reg = 0

336

337 results = measure.measure(

338 mbexp=mbexp,

339 detexp=detexp,

340 sources=sources,

341 fitter=fitter,

342 stamp_size=config['stamp_size'],

343 fwhm_reg=fwhm_reg,

344 rng=rng,

345 )

346

347 return results

348

349

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

351 """

352 calculate and add mfrac to the input result array

353 """

354 if np.any(mfrac > 0):

355

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

357 # this.

358

359 cen, _ = get_integer_center(

360 wcs=exp.getWcs(),

361 bbox=exp.getBBox(),

362 as_double=True,

363 )

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

365

366 res['mfrac'] = measure_weighted_mfrac(

367 mfrac=mfrac,

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

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

370 jac=jac,

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

372 )

373 else:

374 res['mfrac'] = 0

375

376

377def measure_weighted_mfrac(

378 *,

379 mfrac,

380 x,

381 y,

382 jac,

383 fwhm,

384):

385 """

386 Measure a Gaussian-weighted average of an image.

387

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

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

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

391

392 Parameters

393 ----------

394 mfrac : np.ndarray

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

396 x : np.ndarray

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

398 weighted average.

399 y : np.ndarray

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

401 weighted average.

402 box_sizes : np.ndarray

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

404 big enough to hold 2 * `fwhm`.

405 jac: ngmix.Jacobian

406 The jacobian for the data

407 fwhm : float or None

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

409 of 1.2 is used.

410

411 Returns

412 -------

413 mfracs : np.ndarray

414 The weighted averages at each input location.

415 """

416

417 if fwhm is None:

418 fwhm = 1.2

419

420 ny, nx = mfrac.shape

421

422 gauss_wgt = ngmix.GMixModel(

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

424 'gauss',

425 )

426 sigma = ngmix.moments.fwhm_to_sigma(fwhm)

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

428

429 mfracs = []

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

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

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

433

434 xstart = ix - box_rad

435 xend = ix + box_rad + 1

436 ystart = iy - box_rad

437 yend = iy + box_rad + 1

438

439 if xstart < 0:

440 xstart = 0

441 if ystart < 0:

442 ystart = 0

443 if xend > nx:

444 xend = nx

445 if yend > ny:

446 yend = ny

447

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

449 if sub_mfrac.size == 0:

450 mfracs.append(1.0)

451 else:

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

453 this_jac = jac.copy()

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

455

456 obs = ngmix.Observation(

457 image=sub_mfrac,

458 jacobian=this_jac,

459 )

460

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

462

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

464 # Gaussian as the weight.

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

466

467 return np.array(mfracs)

468

469

470def add_ormask(ormask, res):

471 """

472 copy in ormask values using the row, col positions

473 """

474 for i in range(res.size):

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

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

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

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

479

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

481

482

483def add_original_psf(psf_stats, res):

484 """

485 copy in psf results

486 """

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

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

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

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

491

492

493def get_fitter(config, rng=None):

494 """

495 get the fitter based on the 'fitter' input

496 """

497

498 meas_type = config['meas_type']

499

500 if meas_type == 'am':

501 fitter_obj = ngmix.admom.AdmomFitter(rng=rng)

502 guesser = ngmix.guessers.GMixPSFGuesser(

503 rng=rng, ngauss=1, guess_from_moms=True,

504 )

505 fitter = ngmix.runners.Runner(

506 fitter=fitter_obj, guesser=guesser,

507 ntry=2,

508 )

509 # TODO is there a better way?

510 fitter.kind = 'am'

511

512 elif meas_type == 'gauss':

513 fitter = 'gauss'

514 else:

515 fwhm = config['weight']['fwhm']

516

517 if meas_type == 'wmom':

518 fitter = ngmix.gaussmom.GaussMom(fwhm=fwhm)

519 else:

520 # for backwards compatibility with older ngmix that did not

521 # support the fwhm_smooth keyword.

522 sm_kwargs = {}

523 if config['weight']['fwhm_smooth'] > 0:

524 sm_kwargs["fwhm_smooth"] = config['weight']['fwhm_smooth']

525

526 if meas_type == 'ksigma':

527 fitter = ngmix.ksigmamom.KSigmaMom(

528 fwhm=fwhm,

529 **sm_kwargs

530 )

531 elif meas_type == 'pgauss':

532 fitter = ngmix.prepsfmom.PGaussMom(

533 fwhm=fwhm,

534 **sm_kwargs

535 )

536 else:

537 raise ValueError("bad meas_type: '%s'" % meas_type)

538

539 return fitter

540

541

542def combine_ormasks(mbexp, ormasks):

543 """

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

545 versions for each band

546 """

547 if ormasks is None:

548 bands = mbexp.bands

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

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

551 else:

552 for imask, tormask in enumerate(ormasks):

553 if imask == 0:

554 ormask = tormask.copy()

555 else:

556 ormask |= tormask

557

558 return ormask

559

560

561def get_mfrac_mbexp(mbexp, mfrac_mbexp):

562 """

563 set the masked fraction image, averaged over all bands

564

565 Parameters

566 ----------

567 mbexp: lsst.afw.image.MultibandExposure

568 The exposures to process

569

570 Returns

571 -------

572 mfrac: array

573 """

574 wgts = []

575

576 wsum = 0.0

577

578 mfrac = None

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

580 varray = exp.variance.array

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

582

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

584 raise ValueError('no variance are finite')

585

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

587 wgt = 1/var

588 wgts.append(wgt)

589 wsum += wgt

590

591 if mfrac is None:

592 mfrac = wgt * mfrac_exp.image.array

593 else:

594 mfrac += wgt * mfrac_exp.image.array

595

596 mfrac *= 1.0/wsum

597

598 return mfrac, wgts

599

600

601def fit_original_psfs_mbexp(mbexp, rng, wgts):

602 """

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

604 across all bands

605

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

607 """

608

609 assert len(wgts) == len(mbexp)

610 wsum = sum(wgts)

611 if wsum <= 0:

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

613

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

615 guesser = ngmix.guessers.GMixPSFGuesser(

616 rng=rng, ngauss=1, guess_from_moms=True,

617 )

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

619

620 try:

621 g1sum = 0.0

622 g2sum = 0.0

623 Tsum = 0.0

624

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

626 cen, _ = get_integer_center(

627 wcs=exp.getWcs(),

628 bbox=exp.getBBox(),

629 as_double=True,

630 )

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

632

633 psf_im = measure.extract_psf_image(exp, cen)

634

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

636 psf_jacob = jac.copy()

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

638

639 psf_obs = ngmix.Observation(

640 psf_im,

641 jacobian=psf_jacob,

642 )

643 res = runner.go(obs=psf_obs)

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

645 raise BootPSFFailure('failed to fit psf')

646

647 g1, g2 = res['e']

648 T = res['T']

649

650 g1sum += g1*wgt

651 g2sum += g2*wgt

652 Tsum += T*wgt

653

654 g1 = g1sum/wsum

655 g2 = g2sum/wsum

656 T = Tsum/wsum

657

658 flags = 0

659

660 except BootPSFFailure:

661 flags = procflags.PSF_FAILURE

662 g1 = -9999.0

663 g2 = -9999.0

664 T = -9999.0

665

666 return {

667 'flags': flags,

668 'g1': g1,

669 'g2': g2,

670 'T': T,

671 }

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

236 statements