LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
maskStreaks.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 # Copyright 2014 LSST Corporation.
4 #
5 # This product includes software developed by the
6 # LSST Project (http://www.lsst.org/).
7 #
8 # This program is free software: you can redistribute it and/or modify
9 # it under the terms of the GNU General Public License as published by
10 # the Free Software Foundation, either version 3 of the License, or
11 # (at your option) any later version.
12 #
13 # This program is distributed in the hope that it will be useful,
14 # but WITHOUT ANY WARRANTY; without even the implied warranty of
15 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 # GNU General Public License for more details.
17 #
18 # You should have received a copy of the LSST License Statement and
19 # the GNU General Public License along with this program. If not,
20 # see <http://www.lsstcorp.org/LegalNotices/>.
21 #
22 
23 __all__ = ["MaskStreaksConfig", "MaskStreaksTask", "setDetectionMask"]
24 
25 import lsst.pex.config as pexConfig
26 import lsst.pipe.base as pipeBase
27 import lsst.kht
28 from lsst.utils.timer import timeMethod
29 
30 import numpy as np
31 import scipy
32 import textwrap
33 import copy
34 from skimage.feature import canny
35 from sklearn.cluster import KMeans
36 import warnings
37 from dataclasses import dataclass
38 
39 
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
43 
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.
50 
51  Parameters
52  ----------
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
57  give the same result.
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
67  """
68  data = maskedImage.image.array
69  weights = 1 / maskedImage.variance.array
70  mask = maskedImage.getMask()
71 
72  detectionMask = ((mask.array & mask.getPlaneBitMask(detectedPlane)))
73  badPixelMask = mask.getPlaneBitMask(badMaskPlanes)
74  badMask = (mask.array & badPixelMask) > 0
75  fitMask = detectionMask.astype(bool) & ~badMask
76 
77  fitData = np.copy(data)
78  fitData[~fitMask] = 0
79  fitWeights = np.copy(weights)
80  fitWeights[~fitMask] = 0
81 
82  if binning:
83  # Do weighted binning:
84  ymax, xmax = fitData.shape
85  if (ymax % binning == 0) and (xmax % binning == 0) and (not forceSlowBin):
86  # Faster binning method
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)
92  else:
93  # Slower binning method when (image shape mod binsize) != 0
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',
105  bins=(ybins, xbins))
106  binnedDenominator, *_ = scipy.stats.binned_statistic_2d(ymesh.ravel(), xmesh.ravel(),
107  fitWeights.ravel(), statistic='sum',
108  bins=(ybins, xbins))
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]
116  else:
117  outputMask = (fitData * fitWeights**0.5) > detectionThreshold
118 
119  # Clear existing Detected Plane:
120  maskedImage.mask.array &= ~maskedImage.mask.getPlaneBitMask(detectedPlane)
121 
122  # Set Detected Plane with the binary detection mask:
123  maskedImage.mask.array[outputMask] |= maskedImage.mask.getPlaneBitMask(detectedPlane)
124 
125 
126 @dataclass
127 class Line:
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.
131  """
132  rho: float
133  theta: float
134  sigma: float = 0
135 
136 
138  """Collection of `Line` objects.
139 
140  Parameters
141  ----------
142  rhos : np.ndarray
143  Array of `Line` rho parameters
144  thetas : np.ndarray
145  Array of `Line` theta parameters
146  sigmas : np.ndarray (optional)
147  Array of `Line` sigma parameters
148  """
149 
150  def __init__(self, rhos, thetas, sigmas=None):
151  if sigmas is None:
152  sigmas = np.zeros(len(rhos))
153 
154  self._lines_lines = [Line(rho, theta, sigma) for (rho, theta, sigma) in
155  zip(rhos, thetas, sigmas)]
156 
157  def __len__(self):
158  return len(self._lines_lines)
159 
160  def __getitem__(self, index):
161  return self._lines_lines[index]
162 
163  def __iter__(self):
164  return iter(self._lines_lines)
165 
166  def __repr__(self):
167  joinedString = ", ".join(str(line) for line in self._lines_lines)
168  return textwrap.shorten(joinedString, width=160, placeholder="...")
169 
170  @property
171  def rhos(self):
172  return np.array([line.rho for line in self._lines_lines])
173 
174  @property
175  def thetas(self):
176  return np.array([line.theta for line in self._lines_lines])
177 
178  def append(self, newLine):
179  """Add line to current collection of lines.
180 
181  Parameters
182  ----------
183  newLine : `Line`
184  `Line` to add to current collection of lines
185  """
186  self._lines_lines.append(copy.copy(newLine))
187 
188 
190  """Construct and/or fit a model for a linear streak.
191 
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.
200 
201  Parameters
202  ----------
203  data : np.ndarray
204  2d array of data
205  weights : np.ndarray
206  2d array of weights
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
210  out.
211  """
212 
213  def __init__(self, data, weights, line=None):
214  self.datadata = data
215  self.weightsweights = weights
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)
222  self.maskmask = (weights != 0)
223 
224  self._initLine_initLine = line
225  self.setLineMasksetLineMask(line)
226 
227  def setLineMask(self, line):
228  """Set mask around the image region near the line
229 
230  Parameters
231  ----------
232  line : `Line`
233  Parameters of line in the image
234  """
235  if line:
236  # Only fit pixels within 5 sigma of the estimated line
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)
240  self.lineMasklineMask = self.maskmask & m
241  else:
242  self.lineMasklineMask = np.copy(self.maskmask)
243 
244  self.lineMaskSizelineMaskSize = self.lineMasklineMask.sum()
245  self._maskData_maskData = self.datadata[self.lineMasklineMask]
246  self._maskWeights_maskWeights = self.weightsweights[self.lineMasklineMask]
247  self._mxmesh_mxmesh = self._xmesh[self.lineMasklineMask]
248  self._mymesh_mymesh = self._ymesh_ymesh[self.lineMasklineMask]
249 
250  def _makeMaskedProfile(self, line, fitFlux=True):
251  """Construct the line model in the masked region and calculate its
252  derivatives
253 
254  Parameters
255  ----------
256  line : `Line`
257  Parameters of line profile for which to make profile in the masked
258  region
259  fitFlux : bool
260  Fit the amplitude of the line profile to the data
261 
262  Returns
263  -------
264  model : np.ndarray
265  Model in the masked region
266  dModel : np.ndarray
267  Derivative of the model in the masked region
268  """
269  invSigma = line.sigma**-1
270  # Calculate distance between pixels and line
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
276 
277  # Calculate partial derivatives of distance
278  drad = np.pi / 180
279  dDistanceSqdRho = 2 * distance * (-np.ones_like(self._mxmesh_mxmesh))
280  dDistanceSqdTheta = (2 * distance * (-sintheta * self._mxmesh_mxmesh + costheta * self._mymesh_mymesh) * drad)
281 
282  # Use pixel-line distances to make Moffat profile
283  profile = (1 + distanceSquared * invSigma**2)**-2.5
284  dProfile = -2.5 * (1 + distanceSquared * invSigma**2)**-3.5
285 
286  if fitFlux:
287  # Calculate line flux from profile and data
288  flux = ((self._maskWeights_maskWeights * self._maskData_maskData * profile).sum()
289  / (self._maskWeights_maskWeights * profile**2).sum())
290  else:
291  # Approximately normalize the line
292  flux = invSigma**-1
293  if np.isnan(flux):
294  flux = 0
295 
296  model = flux * profile
297 
298  # Calculate model derivatives
299  fluxdProfile = flux * dProfile
300  fluxdProfileInvSigma = fluxdProfile * invSigma**2
301  dModeldRho = fluxdProfileInvSigma * dDistanceSqdRho
302  dModeldTheta = fluxdProfileInvSigma * dDistanceSqdTheta
303  dModeldInvSigma = fluxdProfile * distanceSquared * 2 * invSigma
304 
305  dModel = np.array([dModeldRho, dModeldTheta, dModeldInvSigma])
306  return model, dModel
307 
308  def makeProfile(self, line, fitFlux=True):
309  """Construct the line profile model
310 
311  Parameters
312  ----------
313  line : `Line`
314  Parameters of the line profile to model
315  fitFlux : bool (optional)
316  Fit the amplitude of the line profile to the data
317 
318  Returns
319  -------
320  finalModel : np.ndarray
321  Model for line profile
322  """
323  model, _ = self._makeMaskedProfile_makeMaskedProfile(line, fitFlux=fitFlux)
324  finalModel = np.zeros((self._ymax, self._xmax_xmax), dtype=self._dtype_dtype)
325  finalModel[self.lineMasklineMask] = model
326  return finalModel
327 
328  def _lineChi2(self, line, grad=True):
329  """Construct the chi2 between the data and the model
330 
331  Parameters
332  ----------
333  line : `Line`
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
337 
338  Returns
339  -------
340  reducedChi : float
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
346  """
347  # Calculate chi2
348  model, dModel = self._makeMaskedProfile_makeMaskedProfile(line)
349  chi2 = (self._maskWeights_maskWeights * (self._maskData_maskData - model)**2).sum()
350  if not grad:
351  return chi2.sum() / self.lineMaskSizelineMaskSize
352 
353  # Calculate derivative and Hessian of chi2
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)
356 
357  reducedChi = chi2 / self.lineMaskSizelineMaskSize
358  reducedDChi = derivChi2 / self.lineMaskSizelineMaskSize
359  reducedHessianChi = hessianChi2 / self.lineMaskSizelineMaskSize
360  return reducedChi, reducedDChi, reducedHessianChi
361 
362  def fit(self, dChi2Tol=0.1, maxIter=100):
363  """Perform Newton-Raphson minimization to find line parameters
364 
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
371  needed.
372 
373  Parameters
374  ----------
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.
381 
382  Returns
383  -------
384  outline : np.ndarray
385  Coordinates and inverse width of fit line
386  chi2 : float
387  Reduced Chi2 of model fit to data
388  fitFailure : bool
389  Boolean where `False` corresponds to a successful fit
390  """
391  # Do minimization on inverse of sigma to simplify derivatives:
392  x = np.array([self._initLine_initLine.rho, self._initLine_initLine.theta, self._initLine_initLine.sigma**-1])
393 
394  dChi2 = 1
395  iter = 0
396  oldChi2 = 0
397  fitFailure = False
398 
399  def line_search(c, dx):
400  testx = x - c * dx
401  testLine = Line(testx[0], testx[1], testx[2]**-1)
402  return self._lineChi2_lineChi2(testLine, grad=False)
403 
404  while abs(dChi2) > dChi2Tol:
405  line = Line(x[0], x[1], x[2]**-1)
406  chi2, b, A = self._lineChi2_lineChi2(line)
407  if chi2 == 0:
408  break
409  if not np.isfinite(A).all():
410  # TODO: DM-30797 Add warning here.
411  fitFailure = True
412  break
413  dChi2 = oldChi2 - chi2
414  cholesky = scipy.linalg.cho_factor(A)
415  dx = scipy.linalg.cho_solve(cholesky, b)
416 
417  factor, fmin, _, _ = scipy.optimize.brent(line_search, args=(dx,), full_output=True, tol=0.05)
418  x -= factor * dx
419  if (abs(x[0]) > 1.5 * self._rhoMax_rhoMax) or (iter > maxIter):
420  fitFailure = True
421  break
422  oldChi2 = chi2
423  iter += 1
424 
425  outline = Line(x[0], x[1], abs(x[2])**-1)
426 
427  return outline, chi2, fitFailure
428 
429 
430 class MaskStreaksConfig(pexConfig.Config):
431  """Configuration parameters for `MaskStreaksTask`
432  """
433  minimumKernelHeight = pexConfig.Field(
434  doc="Minimum height of the streak-finding kernel relative to the tallest kernel",
435  dtype=float,
436  default=0.0,
437  )
438  absMinimumKernelHeight = pexConfig.Field(
439  doc="Minimum absolute height of the streak-finding kernel",
440  dtype=float,
441  default=5,
442  )
443  clusterMinimumSize = pexConfig.Field(
444  doc="Minimum size in pixels of detected clusters",
445  dtype=int,
446  default=50,
447  )
448  clusterMinimumDeviation = pexConfig.Field(
449  doc="Allowed deviation (in pixels) from a straight line for a detected "
450  "line",
451  dtype=int,
452  default=2,
453  )
454  delta = pexConfig.Field(
455  doc="Stepsize in angle-radius parameter space",
456  dtype=float,
457  default=0.2,
458  )
459  nSigma = pexConfig.Field(
460  doc="Number of sigmas from center of kernel to include in voting "
461  "procedure",
462  dtype=float,
463  default=2,
464  )
465  rhoBinSize = pexConfig.Field(
466  doc="Binsize in pixels for position parameter rho when finding "
467  "clusters of detected lines",
468  dtype=float,
469  default=30,
470  )
471  thetaBinSize = pexConfig.Field(
472  doc="Binsize in degrees for angle parameter theta when finding "
473  "clusters of detected lines",
474  dtype=float,
475  default=2,
476  )
477  invSigma = pexConfig.Field(
478  doc="Inverse of the Moffat sigma parameter (in units of pixels)"
479  "describing the profile of the streak",
480  dtype=float,
481  default=10.**-1,
482  )
483  footprintThreshold = pexConfig.Field(
484  doc="Threshold at which to determine edge of line, in units of the line"
485  "profile maximum",
486  dtype=float,
487  default=0.01
488  )
489  dChi2Tolerance = pexConfig.Field(
490  doc="Absolute difference in Chi2 between iterations of line profile"
491  "fitting that is acceptable for convergence",
492  dtype=float,
493  default=0.1
494  )
495  detectedMaskPlane = pexConfig.Field(
496  doc="Name of mask with pixels above detection threshold, used for first"
497  "estimate of streak locations",
498  dtype=str,
499  default="DETECTED"
500  )
501  streaksMaskPlane = pexConfig.Field(
502  doc="Name of mask plane holding detected streaks",
503  dtype=str,
504  default="STREAK"
505  )
506 
507 
508 class MaskStreaksTask(pipeBase.Task):
509  """Find streaks or other straight lines in image data.
510 
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
523  (non-binary) image.
524  """
525 
526  ConfigClass = MaskStreaksConfig
527  _DefaultName = "maskStreaks"
528 
529  @timeMethod
530  def find(self, maskedImage):
531  """Find streaks in a masked image
532 
533  Parameters
534  ----------
535  maskedImage : `lsst.afw.image.maskedImage`
536  The image in which to search for streaks.
537 
538  Returns
539  -------
540  result : `lsst.pipe.base.Struct`
541  Result struct with components:
542 
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
547  """
548  mask = maskedImage.getMask()
549  detectionMask = (mask.array & mask.getPlaneBitMask(self.config.detectedMaskPlane))
550 
551  self.edgesedges = self._cannyFilter_cannyFilter(detectionMask)
552  self.lineslines = self._runKHT_runKHT(self.edgesedges)
553 
554  if len(self.lineslines) == 0:
555  lineMask = np.zeros(detectionMask.shape, dtype=bool)
556  fitLines = LineCollection([], [])
557  clusters = LineCollection([], [])
558  else:
559  clusters = self._findClusters_findClusters(self.lineslines)
560  fitLines, lineMask = self._fitProfile_fitProfile(clusters, maskedImage)
561 
562  # The output mask is the intersection of the fit streaks and the image detections
563  outputMask = lineMask & detectionMask.astype(bool)
564 
565  return pipeBase.Struct(
566  lines=fitLines,
567  lineClusters=clusters,
568  originalLines=self.lineslines,
569  mask=outputMask,
570  )
571 
572  @timeMethod
573  def run(self, maskedImage):
574  """Find and mask streaks in a masked image.
575 
576  Finds streaks in the image and modifies maskedImage in place by adding a
577  mask plane with any identified streaks.
578 
579  Parameters
580  ----------
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
584  detected pixels.
585 
586  Returns
587  -------
588  result : `lsst.pipe.base.Struct`
589  Result struct with components:
590 
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
594  """
595  streaks = self.findfind(maskedImage)
596 
597  maskedImage.mask.addMaskPlane(self.config.streaksMaskPlane)
598  maskedImage.mask.array[streaks.mask] |= maskedImage.mask.getPlaneBitMask(self.config.streaksMaskPlane)
599 
600  return pipeBase.Struct(
601  lines=streaks.lines,
602  lineClusters=streaks.lineClusters,
603  originalLines=streaks.originalLines,
604  )
605 
606  def _cannyFilter(self, image):
607  """Apply a canny filter to the data in order to detect edges
608 
609  Parameters
610  ----------
611  image : `np.ndarray`
612  2-d image data on which to run filter
613 
614  Returns
615  -------
616  cannyData : `np.ndarray`
617  2-d image of edges found in input image
618  """
619  filterData = image.astype(int)
620  return canny(filterData, low_threshold=0, high_threshold=1, sigma=0.1)
621 
622  def _runKHT(self, image):
623  """Run Kernel Hough Transform on image.
624 
625  Parameters
626  ----------
627  image : `np.ndarray`
628  2-d image data on which to detect lines
629 
630  Returns
631  -------
632  result : `LineCollection`
633  Collection of detected lines, with their detected rho and theta
634  coordinates.
635  """
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)
640 
641  return LineCollection(lines.rho, lines.theta)
642 
643  def _findClusters(self, lines):
644  """Group lines that are close in parameter space and likely describe
645  the same streak.
646 
647  Parameters
648  ----------
649  lines : `LineCollection`
650  Collection of lines to group into clusters
651 
652  Returns
653  -------
654  result : `LineCollection`
655  Average `Line` for each cluster of `Line`s in the input
656  `LineCollection`
657  """
658  # Scale variables by threshold bin-size variable so that rho and theta
659  # are on the same scale. Since the clustering algorithm below stops when
660  # the standard deviation <= 1, after rescaling each cluster will have a
661  # standard deviation at or below the bin-size.
662  x = lines.rhos / self.config.rhoBinSize
663  y = lines.thetas / self.config.thetaBinSize
664  X = np.array([x, y]).T
665  nClusters = 1
666 
667  # Put line parameters in clusters by starting with all in one, then
668  # subdividing until the parameters of each cluster have std dev=1.
669  # If nClusters == len(lines), each line will have its own 'cluster', so
670  # the standard deviations of each cluster must be zero and the loop
671  # is guaranteed to stop.
672  while True:
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)
678  # Are the rhos and thetas in each cluster all below the threshold?
679  if (clusterStandardDeviations <= 1).all():
680  break
681  nClusters += 1
682 
683  # The cluster centers are final line estimates
684  finalClusters = kmeans.cluster_centers_.T
685 
686  # Rescale variables:
687  finalRhos = finalClusters[0] * self.config.rhoBinSize
688  finalThetas = finalClusters[1] * self.config.thetaBinSize
689  result = LineCollection(finalRhos, finalThetas)
690 
691  return result
692 
693  def _fitProfile(self, lines, maskedImage):
694  """Fit the profile of the streak.
695 
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.
699 
700  Parameters
701  ----------
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.
706 
707  Returns
708  -------
709  lineFits : `LineCollection`
710  Collection of `Line` profiles fit to the data
711  finalMask : `np.ndarray`
712  2d mask array with detected streaks=1.
713  """
714  data = maskedImage.image.array
715  weights = maskedImage.variance.array**-1
716  # Mask out any pixels with non-finite weights
717  weights[~np.isfinite(weights) | ~np.isfinite(data)] = 0
718 
719  lineFits = LineCollection([], [])
720  finalLineMasks = [np.zeros(data.shape, dtype=bool)]
721  for line in lines:
722  line.sigma = self.config.invSigma**-1
723  lineModel = LineProfile(data, weights, line=line)
724  # Skip any lines that do not cover any data (sometimes happens because of chip gaps)
725  if lineModel.lineMaskSize == 0:
726  continue
727 
728  fit, chi2, fitFailure = lineModel.fit(dChi2Tol=self.config.dChi2Tolerance)
729 
730  # Initial estimate should be quite close: fit is deemed unsuccessful if rho or theta
731  # change more than the allowed bin in rho or theta:
732  if ((abs(fit.rho - line.rho) > 2 * self.config.rhoBinSize)
733  or (abs(fit.theta - line.theta) > 2 * self.config.thetaBinSize)):
734  fitFailure = True
735 
736  if fitFailure:
737  continue
738 
739  # Make mask
740  lineModel.setLineMask(fit)
741  finalModel = lineModel.makeProfile(fit)
742  # Take absolute value, as streaks are allowed to be negative
743  finalModelMax = abs(finalModel).max()
744  finalLineMask = abs(finalModel) > self.config.footprintThreshold
745  # Drop this line if the model profile is below the footprint threshold
746  if not finalLineMask.any():
747  continue
748  fit.chi2 = chi2
749  fit.finalModelMax = finalModelMax
750  lineFits.append(fit)
751  finalLineMasks.append(finalLineMask)
752 
753  finalMask = np.array(finalLineMasks).any(axis=0)
754 
755  return lineFits, finalMask
int max
def __init__(self, rhos, thetas, sigmas=None)
Definition: maskStreaks.py:150
def __init__(self, data, weights, line=None)
Definition: maskStreaks.py:213
def makeProfile(self, line, fitFlux=True)
Definition: maskStreaks.py:308
def _makeMaskedProfile(self, line, fitFlux=True)
Definition: maskStreaks.py:250
def _lineChi2(self, line, grad=True)
Definition: maskStreaks.py:328
def fit(self, dChi2Tol=0.1, maxIter=100)
Definition: maskStreaks.py:362
def _fitProfile(self, lines, maskedImage)
Definition: maskStreaks.py:693
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)
Definition: maskStreaks.py:41
Angle abs(Angle const &a)
Definition: Angle.h:106