Coverage for python / lsst / meas / algorithms / objectSizeStarSelector.py: 14%
268 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-22 09:03 +0000
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/>.
22__all__ = ["ObjectSizeStarSelectorConfig", "ObjectSizeStarSelectorTask",
23 "ObjectSizeNoSourcesError", "ObjectSizeNoGoodSourcesError"]
25import sys
27import numpy
28import warnings
29from functools import reduce
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
40afwDisplay.setDefaultMaskTransparency(75)
42_LOG = getLogger(__name__)
45class ObjectSizeNoSourcesError(AlgorithmError):
46 """Raised if no sources were supplied in the input catalog."""
48 def __init__(self) -> None:
49 super().__init__("Input catalog for source selection is empty.")
51 @property
52 def metadata(self) -> dict:
53 return {}
56class ObjectSizeNoGoodSourcesError(AlgorithmError):
57 """Raised if no good sources were available after applying configured
58 selection criteria.
59 """
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.")
65 @property
66 def metadata(self) -> dict:
67 return {
68 "n_input_sources": self._n_available,
69 }
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 )
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)
163class EventHandler:
164 """A class to handle key strokes with matplotlib displays.
165 """
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
175 self.cid = self.axes.figure.canvas.mpl_connect('key_press_event', self)
177 def __call__(self, ev):
178 if ev.inaxes != self.axes:
179 return
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
185 which = numpy.where(dist == min(dist))
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
197def _assignClusters(yvec, centers):
198 """Return a vector of centerIds based on their distance to the centers.
199 """
200 assert len(centers) > 0
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)
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]))
221 minDist[update] = dist[update]
222 clusterId[update] = i
223 numpy.seterr(**oldSettings)
225 return clusterId
228def _kcenters(yvec, nCluster, useMedian=False, widthStdAllowed=0.15):
229 """A classic k-means algorithm, clustering yvec into nCluster clusters
231 Return the set of centres, and the cluster ID for each of the points
233 If useMedian is true, use the median of the cluster as its centre, rather than
234 the traditional mean
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 """
242 assert nCluster > 0
244 mean0 = sorted(yvec)[len(yvec)//10] # guess
245 delta = mean0 * widthStdAllowed * 2.0
246 centers = mean0 + delta * numpy.arange(nCluster)
248 func = numpy.median if useMedian else numpy.mean
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)
256 if numpy.all(clusterId == oclusterId):
257 break
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
267 return centers, clusterId
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 """
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
280 inCluster0 = clusterId == clusterNum
281 yv = yvec[inCluster0]
283 centers[clusterNum] = numpy.median(yv)
284 stdev = numpy.std(yv)
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)]
290 sd = stdev if stdev < stdev_iqr else stdev_iqr
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
298 nMember = sum(clusterId == clusterNum)
299 # 'sd < widthStdAllowed * median' prevents too much rejections
300 if nMember == old_nMember or sd < widthStdAllowed * median:
301 break
303 return clusterId
306def plot(mag, width, centers, clusterId, marker="o", markersize=2, markeredgewidth=0, ltype='-',
307 magType="model", clear=True):
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
316 try:
317 fig
318 except NameError:
319 fig = plt.figure()
320 else:
321 if clear:
322 fig.clf()
324 axes = fig.add_axes((0.1, 0.1, 0.85, 0.80))
326 xmin = sorted(mag)[int(0.05*len(mag))]
327 xmax = sorted(mag)[int(0.95*len(mag))]
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)
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)
338 li = (clusterId == k)
339 axes.plot(mag[li], width[li], marker, markersize=markersize, markeredgewidth=markeredgewidth,
340 color=colors[k % len(colors)])
342 li = (clusterId == -1)
343 axes.plot(mag[li], width[li], marker, markersize=markersize, markeredgewidth=markeredgewidth,
344 color='k')
346 if clear:
347 axes.set_xlabel("Instrumental %s mag" % magType)
348 axes.set_ylabel(r"$\sqrt{(I_{xx} + I_{yy})/2}$")
350 return fig
353@pexConfig.registerConfigurable("objectSize", sourceSelectorRegistry)
354class ObjectSizeStarSelectorTask(BaseSourceSelectorTask):
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 """
359 ConfigClass = ObjectSizeStarSelectorConfig
360 usesMatches = False # selectStars does not use its matches argument
362 def selectSources(self, sourceCat, matches=None, exposure=None):
363 """Return a selection of PSF candidates that represent likely stars.
365 A list of PSF candidates may be used by a PSF fitter to construct a PSF.
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.
378 Returns
379 -------
380 struct : `lsst.pipe.base.Struct`
381 The struct contains the following data:
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
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?
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"]
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()
420 xx[i], xy[i], yy[i] = Ixx, Ixy, Iyy
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)
439 if not numpy.any(good):
440 raise ObjectSizeNoGoodSourcesError(n_input_sources=len(sourceCat))
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
458 with open(pickleFile, "wb") as fd:
459 pickle.dump(mag, fd, -1)
460 pickle.dump(width, fd, -1)
462 centers, clusterId = _kcenters(width, nCluster=4, useMedian=True,
463 widthStdAllowed=self.config.widthStdAllowed)
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
472 clusterId = _improveCluster(width, centers, clusterId,
473 nsigma=self.config.nSigmaClip,
474 widthStdAllowed=self.config.widthStdAllowed)
476 if display and plotMagSize:
477 plot(mag, width, centers, clusterId, marker="x", markersize=3, markeredgewidth=None, clear=False)
479 stellar = (clusterId == 0)
480 #
481 # We know enough to plot, if so requested
482 #
483 frame = 0
485 if fig:
486 if display and displayExposure:
487 disp = afwDisplay.Display(frame=frame)
488 disp.mtv(exposure.getMaskedImage(), title="PSF candidates")
490 global eventHandler
491 eventHandler = EventHandler(fig.get_axes()[0], mag, width,
492 sourceCat.getX()[good], sourceCat.getY()[good], frames=[frame])
494 fig.show()
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"
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).
510 At this prompt, you can continue with almost any key; 'p' enters pdb, and 'h' prints this text
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
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
532 disp.dot("+", source.getX() - mi.getX0(), source.getY() - mi.getY0(), ctype=ctype)
534 # stellar only applies to good==True objects
535 mask = good == True # noqa (numpy bool comparison): E712
536 good[mask] = stellar
538 return Struct(selected=good)