Coverage for python/lsst/images/tests/_minify_for_fixtures.py: 20%

187 statements  

« prev     ^ index     » next       coverage.py v7.16.1, created at 2026-09-22 10:42 +0000

1# This file is part of lsst-images. 

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# Use of this source code is governed by a 3-clause BSD-style 

10# license that can be found in the LICENSE file. 

11 

12"""Minify a real on-disk archive into a small JSON test fixture. 

13 

14Reads a FITS or NDF file via the appropriate input archive, takes a 

15small subset of the in-memory object, and writes JSON via 

16``JsonOutputArchive``. Used to populate the ``as_shipped`` and ``dp1`` / 

17``dp2`` fixture variants under ``tests/data/schemas/`` with 

18derived-from-real test data that exercises the full read path. 

19 

20Per top-level type the subset rule is: 

21 

22 VisitImage Crop the image/mask/variance planes to a small (~16x16) 

23 corner, keeping the real single-instance structures (PSF 

24 such as Piff, detector frames) that synthetic fixtures 

25 cannot reproduce -- the whole point of deriving a fixture 

26 from real data. Homogeneous repeated collections (detector 

27 amplifiers, aperture-correction entries) are trimmed to a 

28 representative few, since one entry exercises the schema as 

29 well as sixteen. The projection's pixel->sky mapping is 

30 replaced by its linear (affine) approximation over the kept 

31 box: a real TAN-SIP WCS serializes as a ~100 KB AST 

32 polynomial dump, but over a 16x16 box it is linear to far 

33 below a pixel, so the affine form is schema-identical and 

34 orders of magnitude smaller. A Piff PSF's field 

35 interpolation is truncated to a low order (the order-4 

36 solution table is ~225 KB; order 0, the field-averaged PSF, 

37 is schema-identical and ~13x smaller). 

38 

39 CellCoadd Crop to a small block of cells (preferring a block that 

40 includes a missing cell so the sparse-grid path is 

41 exercised) and then *morph* that block onto a tiny cell 

42 grid: each cell's planes are decimated from the native 

43 cell size down to a few pixels and re-stitched, and the 

44 PSF kernels are cropped to a small odd window. The grid 

45 topology (number of cells, the missing-cell set, band, 

46 mask schema and provenance shape) is preserved; the pixel 

47 values and WCS are *not* physically meaningful. This is 

48 the "morph cells in place" fallback: it sidesteps the 

49 outer-ring problem (inputs/PSFs that overlap kept cells) 

50 by rebuilding a self-consistent miniature coadd rather 

51 than trying to carve an accurate subset out of the real 

52 one. An accurate per-cell subset would inline several 

53 150x150 planes per cell and produce multi-megabyte JSON, 

54 which defeats the purpose of a fixture. 

55 

56Run interactively (CellCoadd works with just this package installed; 

57VisitImage needs a full Rubin environment so the real PSF can be read):: 

58 

59 python -c " 

60 from lsst.images.tests._minify_for_fixtures import minify 

61 schemas = 'tests/data/schemas' 

62 minify( 

63 'cell_example.fits', 

64 f'{schemas}/cell_coadd/cell_coadd-1.0.0-as_shipped.json', 

65 ) 

66 minify('dp1.fits', f'{schemas}/visit_image/visit_image-1.0.0.dev-dp1.json') 

67 minify('dp2.fits', f'{schemas}/visit_image/visit_image-1.0.0.dev-dp2.json') 

68 " 

69 

70The helper is invoked manually by developers when they have a real 

71on-disk file to derive from; it is not exercised by CI. 

72""" 

73 

74from __future__ import annotations 

75 

76__all__ = ("minify",) 

77 

78import os 

79from collections.abc import Callable 

80from typing import Any, cast 

81 

82import numpy as np 

83 

84from .. import DifferenceImage, VisitImage 

85from .._cell_grid import CellGrid, CellGridBounds, CellIJ, PatchDefinition 

86from .._geom import YX, Box 

87from .._image import Image 

88from .._mask import Mask 

