Coverage for python / lsst / meas / algorithms / maskStreaks.py: 18%

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 

56class LineCollection: 

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 

112class LineProfile: 

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 detectionMask : `np.ndarray`, optional 

135 2-d boolean array where detected pixels are True. 

136 """ 

137 

138 def __init__(self, data, weights, line=None, detectionMask=None): 

139 self.data = data 

140 self.weights = weights 

141 self._ymax, self._xmax = data.shape 

142 self._dtype = data.dtype 

143 xrange = np.arange(self._xmax) - self._xmax / 2. 

144 yrange = np.arange(self._ymax) - self._ymax / 2. 

145 self._rhoMax = ((0.5 * self._ymax)**2 + (0.5 * self._xmax)**2)**0.5 

146 self._xmesh, self._ymesh = np.meshgrid(xrange, yrange) 

147 self.mask = (weights != 0) 

148 

149 self._initLine = line 

150 self.modelFailure = False 

151 self.setLineMask(line, maxStreakWidth=0, nSigmaMask=10, detectionMask=detectionMask) 

152 

153 def getLineXY(self, line): 

154 """Return the pixel coordinates of the ends of the line. 

155 

156 Parameters 

157 ---------- 

158 line : `Line` 

159 Line for which to find the endpoints. 

160 

161 Returns 

162 ------- 

163 boxIntersections : `np.ndarray` 

164 (x, y) coordinates of the start and endpoints of the line. 

165 """ 

166 theta = line.theta * u.deg 

167 # Determine where the line intersects with each edge of the bounding 

168 # box. 

169 # Bottom: 

170 yA = -self._ymax / 2. 

171 xA = (line.rho - yA * np.sin(theta)) / np.cos(theta) 

172 # Left: 

173 xB = -self._xmax / 2. 

174 yB = (line.rho - xB * np.cos(theta)) / np.sin(theta) 

175 # Top: 

176 yC = self._ymax / 2. 

177 xC = (line.rho - yC * np.sin(theta)) / np.cos(theta) 

178 # Right: 

179 xD = self._xmax / 2. 

180 yD = (line.rho - xD * np.cos(theta)) / np.sin(theta) 

181 

182 lineIntersections = np.array([[xA, yA], 

183 [xB, yB], 

184 [xC, yC], 

185 [xD, yD]]) 

186 lineIntersections[:, 0] += self._xmax / 2. 

187 lineIntersections[:, 1] += self._ymax / 2. 

188 

189 # The line will necessarily intersect with exactly two edges of the 

190 # bounding box itself. 

191 inBox = ((lineIntersections[:, 0] >= 0) & (lineIntersections[:, 0] <= self._xmax) 

192 & (lineIntersections[:, 1] >= 0) & (lineIntersections[:, 1] <= self._ymax)) 

193 boxIntersections = lineIntersections[inBox] 

194 

195 return boxIntersections 

196 

197 def setLineMask(self, line, maxStreakWidth, nSigmaMask, logger=None, detectionMask=None): 

198 """Set mask around the image region near the line. 

199 

200 Parameters 

201 ---------- 

202 line : `Line` 

203 Parameters of line in the image. 

204 maxStreakWidth : `float` 

205 Maximum width in pixels of streak mask. 

206 nSigmaMask : `float` 

207 Factor by which to multiply the line's width to get the mask width. 

208 logger : `lsst.utils.logging.LsstLogAdapter`, optional 

209 Logger to use for reporting when maxStreakWidth is reached. 

210 detectionMask : `np.ndarray`, optional 

211 2-d boolean array where detected pixels are True. 

