Coverage for python / lsst / afw / cameraGeom / utils.py: 8%

443 statements  

« prev     ^ index     » next       coverage.py v7.13.5, created at 2026-05-07 01:18 -0700

1# This file is part of afw. 

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""" 

23Support for displaying cameraGeom objects. 

24""" 

25 

26__all__ = ['prepareWcsData', 'plotFocalPlane', 'makeImageFromAmp', 'calcRawCcdBBox', 'makeImageFromCcd', 

27 'FakeImageDataSource', 'ButlerImage', 'rawCallback', 'overlayCcdBoxes', 

28 'showAmp', 'showCcd', 'getCcdInCamBBoxList', 'getCameraImageBBox', 

29 'makeImageFromCamera', 'showCamera', 'makeFocalPlaneWcs', 'findAmp'] 

30 

31import logging 

32import math 

33import numpy 

34 

35import lsst.geom 

36from lsst.afw.fits import FitsError 

37import lsst.afw.geom as afwGeom 

38import lsst.afw.image as afwImage 

39import lsst.afw.math as afwMath 

40import lsst.afw.cameraGeom as afwCameraGeom 

41import lsst.daf.base as dafBase 

42import lsst.pex.exceptions as pexExceptions 

43 

44from ._rotateBBoxBy90 import rotateBBoxBy90 

45from ._assembleImage import assembleAmplifierImage, assembleAmplifierRawImage 

46from ._cameraGeom import FIELD_ANGLE, FOCAL_PLANE 

47from lsst.afw.display.utils import _getDisplayFromDisplayOrFrame 

48from lsst.afw.cameraGeom import DetectorType 

49 

50import lsst.afw.display as afwDisplay 

51import lsst.afw.display.utils as displayUtils 

52 

53_LOG = logging.getLogger(__name__) 

54 

55 

56def prepareWcsData(wcs, amp, isTrimmed=True): 

57 """Put Wcs from an Amp image into CCD coordinates 

58 

59 Parameters 

60 ---------- 

61 wcs : `lsst.afw.geom.SkyWcs` 

62 The WCS object to start from. 

63 amp : `lsst.afw.table.AmpInfoRecord` 

64 Amp object to use 

65 isTrimmed : `bool` 

66 Is the image to which the WCS refers trimmed of non-imaging pixels? 

67 

68 Returns 

69 ------- 

70 ampWcs : `lsst.afw.geom.SkyWcs` 

71 The modified WCS. 

72 """ 

73 if isTrimmed: 

74 ampBox = amp.getRawDataBBox() 

75 else: 

76 ampBox = amp.getRawBBox() 

77 ampCenter = lsst.geom.Point2D(ampBox.getDimensions()/2.0) 

78 wcs = afwGeom.makeFlippedWcs(wcs, amp.getRawFlipX(), amp.getRawFlipY(), ampCenter) 

79 # Shift WCS for trimming 

80 if isTrimmed: 

81 trim_shift = ampBox.getMin() - amp.getBBox().getMin() 

82 wcs = wcs.copyAtShiftedPixelOrigin(lsst.geom.Extent2D(-trim_shift.getX(), -trim_shift.getY())) 

83 # Account for shift of amp data in larger ccd matrix 

84 offset = amp.getRawXYOffset() 

85 return wcs.copyAtShiftedPixelOrigin(lsst.geom.Extent2D(offset)) 

86 

87 

88def plotFocalPlane(camera, fieldSizeDeg_x=0, fieldSizeDeg_y=None, dx=0.1, dy=0.1, figsize=(10., 10.), 

89 useIds=False, showFig=True, savePath=None): 

90 """Make a plot of the focal plane along with a set points that sample 

91 the field of view. 

92 

93 Parameters 

94 ---------- 

95 camera : `lsst.afw.cameraGeom.Camera` 

96 A camera object 

97 fieldSizeDeg_x : `float` 

98 Amount of the field to sample in x in degrees 

99 fieldSizeDeg_y : `float` or `None` 

100 Amount of the field to sample in y in degrees 

101 dx : `float` 

102 Spacing of sample points in x in degrees 

103 dy : `float` 

104 Spacing of sample points in y in degrees 

105 figsize : `tuple` containing two `float` 

106 Matplotlib style tuple indicating the size of the figure in inches 

107 useIds : `bool` 

108 Label detectors by name, not id? 

109 showFig : `bool` 

110 Display the figure on the screen? 

111 savePath : `str` or `None` 

112 If not `None`, save a copy of the figure to this name. 

113 """ 

114 try: 

115 from matplotlib.patches import Polygon 

116 from matplotlib.collections import PatchCollection 

117 import matplotlib.pyplot as plt 

118 except ImportError: 

119 raise ImportError( 

120 "Can't run plotFocalPlane: matplotlib has not been set up") 

121 

122 if fieldSizeDeg_x: 

123 if fieldSizeDeg_y is None: 

124 fieldSizeDeg_y = fieldSizeDeg_x 

125 

126 field_gridx, field_gridy = numpy.meshgrid( 

127 numpy.arange(0., fieldSizeDeg_x + dx, dx) - fieldSizeDeg_x/2., 

128 numpy.arange(0., fieldSizeDeg_y + dy, dy) - fieldSizeDeg_y/2.) 

129 field_gridx, field_gridy = field_gridx.flatten(), field_gridy.flatten() 

130 else: 

131 field_gridx, field_gridy = [], [] 

132 

133 xs = [] 

134 ys = [] 

135 pcolors = [] 

136 

137 # compute focal plane positions corresponding to field angles field_gridx, field_gridy 

138 posFieldAngleList = [lsst.geom.Point2D(x*lsst.geom.radians, y*lsst.geom.radians) 

139 for x, y in zip(field_gridx, field_gridy)] 

140 posFocalPlaneList = camera.transform(posFieldAngleList, FIELD_ANGLE, FOCAL_PLANE) 