89from .._transforms import SkyProjection, TractFrame, Transform 

90from .._transforms._ast import PolyMap 

91from ..cells import CellCoadd, CellField, CellPointSpreadFunction, CoaddProvenance 

92from ..convolution_kernels import ImageBasisConvolutionKernel 

93from ..psfs import PiffWrapper 

94from ..serialization import read_archive 

95from ._creation import make_random_sky_projection 

96 

97# Default morph parameters for CellCoadd. ``CELL_SIZE`` should divide the 

98# native cell size evenly; ``KERNEL_SIZE`` must be odd. ``MAX_INPUTS`` caps 

99# the provenance ``inputs`` table (a real coadd has hundreds of visits); the 

100# full provenance schema is already exercised by the ``coadd_provenance`` 

101# fixture, so here we keep just enough rows to be representative. 

102_CELL_SIZE = 6 

103_KERNEL_SIZE = 5 

104_MAX_INPUTS = 6 

105 

106# Default trim parameters for VisitImage. Amplifiers and aperture-correction 

107# entries are homogeneous collections, so a couple of each cover the schema 

108# just as well as the full set (a real detector has 16 amplifiers and dozens 

109# of aperture corrections). 

110_MAX_AMPLIFIERS = 2 

111_MAX_APERTURE_CORRECTIONS = 2 

112 

113# Field-interpolation order to truncate a Piff PSF to (the solution table of a 

114# real order-4 PixelGrid PSF dominates the fixture at ~225 KB). Order 0 is the 

115# field-averaged PSF; set to `None` to leave the PSF untouched. 

116_PSF_INTERP_ORDER = 0 

117 

118# Maximum permitted deviation (radians) when approximating a projection's 

119# pixel->sky mapping with an affine one. Over a fixture's tiny box the real 

120# mapping is linear well below this, so the fit always succeeds. 

121_PROJECTION_LINEAR_APPROX_TOL = 1e-8 

122 

123 

124def minify(in_path: str, out_path: str) -> None: 

125 """Read a real archive at ``in_path``, take a small subset, and write JSON. 

126 

127 Parameters 

128 ---------- 

129 in_path 

130 Path to a FITS (``.fits`` / ``.fits.gz``) or NDF (``.sdf`` / ``.ndf``) 

131 file to read. 

132 out_path 

133 Path to the output file to write. The parent directory is 

134 created if it does not exist. Can be any supported file format. 

135 File extension controls the output format. 

136 

137 Raises 

138 ------ 

139 ValueError 

140 If the file extension is not recognised. 

141 NotImplementedError 

142 If the top-level type is not one this helper knows how to subset. 

143 """ 

144 obj = read_archive(in_path) 

145 subsetter = _dispatch(obj) 

146 subset = subsetter(obj) 

147 

148 os.makedirs(os.path.dirname(os.path.abspath(out_path)), exist_ok=True) 

149 subset.write(out_path) 

150 

151 

152def _dispatch(image: Any) -> Callable[[Any], Any]: 

153 """Return the relevant subsetter for this image object.""" 

154 match image: 

155 case DifferenceImage(): # this branch needs to go first, as DifferenceImage subclasses VisitImage 155 ↛ 156line 155 didn't jump to line 156 because the pattern on line 155 never matched

156 return _subset_difference_image 

157 case VisitImage(): 157 ↛ 158line 157 didn't jump to line 158 because the pattern on line 157 never matched

158 return _subset_visit_image 

159 case CellCoadd(): 

160 return _subset_cell_coadd 

161 case _: 

162 raise NotImplementedError(f"No minify rule for image of type {type(image)}.") 

163 

164 

165# -- VisitImage ------------------------------------------------------------ 

166 

167 

168def _subset_visit_image[T: VisitImage]( 

169 visit_like_image: T, 

170 *, 

171 size: int = 16, 

172 max_amplifiers: int = _MAX_AMPLIFIERS, 

173 max_aperture_corrections: int = _MAX_APERTURE_CORRECTIONS, 

174 linearize_projection: bool = True, 

175 projection_tol: float = _PROJECTION_LINEAR_APPROX_TOL, 

176 psf_interp_order: int | None = _PSF_INTERP_ORDER, 

177) -> T: 