212 """ 

213 if line: 

214 # Only fit pixels within nSigmaMask * sigma of the estimated line 

215 radtheta = np.deg2rad(line.theta) 

216 distance = (np.cos(radtheta) * self._xmesh + np.sin(radtheta) * self._ymesh - line.rho) 

217 

218 width = 2 * nSigmaMask * line.sigma 

219 if (maxStreakWidth > 0) and (maxStreakWidth < width): 

220 if logger is not None: 

221 logger.info("Line with width %d exceeded maximum configured width of %d", 

222 width, maxStreakWidth) 

223 width = maxStreakWidth 

224 m = (abs(distance) < width/2) 

225 self.lineMask = self.mask & m 

226 if detectionMask is not None: 

227 # Mask out areas where there are no detected pixels. This 

228 # happens when, for example, the streak ends in the middle of 

229 # the image. 

230 lineEnds = self.getLineXY(line) 

231 if lineEnds.size == 0: 

232 if logger is not None: 

233 logger.debug("Calculated line not contained in image bounding box") 

234 self.modelFailure = True 

235 return 

236 xA = lineEnds[0, 0] - self._xmax / 2. 

237 yA = lineEnds[0, 1] - self._ymax / 2. 

238 

239 radtheta = np.deg2rad(line.theta) 

240 costheta = np.cos(radtheta) 

241 sintheta = np.sin(radtheta) 

242 

243 maskDetections = detectionMask[self.lineMask] != 0 

244 distanceFromLineEnd = (- sintheta * self._xmesh[self.lineMask] 

245 + costheta * self._ymesh[self.lineMask] 

246 + sintheta * xA 

247 - costheta * yA) 

248 lineBins = np.arange(distanceFromLineEnd.min(), distanceFromLineEnd.max() + 5.1, 5) 

249 # Get the chi2 of the pixels perpendicular to the streak: 

250 detectionsAlongStreak, _, binnumber = scipy.stats.binned_statistic(distanceFromLineEnd, 

251 maskDetections, 

252 statistic='sum', 

253 bins=lineBins) 

254 countAlongStreak, *_ = scipy.stats.binned_statistic(distanceFromLineEnd, maskDetections, 

255 statistic='count', bins=lineBins) 

256 detectionFraction = detectionsAlongStreak / countAlongStreak 

257 emptyRows = detectionFraction < (0.5 * np.median(detectionFraction[detectionFraction != 0])) 

258 emptyDetections = emptyRows[binnumber - 1] 

259 self.lineMask[self.lineMask] = ~emptyDetections 

260 else: 

261 self.lineMask = np.copy(self.mask) 

262 

263 self._maskData = self.data[self.lineMask] 

264 self._maskWeights = self.weights[self.lineMask] 

265 self._mxmesh = self._xmesh[self.lineMask] 

266 self._mymesh = self._ymesh[self.lineMask] 

267 

268 def _makeMaskedProfile(self, line, fitFlux=True): 

269 """Construct the line model in the masked region and calculate its 

270 derivatives. 

271 

272 Parameters 

273 ---------- 

274 line : `Line` 

275 Parameters of line profile for which to make profile in the masked 

276 region. 

277 fitFlux : `bool` 

278 Fit the amplitude of the line profile to the data. 

279 

280 Returns 

281 ------- 

282 model : `np.ndarray` 

283 Model in the masked region. 

284 dModel : `np.ndarray` 

285 Derivative of the model in the masked region. 

286 """ 

287 invSigma = line.sigma**-1 

288 # Calculate distance between pixels and line 

289 radtheta = np.deg2rad(line.theta) 

290 costheta = np.cos(radtheta) 

291 sintheta = np.sin(radtheta) 

292 distance = (costheta * self._mxmesh + sintheta * self._mymesh - line.rho) 

293 distanceSquared = distance**2 

294 

295 # Calculate partial derivatives of distance 

296 drad = np.pi / 180 

297 dDistanceSqdRho = 2 * distance * (-np.ones_like(self._mxmesh)) 

298 dDistanceSqdTheta = (2 * distance * (-sintheta * self._mxmesh + costheta * self._mymesh) * drad) 

299 

300 # Use pixel-line distances to make Moffat profile 

301 profile = (1 + distanceSquared * invSigma**2)**-2.5 

302 dProfile = -2.5 * (1 + distanceSquared * invSigma**2)**-3.5 

303 

304 if fitFlux: 

305 # Calculate line flux from profile and data 

306 flux = ((self._maskWeights * self._maskData * profile).sum() 

307 / (self._maskWeights * profile**2).sum()) 

308 else: 

309 # Approximately normalize the line 

310 flux = invSigma**-1 

311 if np.isnan(flux): 

312 flux = 0 

313 

314 model = flux * profile 

315 

316 # Calculate model derivatives 

317 fluxdProfile = flux * dProfile 

318 fluxdProfileInvSigma = fluxdProfile * invSigma**2 

319 dModeldRho = fluxdProfileInvSigma * dDistanceSqdRho 

320 dModeldTheta = fluxdProfileInvSigma * dDistanceSqdTheta 

321 dModeldInvSigma = fluxdProfile * distanceSquared * 2 * invSigma 

322 

323 dModel = np.array([dModeldRho, dModeldTheta, dModeldInvSigma]) 

324 return model, dModel 

325 

326 def makeProfile(self, line, fitFlux=True): 

327 """Construct the line profile model. 

