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