178 """Crop a VisitImage's pixel planes to a small corner and trim its 

179 homogeneous collections. 

180 

181 The detector frames are a single structure carried through unchanged by 

182 ``__getitem__``. The detector's amplifiers and the aperture-correction map 

183 are repeated, schema-identical entries, so they are trimmed to a 

184 representative few. The projection's pixel->sky mapping is replaced by its 

185 affine approximation over the kept box (see ``_linear_approx_projection``) 

186 unless ``linearize_projection`` is false. A Piff PSF's field interpolation 

187 is truncated to ``psf_interp_order`` (see ``_simplify_piff_psf``) unless 

188 that is `None`. 

189 """ 

190 bbox = visit_like_image.bbox 

191 y0 = bbox.y.start 

192 x0 = bbox.x.start 

193 y1 = min(y0 + size, bbox.y.stop) 

194 x1 = min(x0 + size, bbox.x.stop) 

195 subset = cast(T, visit_like_image[Box.factory[y0:y1, x0:x1]]) 

196 

197 # ``subset`` is a fresh throwaway object whose detector amplifier list, 

198 # aperture-correction map and PSF are live, mutable components. Trim them 

199 # in place through the public accessors rather than reaching for private 

200 # attributes. 

201 del subset.detector.amplifiers[max_amplifiers:] 

202 aperture_corrections = subset.aperture_corrections 

203 for key in list(aperture_corrections)[max_aperture_corrections:]: 

204 del aperture_corrections[key] 

205 if psf_interp_order is not None and isinstance(subset.psf, PiffWrapper): 

206 _simplify_piff_psf(subset.psf, order=psf_interp_order) 

207 

208 if not linearize_projection or subset.sky_projection is None: 

209 return subset 

210 

211 # The pixel planes carry the projection immutably (there is no public 

212 # setter for it), so install the affine approximation by rebuilding the 

213 # VisitImage from its public components with re-viewed planes. Only the 

214 # image plane's projection is actually serialized, but keeping all three 

215 # consistent avoids surprises. 

216 linear = _linear_approx_projection(subset.sky_projection, subset.image.bbox, tol=projection_tol) 

217 return type(visit_like_image)( 

218 subset.image.view(sky_projection=linear), 

219 mask=subset.mask.view(sky_projection=linear), 

220 variance=subset.variance.view(sky_projection=linear), 

221 sky_projection=linear, 

222 psf=subset.psf, 

223 obs_info=subset.obs_info, 

224 bounds=subset.bounds, 

225 summary_stats=subset.summary_stats, 

226 detector=subset.detector, 

227 photometric_scaling=subset.photometric_scaling, 

228 aperture_corrections=subset.aperture_corrections, 

229 backgrounds=subset.backgrounds, 

230 band=subset.band, 

231 metadata=subset.metadata, 

232 ) 

233 

234 

235def _linear_approx_projection(sky_projection: SkyProjection, bbox: Box, *, tol: float) -> SkyProjection: 

236 """Return a copy of ``sky_projection`` whose pixel->sky mapping is replaced 

237 by its best linear (affine) approximation over ``bbox``. 

238 

239 Real WCS mappings (e.g. TAN-SIP) serialize as large AST polynomial dumps. 

240 Over the small box of a fixture they are linear to far below a pixel, so 

241 an affine approximation is schema-identical but orders of magnitude 

242 smaller. The result carries no FITS approximation (the affine is itself 

243 trivially FITS-representable). 

244 

245 This is written as a self-contained ``sky_projection -> sky_projection`` 

246 transform so it can be promoted to a public 

247 ``SkyProjection.linear_approx(bbox, tol)`` method later with essentially 

248 no change. It assumes a 2-D pixel->sky 

249 mapping. 

250 

251 Parameters 

252 ---------- 

253 sky_projection 

254 The projection to approximate. 

255 bbox 

256 Box (in pixel coordinates) over which the approximation must hold. 

257 tol 

258 Maximum permitted deviation from linearity, as a Cartesian 

259 displacement in the output (sky, radians) coordinates. AST raises 

260 ``RuntimeError`` if no fit within ``tol`` exists. 

261 """ 

