Coverage for python/lsst/meas/algorithms/findGlintTrails.py: 92%
158 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-30 02:00 -0700
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__ = ["FindGlintTrailsConfig", "FindGlintTrailsTask", "FittedLine", "GlintTrailParameters", "fit_line"]
24import collections
25import dataclasses
26import math
28import numpy as np
29import scipy.spatial
31import lsst.afw.table
32import lsst.pex.config
33import lsst.pipe.base
36class FittedLine:
37 """A straight line through ``center`` along the unit vector
38 ``direction``.
40 Parameters
41 ----------
42 center : `numpy.ndarray`
43 ``(x, y)`` of a point on the line: the mean of the points it was fit
44 to.
45 direction : `numpy.ndarray`
46 Unit vector along the line; its sign is arbitrary.
47 """
49 def __init__(self, center, direction):
50 self.center = center
51 self.direction = direction
53 def distance(self, x, y):
54 """Compute the perpendicular distance of points from the line.
56 Parameters
57 ----------
58 x, y : `numpy.ndarray` or `float`
59 Coordinates of the points.
61 Returns
62 -------
63 distance : `numpy.ndarray` or `float`
64 Distance of each point from the line, in the units of ``x`` and
65 ``y``.
66 """
67 return np.abs((x - self.center[0])*self.direction[1] - (y - self.center[1])*self.direction[0])
70def fit_line(x, y):
71 """Fit the line that minimizes the summed squared perpendicular distances
72 of points from it.
74 Unlike a regression of ``y`` on ``x``, this fits lines at any angle,
75 including vertical ones, equally well.
77 Parameters
78 ----------
79 x, y : `numpy.ndarray`
80 Coordinates of at least two points.
82 Returns
83 -------
84 line : `FittedLine`
85 The fitted line.
86 """
87 points = np.column_stack((x, y))
88 center = points.mean(axis=0)
89 direction = np.linalg.svd(points - center)[2][0]
90 return FittedLine(center=center, direction=direction)
93def _ransac_line(x, y, threshold, min_samples, seed, max_trials=100):
94 """Fit a line to points that may include outliers, using RANSAC.
96 Each trial fits a line to ``min_samples`` randomly chosen points. The
97 trial with the most points within ``threshold`` of its line wins, with
98 ties going to the smaller summed squared distance. The line is iteratively
99 refit with new selections of inliers until the results stop changing, so
100 that every inlier is within ``threshold`` pixels of the returned line.
102 Parameters
103 ----------
104 x, y : `numpy.ndarray`
105 Coordinates of the points.
106 threshold : `float`
107 Maximum perpendicular distance of an inlier from the line.
108 min_samples : `int`
109 Number of points each trial line is fit to.
110 seed : `int`
111 Seed for the random choice of points, so that fits are repeatable.
112 max_trials : `int`, optional
113 Maximum number of trial lines.
115 Returns
116 -------
117 line : `FittedLine` or `None`
118 The line fit to the inliers, or `None` if fewer than ``min_samples``
119 points are inliers.
120 inliers : `numpy.ndarray` [`bool`] or `None`
121 Which of the points lie within ``threshold`` of ``line``.
122 """
123 if len(x) < min_samples: 123 ↛ 124line 123 didn't jump to line 124 because the condition on line 123 was never true
124 return None, None
125 rng = np.random.default_rng(seed)
126 best = None
127 best_count = min_samples - 1
128 best_sum = np.inf
129 for _ in range(max_trials):
130 sample = rng.choice(len(x), size=min_samples, replace=False)
131 distance = fit_line(x[sample], y[sample]).distance(x, y)
132 inliers = distance <= threshold
133 count = np.sum(inliers)
134 total = np.sum(distance[inliers]**2)
135 if count > best_count or (count == best_count and total < best_sum):
136 best, best_count, best_sum = inliers, count, total
137 if count == len(x):
138 break
139 if best is None: 139 ↛ 140line 139 didn't jump to line 140 because the condition on line 139 was never true
140 return None, None
141 inliers = best
142 for _ in range(max_trials): 142 ↛ 150line 142 didn't jump to line 150 because the loop on line 142 didn't complete
143 line = fit_line(x[inliers], y[inliers])
144 refit = line.distance(x, y) <= threshold
145 if np.sum(refit) < min_samples: 145 ↛ 146line 145 didn't jump to line 146 because the condition on line 145 was never true
146 return None, None
147 if np.array_equal(refit, inliers):
148 break
149 inliers = refit
150 return line, refit
153class FindGlintTrailsConfig(lsst.pex.config.Config):
154 radius = lsst.pex.config.Field(
155 doc="Radius to search for glint trail candidates from each source (pixels).",
156 dtype=float,
157 default=500,
158 )
159 min_points = lsst.pex.config.Field(
160 doc="Minimum number of points to be considered a possible glint trail.",
161 dtype=int,
162 default=5,
163 check=lambda x: x >= 3,
164 )
165 threshold = lsst.pex.config.Field(
166 doc="Maximum perpendicular distance of a source from the line fit to its trail (pixels)."
167 " Real glints are highly collinear, with a typical deviation from a line of < 1 pixel.",
168 dtype=float,
169 default=1.5,
170 )
171 seed = lsst.pex.config.Field(
172 doc="Random seed for RANSAC fitter, to ensure stable fitting.",
173 dtype=int,
174 default=42,
175 )
176 bad_flags = lsst.pex.config.ListField[str](
177 doc="Do not fit sources that have these flags set.",
178 default=["ip_diffim_DipoleFit_classification",
179 "is_negative",
180 ],
181 )
184@dataclasses.dataclass(frozen=True, kw_only=True)
185class GlintTrailParameters:
186 """Holds values from the line fit to a single glint trail.
188 ``slope`` and ``intercept`` are very large or non-finite for a
189 near-vertical trail; ``angle`` with ``center_x`` and ``center_y`` describe
190 the line at any orientation.
191 """
192 slope: float
193 intercept: float
194 stderr: float # root mean squared perpendicular distance, pixels
195 length: float # pixels
196 angle: float # radians from the +X axis, in (-pi/2, pi/2]
197 center_x: float # pixels; the line passes through (center_x, center_y)
198 center_y: float # pixels
201class FindGlintTrailsTask(lsst.pipe.base.Task):
202 """Find glint trails in a catalog by searching for sources that lie in a
203 line.
205 Notes
206 -----
207 For each source ("anchor") in the input catalog that was not included in
208 an an earlier iteration as part of a trail:
209 * Find all sources within a given radius.
210 * For each pair of anchor and match, identify the other sources that
211 could lie on the same line(s).
212 * Take the longest set of such pairs as a candidate trail.
213 * Fit a line to the identified pairs with the RANSAC algorithm,
214 measuring distances perpendicular to the line so that trails at any
215 angle are found.
216 * Find all sources in the catalog that could lie on that line.
217 * Refit a line to all of the matched sources.
218 * If the line holds at least the minimum number of sources, return the
219 sources that were considered inliers during the fit, and the fit
220 parameters.
221 """
223 ConfigClass = FindGlintTrailsConfig
224 _DefaultName = "findGlintTrails"
226 def run(self, catalog):
227 """Find glint trails in a catalog.
229 Parameters
230 ----------
231 catalog : `lsst.afw.table.SourceCatalog`
232 Catalog to search for glint trails.
234 Returns
235 -------
236 result : `lsst.pipe.base.Struct`
237 Results as a struct with attributes:
239 ``trails``
240 Catalog subsets containing sources in each trail that was found.
241 (`list` [`lsst.afw.table.SourceCatalog`])
242 ``trailed_ids``
243 Ids of all the sources that were included in any fit trail.
244 (`set` [`int`])
245 ``parameters``
246 Parameters of all the trails that were found.
247 (`list` [`GlintTrailParameters`])
248 """
249 good_catalog = self._select_good_sources(catalog)
251 matches = lsst.afw.table.matchXy(good_catalog, self.config.radius)
252 per_id = collections.defaultdict(list)
253 for match in matches:
254 per_id[match.first["id"]].append(match)
255 counts = {id: len(value) for id, value in per_id.items()}
257 trails = []
258 parameters = []
259 trailed_ids = set()
260 # Search starting with the source with the largest number of matches.
261 for id in dict(sorted(counts.items(), key=lambda item: item[1], reverse=True)):
262 # Don't search this point if it was already included in a trail.
263 # The anchor itself is one of the trail's points, but is not one
264 # of its own matches.
265 if counts[id] + 1 < self.config.min_points or id in trailed_ids:
266 continue
268 self.log.debug("id=%d at %.1f,%.1f has %d matches within %d pixels.",
269 id,
270 per_id[id][0].first.getX(),
271 per_id[id][0].first.getY(),
272 counts[id],
273 self.config.radius)
274 if (trail := self._search_one(per_id[id], good_catalog)) is not None:
275 trail, result = trail
276 # Check that we didn't already find this trail.
277 n_new = len(set(trail["id"]).difference(trailed_ids))
278 if n_new > 0: 278 ↛ 261line 278 didn't jump to line 261 because the condition on line 278 was always true
279 self.log.info("Found %.1f pixel length trail with %d points, "
280 "%d not in any other trail (angle=%.2f deg, center=%.1f,%.1f)",
281 result.length, len(trail), n_new, math.degrees(result.angle),
282 result.center_x, result.center_y)
283 trails.append(trail)
284 trailed_ids.update(trail["id"])
285 parameters.append(result)
287 self.log.info("Found %d glint trails containing %d total sources.",
288 len(trails), len(trailed_ids))
289 return lsst.pipe.base.Struct(trails=trails,
290 trailed_ids=trailed_ids,
291 parameters=parameters)
293 def _select_good_sources(self, catalog):
294 """Return sources that could possibly be in a glint trail, i.e. ones
295 that do not have bad flags set.
297 Parameters
298 ----------
299 catalog : `lsst.afw.table.SourceCatalog`
300 Original catalog to be selected from.
302 Returns
303 -------
304 good_catalog : `lsst.afw.table.SourceCatalog`
305 Catalog that has had bad sources removed.
306 """
307 bad = np.zeros(len(catalog), dtype=bool)
308 for flag in self.config.bad_flags:
309 bad |= catalog[flag]
310 return catalog[~bad]
312 def _search_one(self, matches, catalog):
313 """Search one set of matches for a possible trail.
315 Parameters
316 ----------
317 matches : `list` [`lsst.afw.table.Match`]
318 Matches for one anchor source to search for lines.
319 catalog : `lsst.afw.SourceCatalog`
320 Catalog of all sources, to refit lines to.
322 Returns
323 -------
324 trail, result : `tuple` or None
325 If the no trails matching the criteria are found, return None,
326 otherwise return a tuple of the sources in the trail and the
327 trail parameters.
328 """
329 components = collections.defaultdict(list)
330 # Normalized distances from the first record to all the others.
331 xy_deltas = {pair.second["id"]: (pair.second.getX() - pair.first.getX(),
332 pair.second.getY() - pair.first.getY()) for pair in matches}
334 # Find all sets of pairs from this anchor that could lie on a line.
335 for i, (id1, pair1) in enumerate(xy_deltas.items()):
336 distance = math.sqrt(pair1[0]**2 + pair1[1]**2)
337 for j, (id2, pair2) in enumerate(xy_deltas.items()):
338 if i == j:
339 continue
340 delta = abs(pair1[0] * pair2[1] - pair1[1] * pair2[0])
341 # 2x threshold to search more broadly; will be refined later.
342 if delta / distance < 2 * self.config.threshold:
343 components[i].append(j)
345 # There are no lines with at least 3 components.
346 if len(components) == 0:
347 return None
349 longest, value = max(components.items(), key=lambda x: len(x[1]))
350 n_points = len(value)
351 n_points += 2 # to account for the base source and the first pair
352 if n_points < self.config.min_points:
353 return None
355 candidate = [longest] + components[longest]
356 trail, result = self._other_points(n_points, candidate, matches, catalog)
358 if trail is None or len(trail) < self.config.min_points:
359 return None
360 return trail, result
362 def _other_points(self, n_points, indexes, matches, catalog):
363 """Find all catalog records that could lie on this line.
365 Parameters
366 ----------
367 n_points : `int`
368 Number of sources in this candidate trail.
369 indexes : `list` [`int`]
370 Indexes into matches on this candidate trail.
371 matches : `list` [`lsst.afw.table.Match`]
372 Matches for one anchor sources to search for lines.
373 catalog : `lsst.afw.SourceCatalog`
374 Catalog of all sources, to refit lines to.
376 Returns
377 -------
378 trail : `lsst.afw.table.SourceCatalog`
379 Sources that are in the fitted trail.
380 result : `GlintTrailParameters`
381 Parameters of the fitted trail.
382 """
384 def extract(line, inliers, x, y, prefix=""):
385 """Extract values from the fit and log and return them."""
386 x = x[inliers]
387 y = y[inliers]
388 stderr = math.sqrt(np.mean(line.distance(x, y)**2))
389 self.log.debug("%s fit: stderr=%.6f, inliers/total=%d/%d",
390 prefix, stderr, len(x), len(inliers))
391 # Simple O(N^2) search for longest distance; there will never be
392 # enough points in a trail a for "faster" approach to be worth it.
393 length = max(scipy.spatial.distance.pdist(np.column_stack((x, y))), default=0)
394 dx, dy = line.direction
395 # The slope and angle do not depend on the direction's sign.
396 if dx == 0: 396 ↛ 397line 396 didn't jump to line 397 because the condition on line 396 was never true
397 m, b = math.inf, math.nan
398 else:
399 m = dy/dx
400 b = line.center[1] - m*line.center[0]
401 angle = math.atan(m)
402 return GlintTrailParameters(slope=m, intercept=b, stderr=stderr, length=length, angle=angle,
403 center_x=line.center[0], center_y=line.center[1])
405 x = np.empty(n_points)
406 x[0] = matches[0].first.getX()
407 x[1:] = [matches[i].second.getX() for i in indexes]
408 y = np.empty(n_points)
409 y[0] = matches[0].first.getY()
410 y[1:] = [matches[i].second.getY() for i in indexes]
412 line, inliers = _ransac_line(x, y, self.config.threshold, min_samples=2, seed=self.config.seed)
413 if line is None: 413 ↛ 414line 413 didn't jump to line 414 because the condition on line 413 was never true
414 self.log.info("Glint trail interpolation could not find a valid fit.")
415 return None, None
416 extract(line, inliers, x, y, prefix="preliminary")
417 # Reject trails that have too many outliers after the first fit.
418 if (n_inliers := np.sum(inliers)) < self.config.min_points:
419 self.log.debug("Candidate trail rejected with %d < %d points.",
420 n_inliers, self.config.min_points)
421 return None, None
423 # Find all points that are close to this line and refit with them.
424 x = catalog["slot_Centroid_x"]
425 y = catalog["slot_Centroid_y"]
426 # 2x threshold to search more broadly: outlier rejection may change
427 # the line parameters some and we want to grab all candidates here.
428 candidates = line.distance(x, y) < 2*self.config.threshold
429 x = x[candidates]
430 y = y[candidates]
431 # min_samples>2 should make the fit more stable.
432 line, inliers = _ransac_line(x, y, self.config.threshold, min_samples=3, seed=self.config.seed)
433 if line is None: 433 ↛ 434line 433 didn't jump to line 434 because the condition on line 433 was never true
434 self.log.info("Glint trail interpolation could not find a valid fit.")
435 return None, None
436 result = extract(line, inliers, x, y, prefix="final")
438 return catalog[candidates][inliers], result