141 for posFocalPlane in posFocalPlaneList: 

142 xs.append(posFocalPlane.getX()) 

143 ys.append(posFocalPlane.getY()) 

144 dets = camera.findDetectors(posFocalPlane, FOCAL_PLANE) 

145 if len(dets) > 0: 

146 pcolors.append('w') 

147 else: 

148 pcolors.append('k') 

149 

150 colorMap = {DetectorType.SCIENCE: 'b', DetectorType.FOCUS: 'y', 

151 DetectorType.GUIDER: 'g', DetectorType.WAVEFRONT: 'r'} 

152 

153 patches = [] 

154 colors = [] 

155 plt.figure(figsize=figsize) 

156 ax = plt.gca() 

157 xvals = [] 

158 yvals = [] 

159 for det in camera: 

160 corners = [(c.getX(), c.getY()) for c in det.getCorners(FOCAL_PLANE)] 

161 for corner in corners: 

162 xvals.append(corner[0]) 

163 yvals.append(corner[1]) 

164 colors.append(colorMap[det.getType()]) 

165 patches.append(Polygon(corners, closed=True)) 

166 center = det.getOrientation().getFpPosition() 

167 if det.getName() in ['R04_SW0', 'R04_SW1', 'R40_SW0', 'R40_SW1']: 

168 text_rotation = 'vertical' 

169 else: 

170 text_rotation = 'horizontal' 

171 ax.text(center.getX(), center.getY(), det.getId() if useIds else det.getName(), 

172 horizontalalignment='center', rotation=text_rotation, 

173 rotation_mode='anchor', size=6) 

174 patchCollection = PatchCollection(patches, alpha=0.6, facecolor=colors) 

175 ax.add_collection(patchCollection) 

176 ax.scatter(xs, ys, s=10, alpha=.7, linewidths=0., c=pcolors) 

177 ax.set_xlim(min(xvals) - abs(0.1*min(xvals)), 

178 max(xvals) + abs(0.1*max(xvals))) 

179 ax.set_ylim(min(yvals) - abs(0.1*min(yvals)), 

180 max(yvals) + abs(0.1*max(yvals))) 

181 ax.set_xlabel('Focal Plane X (mm)') 

182 ax.set_ylabel('Focal Plane Y (mm)') 

183 if savePath is not None: 

184 plt.savefig(savePath) 

185 if showFig: 

186 plt.show() 

187 

188 