262 transform = sky_projection.pixel_to_sky_transform 

263 mapping = transform._ast_mapping 

264 lbnd = [bbox.x.start, bbox.y.start] 

265 ubnd = [bbox.x.stop, bbox.y.stop] 

266 # linearApprox yields [offsets; Jacobian] as a (1 + n_out, n_in) array on 

267 # both AST backends (astshim returns the flat buffer in the same order, so 

268 # the reshape recovers the same layout the starlink-pyast bridge returns). 

269 fit = np.asarray(mapping.linearApprox(lbnd, ubnd, tol), dtype=float).reshape(3, 2) 

270 offset = fit[0] # (lon0, lat0), radians 

271 jacobian = fit[1:] # jacobian[i, j] = d(out_i) / d(in_j), in = (x, y) 

272 jacobian_inv = np.linalg.inv(jacobian) 

273 forward = _affine_polymap_coeffs(jacobian, offset) 

274 inverse = _affine_polymap_coeffs(jacobian_inv, -jacobian_inv @ offset) 

275 affine = Transform( 

276 transform.in_frame, 

277 transform.out_frame, 

278 PolyMap(forward, inverse), 

279 in_bounds=sky_projection.pixel_bounds, 

280 ) 

281 return SkyProjection(affine) 

282 

283 

284def _affine_polymap_coeffs(matrix: np.ndarray, offset: np.ndarray) -> np.ndarray: 

285 """Build AST ``PolyMap`` coefficients for ``out = matrix @ in + offset``. 

286 

287 Each row is ``[coefficient, output_axis (1-based), power_of_in_1, ...]``; 

288 one constant row plus one row per input per output axis. Returned as a 

289 float array, which is the form both AST backends require. 

290 """ 

291 n = len(offset) 

292 coeffs: list[list[float]] = [] 

293 for i in range(n): 

294 coeffs.append([float(offset[i]), i + 1, *([0] * n)]) 

295 for j in range(n): 

296 powers = [1 if k == j else 0 for k in range(n)] 

297 coeffs.append([float(matrix[i][j]), i + 1, *powers]) 

298 return np.array(coeffs, dtype=float) 

299 

300 

301def _simplify_piff_psf(psf: PiffWrapper, *, order: int) -> None: 

302 """Truncate a Piff PSF's field interpolation to ``order``, in place. 

303 

304 A real Piff PSF interpolates a per-pixel model across the focal plane with 

305 a high-order 2-D polynomial; that solution table dominates the serialized 

306 size (a 25x25 PixelGrid x order-4 polynomial is ~225 KB). Truncating to 

307 ``order`` keeps only the lowest-order field terms -- order 0 is the 

308 field-averaged PSF -- which is schema-identical but far smaller, and needs 

309 no stars or refit (the fitted ``stars`` are already dropped on serialize). 

310 

311 Only ``BasisPolynomial``-interpolated PSFs are handled; anything else (a 

312 higher-order model already at/under ``order``, a non-polynomial interp) is 

313 left untouched. 

314 

315 ``piff`` is imported lazily because it is an optional dependency; this is 

316 only ever reached when the PSF being simplified is itself a Piff PSF. 

317 """ 

318 interp = getattr(psf.piff_psf, "interp", None) 

319 if interp is None or type(interp).__name__ != "BasisPolynomial" or interp.q is None: 

320 return 

321 if order >= max(interp._orders): 

322 return 

323 

324 from piff import BasisPolynomial 

325 

326 # ``q`` has one column per active basis term; the terms are the True cells 

327 # of ``_mask`` in row-major (i, j) order (see BasisPolynomial.basis). Make 

328 # the same ordering for a lower-order interp and copy the shared columns. 

329 def _terms(orders: tuple[int, ...], mask: np.ndarray) -> list[tuple[int, ...]]: 