328 

329 Parameters 

330 ---------- 

331 line : `Line` 

332 Parameters of the line profile to model. 

333 fitFlux : `bool`, optional 

334 Fit the amplitude of the line profile to the data. 

335 

336 Returns 

337 ------- 

338 finalModel : `np.ndarray` 

339 Model for line profile. 

340 """ 

341 model, _ = self._makeMaskedProfile(line, fitFlux=fitFlux) 

342 finalModel = np.zeros((self._ymax, self._xmax), dtype=self._dtype) 

343 finalModel[self.lineMask] = model 

344 return finalModel 

345 

346 def _lineChi2(self, line, grad=True): 

347 """Construct the chi2 between the data and the model. 

348 

349 Parameters 

350 ---------- 

351 line : `Line` 

352 `Line` parameters for which to build model and calculate chi2. 

353 grad : `bool`, optional 

354 Whether or not to return the gradient and hessian. 

355 

356 Returns 

357 ------- 

358 reducedChi : `float` 

359 Reduced chi2 of the model. 

360 reducedDChi : `np.ndarray` 

361 Derivative of the chi2 with respect to rho, theta, invSigma. 

362 reducedHessianChi : `np.ndarray` 

363 Hessian of the chi2 with respect to rho, theta, invSigma. 

364 """ 

365 # Calculate chi2 

366 model, dModel = self._makeMaskedProfile(line) 

367 chi2 = (self._maskWeights * (self._maskData - model)**2).sum() 

368 maskSize = (self._maskWeights != 0).sum() 

369 if not grad: 

370 return chi2.sum() / maskSize 

371 

372 # Calculate derivative and Hessian of chi2 

373 derivChi2 = ((-2 * self._maskWeights * (self._maskData - model))[None, :] * dModel).sum(axis=1) 

374 hessianChi2 = (2 * self._maskWeights * dModel[:, None, :] * dModel[None, :, :]).sum(axis=2) 

375 

376 reducedChi = chi2 / maskSize 

377 reducedDChi = derivChi2 / maskSize 

378 reducedHessianChi = hessianChi2 / maskSize 

379 return reducedChi, reducedDChi, reducedHessianChi 

380 

381 def _rejectOutliers(self, line): 

382 """Reject outlier pixels. 

383 

384 This calculates the chi2/dof in bins of pixels perpendicular to the 

385 streak direction and removes outliers. This is done so that the profile 

386 fitter ignores regions like the area around bright stars. 

387 

388 Parameters 

389 ---------- 

390 line : `Line` 

391 `Line` parameters for which to build model and calculate chi2. 

392 

393 Returns 

394 ------- 

395 nOutliers : `int` 

396 Number of outlier pixels. 

397 """ 

398 model, _ = self._makeMaskedProfile(line) 

399 pixelChi2 = (self._maskWeights * (self._maskData - model)**2) 

400 

401 lineEnds = self.getLineXY(line) 

402 

403 if lineEnds.shape == (2, 2): 

404 xA = lineEnds[0, 0] - self._xmax / 2. 

405 yA = lineEnds[0, 1] - self._ymax / 2. 

406 else: 

407 # Line profile is outside the detector bounding box. Exit outlier rejection. 

408 return 0 

409 

410 radtheta = np.deg2rad(line.theta) 

411 costheta = np.cos(radtheta) 

412 sintheta = np.sin(radtheta) 

413 

414 distanceFromLineEnd = (- sintheta * self._mxmesh + costheta * self._mymesh + sintheta * xA 

415 - costheta * yA) 

416 

417 distanceFromLineEnd = distanceFromLineEnd[self._maskWeights != 0] 

418 nonZeroPixelChi2 = pixelChi2[self._maskWeights != 0] 

419 lineBins = np.arange(distanceFromLineEnd.min(), distanceFromLineEnd.max() + 5.1, 5) 

420 # Get the chi2 of the pixels perpendicular to the streak: 

421 chi2AlongStreak, _, binnumber = scipy.stats.binned_statistic(distanceFromLineEnd, nonZeroPixelChi2, 

422 statistic='sum', bins=lineBins) 

423 countAlongStreak, *_ = scipy.stats.binned_statistic(distanceFromLineEnd, nonZeroPixelChi2, 

424 statistic='count', bins=lineBins) 

425 

426 rChi2AlongStreak = chi2AlongStreak / countAlongStreak 

427 outliers = rChi2AlongStreak > (np.nanmean(rChi2AlongStreak[rChi2AlongStreak != 0]) 

428 + 3 * np.nanstd(rChi2AlongStreak[rChi2AlongStreak != 0])) 

429 

430 outlierPix = outliers[binnumber - 1] 

431 tmpWeights = self._maskWeights[self._maskWeights != 0] 

432 tmpWeights[outlierPix] = 0 

433 self._maskWeights[self._maskWeights != 0] = tmpWeights 

434 

435 nOutliers = outlierPix.sum() 

436 

437 return nOutliers 

438 

439 def fit(self, dChi2Tol=0.1, maxIter=100, log=None): 

440 """Perform Newton-Raphson minimization to find line parameters. 

