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.setLineMask(line, maxStreakWidth=0, nSigmaMask=10, detectionMask=detectionMask) 

151 

152 def getLineXY(self, line): 

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

154 

155 Parameters 

156 ---------- 

157 line : `Line` 

158 Line for which to find the endpoints. 

159 

160 Returns 

161 ------- 

162 boxIntersections : `np.ndarray` 

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

164 """ 

165 theta = line.theta * u.deg 

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

167 # box. 

168 # Bottom: 

169 yA = -self._ymax / 2. 

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

171 # Left: 

172 xB = -self._xmax / 2. 

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

174 # Top: 

175 yC = self._ymax / 2. 

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

177 # Right: 

178 xD = self._xmax / 2. 

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

180 

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

182 [xB, yB], 

183 [xC, yC], 

184 [xD, yD]]) 

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

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

187 

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

189 # bounding box itself. 

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

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

192 boxIntersections = lineIntersections[inBox] 

193 

194 return boxIntersections 

195 

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

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

198 

199 Parameters 

200 ---------- 

201 line : `Line` 

202 Parameters of line in the image. 

203 maxStreakWidth : `float` 

204 Maximum width in pixels of streak mask. 

205 nSigmaMask : `float` 

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

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

208 Logger to use for reporting when maxStreakWidth is reached. 

209 detectionMask : `np.ndarray`, optional 

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

211 """ 

212 if line: 

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

214 radtheta = np.deg2rad(line.theta) 

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

216 

217 width = 2 * nSigmaMask * line.sigma 

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

219 if logger is not None: 

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

221 width, maxStreakWidth) 

222 width = maxStreakWidth 

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

224 self.lineMask = self.mask & m 

225 if detectionMask is not None: 

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

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

228 # the image. 

229 lineEnds = self.getLineXY(line) 

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

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

232 

233 radtheta = np.deg2rad(line.theta) 

234 costheta = np.cos(radtheta) 

235 sintheta = np.sin(radtheta) 

236 

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

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

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

240 + sintheta * xA 

241 - costheta * yA) 

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

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

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

245 maskDetections, 

246 statistic='sum', 

247 bins=lineBins) 

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

249 statistic='count', bins=lineBins) 

250 detectionFraction = detectionsAlongStreak / countAlongStreak 

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

252 emptyDetections = emptyRows[binnumber - 1] 

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

254 else: 

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

256 

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

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

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

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

261 

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

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

264 derivatives. 

265 

266 Parameters 

267 ---------- 

268 line : `Line` 

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

270 region. 

271 fitFlux : `bool` 

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

273 

274 Returns 

275 ------- 

276 model : `np.ndarray` 

277 Model in the masked region. 

278 dModel : `np.ndarray` 

279 Derivative of the model in the masked region. 

280 """ 

281 invSigma = line.sigma**-1 

282 # Calculate distance between pixels and line 

283 radtheta = np.deg2rad(line.theta) 

284 costheta = np.cos(radtheta) 

285 sintheta = np.sin(radtheta) 

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

287 distanceSquared = distance**2 

288 

289 # Calculate partial derivatives of distance 

290 drad = np.pi / 180 

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

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

293 

294 # Use pixel-line distances to make Moffat profile 

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

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

297 

298 if fitFlux: 

299 # Calculate line flux from profile and data 

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

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

302 else: 

303 # Approximately normalize the line 

304 flux = invSigma**-1 

305 if np.isnan(flux): 

306 flux = 0 

307 

308 model = flux * profile 

309 

310 # Calculate model derivatives 

311 fluxdProfile = flux * dProfile 

312 fluxdProfileInvSigma = fluxdProfile * invSigma**2 

313 dModeldRho = fluxdProfileInvSigma * dDistanceSqdRho 

314 dModeldTheta = fluxdProfileInvSigma * dDistanceSqdTheta 

315 dModeldInvSigma = fluxdProfile * distanceSquared * 2 * invSigma 

316 

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

318 return model, dModel 

319 

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

321 """Construct the line profile model. 

322 

323 Parameters 

324 ---------- 

325 line : `Line` 

326 Parameters of the line profile to model. 

327 fitFlux : `bool`, optional 

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

329 

330 Returns 

331 ------- 

332 finalModel : `np.ndarray` 

333 Model for line profile. 

334 """ 

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

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

337 finalModel[self.lineMask] = model 

338 return finalModel 

339 

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

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

342 

343 Parameters 

344 ---------- 

345 line : `Line` 

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

347 grad : `bool`, optional 

348 Whether or not to return the gradient and hessian. 

349 

350 Returns 

351 ------- 

352 reducedChi : `float` 

353 Reduced chi2 of the model. 

354 reducedDChi : `np.ndarray` 

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

356 reducedHessianChi : `np.ndarray` 

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

358 """ 

359 # Calculate chi2 

360 model, dModel = self._makeMaskedProfile(line) 

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

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

363 if not grad: 

364 return chi2.sum() / maskSize 

365 

366 # Calculate derivative and Hessian of chi2 

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

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

369 

370 reducedChi = chi2 / maskSize 