330 grids = np.meshgrid(*[np.arange(o + 1) for o in orders], indexing="ij") 

331 return list(zip(*(grid[mask].tolist() for grid in grids))) 

332 

333 old_terms = _terms(interp._orders, interp._mask) 

334 truncated = BasisPolynomial(order, keys=list(interp._keys)) 

335 new_terms = _terms(truncated._orders, truncated._mask) 

336 column_of = {term: index for index, term in enumerate(old_terms)} 

337 truncated.q = np.ascontiguousarray(interp.q[:, [column_of[term] for term in new_terms]]) 

338 psf.piff_psf.interp = truncated 

339 

340 

341# -- DifferenceImage ------------------------------------------------------- 

342 

343 

344def _subset_difference_image( 

345 difference_image: DifferenceImage, 

346 *, 

347 size: int = 16, 

348 max_amplifiers: int = _MAX_AMPLIFIERS, 

349 max_aperture_corrections: int = _MAX_APERTURE_CORRECTIONS, 

350 linearize_projection: bool = True, 

351 projection_tol: float = _PROJECTION_LINEAR_APPROX_TOL, 

352 psf_interp_order: int | None = _PSF_INTERP_ORDER, 

353) -> DifferenceImage: 

354 """Shrink a difference image. 

355 

356 Most of the shrinking is delegated to `_subset_visit_image`. 

357 

358 Template provenance is shrunk to the first few entries. 

359 

360 Difference kernel basis images are sliced to the innermost pixels and the 

361 number of basis functions is shrunk to the first few. 

362 """ 

363 result = _subset_visit_image( 

364 difference_image, 

365 size=size, 

366 max_amplifiers=max_amplifiers, 

367 max_aperture_corrections=max_aperture_corrections, 

368 linearize_projection=linearize_projection, 

369 projection_tol=projection_tol, 

370 psf_interp_order=psf_interp_order, 

371 ) 

372 if difference_image._kernel is not None: 

373 result.kernel = _subset_difference_kernel(cast(ImageBasisConvolutionKernel, difference_image._kernel)) 

374 result.templates = difference_image.templates[:2] if difference_image.templates is not None else None 

375 return result 

376 

377 

378def _subset_difference_kernel( 

379 kernel: ImageBasisConvolutionKernel, 

380 *, 

381 n_basis_images: int = 2, 

382 basis_radius: int = 1, 

383) -> ImageBasisConvolutionKernel: 

384 kernel_bbox = kernel.kernel_bbox.absolute[ 

385 -basis_radius : basis_radius + 1, -basis_radius : basis_radius + 1 

386 ] 

387 slices = kernel_bbox.slice_within(kernel.kernel_bbox) 

388 return ImageBasisConvolutionKernel( 

389 kernel.basis[:n_basis_images, *slices], 

390 kernel.spatial[:n_basis_images], 

391 ) 

392 

393 

394# -- CellCoadd ------------------------------------------------------------- 

395 

396 

397def _subset_cell_coadd( 

398 cell_coadd: CellCoadd, 

399 *, 

400 cell_size: int = _CELL_SIZE, 

401 kernel_size: int = _KERNEL_SIZE, 

402 max_inputs: int = _MAX_INPUTS, 

403) -> CellCoadd: 

404 """Crop a CellCoadd to a small block of cells and morph it onto a tiny 

405 grid (see the module docstring for the rationale). 

406 """ 

407 if kernel_size % 2 == 0: 

408 raise ValueError(f"kernel_size must be odd, got {kernel_size}.") 

409 

410 # 1. Pick a block of (up to) 2x2 cells, preferring one that contains a 

411 # missing cell so the sparse-grid path is exercised. Falls back to the 

412 # first available block when the coadd is fully dense. 

413 block = cell_coadd[_choose_block_bbox(cell_coadd)] 

414 

415 grid = block.grid 

416 cs = grid.cell_shape 

417 start = block.bounds.subgrid_start 

418 stop = block.bounds.subgrid_stop 

419 n_i = stop.i - start.i 

420 n_j = stop.j - start.j 