441 

442 This method takes advantage of having known derivative and Hessian of 

443 the multivariate function to quickly and efficiently find the minimum. 

444 This is more efficient than the scipy implementation of the Newton- 

445 Raphson method, which doesn't take advantage of the Hessian matrix. The 

446 method here also performs a line search in the direction of the steepest 

447 derivative at each iteration, which reduces the number of iterations 

448 needed. 

449 

450 Parameters 

451 ---------- 

452 dChi2Tol : `float`, optional 

453 Change in Chi2 tolerated for fit convergence. 

454 maxIter : `int`, optional 

455 Maximum number of fit iterations allowed. The fit should converge in 

456 ~10 iterations, depending on the value of dChi2Tol, but this 

457 maximum provides a backup. 

458 log : `lsst.utils.logging.LsstLogAdapter`, optional 

459 Logger to use for reporting more details for fitting failures. 

460 

461 Returns 

462 ------- 

463 outline : `Line` 

464 Coordinates, inverse width, and chi2 of fit line. 

465 fitFailure : `bool` 

466 Boolean where `False` corresponds to a successful fit. 

467 """ 

468 # Do minimization on inverse of sigma to simplify derivatives: 

469 x = np.array([self._initLine.rho, self._initLine.theta, self._initLine.sigma**-1]) 

470 

471 dChi2 = 1 

472 iter = 0 

473 oldChi2 = 0 

474 nOutliers = 1 

475 fitFailure = False 

476 

477 def line_search(c, dx): 

478 testx = x - c * dx 

479 testLine = Line(testx[0], testx[1], testx[2]**-1) 

480 return self._lineChi2(testLine, grad=False) 

481 

482 while (abs(dChi2) > dChi2Tol) or (nOutliers != 0): 

483 line = Line(x[0], x[1], x[2]**-1) 

484 chi2, b, A = self._lineChi2(line) 

485 if chi2 == 0: 

486 break 

487 if not np.isfinite(A).all(): 

488 fitFailure = True 

489 if log is not None: 

490 log.warning("Hessian matrix has non-finite elements.") 

491 break 

492 dChi2 = oldChi2 - chi2 

493 try: 

494 cholesky = scipy.linalg.cho_factor(A) 

495 except np.linalg.LinAlgError: 

496 fitFailure = True 

497 if log is not None: 

498 log.warning("Hessian matrix is not invertible.") 

499 break 

500 dx = scipy.linalg.cho_solve(cholesky, b) 

501 

502 if abs(line_search(1, dx) - chi2) < 1e-12: 

503 # Step size is too small for the brent line search to work well. 

504 # Just use the step fit from the Cholesky solve. 

505 factor = 1 

506 else: 

507 factor, _, _, _ = scipy.optimize.brent(line_search, args=(dx,), full_output=True, tol=0.05) 

508 

509 x -= factor * dx 

510 if (abs(x[0]) > 1.5 * self._rhoMax) or (iter > maxIter): 

511 fitFailure = True 

512 break 

513 oldChi2 = chi2 

514 

515 nOutliers = self._rejectOutliers(line) 

516 iter += 1 

517 

518 outline = Line(x[0], x[1], abs(x[2])**-1, chi2) 

519 return outline, fitFailure 

520 

521 

522class MaskStreaksConfig(pexConfig.Config): 

523 """Configuration parameters for `MaskStreaksTask`. 

