23 __all__ = [
"MaskStreaksConfig",
"MaskStreaksTask",
"setDetectionMask"]
28 from lsst.utils.timer
import timeMethod
34 from skimage.feature
import canny
35 from sklearn.cluster
import KMeans
37 from dataclasses
import dataclass
40 def setDetectionMask(maskedImage, forceSlowBin=False, binning=None, detectedPlane="DETECTED",
41 badMaskPlanes=(
"NO_DATA",
"INTRP",
"BAD",
"SAT",
"EDGE"), detectionThreshold=5):
42 """Make detection mask and set the mask plane
44 Creat a binary image from a masked image by setting all data with signal-to-
45 noise below some threshold to zero, and all data above the threshold to one.
46 If the binning parameter has been set, this procedure will be preceded by a
47 weighted binning of the data in order to smooth the result, after which the
48 result is scaled back to the original dimensions. Set the detection mask
49 plane with this binary image.
53 maskedImage : `lsst.afw.image.maskedImage`
54 Image to be (optionally) binned and converted
55 forceSlowBin : bool (optional)
56 Force usage of slower binning method to check that the two methods
58 binning : int (optional)
59 Number of pixels by which to bin image
60 detectedPlane : str (optional)
61 Name of mask with pixels that were detected above threshold in image
62 badMaskPlanes : set (optional)
63 Names of masks with pixels that are rejected
64 detectionThreshold : float (optional)
65 Boundary in signal-to-noise between non-detections and detections for
66 making a binary image from the original input image
68 data = maskedImage.image.array
69 weights = 1 / maskedImage.variance.array
70 mask = maskedImage.getMask()
72 detectionMask = ((mask.array & mask.getPlaneBitMask(detectedPlane)))
73 badPixelMask = mask.getPlaneBitMask(badMaskPlanes)
74 badMask = (mask.array & badPixelMask) > 0
75 fitMask = detectionMask.astype(bool) & ~badMask
77 fitData = np.copy(data)
79 fitWeights = np.copy(weights)
80 fitWeights[~fitMask] = 0
84 ymax, xmax = fitData.shape
85 if (ymax % binning == 0)
and (xmax % binning == 0)
and (
not forceSlowBin):
87 binNumeratorReshape = (fitData * fitWeights).reshape(ymax // binning, binning,
88 xmax // binning, binning)
89 binDenominatorReshape = fitWeights.reshape(binNumeratorReshape.shape)
90 binnedNumerator = binNumeratorReshape.sum(axis=3).sum(axis=1)
91 binnedDenominator = binDenominatorReshape.sum(axis=3).sum(axis=1)
94 warnings.warn(
'Using slow binning method--consider choosing a binsize that evenly divides '
95 f
'into the image size, so that {ymax} mod binning == 0 '
96 f
'and {xmax} mod binning == 0')
97 xarray = np.arange(xmax)
98 yarray = np.arange(ymax)
99 xmesh, ymesh = np.meshgrid(xarray, yarray)
100 xbins = np.arange(0, xmax + binning, binning)
101 ybins = np.arange(0, ymax + binning, binning)
102 numerator = fitWeights * fitData
103 binnedNumerator, *_ = scipy.stats.binned_statistic_2d(ymesh.ravel(), xmesh.ravel(),
104 numerator.ravel(), statistic=
'sum',
106 binnedDenominator, *_ = scipy.stats.binned_statistic_2d(ymesh.ravel(), xmesh.ravel(),
107 fitWeights.ravel(), statistic=
'sum',
109 binnedData = np.zeros(binnedNumerator.shape)
110 ind = binnedDenominator != 0
111 np.divide(binnedNumerator, binnedDenominator, out=binnedData, where=ind)
112 binnedWeight = binnedDenominator
113 binMask = (binnedData * binnedWeight**0.5) > detectionThreshold
114 tmpOutputMask = binMask.repeat(binning, axis=0)[:ymax]
115 outputMask = tmpOutputMask.repeat(binning, axis=1)[:, :xmax]
117 outputMask = (fitData * fitWeights**0.5) > detectionThreshold
120 maskedImage.mask.array &= ~maskedImage.mask.getPlaneBitMask(detectedPlane)
123 maskedImage.mask.array[outputMask] |= maskedImage.mask.getPlaneBitMask(detectedPlane)
128 """A simple data class to describe a line profile. The parameter `rho`
129 describes the distance from the center of the image, `theta` describes
130 the angle, and `sigma` describes the width of the line.
138 """Collection of `Line` objects.
143 Array of `Line` rho parameters
145 Array of `Line` theta parameters
146 sigmas : np.ndarray (optional)
147 Array of `Line` sigma parameters
152 sigmas = np.zeros(len(rhos))
154 self.
_lines_lines = [
Line(rho, theta, sigma)
for (rho, theta, sigma)
in
155 zip(rhos, thetas, sigmas)]
158 return len(self.
_lines_lines)
161 return self.
_lines_lines[index]
167 joinedString =
", ".join(str(line)
for line
in self.
_lines_lines)
168 return textwrap.shorten(joinedString, width=160, placeholder=
"...")
172 return np.array([line.rho
for line
in self.
_lines_lines])
176 return np.array([line.theta
for line
in self.
_lines_lines])
179 """Add line to current collection of lines.
184 `Line` to add to current collection of lines
190 """Construct and/or fit a model for a linear streak.
192 This assumes a simple model for a streak, in which the streak
193 follows a straight line in pixels space, with a Moffat-shaped profile. The
194 model is fit to data using a Newton-Raphson style minimization algorithm.
195 The initial guess for the line parameters is assumed to be fairly accurate,
196 so only a narrow band of pixels around the initial line estimate is used in
197 fitting the model, which provides a significant speed-up over using all the
198 data. The class can also be used just to construct a model for the data with
199 a line following the given coordinates.
207 line : `Line` (optional)
208 Guess for position of line. Data far from line guess is masked out.
209 Defaults to None, in which case only data with `weights` = 0 is masked
216 self._ymax, self.
_xmax_xmax = data.shape
217 self.
_dtype_dtype = data.dtype
218 xrange = np.arange(self.
_xmax_xmax) - self.
_xmax_xmax / 2.
219 yrange = np.arange(self._ymax) - self._ymax / 2.
220 self.
_rhoMax_rhoMax = ((0.5 * self._ymax)**2 + (0.5 * self.
_xmax_xmax)**2)**0.5
221 self._xmesh, self.
_ymesh_ymesh = np.meshgrid(xrange, yrange)
228 """Set mask around the image region near the line
233 Parameters of line in the image
237 radtheta = np.deg2rad(line.theta)
238 distance = (np.cos(radtheta) * self._xmesh + np.sin(radtheta) * self.
_ymesh_ymesh - line.rho)
239 m = (
abs(distance) < 5 * line.sigma)
250 def _makeMaskedProfile(self, line, fitFlux=True):
251 """Construct the line model in the masked region and calculate its
257 Parameters of line profile for which to make profile in the masked
260 Fit the amplitude of the line profile to the data
265 Model in the masked region
267 Derivative of the model in the masked region
269 invSigma = line.sigma**-1
271 radtheta = np.deg2rad(line.theta)
272 costheta = np.cos(radtheta)
273 sintheta = np.sin(radtheta)
274 distance = (costheta * self.
_mxmesh_mxmesh + sintheta * self.
_mymesh_mymesh - line.rho)
275 distanceSquared = distance**2
279 dDistanceSqdRho = 2 * distance * (-np.ones_like(self.
_mxmesh_mxmesh))
280 dDistanceSqdTheta = (2 * distance * (-sintheta * self.
_mxmesh_mxmesh + costheta * self.
_mymesh_mymesh) * drad)
283 profile = (1 + distanceSquared * invSigma**2)**-2.5
284 dProfile = -2.5 * (1 + distanceSquared * invSigma**2)**-3.5
296 model = flux * profile
299 fluxdProfile = flux * dProfile
300 fluxdProfileInvSigma = fluxdProfile * invSigma**2
301 dModeldRho = fluxdProfileInvSigma * dDistanceSqdRho
302 dModeldTheta = fluxdProfileInvSigma * dDistanceSqdTheta
303 dModeldInvSigma = fluxdProfile * distanceSquared * 2 * invSigma
305 dModel = np.array([dModeldRho, dModeldTheta, dModeldInvSigma])
309 """Construct the line profile model
314 Parameters of the line profile to model
315 fitFlux : bool (optional)
316 Fit the amplitude of the line profile to the data
320 finalModel : np.ndarray
321 Model for line profile
324 finalModel = np.zeros((self._ymax, self.
_xmax_xmax), dtype=self.
_dtype_dtype)
325 finalModel[self.
lineMasklineMask] = model
328 def _lineChi2(self, line, grad=True):
329 """Construct the chi2 between the data and the model
334 `Line` parameters for which to build model and calculate chi2
335 grad : bool (optional)
336 Whether or not to return the gradient and hessian
341 Reduced chi2 of the model
342 reducedDChi : np.ndarray
343 Derivative of the chi2 with respect to rho, theta, invSigma
344 reducedHessianChi : np.ndarray
345 Hessian of the chi2 with respect to rho, theta, invSigma
354 derivChi2 = ((-2 * self.
_maskWeights_maskWeights * (self.
_maskData_maskData - model))[
None, :] * dModel).sum(axis=1)
355 hessianChi2 = (2 * self.
_maskWeights_maskWeights * dModel[:,
None, :] * dModel[
None, :, :]).sum(axis=2)
359 reducedHessianChi = hessianChi2 / self.
lineMaskSizelineMaskSize
360 return reducedChi, reducedDChi, reducedHessianChi
362 def fit(self, dChi2Tol=0.1, maxIter=100):
363 """Perform Newton-Raphson minimization to find line parameters
365 This method takes advantage of having known derivative and Hessian of
366 the multivariate function to quickly and efficiently find the minimum.
367 This is more efficient than the scipy implementation of the Newton-
368 Raphson method, which doesn't take advantage of the Hessian matrix. The
369 method here also performs a line search in the direction of the steepest
370 derivative at each iteration, which reduces the number of iterations
375 dChi2Tol : float (optional)
376 Change in Chi2 tolerated for fit convergence
377 maxIter : int (optional)
378 Maximum number of fit iterations allowed. The fit should converge in
379 ~10 iterations, depending on the value of dChi2Tol, but this
380 maximum provides a backup.
385 Coordinates and inverse width of fit line
387 Reduced Chi2 of model fit to data
389 Boolean where `False` corresponds to a successful fit
399 def line_search(c, dx):
401 testLine =
Line(testx[0], testx[1], testx[2]**-1)
402 return self.
_lineChi2_lineChi2(testLine, grad=
False)
404 while abs(dChi2) > dChi2Tol:
405 line =
Line(x[0], x[1], x[2]**-1)
406 chi2, b, A = self.
_lineChi2_lineChi2(line)
409 if not np.isfinite(A).
all():
413 dChi2 = oldChi2 - chi2
414 cholesky = scipy.linalg.cho_factor(A)
415 dx = scipy.linalg.cho_solve(cholesky, b)
417 factor, fmin, _, _ = scipy.optimize.brent(line_search, args=(dx,), full_output=
True, tol=0.05)
419 if (
abs(x[0]) > 1.5 * self.
_rhoMax_rhoMax)
or (iter > maxIter):
425 outline =
Line(x[0], x[1],
abs(x[2])**-1)
427 return outline, chi2, fitFailure
431 """Configuration parameters for `MaskStreaksTask`
433 minimumKernelHeight = pexConfig.Field(
434 doc=
"Minimum height of the streak-finding kernel relative to the tallest kernel",
438 absMinimumKernelHeight = pexConfig.Field(
439 doc=
"Minimum absolute height of the streak-finding kernel",
443 clusterMinimumSize = pexConfig.Field(
444 doc=
"Minimum size in pixels of detected clusters",
448 clusterMinimumDeviation = pexConfig.Field(
449 doc=
"Allowed deviation (in pixels) from a straight line for a detected "
454 delta = pexConfig.Field(
455 doc=
"Stepsize in angle-radius parameter space",
459 nSigma = pexConfig.Field(
460 doc=
"Number of sigmas from center of kernel to include in voting "
465 rhoBinSize = pexConfig.Field(
466 doc=
"Binsize in pixels for position parameter rho when finding "
467 "clusters of detected lines",
471 thetaBinSize = pexConfig.Field(
472 doc=
"Binsize in degrees for angle parameter theta when finding "
473 "clusters of detected lines",
477 invSigma = pexConfig.Field(
478 doc=
"Inverse of the Moffat sigma parameter (in units of pixels)"
479 "describing the profile of the streak",
483 footprintThreshold = pexConfig.Field(
484 doc=
"Threshold at which to determine edge of line, in units of the line"
489 dChi2Tolerance = pexConfig.Field(
490 doc=
"Absolute difference in Chi2 between iterations of line profile"
491 "fitting that is acceptable for convergence",
495 detectedMaskPlane = pexConfig.Field(
496 doc=
"Name of mask with pixels above detection threshold, used for first"
497 "estimate of streak locations",
501 streaksMaskPlane = pexConfig.Field(
502 doc=
"Name of mask plane holding detected streaks",
509 """Find streaks or other straight lines in image data.
511 Nearby objects passing through the field of view of the telescope leave a
512 bright trail in images. This class uses the Kernel Hough Transform (KHT)
513 (Fernandes and Oliveira, 2007), implemented in `lsst.houghtransform`. The
514 procedure works by taking a binary image, either provided as put or produced
515 from the input data image, using a Canny filter to make an image of the
516 edges in the original image, then running the KHT on the edge image. The KHT
517 identifies clusters of non-zero points, breaks those clusters of points into
518 straight lines, keeps clusters with a size greater than the user-set
519 threshold, then performs a voting procedure to find the best-fit coordinates
520 of any straight lines. Given the results of the KHT algorithm, clusters of
521 lines are identified and grouped (generally these correspond to the two
522 edges of a strea) and a profile is fit to the streak in the original
526 ConfigClass = MaskStreaksConfig
527 _DefaultName =
"maskStreaks"
531 """Find streaks in a masked image
535 maskedImage : `lsst.afw.image.maskedImage`
536 The image in which to search for streaks.
540 result : `lsst.pipe.base.Struct`
541 Result struct with components:
543 - ``originalLines``: lines identified by kernel hough transform
544 - ``lineClusters``: lines grouped into clusters in rho-theta space
545 - ``lines``: final result for lines after line-profile fit
546 - ``mask``: 2-d boolean mask where detected lines are True
548 mask = maskedImage.getMask()
549 detectionMask = (mask.array & mask.getPlaneBitMask(self.config.detectedMaskPlane))
554 if len(self.
lineslines) == 0:
555 lineMask = np.zeros(detectionMask.shape, dtype=bool)
560 fitLines, lineMask = self.
_fitProfile_fitProfile(clusters, maskedImage)
563 outputMask = lineMask & detectionMask.astype(bool)
565 return pipeBase.Struct(
567 lineClusters=clusters,
568 originalLines=self.
lineslines,
573 def run(self, maskedImage):
574 """Find and mask streaks in a masked image.
576 Finds streaks in the image and modifies maskedImage in place by adding a
577 mask plane with any identified streaks.
581 maskedImage : `lsst.afw.image.maskedImage`
582 The image in which to search for streaks. The mask detection plane
583 corresponding to `config.detectedMaskPlane` must be set with the
588 result : `lsst.pipe.base.Struct`
589 Result struct with components:
591 - ``originalLines``: lines identified by kernel hough transform
592 - ``lineClusters``: lines grouped into clusters in rho-theta space
593 - ``lines``: final result for lines after line-profile fit
595 streaks = self.
findfind(maskedImage)
597 maskedImage.mask.addMaskPlane(self.config.streaksMaskPlane)
598 maskedImage.mask.array[streaks.mask] |= maskedImage.mask.getPlaneBitMask(self.config.streaksMaskPlane)
600 return pipeBase.Struct(
602 lineClusters=streaks.lineClusters,
603 originalLines=streaks.originalLines,
606 def _cannyFilter(self, image):
607 """Apply a canny filter to the data in order to detect edges
612 2-d image data on which to run filter
616 cannyData : `np.ndarray`
617 2-d image of edges found in input image
619 filterData = image.astype(int)
620 return canny(filterData, low_threshold=0, high_threshold=1, sigma=0.1)
622 def _runKHT(self, image):
623 """Run Kernel Hough Transform on image.
628 2-d image data on which to detect lines
632 result : `LineCollection`
633 Collection of detected lines, with their detected rho and theta
636 lines = lsst.kht.find_lines(image, self.config.clusterMinimumSize,
637 self.config.clusterMinimumDeviation, self.config.delta,
638 self.config.minimumKernelHeight, self.config.nSigma,
639 self.config.absMinimumKernelHeight)
643 def _findClusters(self, lines):
644 """Group lines that are close in parameter space and likely describe
649 lines : `LineCollection`
650 Collection of lines to group into clusters
654 result : `LineCollection`
655 Average `Line` for each cluster of `Line`s in the input
662 x = lines.rhos / self.config.rhoBinSize
663 y = lines.thetas / self.config.thetaBinSize
664 X = np.array([x, y]).T
673 kmeans = KMeans(n_clusters=nClusters).fit(X)
674 clusterStandardDeviations = np.zeros((nClusters, 2))
675 for c
in range(nClusters):
676 inCluster = X[kmeans.labels_ == c]
677 clusterStandardDeviations[c] = np.std(inCluster, axis=0)
679 if (clusterStandardDeviations <= 1).
all():
684 finalClusters = kmeans.cluster_centers_.T
687 finalRhos = finalClusters[0] * self.config.rhoBinSize
688 finalThetas = finalClusters[1] * self.config.thetaBinSize
693 def _fitProfile(self, lines, maskedImage):
694 """Fit the profile of the streak.
696 Given the initial parameters of detected lines, fit a model for the
697 streak to the original (non-binary image). The assumed model is a
698 straight line with a Moffat profile.
702 lines : `LineCollection`
703 Collection of guesses for `Line`s detected in the image
704 maskedImage : `lsst.afw.image.maskedImage`
705 Original image to be used to fit profile of streak.
709 lineFits : `LineCollection`
710 Collection of `Line` profiles fit to the data
711 finalMask : `np.ndarray`
712 2d mask array with detected streaks=1.
714 data = maskedImage.image.array
715 weights = maskedImage.variance.array**-1
717 weights[~np.isfinite(weights) | ~np.isfinite(data)] = 0
720 finalLineMasks = [np.zeros(data.shape, dtype=bool)]
722 line.sigma = self.config.invSigma**-1
725 if lineModel.lineMaskSize == 0:
728 fit, chi2, fitFailure = lineModel.fit(dChi2Tol=self.config.dChi2Tolerance)
732 if ((
abs(fit.rho - line.rho) > 2 * self.config.rhoBinSize)
733 or (
abs(fit.theta - line.theta) > 2 * self.config.thetaBinSize)):
740 lineModel.setLineMask(fit)
741 finalModel = lineModel.makeProfile(fit)
743 finalModelMax =
abs(finalModel).
max()
744 finalLineMask =
abs(finalModel) > self.config.footprintThreshold
746 if not finalLineMask.any():
749 fit.finalModelMax = finalModelMax
751 finalLineMasks.append(finalLineMask)
753 finalMask = np.array(finalLineMasks).
any(axis=0)
755 return lineFits, finalMask
def __init__(self, rhos, thetas, sigmas=None)
def append(self, newLine)
def __getitem__(self, index)
def __init__(self, data, weights, line=None)
def makeProfile(self, line, fitFlux=True)
def _makeMaskedProfile(self, line, fitFlux=True)
def _lineChi2(self, line, grad=True)
def fit(self, dChi2Tol=0.1, maxIter=100)
def setLineMask(self, line)
def run(self, maskedImage)
def find(self, maskedImage)
def _fitProfile(self, lines, maskedImage)
def _findClusters(self, lines)
def _cannyFilter(self, image)
bool any(CoordinateExpr< N > const &expr) noexcept
Return true if any elements are true.
bool all(CoordinateExpr< N > const &expr) noexcept
Return true if all elements are true.
def setDetectionMask(maskedImage, forceSlowBin=False, binning=None, detectedPlane="DETECTED", badMaskPlanes=("NO_DATA", "INTRP", "BAD", "SAT", "EDGE"), detectionThreshold=5)
Angle abs(Angle const &a)