371 reducedDChi = derivChi2 / maskSize 

372 reducedHessianChi = hessianChi2 / maskSize 

373 return reducedChi, reducedDChi, reducedHessianChi 

374 

375 def _rejectOutliers(self, line): 

376 """Reject outlier pixels. 

377 

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

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

380 fitter ignores regions like the area around bright stars. 

381 

382 Parameters 

383 ---------- 

384 line : `Line` 

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

386 

387 Returns 

388 ------- 

389 nOutliers : `int` 

390 Number of outlier pixels. 

391 """ 

392 model, _ = self._makeMaskedProfile(line) 

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

394 

395 lineEnds = self.getLineXY(line) 

396 

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

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

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

400 else: 

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

402 return 0 

403 

404 radtheta = np.deg2rad(line.theta) 

405 costheta = np.cos(radtheta) 

406 sintheta = np.sin(radtheta) 

407 

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

409 - costheta * yA) 

410 

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

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

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

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

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

416 statistic='sum', bins=lineBins) 

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

418 statistic='count', bins=lineBins) 

419 

420 rChi2AlongStreak = chi2AlongStreak / countAlongStreak 

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

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

423 

424 outlierPix = outliers[binnumber - 1] 

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

426 tmpWeights[outlierPix] = 0 

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

428 

429 nOutliers = outlierPix.sum() 

430 

431 return nOutliers 

432 

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

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

435 

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

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

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

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

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

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

442 needed. 

443 

444 Parameters 

445 ---------- 

446 dChi2Tol : `float`, optional 

447 Change in Chi2 tolerated for fit convergence. 

448 maxIter : `int`, optional 

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

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

451 maximum provides a backup. 

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

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

454 

455 Returns 

456 ------- 

457 outline : `Line` 

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

459 fitFailure : `bool` 

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

461 """ 

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

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

464 

465 dChi2 = 1 

466 iter = 0 

467 oldChi2 = 0 

468 nOutliers = 1 

469 fitFailure = False 

470 

471 def line_search(c, dx): 

472 testx = x - c * dx 

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

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

475 

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

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

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

479 if chi2 == 0: 

480 break 

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

482 fitFailure = True 

483 if log is not None: 

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

485 break 

486 dChi2 = oldChi2 - chi2 

487 try: 

488 cholesky = scipy.linalg.cho_factor(A) 

489 except np.linalg.LinAlgError: 

490 fitFailure = True 

491 if log is not None: 

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

493 break 

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

495 

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

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

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

499 factor = 1 

500 else: 

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

502 

503 x -= factor * dx 

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

505 fitFailure = True 

506 break 

507 oldChi2 = chi2 

508 

509 nOutliers = self._rejectOutliers(line) 

510 iter += 1 

511 

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

513 return outline, fitFailure 

514 

515 

516class MaskStreaksConfig(pexConfig.Config): 

517 """Configuration parameters for `MaskStreaksTask`. 

518 """ 

519 

520 minimumKernelHeight = pexConfig.Field( 

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

522 dtype=float, 

523 default=0.0, 

524 ) 

525 absMinimumKernelHeight = pexConfig.Field( 

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

527 dtype=float, 

528 default=5, 

529 ) 

530 clusterMinimumSize = pexConfig.Field( 

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

532 dtype=int, 

533 default=50, 

534 ) 

535 clusterMinimumDeviation = pexConfig.Field( 

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

537 "line", 

538 dtype=int, 

539 default=2, 

540 ) 

541 delta = pexConfig.Field( 

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

543 dtype=float, 

544 default=0.2, 

545 ) 

546 nSigma = pexConfig.Field( 

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

548 "procedure", 

549 dtype=float, 

550 default=2, 

551 ) 

552 nSigmaMask = pexConfig.Field( 

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

554 dtype=float, 

555 default=5, 

556 ) 

557 rhoBinSize = pexConfig.Field( 

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

559 "clusters of detected lines", 

560 dtype=float, 

561 default=40, 

562 ) 

563 thetaBinSize = pexConfig.Field( 

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

565 "clusters of detected lines", 

566 dtype=float, 

567 default=2, 

568 ) 

569 invSigma = pexConfig.Field( 

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

571 "describing the profile of the streak", 

572 dtype=float, 

573 default=10.**-1, 

574 ) 

575 footprintThreshold = pexConfig.Field( 

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

577 "nanoJanskys", 

578 dtype=float, 

579 default=0.01 

580 ) 

581 dChi2Tolerance = pexConfig.Field( 

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

583 "fitting that is acceptable for convergence", 

584 dtype=float, 

585 default=0.1 

586 ) 

587 maxFitIter = pexConfig.Field( 

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

589 "acceptable for convergence", 

590 dtype=int, 

591 default=100 

592 ) 

593 detectedMaskPlane = pexConfig.Field( 

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

595 "estimate of streak locations", 

596 dtype=str, 

597 default="DETECTED" 

598 ) 

599 onlyMaskDetected = pexConfig.Field( 

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

601 "overlaps with the detection mask."), 

602 dtype=bool, 

603 default=True, 

604 ) 

605 streaksMaskPlane = pexConfig.Field( 

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

607 dtype=str, 

608 default="STREAK" 

609 ) 

610 badMaskPlanes = pexConfig.ListField( 

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

612 " detection"), 

613 dtype=str, 

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

615 ) 

616 maxStreakWidth = pexConfig.Field( 

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

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

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

620 "Set to 0 to disable.", 

621 dtype=float, 

622 default=0., 

623 ) 

624 saturatedDetectionsDilation = pexConfig.Field( 

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

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

627 dtype=int, 

628 default=250, 

629 ) 

630 

631 

632class MaskStreaksTask(pipeBase.Task): 

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

634 

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

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

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

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

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

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

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

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

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

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

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

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

647 (non-binary) image. 

648 """ 