524 """ 

525 

526 minimumKernelHeight = pexConfig.Field( 

527 doc="Minimum height of the streak-finding kernel relative to the tallest kernel", 

528 dtype=float, 

529 default=0.0, 

530 ) 

531 absMinimumKernelHeight = pexConfig.Field( 

532 doc="Minimum absolute height of the streak-finding kernel", 

533 dtype=float, 

534 default=5, 

535 ) 

536 clusterMinimumSize = pexConfig.Field( 

537 doc="Minimum size in pixels of detected clusters", 

538 dtype=int, 

539 default=50, 

540 ) 

541 clusterMinimumDeviation = pexConfig.Field( 

542 doc="Allowed deviation (in pixels) from a straight line for a detected " 

543 "line", 

544 dtype=int, 

545 default=2, 

546 ) 

547 delta = pexConfig.Field( 

548 doc="Stepsize in angle-radius parameter space", 

549 dtype=float, 

550 default=0.2, 

551 ) 

552 nSigma = pexConfig.Field( 

553 doc="Number of sigmas from center of kernel to include in voting " 

554 "procedure", 

555 dtype=float, 

556 default=2, 

557 ) 

558 nSigmaMask = pexConfig.Field( 

559 doc="Number of sigmas from center of kernel to mask", 

560 dtype=float, 

561 default=5, 

562 ) 

563 rhoBinSize = pexConfig.Field( 

564 doc="Binsize in pixels for position parameter rho when finding " 

565 "clusters of detected lines", 

566 dtype=float, 

567 default=40, 

568 ) 

569 thetaBinSize = pexConfig.Field( 

570 doc="Binsize in degrees for angle parameter theta when finding " 

571 "clusters of detected lines", 

572 dtype=float, 

573 default=2, 

574 ) 

575 invSigma = pexConfig.Field( 

576 doc="Inverse of the Moffat sigma parameter (in units of pixels)" 

577 "describing the profile of the streak", 

578 dtype=float, 

579 default=10.**-1, 

580 ) 

581 footprintThreshold = pexConfig.Field( 

582 doc="Threshold at which to determine edge of line, in units of " 

583 "nanoJanskys", 

584 dtype=float, 

585 default=0.01 

586 ) 

587 dChi2Tolerance = pexConfig.Field( 

588 doc="Absolute difference in Chi2 between iterations of line profile" 

589 "fitting that is acceptable for convergence", 

590 dtype=float, 

591 default=0.1 

592 ) 

593 maxFitIter = pexConfig.Field( 

594 doc="Maximum number of line profile fitting iterations that is " 

595 "acceptable for convergence", 

596 dtype=int, 

597 default=100 

598 ) 

599 detectedMaskPlane = pexConfig.Field( 

600 doc="Name of mask with pixels above detection threshold, used for first" 

601 "estimate of streak locations", 

602 dtype=str, 

603 default="DETECTED" 

604 ) 

605 onlyMaskDetected = pexConfig.Field( 

606 doc=("If true, only propagate the part of the streak mask that " 

607 "overlaps with the detection mask."), 

608 dtype=bool, 

609 default=True, 

610 ) 

611 streaksMaskPlane = pexConfig.Field( 

612 doc="Name of mask plane holding detected streaks", 

613 dtype=str, 

614 default="STREAK" 

615 ) 

616 badMaskPlanes = pexConfig.ListField( 

617 doc=("Names of mask plane regions to ignore entirely when doing streak" 

618 " detection"), 

619 dtype=str, 

620 default=("NO_DATA", "INTRP", "BAD", "SAT", "EDGE"), 

621 ) 

622 maxStreakWidth = pexConfig.Field( 

623 doc="Maximum width in pixels to allow for masking a streak." 

624 "The fit streak parameters will not be modified, and a warning will " 

625 "be issued if the fitted width is larger than this value." 

626 "Set to 0 to disable.", 

627 dtype=float, 

628 default=0., 

629 ) 

630 saturatedDetectionsDilation = pexConfig.Field( 

631 doc="Mask out the region around saturated detections by dilating the " 

632 "existing mask by this number of pixels.", 

633 dtype=int, 

634 default=250, 

635 ) 

636 

637 

638class MaskStreaksTask(pipeBase.Task): 

639 """Find streaks or other straight lines in image data. 

640 

641 Nearby objects passing through the field of view of the telescope leave a 

642 bright trail in images. This class uses the Kernel Hough Transform (KHT) 

643 (Fernandes and Oliveira, 2007), implemented in `lsst.houghtransform`. The 

644 procedure works by taking a binary image, either provided as put or produced 

