Loading [MathJax]/extensions/tex2jax.js
LSST Applications g0fba68d861+c32a75abc2,g1ec0fe41b4+3ea9d11450,g1fd858c14a+3f2e8fb4ee,g2440f9efcc+8c5ae1fdc5,g2d6b4b5a25+53e3e6ec14,g35bb328faa+8c5ae1fdc5,g4a4af6cd76+de8e3777d4,g4d2262a081+2ada2482a1,g53246c7159+8c5ae1fdc5,g55585698de+c901811295,g56a49b3a55+f42f167f24,g60b5630c4e+c901811295,g67b6fd64d1+3fc8cb0b9e,g78460c75b0+7e33a9eb6d,g786e29fd12+668abc6043,g8352419a5c+8c5ae1fdc5,g8852436030+23a2163ad3,g89139ef638+3fc8cb0b9e,g94187f82dc+c901811295,g989de1cb63+3fc8cb0b9e,g9d31334357+c901811295,g9f33ca652e+90b8cadf48,gabe3b4be73+8856018cbb,gabf8522325+977d9fabaf,gae5fe98123+65ba0d7b48,gb1101e3267+d2a4d34c5e,gb89ab40317+3fc8cb0b9e,gc91f06edcd+2f65d60c96,gcf25f946ba+23a2163ad3,gd6cbbdb0b4+1cc2750d2e,gdb1c4ca869+be65c9c1d7,gde0f65d7ad+b3e1ecbd0f,ge278dab8ac+6b863515ed,ge410e46f29+3fc8cb0b9e,gf35d7ec915+97dd712d81,gf5e32f922b+8c5ae1fdc5,gf618743f1b+c6bc2b68fc,gf67bdafdda+3fc8cb0b9e,w.2025.18
LSST Data Management Base Package
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
objectSizeStarSelector.py
Go to the documentation of this file.
1# This file is part of meas_algorithms.
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__ = ["ObjectSizeStarSelectorConfig", "ObjectSizeStarSelectorTask",
23 "ObjectSizeNoSourcesError", "ObjectSizeNoGoodSourcesError"]
24
25import sys
26
27import numpy
28import warnings
29from functools import reduce
30
31from lsst.utils.logging import getLogger
32from lsst.pipe.base import Struct, AlgorithmError
33import lsst.geom
34from lsst.afw.cameraGeom import PIXELS, TAN_PIXELS
35import lsst.afw.geom as afwGeom
36import lsst.pex.config as pexConfig
37import lsst.afw.display as afwDisplay
38from .sourceSelector import BaseSourceSelectorTask, sourceSelectorRegistry
39
40afwDisplay.setDefaultMaskTransparency(75)
41
42_LOG = getLogger(__name__)
43
44
45class ObjectSizeNoSourcesError(AlgorithmError):
46 """Raised if no sources were supplied in the input catalog."""
47
48 def __init__(self) -> None:
49 super().__init__("Input catalog for source selection is empty.")
50
51 @property
52 def metadata(self) -> dict:
53 return {}
54
55
56class ObjectSizeNoGoodSourcesError(AlgorithmError):
57 """Raised if no good sources were available after applying configured
58 selection criteria.
59 """
60
61 def __init__(self, n_input_sources) -> None:
62 self._n_available = n_input_sources
63 super().__init__(f"No good star candidates out of {n_input_sources} input sources.")
64
65 @property
66 def metadata(self) -> dict:
67 return {
68 "n_input_sources": self._n_available,
69 }
70
71
72class ObjectSizeStarSelectorConfig(BaseSourceSelectorTask.ConfigClass):
73 doFluxLimit = pexConfig.Field(
74 doc="Apply flux limit to Psf Candidate selection?",
75 dtype=bool,
76 default=False,
77 )
78 fluxMin = pexConfig.Field(
79 doc="Minimum flux value for good Psf Candidates.",
80 dtype=float,
81 default=12500.0,
82 check=lambda x: x >= 0.0,
83 )
84 fluxMax = pexConfig.Field(
85 doc="Maximum flux value for good Psf Candidates (ignored if == 0).",
86 dtype=float,
87 default=0.0,
88 check=lambda x: x >= 0.0,
89 )
90 doSignalToNoiseLimit = pexConfig.Field(
91 doc="Apply signal-to-noise (i.e. flux/fluxErr) limit to Psf Candidate selection?",
92 dtype=bool,
93 default=True,
94 )
95 # Note that the current default is conditioned on the detection thresholds
96 # set in the characterizeImage setDefaults function for the measurePsf
97 # stage.
98 signalToNoiseMin = pexConfig.Field(
99 doc="Minimum signal-to-noise for good Psf Candidates "
100 "(value should take into consideration the detection thresholds "
101 "set for the catalog of interest).",
102 dtype=float,
103 default=50.0,
104 check=lambda x: x >= 0.0,
105 )
106 signalToNoiseMax = pexConfig.Field(
107 doc="Maximum signal-to-noise for good Psf Candidates (ignored if == 0).",
108 dtype=float,
109 default=0.0,
110 check=lambda x: x >= 0.0,
111 )
112 widthMin = pexConfig.Field(
113 doc="Minimum width to include in histogram.",
114 dtype=float,
115 default=0.9,
116 check=lambda x: x >= 0.0,
117 )
118 widthMax = pexConfig.Field(
119 doc="Maximum width to include in histogram.",
120 dtype=float,
121 default=10.0,
122 check=lambda x: x >= 0.0,
123 )
124 sourceFluxField = pexConfig.Field(
125 doc="Name of field in Source to use for flux measurement.",
126 dtype=str,
127 default="base_PsfFlux_instFlux",
128 )
129 widthStdAllowed = pexConfig.Field(
130 doc="Standard deviation of width allowed to be interpreted as good stars.",
131 dtype=float,
132 default=0.15,
133 check=lambda x: x >= 0.0,
134 )
135 nSigmaClip = pexConfig.Field(
136 doc="Keep objects within this many sigma of cluster 0's median.",
137 dtype=float,
138 default=2.0,
139 check=lambda x: x >= 0.0,
140 )
141 badFlags = pexConfig.ListField(
142 doc="List of flags which cause a source to be rejected as bad.",
143 dtype=str,
144 default=[
145 "base_PixelFlags_flag_edge",
146 "base_PixelFlags_flag_nodata",
147 "base_PixelFlags_flag_interpolatedCenter",
148 "base_PixelFlags_flag_saturatedCenter",
149 "base_PixelFlags_flag_crCenter",
150 "base_PixelFlags_flag_bad",
151 "base_PixelFlags_flag_interpolated",
152 "slot_Centroid_flag",
153 ],
154 )
155
156 def validate(self):
157 BaseSourceSelectorTask.ConfigClass.validate(self)
158 if self.widthMin > self.widthMax:
159 msg = f"widthMin ({self.widthMin}) > widthMax ({self.widthMax})"
160 raise pexConfig.FieldValidationError(ObjectSizeStarSelectorConfig.widthMin, self, msg)
161
162
164 """A class to handle key strokes with matplotlib displays.
165 """
166
167 def __init__(self, axes, xs, ys, x, y, frames=[0]):
168 self.axes = axes
169 self.xs = xs
170 self.ys = ys
171 self.x = x
172 self.y = y
173 self.frames = frames
174
175 self.cid = self.axes.figure.canvas.mpl_connect('key_press_event', self)
176
177 def __call__(self, ev):
178 if ev.inaxes != self.axes:
179 return
180
181 if ev.key and ev.key in ("p"):
182 dist = numpy.hypot(self.xs - ev.xdata, self.ys - ev.ydata)
183 dist[numpy.where(numpy.isnan(dist))] = 1e30
184
185 which = numpy.where(dist == min(dist))
186
187 x = self.x[which][0]
188 y = self.y[which][0]
189 for frame in self.frames:
190 disp = afwDisplay.Display(frame=frame)
191 disp.pan(x, y)
192 disp.flush()
193 else:
194 pass
195
196
197def _assignClusters(yvec, centers):
198 """Return a vector of centerIds based on their distance to the centers.
199 """
200 assert len(centers) > 0
201
202 minDist = numpy.nan*numpy.ones_like(yvec)
203 clusterId = numpy.empty_like(yvec)
204 clusterId.dtype = int # zeros_like(..., dtype=int) isn't in numpy 1.5
205 dbl = _LOG.getChild("_assignClusters")
206 dbl.setLevel(dbl.INFO)
207
208 # Make sure we are logging aall numpy warnings...
209 oldSettings = numpy.seterr(all="warn")
210 with warnings.catch_warnings(record=True) as w:
211 warnings.simplefilter("always")
212 for i, mean in enumerate(centers):
213 dist = abs(yvec - mean)
214 if i == 0:
215 update = dist == dist # True for all points
216 else:
217 update = dist < minDist
218 if w: # Only do if w is not empty i.e. contains a warning message
219 dbl.trace(str(w[-1]))
220
221 minDist[update] = dist[update]
222 clusterId[update] = i
223 numpy.seterr(**oldSettings)
224
225 return clusterId
226
227
228def _kcenters(yvec, nCluster, useMedian=False, widthStdAllowed=0.15):
229 """A classic k-means algorithm, clustering yvec into nCluster clusters
230
231 Return the set of centres, and the cluster ID for each of the points
232
233 If useMedian is true, use the median of the cluster as its centre, rather than
234 the traditional mean
235
236 Serge Monkewitz points out that there other (maybe smarter) ways of seeding the means:
237 "e.g. why not use the Forgy or random partition initialization methods"
238 however, the approach adopted here seems to work well for the particular sorts of things
239 we're clustering in this application
240 """
241
242 assert nCluster > 0
243
244 mean0 = sorted(yvec)[len(yvec)//10] # guess
245 delta = mean0 * widthStdAllowed * 2.0
246 centers = mean0 + delta * numpy.arange(nCluster)
247
248 func = numpy.median if useMedian else numpy.mean
249
250 clusterId = numpy.zeros_like(yvec) - 1 # which cluster the points are assigned to
251 clusterId.dtype = int # zeros_like(..., dtype=int) isn't in numpy 1.5
252 while True:
253 oclusterId = clusterId
254 clusterId = _assignClusters(yvec, centers)
255
256 if numpy.all(clusterId == oclusterId):
257 break
258
259 for i in range(nCluster):
260 # Only compute func if some points are available; otherwise, default to NaN.
261 pointsInCluster = (clusterId == i)
262 if numpy.any(pointsInCluster):
263 centers[i] = func(yvec[pointsInCluster])
264 else:
265 centers[i] = numpy.nan
266
267 return centers, clusterId
268
269
270def _improveCluster(yvec, centers, clusterId, nsigma=2.0, nIteration=10, clusterNum=0, widthStdAllowed=0.15):
271 """Improve our estimate of one of the clusters (clusterNum) by sigma-clipping around its median.
272 """
273
274 nMember = sum(clusterId == clusterNum)
275 if nMember < 5: # can't compute meaningful interquartile range, so no chance of improvement
276 return clusterId
277 for iter in range(nIteration):
278 old_nMember = nMember
279
280 inCluster0 = clusterId == clusterNum
281 yv = yvec[inCluster0]
282
283 centers[clusterNum] = numpy.median(yv)
284 stdev = numpy.std(yv)
285
286 syv = sorted(yv)
287 stdev_iqr = 0.741*(syv[int(0.75*nMember)] - syv[int(0.25*nMember)])
288 median = syv[int(0.5*nMember)]
289
290 sd = stdev if stdev < stdev_iqr else stdev_iqr
291
292 if False:
293 print("sigma(iqr) = %.3f, sigma = %.3f" % (stdev_iqr, numpy.std(yv)))
294 newCluster0 = abs(yvec - centers[clusterNum]) < nsigma*sd
295 clusterId[numpy.logical_and(inCluster0, newCluster0)] = clusterNum
296 clusterId[numpy.logical_and(inCluster0, numpy.logical_not(newCluster0))] = -1
297
298 nMember = sum(clusterId == clusterNum)
299 # 'sd < widthStdAllowed * median' prevents too much rejections
300 if nMember == old_nMember or sd < widthStdAllowed * median:
301 break
302
303 return clusterId
304
305
306def plot(mag, width, centers, clusterId, marker="o", markersize=2, markeredgewidth=0, ltype='-',
307 magType="model", clear=True):
308
309 log = _LOG.getChild("plot")
310 try:
311 import matplotlib.pyplot as plt
312 except ImportError as e:
313 log.warning("Unable to import matplotlib: %s", e)
314 return
315
316 try:
317 fig
318 except NameError:
319 fig = plt.figure()
320 else:
321 if clear:
322 fig.clf()
323
324 axes = fig.add_axes((0.1, 0.1, 0.85, 0.80))
325
326 xmin = sorted(mag)[int(0.05*len(mag))]
327 xmax = sorted(mag)[int(0.95*len(mag))]
328
329 axes.set_xlim(-17.5, -13)
330 axes.set_xlim(xmin - 0.1*(xmax - xmin), xmax + 0.1*(xmax - xmin))
331 axes.set_ylim(0, 10)
332
333 colors = ["r", "g", "b", "c", "m", "k", ]
334 for k, mean in enumerate(centers):
335 if k == 0:
336 axes.plot(axes.get_xlim(), (mean, mean,), "k%s" % ltype)
337
338 li = (clusterId == k)
339 axes.plot(mag[li], width[li], marker, markersize=markersize, markeredgewidth=markeredgewidth,
340 color=colors[k % len(colors)])
341
342 li = (clusterId == -1)
343 axes.plot(mag[li], width[li], marker, markersize=markersize, markeredgewidth=markeredgewidth,
344 color='k')
345
346 if clear:
347 axes.set_xlabel("Instrumental %s mag" % magType)
348 axes.set_ylabel(r"$\sqrt{(I_{xx} + I_{yy})/2}$")
349
350 return fig
351
352
353@pexConfig.registerConfigurable("objectSize", sourceSelectorRegistry)
355 """A star selector that looks for a cluster of small objects in a
356 size-magnitude plot, typically to select stars for PSF estimation.
357 """
358
359 ConfigClass = ObjectSizeStarSelectorConfig
360 usesMatches = False # selectStars does not use its matches argument
361
362 def selectSources(self, sourceCat, matches=None, exposure=None):
363 """Return a selection of PSF candidates that represent likely stars.
364
365 A list of PSF candidates may be used by a PSF fitter to construct a PSF.
366
367 Parameters
368 ----------
369 sourceCat : `lsst.afw.table.SourceCatalog`
370 Catalog of sources to select from.
371 This catalog must be contiguous in memory.
372 matches : `list` of `lsst.afw.table.ReferenceMatch` or None
373 Ignored in this SourceSelector.
374 exposure : `lsst.afw.image.Exposure` or None
375 The exposure the catalog was built from; used to get the detector
376 to transform to TanPix, and for debug display.
377
378 Returns
379 -------
380 struct : `lsst.pipe.base.Struct`
381 The struct contains the following data:
382
383 ``selected``
384 Boolean array of sources that were selected, same length as
385 sourceCat. (`numpy.ndarray` of `bool`)
386 """
387 if len(sourceCat) == 0:
388 raise ObjectSizeNoSourcesError
389
390 import lsstDebug
391 display = lsstDebug.Info(__name__).display
392 displayExposure = lsstDebug.Info(__name__).displayExposure # display the Exposure + spatialCells
393 plotMagSize = lsstDebug.Info(__name__).plotMagSize # display the magnitude-size relation
394 dumpData = lsstDebug.Info(__name__).dumpData # dump data to pickle file?
395
396 detector = None
397 pixToTanPix = None
398 if exposure:
399 detector = exposure.getDetector()
400 if detector:
401 pixToTanPix = detector.getTransform(PIXELS, TAN_PIXELS)
402 #
403 # Look at the distribution of stars in the magnitude-size plane
404 #
405 flux = sourceCat[self.config.sourceFluxField]
406 fluxErr = sourceCat[self.config.sourceFluxField + "Err"]
407
408 xx = numpy.empty(len(sourceCat))
409 xy = numpy.empty_like(xx)
410 yy = numpy.empty_like(xx)
411 for i, source in enumerate(sourceCat):
412 Ixx, Ixy, Iyy = source.getIxx(), source.getIxy(), source.getIyy()
413 if pixToTanPix:
414 p = lsst.geom.Point2D(source.getX(), source.getY())
415 linTransform = afwGeom.linearizeTransform(pixToTanPix, p).getLinear()
416 m = afwGeom.Quadrupole(Ixx, Iyy, Ixy)
417 m.transform(linTransform)
418 Ixx, Iyy, Ixy = m.getIxx(), m.getIyy(), m.getIxy()
419
420 xx[i], xy[i], yy[i] = Ixx, Ixy, Iyy
421
422 width = numpy.sqrt(0.5*(xx + yy))
423 with numpy.errstate(invalid="ignore"): # suppress NAN warnings
424 bad = reduce(lambda x, y: numpy.logical_or(x, sourceCat[y]), self.config.badFlags, False)
425 bad = numpy.logical_or(bad, numpy.logical_not(numpy.isfinite(width)))
426 bad = numpy.logical_or(bad, numpy.logical_not(numpy.isfinite(flux)))
427 if self.config.doFluxLimit:
428 bad = numpy.logical_or(bad, flux < self.config.fluxMin)
429 if self.config.fluxMax > 0:
430 bad = numpy.logical_or(bad, flux > self.config.fluxMax)
431 if self.config.doSignalToNoiseLimit:
432 bad = numpy.logical_or(bad, flux/fluxErr < self.config.signalToNoiseMin)
433 if self.config.signalToNoiseMax > 0:
434 bad = numpy.logical_or(bad, flux/fluxErr > self.config.signalToNoiseMax)
435 bad = numpy.logical_or(bad, width < self.config.widthMin)
436 bad = numpy.logical_or(bad, width > self.config.widthMax)
437 good = numpy.logical_not(bad)
438
439 if not numpy.any(good):
440 raise ObjectSizeNoGoodSourcesError(n_input_sources=len(sourceCat))
441
442 mag = -2.5*numpy.log10(flux[good])
443 width = width[good]
444 #
445 # Look for the maximum in the size histogram, then search upwards for the minimum that separates
446 # the initial peak (of, we presume, stars) from the galaxies
447 #
448 if dumpData:
449 import os
450 import pickle as pickle
451 _ii = 0
452 while True:
453 pickleFile = os.path.expanduser(os.path.join("~", "widths-%d.pkl" % _ii))
454 if not os.path.exists(pickleFile):
455 break
456 _ii += 1
457
458 with open(pickleFile, "wb") as fd:
459 pickle.dump(mag, fd, -1)
460 pickle.dump(width, fd, -1)
461
462 centers, clusterId = _kcenters(width, nCluster=4, useMedian=True,
463 widthStdAllowed=self.config.widthStdAllowed)
464
465 if display and plotMagSize:
466 fig = plot(mag, width, centers, clusterId,
467 magType=self.config.sourceFluxField.split(".")[-1].title(),
468 marker="+", markersize=3, markeredgewidth=None, ltype=':', clear=True)
469 else:
470 fig = None
471
472 clusterId = _improveCluster(width, centers, clusterId,
473 nsigma=self.config.nSigmaClip,
474 widthStdAllowed=self.config.widthStdAllowed)
475
476 if display and plotMagSize:
477 plot(mag, width, centers, clusterId, marker="x", markersize=3, markeredgewidth=None, clear=False)
478
479 stellar = (clusterId == 0)
480 #
481 # We know enough to plot, if so requested
482 #
483 frame = 0
484
485 if fig:
486 if display and displayExposure:
487 disp = afwDisplay.Display(frame=frame)
488 disp.mtv(exposure.getMaskedImage(), title="PSF candidates")
489
490 global eventHandler
491 eventHandler = EventHandler(fig.get_axes()[0], mag, width,
492 sourceCat.getX()[good], sourceCat.getY()[good], frames=[frame])
493
494 fig.show()
495
496 while True:
497 try:
498 reply = input("continue? [c h(elp) q(uit) p(db)] ").strip()
499 except EOFError:
500 reply = None
501 if not reply:
502 reply = "c"
503
504 if reply:
505 if reply[0] == "h":
506 print("""\
507 We cluster the points; red are the stellar candidates and the other colours are other clusters.
508 Points labelled + are rejects from the cluster (only for cluster 0).
509
510 At this prompt, you can continue with almost any key; 'p' enters pdb, and 'h' prints this text
511
512 If displayExposure is true, you can put the cursor on a point and hit 'p' to see it in the
513 image display.
514 """)
515 elif reply[0] == "p":
516 import pdb
517 pdb.set_trace()
518 elif reply[0] == 'q':
519 sys.exit(1)
520 else:
521 break
522
523 if display and displayExposure:
524 mi = exposure.getMaskedImage()
525 with disp.Buffering():
526 for i, source in enumerate(sourceCat):
527 if good[i]:
528 ctype = afwDisplay.GREEN # star candidate
529 else:
530 ctype = afwDisplay.RED # not star
531
532 disp.dot("+", source.getX() - mi.getX0(), source.getY() - mi.getY0(), ctype=ctype)
533
534 # stellar only applies to good==True objects
535 mask = good == True # noqa (numpy bool comparison): E712
536 good[mask] = stellar
537
538 return Struct(selected=good)
An ellipse core with quadrupole moments as parameters.
Definition Quadrupole.h:47
lsst::geom::AffineTransform linearizeTransform(TransformPoint2ToPoint2 const &original, lsst::geom::Point2D const &inPoint)
Approximate a Transform by its local linearization.
_improveCluster(yvec, centers, clusterId, nsigma=2.0, nIteration=10, clusterNum=0, widthStdAllowed=0.15)
_kcenters(yvec, nCluster, useMedian=False, widthStdAllowed=0.15)
plot(mag, width, centers, clusterId, marker="o", markersize=2, markeredgewidth=0, ltype='-', magType="model", clear=True)