LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
quickFrameMeasurement.py
Go to the documentation of this file.
1# This file is part of pipe_tasks.
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
22__all__ = ["QuickFrameMeasurementTaskConfig", "QuickFrameMeasurementTask"]
23
24import numpy as np
25import scipy.ndimage as ndImage
26
27import lsst.afw.detection as afwDetect
28import lsst.afw.table as afwTable
29import lsst.geom as geom
30import lsst.meas.base as measBase
31import lsst.daf.base as dafBase
32import lsst.pipe.base as pipeBase
33import lsst.pex.config as pexConfig
34from lsst.meas.base import MeasurementError
35from lsst.meas.algorithms.installGaussianPsf import InstallGaussianPsfTask
36
37
38class QuickFrameMeasurementTaskConfig(pexConfig.Config):
39 """Config class for the QuickFrameMeasurementTask.
40 """
41 installPsf = pexConfig.ConfigurableField(
42 target=InstallGaussianPsfTask,
43 doc="Task for installing an initial PSF",
44 )
45 maxNonRoundness = pexConfig.Field(
46 dtype=float,
47 doc="Ratio of xx to yy (or vice versa) above which to cut, in order to exclude spectra",
48 default=5.,
49 )
50 maxExtendedness = pexConfig.Field(
51 dtype=float,
52 doc="Max absolute value of xx and yy above which to cut, in order to exclude large/things",
53 default=100,
54 )
55 doExtendednessCut = pexConfig.Field(
56 dtype=bool,
57 doc="Apply the extendeness cut, as definted by maxExtendedness",
58 default=False,
59 )
60 centroidPixelPercentile = pexConfig.Field(
61 dtype=float,
62 doc="The image's percentile value which the centroid must be greater than to pass the final peak"
63 " check. Ignored if doCheckCentroidPixelValue is False",
64 default=90,
65 )
66 doCheckCentroidPixelValue = pexConfig.Field(
67 dtype=bool,
68 doc="Check that the centroid found is actually in the centroidPixelPercentile percentile of the"
69 " image? Set to False for donut images.",
70 default=True,
71 )
72 initialPsfWidth = pexConfig.Field(
73 dtype=float,
74 doc="Guess at the initial PSF FWHM in pixels.",
75 default=10,
76 )
77 nSigmaDetection = pexConfig.Field(
78 dtype=float,
79 doc="Number of sigma for the detection limit.",
80 default=20,
81 )
82 nPixMinDetection = pexConfig.Field(
83 dtype=int,
84 doc="Minimum number of pixels in a detected source.",
85 default=10,
86 )
87 donutDiameter = pexConfig.Field(
88 dtype=int,
89 doc="The expected diameter of donuts in a donut image, in pixels.",
90 default=400,
91 )
92
93 def setDefaults(self):
94 super().setDefaults()
95 self.installPsf.fwhm = self.initialPsfWidth
96
97
98class QuickFrameMeasurementTask(pipeBase.Task):
99 """WARNING: An experimental new task with changable API! Do not rely on yet!
100
101 This task finds the centroid of the brightest source in a given CCD-image
102 and returns its centroid and a rough estimate of the seeing/PSF.
103
104 It is designed for speed, such that it can be used in observing scripts
105 to provide pointing offsets, allowing subsequent pointings to place
106 a source at an exact pixel position.
107
108 The approach taken here is deliberately sub-optimal in the detection and
109 measurement sense, with all optimisation being done for speed and robustness
110 of the result.
111
112 A small set of unit tests exist for this task, which run automatically
113 if afwdata is setup. These, however, are stricky unit tests, and will not
114 catch algorithmic regressions. TODO: DM-29038 exists to merge a regression
115 real test which runs against 1,000 LATISS images, but is therefore slow
116 and requires access to the data.
117
118 Parameters
119 ----------
120 config : `lsst.pipe.tasks.quickFrameMeasurement.QuickFrameMeasurementTaskConfig`
121 Configuration class for the QuickFrameMeasurementTask.
122 display : `lsst.afw.display.Display`, optional
123 The display to use for showing the images, detections and centroids.
124
125 Returns
126 -------
127 result : `lsst.pipe.base.Struct`
128 Return strucure containing whether the task was successful, the main
129 source's centroid, its the aperture fluxes, the ixx and iyy of the
130 source, and the median ixx, iyy of the detections in the exposure.
131 See run() method for further details.
132
133 Raises
134 ------
135 This task should *never* raise, as the run() method is enclosed in an
136 except Exception block, so that it will never fail during observing.
137 Failure modes should be limited to returning a return Struct() with the same
138 structure as the success case, with all value set to np.nan but with
139 result.success=False.
140 """
141 ConfigClass = QuickFrameMeasurementTaskConfig
142 _DefaultName = 'quickFrameMeasurementTask'
143
144 def __init__(self, config, *, display=None, **kwargs):
145 super().__init__(config=config, **kwargs)
146 self.makeSubtask("installPsf")
147
148 self.display = None
149 if display:
150 self.display = display
151
152 self.centroidName = "base_SdssCentroid"
153 self.shapeName = "base_SdssShape"
154 self.schema = afwTable.SourceTable.makeMinimalSchema()
155 self.schema.getAliasMap().set("slot_Centroid", self.centroidName)
156 self.schema.getAliasMap().set("slot_Shape", self.shapeName)
157 self.control = measBase.SdssCentroidControl()
158 self.control.maxDistToPeak = -1
159 self.centroider = measBase.SdssCentroidAlgorithm(self.control, self.centroidName, self.schema)
160 self.sdssShape = measBase.SdssShapeControl()
161 self.shaper = measBase.SdssShapeAlgorithm(self.sdssShape, self.shapeName, self.schema)
162 self.apFluxControl = measBase.ApertureFluxControl()
164 self.apFluxer = measBase.CircularApertureFluxAlgorithm(self.apFluxControl, "aperFlux",
165 self.schema, md)
166
167 self.table = afwTable.SourceTable.make(self.schema) # make sure to call this last!
168
169 @staticmethod
170 def detectObjectsInExp(exp, nSigma, nPixMin, grow=0):
171 """Run a very basic but fast threshold-based object detection on an exposure
172 Return the footPrintSet for the objects in a postISR exposure.
173
174 Parameters
175 ----------
176 exp : `lsst.afw.image.Exposure`
177 Image in which to detect objects.
178 nSigma : `float`
179 nSigma above image's stddev at which to set the detection threshold.
180 nPixMin : `int`
181 Minimum number of pixels for detection.
182 grow : `int`
183 Grow the detected footprint by this many pixels.
184
185 Returns
186 -------
187 footPrintSet : `lsst.afw.detection.FootprintSet`
188 FootprintSet containing the detections.
189 """
190 threshold = afwDetect.Threshold(nSigma, afwDetect.Threshold.STDEV)
191 footPrintSet = afwDetect.FootprintSet(exp.getMaskedImage(), threshold, "DETECTED", nPixMin)
192 if grow > 0:
193 isotropic = True
194 footPrintSet = afwDetect.FootprintSet(footPrintSet, grow, isotropic)
195 return footPrintSet
196
197 @staticmethod
198 def checkResult(exp, centroid, srcNum, percentile):
199 """Perform a final check that centroid location is actually bright.
200
201 Parameters
202 ----------
203 exp : `lsst.afw.image.Exposure`
204 The exposure on which to operate
205 centroid : `tuple` of `float`
206 Location of the centroid in pixel coordinates
207 scrNum : `int`
208 Number of the source in the source catalog. Only used if the check
209 is failed, for debug purposes.
210 percentile : `float`
211 Image's percentile above which the pixel containing the centroid
212 must be in order to pass the check.
213
214 Raises
215 ------
216 ValueError
217 Raised if the centroid's pixel is not above the percentile threshold
218 """
219 threshold = np.percentile(exp.image.array, percentile)
220 pixelValue = exp.image[centroid]
221 if pixelValue < threshold:
222 msg = (f"Final centroid pixel value check failed: srcNum {srcNum} at {centroid}"
223 f" has central pixel = {pixelValue:3f} <"
224 f" {percentile} percentile of image = {threshold:3f}")
225 raise ValueError(msg)
226 return
227
228 @staticmethod
229 def _calcMedianXxYy(objData):
230 """Return the median ixx and iyy for object in the image.
231 """
232 medianXx = np.nanmedian([element['xx'] for element in objData.values()])
233 medianYy = np.nanmedian([element['yy'] for element in objData.values()])
234 return medianXx, medianYy
235
236 @staticmethod
237 def _getCenterOfMass(exp, nominalCentroid, boxSize):
238 """Get the centre of mass around a point in the image.
239
240 Parameters
241 ----------
242 exp : `lsst.afw.image.Exposure`
243 The exposure in question.
244 nominalCentroid : `tuple` of `float`
245 Nominal location of the centroid in pixel coordinates.
246 boxSize : `int`
247 The size of the box around the nominalCentroid in which to measure
248 the centre of mass.
249
250 Returns
251 -------
252 com : `tuple` of `float`
253 The locaiton of the centre of mass of the brightest source in pixel
254 coordinates.
255 """
256 centroidPoint = geom.Point2I(nominalCentroid)
257 extent = geom.Extent2I(1, 1)
258 bbox = geom.Box2I(centroidPoint, extent)
259 bbox = bbox.dilatedBy(int(boxSize//2))
260 bbox = bbox.clippedTo(exp.getBBox())
261 data = exp[bbox].image.array
262 xy0 = exp[bbox].getXY0()
263
264 peak = ndImage.center_of_mass(data)
265 peak = (peak[1], peak[0]) # numpy coords returned
266 com = geom.Point2D(xy0)
267 com.shift(geom.Extent2D(*peak))
268 return (com[0], com[1])
269
270 def _calcBrightestObjSrcNum(self, objData):
271 """Find the brightest source which passes the cuts among the sources.
272
273 Parameters
274 ----------
275 objData : `dict` of `dict`
276 Dictionary, keyed by source number, containing the measurements.
277
278 Returns
279 -------
280 srcNum : `int`
281 The source number of the brightest source which passes the cuts.
282 """
283 max70, max70srcNum = -1, -1
284 max25, max25srcNum = -1, -1
285
286 for srcNum in sorted(objData.keys()): # srcNum not contiguous so don't use a list comp
287 # skip flag used rather than continue statements so we have all the
288 # metrics computed for debug purposes as this task is whack-a-mole
289 skip = False
290 xx = objData[srcNum]['xx']
291 yy = objData[srcNum]['yy']
292
293 xx = max(xx, 1e-9) # need to protect against division by zero
294 yy = max(yy, 1e-9) # because we don't `continue` on zero moments
295
296 if self.config.doExtendednessCut:
297 if xx > self.config.maxExtendedness or yy > self.config.maxExtendedness:
298 skip = True
299
300 nonRoundness = xx/yy
301 nonRoundness = max(nonRoundness, 1/nonRoundness)
302 if nonRoundness > self.config.maxNonRoundness:
303 skip = True
304
305 if self.log.isEnabledFor(self.log.DEBUG):
306 text = f"src {srcNum}: {objData[srcNum]['xCentroid']:.0f}, {objData[srcNum]['yCentroid']:.0f}"
307 text += f" - xx={xx:.1f}, yy={yy:.1f}, nonRound={nonRoundness:.1f}"
308 text += f" - ap70={objData[srcNum]['apFlux70']:,.0f}"
309 text += f" - ap25={objData[srcNum]['apFlux25']:,.0f}"
310 text += f" - skip={skip}"
311 self.log.debug(text)
312
313 if skip:
314 continue
315
316 ap70 = objData[srcNum]['apFlux70']
317 ap25 = objData[srcNum]['apFlux25']
318 if ap70 > max70:
319 max70 = ap70
320 max70srcNum = srcNum
321 if ap25 > max25:
322 max25 = ap25
323 max25srcNum = srcNum
324 if max70srcNum != max25srcNum:
325 self.log.warning("WARNING! Max apFlux70 for different object than with max apFlux25")
326
327 if max70srcNum >= 0: # starts as -1, return None if nothing is acceptable
328 return max70srcNum
329 return None
330
331 def _measureFp(self, fp, exp):
332 """Run the measurements on a footprint.
333
334 Parameters
335 ----------
336 fp : `lsst.afw.detection.Footprint`
337 The footprint to measure.
338 exp : `lsst.afw.image.Exposure`
339 The footprint's parent exposure.
340
341 Returns
342 -------
343 src : `lsst.afw.table.SourceRecord`
344 The source record containing the measurements.
345 """
346 src = self.table.makeRecord()
347 src.setFootprint(fp)
348 self.centroider.measure(src, exp)
349 self.shaper.measure(src, exp)
350 self.apFluxer.measure(src, exp)
351 return src
352
353 def _getDataFromSrcRecord(self, src):
354 """Extract the shapes and centroids from a source record.
355
356 Parameters
357 ----------
358 src : `lsst.afw.table.SourceRecord`
359 The source record from which to extract the measurements.
360
361 Returns
362 -------
363 srcData : `lsst.pipe.base.Struct`
364 The struct containing the extracted measurements.
365 """
366 pScale = self.plateScale
367 xx = np.sqrt(src['base_SdssShape_xx'])*2.355*pScale # 2.355 for FWHM, pScale for platescale from exp
368 yy = np.sqrt(src['base_SdssShape_yy'])*2.355*pScale
369 xCentroid = src['base_SdssCentroid_x']
370 yCentroid = src['base_SdssCentroid_y']
371 # apFluxes available: 70, 50, 35, 25, 17, 12 9, 6, 4.5, 3
372 apFlux70 = src['aperFlux_70_0_instFlux']
373 apFlux25 = src['aperFlux_25_0_instFlux']
374 return pipeBase.Struct(xx=xx,
375 yy=yy,
376 xCentroid=xCentroid,
377 yCentroid=yCentroid,
378 apFlux70=apFlux70,
379 apFlux25=apFlux25)
380
381 @staticmethod
383 """Get the shape, centroid and flux from a footprint.
384
385 Parameters
386 ----------
387 fp : `lsst.afw.detection.Footprint`
388 The footprint to measure.
389 exp : `lsst.afw.image.Exposure`
390 The footprint's parent exposure.
391
392 Returns
393 -------
394 srcData : `lsst.pipe.base.Struct`
395 The struct containing the extracted measurements.
396 """
397 xx = fp.getShape().getIxx()
398 yy = fp.getShape().getIyy()
399 xCentroid, yCentroid = fp.getCentroid()
400 apFlux70 = np.sum(exp[fp.getBBox()].image.array)
401 apFlux25 = np.sum(exp[fp.getBBox()].image.array)
402 return pipeBase.Struct(xx=xx,
403 yy=yy,
404 xCentroid=xCentroid,
405 yCentroid=yCentroid,
406 apFlux70=apFlux70,
407 apFlux25=apFlux25)
408
409 @staticmethod
410 def _measurementResultToDict(measurementResult):
411 """Convenience function to repackage measurement results to a dict.
412
413 Parameters
414 ----------
415 measurementResult : `lsst.afw.table.SourceRecord`
416 The source record to convert to a dict.
417
418 Returns
419 -------
420 objData : `dict`
421 The dict containing the extracted data.
422 """
423 objData = {}
424 objData['xx'] = measurementResult.xx
425 objData['yy'] = measurementResult.yy
426 objData['xCentroid'] = measurementResult.xCentroid
427 objData['yCentroid'] = measurementResult.yCentroid
428 objData['apFlux70'] = measurementResult.apFlux70
429 objData['apFlux25'] = measurementResult.apFlux25
430 return objData
431
432 @staticmethod
434 """Make the default/template return struct, with defaults to False/nan.
435
436 Returns
437 -------
438 objData : `lsst.pipe.base.Struct`
439 The default template return structure.
440 """
441 result = pipeBase.Struct()
442 result.success = False
443 result.brightestObjCentroid = (np.nan, np.nan)
444 result.brightestObjCentroidCofM = None
445 result.brightestObj_xXyY = (np.nan, np.nan)
446 result.brightestObjApFlux70 = np.nan
447 result.brightestObjApFlux25 = np.nan
448 result.medianXxYy = (np.nan, np.nan)
449 return result
450
451 def run(self, exp, *, donutDiameter=None, doDisplay=False):
452 """Calculate position, flux and shape of the brightest star in an image.
453
454 Given an an assembled (and at least minimally ISRed exposure),
455 quickly and robustly calculate the centroid of the
456 brightest star in the image.
457
458 Parameters
459 ----------
460 exp : `lsst.afw.image.Exposure`
461 The exposure in which to find and measure the brightest star.
462 donutDiameter : `int` or `float`, optional
463 The expected diameter of donuts in pixels for use in the centre of
464 mass centroid measurement. If None is provided, the config option
465 is used.
466 doDisplay : `bool`
467 Display the image and found sources. A diplay object must have
468 been passed to the task constructor.
469
470 Returns
471 -------
472 result : `lsst.pipe.base.Struct`
473 Struct containing:
474 Whether the task ran successfully and found the object (bool)
475 The object's centroid (float, float)
476 The object's ixx, iyy (float, float)
477 The object's 70 pixel aperture flux (float)
478 The object's 25 pixel aperture flux (float)
479 The images's median ixx, iyy (float, float)
480 If unsuccessful, the success field is False and all other results
481 are np.nan of the expected shape.
482
483 Notes
484 -----
485 Because of this task's involvement in observing scripts, the run method
486 should *never* raise. Failure modes are noted by returning a Struct with
487 the same structure as the success case, with all value set to np.nan and
488 result.success=False.
489 """
490 try:
491 result = self._run(exp=exp, donutDiameter=donutDiameter, doDisplay=doDisplay)
492 return result
493 except Exception as e:
494 self.log.warning("Failed to find main source centroid %s", e)
495 result = self._makeEmptyReturnStruct()
496 return result
497
498 def _run(self, exp, *, donutDiameter=None, doDisplay=False):
499 """The actual run method, called by run()
500
501 Behaviour is documented in detail in the main run().
502 """
503 if donutDiameter is None:
504 donutDiameter = self.config.donutDiameter
505
506 self.plateScale = exp.getWcs().getPixelScale().asArcseconds()
507 median = np.nanmedian(exp.image.array)
508 exp.image -= median # is put back later
509 self.installPsf.run(exp)
510 sources = self.detectObjectsInExp(exp, nSigma=self.config.nSigmaDetection,
511 nPixMin=self.config.nPixMinDetection)
512
513 if doDisplay:
514 if self.display is None:
515 raise RuntimeError("Display failed as no display provided during init()")
516 self.display.mtv(exp)
517
518 fpSet = sources.getFootprints()
519 self.log.info("Found %d sources in exposure", len(fpSet))
520
521 objData = {}
522 nMeasured = 0
523
524 for srcNum, fp in enumerate(fpSet):
525 try:
526 src = self._measureFp(fp, exp)
527 result = self._getDataFromSrcRecord(src)
528 except MeasurementError:
529 try:
530 # gets shape and centroid from footprint
531 result = self._getDataFromFootprintOnly(fp, exp)
532 except MeasurementError as e:
533 self.log.info("Skipped measuring source %s: %s", srcNum, e)
534 continue
535 objData[srcNum] = self._measurementResultToDict(result)
536 nMeasured += 1
537
538 self.log.info("Measured %d of %d sources in exposure", nMeasured, len(fpSet))
539
540 medianXxYy = self._calcMedianXxYy(objData)
541
542 brightestObjSrcNum = self._calcBrightestObjSrcNum(objData)
543 if brightestObjSrcNum is None:
544 raise RuntimeError("No sources in image passed cuts")
545
546 x = objData[brightestObjSrcNum]['xCentroid']
547 y = objData[brightestObjSrcNum]['yCentroid']
548 brightestObjCentroid = (x, y)
549 xx = objData[brightestObjSrcNum]['xx']
550 yy = objData[brightestObjSrcNum]['yy']
551 brightestObjApFlux70 = objData[brightestObjSrcNum]['apFlux70']
552 brightestObjApFlux25 = objData[brightestObjSrcNum]['apFlux25']
553
554 exp.image += median # put background back in
555 if self.config.doCheckCentroidPixelValue:
556 self.checkResult(exp, brightestObjCentroid, brightestObjSrcNum,
557 self.config.centroidPixelPercentile)
558
559 boxSize = donutDiameter * 1.3 # allow some slack, as cutting off side of donut is very bad
560 centreOfMass = self._getCenterOfMass(exp, brightestObjCentroid, boxSize)
561
562 result = self._makeEmptyReturnStruct()
563 result.success = True
564 result.brightestObjCentroid = brightestObjCentroid
565 result.brightestObj_xXyY = (xx, yy)
566 result.brightestObjApFlux70 = brightestObjApFlux70
567 result.brightestObjApFlux25 = brightestObjApFlux25
568 result.medianXxYy = medianXxYy
569 result.brightestObjCentroidCofM = centreOfMass
570
571 return result
int max
Class for storing generic metadata.
Definition PropertySet.h:66
An integer coordinate rectangle.
Definition Box.h:55
_run(self, exp, *donutDiameter=None, doDisplay=False)
daf::base::PropertySet * set
Definition fits.cc:931