Coverage for python / lsst / meas / extensions / trailedSources / NaivePlugin.py: 14%
208 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-22 09:05 +0000
1#
2# This file is part of meas_extensions_trailedSources.
3#
4# Developed for the LSST Data Management System.
5# This product includes software developed by the LSST Project
6# (http://www.lsst.org).
7# See the COPYRIGHT file at the top-level directory of this distribution
8# for details of code ownership.
9#
10# This program is free software: you can redistribute it and/or modify
11# it under the terms of the GNU General Public License as published by
12# the Free Software Foundation, either version 3 of the License, or
13# (at your option) any later version.
14#
15# This program is distributed in the hope that it will be useful,
16# but WITHOUT ANY WARRANTY; without even the implied warranty of
17# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18# GNU General Public License for more details.
19#
20# You should have received a copy of the GNU General Public License
21# along with this program. If not, see <http://www.gnu.org/licenses/>.
22#
24import logging
25import numpy as np
26import scipy.optimize as sciOpt
27from scipy.special import erf
28from math import sqrt
30from lsst.geom import Point2D, Point2I
31from lsst.meas.base.pluginRegistry import register
32from lsst.meas.base import SingleFramePlugin, SingleFramePluginConfig
33from lsst.meas.base import FlagHandler, FlagDefinitionList
34import lsst.pex.config
36from ._trailedSources import VeresModel
37from .utils import getMeasurementCutout
39__all__ = ("SingleFrameNaiveTrailConfig", "SingleFrameNaiveTrailPlugin")
42class SingleFrameNaiveTrailConfig(SingleFramePluginConfig):
43 """Config class for SingleFrameNaiveTrailPlugin.
44 """
45 maxFlux = lsst.pex.config.Field(
46 dtype=float,
47 default=1e10,
48 doc="Maximum calculated model flux before falling back on aperture flux."
49 )
52@register("ext_trailedSources_Naive")
53class SingleFrameNaiveTrailPlugin(SingleFramePlugin):
54 """Naive trailed source measurement plugin
56 Measures the length, angle from +x-axis, and end points of an extended
57 source using the second moments.
59 Parameters
60 ----------
61 config: `SingleFrameNaiveTrailConfig`
62 Plugin configuration.
63 name: `str`
64 Plugin name.
65 schema: `lsst.afw.table.Schema`
66 Schema for the output catalog.
67 metadata: `lsst.daf.base.PropertySet`
68 Metadata to be attached to output catalog.
70 Notes
71 -----
72 This measurement plugin aims to utilize the already measured adaptive
73 second moments to naively estimate the length and angle, and thus
74 end-points, of a fast-moving, trailed source. The length is solved for via
75 finding the root of the difference between the numerical (stack computed)
76 and the analytic adaptive second moments. The angle, theta, from the x-axis
77 is also computed via adaptive moments: theta = arctan(2*Ixy/(Ixx - Iyy))/2.
78 The end points of the trail are then given by (xc +/- (length/2)*cos(theta)
79 and yc +/- (length/2)*sin(theta)), with xc and yc being the centroid
80 coordinates.
82 See also
83 --------
84 lsst.meas.base.SingleFramePlugin
85 """
87 ConfigClass = SingleFrameNaiveTrailConfig
89 @classmethod
90 def getExecutionOrder(cls):
91 # Needs centroids, shape, and flux measurements.
92 # VeresPlugin is run after, which requires image data.
93 return cls.APCORR_ORDER + 0.1
95 def __init__(self, config, name, schema, metadata, logName=None):
96 if logName is None:
97 logName = __name__
98 super().__init__(config, name, schema, metadata, logName=logName)
100 # Measurement Keys
101 self.keyRa = schema.addField(name + "_ra", type="D", doc="Trail centroid right ascension.")
102 self.keyDec = schema.addField(name + "_dec", type="D", doc="Trail centroid declination.")
103 self.keyX0 = schema.addField(name + "_x0", type="D", doc="Trail head X coordinate.", units="pixel")
104 self.keyY0 = schema.addField(name + "_y0", type="D", doc="Trail head Y coordinate.", units="pixel")
105 self.keyX1 = schema.addField(name + "_x1", type="D", doc="Trail tail X coordinate.", units="pixel")
106 self.keyY1 = schema.addField(name + "_y1", type="D", doc="Trail tail Y coordinate.", units="pixel")
107 self.keyFlux = schema.addField(name + "_flux", type="D", doc="Trailed source flux.", units="count")
108 self.keyLength = schema.addField(name + "_length", type="D", doc="Trail length.", units="pixel")
109 self.keyAngle = schema.addField(name + "_angle", type="D", doc="Angle measured from +x-axis.")
111 # Measurement Error Keys
112 self.keyX0Err = schema.addField(name + "_x0Err", type="D",
113 doc="Trail head X coordinate error.", units="pixel")
114 self.keyY0Err = schema.addField(name + "_y0Err", type="D",
115 doc="Trail head Y coordinate error.", units="pixel")
116 self.keyX1Err = schema.addField(name + "_x1Err", type="D",
117 doc="Trail tail X coordinate error.", units="pixel")
118 self.keyY1Err = schema.addField(name + "_y1Err", type="D",
119 doc="Trail tail Y coordinate error.", units="pixel")
120 self.keyFluxErr = schema.addField(name + "_fluxErr", type="D",
121 doc="Trail flux error.", units="count")
122 self.keyLengthErr = schema.addField(name + "_lengthErr", type="D",
123 doc="Trail length error.", units="pixel")
124 self.keyAngleErr = schema.addField(name + "_angleErr", type="D", doc="Trail angle error.")
126 flagDefs = FlagDefinitionList()
127 self.FAILURE = flagDefs.addFailureFlag("No trailed-source measured")
128 self.NO_FLUX = flagDefs.add("flag_noFlux", "No suitable prior flux measurement")
129 self.NO_CONVERGE = flagDefs.add("flag_noConverge", "The root finder did not converge")
130 self.NO_SIGMA = flagDefs.add("flag_noSigma", "No PSF width (sigma)")
131 self.EDGE = flagDefs.add("flag_edge", "Trail contains edge pixels")
132 self.OFFIMAGE = flagDefs.add("flag_off_image", "Trail extends off image")
133 self.NAN = flagDefs.add("flag_nan", "One or more trail coordinates are missing")
134 self.SUSPECT_LONG_TRAIL = flagDefs.add("flag_suspect_long_trail",
135 "Trail length is greater than three times the psf radius")
136 self.SHAPE = flagDefs.add("flag_shape", "Shape flag is set, trail length not calculated")
137 self.flagHandler = FlagHandler.addFields(schema, name, flagDefs)
139 self.log = logging.getLogger(self.logName)
141 def measure(self, measRecord, exposure):
142 """Run the Naive trailed source measurement algorithm.
144 Parameters
145 ----------
146 measRecord : `lsst.afw.table.SourceRecord`
147 Record describing the object being measured.
148 exposure : `lsst.afw.image.Exposure`
149 Pixel data to be measured.
151 See also
152 --------
153 lsst.meas.base.SingleFramePlugin.measure
154 """
155 if measRecord.getShapeFlag():
156 self.log.debug("Shape flag is set for measRecord: %s. Trail measurement "
157 "will not be made. All trail values will be set to nan.", measRecord.getId())
158 self.flagHandler.setValue(measRecord, self.FAILURE.number, True)
159 self.flagHandler.setValue(measRecord, self.SHAPE.number, True)
160 return
162 xc = measRecord["slot_Shape_x"]
163 yc = measRecord["slot_Shape_y"]
164 if not np.isfinite(xc) or not np.isfinite(yc):
165 self.flagHandler.setValue(measRecord, self.SAFE_CENTROID.number, True)
166 self.flagHandler.setValue(measRecord, self.FAILURE.number, True)
167 return
168 ra, dec = self.computeRaDec(exposure, xc, yc)
170 # Transform the second-moments to semi-major and minor axes
171 Ixx, Iyy, Ixy = measRecord.getShape().getParameterVector()
172 xmy = Ixx - Iyy
173 xpy = Ixx + Iyy
174 xmy2 = xmy*xmy
175 xy2 = Ixy*Ixy
176 a2 = 0.5 * (xpy + sqrt(xmy2 + 4.0*xy2))
177 b2 = 0.5 * (xpy - sqrt(xmy2 + 4.0*xy2))
179 # Measure the trail length
180 length, gradLength, results = self.findLength(a2, b2)
181 if not results.converged:
182 self.log.info("Results not converged: %s", results.flag)
183 self.flagHandler.setValue(measRecord, self.NO_CONVERGE.number, True)
184 self.flagHandler.setValue(measRecord, self.FAILURE.number, True)
185 return
187 # Compute the angle of the trail from the x-axis
188 theta = 0.5 * np.arctan2(2.0 * Ixy, xmy)
190 # Get end-points of the trail (there is a degeneracy here)
191 radius = length/2.0 # Trail 'radius'
192 dydtheta = radius*np.cos(theta)
193 dxdtheta = radius*np.sin(theta)
194 x0 = xc - dydtheta
195 y0 = yc - dxdtheta
196 x1 = xc + dydtheta
197 y1 = yc + dxdtheta
199 self.check_trail(measRecord, exposure, x0, y0, x1, y1, length)
201 # Get a cutout of the object from the exposure
202 cutout = getMeasurementCutout(measRecord, exposure)
204 # Compute flux assuming fixed parameters for VeresModel
205 params = np.array([xc, yc, 1.0, length, theta]) # Flux = 1.0
206 model = VeresModel(cutout)
207 flux, gradFlux = model.computeFluxWithGradient(params)
209 # Fall back to aperture flux
210 if (not np.isfinite(flux)) | (np.abs(flux) > self.config.maxFlux):
211 if np.isfinite(measRecord.getApInstFlux()):
212 flux = measRecord.getApInstFlux()
213 else:
214 self.flagHandler.setValue(measRecord, self.NO_FLUX.number, True)
215 self.flagHandler.setValue(measRecord, self.FAILURE.number, True)
216 return
218 # Propogate errors from second moments and centroid
219 IxxErr2, IyyErr2, IxyErr2 = np.diag(measRecord.getShapeErr())
221 # SdssShape does not produce centroid errors. The
222 # Slot centroid errors will suffice for now.
223 xcErr2, ycErr2 = np.diag(measRecord.getCentroidErr())
225 # Error in length
226 desc = sqrt(xmy2 + 4.0*xy2) # Descriminant^1/2 of EV equation
227 da2dIxx = 0.5*(1.0 + (xmy/desc))
228 da2dIyy = 0.5*(1.0 - (xmy/desc))
229 da2dIxy = 2.0*Ixy / desc
230 a2Err2 = IxxErr2*da2dIxx*da2dIxx + IyyErr2*da2dIyy*da2dIyy + IxyErr2*da2dIxy*da2dIxy
231 b2Err2 = IxxErr2*da2dIyy*da2dIyy + IyyErr2*da2dIxx*da2dIxx + IxyErr2*da2dIxy*da2dIxy
232 dLda2, dLdb2 = gradLength
233 lengthErr = np.sqrt(dLda2*dLda2*a2Err2 + dLdb2*dLdb2*b2Err2)
235 # Error in theta
236 dThetadIxx = -Ixy / (xmy2 + 4.0*xy2) # dThetadIxx = -dThetadIyy
237 dThetadIxy = xmy / (xmy2 + 4.0*xy2)
238 thetaErr = sqrt(dThetadIxx*dThetadIxx*(IxxErr2 + IyyErr2) + dThetadIxy*dThetadIxy*IxyErr2)
240 # Error in flux
241 dFdxc, dFdyc, _, dFdL, dFdTheta = gradFlux
242 fluxErr = sqrt(dFdL*dFdL*lengthErr*lengthErr + dFdTheta*dFdTheta*thetaErr*thetaErr
243 + dFdxc*dFdxc*xcErr2 + dFdyc*dFdyc*ycErr2)
245 # Errors in end-points
246 dxdradius = np.cos(theta)
247 dydradius = np.sin(theta)
248 radiusErr2 = lengthErr*lengthErr/4.0
249 xErr2 = sqrt(xcErr2 + radiusErr2*dxdradius*dxdradius + thetaErr*thetaErr*dxdtheta*dxdtheta)
250 yErr2 = sqrt(ycErr2 + radiusErr2*dydradius*dydradius + thetaErr*thetaErr*dydtheta*dydtheta)
251 x0Err = sqrt(xErr2) # Same for x1
252 y0Err = sqrt(yErr2) # Same for y1
254 # Set flags
255 measRecord.set(self.keyRa, ra)
256 measRecord.set(self.keyDec, dec)
257 measRecord.set(self.keyX0, x0)
258 measRecord.set(self.keyY0, y0)
259 measRecord.set(self.keyX1, x1)
260 measRecord.set(self.keyY1, y1)
261 measRecord.set(self.keyFlux, flux)
262 measRecord.set(self.keyLength, length)
263 measRecord.set(self.keyAngle, theta)
264 measRecord.set(self.keyX0Err, x0Err)
265 measRecord.set(self.keyY0Err, y0Err)
266 measRecord.set(self.keyX1Err, x0Err)
267 measRecord.set(self.keyY1Err, y0Err)
268 measRecord.set(self.keyFluxErr, fluxErr)
269 measRecord.set(self.keyLengthErr, lengthErr)
270 measRecord.set(self.keyAngleErr, thetaErr)
272 def check_trail(self, measRecord, exposure, x0, y0, x1, y1, length):
273 """ Set flags for edge pixels, off chip, and nan trail coordinates and
274 flag if trail length is three times larger than psf.
276 Check if the coordinates of the beginning and ending of the trail fall
277 inside the exposures bounding box. If not, set the off_chip flag.
278 If the beginning or ending falls within a pixel marked as edge, set the
279 edge flag. If any of the coordinates happens to fall on a nan, then
280 set the nan flag.
281 Additionally, check if the trail is three times larger than the psf. If
282 so, set the suspect trail flag.
284 Parameters
285 ----------
286 measRecord: `lsst.afw.MeasurementRecord`
287 Record describing the object being measured.
288 exposure: `lsst.afw.Exposure`
289 Pixel data to be measured.
291 x0: `float`
292 x coordinate of the beginning of the trail.
293 y0: `float`
294 y coordinate of the beginning of the trail.
295 x1: `float`
296 x coordinate of the end of the trail.
297 y1: `float`
298 y coordinate of the end of the trail.
299 """
300 x_coords = [x0, x1]
301 y_coords = [y0, y1]
303 # Check if one of the end points of the trail sources is nan. If so,
304 # set the trailed source nan flag.
305 if np.isnan(x_coords).any() or np.isnan(y_coords).any():
306 self.flagHandler.setValue(measRecord, self.NAN.number, True)
307 x_coords = [x for x in x_coords if not np.isnan(x)]
308 y_coords = [y for y in y_coords if not np.isnan(y)]
310 # Check if the non-nan coordinates are within the bounding box
311 if not (all(exposure.getBBox().beginX <= x <= exposure.getBBox().endX for x in x_coords)
312 and all(exposure.getBBox().beginY <= y <= exposure.getBBox().endY for y in y_coords)):
313 self.flagHandler.setValue(measRecord, self.EDGE.number, True)
314 self.flagHandler.setValue(measRecord, self.OFFIMAGE.number, True)
315 else:
316 # Check if edge is set for any of the pixel pairs. Do not
317 # check any that have a nan.
318 for (x_val, y_val) in zip(x_coords, y_coords):
319 if x_val is not np.nan and y_val is not np.nan:
320 if exposure.mask[Point2I(int(x_val),
321 int(y_val))] & exposure.mask.getPlaneBitMask('EDGE') != 0:
322 self.flagHandler.setValue(measRecord, self.EDGE.number, True)
323 # Check whether trail extends off the edge of the exposure. Allows nans
324 # as their location
325 elif not (all(exposure.getBBox().beginX <= x <= exposure.getBBox().endX for x in x_coords)
326 and all(exposure.getBBox().beginY <= y <= exposure.getBBox().endY for y in y_coords)):
327 self.flagHandler.setValue(measRecord, self.EDGE.number, True)
328 self.flagHandler.setValue(measRecord, self.OFFIMAGE.number, True)
329 else:
330 # Check whether the beginning or end point of the trail has the
331 # edge flag set. The end points are not whole pixel values, so
332 # the pixel value must be rounded.
333 if exposure.mask[Point2I(int(x0), int(y0))] and exposure.mask[Point2I(int(x1), int(y1))]:
334 if ((exposure.mask[Point2I(int(x0), int(y0))] & exposure.mask.getPlaneBitMask('EDGE') != 0)
335 or (exposure.mask[Point2I(int(x1), int(y1))]
336 & exposure.mask.getPlaneBitMask('EDGE') != 0)):
337 self.flagHandler.setValue(measRecord, self.EDGE.number, True)
339 psfShape = exposure.psf.computeShape(exposure.getBBox().getCenter())
340 psfRadius = psfShape.getDeterminantRadius()
342 if length > psfRadius*3.0:
343 self.flagHandler.setValue(measRecord, self.SUSPECT_LONG_TRAIL.number, True)
345 def fail(self, measRecord, error=None):
346 """Record failure
348 See also
349 --------
350 lsst.meas.base.SingleFramePlugin.fail
351 """
352 if error is None:
353 self.flagHandler.handleFailure(measRecord)
354 else:
355 self.flagHandler.handleFailure(measRecord, error.cpp)
357 @staticmethod
358 def _computeSecondMomentDiff(z, c):
359 """Compute difference of the numerical and analytic second moments.
361 Parameters
362 ----------
363 z : `float`
364 Proportional to the length of the trail. (see notes)
365 c : `float`
366 Constant (see notes)
368 Returns
369 -------
370 diff : `float`
371 Difference in numerical and analytic second moments.
373 Notes
374 -----
375 This is a simplified expression for the difference between the stack
376 computed adaptive second-moment and the analytic solution. The variable
377 z is proportional to the length such that length=2*z*sqrt(2*(Ixx+Iyy)),
378 and c is a constant (c = 4*Ixx/((Ixx+Iyy)*sqrt(pi))). Both have been
379 defined to avoid unnecessary floating-point operations in the root
380 finder.
381 """
383 diff = erf(z) - c*z*np.exp(-z*z)
384 return diff
386 @classmethod
387 def findLength(cls, Ixx, Iyy):
388 """Find the length of a trail, given adaptive second-moments.
390 Uses a root finder to compute the length of a trail corresponding to
391 the adaptive second-moments computed by previous measurements
392 (ie. SdssShape).
394 Parameters
395 ----------
396 Ixx : `float`
397 Adaptive second-moment along x-axis.
398 Iyy : `float`
399 Adaptive second-moment along y-axis.
401 Returns
402 -------
403 length : `float`
404 Length of the trail.
405 results : `scipy.optimize.RootResults`
406 Contains messages about convergence from the root finder.
407 """
409 xpy = Ixx + Iyy
410 c = 4.0*Ixx/(xpy*np.sqrt(np.pi))
412 # Given a 'c' in (c_min, c_max], the root is contained in (0,1].
413 # c_min is given by the case: Ixx == Iyy, ie. a point source.
414 # c_max is given by the limit Ixx >> Iyy.
415 # Empirically, 0.001 is a suitable lower bound, assuming Ixx > Iyy.
416 z, results = sciOpt.brentq(lambda z: cls._computeSecondMomentDiff(z, c),
417 0.001, 1.0, full_output=True)
419 length = 2.0*z*np.sqrt(2.0*xpy)
420 gradLength = cls._gradFindLength(Ixx, Iyy, z, c)
421 return length, gradLength, results
423 @staticmethod
424 def _gradFindLength(Ixx, Iyy, z, c):
425 """Compute the gradient of the findLength function.
426 """
427 spi = np.sqrt(np.pi)
428 xpy = Ixx+Iyy
429 xpy2 = xpy*xpy
430 enz2 = np.exp(-z*z)
431 sxpy = np.sqrt(xpy)
433 fac = 4.0 / (spi*xpy2)
434 dcdIxx = Iyy*fac
435 dcdIyy = -Ixx*fac
437 # Derivatives of the _computeMomentsDiff function
438 dfdc = z*enz2
439 dzdf = spi / (enz2*(spi*c*(2.0*z*z - 1.0) + 2.0)) # inverse of dfdz
441 dLdz = 2.0*np.sqrt(2.0)*sxpy
442 pLpIxx = np.sqrt(2.0)*z / sxpy # Same as pLpIyy
444 dLdc = dLdz*dzdf*dfdc
445 dLdIxx = dLdc*dcdIxx + pLpIxx
446 dLdIyy = dLdc*dcdIyy + pLpIxx
447 return dLdIxx, dLdIyy
449 @staticmethod
450 def computeLength(Ixx, Iyy):
451 """Compute the length of a trail, given unweighted second-moments.
452 """
453 denom = np.sqrt(Ixx - 2.0*Iyy)
455 length = np.sqrt(6.0)*denom
457 dLdIxx = np.sqrt(1.5) / denom
458 dLdIyy = -np.sqrt(6.0) / denom
459 return length, (dLdIxx, dLdIyy)
461 @staticmethod
462 def computeRaDec(exposure, x, y):
463 """Convert pixel coordinates to RA and Dec.
465 Parameters
466 ----------
467 exposure : `lsst.afw.image.ExposureF`
468 Exposure object containing the WCS.
469 x : `float`
470 x coordinate of the trail centroid
471 y : `float`
472 y coodinate of the trail centroid
474 Returns
475 -------
476 ra : `float`
477 Right ascension.
478 dec : `float`
479 Declination.
480 """
482 wcs = exposure.getWcs()
483 center = wcs.pixelToSky(Point2D(x, y))
484 ra = center.getRa().asDegrees()
485 dec = center.getDec().asDegrees()
486 return ra, dec