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
maskStreaks.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__ = ["MaskStreaksConfig", "MaskStreaksTask"]
23
24import numpy as np
25import scipy
26import textwrap
27import copy
28from skimage.feature import canny
29from sklearn.cluster import KMeans
30from dataclasses import dataclass
31import astropy.units as u
32
33from lsst.afw.geom import SpanSet
34import lsst.pex.config as pexConfig
35import lsst.pipe.base as pipeBase
36import lsst.kht
37from lsst.utils.timer import timeMethod
38
39
40@dataclass
41class Line:
42 """A simple data class to describe a line profile. The parameter `rho`
43 describes the distance from the center of the image, `theta` describes
44 the angle, `sigma` describes the width of the line, `reducedChi2` gives
45 the reduced chi2 of the fit, and `modelMaximum` gives the peak value of the
46 fit line profile.
47 """
48
49 rho: float
50 theta: float
51 sigma: float = 0
52 reducedChi2: float = 0
53 modelMaximum: float = 0
54
55
57 """Collection of `Line` objects.
58
59 Parameters
60 ----------
61 rhos : `np.ndarray`
62 Array of `Line` rho parameters.
63 thetas : `np.ndarray`
64 Array of `Line` theta parameters.
65 sigmas : `np.ndarray`, optional
66 Array of `Line` sigma parameters.
67 """
68
69 def __init__(self, rhos, thetas, sigmas=None):
70 if sigmas is None:
71 sigmas = np.zeros(len(rhos))
72
73 self._lines = [Line(rho, theta, sigma) for (rho, theta, sigma) in
74 zip(rhos, thetas, sigmas)]
75
76 def __len__(self):
77 return len(self._lines)
78
79 def __getitem__(self, index):
80 return self._lines[index]
81
82 def __iter__(self):
83 return iter(self._lines)
84
85 def __repr__(self):
86 joinedString = ", ".join(str(line) for line in self._lines)
87 return textwrap.shorten(joinedString, width=160, placeholder="...")
88
89 @property
90 def rhos(self):
91 return np.array([line.rho for line in self._lines])
92
93 @property
94 def thetas(self):
95 return np.array([line.theta for line in self._lines])
96
97 @property
98 def sigmas(self):
99 return np.array([line.sigma for line in self._lines])
100
101 def append(self, newLine):
102 """Add line to current collection of lines.
103
104 Parameters
105 ----------
106 newLine : `Line`
107 `Line` to add to current collection of lines
108 """
109 self._lines.append(copy.copy(newLine))
110
111
113 """Construct and/or fit a model for a linear streak.
114
115 This assumes a simple model for a streak, in which the streak
116 follows a straight line in pixels space, with a Moffat-shaped profile. The
117 model is fit to data using a Newton-Raphson style minimization algorithm.
118 The initial guess for the line parameters is assumed to be fairly accurate,
119 so only a narrow band of pixels around the initial line estimate is used in
120 fitting the model, which provides a significant speed-up over using all the
121 data. The class can also be used just to construct a model for the data with
122 a line following the given coordinates.
123
124 Parameters
125 ----------
126 data : `np.ndarray`
127 2d array of data.
128 weights : `np.ndarray`
129 2d array of weights.
130 line : `Line`, optional
131 Guess for position of line. Data far from line guess is masked out.
132 Defaults to None, in which case only data with `weights` = 0 is masked
133 out.
134 """
135
136 def __init__(self, data, weights, line=None):
137 self.data = data
138 self.weights = weights
139 self._ymax, self._xmax = data.shape
140 self._dtype = data.dtype
141 xrange = np.arange(self._xmax) - self._xmax / 2.
142 yrange = np.arange(self._ymax) - self._ymax / 2.
143 self._rhoMax = ((0.5 * self._ymax)**2 + (0.5 * self._xmax)**2)**0.5
144 self._xmesh, self._ymesh = np.meshgrid(xrange, yrange)
145 self.mask = (weights != 0)
146
147 self._initLine = line
148 self.setLineMask(line, maxStreakWidth=0, nSigmaMask=5)
149
150 def getLineXY(self, line):
151 """Return the pixel coordinates of the ends of the line.
152
153 Parameters
154 ----------
155 line : `Line`
156 Line for which to find the endpoints.
157
158 Returns
159 -------
160 boxIntersections : `np.ndarray`
161 (x, y) coordinates of the start and endpoints of the line.
162 """
163 theta = line.theta * u.deg
164 # Determine where the line intersects with each edge of the bounding
165 # box.
166 # Bottom:
167 yA = -self._ymax / 2.
168 xA = (line.rho - yA * np.sin(theta)) / np.cos(theta)
169 # Left:
170 xB = -self._xmax / 2.
171 yB = (line.rho - xB * np.cos(theta)) / np.sin(theta)
172 # Top:
173 yC = self._ymax / 2.
174 xC = (line.rho - yC * np.sin(theta)) / np.cos(theta)
175 # Right:
176 xD = self._xmax / 2.
177 yD = (line.rho - xD * np.cos(theta)) / np.sin(theta)
178
179 lineIntersections = np.array([[xA, yA],
180 [xB, yB],
181 [xC, yC],
182 [xD, yD]])
183 lineIntersections[:, 0] += self._xmax / 2.
184 lineIntersections[:, 1] += self._ymax / 2.
185
186 # The line will necessarily intersect with exactly two edges of the
187 # bounding box itself.
188 inBox = ((lineIntersections[:, 0] >= 0) & (lineIntersections[:, 0] <= self._xmax)
189 & (lineIntersections[:, 1] >= 0) & (lineIntersections[:, 1] <= self._ymax))
190 boxIntersections = lineIntersections[inBox]
191
192 return boxIntersections
193
194 def setLineMask(self, line, maxStreakWidth, nSigmaMask, logger=None):
195 """Set mask around the image region near the line.
196
197 Parameters
198 ----------
199 line : `Line`
200 Parameters of line in the image.
201 """
202 if line:
203 # Only fit pixels within nSigmaMask * sigma of the estimated line
204 radtheta = np.deg2rad(line.theta)
205 distance = (np.cos(radtheta) * self._xmesh + np.sin(radtheta) * self._ymesh - line.rho)
206
207 width = 2 * nSigmaMask * line.sigma
208 if (maxStreakWidth > 0) and (maxStreakWidth < width):
209 if logger is not None:
210 logger.info("Line with width %d exceeded maximum configured width of %d",
211 width, maxStreakWidth)
212 width = maxStreakWidth
213 m = (abs(distance) < width/2)
214 self.lineMask = self.mask & m
215 else:
216 self.lineMask = np.copy(self.mask)
217
218 self.lineMaskSize = self.lineMask.sum()
219 self._maskData = self.data[self.lineMask]
220 self._maskWeights = self.weights[self.lineMask]
221 self._mxmesh = self._xmesh[self.lineMask]
222 self._mymesh = self._ymesh[self.lineMask]
223
224 def _makeMaskedProfile(self, line, fitFlux=True):
225 """Construct the line model in the masked region and calculate its
226 derivatives.
227
228 Parameters
229 ----------
230 line : `Line`
231 Parameters of line profile for which to make profile in the masked
232 region.
233 fitFlux : `bool`
234 Fit the amplitude of the line profile to the data.
235
236 Returns
237 -------
238 model : `np.ndarray`
239 Model in the masked region.
240 dModel : `np.ndarray`
241 Derivative of the model in the masked region.
242 """
243 invSigma = line.sigma**-1
244 # Calculate distance between pixels and line
245 radtheta = np.deg2rad(line.theta)
246 costheta = np.cos(radtheta)
247 sintheta = np.sin(radtheta)
248 distance = (costheta * self._mxmesh + sintheta * self._mymesh - line.rho)
249 distanceSquared = distance**2
250
251 # Calculate partial derivatives of distance
252 drad = np.pi / 180
253 dDistanceSqdRho = 2 * distance * (-np.ones_like(self._mxmesh))
254 dDistanceSqdTheta = (2 * distance * (-sintheta * self._mxmesh + costheta * self._mymesh) * drad)
255
256 # Use pixel-line distances to make Moffat profile
257 profile = (1 + distanceSquared * invSigma**2)**-2.5
258 dProfile = -2.5 * (1 + distanceSquared * invSigma**2)**-3.5
259
260 if fitFlux:
261 # Calculate line flux from profile and data
262 flux = ((self._maskWeights * self._maskData * profile).sum()
263 / (self._maskWeights * profile**2).sum())
264 else:
265 # Approximately normalize the line
266 flux = invSigma**-1
267 if np.isnan(flux):
268 flux = 0
269
270 model = flux * profile
271
272 # Calculate model derivatives
273 fluxdProfile = flux * dProfile
274 fluxdProfileInvSigma = fluxdProfile * invSigma**2
275 dModeldRho = fluxdProfileInvSigma * dDistanceSqdRho
276 dModeldTheta = fluxdProfileInvSigma * dDistanceSqdTheta
277 dModeldInvSigma = fluxdProfile * distanceSquared * 2 * invSigma
278
279 dModel = np.array([dModeldRho, dModeldTheta, dModeldInvSigma])
280 return model, dModel
281
282 def makeProfile(self, line, fitFlux=True):
283 """Construct the line profile model.
284
285 Parameters
286 ----------
287 line : `Line`
288 Parameters of the line profile to model.
289 fitFlux : `bool`, optional
290 Fit the amplitude of the line profile to the data.
291
292 Returns
293 -------
294 finalModel : `np.ndarray`
295 Model for line profile.
296 """
297 model, _ = self._makeMaskedProfile(line, fitFlux=fitFlux)
298 finalModel = np.zeros((self._ymax, self._xmax), dtype=self._dtype)
299 finalModel[self.lineMask] = model
300 return finalModel
301
302 def _lineChi2(self, line, grad=True):
303 """Construct the chi2 between the data and the model.
304
305 Parameters
306 ----------
307 line : `Line`
308 `Line` parameters for which to build model and calculate chi2.
309 grad : `bool`, optional
310 Whether or not to return the gradient and hessian.
311
312 Returns
313 -------
314 reducedChi : `float`
315 Reduced chi2 of the model.
316 reducedDChi : `np.ndarray`
317 Derivative of the chi2 with respect to rho, theta, invSigma.
318 reducedHessianChi : `np.ndarray`
319 Hessian of the chi2 with respect to rho, theta, invSigma.
320 """
321 # Calculate chi2
322 model, dModel = self._makeMaskedProfile(line)
323 chi2 = (self._maskWeights * (self._maskData - model)**2).sum()
324 if not grad:
325 return chi2.sum() / self.lineMaskSize
326
327 # Calculate derivative and Hessian of chi2
328 derivChi2 = ((-2 * self._maskWeights * (self._maskData - model))[None, :] * dModel).sum(axis=1)
329 hessianChi2 = (2 * self._maskWeights * dModel[:, None, :] * dModel[None, :, :]).sum(axis=2)
330
331 reducedChi = chi2 / self.lineMaskSize
332 reducedDChi = derivChi2 / self.lineMaskSize
333 reducedHessianChi = hessianChi2 / self.lineMaskSize
334 return reducedChi, reducedDChi, reducedHessianChi
335
336 def fit(self, dChi2Tol=0.1, maxIter=100, log=None):
337 """Perform Newton-Raphson minimization to find line parameters.
338
339 This method takes advantage of having known derivative and Hessian of
340 the multivariate function to quickly and efficiently find the minimum.
341 This is more efficient than the scipy implementation of the Newton-
342 Raphson method, which doesn't take advantage of the Hessian matrix. The
343 method here also performs a line search in the direction of the steepest
344 derivative at each iteration, which reduces the number of iterations
345 needed.
346
347 Parameters
348 ----------
349 dChi2Tol : `float`, optional
350 Change in Chi2 tolerated for fit convergence.
351 maxIter : `int`, optional
352 Maximum number of fit iterations allowed. The fit should converge in
353 ~10 iterations, depending on the value of dChi2Tol, but this
354 maximum provides a backup.
355 log : `lsst.utils.logging.LsstLogAdapter`, optional
356 Logger to use for reporting more details for fitting failures.
357
358 Returns
359 -------
360 outline : `Line`
361 Coordinates, inverse width, and chi2 of fit line.
362 fitFailure : `bool`
363 Boolean where `False` corresponds to a successful fit.
364 """
365 # Do minimization on inverse of sigma to simplify derivatives:
366 x = np.array([self._initLine.rho, self._initLine.theta, self._initLine.sigma**-1])
367
368 dChi2 = 1
369 iter = 0
370 oldChi2 = 0
371 fitFailure = False
372
373 def line_search(c, dx):
374 testx = x - c * dx
375 testLine = Line(testx[0], testx[1], testx[2]**-1)
376 return self._lineChi2(testLine, grad=False)
377
378 while abs(dChi2) > dChi2Tol:
379 line = Line(x[0], x[1], x[2]**-1)
380 chi2, b, A = self._lineChi2(line)
381 if chi2 == 0:
382 break
383 if not np.isfinite(A).all():
384 fitFailure = True
385 if log is not None:
386 log.warning("Hessian matrix has non-finite elements.")
387 break
388 dChi2 = oldChi2 - chi2
389 try:
390 cholesky = scipy.linalg.cho_factor(A)
391 except np.linalg.LinAlgError:
392 fitFailure = True
393 if log is not None:
394 log.warning("Hessian matrix is not invertible.")
395 break
396 dx = scipy.linalg.cho_solve(cholesky, b)
397
398 factor, fmin, _, _ = scipy.optimize.brent(line_search, args=(dx,), full_output=True, tol=0.05)
399 x -= factor * dx
400 if (abs(x[0]) > 1.5 * self._rhoMax) or (iter > maxIter):
401 fitFailure = True
402 break
403 oldChi2 = chi2
404 iter += 1
405
406 outline = Line(x[0], x[1], abs(x[2])**-1, chi2)
407
408 return outline, fitFailure
409
410
411class MaskStreaksConfig(pexConfig.Config):
412 """Configuration parameters for `MaskStreaksTask`.
413 """
414
415 minimumKernelHeight = pexConfig.Field(
416 doc="Minimum height of the streak-finding kernel relative to the tallest kernel",
417 dtype=float,
418 default=0.0,
419 )
420 absMinimumKernelHeight = pexConfig.Field(
421 doc="Minimum absolute height of the streak-finding kernel",
422 dtype=float,
423 default=5,
424 )
425 clusterMinimumSize = pexConfig.Field(
426 doc="Minimum size in pixels of detected clusters",
427 dtype=int,
428 default=50,
429 )
430 clusterMinimumDeviation = pexConfig.Field(
431 doc="Allowed deviation (in pixels) from a straight line for a detected "
432 "line",
433 dtype=int,
434 default=2,
435 )
436 delta = pexConfig.Field(
437 doc="Stepsize in angle-radius parameter space",
438 dtype=float,
439 default=0.2,
440 )
441 nSigma = pexConfig.Field(
442 doc="Number of sigmas from center of kernel to include in voting "
443 "procedure",
444 dtype=float,
445 default=2,
446 )
447 nSigmaMask = pexConfig.Field(
448 doc="Number of sigmas from center of kernel to mask",
449 dtype=float,
450 default=5,
451 )
452 rhoBinSize = pexConfig.Field(
453 doc="Binsize in pixels for position parameter rho when finding "
454 "clusters of detected lines",
455 dtype=float,
456 default=30,
457 )
458 thetaBinSize = pexConfig.Field(
459 doc="Binsize in degrees for angle parameter theta when finding "
460 "clusters of detected lines",
461 dtype=float,
462 default=2,
463 )
464 invSigma = pexConfig.Field(
465 doc="Inverse of the Moffat sigma parameter (in units of pixels)"
466 "describing the profile of the streak",
467 dtype=float,
468 default=10.**-1,
469 )
470 footprintThreshold = pexConfig.Field(
471 doc="Threshold at which to determine edge of line, in units of "
472 "nanoJanskys",
473 dtype=float,
474 default=0.01
475 )
476 dChi2Tolerance = pexConfig.Field(
477 doc="Absolute difference in Chi2 between iterations of line profile"
478 "fitting that is acceptable for convergence",
479 dtype=float,
480 default=0.1
481 )
482 maxFitIter = pexConfig.Field(
483 doc="Maximum number of line profile fitting iterations that is "
484 "acceptable for convergence",
485 dtype=int,
486 default=100
487 )
488 detectedMaskPlane = pexConfig.Field(
489 doc="Name of mask with pixels above detection threshold, used for first"
490 "estimate of streak locations",
491 dtype=str,
492 default="DETECTED"
493 )
494 onlyMaskDetected = pexConfig.Field(
495 doc=("If true, only propagate the part of the streak mask that "
496 "overlaps with the detection mask."),
497 dtype=bool,
498 default=True,
499 )
500 streaksMaskPlane = pexConfig.Field(
501 doc="Name of mask plane holding detected streaks",
502 dtype=str,
503 default="STREAK"
504 )
505 badMaskPlanes = pexConfig.ListField(
506 doc=("Names of mask plane regions to ignore entirely when doing streak"
507 " detection"),
508 dtype=str,
509 default=("NO_DATA", "INTRP", "BAD", "SAT", "EDGE"),
510 )
511 maxStreakWidth = pexConfig.Field(
512 doc="Maximum width in pixels to allow for masking a streak."
513 "The fit streak parameters will not be modified, and a warning will "
514 "be issued if the fitted width is larger than this value."
515 "Set to 0 to disable.",
516 dtype=float,
517 default=0.,
518 )
519
520
521class MaskStreaksTask(pipeBase.Task):
522 """Find streaks or other straight lines in image data.
523
524 Nearby objects passing through the field of view of the telescope leave a
525 bright trail in images. This class uses the Kernel Hough Transform (KHT)
526 (Fernandes and Oliveira, 2007), implemented in `lsst.houghtransform`. The
527 procedure works by taking a binary image, either provided as put or produced
528 from the input data image, using a Canny filter to make an image of the
529 edges in the original image, then running the KHT on the edge image. The KHT
530 identifies clusters of non-zero points, breaks those clusters of points into
531 straight lines, keeps clusters with a size greater than the user-set
532 threshold, then performs a voting procedure to find the best-fit coordinates
533 of any straight lines. Given the results of the KHT algorithm, clusters of
534 lines are identified and grouped (generally these correspond to the two
535 edges of a strea) and a profile is fit to the streak in the original
536 (non-binary) image.
537 """
538
539 ConfigClass = MaskStreaksConfig
540 _DefaultName = "maskStreaks"
541
542 @timeMethod
543 def find(self, maskedImage):
544 """Find streaks in a masked image.
545
546 Parameters
547 ----------
548 maskedImage : `lsst.afw.image.maskedImage`
549 The image in which to search for streaks.
550
551 Returns
552 -------
553 result : `lsst.pipe.base.Struct`
554 Results as a struct with attributes:
555
556 ``originalLines``
557 Lines identified by kernel hough transform.
558 ``lineClusters``
559 Lines grouped into clusters in rho-theta space.
560 ``lines``
561 Final result for lines after line-profile fit.
562 ``mask``
563 2-d boolean mask where detected lines are True.
564 """
565 mask = maskedImage.mask
566 detectionMask = (mask.array & mask.getPlaneBitMask(self.config.detectedMaskPlane))
567
568 initEdges = self._cannyFilter(detectionMask)
569 # Ignore regions with known bad masks, adding a one-pixel buffer around
570 # each to ensure that the edges of bad regions are also ignored.
571 badPixelMask = mask.getPlaneBitMask(self.config.badMaskPlanes)
572 badMaskSpanSet = SpanSet.fromMask(mask, badPixelMask).split()
573 for sset in badMaskSpanSet:
574 sset_dilated = sset.dilated(1)
575 sset_dilated.clippedTo(mask.getBBox()).setMask(mask, mask.getPlaneBitMask("BAD"))
576 dilatedBadMask = (mask.array & badPixelMask) > 0
577 self.edges = initEdges & ~dilatedBadMask
578 self.lines = self._runKHT(self.edges)
579
580 if len(self.lines) == 0:
581 lineMask = np.zeros(detectionMask.shape, dtype=bool)
582 fitLines = LineCollection([], [])
583 clusters = LineCollection([], [])
584 else:
585 clusters = self._findClusters(self.lines)
586 fitLines, lineMask = self._fitProfile(clusters, maskedImage)
587
588 if self.config.onlyMaskDetected:
589 # The output mask is the intersection of the fit streaks and the image detections
590 lineMask &= detectionMask.astype(bool)
591
592 return pipeBase.Struct(
593 lines=fitLines,
594 lineClusters=clusters,
595 originalLines=self.lines,
596 mask=lineMask,
597 )
598
599 @timeMethod
600 def run(self, maskedImage):
601 """Find and mask streaks in a masked image.
602
603 Finds streaks in the image and modifies maskedImage in place by adding a
604 mask plane with any identified streaks.
605
606 Parameters
607 ----------
608 maskedImage : `lsst.afw.image.Exposure` or `lsst.afw.image.maskedImage`
609 The image in which to search for streaks. The mask detection plane
610 corresponding to `config.detectedMaskPlane` must be set with the
611 detected pixels. The mask will have a plane added with any detected
612 streaks, and with the mask plane name set by
613 self.config.streaksMaskPlane.
614
615 Returns
616 -------
617 result : `lsst.pipe.base.Struct`
618 Results as a struct with attributes:
619
620 ``originalLines``
621 Lines identified by kernel hough transform.
622 ``lineClusters``
623 Lines grouped into clusters in rho-theta space.
624 ``lines``
625 Final result for lines after line-profile fit.
626 """
627 streaks = self.find(maskedImage)
628
629 if (self.config.streaksMaskPlane != "STREAK") and \
630 (self.config.streaksMaskPlane not in maskedImage.mask.getMaskPlaneDict()):
631 maskedImage.mask.addMaskPlane(self.config.streaksMaskPlane)
632 maskedImage.mask.array[streaks.mask] |= maskedImage.mask.getPlaneBitMask(self.config.streaksMaskPlane)
633
634 return pipeBase.Struct(
635 lines=streaks.lines,
636 lineClusters=streaks.lineClusters,
637 originalLines=streaks.originalLines,
638 )
639
640 def _cannyFilter(self, image):
641 """Apply a canny filter to the data in order to detect edges.
642
643 Parameters
644 ----------
645 image : `np.ndarray`
646 2-d image data on which to run filter.
647
648 Returns
649 -------
650 cannyData : `np.ndarray`
651 2-d image of edges found in input image.
652 """
653 # Ensure that the pixels are zero or one. Change the datatype to
654 # np.float64 to be compatible with the Canny filter routine.
655 filterData = (image > 0).astype(np.float64)
656 return canny(filterData, use_quantiles=True, sigma=0.1)
657
658 def _runKHT(self, image):
659 """Run Kernel Hough Transform on image.
660
661 Parameters
662 ----------
663 image : `np.ndarray`
664 2-d image data on which to detect lines.
665
666 Returns
667 -------
668 result : `LineCollection`
669 Collection of detected lines, with their detected rho and theta
670 coordinates.
671 """
672 lines = lsst.kht.find_lines(image, self.config.clusterMinimumSize,
673 self.config.clusterMinimumDeviation, self.config.delta,
674 self.config.minimumKernelHeight, self.config.nSigma,
675 self.config.absMinimumKernelHeight)
676 self.log.info("The Kernel Hough Transform detected %s line(s)", len(lines))
677
678 return LineCollection(lines.rho, lines.theta)
679
680 def _findClusters(self, lines):
681 """Group lines that are close in parameter space and likely describe
682 the same streak.
683
684 Parameters
685 ----------
686 lines : `LineCollection`
687 Collection of lines to group into clusters.
688
689 Returns
690 -------
691 result : `LineCollection`
692 Average `Line` for each cluster of `Line`s in the input
693 `LineCollection`.
694 """
695 # Scale variables by threshold bin-size variable so that rho and theta
696 # are on the same scale. Since the clustering algorithm below stops when
697 # the standard deviation <= 1, after rescaling each cluster will have a
698 # standard deviation at or below the bin-size.
699 x = lines.rhos / self.config.rhoBinSize
700 y = lines.thetas / self.config.thetaBinSize
701 X = np.array([x, y]).T
702 nClusters = 1
703
704 # Put line parameters in clusters by starting with all in one, then
705 # subdividing until the parameters of each cluster have std dev=1.
706 # If nClusters == len(lines), each line will have its own 'cluster', so
707 # the standard deviations of each cluster must be zero and the loop
708 # is guaranteed to stop.
709 while True:
710 kmeans = KMeans(n_clusters=nClusters, n_init='auto').fit(X)
711 clusterStandardDeviations = np.zeros((nClusters, 2))
712 for c in range(nClusters):
713 inCluster = X[kmeans.labels_ == c]
714 clusterStandardDeviations[c] = np.std(inCluster, axis=0)
715 # Are the rhos and thetas in each cluster all below the threshold?
716 if (clusterStandardDeviations <= 1).all():
717 break
718 nClusters += 1
719
720 # The cluster centers are final line estimates
721 finalClusters = kmeans.cluster_centers_.T
722
723 # Rescale variables:
724 finalRhos = finalClusters[0] * self.config.rhoBinSize
725 finalThetas = finalClusters[1] * self.config.thetaBinSize
726 result = LineCollection(finalRhos, finalThetas)
727 self.log.info("Lines were grouped into %s potential streak(s)", len(finalRhos))
728
729 return result
730
731 def _fitProfile(self, lines, maskedImage):
732 """Fit the profile of the streak.
733
734 Given the initial parameters of detected lines, fit a model for the
735 streak to the original (non-binary image). The assumed model is a
736 straight line with a Moffat profile.
737
738 Parameters
739 ----------
740 lines : `LineCollection`
741 Collection of guesses for `Line`s detected in the image.
742 maskedImage : `lsst.afw.image.maskedImage`
743 Original image to be used to fit profile of streak.
744
745 Returns
746 -------
747 lineFits : `LineCollection`
748 Collection of `Line` profiles fit to the data.
749 finalMask : `np.ndarray`
750 2d mask array with detected streaks=1.
751 """
752 data = maskedImage.image.array
753 weights = maskedImage.variance.array**-1
754 mask = maskedImage.mask
755 badPixelMask = mask.getPlaneBitMask(self.config.badMaskPlanes)
756 badMask = (mask.array & badPixelMask) > 0
757 # Mask out any pixels with non-finite weights
758 weights[~np.isfinite(weights) | ~np.isfinite(data)] = 0
759 weights[badMask] = 0
760
761 lineFits = LineCollection([], [])
762 finalLineMasks = [np.zeros(data.shape, dtype=bool)]
763 nFinalLines = 0
764 nFitFailures = 0
765 for line in lines:
766 line.sigma = self.config.invSigma**-1
767 lineModel = LineProfile(data, weights, line=line)
768 # Skip any lines that do not cover any data (sometimes happens because of chip gaps)
769 if lineModel.lineMaskSize == 0:
770 continue
771
772 fit, fitFailure = lineModel.fit(dChi2Tol=self.config.dChi2Tolerance, log=self.log,
773 maxIter=self.config.maxFitIter)
774
775 # Initial estimate should be quite close: fit is deemed unsuccessful if rho or theta
776 # change more than the allowed bin in rho or theta:
777 if ((abs(fit.rho - line.rho) > 2 * self.config.rhoBinSize)
778 or (abs(fit.theta - line.theta) > 2 * self.config.thetaBinSize)):
779 fitFailure = True
780 self.log.debug("Streak fit moved too far from initial estimate. Line will be dropped.")
781
782 if fitFailure:
783 nFitFailures += 1
784 continue
785
786 # Make mask
787 lineModel.setLineMask(fit, self.config.maxStreakWidth, self.config.nSigmaMask, logger=self.log)
788 finalModel = lineModel.makeProfile(fit)
789 # Take absolute value, as streaks are allowed to be negative
790 finalModelMax = abs(finalModel).max()
791 finalLineMask = abs(finalModel) > self.config.footprintThreshold
792 # Drop this line if the model profile is below the footprint
793 # threshold
794 if not finalLineMask.any():
795 continue
796 fit.modelMaximum = finalModelMax
797 lineFits.append(fit)
798 finalLineMasks.append(finalLineMask)
799 nFinalLines += 1
800
801 if nFitFailures > 0:
802 self.log.info("Streak profile could not be fit for %d out of %d detected lines.", nFitFailures,
803 len(lines))
804 finalMask = np.array(finalLineMasks).any(axis=0)
805
806 return lineFits, finalMask
__init__(self, rhos, thetas, sigmas=None)
setLineMask(self, line, maxStreakWidth, nSigmaMask, logger=None)
__init__(self, data, weights, line=None)
Definition __init__.py:1