645 from the input data image, using a Canny filter to make an image of the 

646 edges in the original image, then running the KHT on the edge image. The KHT 

647 identifies clusters of non-zero points, breaks those clusters of points into 

648 straight lines, keeps clusters with a size greater than the user-set 

649 threshold, then performs a voting procedure to find the best-fit coordinates 

650 of any straight lines. Given the results of the KHT algorithm, clusters of 

651 lines are identified and grouped (generally these correspond to the two 

652 edges of a strea) and a profile is fit to the streak in the original 

653 (non-binary) image. 

654 """ 

655 

656 ConfigClass = MaskStreaksConfig 

657 _DefaultName = "maskStreaks" 

658 

659 @timeMethod 

660 def find(self, maskedImage): 

661 """Find streaks in a masked image. 

662 

663 Parameters 

664 ---------- 

665 maskedImage : `lsst.afw.image.maskedImage` 

666 The image in which to search for streaks. 

667 

668 Returns 

669 ------- 

670 result : `lsst.pipe.base.Struct` 

671 Results as a struct with attributes: 

672 

673 ``originalLines`` 

674 Lines identified by kernel hough transform. 

675 ``lineClusters`` 

676 Lines grouped into clusters in rho-theta space. 

677 ``lines`` 

678 Final result for lines after line-profile fit. 

679 ``mask`` 

680 2-d boolean mask where detected lines are True. 

681 """ 

682 mask = maskedImage.mask 

683 detectionMask = (mask.array & mask.getPlaneBitMask(self.config.detectedMaskPlane)) 

684 

685 initEdges = self._cannyFilter(detectionMask) 

686 # Ignore regions with known bad masks, adding a one-pixel buffer around 

687 # each to ensure that the edges of bad regions are also ignored. 

688 ignoreMask = mask.clone() 

689 

690 badPixelMask = mask.getPlaneBitMask(self.config.badMaskPlanes) 

691 badMaskSpanSet = SpanSet.fromMask(mask, badPixelMask).split() 

692 for sset in badMaskSpanSet: 

693 sset_dilated = sset.dilated(1) 

694 sset_dilated.clippedTo( 

695 ignoreMask.getBBox()).setMask(ignoreMask, ignoreMask.getPlaneBitMask("BAD")) 

696 

697 # TODO: DM-52769, replace this with a model for the diffraction spikes 

698 # around bright stars once DM-52541 is done. 

699 if self.config.saturatedDetectionsDilation: 

700 # Dilate spansets that are both detected and saturated mask by a lot more: 

701 satMask = mask.getPlaneBitMask("SAT") 

702 satMask = (mask.array & mask.getPlaneBitMask("SAT")) 

703 satDetMask = (satMask != 0) & (detectionMask != 0) 

704 satDetIm = lsst.afw.image.Mask(satDetMask.astype(np.int32)) 

705 satSpanSet = SpanSet.fromMask(satDetIm, 1).split() 

706 for sset in satSpanSet: 

707 sset_dilated = sset.dilated(self.config.saturatedDetectionsDilation) 

708 sset_dilated.clippedTo( 

709 ignoreMask.getBBox()).setMask(ignoreMask, ignoreMask.getPlaneBitMask("BAD")) 

710 

711 dilatedBadMask = (ignoreMask.array & badPixelMask) > 0 

712 self.edges = initEdges & ~dilatedBadMask 

713 self.lines = self._runKHT(self.edges) 

714 

715 if len(self.lines) == 0: 

716 lineMask = np.zeros(detectionMask.shape, dtype=bool) 

717 fitLines = LineCollection([], []) 

718 clusters = LineCollection([], []) 

719 else: 

720 clusters = self._findClusters(self.lines) 

721 fitLines, lineMask = self._fitProfile(clusters, maskedImage, detectionMask=detectionMask) 

722 

723 if self.config.onlyMaskDetected: 

724 # The output mask is the intersection of the fit streaks and the image detections 

725 lineMask &= detectionMask.astype(bool) 

726 

727 return pipeBase.Struct( 

728 lines=fitLines, 

729 lineClusters=clusters, 

730 originalLines=self.lines, 

731 mask=lineMask, 

732 ) 

733 

734 @timeMethod 

735 def run(self, maskedImage): 

736 """Find and mask streaks in a masked image. 

737 

738 Finds streaks in the image and modifies maskedImage in place by adding a 

739 mask plane with any identified streaks. 

740 

741 Parameters 

742 ---------- 

