Coverage for python/lsst/images/_mask.py: 73%
533 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-13 09:31 +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.
12from __future__ import annotations
14__all__ = (
15 "Mask",
16 "MaskPlane",
17 "MaskPlaneBit",
18 "MaskSchema",
19 "MaskSerializationModel",
20 "get_legacy_deep_coadd_mask_planes",
21 "get_legacy_difference_image_mask_planes",
22 "get_legacy_non_cell_coadd_mask_planes",
23 "get_legacy_optional_mask_planes",
24 "get_legacy_template_mask_planes",
25 "get_legacy_visit_image_mask_planes",
26)
28import dataclasses
29import math
30from collections.abc import Callable, Iterable, Iterator, Mapping, Sequence, Set
31from types import EllipsisType
32from typing import TYPE_CHECKING, Any, ClassVar, cast
34import astropy.io.fits
35import astropy.wcs
36import numpy as np
37import numpy.typing as npt
38import pydantic
40from lsst.resources import ResourcePath, ResourcePathExpression
42from . import fits
43from ._generalized_image import GeneralizedImage
44from ._geom import YX, Box, NoOverlapError
45from ._transforms import Frame, SkyProjection, SkyProjectionSerializationModel
46from .describe import DescribableMixin, DescribeOptions, FieldRole, Report, ReportField, ReportTable
47from .serialization import (
48 ArchiveReadError,
49 ArchiveTree,
50 ArrayReferenceModel,
51 InlineArrayModel,
52 InputArchive,
53 IntegerType,
54 InvalidParameterError,
55 MetadataValue,
56 NumberType,
57 OutputArchive,
58 is_integer,
59 no_header_updates,
60)
61from .utils import is_none
63if TYPE_CHECKING:
64 try:
65 from lsst.afw.image import Mask as LegacyMask
66 except ImportError:
67 type LegacyMask = Any # type: ignore[no-redef]
70@dataclasses.dataclass(frozen=True)
71class MaskPlane:
72 """Name and description of a single plane in a mask array."""
74 name: str
75 """Unique name for the mask plane (`str`)."""
77 description: str
78 """Human-readable documentation for the mask plane (`str`)."""
80 @classmethod
81 def read_legacy(cls, header: astropy.io.fits.Header, *, strip: bool = True) -> dict[str, int]:
82 """Read mask plane descriptions written by
83 `lsst.afw.image.Mask.writeFits`.
85 Parameters
86 ----------
87 header
88 FITS header.
89 strip
90 If `True` (default), delete the ``MP_`` cards from the header after
91 reading them, as appropriate when the mask is being reinterpreted
92 for new code only. If `False`, leave them in place so they can be
93 propagated for backwards compatibility (re-indexed to the new
94 schema by the caller).
96 Returns
97 -------
98 `dict` [`str`, `int`]
99 A dictionary mapping mask plane name to integer bit index.
100 """
101 result: dict[str, int] = {}
102 for card in list(header.cards):
103 if card.keyword.startswith("MP_"):
104 result[card.keyword.removeprefix("MP_")] = card.value
105 if strip:
106 del header[card.keyword]
107 return result
110@dataclasses.dataclass(frozen=True)
111class MaskPlaneBit:
112 """The nested array index and mask value associated with a single mask
113 plane.
114 """
116 index: int
117 """Index into the last dimension of the mask array where this plane's bit
118 is stored.
119 """
121 mask: np.integer
122 """Bitmask that selects just this plane's bit from a mask array value
123 (`numpy.integer`).
124 """
126 @classmethod
127 def compute(cls, overall_index: int, stride: int, mask_type: type[np.integer]) -> MaskPlaneBit:
128 """Construct a `MaskPlaneBit` from the overall index of a plane in a
129 `MaskSchema` and the stride (number of bits per mask array element).
131 Parameters
132 ----------
133 overall_index
134 Index of the plane across the whole schema.
135 stride
136 Number of mask bits per array element.
137 mask_type
138 Integer dtype of the mask array elements.
139 """
140 index, bit = divmod(overall_index, stride)
141 return cls(index, mask_type(1 << bit))
143 def check(self, value: np.ndarray) -> bool:
144 """Test if this bit is set on a single `Mask` pixel value.
146 Parameters
147 ----------
148 value
149 A 1-d array of length `MaskSchema.mask_size`, representing a
150 single pixel in a `Mask`.
151 """
152 return bool(value[self.index] & self.mask)
155class MaskSchema(DescribableMixin):
156 """A schema for a bit-packed mask array.
158 Parameters
159 ----------
160 planes
161 Iterable of `MaskPlane` instances that define the schema. `None`
162 values may be included to reserve bits for future use.
163 dtype
164 The numpy data type of the mask arrays that use this schema.
166 Notes
167 -----
168 A `MaskSchema` is a collection of mask planes, which each correspond to a
169 single bit in a mask array. Mask schemas are immutable and associated with
170 a particular array data type, allowing them to safely precompute the index
171 and bitmask for each plane.
173 `MaskSchema` indexing is by integer (the overall index of a plane in the
174 schema). The `descriptions` attribute may be indexed by plane name to get
175 the description for that plane, and the `bitmask` method can be used to
176 obtain an array that can be used to select one or more planes by name in
177 a mask array that uses this schema.
179 If no mask planes are provided, a `None` placeholder is automatically
180 added.
181 """
183 def __init__(self, planes: Iterable[MaskPlane | None], dtype: npt.DTypeLike = np.uint8) -> None:
184 self._planes: tuple[MaskPlane | None, ...] = tuple(planes) or (None,)
185 self._dtype = cast(np.dtype[np.integer], np.dtype(dtype))
186 stride = self.bits_per_element(self._dtype)
187 self._descriptions = {plane.name: plane.description for plane in self._planes if plane is not None}
188 self._mask_size = math.ceil(len(self._planes) / stride)
189 self._bits: dict[str, MaskPlaneBit] = {
190 plane.name: MaskPlaneBit.compute(n, stride, self._dtype.type)
191 for n, plane in enumerate(self._planes)
192 if plane is not None
193 }
195 @staticmethod
196 def bits_per_element(dtype: npt.DTypeLike) -> int:
197 """Return the number of mask bits per array element for the given
198 data type.
200 Parameters
201 ----------
202 dtype
203 Data type of the mask array elements.
204 """
205 dtype = np.dtype(dtype)
206 match dtype.kind:
207 case "u":
208 return dtype.itemsize * 8
209 case "i":
210 return dtype.itemsize * 8 - 1
211 case _:
212 raise TypeError(f"dtype for masks must be an integer; got {dtype} with kind={dtype.kind}.")
214 def __iter__(self) -> Iterator[MaskPlane | None]:
215 return iter(self._planes)
217 def __len__(self) -> int:
218 return len(self._planes)
220 def __contains__(self, plane: str | MaskPlane) -> bool:
221 return getattr(plane, "name", plane) in self.names
223 def __getitem__(self, i: int) -> MaskPlane | None:
224 return self._planes[i]
226 def _describe(
227 self,
228 options: DescribeOptions = DescribeOptions(),
229 /,
230 *,
231 counts: Mapping[str, int] | None = None,
232 ) -> Report:
233 """Return a `Report` describing this mask schema.
235 Parameters
236 ----------
237 options : `DescribeOptions`, optional
238 Rendering options.
239 counts : `~collections.abc.Mapping` [`str`, `int`], optional
240 Number of set pixels per plane name. When provided, the mask
241 planes table gains a ``"Set pixels"`` column.
242 """
243 # The table lists every plane, so the plane count is REPR_ONLY: repr
244 # needs it to round-trip, but next to the table it would just be noise.
245 fields = [
246 ReportField(
247 label="planes",
248 value=f"<{len(self._planes)} planes>",
249 repr_value=repr(list(self._planes)),
250 positional=True,
251 role=FieldRole.REPR_ONLY,
252 ),
253 ReportField(label="dtype", value=str(self._dtype), repr_value=repr(self._dtype)),
254 ]
255 if options.brief:
256 return Report(type_name="MaskSchema", fields=fields)
257 columns = ["Bit", "Index", "Mask", "Name", "Description"]
258 if counts is not None:
259 columns.append("Set pixels")
260 rows = []
261 for n, plane in enumerate(self._planes):
262 if plane is None:
263 continue
264 row: list[Any] = [
265 n,
266 self._bits[plane.name].index,
267 hex(self._bits[plane.name].mask),
268 plane.name,
269 plane.description,
270 ]
271 if counts is not None:
272 row.append(counts.get(plane.name, 0))
273 rows.append(row)
274 return Report(
275 type_name="MaskSchema",
276 fields=fields,
277 tables=[
278 ReportTable(
279 title="Mask planes",
280 columns=columns,
281 rows=rows,
282 )
283 ],
284 )
286 def __eq__(self, other: object) -> bool:
287 if isinstance(other, MaskSchema): 287 ↛ 289line 287 didn't jump to line 289 because the condition on line 287 was always true
288 return self._planes == other._planes and self._dtype == other._dtype
289 return False
291 @property
292 def dtype(self) -> np.dtype:
293 """The numpy data type of the mask arrays that use this schema."""
294 return self._dtype
296 @property
297 def mask_size(self) -> int:
298 """The number of elements in the last dimension of any mask array that
299 uses this schema.
300 """
301 return self._mask_size
303 @property
304 def names(self) -> Set[str]:
305 """The names of the mask planes, in bit order."""
306 return self._bits.keys()
308 @property
309 def descriptions(self) -> Mapping[str, str]:
310 """A mapping from plane name to description."""
311 return self._descriptions
313 def bit(self, plane: str) -> MaskPlaneBit:
314 """Return the last array index and mask for the given mask plane.
316 Parameters
317 ----------
318 plane
319 Name of the mask plane.
320 """
321 return self._bits[plane]
323 def bitmask(self, *planes: str) -> np.ndarray:
324 """Return a 1-d mask array that represents the union (i.e. bitwise OR)
325 of the planes with the given names.
327 Parameters
328 ----------
329 *planes
330 Mask plane names.
332 Returns
333 -------
334 numpy.ndarray
335 A 1-d array with shape ``(mask_size,)``.
336 """
337 result = np.zeros(self.mask_size, dtype=self._dtype)
338 for plane in planes:
339 bit = self._bits[plane]
340 result[bit.index] |= bit.mask
341 return result
343 def split(self, dtype: npt.DTypeLike) -> list[MaskSchema]:
344 """Split the schema into an equivalent series of schemas that each
345 have a `mask_size` of ``1``, dropping all `None` placeholders.
347 Parameters
348 ----------
349 dtype
350 Data type of the new mask pixels.
352 Returns
353 -------
354 `list` [`MaskSchema`]
355 A list of mask schemas that together include all planes in
356 ``self`` and have `mask_size` equal to ``1``. If there are no
357 mask planes (only `None` placeholders) in ``self``, a single mask
358 schema with a `None` placeholder is returned; otherwise `None`
359 placeholders are returned.
360 """
361 dtype = np.dtype(dtype)
362 planes: list[MaskPlane] = []
363 schemas: list[MaskSchema] = []
364 n_planes_per_schema = self.bits_per_element(dtype)
365 for plane in self._planes:
366 if plane is not None:
367 planes.append(plane)
368 if len(planes) == n_planes_per_schema:
369 schemas.append(MaskSchema(planes, dtype=dtype))
370 planes.clear()
371 if planes: 371 ↛ 373line 371 didn't jump to line 373 because the condition on line 371 was always true
372 schemas.append(MaskSchema(planes, dtype=dtype))
373 if not schemas: 373 ↛ 374line 373 didn't jump to line 374 because the condition on line 373 was never true
374 schemas.append(MaskSchema([None], dtype=dtype))
375 return schemas
377 def update_header(self, header: astropy.io.fits.Header) -> None:
378 """Add a description of this mask schema to a FITS header.
380 Parameters
381 ----------
382 header
383 FITS header to add the mask schema description to.
384 """
385 for n, plane in enumerate(self):
386 if plane is not None:
387 bit = self.bit(plane.name)
388 if bit.index != 0: 388 ↛ 389line 388 didn't jump to line 389 because the condition on line 388 was never true
389 raise TypeError("Only mask schemas with mask_size==1 can be described in FITS.")
390 header.set(f"MSKN{n:04d}", plane.name, f"Name for mask plane {n}.")
391 header.set(f"MSKM{n:04d}", bit.mask, f"Bitmask for plane n={n}; always 1<<n.")
392 # We don't add a comment to the description card, because it's
393 # likely to overrun a single card and get the CONTINUE
394 # treatment. That will cause Astropy to warn about the comment
395 # being truncated and that's worse than just leaving it
396 # unexplained; it's pretty obvious from context what it is.
397 header.set(f"MSKD{n:04d}", plane.description)
399 def strip_header(self, header: astropy.io.fits.Header) -> None:
400 """Remove all header cards added by `update_header`.
402 Parameters
403 ----------
404 header
405 FITS header to remove the mask schema cards from.
406 """
407 for n, plane in enumerate(self):
408 if plane is not None: 408 ↛ 407line 408 didn't jump to line 407 because the condition on line 408 was always true
409 header.remove(f"MSKN{n:04d}", ignore_missing=True)
410 header.remove(f"MSKM{n:04d}", ignore_missing=True)
411 header.remove(f"MSKD{n:04d}", ignore_missing=True)
413 @classmethod
414 def from_fits_header(cls, header: astropy.io.fits.Header, dtype: npt.DTypeLike = np.uint8) -> MaskSchema:
415 """Reconstruct a schema from the ``MSKN``/``MSKD`` cards written by
416 `update_header`.
418 Parameters
419 ----------
420 header
421 FITS header containing ``MSKN{n:04d}`` plane-name cards and
422 ``MSKD{n:04d}`` description cards.
423 dtype
424 Data type of the mask arrays that will use this schema. The cards
425 describe a ``mask_size==1`` serialized form and do not record the
426 in-memory dtype, so the caller must supply it; it defaults to the
427 same ``uint8`` used by the `Mask` constructor.
429 Returns
430 -------
431 `MaskSchema`
432 Schema whose planes are ordered by their ``MSKN`` index, with
433 `None` placeholders inserted for any gaps in that numbering.
435 Raises
436 ------
437 ValueError
438 Raised if the header contains no ``MSKN`` cards.
439 """
440 planes_by_index: dict[int, MaskPlane] = {}
441 for card in header.cards:
442 if card.keyword.startswith("MSKN"):
443 n = int(card.keyword.removeprefix("MSKN"))
444 planes_by_index[n] = MaskPlane(card.value, header.get(f"MSKD{n:04d}", ""))
445 if not planes_by_index:
446 raise ValueError("Header has no MSKN cards describing a mask schema.")
447 planes = [planes_by_index.get(n) for n in range(max(planes_by_index) + 1)]
448 return cls(planes, dtype=dtype)
450 def interpret(self, value: np.ndarray) -> list[str]:
451 """Return the names of the mask planes that are set in the given
452 pixel value.
454 Parameters
455 ----------
456 value
457 A 1-d array of length `mask_size`, representing a single pixel in
458 a `Mask`.
459 """
460 return [name for name, bit in self._bits.items() if bit.check(value)]
463class Mask(GeneralizedImage):
464 """A 2-d bitmask image backed by a 3-d byte array.
466 Parameters
467 ----------
468 array_or_fill
469 Array or fill value for the mask. If a fill value, ``bbox`` or
470 ``shape`` must be provided.
471 schema
472 Schema that defines the planes and their bit assignments.
473 bbox
474 Bounding box for the mask. This sets the shape of the first two
475 dimensions of the array.
476 yx0
477 Logical coordinates of the first pixel in the array, ordered ``y``,
478 ``x`` (unless an `XY` instance is passed). Ignored if
479 ``bbox`` is provided. Defaults to zeros.
480 shape
481 Leading dimensions of the array, ordered ``y``, ``x`` (unless an `XY`
482 instance is passed). Only needed if ``array_or_fill`` is not an
483 array and ``bbox`` is not provided. Like the bbox, this does not
484 include the last dimension of the array.
485 sky_projection
486 Projection that maps the pixel grid to the sky.
487 metadata
488 Arbitrary flexible metadata to associate with the mask.
490 Notes
491 -----
492 Indexing the `array` attribute of a `Mask` does not take into account its
493 ``yx0`` offset, but accessing a subimage mask by indexing a `Mask` with
494 a `Box` does, and the `bbox` of the subimage is set to match its location
495 within the original mask.
497 A mask's ``bbox`` corresponds to the leading dimensions of its backing
498 `numpy.ndarray`, while the last dimension's size is always equal to the
499 `~MaskSchema.mask_size` of its schema, since a schema can in general
500 require multiple array elements to represent all of its planes.
501 """
503 def __init__(
504 self,
505 array_or_fill: np.ndarray | int = 0,
506 /,
507 *,
508 schema: MaskSchema,
509 bbox: Box | None = None,
510 yx0: Sequence[int] | None = None,
511 shape: Sequence[int] | None = None,
512 sky_projection: SkyProjection | None = None,
513 metadata: dict[str, MetadataValue] | None = None,
514 ) -> None:
515 super().__init__(metadata)
516 if shape is not None:
517 shape = tuple(shape)
518 if isinstance(array_or_fill, np.ndarray):
519 array = np.array(array_or_fill, dtype=schema.dtype, copy=None)
520 if array.ndim != 3:
521 raise ValueError("Mask array must be 3-d.")
522 if bbox is None:
523 bbox = Box.from_shape(array.shape[:-1], start=yx0)
524 elif bbox.shape + (schema.mask_size,) != array.shape:
525 raise ValueError(
526 f"Explicit bbox shape {bbox.shape} and schema of size {schema.mask_size} do not "
527 f"match array with shape {array.shape}."
528 )
529 if shape is not None and shape + (schema.mask_size,) != array.shape:
530 raise ValueError(
531 f"Explicit shape {shape} and schema of size {schema.mask_size} do "
532 f"not match array with shape {array.shape}."
533 )
535 else:
536 if bbox is None:
537 if shape is None:
538 raise TypeError("No bbox, size, or array provided.")
539 bbox = Box.from_shape(shape, start=yx0)
540 array = np.full(bbox.shape + (schema.mask_size,), array_or_fill, dtype=schema.dtype)
541 self._array = array
542 self._bbox: Box = bbox
543 self._schema: MaskSchema = schema
544 self._sky_projection = sky_projection
546 @property
547 def array(self) -> np.ndarray:
548 """The low-level array (`numpy.ndarray`).
550 Assigning to this attribute modifies the existing array in place; the
551 bounding box and underlying data pointer are never changed.
552 """
553 return self._array
555 @array.setter
556 def array(self, value: np.ndarray | int) -> None:
557 self._array[:, :] = value
559 @property
560 def schema(self) -> MaskSchema:
561 """Schema that defines the planes and their bit assignments
562 (`MaskSchema`).
563 """
564 return self._schema
566 @property
567 def bbox(self) -> Box:
568 """2-d bounding box of the mask (`Box`).
570 This sets the shape of the first two dimensions of the array.
571 """
572 return self._bbox
574 @property
575 def sky_projection(self) -> SkyProjection[Any] | None:
576 """The projection that maps this mask's pixel grid to the sky
577 (`SkyProjection` | `None`).
579 Notes
580 -----
581 The pixel coordinates used by this projection account for the bounding
582 box ``start`` (i.e. ``yx0``); they are not just array indices.
583 """
584 return self._sky_projection
586 def __getitem__(self, bbox: Box | EllipsisType) -> Mask:
587 bbox, indices = self._handle_getitem_args(bbox)
588 return self._transfer_metadata(
589 Mask(
590 self.array[indices + (slice(None),)],
591 bbox=bbox,
592 schema=self.schema,
593 sky_projection=self._sky_projection,
594 ),
595 bbox=bbox,
596 )
598 def __setitem__(self, bbox: Box | EllipsisType, value: Mask) -> None:
599 subview = self[bbox]
600 subview.clear()
601 subview.update(value)
603 def _describe(self, options: DescribeOptions = DescribeOptions(), /) -> Report:
604 """Return a `Report` describing this mask.
606 Parameters
607 ----------
608 options : `DescribeOptions`, optional
609 Rendering options. ``"bbox"`` and ``"sky_projection"`` are
610 recognized in `DescribeOptions.exclude`, and
611 `DescribeOptions.detail` adds per-plane set-pixel counts to the
612 schema table, which scans the pixel data.
613 """
614 summary = f"Mask({self.bbox!s}, {list(self.schema.names)})"
615 fields = [
616 ReportField(
617 label="array",
618 value="<array>",
619 repr_value="...",
620 positional=True,
621 role=FieldRole.REPR_ONLY,
622 ),
623 ]
624 if "bbox" not in options.exclude:
625 fields.append(ReportField(label="bbox", value=self.bbox, repr_value=repr(self.bbox)))
626 # The schema is rendered as a child below, so repr is the only place
627 # this field needs to appear.
628 fields.append(
629 ReportField(
630 label="schema",
631 value=self.schema,
632 repr_value=repr(self.schema),
633 role=FieldRole.REPR_ONLY,
634 )
635 )
636 if options.brief:
637 return Report(type_name="Mask", summary=summary, fields=fields)
638 child = options.for_child()
639 if options.detail:
640 counts = {name: int(np.count_nonzero(self.get(name))) for name in self.schema.names}
641 schema_report = self.schema._describe(child, counts=counts)
642 else:
643 schema_report = self.schema._describe(child)
644 children: dict[str, Report] = {"schema": schema_report}
645 if "sky_projection" not in options.exclude and self._sky_projection is not None: 645 ↛ 646line 645 didn't jump to line 646 because the condition on line 645 was never true
646 children["sky_projection"] = self._sky_projection._describe(child, bbox=self._bbox)
647 return Report(
648 type_name="Mask",
649 summary=summary,
650 fields=fields,
651 children=children,
652 )
654 def __eq__(self, other: object) -> bool:
655 if not isinstance(other, Mask):
656 return NotImplemented
657 return (
658 self._bbox == other._bbox
659 and self._schema == other._schema
660 and np.array_equal(self._array, other._array, equal_nan=True)
661 )
663 def copy(self) -> Mask:
664 """Deep-copy the mask and metadata."""
665 return self._transfer_metadata(
666 Mask(
667 self._array.copy(), bbox=self._bbox, schema=self._schema, sky_projection=self._sky_projection
668 ),
669 copy=True,
670 )
672 def view(
673 self,
674 *,
675 schema: MaskSchema | EllipsisType = ...,
676 sky_projection: SkyProjection | None | EllipsisType = ...,
677 yx0: Sequence[int] | EllipsisType = ...,
678 ) -> Mask:
679 """Make a view of the mask, with optional updates.
681 Parameters
682 ----------
683 schema
684 Replacement schema; defaults to the current schema.
685 sky_projection
686 Replacement sky projection; defaults to the current one.
687 yx0
688 Replacement origin of the mask; defaults to the current origin.
690 Notes
691 -----
692 This can only be used to make changes to schema descriptions; plane
693 names must remain the same (in the same order).
694 """
695 if schema is ...: 695 ↛ 698line 695 didn't jump to line 698 because the condition on line 695 was always true
696 schema = self._schema
697 else:
698 if list(schema.names) != list(self.schema.names):
699 raise ValueError("Cannot create a mask view with a schema with different names.")
700 if sky_projection is ...: 700 ↛ 701line 700 didn't jump to line 701 because the condition on line 700 was never true
701 sky_projection = self._sky_projection
702 if yx0 is ...: 702 ↛ 704line 702 didn't jump to line 704 because the condition on line 702 was always true
703 yx0 = self._bbox.start
704 return self._transfer_metadata(
705 Mask(self._array, yx0=yx0, schema=schema, sky_projection=sky_projection)
706 )
708 def update(self, other: Mask) -> None:
709 """Update ``self`` to include all common mask values set in ``other``.
711 Parameters
712 ----------
713 other
714 Mask whose set bits are merged into ``self``.
716 Notes
717 -----
718 This only operates on the intersection of the two mask bounding boxes
719 and the mask planes that are present in both. Mask bits are only set,
720 not cleared (i.e. this uses ``|=`` updates, not ``=`` assignments).
721 """
722 lhs = self
723 rhs = other
724 if other.bbox != self.bbox: 724 ↛ 725line 724 didn't jump to line 725 because the condition on line 724 was never true
725 try:
726 bbox = self.bbox.intersection(other.bbox)
727 except NoOverlapError:
728 return
729 lhs = self[bbox]
730 rhs = other[bbox]
731 for name in self.schema.names & other.schema.names:
732 lhs.set(name, rhs.get(name))
734 def get(self, plane: str) -> np.ndarray:
735 """Return a 2-d boolean array for the given mask plane.
737 Parameters
738 ----------
739 plane
740 Name of the mask plane.
742 Returns
743 -------
744 numpy.ndarray
745 A 2-d boolean array with the same shape as `bbox` that is `True`
746 where the bit for ``plane`` is set and `False` elsewhere.
747 """
748 bit = self.schema.bit(plane)
749 return (self._array[..., bit.index] & bit.mask).astype(bool)
751 def compare(self, other: Mask) -> dict[str, tuple[int, int]]:
752 """Return a plane-by-plane comparison with another mask.
754 Parameters
755 ----------
756 other
757 The mask to compare against.
759 Returns
760 -------
761 `dict` [`str`, `tuple` [`int`, `int`]]
762 Dictionary mask planes as keys, where the values are:
764 - the number of pixels with that plane set in ``self`` but not
765 ``other``;
766 - the number of pixels with that plane set in ``other`` but not
767 ``self``.
769 Mask planes where the images have the same pixels set are not
770 included, so the result is empty for identical masks. Mask planes
771 present in only one operand are treated as though they were unset
772 for all pixels in the other operand.
774 Notes
775 -----
776 Two masks with different schemas can have an empty `compare` result
777 while still not satisfying an equality check, as long as the planes
778 they have in common have the same pixels set and any planes present in
779 only one operand have no pixels set.
780 """
781 if self.bbox != other.bbox:
782 raise ValueError("masks must have the same bounding box to compute a difference")
783 empty = np.zeros(self._array.shape[:-1], dtype=bool)
784 result: dict[str, tuple[int, int]] = {}
785 for name in sorted(self.schema.names | other.schema.names):
786 a = self.get(name) if name in self.schema.names else empty
787 b = other.get(name) if name in other.schema.names else empty
788 added = int(np.count_nonzero(a & ~b))
789 removed = int(np.count_nonzero(b & ~a))
790 if added or removed:
791 result[name] = (added, removed)
792 return result
794 def set(self, plane: str, boolean_mask: np.ndarray | EllipsisType = ...) -> None:
795 """Set a mask plane.
797 Parameters
798 ----------
799 plane
800 Name of the mask plane to set.
801 boolean_mask
802 A 2-d boolean array with the same shape as `bbox` that is `True`
803 where the bit for ``plane`` should be set and `False` where it
804 should be left unchanged (*not* set to zero). May be ``...`` to
805 set the bit everywhere.
806 """
807 bit = self.schema.bit(plane)
808 if boolean_mask is not ...: 808 ↛ 810line 808 didn't jump to line 810 because the condition on line 808 was always true
809 boolean_mask = boolean_mask.astype(bool)
810 self._array[boolean_mask, bit.index] |= bit.mask
812 def clear(self, plane: str | None = None, boolean_mask: np.ndarray | EllipsisType = ...) -> None:
813 """Clear one or more mask planes.
815 Parameters
816 ----------
817 plane
818 Name of the mask plane to set. If `None` all mask planes are
819 cleared.
820 boolean_mask
821 A 2-d boolean array with the same shape as `bbox` that is `True`
822 where the bit for ``plane`` should be cleared and `False` where it
823 should be left unchanged. May be ``...`` to clear the bit
824 everywhere.
825 """
826 if boolean_mask is not ...: 826 ↛ 827line 826 didn't jump to line 827 because the condition on line 826 was never true
827 boolean_mask = boolean_mask.astype(bool)
828 if plane is None: 828 ↛ 831line 828 didn't jump to line 831 because the condition on line 828 was always true
829 self._array[boolean_mask, :] = 0
830 else:
831 bit = self.schema.bit(plane)
832 self._array[boolean_mask, bit.index] &= ~bit.mask
834 def add_plane(self, name: str, description: str) -> Mask:
835 """Return a new mask with one additional mask plane.
837 This is a convenience wrapper around `add_planes` for the common case
838 of adding a single plane.
840 Parameters
841 ----------
842 name
843 Unique name for the new mask plane.
844 description
845 Human-readable documentation for the new mask plane.
847 Returns
848 -------
849 `Mask`
850 A new mask whose schema includes the new plane; see `add_planes`
851 for the reallocation and view semantics.
853 Raises
854 ------
855 ValueError
856 Raised if a plane named ``name`` already exists.
857 """
858 return self.add_planes([MaskPlane(name, description)])
860 def add_planes(self, planes: Iterable[MaskPlane | None], *, drop: Iterable[str] = ()) -> Mask:
861 """Return a new mask with planes added and/or dropped.
863 Parameters
864 ----------
865 planes
866 New mask planes to append, in order, after the planes retained
867 from this mask. `None` entries reserve unused bits (placeholders),
868 exactly as in `MaskSchema`.
869 drop
870 Names of existing planes to remove from the schema.
872 Returns
873 -------
874 `Mask`
875 A new mask with the updated schema. Retained planes keep their
876 pixel values (copied by name); newly added planes start cleared.
878 Raises
879 ------
880 ValueError
881 Raised if a name in ``drop`` is not an existing plane, or if a
882 plane in ``planes`` collides with a retained plane name.
884 Notes
885 -----
886 Adding or dropping planes always reallocates the backing array and
887 returns a new `Mask`; this mask is left unchanged and any views or
888 subimages of it continue to refer to the original array with the
889 original schema. This is deliberate: there is no way to update the
890 schema of an existing view, and a stale view must never set bits that
891 its now-outdated schema regards as unused. Dropping a plane compacts
892 the schema, so planes after it are reassigned to lower bits and the
893 pixel values are repacked by plane name to match.
894 """
895 drop_set = set(drop)
896 if unknown := drop_set - set(self._schema.names):
897 raise ValueError(f"Cannot drop mask planes that do not exist: {sorted(unknown)}.")
898 retained = [plane for plane in self._schema if plane is None or plane.name not in drop_set]
899 names = {plane.name for plane in retained if plane is not None}
900 new_planes = list(planes)
901 for plane in new_planes:
902 if plane is None:
903 continue
904 if plane.name in names:
905 raise ValueError(f"Mask plane {plane.name!r} already exists.")
906 names.add(plane.name)
907 new_schema = MaskSchema([*retained, *new_planes], dtype=self._schema.dtype)
908 result = Mask(0, schema=new_schema, bbox=self._bbox, sky_projection=self._sky_projection)
909 # The retained planes are exactly the names common to both schemas, and
910 # ``result`` starts cleared and shares this mask's bbox, so ``update``
911 # transfers their pixel values (and nothing else) by name.
912 result.update(self)
913 return self._transfer_metadata(result, copy=True)
915 def serialize[P: pydantic.BaseModel](
916 self,
917 archive: OutputArchive[P],
918 *,
919 update_header: Callable[[astropy.io.fits.Header], None] = no_header_updates,
920 save_projection: bool = True,
921 add_offset_wcs: str | None = "A",
922 tile_shape: tuple[int, ...] | None = None,
923 options_name: str | None = None,
924 ) -> MaskSerializationModel[P]:
925 """Serialize the mask to an output archive.
927 Parameters
928 ----------
929 archive
930 Archive to write to.
931 update_header
932 A callback that will be given the FITS header for the HDU
933 containing this mask in order to add keys to it. This callback
934 may be provided but will not be called if the output format is not
935 FITS. As multiple HDUs may be added, this function may be called
936 multiple times.
937 save_projection
938 If `True`, save the `SkyProjection` attached to the image, if there
939 is one. This does not affect whether a FITS WCS corresponding to
940 the projection is written (it always is, if available, and if
941 ``add_offset_wcs`` is not ``" "``).
942 add_offset_wcs
943 A FITS WCS single-character suffix to use when adding a linear
944 WCS that maps the FITS array to the logical pixel coordinates
945 defined by ``bbox.start`` / ``yx0``. Set to `None` to not write
946 this WCS. If this is set to ``" "``, it will prevent the
947 `SkyProjection` from being saved as a FITS WCS.
948 tile_shape
949 The recommended shape of each tile, if the archive will save
950 the array in distinct tiles for faster subarray retrieval.
951 This is a hint; archives are not required to use this value.
952 options_name
953 Use this name to look up archive options.
954 """
955 if _archive_prefers_native_mask_arrays(archive):
956 # HDS presents array dimensions in Fortran order, which is the
957 # reverse of the h5py dataset shape. Store the in-memory trailing
958 # mask-byte axis first in HDF5 so Starlink tools see HDS axes
959 # (x, y, byte), without changing the bit packing within a pixel.
960 array_model = archive.add_array(
961 np.moveaxis(self._array, -1, 0),
962 update_header=update_header,
963 tile_shape=tile_shape,
964 options_name=options_name,
965 )
966 if not isinstance(array_model, ArrayReferenceModel): 966 ↛ 967line 966 didn't jump to line 967 because the condition on line 966 was never true
967 raise RuntimeError("Native mask arrays require reference array storage.")
968 array_model.shape = list(self._array.shape)
969 data: list[ArrayReferenceModel | InlineArrayModel] = [array_model]
970 else:
971 data = []
972 for schema_2d in self.schema.split(np.int32):
973 mask_2d = Mask(0, bbox=self.bbox, schema=schema_2d, sky_projection=self._sky_projection)
974 mask_2d.update(self)
975 data.append(
976 mask_2d._serialize_2d(
977 archive,
978 update_header=update_header,
979 add_offset_wcs=add_offset_wcs,
980 tile_shape=tile_shape,
981 options_name=options_name,
982 )
983 )
984 serialized_projection: SkyProjectionSerializationModel[P] | None = None
985 if save_projection and self.sky_projection is not None: 985 ↛ 986line 985 didn't jump to line 986 because the condition on line 985 was never true
986 serialized_projection = archive.serialize_direct("sky_projection", self.sky_projection.serialize)
987 serialized_dtype = NumberType.from_numpy(self.schema.dtype)
988 assert is_integer(serialized_dtype), "Mask dtypes should always be integers."
989 return MaskSerializationModel.model_construct(
990 data=data,
991 yx0=list(self.bbox.start),
992 planes=list(self.schema),
993 dtype=serialized_dtype,
994 sky_projection=serialized_projection,
995 metadata=self._metadata,
996 )
998 def _serialize_2d[P: pydantic.BaseModel](
999 self,
1000 archive: OutputArchive[P],
1001 *,
1002 update_header: Callable[[astropy.io.fits.Header], None] = no_header_updates,
1003 add_offset_wcs: str | None = "A",
1004 tile_shape: tuple[int, ...] | None = None,
1005 options_name: str | None = None,
1006 ) -> ArrayReferenceModel | InlineArrayModel:
1007 def _update_header(header: astropy.io.fits.Header) -> None:
1008 update_header(header)
1009 self.schema.update_header(header)
1010 if self.sky_projection is not None and add_offset_wcs != " ":
1011 if self.fits_wcs:
1012 header.update(self.fits_wcs.to_header(relax=True))
1013 if add_offset_wcs is not None: 1013 ↛ exitline 1013 didn't return from function '_update_header' because the condition on line 1013 was always true
1014 fits.add_offset_wcs(header, x=self.bbox.x.start, y=self.bbox.y.start, key=add_offset_wcs)
1016 assert self.array.shape[2] == 1, "Mask should be split before calling this method."
1017 return archive.add_array(
1018 self._array[:, :, 0],
1019 update_header=_update_header,
1020 tile_shape=tile_shape,
1021 options_name=options_name,
1022 )
1024 @staticmethod
1025 def _get_archive_tree_type[P: pydantic.BaseModel](
1026 pointer_type: type[P],
1027 ) -> type[MaskSerializationModel[P]]:
1028 """Return the serialization model type for this object for an archive
1029 type that uses the given pointer type.
1030 """
1031 return MaskSerializationModel[pointer_type] # type: ignore
1033 _archive_default_name: ClassVar[str] = "mask"
1034 """The name this object should be serialized with when written as the
1035 top-level object.
1036 """
1038 @staticmethod
1039 def from_legacy(
1040 legacy: Any,
1041 plane_map: Mapping[str, MaskPlane] | None = None,
1042 *,
1043 sky_projection: SkyProjection[Any] | None = None,
1044 ) -> Mask:
1045 """Convert from an `lsst.afw.image.Mask` instance.
1047 Parameters
1048 ----------
1049 legacy
1050 An `lsst.afw.image.Mask` instance. This will not share pixel
1051 data with the new object.
1052 plane_map
1053 A mapping from legacy mask plane name to the new plane name and
1054 description. If not provided, the right legacy mask plane will be
1055 guessed, but this can depend on which mask planes the legacy
1056 mask actually has set. Planes from
1057 `get_legacy_optional_mask_planes` are added when the legacy mask
1058 defines them and has pixels set in them.
1059 sky_projection
1060 Projection from pixels to sky.
1062 Returns
1063 -------
1064 `Mask`
1065 Converted mask.
1066 """
1067 return Mask._from_legacy_array(
1068 legacy.array,
1069 legacy.getMaskPlaneDict(),
1070 yx0=YX(y=legacy.getY0(), x=legacy.getX0()),
1071 plane_map=plane_map,
1072 sky_projection=sky_projection,
1073 )
1075 def to_legacy(self, plane_map: Mapping[str, MaskPlane] | None = None) -> Any:
1076 """Convert to an `lsst.afw.image.Mask` instance.
1078 The pixel data will not be shared between the two objects.
1080 Parameters
1081 ----------
1082 plane_map
1083 A mapping from legacy mask plane name to the new plane name and
1084 description. Schema planes the map does not name are written
1085 under their own names. If not provided, every schema plane is
1086 written under its own name.
1088 Returns
1089 -------
1090 `lsst.afw.image.Mask`
1091 Converted mask.
1092 """
1093 import lsst.afw.image
1094 import lsst.geom
1096 result = lsst.afw.image.Mask(self.bbox.to_legacy())
1097 if plane_map is None:
1098 plane_map = {plane.name: plane for plane in self.schema if plane is not None}
1100 # Add a legacy plane to ``result`` and return its bitmask.
1101 def add_legacy_plane(old_name: str) -> int:
1102 return _legacy_bitmask(result.addMaskPlane(old_name))
1104 written: set[str] = set()
1105 for old_name, new_plane in plane_map.items():
1106 old_bitmask = add_legacy_plane(old_name)
1107 if new_plane in self.schema:
1108 result.array[self.get(new_plane.name)] |= old_bitmask
1109 written.add(new_plane.name)
1110 # Write schema planes the map does not name under their own names.
1111 for plane in self.schema:
1112 if plane is None or plane.name in written:
1113 continue
1114 result.array[self.get(plane.name)] |= add_legacy_plane(plane.name)
1115 return result
1117 @staticmethod
1118 def _from_legacy_array(
1119 array2d: np.ndarray,
1120 old_planes: Mapping[str, int],
1121 *,
1122 yx0: YX[int],
1123 plane_map: Mapping[str, MaskPlane] | None = None,
1124 sky_projection: SkyProjection | None = None,
1125 ) -> Mask:
1126 plane_map = _resolve_legacy_plane_map(plane_map, old_planes, array2d)
1127 planes: list[MaskPlane] = list(plane_map.values())
1128 new_name_to_old_bitmask: dict[str, int] = {}
1129 for old_name, old_bit in old_planes.items():
1130 old_bitmask = _legacy_bitmask(old_bit)
1131 if new_plane := plane_map.get(old_name):
1132 # Already added to 'planes' at initialization.
1133 new_name_to_old_bitmask[new_plane.name] = old_bitmask
1134 else:
1135 if n_orphaned := np.count_nonzero(array2d & old_bitmask): 1135 ↛ 1136line 1135 didn't jump to line 1136 because the condition on line 1135 was never true
1136 raise RuntimeError(
1137 f"Legacy mask plane {old_name!r} is not remapped, "
1138 f"but {n_orphaned} pixels have this bit set."
1139 )
1140 schema = MaskSchema(planes)
1141 mask = Mask(0, schema=schema, yx0=yx0, shape=array2d.shape, sky_projection=sky_projection)
1142 for new_name, old_bitmask in new_name_to_old_bitmask.items():
1143 mask.set(new_name, array2d & old_bitmask)
1144 return mask
1146 @staticmethod
1147 def read_legacy(
1148 uri: ResourcePathExpression,
1149 *,
1150 plane_map: Mapping[str, MaskPlane] | None = None,
1151 ext: str | int = 1,
1152 fits_wcs_frame: Frame | None = None,
1153 ) -> Mask:
1154 """Read a FITS file written by `lsst.afw.image.Mask.writeFits`.
1156 Parameters
1157 ----------
1158 uri
1159 URI or file name.
1160 plane_map
1161 A mapping from legacy mask plane name to the new plane name and
1162 description. If not provided, the right legacy mask plane will be
1163 guessed, but this can depend on which mask planes the legacy
1164 mask actually has set. Planes from
1165 `get_legacy_optional_mask_planes` are added when the legacy mask
1166 defines them and has pixels set in them.
1167 ext
1168 Name or index of the FITS HDU to read.
1169 fits_wcs_frame
1170 If not `None` and the HDU containing the mask has a FITS WCS,
1171 attach a `SkyProjection` to the returned mask by converting that
1172 WCS.
1174 Returns
1175 -------
1176 `Mask`
1177 Mask read from the file.
1178 """
1179 opaque_metadata = fits.FitsOpaqueMetadata()
1180 fs, fspath = ResourcePath(uri).to_fsspec()
1181 with fs.open(fspath) as stream, astropy.io.fits.open(stream) as hdu_list:
1182 native_metadata = opaque_metadata.extract_legacy_primary_header(hdu_list[0].header)
1183 result = Mask._read_legacy_hdu(
1184 hdu_list[ext], opaque_metadata, plane_map=plane_map, fits_wcs_frame=fits_wcs_frame
1185 )
1186 result._opaque_metadata = opaque_metadata
1187 result._metadata = native_metadata
1188 return result
1190 @staticmethod
1191 def _read_legacy_hdu(
1192 hdu: astropy.io.fits.ImageHDU | astropy.io.fits.CompImageHDU | astropy.io.fits.BinTableHDU,
1193 opaque_metadata: fits.FitsOpaqueMetadata,
1194 plane_map: Mapping[str, MaskPlane] | None = None,
1195 fits_wcs_frame: Frame | None = None,
1196 strip_legacy_planes: bool = True,
1197 ) -> Mask:
1198 if isinstance(hdu, astropy.io.fits.BinTableHDU): 1198 ↛ 1199line 1198 didn't jump to line 1199 because the condition on line 1198 was never true
1199 hdu = astropy.io.fits.CompImageHDU(bintable=hdu)
1200 yx0 = fits.read_yx0(hdu.header)
1201 hdu.header.remove("LTV1", ignore_missing=True)
1202 hdu.header.remove("LTV2", ignore_missing=True)
1203 sky_projection: SkyProjection | None = None
1204 if fits_wcs_frame is not None: 1204 ↛ 1205line 1204 didn't jump to line 1205 because the condition on line 1204 was never true
1205 try:
1206 fits_wcs = astropy.wcs.WCS(hdu.header)
1207 except KeyError:
1208 pass
1209 else:
1210 sky_projection = SkyProjection.from_fits_wcs(
1211 fits_wcs, pixel_frame=fits_wcs_frame, x0=yx0.x, y0=yx0.y
1212 )
1213 if any(card.keyword.startswith("MSKN") for card in hdu.header.cards):
1214 # New ``lsst.images`` form: plane definitions are self-describing
1215 # via MSKN/MSKM/MSKD cards, so no plane_map is needed. The on-disk
1216 # array packs every plane into one element; ``set`` repacks each
1217 # plane into the (default uint8) in-memory layout by name.
1218 schema = MaskSchema.from_fits_header(hdu.header)
1219 mask = Mask(0, schema=schema, yx0=yx0, shape=hdu.data.shape, sky_projection=sky_projection)
1220 for n, plane in enumerate(schema):
1221 if plane is not None: 1221 ↛ 1220line 1221 didn't jump to line 1220 because the condition on line 1221 was always true
1222 mask.set(plane.name, hdu.data & hdu.header.get(f"MSKM{n:04d}", 1 << n))
1223 schema.strip_header(hdu.header)
1224 else:
1225 # Legacy ``lsst.afw.image`` form: bit indices in MP_* cards are
1226 # mapped to new planes via ``plane_map``.
1227 old_planes = MaskPlane.read_legacy(hdu.header, strip=strip_legacy_planes)
1228 resolved_map = _resolve_legacy_plane_map(plane_map, old_planes, hdu.data)
1229 mask = Mask._from_legacy_array(
1230 hdu.data, old_planes, yx0=yx0, plane_map=resolved_map, sky_projection=sky_projection
1231 )
1232 if not strip_legacy_planes:
1233 # Keep the MP_ cards for backwards compatibility, but re-index
1234 # them to the (reshuffled) positions of the new schema so a
1235 # legacy reader sees each plane at the bit it is actually
1236 # packed into on disk.
1237 _reindex_legacy_plane_cards(hdu.header, old_planes, resolved_map, mask.schema)
1238 fits.strip_wcs_cards(hdu.header)
1239 hdu.header.strip()
1240 hdu.header.remove("EXTTYPE", ignore_missing=True)
1241 hdu.header.remove("INHERIT", ignore_missing=True)
1242 # afw set BUNIT on masks because of limitations in how FITS
1243 # metadata is handled there.
1244 hdu.header.remove("BUNIT", ignore_missing=True)
1245 opaque_metadata.add_header(hdu.header)
1246 return mask
1249class MaskSerializationModel[P: pydantic.BaseModel](ArchiveTree):
1250 """Pydantic model used to represent the serialized form of a `.Mask`."""
1252 SCHEMA_NAME: ClassVar[str] = "mask"
1253 SCHEMA_VERSION: ClassVar[str] = "1.0.0"
1254 MIN_READ_VERSION: ClassVar[int] = 1
1255 PUBLIC_TYPE: ClassVar[type] = Mask
1257 data: list[ArrayReferenceModel | InlineArrayModel] = pydantic.Field(
1258 description="References to pixel data."
1259 )
1260 yx0: list[int] = pydantic.Field(
1261 description="Coordinate of the first pixels in the array, ordered (y, x)."
1262 )
1263 planes: list[MaskPlane | None] = pydantic.Field(description="Definitions of the bitplanes in the mask.")
1264 dtype: IntegerType = pydantic.Field(description="Data type of the in-memory mask.")
1265 sky_projection: SkyProjectionSerializationModel[P] | None = pydantic.Field(
1266 default=None,
1267 exclude_if=is_none,
1268 description="Projection that maps the logical pixel grid onto the sky.",
1269 )
1271 @property
1272 def bbox(self) -> Box:
1273 """The 2-d bounding box of the mask."""
1274 shape = self.data[0].shape
1275 if len(shape) == 3:
1276 shape = shape[:2]
1277 return Box.from_shape(shape, start=self.yx0)
1279 def deserialize(
1280 self,
1281 archive: InputArchive[Any],
1282 *,
1283 bbox: Box | None = None,
1284 strip_header: Callable[[astropy.io.fits.Header], None] = no_header_updates,
1285 **kwargs: Any,
1286 ) -> Mask:
1287 """Deserialize a mask from an input archive.
1289 Parameters
1290 ----------
1291 archive
1292 Archive to read from.
1293 bbox
1294 Bounding box of a subimage to read instead.
1295 strip_header
1296 A callable that strips out any FITS header cards added by the
1297 ``update_header`` argument in the corresponding call to
1298 `Mask.serialize`.
1299 **kwargs
1300 Unsupported keyword arguments are accepted only to provide better
1301 error messages (raising `serialization.InvalidParameterError`).
1302 """
1303 if kwargs: 1303 ↛ 1304line 1303 didn't jump to line 1304 because the condition on line 1303 was never true
1304 raise InvalidParameterError(f"Unrecognized parameters for Mask: {set(kwargs.keys())}.")
1306 def strip_header_and_legacy_planes(header: astropy.io.fits.Header) -> None:
1307 # The authoritative schema comes from the serialized tree, so drop
1308 # any legacy MP_* cards (written only for afw compatibility in the
1309 # legacy-cutout scenario) rather than carrying them as opaque
1310 # metadata, where they could drift out of sync or be re-propagated.
1311 strip_header(header)
1312 _strip_legacy_plane_cards(header)
1314 slices: tuple[slice, ...] | EllipsisType = ...
1315 if bbox is not None:
1316 slices = bbox.slice_within(self.bbox)
1317 else:
1318 bbox = self.bbox
1319 if not is_integer(self.dtype): 1319 ↛ 1320line 1319 didn't jump to line 1320 because the condition on line 1319 was never true
1320 raise ArchiveReadError(f"Mask array has a non-integer dtype: {self.dtype}.")
1321 schema = MaskSchema(self.planes, dtype=self.dtype.to_numpy())
1322 sky_projection = self.sky_projection.deserialize(archive) if self.sky_projection is not None else None
1323 if len(self.data) == 1 and tuple(self.data[0].shape) == tuple(self.bbox.shape) + (schema.mask_size,):
1324 storage_slices = slices if slices is ... else (slice(None),) + slices
1325 array = archive.get_array(
1326 self.data[0], strip_header=strip_header_and_legacy_planes, slices=storage_slices
1327 )
1328 array = np.moveaxis(array, 0, -1)
1329 return Mask(array, schema=schema, bbox=bbox, sky_projection=sky_projection)._finish_deserialize(
1330 self
1331 )
1332 result = Mask(0, schema=schema, bbox=bbox, sky_projection=sky_projection)
1333 schemas_2d = schema.split(np.int32)
1334 if len(schemas_2d) != len(self.data): 1334 ↛ 1335line 1334 didn't jump to line 1335 because the condition on line 1334 was never true
1335 raise ArchiveReadError(
1336 f"Number of mask arrays ({len(self.data)}) does not match expectation ({len(schemas_2d)})."
1337 )
1338 for array_model, schema_2d in zip(self.data, schemas_2d):
1339 mask_2d = self._deserialize_2d(
1340 array_model,
1341 schema_2d,
1342 bbox.start,
1343 archive,
1344 strip_header=strip_header_and_legacy_planes,
1345 slices=slices,
1346 )
1347 result.update(mask_2d)
1348 return result._finish_deserialize(self)
1350 @staticmethod
1351 def _deserialize_2d(
1352 ref: ArrayReferenceModel | InlineArrayModel,
1353 schema_2d: MaskSchema,
1354 yx0: Sequence[int],
1355 archive: InputArchive[Any],
1356 *,
1357 slices: tuple[slice, ...] | EllipsisType = ...,
1358 strip_header: Callable[[astropy.io.fits.Header], None] = no_header_updates,
1359 ) -> Mask:
1360 def _strip_header(header: astropy.io.fits.Header) -> None:
1361 strip_header(header)
1362 schema_2d.strip_header(header)
1363 fits.strip_wcs_cards(header)
1365 array_2d = archive.get_array(ref, strip_header=_strip_header, slices=slices)
1366 return Mask(array_2d[:, :, np.newaxis], schema=schema_2d, yx0=yx0)
1368 def deserialize_component(self, component: str, archive: InputArchive[Any], **kwargs: Any) -> Any:
1369 if kwargs:
1370 raise InvalidParameterError(f"Unsupported parameters for Mask components: {set(kwargs.keys())}.")
1371 return super().deserialize_component(component, archive)
1374def _archive_prefers_native_mask_arrays(archive: OutputArchive[Any]) -> bool:
1375 """Return whether an archive wants masks in their native 3-D layout."""
1376 current: Any = archive
1377 while current is not None:
1378 if getattr(current, "_prefer_native_mask_arrays", False):
1379 return True
1380 current = getattr(current, "_parent", None)
1381 return False
1384def get_legacy_visit_image_mask_planes() -> dict[str, MaskPlane]:
1385 """Return a mapping from legacy mask plane name to `MaskPlane` instance
1386 for LSST visit images, c. DP2.
1387 """
1388 return {
1389 "BAD": MaskPlane("BAD", "Bad pixel in the instrument, including bad amplifiers."),
1390 "SAT": MaskPlane(
1391 "SATURATED", "Pixel was saturated or affected by saturation in a neighboring pixel."
1392 ),
1393 "INTRP": MaskPlane("INTERPOLATED", "Original pixel value was interpolated."),
1394 "CR": MaskPlane("COSMIC_RAY", "A cosmic ray affected this pixel."),
1395 "EDGE": MaskPlane(
1396 "DETECTION_EDGE",
1397 "Pixel was too close to the edge to be considered for detection, "
1398 "due to the finite size of the detection kernel.",
1399 ),
1400 "DETECTED": MaskPlane("DETECTED", "Pixel was part of a detected source."),
1401 "SUSPECT": MaskPlane("SUSPECT", "Pixel was close to the saturation level. "),
1402 "NO_DATA": MaskPlane("NO_DATA", "No data was available for this pixel."),
1403 "VIGNETTED": MaskPlane("VIGNETTED", "Pixel was vignetted by the optics."),
1404 "PARTLY_VIGNETTED": MaskPlane("PARTLY_VIGNETTED", "Pixel was partly vignetted by the optics."),
1405 "CROSSTALK": MaskPlane("CROSSTALK", "Pixel was affected by crosstalk and corrected accordingly."),
1406 "ITL_DIP": MaskPlane(
1407 "ITL_DIP", "Pixel was affected by a dark vertical trail from a bright source, on an ITL CCD."
1408 ),
1409 "NOT_DEBLENDED": MaskPlane(
1410 "NOT_DEBLENDED",
1411 "Pixel belonged to a detection that was not deblended, usually due to size limits.",
1412 ),
1413 "SPIKE": MaskPlane(
1414 "SPIKE", "Pixel is in the neighborhood of a diffraction spike from a bright star."
1415 ),
1416 "UNMASKEDNAN": MaskPlane("UNMASKED_NAN", "Pixel was found to be NaN unexpectedly."),
1417 }
1420def get_legacy_difference_image_mask_planes() -> dict[str, MaskPlane]:
1421 """Return a mapping from legacy mask plane name to `MaskPlane` instance
1422 for LSST difference images, c. DP2.
1423 """
1424 result = get_legacy_visit_image_mask_planes()
1425 result["DETECTED_NEGATIVE"] = MaskPlane(
1426 "DETECTED_NEGATIVE", "Pixel was part of a detected source with negative flux."
1427 )
1428 result["SAT_TEMPLATE"] = MaskPlane("SAT_TEMPLATE", "Template pixel was saturated.")
1429 result["HIGH_VARIANCE"] = MaskPlane(
1430 "HIGH_VARIANCE", "Template pixel had fewer-than-usual input epochs and hence high noise."
1431 )
1432 result["STREAK"] = MaskPlane(
1433 "STREAK", "An extended streak (probably an artificial satellite) affected this pixel."
1434 )
1435 return result
1438def get_legacy_deep_coadd_mask_planes() -> dict[str, MaskPlane]:
1439 """Return a mapping from legacy mask plane name to `MaskPlane` instance
1440 for LSST deep coadds, c. DP2.
1441 """
1442 return {
1443 "NO_DATA": MaskPlane("NO_DATA", "No data was available for this pixel."),
1444 "INTRP": MaskPlane("INTERPOLATED", "Pixel value is the result of interpolating nearby good pixels."),
1445 "CR": MaskPlane(
1446 "COSMIC_RAY",
1447 "A cosmic ray affected this pixel on at least one input image (and was interpolated).",
1448 ),
1449 "SAT": MaskPlane(
1450 "SATURATED",
1451 "More than 10% of the potential input visits had a saturated pixel at this location "
1452 "('potential' because saturated pixel values are not actually propagated to the coadd). "
1453 "SATURATED always implies REJECTED, and is often a reason for NO_DATA.",
1454 ),
1455 "EDGE": MaskPlane(
1456 "DETECTION_EDGE",
1457 "Pixel was too close to the edge of the patch to be considered for detection, "
1458 "due to the finite size of the detection kernel.",
1459 ),
1460 "CLIPPED": MaskPlane(
1461 "CLIPPED",
1462 "Region was identified as a probable artifact when comparing multiple single-visit warps. "
1463 "CLIPPED always implies REJECTED.",
1464 ),
1465 "REJECTED": MaskPlane(
1466 "REJECTED",
1467 "At least one input visit was left out of the coadd for this pixel due to masking. "
1468 "REJECTED always implies INEXACT_PSF.",
1469 ),
1470 "DETECTED": MaskPlane("DETECTED", "Pixel was part of a detected source."),
1471 "INEXACT_PSF": MaskPlane(
1472 "INEXACT_PSF",
1473 "The set of visits contributing to this pixel differs from the set of visits "
1474 "contributing to the PSF model for its cell.",
1475 ),
1476 }
1479def get_legacy_non_cell_coadd_mask_planes() -> dict[str, MaskPlane]:
1480 """Return a mapping from legacy mask plane name to `MaskPlane` instance
1481 for LSST non-cell coadds such as ``template_coadd`` in DP2, and all
1482 DP1 coadds.
1484 These coadds carry the visit-level planes propagated from their input
1485 warps in addition to the coadd-specific planes, and flag chip edges with
1486 ``SENSOR_EDGE`` (cell coadds use ``CELL_EDGE`` instead).
1487 """
1488 result = get_legacy_deep_coadd_mask_planes()
1489 result["BAD"] = MaskPlane("BAD", "Bad pixel in the instrument, including bad amplifiers.")
1490 result["SUSPECT"] = MaskPlane("SUSPECT", "Pixel was close to the saturation level.")
1491 result["CROSSTALK"] = MaskPlane("CROSSTALK", "Pixel was affected by crosstalk and corrected accordingly.")
1492 result["DETECTED_NEGATIVE"] = MaskPlane(
1493 "DETECTED_NEGATIVE", "Pixel was part of a detected source with negative flux."
1494 )
1495 result["NOT_DEBLENDED"] = MaskPlane(
1496 "NOT_DEBLENDED",
1497 "Pixel belonged to a detection that was not deblended, usually due to size limits.",
1498 )
1499 result["UNMASKEDNAN"] = MaskPlane("UNMASKED_NAN", "Pixel was found to be NaN unexpectedly.")
1500 result["SENSOR_EDGE"] = MaskPlane(
1501 "SENSOR_EDGE",
1502 "Pixel is near the edge of a contributing sensor/chip, so the coadd PSF is discontinuous there.",
1503 )
1504 return result
1507def get_legacy_template_mask_planes() -> dict[str, MaskPlane]:
1508 """Return a mapping from legacy mask plane name to `MaskPlane` instance
1509 for image differencing template images.
1511 Returns
1512 -------
1513 plane_map : `dict` [`str`, `MaskPlane`]
1514 Mapping from legacy mask plane name to the plane it becomes.
1515 """
1516 result = get_legacy_non_cell_coadd_mask_planes()
1517 result["HIGH_VARIANCE"] = MaskPlane(
1518 "HIGH_VARIANCE", "Template pixel had fewer-than-usual input epochs and hence high noise."
1519 )
1520 return result
1523def get_legacy_optional_mask_planes() -> dict[str, MaskPlane]:
1524 """Return a mapping from legacy mask plane name to `MaskPlane` instance
1525 for the planes that only some images carry.
1527 Source injection defines these planes, but otherwise they are not used.
1529 Returns
1530 -------
1531 plane_map : `dict` [`str`, `MaskPlane`]
1532 Mapping from legacy mask plane name to the plane it becomes.
1533 """
1534 return {
1535 "INJECTED": MaskPlane(
1536 "INJECTED", "Pixel was affected by a synthetic source injected into the image."
1537 ),
1538 "INJECTED_CORE": MaskPlane(
1539 "INJECTED_CORE", "Pixel is in the core of a synthetic source injected into the image."
1540 ),
1541 "INJECTED_TEMPLATE": MaskPlane(
1542 "INJECTED_TEMPLATE", "Pixel was affected by a synthetic source injected into the template."
1543 ),
1544 "INJECTED_CORE_TEMPLATE": MaskPlane(
1545 "INJECTED_CORE_TEMPLATE",
1546 "Pixel is in the core of a synthetic source injected into the template.",
1547 ),
1548 }
1551def _legacy_bitmask(old_bit: int) -> int:
1552 """Return the bitmask that selects a legacy mask plane's bit."""
1553 old_bitmask = 1 << old_bit
1554 if old_bitmask == 2147483648: 1554 ↛ 1557line 1554 didn't jump to line 1557 because the condition on line 1554 was never true
1555 # afw uses int32 masks, but relies on overflow wrapping, which
1556 # numpy doesn't like.
1557 old_bitmask = -2147483648
1558 return old_bitmask
1561def _resolve_legacy_plane_map(
1562 plane_map: Mapping[str, MaskPlane] | None,
1563 old_planes: Mapping[str, int],
1564 array: np.ndarray,
1565) -> dict[str, MaskPlane]:
1566 """Return the plane map to convert a legacy mask array with.
1568 Parameters
1569 ----------
1570 plane_map
1571 Map given by the caller, or `None` to guess one from ``old_planes``.
1572 old_planes
1573 Mapping from legacy mask plane name to its bit index.
1574 array
1575 Legacy mask array, read to see which optional planes are in use.
1577 Returns
1578 -------
1579 plane_map
1580 The given (or guessed) map, plus the planes of
1581 `get_legacy_optional_mask_planes` that ``array`` has pixels set in.
1582 """
1583 resolved = dict(plane_map) if plane_map is not None else _guess_legacy_plane_map(old_planes)
1584 for old_name, plane in get_legacy_optional_mask_planes().items():
1585 if old_name in resolved or (old_bit := old_planes.get(old_name)) is None:
1586 continue
1587 if np.any(array & _legacy_bitmask(old_bit)):
1588 resolved[old_name] = plane
1589 return resolved
1592def _guess_legacy_plane_map(old_planes: Mapping[str, int]) -> dict[str, MaskPlane]:
1593 """Guess which of the ``get_legacy_*_plane_map`` created the given mask
1594 plane dictionary and call it.
1595 """
1596 if "SAT_TEMPLATE" in old_planes: 1596 ↛ 1597line 1596 didn't jump to line 1597 because the condition on line 1596 was never true
1597 return get_legacy_difference_image_mask_planes()
1598 if "INEXACT_PSF" in old_planes:
1599 # Both cell and non-cell coadds have INEXACT_PSF, but only non-cell
1600 # (assemble_coadd) coadds flag chip edges with SENSOR_EDGE; cell coadds
1601 # use CELL_EDGE.
1602 if "SENSOR_EDGE" in old_planes:
1603 if "HIGH_VARIANCE" in old_planes:
1604 return get_legacy_template_mask_planes()
1605 return get_legacy_non_cell_coadd_mask_planes()
1606 return get_legacy_deep_coadd_mask_planes()
1607 return get_legacy_visit_image_mask_planes()
1610def _reindex_legacy_plane_cards(
1611 header: astropy.io.fits.Header,
1612 old_planes: Mapping[str, int],
1613 plane_map: Mapping[str, MaskPlane],
1614 schema: MaskSchema,
1615) -> None:
1616 """Rewrite retained legacy ``MP_`` cards in place to match a reshuffled
1617 schema.
1619 Parameters
1620 ----------
1621 header
1622 Header whose ``MP_`` cards are updated in place.
1623 old_planes
1624 Mapping from legacy mask plane name to its original (on-disk) bit
1625 index, as returned by `MaskPlane.read_legacy`.
1626 plane_map
1627 Mapping from legacy mask plane name to the `MaskPlane` it was remapped
1628 to in ``schema``.
1629 schema
1630 The reconstructed schema that defines the new bit positions.
1632 Notes
1633 -----
1634 Each ``MP_<legacy name>`` card is set to the index that its remapped plane
1635 occupies in ``schema`` (equivalently, the ``MSKN`` index written on
1636 serialization). Cards for legacy planes that are not represented in the
1637 new schema are removed, since they no longer correspond to any stored bit.
1638 Legacy masks have at most 31 planes, so every plane maps to a single bit in
1639 one on-disk element and the index is unambiguous.
1640 """
1641 new_index = {plane.name: n for n, plane in enumerate(schema) if plane is not None}
1642 for old_name in old_planes:
1643 keyword = f"MP_{old_name}"
1644 new_plane = plane_map.get(old_name)
1645 if new_plane is not None and (index := new_index.get(new_plane.name)) is not None:
1646 header[keyword] = index
1647 else:
1648 del header[keyword]
1651def _strip_legacy_plane_cards(header: astropy.io.fits.Header) -> None:
1652 """Remove all legacy ``MP_*`` mask-plane cards from a FITS header.
1654 These are written only so that legacy tooling can read masks reconstructed
1655 from legacy cutouts; the ``lsst.images`` reader uses the serialized schema
1656 instead, so it strips them rather than carrying them as opaque metadata.
1657 """
1658 for keyword in [card.keyword for card in header.cards if card.keyword.startswith("MP_")]:
1659 del header[keyword]