Coverage for metadetect / lsst / measure.py: 18%
301 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-21 00:03 -0800
1import logging
2import warnings
3import numpy as np
4import esutil as eu
5import ngmix
7import lsst.afw.table as afw_table
8from lsst.meas.algorithms import SourceDetectionTask
9from lsst.meas.algorithms import SourceDetectionConfig as OriginalSourceDetectionConfig
10from lsst.meas.deblender import SourceDeblendTask, SourceDeblendConfig
11from lsst.meas.base import (
12 SingleFrameMeasurementConfig,
13 SingleFrameMeasurementTask,
14)
15from lsst.pex.config import Config, ConfigurableField, Field
16from lsst.pipe.base import Task
17import lsst.afw.image as afw_image
18import lsst.geom as geom
19from lsst.pex.exceptions import (
20 InvalidParameterError,
21 LengthError,
22)
24from ..procflags import (
25 EDGE_HIT, ZERO_WEIGHTS, CENTROID_FAILURE, NO_ATTEMPT,
26)
27from ..fitting import (
28 fit_mbobs_gauss,
29 fit_mbobs_wavg,
30 get_wavg_output_struct,
31)
33from . import util
34from .util import ContextNoiseReplacer, get_jacobian
35from . import vis
36from .defaults import DEFAULT_THRESH
38warnings.filterwarnings('ignore', category=FutureWarning)
40LOG = logging.getLogger('lsst_measure')
43class SourceDetectionConfig(OriginalSourceDetectionConfig):
44 @property
45 def thresh(self):
46 return self.thresholdValue
48 @thresh.setter
49 def thresh(self, value):
50 self.thresholdValue = value
53SourceDetectionTask.ConfigClass = SourceDetectionConfig
56class DetectAndDeblendConfig(Config):
57 meas = ConfigurableField[SingleFrameMeasurementConfig](
58 doc="Measurement config",
59 target=SingleFrameMeasurementTask,
60 )
62 detect = ConfigurableField[SourceDetectionConfig](
63 doc="Detection config",
64 target=SourceDetectionTask
65 )
67 deblend = ConfigurableField[SourceDeblendConfig](
68 doc="Deblend config",
69 target=SourceDeblendTask
70 )
72 seed = Field[int](
73 doc="Random rng seed",
74 default=42,
75 )
77 def setDefaults(self):
78 super().setDefaults()
80 # defaults for measurement config
81 self.meas.plugins.names = [
82 "base_SdssCentroid",
83 "base_PsfFlux",
84 "base_SkyCoord",
85 ]
87 # set these slots to none because we aren't running these algorithms
88 self.meas.slots.apFlux = None
89 self.meas.slots.gaussianFlux = None
90 self.meas.slots.calibFlux = None
91 self.meas.slots.modelFlux = None
93 # goes with SdssShape above
94 self.meas.slots.shape = None
96 # fix odd issue where it things things are near the edge
97 self.meas.plugins['base_SdssCentroid'].binmax = 1
99 # defaults for detection config
100 # DM does not have config default stability. Set all of them explicitly
101 self.detect.minPixels = 1
102 self.detect.isotropicGrow = True
103 self.detect.combinedGrow = True
104 self.detect.nSigmaToGrow = 2.4
105 self.detect.returnOriginalFootprints = False
106 self.detect.includeThresholdMultiplier = 1.0
107 self.detect.thresholdPolarity = "positive"
108 self.detect.adjustBackground = 0.0
109 self.detect.doApplyFlatBackgroundRatio = False
110 # these are ignored since we are doing reEstimateBackground = False
111 # detection_config.background
112 # detection_config.tempLocalBackground
113 # detection_config.doTempLocalBackground
114 # detection_config.tempWideBackground
115 # detection_config.doTempWideBackground
116 self.detect.nPeaksMaxSimple = 1
117 self.detect.nSigmaForKernel = 7.0
118 self.detect.excludeMaskPlanes = util.get_detection_mask()
120 # the defaults changed from from stdev to pixel_std but
121 # we don't want that
123 self.detect.thresholdType = "stdev"
124 # our changes from defaults
125 self.detect.reEstimateBackground = False
127 self.detect.thresholdValue = DEFAULT_THRESH
129 # these will be ignored when finding the image standard deviation
130 self.detect.statsMask = util.get_stats_mask()
132 # deblend config
133 # these tasks must use the same schema and all be constructed before any
134 # other tasks using the same schema are run because schema is modified in
135 # place by tasks, and the constructor does a check that fails if we do this
136 # afterward
137 self.deblend.maxFootprintArea = 0
140class DetectAndDeblendTask(Task):
141 ConfigClass = DetectAndDeblendConfig
142 _DefaultName = "detect_and_deblend"
144 def __init__(self, *args, **kwargs):
145 super().__init__(*args, **kwargs)
146 schema = afw_table.SourceTable.makeMinimalSchema()
147 self.makeSubtask("meas", schema=schema)
148 self.makeSubtask("detect", schema=schema)
149 self.makeSubtask("deblend", schema=schema)
150 self.rng = np.random.RandomState(seed=self.config.seed)
152 def run(self, mbexp, show=False):
153 if len(mbexp.singles) > 1:
154 detexp = util.coadd_exposures(mbexp.singles)
155 if detexp is None:
156 return [], None
157 else:
158 detexp = mbexp.singles[0]
160 # background measurement within the detection code requires ExposureF
161 if not isinstance(detexp, afw_image.ExposureF):
162 detexp = afw_image.ExposureF(detexp, deep=True)
164 schema = self.deblend.schema # should be the same for all tasks
165 table = afw_table.SourceTable.make(schema)
166 result = self.detect.run(table, detexp)
168 if result is not None:
169 sources = result.sources
170 self.deblend.run(detexp, sources)
172 with ContextNoiseReplacer(
173 detexp, sources, self.rng, config=self.meas.config.noiseReplacer
174 ) as replacer:
176 for source in sources:
178 if source.get('deblend_nChild') != 0:
179 continue
181 source_id = source.getId()
183 with replacer.sourceInserted(source_id):
184 self.meas.callMeasure(source, detexp)
186 else:
187 sources = []
189 if show:
190 vis.show_exp(detexp, use_mpl=True, sources=sources)
192 return sources, detexp
195def detect_and_deblend(
196 mbexp,
197 rng=None,
198 thresh=DEFAULT_THRESH,
199 show=False,
200 config=None,
201):
202 """
203 run detection and deblending of peaks, as well as basic measurments such as
204 centroid. The SDSS deblender is run in order to split footprints.
206 We must combine detection and deblending in the same function because the
207 schema gets modified in place, which means we must construct the deblend
208 task at the same time as the detect task
210 Parameters
211 ----------
212 mbexp: lsst.afw.image.MultibandExposure
213 The exposures to process
214 rng: np.random.RandomState
215 Random number generator for noise replacer
216 thresh: float, optional
217 The detection threshold in units of the sky noise
218 show: bool, optional
219 If set to True, show images
221 Returns
222 -------
223 sources, detexp
224 The sources and the detection exposure
225 """
226 config_override = config if config is not None else {}
227 if thresh:
228 if 'detect' not in config_override:
229 config_override['detect'] = {}
230 config_override['detect']['thresholdValue'] = thresh
232 config = DetectAndDeblendConfig()
233 config.setDefaults()
235 util.override_config(config, config_override)
237 config.freeze()
238 config.validate()
239 task = DetectAndDeblendTask(config=config)
240 if rng is not None:
241 task.rng = rng
242 return task.run(mbexp, show)
245def measure(
246 mbexp,
247 detexp,
248 sources,
249 fitter,
250 stamp_size,
251 fwhm_reg=0,
252 rng=None,
253):
254 """
255 run measurements on the input exposure, given the input measurement task,
256 list of sources, and fitter.
258 Parameters
259 ----------
260 mbexp: lsst.afw.image.MultibandExposure
261 The exposure on which to detect and measure
262 detexp: lsst.afw.image.Exposure*
263 The detection exposure, used for bmask info
264 sources: list of sources
265 From a detection task
266 fitter: e.g. ngmix.gaussmom.GaussMom or ngmix.ksigmamom.PGaussMom
267 For calculating moments. Can also be string 'gauss'
268 stamp_size: int
269 Size for postage stamps
270 fwhm_reg: float, optional
271 Optional regularization for calculating shapes. The fwhm is converted
272 to T and T+Treg is used in the denominator
273 rng: np.random.RandomState
274 Random number generator, only used when fitter == 'gauss'
276 Returns
277 -------
278 ndarray with results or None
279 """
281 if len(sources) == 0:
282 return None
284 if fitter == 'gauss':
285 fitter_kind = fitter
286 else:
287 fitter_kind = fitter.kind
289 nband = len(mbexp.bands)
290 exp_bbox = mbexp.getBBox()
291 wcs = mbexp.singles[0].getWcs()
292 results = []
294 # bmasks will be different within the loop below due to the replacer
295 bmasks = get_bmasks(sources=sources, exposure=detexp)
297 for i, source in enumerate(sources):
299 if source.get('deblend_nChild') != 0:
300 continue
302 bmask = bmasks[i]
304 flags = 0
305 try:
306 mbobs = _get_stamp_mbobs(
307 mbexp=mbexp, source=source, stamp_size=stamp_size,
308 )
310 # TODO do something with bmask_flags?
311 if fitter == 'gauss':
312 psf_fitter = ngmix.admom.AdmomFitter(rng=rng)
313 psf_guesser = ngmix.guessers.GMixPSFGuesser(
314 rng=rng, ngauss=1, guess_from_moms=True,
315 )
316 psf_runner = ngmix.runners.PSFRunner(
317 fitter=psf_fitter, guesser=psf_guesser, ntry=4,
318 )
320 this_res = fit_mbobs_gauss(
321 mbobs=mbobs,
322 bmask_flags=0,
323 psf_runner=psf_runner,
324 rng=rng,
325 )
326 else:
327 this_res = fit_mbobs_wavg(
328 mbobs=mbobs,
329 fitter=fitter,
330 bmask_flags=0,
331 fwhm_reg=fwhm_reg,
332 )
333 except LengthError as err:
334 # This is raised when a bbox hits an edge
335 LOG.debug('%s', err)
336 flags = EDGE_HIT
337 except AllZeroWeightError as err:
338 # failure creating some observation due to zero weights
339 LOG.info('%s', err)
340 flags = ZERO_WEIGHTS
341 except CentroidFailError as err:
342 # failure in the center finding
343 LOG.info(str(err))
344 flags = CENTROID_FAILURE
346 if flags != 0:
347 this_res = get_wavg_output_struct(nband=nband, model=fitter_kind)
348 this_res[fitter_kind + '_flags'] = flags
350 res = get_output(
351 wcs=wcs, source=source, res=this_res,
352 bmask=bmask, stamp_size=stamp_size, exp_bbox=exp_bbox,
353 )
355 results.append(res)
357 if len(results) > 0:
358 results = eu.numpy_util.combine_arrlist(results)
359 else:
360 results = None
362 return results
365def get_bmasks(sources, exposure):
366 """
367 get a list of all the bmasks for the sources
369 Parameters
370 ----------
371 sources: lsst.afw.table.SourceCatalog
372 The sources
373 exposure: lsst.afw.image.ExposureF
374 The exposure
376 Returns
377 -------
378 list of bmask values
379 """
380 bmasks = []
381 for source in sources:
382 bmask = get_bmask(source=source, exposure=exposure)
383 bmasks.append(bmask)
384 return bmasks
387def get_bmask(source, exposure):
388 """
389 get bmask based on original peak position
391 Parameters
392 ----------
393 sources: lsst.afw.table.SourceRecord
394 The sources
395 exposure: lsst.afw.image.ExposureF
396 The exposure
398 Returns
399 -------
400 bmask value
401 """
402 peak = source.getFootprint().getPeaks()[0]
403 orig_cen = peak.getI()
404 maskval = exposure.mask[orig_cen]
405 return maskval
408def extract_obs(exp, source):
409 """
410 convert an image object into an ngmix.Observation, including
411 a psf observation
413 parameters
414 ----------
415 imobj: lsst.afw.image.ExposureF
416 The exposure
417 source: lsst.afw.table.SourceRecord
418 The source record
420 returns
421 --------
422 obs: ngmix.Observation
423 The Observation unless all the weight are zero, in which
424 case AllZeroWeightError is raised
425 """
427 im = exp.image.array
429 wt = _extract_weight(exp)
430 if np.all(wt <= 0):
431 raise AllZeroWeightError('all weights <= 0')
433 bmask = exp.mask.array
434 jacob = _extract_jacobian_at_source(
435 exp=exp,
436 source=source,
437 )
439 orig_cen = source.getCentroid()
441 psf_im = extract_psf_image(exposure=exp, orig_cen=orig_cen)
443 # fake the psf pixel noise
444 psf_err = psf_im.max()*0.0001
445 psf_wt = psf_im*0 + 1.0/psf_err**2
447 # use canonical center for the psf
448 psf_cen = (np.array(psf_im.shape)-1.0)/2.0
449 psf_jacob = jacob.copy()
450 psf_jacob.set_cen(row=psf_cen[0], col=psf_cen[1])
452 # we will have need of the bit names which we can only
453 # get from the mask object
454 # this is sort of monkey patching, but I'm not sure of
455 # a better solution
457 meta = {'orig_cen': orig_cen}
459 psf_obs = ngmix.Observation(
460 psf_im,
461 weight=psf_wt,
462 jacobian=psf_jacob,
463 )
464 obs = ngmix.Observation(
465 im,
466 weight=wt,
467 bmask=bmask,
468 jacobian=jacob,
469 psf=psf_obs,
470 meta=meta,
471 )
473 return obs
476def _get_stamp_mbobs(mbexp, source, stamp_size, clip=False):
477 """
478 Get a postage stamp MultibandExposure
480 Parameters
481 ----------
482 mbexp: lsst.afw.image.MultibandExposure
483 The exposures
484 source: lsst.afw.table.SourceRecord
485 The source for which to get the stamp
486 stamp_size: int
487 If sent, a bounding box is created with about this size rather than
488 using the footprint bounding box. Typically the returned size is
489 stamp_size + 1
490 clip: bool, optional
491 If set to True, clip the bbox to fit into the exposure.
493 If clip is False and the bbox does not fit, a
494 lsst.pex.exceptions.LengthError is raised
496 Only relevant if stamp_size is sent. Default False
498 Returns
499 -------
500 lsst.afw.image.ExposureF
501 """
503 bbox = _get_bbox(mbexp, source, stamp_size, clip=clip)
505 mbobs = ngmix.MultiBandObsList()
506 for band in mbexp.bands:
508 subexp = mbexp[band][bbox]
509 obs = extract_obs(
510 exp=subexp,
511 source=source,
512 )
514 obslist = ngmix.ObsList(meta={'band': band})
515 obslist.append(obs)
516 mbobs.append(obslist)
518 return mbobs
521def _get_bbox(mbexp, source, stamp_size, clip=False):
522 """
523 Get a bounding box at the location of the specified source.
525 Parameters
526 ----------
527 mbexp: lsst.afw.image.MultibandExposure
528 The exposures
529 source: lsst.afw.table.SourceRecord
530 The source for which to get the stamp
531 stamp_size: int
532 If sent, a bounding box is created with about this size rather than
533 using the footprint bounding box. Typically the returned size is
534 stamp_size + 1
535 clip: bool, optional
536 If set to True, clip the bbox to fit into the exposure.
538 If clip is False and the bbox does not fit, a
539 lsst.pex.exceptions.LengthError is raised
541 Only relevant if stamp_size is sent. Default False
543 Returns
544 -------
545 lsst.geom.Box2I
546 """
548 fp = source.getFootprint()
549 peak = fp.getPeaks()[0]
551 x_peak, y_peak = peak.getIx(), peak.getIy()
553 bbox = geom.Box2I(
554 geom.Point2I(x_peak, y_peak),
555 geom.Extent2I(1, 1),
556 )
557 bbox.grow(stamp_size // 2)
559 exp_bbox = mbexp.getBBox()
560 if clip:
561 bbox.clip(exp_bbox)
562 else:
563 if not exp_bbox.contains(bbox):
564 source_id = source.getId()
565 raise LengthError(
566 f'requested stamp size {stamp_size} for source '
567 f'{source_id} does not fit into the exposoure. '
568 f'Use clip=True to clip the bbox to fit'
569 )
571 return bbox
574def extract_psf_image(exposure, orig_cen):
575 """
576 get the psf associated with this image.
578 coadded psfs from DM are generally not square, but the coadd in cells code
579 makes them so. We will assert they are square and odd dimensions
581 Parameters
582 ----------
583 exposure: lsst.afw.image.ExposureF
584 The exposure data
585 orig_cen: lsst.geom.Point2D
586 The location at which to draw the image
588 Returns
589 -------
590 ndarray
591 """
592 try:
593 psfobj = exposure.getPsf()
594 psfim = psfobj.computeKernelImage(orig_cen).array
595 except InvalidParameterError:
596 raise MissingDataError("could not reconstruct PSF")
598 psfim = np.asarray(psfim, dtype='f4')
600 shape = psfim.shape
601 assert shape[0] == shape[1], 'require square psf images'
602 assert shape[0] % 2 != 0, 'require odd psf images'
604 return psfim
607def _extract_psf_image_fix(exposure, orig_cen):
608 """
609 get the psf associated with this image
611 coadded psfs are generally not square, so we will
612 trim it to be square and preserve the center to
613 be at the new canonical center
615 TODO: should we really trim the psf to be even? will this
616 cause a shift due being off-center?
617 """
618 try:
619 psfobj = exposure.getPsf()
620 psfim = psfobj.computeKernelImage(orig_cen).array
621 except InvalidParameterError:
622 raise MissingDataError("could not reconstruct PSF")
624 psfim = np.asarray(psfim, dtype='f4')
626 psfim = util.trim_odd_image(psfim)
627 return psfim
630def _extract_weight(exp):
631 """
632 TODO get the estimated sky variance rather than this hack
633 TODO should we zero out other bits?
635 extract a weight map
637 Areas with NO_DATA will get zero weight.
639 Because the variance map includes the full poisson variance, which
640 depends on the signal, we instead extract the median of the parts of
641 the image without NO_DATA set
643 parameters
644 ----------
645 exp: sub exposure object
646 """
648 # TODO implement bit checking
649 var_image = exp.variance.array
651 weight = var_image.copy()
653 weight[:, :] = 0
655 wuse = np.where(var_image > 0)
657 if wuse[0].size > 0:
658 medvar = np.median(var_image[wuse])
659 weight[:, :] = 1.0/medvar
660 else:
661 print(' weight is all zero, found '
662 'none that passed cuts')
664 return weight
667def _extract_jacobian_at_source(exp, source):
668 """
669 extract an ngmix.Jacobian from the image object
670 and object record
672 exp: lsst.afw.image.ExposureF
673 An exposure object
674 source: lsst.afw.table.SourceRecord
675 The source record created during detection
677 returns
678 --------
679 Jacobian: ngmix.Jacobian
680 The local jacobian
681 """
683 orig_cen = exp.getWcs().skyToPixel(source.getCoord())
685 if np.isnan(orig_cen.getY()):
686 LOG.info('falling back on integer location')
687 # fall back to integer pixel location
688 peak = source.getFootprint().getPeaks()[0]
689 orig_cen_i = peak.getI()
690 orig_cen = geom.Point2D(
691 x=orig_cen_i.getX(),
692 y=orig_cen_i.getY(),
693 )
695 return get_jacobian(exp, orig_cen)
698def get_output_dtype():
700 dt = [
701 ('stamp_size', 'i4'),
702 ('row0', 'i4'), # bbox row start
703 ('col0', 'i4'), # bbox col start
704 ('row', 'f4'), # row in image. Use row0 to get to global pixel coords
705 ('col', 'f4'), # col in image. Use col0 to get to global pixel coords
706 ('row_diff', 'f4'), # difference from peak location
707 ('col_diff', 'f4'), # difference from peak location
708 ('ra', 'f8'),
709 ('dec', 'f8'),
711 ('psfrec_flags', 'i4'), # psfrec is the original psf
712 ('psfrec_g', 'f8', 2),
713 ('psfrec_T', 'f8'),
715 # values from .mask of input exposures
716 ('bmask', 'i4'),
717 # values for ormask across all input exposures to coadd
718 ('ormask', 'i4'),
719 # fraction of images going into a pixel that were masked
720 ('mfrac', 'f4'),
721 ]
723 return dt
726def get_output_struct(res):
727 """
728 get the output struct
730 Parameters
731 ----------
732 res: ndarray
733 The result from running metadetect.fitting.fit_mbobs_wavg
735 Returns
736 -------
737 ndarray
738 Has the fields from res, with new fields added, see get_output_dtype
739 """
740 dt = get_output_dtype()
741 output = eu.numpy_util.add_fields(res, dt)
743 for subdt in dt:
744 name = subdt[0]
745 dtype = subdt[1]
747 if 'flags' in name:
748 output[name] = NO_ATTEMPT
749 elif name in ('bmask', 'ormask'):
750 output[name] = 0
751 elif dtype[0] == 'i':
752 output[name] = -9999
753 else:
754 output[name] = np.nan
756 return output
759def get_output(wcs, source, res, bmask, stamp_size, exp_bbox):
760 """
761 get the output structure, copying in results
763 The following fields are not set:
764 psfrec_flags, psfrec_g, psfrec_T
765 mfrac
767 Parameters
768 ----------
769 wcs: a stack wcs
770 The wcs with which to determine the ra, dec
771 res: ndarray
772 The result from running metadetect.fitting.fit_mbobs_wavg
773 bmask: int
774 The bmask value at the location of this object
775 stamp_size: int
776 The stamp size used for the measurement
777 exp_bbox: lsst.geom.Box2I
778 The bounding box used for measurement
780 Returns
781 -------
782 ndarray
783 Has the fields from res, with new fields added, see get_output_dtype
784 """
785 output = get_output_struct(res)
787 orig_cen = source.getCentroid()
789 skypos = wcs.pixelToSky(orig_cen)
791 peak = source.getFootprint().getPeaks()[0]
792 peak_loc = peak.getI()
794 if np.isnan(orig_cen.getY()):
795 orig_cen = peak.getCentroid()
796 cen_offset = geom.Point2D(np.nan, np.nan)
797 else:
798 cen_offset = geom.Point2D(
799 orig_cen.getX() - peak_loc.getX(),
800 orig_cen.getY() - peak_loc.getY(),
801 )
803 output['stamp_size'] = stamp_size
804 output['row0'] = exp_bbox.getBeginY()
805 output['col0'] = exp_bbox.getBeginX()
806 output['row'] = orig_cen.getY()
807 output['col'] = orig_cen.getX()
808 output['row_diff'] = cen_offset.getY()
809 output['col_diff'] = cen_offset.getX()
811 output['ra'] = skypos.getRa().asDegrees()
812 output['dec'] = skypos.getDec().asDegrees()
814 # remove DETECTED bit, it is just clutter since all detected
815 # objects have this bit set
816 detected = afw_image.Mask.getPlaneBitMask('DETECTED')
817 output['bmask'] = bmask & ~detected
819 return output
822class MissingDataError(Exception):
823 """
824 Some number was out of range
825 """
827 def __init__(self, value):
828 super().__init__(value)
829 self.value = value
831 def __str__(self):
832 return repr(self.value)
835class AllZeroWeightError(Exception):
836 """
837 Some number was out of range
838 """
840 def __init__(self, value):
841 super().__init__(value)
842 self.value = value
844 def __str__(self):
845 return repr(self.value)
848class CentroidFailError(Exception):
849 """
850 Some number was out of range
851 """
853 def __init__(self, value):
854 super().__init__(value)
855 self.value = value
857 def __str__(self):
858 return repr(self.value)