Coverage for python / lsst / meas / extensions / psfex / psfexStarSelector.py: 0%
198 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-30 09:04 +0000
1# This file is part of meas_extensions_psfex.
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__ = ("PsfexStarSelectorConfig", "PsfexStarSelectorTask")
24import re
25import sys
27import numpy as np
28try:
29 import matplotlib.pyplot as plt
30 fig = None
31except ImportError:
32 plt = None
34from lsst.pipe.base import Struct
35import lsst.pex.config as pexConfig
36import lsst.afw.display as afwDisplay
37from lsst.meas.algorithms import BaseSourceSelectorTask, sourceSelectorRegistry
38from . import psfexLib
39from .psfex import compute_fwhmrange
42class PsfexStarSelectorConfig(BaseSourceSelectorTask.ConfigClass):
43 fluxName = pexConfig.Field(
44 dtype=str,
45 doc="Name of photometric flux key ",
46 default="base_PsfFlux",
47 )
48 fluxErrName = pexConfig.Field(
49 dtype=str,
50 doc="Name of phot. flux err. key",
51 default="",
52 )
53 minFwhm = pexConfig.Field(
54 dtype=float,
55 doc="Maximum allowed FWHM ",
56 default=2,
57 )
58 maxFwhm = pexConfig.Field(
59 dtype=float,
60 doc="Minimum allowed FWHM ",
61 default=10,
62 )
63 maxFwhmVariability = pexConfig.Field(
64 dtype=float,
65 doc="Allowed FWHM variability (1.0 = 100%)",
66 default=0.2,
67 )
68 maxellip = pexConfig.Field(
69 dtype=float,
70 doc="Maximum (A-B)/(A+B) ",
71 default=0.3,
72 check=lambda x: x >= 0.0,
73 )
74 minsn = pexConfig.Field(
75 dtype=float,
76 doc="Minimum S/N for candidates",
77 default=100,
78 check=lambda x: x >= 0.0,
79 )
81 def validate(self):
82 pexConfig.Config.validate(self)
84 if self.fluxErrName == "":
85 self.fluxErrName = self.fluxName + ".err"
86 elif self.fluxErrName != self.fluxName + ".err":
87 msg = f"fluxErrName ({self.fluxErrName}) doesn't correspond to fluxName ({self.fluxName})"
88 raise pexConfig.FieldValidationError(PsfexStarSelectorConfig.fluxName, self, msg)
90 if self.minFwhm > self.maxFwhm:
91 raise pexConfig.FieldValidationError(PsfexStarSelectorConfig.minFwhm, self,
92 f"minFwhm ({self.minFwhm}) > maxFwhm ({self.maxFwhm})")
94 def setDefaults(self):
95 super().setDefaults()
96 self.badFlags = [
97 "base_PixelFlags_flag_edge",
98 "base_PixelFlags_flag_nodata",
99 "base_PixelFlags_flag_saturatedCenter",
100 "base_PixelFlags_flag_crCenter",
101 "base_PixelFlags_flag_bad",
102 "base_PixelFlags_flag_suspectCenter",
103 "base_PsfFlux_flag",
104 # "parent", # actually this is a test on deblend_nChild
105 ]
108class EventHandler():
109 """A class to handle key strokes with matplotlib displays
110 """
112 def __init__(self, axes, xs, ys, x, y, frames=[0]):
113 self.axes = axes
114 self.xs = xs
115 self.ys = ys
116 self.x = x
117 self.y = y
118 self.frames = frames
120 self.cid = self.axes.figure.canvas.mpl_connect('key_press_event', self)
122 def __call__(self, ev):
123 if ev.inaxes != self.axes:
124 return
126 if ev.key and ev.key in ("p"):
127 dist = np.hypot(self.xs - ev.xdata, self.ys - ev.ydata)
128 dist[np.where(np.isnan(dist))] = 1e30
130 which = np.where(dist == min(dist))
132 x = self.x[which][0]
133 y = self.y[which][0]
134 for frame in self.frames:
135 disp = afwDisplay.Display(frame=frame)
136 disp.pan(x, y)
137 disp.flush()
138 else:
139 pass
142def plot(mag, width, centers, clusterId, marker="o", markersize=2, markeredgewidth=0, ltype='-',
143 clear=True):
145 global fig
146 if not fig:
147 fig = plt.figure()
148 newFig = True
149 else:
150 newFig = False
151 if clear:
152 fig.clf()
154 axes = fig.add_axes((0.1, 0.1, 0.85, 0.80))
156 xmin = sorted(mag)[int(0.05*len(mag))]
157 xmax = sorted(mag)[int(0.95*len(mag))]
159 axes.set_xlim(-17.5, -13)
160 axes.set_xlim(xmin - 0.1*(xmax - xmin), xmax + 0.1*(xmax - xmin))
161 axes.set_ylim(0, 10)
163 colors = ["r", "g", "b", "c", "m", "k", ]
164 for k, mean in enumerate(centers):
165 if k == 0:
166 axes.plot(axes.get_xlim(), (mean, mean,), "k%s" % ltype)
168 ll = (clusterId == k)
169 axes.plot(mag[ll], width[ll], marker, markersize=markersize, markeredgewidth=markeredgewidth,
170 color=colors[k%len(colors)])
172 ll = (clusterId == -1)
173 axes.plot(mag[ll], width[ll], marker, markersize=markersize, markeredgewidth=markeredgewidth,
174 color='k')
176 if newFig:
177 axes.set_xlabel("model")
178 axes.set_ylabel(r"$\sqrt{I_{xx} + I_{yy}}$")
180 return fig
182## \addtogroup LSST_task_documentation
183## \{
184## \page PsfexStarSelectorTask
185## \ref PsfexStarSelectorTask_ "PsfexStarSelectorTask"
186## \copybrief PsfexStarSelectorTask
187## \}
190@pexConfig.registerConfigurable("psfex", sourceSelectorRegistry)
191class PsfexStarSelectorTask(BaseSourceSelectorTask):
192 r"""A star selector whose algorithm is not yet documented.
194 @anchor PsfexStarSelectorTask_
196 @section meas_extensions_psfex_psfexStarSelectorStarSelector_Contents Contents
198 - @ref meas_extensions_psfex_psfexStarSelectorStarSelector_Purpose
199 - @ref meas_extensions_psfex_psfexStarSelectorStarSelector_Initialize
200 - @ref meas_extensions_psfex_psfexStarSelectorStarSelector_IO
201 - @ref meas_extensions_psfex_psfexStarSelectorStarSelector_Config
202 - @ref meas_extensions_psfex_psfexStarSelectorStarSelector_Debug
204 @section meas_extensions_psfex_psfexStarSelectorStarSelector_Purpose Description
206 A star selector whose algorithm is not yet documented
208 @section meas_extensions_psfex_psfexStarSelectorStarSelector_Initialize Task initialisation
210 @copydoc \_\_init\_\_
212 @section meas_extensions_psfex_psfexStarSelectorStarSelector_IO Invoking the Task
214 Like all star selectors, the main method is `run`.
216 @section meas_extensions_psfex_psfexStarSelectorStarSelector_Config Configuration parameters
218 See @ref PsfexStarSelectorConfig
220 @section meas_extensions_psfex_psfexStarSelectorStarSelector_Debug Debug variables
222 PsfexStarSelectorTask has a debug dictionary with the following keys:
223 <dl>
224 <dt>display
225 <dd>bool; if True display debug information
226 <dt>displayExposure
227 <dd>bool; if True display the exposure and spatial cells
228 <dt>plotFwhmHistogram
229 <dd>bool; if True plot histogram of FWHM
230 <dt>plotFlags
231 <dd>bool: if True plot the sources coloured by their flags
232 <dt>plotRejection
233 <dd>bool; if True plot why sources are rejected
234 </dl>
236 For example, put something like:
237 @code{.py}
238 import lsstDebug
239 def DebugInfo(name):
240 di = lsstDebug.getInfo(name) # N.b. lsstDebug.Info(name) would call us recursively
241 if name.endswith("objectSizeStarSelector"):
242 di.display = True
243 di.displayExposure = True
244 di.plotFwhmHistogram = True
246 return di
248 lsstDebug.Info = DebugInfo
249 @endcode
250 into your `debug.py` file and run your task with the `--debug` flag.
251 """ # noqa: W505
252 ConfigClass = PsfexStarSelectorConfig
253 usesMatches = False # selectStars does not use its matches argument
255 def selectSources(self, sourceCat, matches=None, exposure=None):
256 """Return a selection of psf-like objects.
258 Parameters
259 ----------
260 sourceCat : `lsst.afw.table.SourceCatalog`
261 Catalog of sources to select from.
262 This catalog must be contiguous in memory.
263 matches : `list` of `lsst.afw.table.ReferenceMatch` or None
264 Ignored by this source selector.
265 exposure : `lsst.afw.image.Exposure` or None
266 The exposure the catalog was built from; used for debug display.
268 Return
269 ------
270 struct : `lsst.pipe.base.Struct`
271 The struct contains the following data:
273 - selected : `numpy.ndarray` of `bool``
274 Boolean array of sources that were selected, same length as
275 sourceCat.
276 """
277 import lsstDebug
278 display = lsstDebug.Info(__name__).display
280 displayExposure = display and \
281 lsstDebug.Info(__name__).displayExposure # display the Exposure + spatialCells
282 plotFwhmHistogram = display and plt and \
283 lsstDebug.Info(__name__).plotFwhmHistogram # Plot histogram of FWHM
284 plotFlags = display and plt and \
285 lsstDebug.Info(__name__).plotFlags # Plot the sources coloured by their flags
286 plotRejection = display and plt and \
287 lsstDebug.Info(__name__).plotRejection # Plot why sources are rejected
288 afwDisplay.setDefaultMaskTransparency(75)
290 fluxName = self.config.fluxName
291 fluxErrName = self.config.fluxErrName
292 minFwhm = self.config.minFwhm
293 maxFwhm = self.config.maxFwhm
294 maxFwhmVariability = self.config.maxFwhmVariability
295 maxellip = self.config.maxellip
296 minsn = self.config.minsn
298 maxelong = (maxellip + 1.0)/(1.0 - maxellip) if maxellip < 1.0 else 100
300 # Unpack the catalogue
301 shape = sourceCat.getShapeDefinition()
302 ixx = sourceCat.get("%s.xx" % shape)
303 iyy = sourceCat.get("%s.yy" % shape)
305 fwhm = 2*np.sqrt(2*np.log(2))*np.sqrt(0.5*(ixx + iyy))
306 elong = 0.5*(ixx - iyy)/(ixx + iyy)
308 flux = sourceCat.get(fluxName)
309 fluxErr = sourceCat.get(fluxErrName)
310 sn = flux/np.where(fluxErr > 0, fluxErr, 1)
311 sn[fluxErr <= 0] = -psfexLib.BIG
313 flags = 0x0
314 for i, f in enumerate(self.config.badFlags):
315 flags = np.bitwise_or(flags, np.where(sourceCat.get(f), 1 << i, 0))
316 #
317 # Estimate the acceptable range of source widths
318 #
319 good = np.logical_and(sn > minsn, np.logical_not(flags))
320 good = np.logical_and(good, elong < maxelong)
321 good = np.logical_and(good, fwhm >= minFwhm)
322 good = np.logical_and(good, fwhm < maxFwhm)
324 fwhmMode, fwhmMin, fwhmMax = compute_fwhmrange(fwhm[good], maxFwhmVariability, minFwhm, maxFwhm,
325 plot=dict(fwhmHistogram=plotFwhmHistogram))
327 # -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
328 #
329 # Here's select_candidates
330 #
331 # ---- Apply some selection over flags, fluxes...
333 bad = (flags != 0)
334 # set.setBadFlags(int(sum(bad)))
336 if plotRejection:
337 selectionVectors = []
338 selectionVectors.append((bad, "flags %d" % sum(bad)))
340 dbad = sn < minsn
341 # set.setBadSN(int(sum(dbad)))
342 bad = np.logical_or(bad, dbad)
343 if plotRejection:
344 selectionVectors.append((dbad, "S/N %d" % sum(dbad)))
346 dbad = fwhm < fwhmMin
347 # set.setBadFrmin(int(sum(dbad)))
348 bad = np.logical_or(bad, dbad)
349 if plotRejection:
350 selectionVectors.append((dbad, "fwhmMin %d" % sum(dbad)))
352 dbad = fwhm > fwhmMax
353 # set.setBadFrmax(int(sum(dbad)))
354 bad = np.logical_or(bad, dbad)
355 if plotRejection:
356 selectionVectors.append((dbad, "fwhmMax %d" % sum(dbad)))
358 dbad = elong > maxelong
359 # set.setBadElong(int(sum(dbad)))
360 bad = np.logical_or(bad, dbad)
361 if plotRejection:
362 selectionVectors.append((dbad, "elong %d" % sum(dbad)))
364 good = np.logical_not(bad)
365 #
366 # We know enough to plot, if so requested
367 #
368 frame = 0
369 if displayExposure:
370 mi = exposure.getMaskedImage()
371 disp = afwDisplay.Display(frame=frame)
372 disp.mtv(mi, title="PSF candidates")
374 with disp.Buffering():
375 for i, source in enumerate(sourceCat):
376 if good[i]:
377 ctype = afwDisplay.GREEN # star candidate
378 else:
379 ctype = afwDisplay.RED # not star
381 disp.dot("+", source.getX() - mi.getX0(), source.getY() - mi.getY0(),
382 ctype=ctype)
384 if plotFlags or plotRejection:
385 imag = -2.5*np.log10(flux)
386 plt.clf()
388 alpha = 0.5
389 if plotFlags:
390 isSet = np.where(flags == 0x0)[0]
391 plt.plot(imag[isSet], fwhm[isSet], 'o', alpha=alpha, label="good")
393 for i, f in enumerate(self.config.badFlags):
394 mask = 1 << i
395 isSet = np.where(np.bitwise_and(flags, mask))[0]
396 if isSet.any():
397 if np.isfinite(imag[isSet] + fwhm[isSet]).any():
398 label = re.sub(r"\_flag", "",
399 re.sub(r"^base\_", "",
400 re.sub(r"^.*base\_PixelFlags\_flag\_", "", f)))
401 plt.plot(imag[isSet], fwhm[isSet], 'o', alpha=alpha, label=label)
402 else:
403 for bad, label in selectionVectors:
404 plt.plot(imag[bad], fwhm[bad], 'o', alpha=alpha, label=label)
406 plt.plot(imag[good], fwhm[good], 'o', color="black", label="selected")
407 [plt.axhline(_, color='red') for _ in [fwhmMin, fwhmMax]]
408 plt.xlim(np.median(imag[good]) + 5*np.array([-1, 1]))
409 plt.ylim(fwhm[np.where(np.isfinite(fwhm + imag))].min(), 2*fwhmMax)
410 plt.legend(loc=2)
411 plt.xlabel("Instrumental %s Magnitude" % fluxName.split(".")[-1].title())
412 plt.ylabel("fwhm")
413 title = "PSFEX Star Selection"
414 plt.title("%s %d selected" % (title, sum(good)))
416 if displayExposure:
417 global eventHandler
418 eventHandler = EventHandler(plt.axes(), imag, fwhm, sourceCat.getX(), sourceCat.getY(),
419 frames=[frame])
421 if plotFlags or plotRejection:
422 while True:
423 try:
424 reply = input("continue? [y[es] h(elp) p(db) q(uit)] ").strip()
425 except EOFError:
426 reply = "y"
428 if not reply:
429 reply = "y"
431 if reply[0] == "h":
432 print("""\
433At this prompt, you can continue with almost any key; 'p' enters pdb,
434 'q' returns to the shell, and
435 'h' prints this text
436""", end=' ')
438 if displayExposure:
439 print("""
440If you put the cursor on a point in the matplotlib scatter plot and hit 'p' you'll see it in ds9.""")
441 elif reply[0] == "p":
442 import pdb
443 pdb.set_trace()
444 elif reply[0] == 'q':
445 sys.exit(1)
446 else:
447 break
449 return Struct(selected=good)