649 

650 ConfigClass = MaskStreaksConfig 

651 _DefaultName = "maskStreaks" 

652 

653 @timeMethod 

654 def find(self, maskedImage): 

655 """Find streaks in a masked image. 

656 

657 Parameters 

658 ---------- 

659 maskedImage : `lsst.afw.image.maskedImage` 

660 The image in which to search for streaks. 

661 

662 Returns 

663 ------- 

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

665 Results as a struct with attributes: 

666 

667 ``originalLines`` 

668 Lines identified by kernel hough transform. 

669 ``lineClusters`` 

670 Lines grouped into clusters in rho-theta space. 

671 ``lines`` 

672 Final result for lines after line-profile fit. 

673 ``mask`` 

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

675 """ 

676 mask = maskedImage.mask 

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

678 

679 initEdges = self._cannyFilter(detectionMask) 

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

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

682 ignoreMask = mask.clone() 

683 

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

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

686 for sset in badMaskSpanSet: 

687 sset_dilated = sset.dilated(1) 

688 sset_dilated.clippedTo( 

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

690 

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

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

693 if self.config.saturatedDetectionsDilation: 

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

695 satMask = mask.getPlaneBitMask("SAT") 

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

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

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

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

700 for sset in satSpanSet: 

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

702 sset_dilated.clippedTo( 

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

704 

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

706 self.edges = initEdges & ~dilatedBadMask 

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

708 

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

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

711 fitLines = LineCollection([], []) 

712 clusters = LineCollection([], []) 

713 else: 

714 clusters = self._findClusters(self.lines) 

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

716 

717 if self.config.onlyMaskDetected: 

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

719 lineMask &= detectionMask.astype(bool) 

720 

721 return pipeBase.Struct( 

722 lines=fitLines, 

723 lineClusters=clusters, 

724 originalLines=self.lines, 

725 mask=lineMask, 

726 ) 

727 

728 @timeMethod 

729 def run(self, maskedImage): 

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

731 

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

733 mask plane with any identified streaks. 

734 

735 Parameters 

736 ---------- 

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

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

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

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

741 streaks, and with the mask plane name set by 

742 self.config.streaksMaskPlane. 

743 

744 Returns 

745 ------- 

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

747 Results as a struct with attributes: 

748 

749 ``originalLines`` 

750 Lines identified by kernel hough transform. 

751 ``lineClusters`` 

752 Lines grouped into clusters in rho-theta space. 

753 ``lines`` 

754 Final result for lines after line-profile fit. 

755 """ 

756 streaks = self.find(maskedImage) 

757 

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

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

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

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

762 

763 return pipeBase.Struct( 

764 lines=streaks.lines, 

765 lineClusters=streaks.lineClusters, 

766 originalLines=streaks.originalLines, 

767 ) 

768 

769 def _cannyFilter(self, image): 

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

771 

772 Parameters 

773 ---------- 

774 image : `np.ndarray` 

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

776 

777 Returns 

778 ------- 

779 cannyData : `np.ndarray` 

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

781 """ 

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

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

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

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

786 

787 def _runKHT(self, image): 

788 """Run Kernel Hough Transform on image. 

789 

790 Parameters 

791 ---------- 

792 image : `np.ndarray` 

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

794 

795 Returns 

796 ------- 

797 result : `LineCollection` 

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

799 coordinates. 

800 """ 

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

802 self.config.clusterMinimumDeviation, self.config.delta, 

803 self.config.minimumKernelHeight, self.config.nSigma, 

804 self.config.absMinimumKernelHeight) 

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

806 

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

808 

809 def _findClusters(self, lines): 

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

811 the same streak. 

812 

813 Parameters 

814 ---------- 

815 lines : `LineCollection` 

816 Collection of lines to group into clusters. 

817 

818 Returns 

819 ------- 

820 result : `LineCollection` 

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

822 `LineCollection`. 

823 """ 

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

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

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

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

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

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

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

831 nClusters = 1 

832 

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

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

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

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

837 # is guaranteed to stop. 

838 while True: 

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

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

841 for c in range(nClusters): 

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

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

844 # Are the rhos and thetas in each cluster all below the threshold? 

845 if (clusterStandardDeviations <= 1).all(): 

846 break 

847 nClusters += 1 

848 

849 # The cluster centers are final line estimates 

850 finalClusters = kmeans.cluster_centers_.T