189def makeImageFromAmp(amp, imValue=None, imageFactory=afwImage.ImageU, markSize=10, markValue=0, 

190 scaleGain=lambda gain: (gain*1000)//10): 

191 """Make an image from an amp object. 

192 

193 Since images are integer images by default, the gain needs to be scaled to 

194 give enough dynamic range to see variation from amp to amp. 

195 The scaling algorithm is assignable. 

196 

197 Parameters 

198 ---------- 

199 amp : `lsst.afw.table.AmpInfoRecord` 

200 Amp record to use for constructing the raw amp image. 

201 imValue : `float` or `None` 

202 Value to assign to the constructed image, or scaleGain(gain) if `None`. 

203 imageFactory : callable like `lsst.afw.image.Image` 

204 Type of image to construct. 

205 markSize : `float` 

206 Size of mark at read corner in pixels. 

207 markValue : `float` 

208 Value of pixels in the read corner mark. 

209 scaleGain : callable 

210 The function by which to scale the gain (must take a single argument). 

211 

212 Returns 

213 ------- 

214 ampImage : `lsst.afw.image` 

215 An untrimmed amp image, of the type produced by ``imageFactory``. 

216 """ 

217 bbox = amp.getRawBBox() 

218 dbbox = amp.getRawDataBBox() 

219 img = imageFactory(bbox) 

220 if imValue is None: 

221 img.set(int(scaleGain(amp.getGain()))) 

222 else: 

223 img.set(imValue) 

224 # Set the first pixel read to a different value 

225 markbbox = lsst.geom.Box2I() 

226 if amp.getReadoutCorner() == afwCameraGeom.ReadoutCorner.LL: 

227 markbbox.include(dbbox.getMin()) 

228 markbbox.include(dbbox.getMin() + lsst.geom.Extent2I(markSize, markSize)) 

229 elif amp.getReadoutCorner() == afwCameraGeom.ReadoutCorner.LR: 

230 cornerPoint = lsst.geom.Point2I(dbbox.getMaxX(), dbbox.getMinY()) 

231 markbbox.include(cornerPoint) 

232 markbbox.include(cornerPoint + lsst.geom.Extent2I(-markSize, markSize)) 

233 elif amp.getReadoutCorner() == afwCameraGeom.ReadoutCorner.UR: 

234 cornerPoint = lsst.geom.Point2I(dbbox.getMax()) 

235 markbbox.include(cornerPoint) 

236 markbbox.include(cornerPoint + lsst.geom.Extent2I(-markSize, -markSize)) 

237 elif amp.getReadoutCorner() == afwCameraGeom.ReadoutCorner.UL: 

238 cornerPoint = lsst.geom.Point2I(dbbox.getMinX(), dbbox.getMaxY()) 

239 markbbox.include(cornerPoint) 

240 markbbox.include(cornerPoint + lsst.geom.Extent2I(markSize, -markSize)) 

241 else: 

242 raise RuntimeError("Could not set readout corner") 

243 mimg = imageFactory(img, markbbox) 

244 mimg.set(markValue) 

245 return img 

246 

247 

248def calcRawCcdBBox(ccd): 

249 """Calculate the raw ccd bounding box. 

250 

251 Parameters 

252 ---------- 

253 ccd : `lsst.afw.cameraGeom.Detector` 

254 Detector for which to calculate the un-trimmed bounding box. 

255 

256 Returns 

257 ------- 

258 bbox : `lsst.geom.Box2I` or `None` 

259 Bounding box of the un-trimmed Detector, or `None` if there is not enough 

260 information to calculate raw BBox. 

261 """ 

262 bbox = lsst.geom.Box2I() 

263 for amp in ccd: 

264 tbbox = amp.getRawBBox() 

265 tbbox.shift(amp.getRawXYOffset()) 

266 bbox.include(tbbox) 

267 return bbox 

268 

269 

270def makeImageFromCcd(ccd, isTrimmed=True, showAmpGain=True, imageFactory=afwImage.ImageU, rcMarkSize=10, 

271 binSize=1): 

272 """Make an Image of a CCD. 

273 

274 Parameters 

275 ---------- 

276 ccd : `lsst.afw.cameraGeom.Detector` 

277 Detector to use in making the image. 

278 isTrimmed : `bool` 

279 Assemble a trimmed Detector image. 

280 showAmpGain : `bool` 

281 Use the per-amp gain to color the pixels in the image? 

282 imageFactory : callable like `lsst.afw.image.Image` 

283 Image type to generate. 

284 rcMarkSize : `float` 

285 Size of the mark to make in the amp images at the read corner. 

286 binSize : `int` 

287 Bin the image by this factor in both dimensions. 

288 

289 Returns 

290 ------- 

291 image : `lsst.afw.image.Image` 

292 Image of the Detector (type returned by ``imageFactory``). 

293 """ 

294 ampImages = [] 

295 index = 0 

296 if isTrimmed: 

297 bbox = ccd.getBBox() 

298 else: 

299 bbox = calcRawCcdBBox(ccd) 

300 for amp in ccd: 

301 if showAmpGain: 

302 ampImages.append(makeImageFromAmp( 

303 amp, imageFactory=imageFactory, markSize=rcMarkSize)) 

304 else: 

305 ampImages.append(makeImageFromAmp(amp, imValue=(index + 1)*1000, 

306 imageFactory=imageFactory, markSize=rcMarkSize)) 

307 index += 1 

308 

309 if len(ampImages) > 0: 

310 ccdImage = imageFactory(bbox) 

311 for ampImage, amp in zip(ampImages, ccd): 

312 if isTrimmed: 

313 assembleAmplifierImage(ccdImage, ampImage, amp) 

314 else: 

315 assembleAmplifierRawImage(ccdImage, ampImage, amp) 

316 else: 

317 if not isTrimmed: 

318 raise RuntimeError( 

319 "Cannot create untrimmed CCD without amps with raw information") 

320 ccdImage = imageFactory(ccd.getBBox()) 

321 ccdImage = afwMath.binImage(ccdImage, binSize) 

322 return ccdImage 

323 

324 

325class FakeImageDataSource: 

326 """A class to retrieve synthetic images for display by the show* methods 

327 

328 Parameters 

329 ---------- 

330 isTrimmed : `bool` 

331 Should amps be trimmed? 

332 verbose : `bool` 

333 Be chatty? 

334 background : `float` 

335 The value of any pixels that lie outside the CCDs. 

336 showAmpGain : `bool` 

337 Color the amp segments with the gain of the amp? 

338 markSize : `float` 

339 Size of the side of the box used to mark the read corner. 

340 markValue : `float` 

341 Value to assign the read corner mark. 

342 ampImValue : `float` or `None` 

343 Value to assign to amps; scaleGain(gain) is used if `None`. 

344 scaleGain : callable 

345 Function to scale the gain by. 

346 """ 

347 def __init__(self, isTrimmed=True, verbose=False, background=numpy.nan, 

348 showAmpGain=True, markSize=10, markValue=0, 

349 ampImValue=None, scaleGain=lambda gain: (gain*1000)//10): 

350 self.isTrimmed = isTrimmed 

351 self.verbose = verbose 

352 self.background = background 

353 self.showAmpGain = showAmpGain 

354 self.markSize = markSize 

355 self.markValue = markValue 

356 self.ampImValue = ampImValue 

357 self.scaleGain = scaleGain 

358 

359 def getCcdImage(self, det, imageFactory, binSize): 

360 """Return a CCD image for the detector and the (possibly updated) Detector. 

361 

362 Parameters 

363 ---------- 

364 det : `lsst.afw.cameraGeom.Detector` 

365 Detector to use for making the image. 

366 imageFactory : callable like `lsst.afw.image.Image` 

367 Image constructor for making the image. 

368 binSize : `int` 

369 Bin the image by this factor in both dimensions. 

370 

371 Returns 

372 ------- 

373 ccdImage : `lsst.afw.image.Image` 

374 The constructed image. 

375 """ 

376 ccdImage = makeImageFromCcd(det, isTrimmed=self.isTrimmed, showAmpGain=self.showAmpGain, 

377 imageFactory=imageFactory, binSize=binSize) 

378 return afwMath.rotateImageBy90(ccdImage, det.getOrientation().getNQuarter()), det 

379 

380 def getAmpImage(self, amp, imageFactory): 

381 """Return an amp segment image. 

382 

383 Parameters 

384 ---------- 

385 amp : `lsst.afw.table.AmpInfoTable` 

386 AmpInfoTable for this amp. 

387 imageFactory : callable like `lsst.afw.image.Image` 

388 Image constructor for making the image. 

389 

390 Returns 

391 ------- 

392 ampImage : `lsst.afw.image.Image` 

393 The constructed image. 

394 """ 

395 ampImage = makeImageFromAmp(amp, imValue=self.ampImValue, imageFactory=imageFactory, 

396 markSize=self.markSize, markValue=self.markValue, 

397 scaleGain=self.scaleGain) 

398 if self.isTrimmed: 

399 ampImage = ampImage.Factory(ampImage, amp.getRawDataBBox()) 

400 return ampImage 

401 

402 

403class ButlerImage(FakeImageDataSource): 

404 """A class to return an Image of a given Ccd using the butler. 

405 

406 Parameters 

407 ---------- 

408 butler : `lsst.daf.butler.Butler` or `None` 

409 The butler to use. If `None`, an empty image is returned. Assumes that 

410 the instrument was specified during butler construction or is included 

411 in the ``kwargs`` parameter. 

412 type : `str` 

413 The type of image to read (e.g. raw, bias, flat, calexp). 

414 isTrimmed : `bool` 

415 If true, the showCamera command expects to be given trimmed images. 

416 verbose : `bool` 

417 Be chatty (in particular, log any error messages from the butler)? 

418 background : `float` 

419 The value of any pixels that lie outside the CCDs. 

420 callback : callable 

421 A function called with (image, detector, butler) for every image, which 

422 returns the image to be displayed (e.g. rawCallback). The image must 

423 be of the correct size, allowing for the value of isTrimmed. 

424 *args : `list` 

425 Passed to the base class constructor. 

426 **kwargs : `dict` 

427 Passed to the butler. 

428 

429 Notes 

430 ----- 

431 You can define a short named function as a callback:: 

432 

433 def callback(im, ccd, imageSource): 

434 return cameraGeom.utils.rawCallback(im, ccd, imageSource, correctGain=True) 

435 """ 

436 def __init__(self, butler=None, type="raw", 

437 isTrimmed=True, verbose=False, background=numpy.nan, 

438 callback=None, *args, **kwargs): 

439 super().__init__(*args) 

440 self.isTrimmed = isTrimmed 

441 self.type = type 

442 self.butler = butler 

443 self.kwargs = kwargs 

444 self.isRaw = False 

445 self.background = background 

446 self.verbose = verbose 

447 self.callback = callback 

448 

449 def _prepareImage(self, ccd, im, binSize, allowRotate=True): 

450 if binSize > 1: 

451 im = afwMath.binImage(im, binSize) 

452 

453 if allowRotate: 

454 im = afwMath.rotateImageBy90( 

455 im, ccd.getOrientation().getNQuarter()) 

456 

457 return im 

458 

459 def getCcdImage(self, ccd, imageFactory=afwImage.ImageF, binSize=1, asMaskedImage=False): 

460 """Return an image of the specified ccd, and also the (possibly updated) ccd""" 

461 

462 log = _LOG.getChild("ButlerImage") 

463 

464 if self.isTrimmed: 

465 bbox = ccd.getBBox() 

466 else: 

467 bbox = calcRawCcdBBox(ccd) 

468 

469 im = None 

470 if self.butler is not None: 

471 err = None 

472 try: 

473 im = self.butler.get(self.type, detector=ccd.getId(), **self.kwargs) 

474 except FitsError as e: 

475 err = IOError(e.args[0].split('\n')[0]) # It's a very chatty error 

476 except Exception as e: # try a different dataId 

477 err = e 

478 else: 

479 ccd = im.getDetector() # possibly modified by assembleCcdTask 

480 

481 if im: 

482 if asMaskedImage: 

483 im = im.getMaskedImage() 

484 else: 

485 im = im.getMaskedImage().getImage() 

486 else: 

487 if self.verbose: 

488 # Lost by jupyterlab. 

489 print(f"Reading {ccd.getId()}: {err}") 

490 

491 log.warning("Reading %s: %s", ccd.getId(), err) 

492 

493 if im is None: 

494 return self._prepareImage(ccd, imageFactory(*bbox.getDimensions()), binSize), ccd 

495 

496 if self.type == "raw": 

497 if hasattr(im, 'convertF'): 

498 im = im.convertF() 

499 if False and self.callback is None: # we need to trim the raw image 

500 self.callback = rawCallback 

501 

502 allowRotate = True 

503 if self.callback: 

504 try: 

505 im = self.callback(im, ccd, imageSource=self) 

506 except Exception: 

507 if self.verbose: 

508 log.exception("callback failed.") 

509 im = imageFactory(*bbox.getDimensions()) 

510 else: 

511 allowRotate = False # the callback was responsible for any rotations 

512 

513 return self._prepareImage(ccd, im, binSize, allowRotate=allowRotate), ccd 

514 

515 

516def rawCallback(im, ccd=None, imageSource=None, 

517 correctGain=False, subtractBias=False, convertToFloat=False, obeyNQuarter=True): 

518 """A callback function that may or may not subtract bias/correct gain/trim 

519 a raw image. 

520 

521 Parameters 

522 ---------- 

523 im : `lsst.afw.image.Image` or `lsst.afw.image.MaskedImage` or `lsst.afw.image.Exposure` 

524 An image of a chip, ready to be binned and maybe rotated. 

525 ccd : `lsst.afw.cameraGeom.Detector` or `None` 

526 The Detector; if `None` assume that im is an exposure and extract its Detector. 

527 imageSource : `FakeImageDataSource` or `None` 

528 Source to get ccd images. Must have a `getCcdImage()` method. 

529 correctGain : `bool` 

530 Correct each amplifier for its gain? 

531 subtractBias : `bool` 

532 Subtract the bias from each amplifier? 

533 convertToFloat : `bool` 

534 Convert ``im`` to floating point if possible. 

535 obeyNQuarter : `bool` 

536 Obey nQuarter from the Detector (default: True) 

537 

538 Returns 

539 ------- 

540 image : `lsst.afw.image.Image` like 

541 The constructed image (type returned by ``im.Factory``). 

542 

543 Notes 

544 ----- 

545 If imageSource is derived from ButlerImage, imageSource.butler is available. 

546 """ 

547 if ccd is None: 

548 ccd = im.getDetector() 

549 if hasattr(im, "getMaskedImage"): 

550 im = im.getMaskedImage() 

551 if convertToFloat and hasattr(im, "convertF"): 

552 im = im.convertF() 

553 

554 isTrimmed = imageSource.isTrimmed 

555 if isTrimmed: 

556 bbox = ccd.getBBox() 

557 else: 

558 bbox = calcRawCcdBBox(ccd) 

559 

560 ampImages = [] 

561 for a in ccd: 

562 if isTrimmed: 

563 data = im[a.getRawDataBBox()] 

564 else: 

565 data = im 

566 

567 if subtractBias: 

568 bias = im[a.getRawHorizontalOverscanBBox()] 

569 data -= afwMath.makeStatistics(bias, afwMath.MEANCLIP).getValue() 

570 if correctGain: 

571 data *= a.getGain() 

572 

573 ampImages.append(data) 

574 

575 ccdImage = im.Factory(bbox) 

576 for ampImage, amp in zip(ampImages, ccd): 

577 if isTrimmed: 

578 assembleAmplifierImage(ccdImage, ampImage, amp) 

579 else: 

580 assembleAmplifierRawImage(ccdImage, ampImage, amp) 

581 

582 if obeyNQuarter: 

583 nQuarter = ccd.getOrientation().getNQuarter() 

584 ccdImage = afwMath.rotateImageBy90(ccdImage, nQuarter) 

585 

586 return ccdImage 

587 

588 

589def overlayCcdBoxes(ccd, untrimmedCcdBbox=None, nQuarter=0, 

590 isTrimmed=False, ccdOrigin=(0, 0), display=None, binSize=1): 

591 """Overlay bounding boxes on an image display. 

592 

593 Parameters 

594 ---------- 

595 ccd : `lsst.afw.cameraGeom.Detector` 

596 Detector to iterate for the amp bounding boxes. 

597 untrimmedCcdBbox : `lsst.geom.Box2I` or `None` 

598 Bounding box of the un-trimmed Detector. 

599 nQuarter : `int` 

600 number of 90 degree rotations to apply to the bounding boxes (used for rotated chips). 

601 isTrimmed : `bool` 

602 Is the Detector image over which the boxes are layed trimmed? 

603 ccdOrigin : `tuple` of `float` 

604 Detector origin relative to the parent origin if in a larger pixel grid. 

605 display : `lsst.afw.display.Display` 

606 Image display to display on. 

607 binSize : `int` 

608 Bin the image by this factor in both dimensions. 

609 

610 Notes 

611 ----- 

612 The colours are: 

613 - Entire detector GREEN 

614 - All data for amp GREEN 

615 - HorizontalPrescan YELLOW 

616 - HorizontalOverscan RED 

617 - Data BLUE 

618 - VerticalOverscan MAGENTA 

619 - VerticalOverscan MAGENTA 

620 """ 

621 if not display: # should be second parameter, and not defaulted!! 

622 raise RuntimeError("Please specify a display") 

623 

624 if untrimmedCcdBbox is None: 

625 if isTrimmed: 

626 untrimmedCcdBbox = ccd.getBBox() 

627 else: 

628 untrimmedCcdBbox = lsst.geom.Box2I() 

629 for a in ccd.getAmplifiers(): 

630 bbox = a.getRawBBox() 

631 untrimmedCcdBbox.include(bbox) 

632 

633 with display.Buffering(): 

634 ccdDim = untrimmedCcdBbox.getDimensions() 

635 ccdBbox = rotateBBoxBy90(untrimmedCcdBbox, nQuarter, ccdDim) 

636 for amp in ccd: 

637 if isTrimmed: 

638 ampbbox = amp.getBBox() 

639 else: 

640 ampbbox = amp.getRawBBox() 

641 if nQuarter != 0: 

642 ampbbox = rotateBBoxBy90(ampbbox, nQuarter, ccdDim) 

643 

644 displayUtils.drawBBox(ampbbox, origin=ccdOrigin, borderWidth=0.49, 

645 display=display, bin=binSize) 

646 

647 if not isTrimmed: 

648 for bbox, ctype in ((amp.getRawHorizontalOverscanBBox(), afwDisplay.RED), 

649 (amp.getRawDataBBox(), afwDisplay.BLUE), 

650 (amp.getRawVerticalOverscanBBox(), 

651 afwDisplay.MAGENTA), 

652 (amp.getRawPrescanBBox(), afwDisplay.YELLOW)): 

653 if nQuarter != 0: 

654 bbox = rotateBBoxBy90(bbox, nQuarter, ccdDim) 

655 displayUtils.drawBBox(bbox, origin=ccdOrigin, borderWidth=0.49, ctype=ctype, 

656 display=display, bin=binSize) 

657 # Label each Amp 

658 xc, yc = ((ampbbox.getMin()[0] + ampbbox.getMax()[0])//2, 

659 (ampbbox.getMin()[1] + ampbbox.getMax()[1])//2) 

660 # 

661 # Rotate the amp labels too 

662 # 

663 if nQuarter == 0: 

664 c, s = 1, 0 

665 elif nQuarter == 1: 

666 c, s = 0, -1 

667 elif nQuarter == 2: 

668 c, s = -1, 0 

669 elif nQuarter == 3: 

670 c, s = 0, 1 

671 c, s = 1, 0 

672 ccdHeight = ccdBbox.getHeight() 

673 ccdWidth = ccdBbox.getWidth() 

674 xc -= 0.5*ccdHeight 

675 yc -= 0.5*ccdWidth 

676 

677 xc, yc = 0.5*ccdHeight + c*xc + s*yc, 0.5*ccdWidth + -s*xc + c*yc 

678 

679 if ccdOrigin: 

680 xc += ccdOrigin[0] 

681 yc += ccdOrigin[1] 

682 display.dot(str(amp.getName()), xc/binSize, 

683 yc/binSize, textAngle=nQuarter*90) 

684 

685 displayUtils.drawBBox(ccdBbox, origin=ccdOrigin, 

686 borderWidth=0.49, ctype=afwDisplay.MAGENTA, display=display, bin=binSize) 

687 

688 

689def showAmp(amp, imageSource=FakeImageDataSource(isTrimmed=False), display=None, overlay=True, 

690 imageFactory=afwImage.ImageU): 

691 """Show an amp in an image display. 

692 

693 Parameters 

694 ---------- 

695 amp : `lsst.afw.tables.AmpInfoRecord` 

696 Amp record to use in display. 

697 imageSource : `FakeImageDataSource` or `None` 

698 Source for getting the amp image. Must have a ``getAmpImage()`` method. 

699 display : `lsst.afw.display.Display` 

700 Image display to use. 

701 overlay : `bool` 

702 Overlay bounding boxes? 

703 imageFactory : callable like `lsst.afw.image.Image` 

704 Type of image to display (only used if ampImage is `None`). 

705 """ 

706 if not display: 

707 display = _getDisplayFromDisplayOrFrame(display) 

708 

709 ampImage = imageSource.getAmpImage(amp, imageFactory=imageFactory) 

710 ampImSize = ampImage.getDimensions() 

711 title = amp.getName() 

712 display.mtv(ampImage, title=title) 

713 if overlay: 

714 with display.Buffering(): 

715 if ampImSize == amp.getRawBBox().getDimensions(): 

716 bboxes = [(amp.getRawBBox(), 0.49, afwDisplay.GREEN), ] 

717 xy0 = bboxes[0][0].getMin() 

718 bboxes.append( 

719 (amp.getRawHorizontalOverscanBBox(), 0.49, afwDisplay.RED)) 

720 bboxes.append((amp.getRawDataBBox(), 0.49, afwDisplay.BLUE)) 

721 bboxes.append((amp.getRawPrescanBBox(), 

722 0.49, afwDisplay.YELLOW)) 

723 bboxes.append((amp.getRawVerticalOverscanBBox(), 

724 0.49, afwDisplay.MAGENTA)) 

725 else: 

726 bboxes = [(amp.getBBox(), 0.49, None), ] 

727 xy0 = bboxes[0][0].getMin() 

728 

729 for bbox, borderWidth, ctype in bboxes: 

730 if bbox.isEmpty(): 

731 continue 

732 bbox = lsst.geom.Box2I(bbox) 

733 bbox.shift(-lsst.geom.ExtentI(xy0)) 

734 displayUtils.drawBBox( 

735 bbox, borderWidth=borderWidth, ctype=ctype, display=display) 

736 

737 

738def showCcd(ccd, imageSource=FakeImageDataSource(), display=None, overlay=True, 

739 imageFactory=afwImage.ImageF, binSize=1, inCameraCoords=False): 

740 """Show a CCD on display. 

741 

742 Parameters 

743 ---------- 

744 ccd : `lsst.afw.cameraGeom.Detector` 

745 Detector to use in display. 

746 imageSource : `FakeImageDataSource` or `None` 

747 Source to get ccd images. Must have a ``getCcdImage()`` method. 

748 display : `lsst.afw.display.Display` 

749 image display to use. 

750 overlay : `bool` 

751 Show amp bounding boxes on the displayed image? 

752 imageFactory : callable like `lsst.afw.image.Image` 

753 The image factory to use in generating the images. 

754 binSize : `int` 

755 Bin the image by this factor in both dimensions. 

756 inCameraCoords : `bool` 

757 Show the Detector in camera coordinates? 

758 """ 

759 display = _getDisplayFromDisplayOrFrame(display) 

760 

761 ccdOrigin = lsst.geom.Point2I(0, 0) 

762 nQuarter = 0 

763 ccdImage, ccd = imageSource.getCcdImage( 

764 ccd, imageFactory=imageFactory, binSize=binSize) 

765 

766 ccdBbox = ccdImage.getBBox() 

767 if ccdBbox.getDimensions() == ccd.getBBox().getDimensions(): 

768 isTrimmed = True 

769 else: 

770 isTrimmed = False 

771 

772 if inCameraCoords: 

773 nQuarter = ccd.getOrientation().getNQuarter() 

774 ccdImage = afwMath.rotateImageBy90(ccdImage, nQuarter) 

775 title = "{} [{}]".format(ccd.getName(), ccd.getId()) 

776 if isTrimmed: 

777 title += "(trimmed)" 

778 

779 if display: 

780 display.mtv(ccdImage, title=title) 

781 

782 if overlay: 

783 overlayCcdBoxes(ccd, ccdBbox, nQuarter, isTrimmed, 

784 ccdOrigin, display, binSize) 

785 

786 return ccdImage 

787 

788 

789def getCcdInCamBBoxList(ccdList, binSize, pixelSize_o, origin): 

790 """Get the bounding boxes of a list of Detectors within a camera sized pixel grid 

791 

792 Parameters 

793 ---------- 

794 ccdList : `lsst.afw.cameraGeom.Detector` 

795 List of Detector. 

796 binSize : `int` 

797 Bin the image by this factor in both dimensions. 

798 pixelSize_o : `float` 

799 Size of the pixel in mm. 

800 origin : `int` 

801 Origin of the camera pixel grid in pixels. 

802 

803 Returns 

804 ------- 

805 boxList : `list` [`lsst.geom.Box2I`] 

806 A list of bounding boxes in camera pixel coordinates. 

807 """ 

808 boxList = [] 

809 for ccd in ccdList: 

810 if not pixelSize_o == ccd.getPixelSize(): 

811 raise RuntimeError( 

812 "Cameras with detectors with different pixel scales are not currently supported") 

813 

814 dbbox = lsst.geom.Box2D() 

815 for corner in ccd.getCorners(FOCAL_PLANE): 

816 dbbox.include(corner) 

817 llc = dbbox.getMin() 

818 nQuarter = ccd.getOrientation().getNQuarter() 

819 cbbox = ccd.getBBox() 

820 ex = cbbox.getDimensions().getX()//binSize 

821 ey = cbbox.getDimensions().getY()//binSize 

822 bbox = lsst.geom.Box2I( 

823 cbbox.getMin(), lsst.geom.Extent2I(int(ex), int(ey))) 

824 bbox = rotateBBoxBy90(bbox, nQuarter, bbox.getDimensions()) 

825 bbox.shift(lsst.geom.Extent2I(int(llc.getX()//pixelSize_o.getX()/binSize), 

826 int(llc.getY()//pixelSize_o.getY()/binSize))) 

827 bbox.shift(lsst.geom.Extent2I(-int(origin.getX()//binSize), 

828 -int(origin.getY())//binSize)) 

829 boxList.append(bbox) 

830 return boxList 

831 

832 

833def getCameraImageBBox(camBbox, pixelSize, bufferSize): 

834 """Get the bounding box of a camera sized image in pixels 

835 

836 Parameters 

837 ---------- 

838 camBbox : `lsst.geom.Box2D` 

839 Camera bounding box in focal plane coordinates (mm). 

840 pixelSize : `float` 

841 Size of a detector pixel in mm. 

842 bufferSize : `int` 

843 Buffer around edge of image in pixels. 

844 

845 Returns 

846 ------- 

847 box : `lsst.geom.Box2I` 

848 The resulting bounding box. 

849 """ 

850 pixMin = lsst.geom.Point2I(int(camBbox.getMinX()//pixelSize.getX()), 

851 int(camBbox.getMinY()//pixelSize.getY())) 

852 pixMax = lsst.geom.Point2I(int(camBbox.getMaxX()//pixelSize.getX()), 

853 int(camBbox.getMaxY()//pixelSize.getY())) 

854 retBox = lsst.geom.Box2I(pixMin, pixMax) 

855 retBox.grow(bufferSize) 

856 return retBox 

857 

858 

859def makeImageFromCamera(camera, detectorNameList=None, background=numpy.nan, bufferSize=10, 

860 imageSource=FakeImageDataSource(), imageFactory=afwImage.ImageU, binSize=1): 

861 """Make an Image of a Camera. 

862 

863 Put each detector's image in the correct location and orientation on the 

864 focal plane. The input images can be binned to an integer fraction of their 

865 original bboxes. 

866 

867 Parameters 

868 ---------- 

869 camera : `lsst.afw.cameraGeom.Camera` 

870 Camera object to use to make the image. 

871 detectorNameList : `list` [`str`] 

872 List of detector names from ``camera`` to use in building the image. 

873 Use all Detectors if `None`. 

874 background : `float` 

875 Value to use where there is no Detector. 

876 bufferSize : `int` 

877 Size of border in binned pixels to make around the camera image. 

878 imageSource : `FakeImageDataSource` or `None` 

879 Source to get ccd images. Must have a ``getCcdImage()`` method. 

880 imageFactory : callable like `lsst.afw.image.Image` 

881 Type of image to build. 

882 binSize : `int` 

883 Bin the image by this factor in both dimensions. 

884 

885 Returns 

886 ------- 

887 image : `lsst.afw.image.Image` 

888 Image of the entire camera. 

889 """ 

890 log = _LOG.getChild("makeImageFromCamera") 

891 

892 if detectorNameList is None: 

893 ccdList = camera 

894 else: 

895 ccdList = [camera[name] for name in detectorNameList] 

896 

897 if detectorNameList is None: 

898 camBbox = camera.getFpBBox() 

899 else: 

900 camBbox = lsst.geom.Box2D() 

901 for detName in detectorNameList: 

902 for corner in camera[detName].getCorners(FOCAL_PLANE): 

903 camBbox.include(corner) 

904 

905 pixelSize_o = camera[next(camera.getNameIter())].getPixelSize() 

906 camBbox = getCameraImageBBox(camBbox, pixelSize_o, bufferSize*binSize) 

907 origin = camBbox.getMin() 

908 

909 camIm = imageFactory(int(math.ceil(camBbox.getDimensions().getX()/binSize)), 

910 int(math.ceil(camBbox.getDimensions().getY()/binSize))) 

911 camIm[:] = imageSource.background 

912 

913 assert imageSource.isTrimmed, "isTrimmed is False isn't supported by getCcdInCamBBoxList" 

914 

915 boxList = getCcdInCamBBoxList(ccdList, binSize, pixelSize_o, origin) 

916 for det, bbox in zip(ccdList, boxList): 

917 im = imageSource.getCcdImage(det, imageFactory, binSize)[0] 

918 if im is None: 

919 continue 

920 

921 imView = camIm.Factory(camIm, bbox, afwImage.LOCAL) 

922 try: 

923 imView[:] = im 

924 except pexExceptions.LengthError: 

925 log.exception("Unable to fit image for detector \"%s\" into image of camera.", det.getName()) 

926 

927 return camIm 

928 

929 

930def showCamera(camera, imageSource=FakeImageDataSource(), imageFactory=afwImage.ImageF, 

931 detectorNameList=None, binSize=10, bufferSize=10, overlay=True, title="", 

932 showWcs=None, ctype=afwDisplay.GREEN, textSize=1.25, originAtCenter=True, display=None, 

933 **kwargs): 

934 """Show a Camera on display, with the specified display. 

935 

936 The rotation of the sensors is snapped to the nearest multiple of 90 deg. 

937 Also note that the pixel size is constant over the image array. The lower 

938 left corner (LLC) of each sensor amp is snapped to the LLC of the pixel 

939 containing the LLC of the image. 

940 

941 Parameters 

942 ---------- 

943 camera : `lsst.afw.cameraGeom.Camera` 

944 Camera object to use to make the image. 

945 imageSource : `FakeImageDataSource` or `None` 

946 Source to get ccd images. Must have a ``getCcdImage()`` method. 

947 imageFactory : `lsst.afw.image.Image` 

948 Type of image to make 

949 detectorNameList : `list` [`str`] or `None` 

950 List of detector names from `camera` to use in building the image. 

951 Use all Detectors if `None`. 

952 binSize : `int` 

953 Bin the image by this factor in both dimensions. 

954 bufferSize : `int` 

955 Size of border in binned pixels to make around the camera image. 

956 overlay : `bool` 

957 Overlay Detector IDs and boundaries? 

958 title : `str` 

959 Title to use in display. 

960 showWcs : `bool` 

961 Include a WCS in the display? 

962 ctype : `lsst.afw.display.COLOR` or `str` 

963 Color to use when drawing Detector boundaries. 

964 textSize : `float` 

965 Size of detector labels 

966 originAtCenter : `bool` 

967 Put origin of the camera WCS at the center of the image? 

968 If `False`, the origin will be at the lower left. 

969 display : `lsst.afw.display` 

970 Image display on which to display. 

971 **kwargs : 

972 All remaining keyword arguments are passed to makeImageFromCamera 

973 

974 Returns 

975 ------- 

976 image : `lsst.afw.image.Image` 

977 The mosaic image. 

978 """ 

979 display = _getDisplayFromDisplayOrFrame(display) 

980 

981 if binSize < 1: 

982 binSize = 1 

983 cameraImage = makeImageFromCamera(camera, detectorNameList=detectorNameList, bufferSize=bufferSize, 

984 imageSource=imageSource, imageFactory=imageFactory, binSize=binSize, 

985 **kwargs) 

986 

987 if detectorNameList is None: 

988 ccdList = [camera[name] for name in camera.getNameIter()] 

989 else: 

990 ccdList = [camera[name] for name in detectorNameList] 

991 

992 if detectorNameList is None: 

993 camBbox = camera.getFpBBox() 

994 else: 

995 camBbox = lsst.geom.Box2D() 

996 for detName in detectorNameList: 

997 for corner in camera[detName].getCorners(FOCAL_PLANE): 

998 camBbox.include(corner) 

999 pixelSize = ccdList[0].getPixelSize() 

1000 

1001 if showWcs: 

1002 if originAtCenter: 

1003 wcsReferencePixel = lsst.geom.Box2D( 

1004 cameraImage.getBBox()).getCenter() 

1005 else: 

1006 wcsReferencePixel = lsst.geom.Point2I(0, 0) 

1007 wcs = makeFocalPlaneWcs(pixelSize*binSize, wcsReferencePixel) 

1008 else: 

1009 wcs = None 

1010 

1011 if display: 

1012 if title == "": 

1013 title = camera.getName() 

1014 display.mtv(cameraImage, title=title, wcs=wcs) 

1015 

1016 if overlay: 

1017 with display.Buffering(): 

1018 camBbox = getCameraImageBBox( 

1019 camBbox, pixelSize, bufferSize*binSize) 

1020 bboxList = getCcdInCamBBoxList( 

1021 ccdList, binSize, pixelSize, camBbox.getMin()) 

1022 for bbox, ccd in zip(bboxList, ccdList): 

1023 nQuarter = ccd.getOrientation().getNQuarter() 

1024 # borderWidth to 0.5 to align with the outside edge of the 

1025 # pixel 

1026 displayUtils.drawBBox( 

1027 bbox, borderWidth=0.5, ctype=ctype, display=display) 

1028 dims = bbox.getDimensions() 

1029 ccdLabel = "{}\n[{}]".format(ccd.getName(), ccd.getId()) 

1030 display.dot(ccdLabel, bbox.getMinX() + dims.getX()/2, bbox.getMinY() + dims.getY()/2, 

1031 ctype=ctype, size=textSize, textAngle=nQuarter*90) 

1032 

1033 return cameraImage 

1034 

1035 

1036def makeFocalPlaneWcs(pixelSize, referencePixel): 

1037 """Make a WCS for the focal plane geometry 

1038 (i.e. one that returns positions in "mm") 

1039 

1040 Parameters 

1041 ---------- 

1042 pixelSize : `float` 

1043 Size of the image pixels in physical units 

1044 referencePixel : `lsst.geom.Point2D` 

1045 Pixel for origin of WCS 

1046 

1047 Returns 

1048 ------- 

1049 `lsst.afw.geom.Wcs` 

1050 Wcs object for mapping between pixels and focal plane. 

1051 """ 

1052 md = dafBase.PropertySet() 

1053 if referencePixel is None: 

1054 referencePixel = lsst.geom.PointD(0, 0) 

1055 for i in range(2): 

1056 md.set("CRPIX%d"%(i + 1), referencePixel[i]) 

1057 md.set("CRVAL%d"%(i + 1), 0.) 

1058 md.set("CDELT1", pixelSize[0]) 

1059 md.set("CDELT2", pixelSize[1]) 

1060 md.set("CTYPE1", "CAMERA_X") 

1061 md.set("CTYPE2", "CAMERA_Y") 

1062 md.set("CUNIT1", "mm") 

1063 md.set("CUNIT2", "mm") 

1064 

1065 return afwGeom.makeSkyWcs(md) 

1066 

1067 

1068def findAmp(ccd, pixelPosition): 

1069 """Find the Amp with the specified pixel position within the composite 

1070 

1071 Parameters 

1072 ---------- 

1073 ccd : `lsst.afw.cameraGeom.Detector` 

1074 Detector to look in. 

1075 pixelPosition : `lsst.geom.Point2I` 

1076 The pixel position to find the amp for. 

1077 

1078 Returns 

1079 ------- 

1080 `lsst.afw.table.AmpInfoCatalog` 

1081 Amp record in which ``pixelPosition`` falls or `None` if no Amp found. 

1082 """ 

1083 for amp in ccd: 

1084 if amp.getBBox().contains(pixelPosition): 

1085 return amp 

1086 

1087 return None