LSST Applications g2079a07aa2+86d27d4dc4,g2305ad1205+a659bff248,g2bbee38e9b+3c60f8fe34,g337abbeb29+3c60f8fe34,g33d1c0ed96+3c60f8fe34,g3502564af9+d77d6d1350,g3a166c0a6a+3c60f8fe34,g487adcacf7+25d9892218,g4be5004598+d77d6d1350,g50ff169b8f+96c6868917,g52b1c1532d+585e252eca,g591dd9f2cf+4d81263f9a,g5cd07815a0+980d2b1c3b,g607f77f49a+d77d6d1350,g858d7b2824+d77d6d1350,g88963caddf+83e433e629,g99cad8db69+a4d3c48eeb,g9ddcbc5298+9a081db1e4,ga1e77700b3+bcf1af89ad,ga57fefb910+9a39d7b2d7,gae0086650b+585e252eca,gb065fddaf9+4f9fd82a2c,gb0e22166c9+60f28cb32d,gb363559e06+d84b1d3d07,gb3b7280ab2+4563d032e1,gb4b16eec92+babe958938,gba4ed39666+c2a2e4ac27,gbb8dafda3b+ed6854b564,gc120e1dc64+b72d212f87,gc28159a63d+3c60f8fe34,gc3e9b769f7+921dbcd359,gcf0d15dbbd+9a39d7b2d7,gdaeeff99f8+f9a426f77a,gddc38dedce+585e252eca,ge79ae78c31+3c60f8fe34,w.2024.21
LSST Data Management Base Package
Loading...
Searching...
No Matches
dipoleMeasurement.py
Go to the documentation of this file.
1# This file is part of ip_diffim.
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
22import numpy as np
23
24import lsst.afw.image as afwImage
25import lsst.geom as geom
26import lsst.pex.config as pexConfig
27import lsst.meas.deblender.baseline as deblendBaseline
28from lsst.meas.base.pluginRegistry import register
29from lsst.meas.base import SingleFrameMeasurementTask, SingleFrameMeasurementConfig, \
30 SingleFramePluginConfig, SingleFramePlugin
31import lsst.afw.display as afwDisplay
32from lsst.utils.logging import getLogger
33
34__all__ = ("DipoleMeasurementConfig", "DipoleMeasurementTask", "DipoleAnalysis", "DipoleDeblender",
35 "SourceFlagChecker", "ClassificationDipoleConfig", "ClassificationDipolePlugin")
36
37
39 """Configuration for classification of detected diaSources as dipole or not"""
40 minSn = pexConfig.Field(
41 doc="Minimum quadrature sum of positive+negative lobe S/N to be considered a dipole",
42 dtype=float, default=np.sqrt(2) * 5.0,
43 )
44 maxFluxRatio = pexConfig.Field(
45 doc="Maximum flux ratio in either lobe to be considered a dipole",
46 dtype=float, default=0.65
47 )
48
49
50@register("ip_diffim_ClassificationDipole")
52 """A plugin to classify whether a diaSource is a dipole.
53 """
54
55 ConfigClass = ClassificationDipoleConfig
56
57 @classmethod
59 """
60 Returns
61 -------
62 result : `callable`
63 """
64 return cls.APCORR_ORDER
65
66 def __init__(self, config, name, schema, metadata):
67 SingleFramePlugin.__init__(self, config, name, schema, metadata)
69 self.keyProbability = schema.addField(name + "_value", type="D",
70 doc="Set to 1 for dipoles, else 0.")
71 self.keyFlag = schema.addField(name + "_flag", type="Flag", doc="Set to 1 for any fatal failure.")
72
73 def measure(self, measRecord, exposure):
74 passesSn = self.dipoleAnalysis.getSn(measRecord) > self.config.minSn
75 negFlux = np.abs(measRecord.get("ip_diffim_PsfDipoleFlux_neg_instFlux"))
76 negFluxFlag = measRecord.get("ip_diffim_PsfDipoleFlux_neg_flag")
77 posFlux = np.abs(measRecord.get("ip_diffim_PsfDipoleFlux_pos_instFlux"))
78 posFluxFlag = measRecord.get("ip_diffim_PsfDipoleFlux_pos_flag")
79
80 if negFluxFlag or posFluxFlag:
81 self.failfail(measRecord)
82 # continue on to classify
83
84 totalFlux = negFlux + posFlux
85
86 # If negFlux or posFlux are NaN, these evaluate to False
87 passesFluxNeg = (negFlux / totalFlux) < self.config.maxFluxRatio
88 passesFluxPos = (posFlux / totalFlux) < self.config.maxFluxRatio
89 if (passesSn and passesFluxPos and passesFluxNeg):
90 val = 1.0
91 else:
92 val = 0.0
93
94 measRecord.set(self.keyProbability, val)
95
96 def fail(self, measRecord, error=None):
97 measRecord.set(self.keyFlag, True)
98
99
101 """Measurement of detected diaSources as dipoles"""
102
103 def setDefaults(self):
104 SingleFrameMeasurementConfig.setDefaults(self)
105 self.pluginsplugins = ["base_CircularApertureFlux",
106 "base_PixelFlags",
107 "base_SkyCoord",
108 "base_PsfFlux",
109 "ip_diffim_PsfDipoleFlux",
110 "ip_diffim_ClassificationDipole",
111 ]
112
113 self.slots.calibFlux = None
114 self.slots.modelFlux = None
115 self.slots.gaussianFlux = None
116 self.slots.shape = None
117 self.slots.centroid = "ip_diffim_PsfDipoleFlux"
119
120
122 """Measurement of Sources, specifically ones from difference images, for characterization as dipoles
123
124 Parameters
125 ----------
126 sources : 'lsst.afw.table.SourceCatalog'
127 Sources that will be measured
128 badFlags : `list` of `dict`
129 A list of flags that will be used to determine if there was a measurement problem
130
131 """
132 ConfigClass = DipoleMeasurementConfig
133 _DefaultName = "dipoleMeasurement"
134
135
136
139
140class SourceFlagChecker(object):
141 """Functor class to check whether a diaSource has flags set that should cause it to be labeled bad."""
142
143 def __init__(self, sources, badFlags=None):
144 self.badFlags = ['base_PixelFlags_flag_edge', 'base_PixelFlags_flag_interpolatedCenter',
145 'base_PixelFlags_flag_saturatedCenter']
146 if badFlags is not None:
147 for flag in badFlags:
148 self.badFlags.append(flag)
149 self.keys = [sources.getSchema().find(name).key for name in self.badFlags]
150 self.keys.append(sources.table.getCentroidFlagKey())
151
152 def __call__(self, source):
153 """Call the source flag checker on a single Source
154
155 Parameters
156 ----------
157 source :
158 Source that will be examined
159 """
160 for k in self.keys:
161 if source.get(k):
162 return False
163 return True
164
165
166class DipoleAnalysis(object):
167 """Functor class that provides (S/N, position, orientation) of measured dipoles"""
168
169 def __init__(self):
170 pass
171
172 def __call__(self, source):
173 """Parse information returned from dipole measurement
174
175 Parameters
176 ----------
177 source : `lsst.afw.table.SourceRecord`
178 The source that will be examined"""
179 return self.getSn(source), self.getCentroid(source), self.getOrientation(source)
180
181 def getSn(self, source):
182 """Get the total signal-to-noise of the dipole; total S/N is from positive and negative lobe
183
184 Parameters
185 ----------
186 source : `lsst.afw.table.SourceRecord`
187 The source that will be examined"""
188
189 posflux = source.get("ip_diffim_PsfDipoleFlux_pos_instFlux")
190 posfluxErr = source.get("ip_diffim_PsfDipoleFlux_pos_instFluxErr")
191 negflux = source.get("ip_diffim_PsfDipoleFlux_neg_instFlux")
192 negfluxErr = source.get("ip_diffim_PsfDipoleFlux_neg_instFluxErr")
193
194 # Not a dipole!
195 if (posflux < 0) is (negflux < 0):
196 return 0
197
198 return np.sqrt((posflux/posfluxErr)**2 + (negflux/negfluxErr)**2)
199
200 def getCentroid(self, source):
201 """Get the centroid of the dipole; average of positive and negative lobe
202
203 Parameters
204 ----------
205 source : `lsst.afw.table.SourceRecord`
206 The source that will be examined"""
207
208 negCenX = source.get("ip_diffim_PsfDipoleFlux_neg_centroid_x")
209 negCenY = source.get("ip_diffim_PsfDipoleFlux_neg_centroid_y")
210 posCenX = source.get("ip_diffim_PsfDipoleFlux_pos_centroid_x")
211 posCenY = source.get("ip_diffim_PsfDipoleFlux_pos_centroid_y")
212 if (np.isinf(negCenX) or np.isinf(negCenY) or np.isinf(posCenX) or np.isinf(posCenY)):
213 return None
214
215 center = geom.Point2D(0.5*(negCenX+posCenX),
216 0.5*(negCenY+posCenY))
217 return center
218
219 def getOrientation(self, source):
220 """Calculate the orientation of dipole; vector from negative to positive lobe
221
222 Parameters
223 ----------
224 source : `lsst.afw.table.SourceRecord`
225 The source that will be examined"""
226
227 negCenX = source.get("ip_diffim_PsfDipoleFlux_neg_centroid_x")
228 negCenY = source.get("ip_diffim_PsfDipoleFlux_neg_centroid_y")
229 posCenX = source.get("ip_diffim_PsfDipoleFlux_pos_centroid_x")
230 posCenY = source.get("ip_diffim_PsfDipoleFlux_pos_centroid_y")
231 if (np.isinf(negCenX) or np.isinf(negCenY) or np.isinf(posCenX) or np.isinf(posCenY)):
232 return None
233
234 dx, dy = posCenX-negCenX, posCenY-negCenY
235 angle = geom.Angle(np.arctan2(dx, dy), geom.radians)
236 return angle
237
238 def displayDipoles(self, exposure, sources):
239 """Display debugging information on the detected dipoles
240
241 Parameters
242 ----------
243 exposure : `lsst.afw.image.Exposure`
244 Image the dipoles were measured on
245 sources : `lsst.afw.table.SourceCatalog`
246 The set of diaSources that were measured"""
247
248 import lsstDebug
249 display = lsstDebug.Info(__name__).display
250 displayDiaSources = lsstDebug.Info(__name__).displayDiaSources
251 maskTransparency = lsstDebug.Info(__name__).maskTransparency
252 if not maskTransparency:
253 maskTransparency = 90
254 disp = afwDisplay.Display(frame=lsstDebug.frame)
255 disp.setMaskTransparency(maskTransparency)
256 disp.mtv(exposure)
257
258 if display and displayDiaSources:
259 with disp.Buffering():
260 for source in sources:
261 cenX = source.get("ipdiffim_DipolePsfFlux_x")
262 cenY = source.get("ipdiffim_DipolePsfFlux_y")
263 if np.isinf(cenX) or np.isinf(cenY):
264 cenX, cenY = source.getCentroid()
265
266 isdipole = source.get("ip_diffim_ClassificationDipole_value")
267 if isdipole and np.isfinite(isdipole):
268 # Dipole
269 ctype = afwDisplay.GREEN
270 else:
271 # Not dipole
272 ctype = afwDisplay.RED
273
274 disp.dot("o", cenX, cenY, size=2, ctype=ctype)
275
276 negCenX = source.get("ip_diffim_PsfDipoleFlux_neg_centroid_x")
277 negCenY = source.get("ip_diffim_PsfDipoleFlux_neg_centroid_y")
278 posCenX = source.get("ip_diffim_PsfDipoleFlux_pos_centroid_x")
279 posCenY = source.get("ip_diffim_PsfDipoleFlux_pos_centroid_y")
280 if (np.isinf(negCenX) or np.isinf(negCenY) or np.isinf(posCenX) or np.isinf(posCenY)):
281 continue
282
283 disp.line([(negCenX, negCenY), (posCenX, posCenY)], ctype=afwDisplay.YELLOW)
284
285 lsstDebug.frame += 1
286
287
288class DipoleDeblender(object):
289 """Functor to deblend a source as a dipole, and return a new source with deblended footprints.
290
291 This necessarily overrides some of the functionality from
292 meas_algorithms/python/lsst/meas/algorithms/deblend.py since we
293 need a single source that contains the blended peaks, not
294 multiple children sources. This directly calls the core
295 deblending code deblendBaseline.deblend (optionally _fitPsf for
296 debugging).
297
298 Not actively being used, but there is a unit test for it in
299 dipoleAlgorithm.py.
300 """
301
302 def __init__(self):
303 # Set up defaults to send to deblender
304
305 # Always deblend as Psf
306 self.psfChisqCut1 = self.psfChisqCut2 = self.psfChisqCut2b = np.inf
307 self.log = getLogger('lsst.ip.diffim.DipoleDeblender')
308 self.sigma2fwhm = 2. * np.sqrt(2. * np.log(2.))
309
310 def __call__(self, source, exposure):
311 fp = source.getFootprint()
312 peaks = fp.getPeaks()
313 peaksF = [pk.getF() for pk in peaks]
314 fbb = fp.getBBox()
315 fmask = afwImage.Mask(fbb)
316 fmask.setXY0(fbb.getMinX(), fbb.getMinY())
317 fp.spans.setMask(fmask, 1)
318
319 psf = exposure.getPsf()
320 psfSigPix = psf.computeShape(psf.getAveragePosition()).getDeterminantRadius()
321 psfFwhmPix = psfSigPix * self.sigma2fwhm
322 subimage = afwImage.ExposureF(exposure, bbox=fbb, deep=True)
323 cpsf = deblendBaseline.CachingPsf(psf)
324
325 # if fewer than 2 peaks, just return a copy of the source
326 if len(peaks) < 2:
327 return source.getTable().copyRecord(source)
328
329 # make sure you only deblend 2 peaks; take the brighest and faintest
330 speaks = [(p.getPeakValue(), p) for p in peaks]
331 speaks.sort()
332 dpeaks = [speaks[0][1], speaks[-1][1]]
333
334 # and only set these peaks in the footprint (peaks is mutable)
335 peaks.clear()
336 for peak in dpeaks:
337 peaks.append(peak)
338
339 if True:
340 # Call top-level deblend task
341 fpres = deblendBaseline.deblend(fp, exposure.getMaskedImage(), psf, psfFwhmPix,
342 log=self.log,
343 psfChisqCut1=self.psfChisqCut1,
344 psfChisqCut2=self.psfChisqCut2,
345 psfChisqCut2b=self.psfChisqCut2b)
346 else:
347 # Call lower-level _fit_psf task
348
349 # Prepare results structure
350 fpres = deblendBaseline.DeblenderResult(fp, exposure.getMaskedImage(), psf, psfFwhmPix, self.log)
351
352 for pki, (pk, pkres, pkF) in enumerate(zip(dpeaks, fpres.deblendedParents[0].peaks, peaksF)):
353 self.log.debug('Peak %i', pki)
354 deblendBaseline._fitPsf(fp, fmask, pk, pkF, pkres, fbb, dpeaks, peaksF, self.log,
355 cpsf, psfFwhmPix,
356 subimage.image,
357 subimage.variance,
359
360 deblendedSource = source.getTable().copyRecord(source)
361 deblendedSource.setParent(source.getId())
362 peakList = deblendedSource.getFootprint().getPeaks()
363 peakList.clear()
364
365 for i, peak in enumerate(fpres.deblendedParents[0].peaks):
366 if peak.psfFitFlux > 0:
367 suffix = "pos"
368 else:
369 suffix = "neg"
370 c = peak.psfFitCenter
371 self.log.info("deblended.centroid.dipole.psf.%s %f %f",
372 suffix, c[0], c[1])
373 self.log.info("deblended.chi2dof.dipole.%s %f",
374 suffix, peak.psfFitChisq / peak.psfFitDof)
375 self.log.info("deblended.flux.dipole.psf.%s %f",
376 suffix, peak.psfFitFlux * np.sum(peak.templateImage.array))
377 peakList.append(peak.peak)
378 return deblendedSource
Represent a 2-dimensional array of bitmask pixels.
Definition Mask.h:77
A class representing an angle.
Definition Angle.h:128
fail(self, measRecord, error=None)