743 maskedImage : `lsst.afw.image.Exposure` or `lsst.afw.image.maskedImage` 

744 The image in which to search for streaks. The mask detection plane 

745 corresponding to `config.detectedMaskPlane` must be set with the 

746 detected pixels. The mask will have a plane added with any detected 

747 streaks, and with the mask plane name set by 

748 self.config.streaksMaskPlane. 

749 

750 Returns 

751 ------- 

752 result : `lsst.pipe.base.Struct` 

753 Results as a struct with attributes: 

754 

755 ``originalLines`` 

756 Lines identified by kernel hough transform. 

757 ``lineClusters`` 

758 Lines grouped into clusters in rho-theta space. 

759 ``lines`` 

760 Final result for lines after line-profile fit. 

761 """ 

762 streaks = self.find(maskedImage) 

763 

764 if (self.config.streaksMaskPlane != "STREAK") and \ 

765 (self.config.streaksMaskPlane not in maskedImage.mask.getMaskPlaneDict()): 

766 maskedImage.mask.addMaskPlane(self.config.streaksMaskPlane) 

767 maskedImage.mask.array[streaks.mask] |= maskedImage.mask.getPlaneBitMask(self.config.streaksMaskPlane) 

768 

769 return pipeBase.Struct( 

770 lines=streaks.lines, 

771 lineClusters=streaks.lineClusters, 

772 originalLines=streaks.originalLines, 

773 ) 

774 

775 def _cannyFilter(self, image): 

776 """Apply a canny filter to the data in order to detect edges. 

777 

778 Parameters 

779 ---------- 

780 image : `np.ndarray` 

781 2-d image data on which to run filter. 

782 

783 Returns 

784 ------- 

785 cannyData : `np.ndarray` 

786 2-d image of edges found in input image. 

787 """ 

788 # Ensure that the pixels are zero or one. Change the datatype to 

789 # np.float64 to be compatible with the Canny filter routine. 

790 filterData = (image > 0).astype(np.float64) 

791 return canny(filterData, use_quantiles=True, sigma=0.1) 

792 

793 def _runKHT(self, image): 

794 """Run Kernel Hough Transform on image. 

795 

796 Parameters 

797 ---------- 

798 image : `np.ndarray` 

799 2-d image data on which to detect lines. 

800 

801 Returns 

802 ------- 

803 result : `LineCollection` 

804 Collection of detected lines, with their detected rho and theta 

805 coordinates. 

806 """ 

807 lines = lsst.kht.find_lines(image, self.config.clusterMinimumSize, 

808 self.config.clusterMinimumDeviation, self.config.delta, 

809 self.config.minimumKernelHeight, self.config.nSigma, 

810 self.config.absMinimumKernelHeight) 

811 self.log.info("The Kernel Hough Transform detected %s line(s)", len(lines)) 

812 

813 return LineCollection(lines.rho, lines.theta) 

814 

815 def _findClusters(self, lines): 

816 """Group lines that are close in parameter space and likely describe 

817 the same streak. 

818 

819 Parameters 

820 ---------- 

821 lines : `LineCollection` 

822 Collection of lines to group into clusters. 

823 

824 Returns 

825 ------- 

826 result : `LineCollection` 

827 Average `Line` for each cluster of `Line`s in the input 

828 `LineCollection`. 

829 """ 

830 # Scale variables by threshold bin-size variable so that rho and theta 

831 # are on the same scale. Since the clustering algorithm below stops when 

832 # the standard deviation <= 1, after rescaling each cluster will have a 

833 # standard deviation at or below the bin-size. 

834 x = lines.rhos / self.config.rhoBinSize 

835 y = lines.thetas / self.config.thetaBinSize 

836 X = np.array([x, y]).T 

837 nClusters = 1 

838 

839 # Put line parameters in clusters by starting with all in one, then 

840 # subdividing until the parameters of each cluster have std dev=1. 

841 # If nClusters == len(lines), each line will have its own 'cluster', so 

842 # the standard deviations of each cluster must be zero and the loop 

843 # is guaranteed to stop. 

844 while True: 

845 kmeans = KMeans(n_clusters=nClusters, n_init='auto').fit(X) 

846 clusterStandardDeviations = np.zeros((nClusters, 2)) 

847 for c in range(nClusters): 

848 inCluster = X[kmeans.labels_ == c] 

849 clusterStandardDeviations[c] = np.std(inCluster, axis=0)