421 

422 # 2. Build a tiny full-patch grid with the same cell *count* as the 

423 # original patch but ``cell_size`` pixels per cell, anchored at (0, 0). 

424 full_shape = grid.grid_size 

425 new_grid = CellGrid( 

426 bbox=Box.factory[0 : full_shape.i * cell_size, 0 : full_shape.j * cell_size], 

427 cell_shape=YX(y=cell_size, x=cell_size), 

428 ) 

429 new_block_bbox = _scale_box_to_grid(block.bbox, grid, cell_size) 

430 new_bounds = CellGridBounds(grid=new_grid, bbox=new_block_bbox, missing=block.bounds.missing) 

431 

432 # 3. Decimate each plane. Because the block's planes tile the kept cells 

433 # contiguously, a uniform stride that maps one native cell onto 

434 # ``cell_size`` samples is equivalent to per-cell decimation. 

435 step_y = max(1, cs.y // cell_size) 

436 step_x = max(1, cs.x // cell_size) 

437 ny = n_i * cell_size 

438 nx = n_j * cell_size 

439 

440 def shrink2d(array: np.ndarray) -> np.ndarray: 

441 return np.ascontiguousarray(array[::step_y, ::step_x][:ny, :nx]) 

442 

443 def shrink3d(array: np.ndarray) -> np.ndarray: 

444 return np.ascontiguousarray(array[::step_y, ::step_x, :][:ny, :nx, :]) 

445 

446 # 4. Synthetic-but-valid sky_projection over the tiny tract frame. 

447 rng = np.random.default_rng(0) 

448 tract_frame = TractFrame(skymap=cell_coadd.skymap, tract=cell_coadd.tract, bbox=new_grid.bbox) 

449 sky_projection = make_random_sky_projection(rng, tract_frame, new_block_bbox) 

450 

451 unit = cell_coadd.unit 

452 image = Image(shrink2d(block.image.array), bbox=new_block_bbox, unit=unit, sky_projection=sky_projection) 

453 mask = Mask(shrink3d(block.mask.array), schema=block.mask.schema, bbox=new_block_bbox) 

454 variance = Image(shrink2d(block.variance.array), bbox=new_block_bbox, unit=unit**2) 

455 mask_fractions = { 

456 name: Image(shrink2d(plane.array), bbox=new_block_bbox) 

457 for name, plane in block.mask_fractions.items() 

458 } 

459 noise_realizations = [ 

460 Image(shrink2d(plane.array), bbox=new_block_bbox) for plane in block.noise_realizations 

461 ] 

462 

463 # 5. Crop the PSF kernels to a small odd window about their centre, 

464 # keeping the (n_i, n_j) per-cell structure and NaN-for-missing cells. 

465 psf_array = block.psf._array 

466 ky, kx = psf_array.shape[2:] 

467 half = kernel_size // 2 

468 cy, cx = ky // 2, kx // 2 

469 psf_array = np.ascontiguousarray(psf_array[:, :, cy - half : cy + half + 1, cx - half : cx + half + 1]) 

470 psf = CellPointSpreadFunction(psf_array, bounds=new_bounds) 

471 

472 # 6. Patch geometry scaled onto the tiny grid; provenance and backgrounds 

473 # are reused as-is (provenance is cell-indexed and already subset). 

474 patch = PatchDefinition( 

475 id=block.patch.id, 

476 index=block.patch.index, 

477 inner_bbox=_scale_box_to_grid(block.patch.inner_bbox, grid, cell_size), 

478 cells=new_grid, 

479 ) 

480 

481 provenance = block._provenance 

482 if provenance is not None: 

483 provenance = _trim_provenance(provenance, max_inputs=max_inputs) 

484 

485 # Aperture corrections are not subset when CellCoadd is subset with a 

486 # bounding box, because they're always tiny. But that makes setting 

487 # up a consistent new grid for them tricky. 

488 aperture_corrections = {} 

489 new_apcorr_bounds = None 

490 for i, (name, field) in enumerate(block.aperture_corrections.items()): 

491 if new_apcorr_bounds is None: 

492 new_apcorr_bounds = CellGridBounds( 

493 grid=new_grid, 

494 bbox=_scale_box_to_grid(field.bounds.bbox, grid, cell_size), 

495 missing=cell_coadd.bounds.missing, 

496 ) 

497 aperture_corrections[name] = CellField(new_apcorr_bounds, field._array) 

498 if i >= 2: 

499 break 

500 

501 return CellCoadd( 

502 image, 

503 mask=mask, 

504 variance=variance, 

505 mask_fractions=mask_fractions, 

506 noise_realizations=noise_realizations, 

507 sky_projection=sky_projection, 

508 band=block.band, 

509 psf=psf, 

510 patch=patch, 

511 provenance=provenance, 

512 backgrounds=block._backgrounds, 

513 aperture_corrections=aperture_corrections, 

514 ) 

515 

516 

517def _trim_provenance(provenance: CoaddProvenance, *, max_inputs: int) -> CoaddProvenance: 

518 """Cap the provenance ``inputs`` table to ``max_inputs`` rows and drop any 

519 contributions that reference the removed inputs. 

520 

521 The two-table structure, polygon arrays and string dictionary-compression 

522 paths are all preserved; only the number of contributing visits shrinks. 

523 """ 

524 inputs = provenance.inputs 

525 if len(inputs) <= max_inputs: 

526 return provenance 

527 kept_inputs = inputs[:max_inputs] 

528 keys = {(str(row["instrument"]), int(row["visit"]), int(row["detector"])) for row in kept_inputs} 

529 contributions = provenance.contributions 

530 mask = np.array( 

531 [ 

532 (str(instrument), int(visit), int(detector)) in keys 

533 for instrument, visit, detector in zip( 

534 contributions["instrument"], contributions["visit"], contributions["detector"] 

535 ) 

536 ], 

537 dtype=bool, 

538 ) 

539 return CoaddProvenance(inputs=kept_inputs, contributions=contributions[mask]) 

540 

541 

542def _choose_block_bbox(cell_coadd: CellCoadd) -> Box: 

543 """Return the pixel bbox of a (up to) 2x2 block of cells to keep. 

544 

545 Prefers a block containing a missing cell; otherwise the block anchored at 

546 the start of the populated region. Never raises if there is no missing 

547 cell. 

548 """ 

549 bounds = cell_coadd.bounds 

550 grid = bounds.grid 

551 start = bounds.subgrid_start 

552 stop = bounds.subgrid_stop 

553 span_i = min(2, stop.i - start.i) 

554 span_j = min(2, stop.j - start.j) 

555 

556 target = next(iter(sorted(bounds.missing)), None) 

557 if target is not None: 

558 # Anchor the block so it includes the missing cell, clamped to the 

559 # populated index range. 

560 i0 = min(max(target.i, start.i), stop.i - span_i) 

561 j0 = min(max(target.j, start.j), stop.j - span_j) 

562 else: 

563 i0 = start.i 

564 j0 = start.j 

565 

566 lo = grid.bbox_of(CellIJ(i=i0, j=j0)) 

567 hi = grid.bbox_of(CellIJ(i=i0 + span_i - 1, j=j0 + span_j - 1)) 

568 return Box.factory[lo.y.start : hi.y.stop, lo.x.start : hi.x.stop] 

569 

570 

571def _scale_box_to_grid(box: Box, grid: CellGrid, cell_size: int) -> Box: 

572 """Map a grid-aligned box onto a grid with ``cell_size`` pixels per cell, 

573 anchored at the origin. 

574 """ 

575 cs = grid.cell_shape 

576 s = grid.bbox.start 

577 iy0 = (box.y.start - s.y) // cs.y 

578 iy1 = (box.y.stop - s.y) // cs.y 

579 ix0 = (box.x.start - s.x) // cs.x 

580 ix1 = (box.x.stop - s.x) // cs.x 

581 return Box.factory[iy0 * cell_size : iy1 * cell_size, ix0 * cell_size : ix1 * cell_size]