Coverage for python/lsst/scarlet/lite/image.py: 19%
443 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-07 01:40 -0800
1# This file is part of scarlet_lite.
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/>.
22from __future__ import annotations
24import operator
25from typing import Any, Callable, Sequence, cast
27import numpy as np
28from numpy.typing import DTypeLike
30from .bbox import Box
31from .utils import ScalarLike, ScalarTypes
33__all__ = ["Image", "MismatchedBoxError", "MismatchedBandsError"]
36class MismatchedBandsError(Exception):
37 """Attempt to compare images with different bands"""
40class MismatchedBoxError(Exception):
41 """Attempt to compare images in different bounding boxes"""
44def get_dtypes(*data: np.ndarray | Image | ScalarLike) -> list[DTypeLike]:
45 """Get a list of dtypes from a list of arrays, images, or scalars
47 Parameters
48 ----------
49 data:
50 The arrays to use for calculating the dtype
52 Returns
53 -------
54 result:
55 A list of datatypes.
56 """
57 dtypes: list[DTypeLike] = [None] * len(data)
58 for d, element in enumerate(data):
59 if hasattr(element, "dtype"):
60 dtypes[d] = cast(np.ndarray, element).dtype
61 else:
62 dtypes[d] = np.dtype(type(element))
63 return dtypes
66def get_combined_dtype(*data: np.ndarray | Image | ScalarLike) -> DTypeLike:
67 """Get the combined dtype for a collection of arrays to prevent loss
68 of precision.
70 Parameters
71 ----------
72 data:
73 The arrays to use for calculating the dtype
75 Returns
76 -------
77 result: np.dtype
78 The resulting dtype.
79 """
80 dtypes = get_dtypes(*data)
81 return max(dtypes) # type: ignore
84class Image:
85 """A numpy array with an origin and (optional) bands
87 This class contains a 2D numpy array with the addition of an
88 origin (``yx0``) and an optional first index (``bands``) that
89 allows an immutable named index to be used.
91 Notes
92 -----
93 One of the main limitations of using numpy arrays to store image data
94 is the lack of an ``origin`` attribute that allows an array to retain
95 knowledge of it's location in a larger scene.
96 For example, if a numpy array ``x`` is sliced, eg. ``x[10:20, 30:40]``
97 the result will be a new ``10x10`` numpy array that has no meta
98 data to inform the user that it was sliced from a larger image.
99 In addition, astrophysical images are also multi-band data cubes,
100 with a 2D image in each band (in fact this is the simplifying
101 assumption that distinguishes scarlet lite from scarlet main).
102 However, the ordering of the bands during processing might differ from
103 the ordering of the bands to display multiband data.
104 So a mechanism was also desired to simplify the sorting and index of
105 an image by band name.
107 Thus, scarlet lite creates a numpy-array like class with the additional
108 ``bands`` and ``yx0`` attributes to keep track of the bands contained
109 in an array and the origin of that array (we specify ``yx0`` as opposed
110 to ``xy0`` to be consistent with the default numpy/C++ ``(y, x)``
111 ordering of arrays as opposed to the traditional cartesian ``(x, y)``
112 ordering used in astronomy and other modules in the science pipelines.
113 While this may be a small source of confusion for the user,
114 it is consistent with the ordering in the original scarlet package and
115 ensures the consistency of scarlet lite images and python index slicing.
117 Examples
118 --------
120 The easiest way to create a new image is to use ``Image(numpy_array)``,
121 for example
123 >>> import numpy as np
124 >>> from lsst.scarlet.lite import Image
125 >>>
126 >>> x = np.arange(12).reshape(3, 4)
127 >>> image = Image(x)
128 >>> print(image)
129 Image:
130 [[ 0 1 2 3]
131 [ 4 5 6 7]
132 [ 8 9 10 11]]
133 bands=()
134 bbox=Box(shape=(3, 4), origin=(0, 0))
136 This will create a single band :py:class:`~lsst.scarlet.lite.Image` with
137 origin ``(0, 0)``.
138 To create a multi-band image the input array must have 3 dimensions and
139 the ``bands`` property must be specified:
141 >>> x = np.arange(24).reshape(2, 3, 4)
142 >>> image = Image(x, bands=("i", "z"))
143 >>> print(image)
144 Image:
145 [[[ 0 1 2 3]
146 [ 4 5 6 7]
147 [ 8 9 10 11]]
148 <BLANKLINE>
149 [[12 13 14 15]
150 [16 17 18 19]
151 [20 21 22 23]]]
152 bands=('i', 'z')
153 bbox=Box(shape=(3, 4), origin=(0, 0))
155 It is also possible to create an empty single-band image using the
156 ``from_box`` static method:
158 >>> from lsst.scarlet.lite import Box
159 >>> image = Image.from_box(Box((3, 4), (100, 120)))
160 >>> print(image)
161 Image:
162 [[0. 0. 0. 0.]
163 [0. 0. 0. 0.]
164 [0. 0. 0. 0.]]
165 bands=()
166 bbox=Box(shape=(3, 4), origin=(100, 120))
168 Similarly, an empty multi-band image can be created by passing a tuple
169 of ``bands``:
171 >>> image = Image.from_box(Box((3, 4)), bands=("r", "i"))
172 >>> print(image)
173 Image:
174 [[[0. 0. 0. 0.]
175 [0. 0. 0. 0.]
176 [0. 0. 0. 0.]]
177 <BLANKLINE>
178 [[0. 0. 0. 0.]
179 [0. 0. 0. 0.]
180 [0. 0. 0. 0.]]]
181 bands=('r', 'i')
182 bbox=Box(shape=(3, 4), origin=(0, 0))
184 To select a sub-image use a ``Box`` to select a spatial region in either a
185 single-band or multi-band image:
187 >>> x = np.arange(60).reshape(3, 4, 5)
188 >>> image = Image(x, bands=("g", "r", "i"), yx0=(20, 30))
189 >>> bbox = Box((2, 2), (21, 32))
190 >>> print(image[bbox])
191 Image:
192 [[[ 7 8]
193 [12 13]]
194 <BLANKLINE>
195 [[27 28]
196 [32 33]]
197 <BLANKLINE>
198 [[47 48]
199 [52 53]]]
200 bands=('g', 'r', 'i')
201 bbox=Box(shape=(2, 2), origin=(21, 32))
204 To select a single-band image from a multi-band image,
205 pass the name of the band as an index:
207 >>> print(image["r"])
208 Image:
209 [[20 21 22 23 24]
210 [25 26 27 28 29]
211 [30 31 32 33 34]
212 [35 36 37 38 39]]
213 bands=()
214 bbox=Box(shape=(4, 5), origin=(20, 30))
216 Multi-band images can also be sliced in the spatial dimension, for example
218 >>> print(image["g":"r"])
219 Image:
220 [[[ 0 1 2 3 4]
221 [ 5 6 7 8 9]
222 [10 11 12 13 14]
223 [15 16 17 18 19]]
224 <BLANKLINE>
225 [[20 21 22 23 24]
226 [25 26 27 28 29]
227 [30 31 32 33 34]
228 [35 36 37 38 39]]]
229 bands=('g', 'r')
230 bbox=Box(shape=(4, 5), origin=(20, 30))
232 and
234 >>> print(image["r":"r"])
235 Image:
236 [[[20 21 22 23 24]
237 [25 26 27 28 29]
238 [30 31 32 33 34]
239 [35 36 37 38 39]]]
240 bands=('r',)
241 bbox=Box(shape=(4, 5), origin=(20, 30))
243 both extract a slice of a multi-band image.
245 .. warning::
246 Unlike numerical indices, where ``slice(x, y)`` will select the
247 subset of an array from ``x`` to ``y-1`` (excluding ``y``),
248 a spectral slice of an ``Image`` will return the image slice
249 including band ``y``.
251 It is also possible to change the order or index a subset of bands
252 in an image. For example:
254 >>> print(image[("r", "g", "i")])
255 Image:
256 [[[20 21 22 23 24]
257 [25 26 27 28 29]
258 [30 31 32 33 34]
259 [35 36 37 38 39]]
260 <BLANKLINE>
261 [[ 0 1 2 3 4]
262 [ 5 6 7 8 9]
263 [10 11 12 13 14]
264 [15 16 17 18 19]]
265 <BLANKLINE>
266 [[40 41 42 43 44]
267 [45 46 47 48 49]
268 [50 51 52 53 54]
269 [55 56 57 58 59]]]
270 bands=('r', 'g', 'i')
271 bbox=Box(shape=(4, 5), origin=(20, 30))
274 will return a new image with the bands re-ordered.
276 Images can be combined using the standard arithmetic operations similar to
277 numpy arrays, including ``+, -, *, /, **`` etc, however, if two images are
278 combined with different bounding boxes, the _union_ of the two
279 boxes is used for the result. For example:
281 >>> image1 = Image(np.ones((2, 3, 4)), bands=tuple("gr"))
282 >>> image2 = Image(np.ones((2, 3, 4)), bands=tuple("gr"), yx0=(2, 3))
283 >>> result = image1 + image2
284 >>> print(result)
285 Image:
286 [[[1. 1. 1. 1. 0. 0. 0.]
287 [1. 1. 1. 1. 0. 0. 0.]
288 [1. 1. 1. 2. 1. 1. 1.]
289 [0. 0. 0. 1. 1. 1. 1.]
290 [0. 0. 0. 1. 1. 1. 1.]]
291 <BLANKLINE>
292 [[1. 1. 1. 1. 0. 0. 0.]
293 [1. 1. 1. 1. 0. 0. 0.]
294 [1. 1. 1. 2. 1. 1. 1.]
295 [0. 0. 0. 1. 1. 1. 1.]
296 [0. 0. 0. 1. 1. 1. 1.]]]
297 bands=('g', 'r')
298 bbox=Box(shape=(5, 7), origin=(0, 0))
300 If instead you want to additively ``insert`` image 1 into image 2,
301 so that they have the same bounding box as image 2, use
303 >>> _ = image2.insert(image1)
304 >>> print(image2)
305 Image:
306 [[[2. 1. 1. 1.]
307 [1. 1. 1. 1.]
308 [1. 1. 1. 1.]]
309 <BLANKLINE>
310 [[2. 1. 1. 1.]
311 [1. 1. 1. 1.]
312 [1. 1. 1. 1.]]]
313 bands=('g', 'r')
314 bbox=Box(shape=(3, 4), origin=(2, 3))
316 To insert an image using a different operation use
318 >>> from operator import truediv
319 >>> _ = image2.insert(image1, truediv)
320 >>> print(image2)
321 Image:
322 [[[2. 1. 1. 1.]
323 [1. 1. 1. 1.]
324 [1. 1. 1. 1.]]
325 <BLANKLINE>
326 [[2. 1. 1. 1.]
327 [1. 1. 1. 1.]
328 [1. 1. 1. 1.]]]
329 bands=('g', 'r')
330 bbox=Box(shape=(3, 4), origin=(2, 3))
333 However, depending on the operation you may get unexpected results
334 since now there could be ``NaN`` and ``inf`` values due to the zeros
335 in the non-overlapping regions.
336 Instead, to select only the overlap region one can use
338 >>> result = image1 / image2
339 >>> print(result[image1.bbox & image2.bbox])
340 Image:
341 [[[0.5]]
342 <BLANKLINE>
343 [[0.5]]]
344 bands=('g', 'r')
345 bbox=Box(shape=(1, 1), origin=(2, 3))
348 Parameters
349 ----------
350 data:
351 The array data for the image.
352 bands:
353 The bands coving the image.
354 yx0:
355 The (y, x) offset for the lower left of the image.
356 """
358 def __init__(
359 self,
360 data: np.ndarray,
361 bands: Sequence | None = None,
362 yx0: tuple[int, int] | None = None,
363 ):
364 if bands is None or len(bands) == 0:
365 # Using an empty tuple for the bands will result in a 2D image
366 bands = ()
367 assert data.ndim == 2
368 else:
369 bands = tuple(bands)
370 assert data.ndim == 3
371 if data.shape[0] != len(bands):
372 raise ValueError(f"Array has spectral size {data.shape[0]}, but {bands} bands")
373 if yx0 is None:
374 yx0 = (0, 0)
375 self._data = data
376 self._yx0 = yx0
377 self._bands = bands
379 @staticmethod
380 def from_box(bbox: Box, bands: tuple | None = None, dtype: DTypeLike = float) -> Image:
381 """Initialize an empty image from a bounding Box and optional bands
383 Parameters
384 ----------
385 bbox:
386 The bounding box that contains the image.
387 bands:
388 The bands for the image.
389 If bands is `None` then a 2D image is created.
390 dtype:
391 The numpy dtype of the image.
393 Returns
394 -------
395 image:
396 An empty image contained in ``bbox`` with ``bands`` bands.
397 """
398 if bands is not None and len(bands) > 0:
399 shape = (len(bands),) + bbox.shape
400 else:
401 shape = bbox.shape
402 data = np.zeros(shape, dtype=dtype)
403 return Image(data, bands=bands, yx0=cast(tuple[int, int], bbox.origin))
405 @property
406 def shape(self) -> tuple[int, ...]:
407 """The shape of the image.
409 This includes the spectral dimension, if there is one.
410 """
411 return self._data.shape
413 @property
414 def dtype(self) -> DTypeLike:
415 """The numpy dtype of the image."""
416 return self._data.dtype
418 @property
419 def bands(self) -> tuple:
420 """The bands used in the image."""
421 return self._bands
423 @property
424 def n_bands(self) -> int:
425 """Number of bands in the image.
427 If `n_bands == 0` then the image is 2D and does not have a spectral
428 dimension.
429 """
430 return len(self._bands)
432 @property
433 def is_multiband(self) -> bool:
434 """Whether or not the image has a spectral dimension."""
435 return self.n_bands > 0
437 @property
438 def height(self) -> int:
439 """Height of the image."""
440 return self.shape[-2]
442 @property
443 def width(self) -> int:
444 """Width of the image."""
445 return self.shape[-1]
447 @property
448 def yx0(self) -> tuple[int, int]:
449 """Origin of the image, in numpy/C++ y,x ordering."""
450 return self._yx0
452 @property
453 def y0(self) -> int:
454 """location of the y-offset."""
455 return self._yx0[0]
457 @property
458 def x0(self) -> int:
459 """Location of the x-offset."""
460 return self._yx0[1]
462 @property
463 def bbox(self) -> Box:
464 """Bounding box for the special dimensions in the image."""
465 return Box(self.shape[-2:], self._yx0)
467 @property
468 def data(self) -> np.ndarray:
469 """The image viewed as a numpy array."""
470 return self._data
472 @property
473 def ndim(self) -> int:
474 """Number of dimensions in the image."""
475 return self._data.ndim
477 def spectral_indices(self, bands: Sequence | slice) -> tuple[int, ...] | slice:
478 """The indices to extract each band in `bands` in order from the image
480 This converts a band name, or list of band names,
481 into numerical indices that can be used to slice the internal numpy
482 `data` array.
484 Parameters
485 ---------
486 bands:
487 If `bands` is a list of band names, then the result will be an
488 index corresponding to each band, in order.
489 If `bands` is a slice, then the ``start`` and ``stop`` properties
490 should be band names, and the result will be a slice with the
491 appropriate indices to start at `bands.start` and end at
492 `bands.stop`.
494 Returns
495 -------
496 band_indices:
497 Tuple of indices for each band in this image.
498 """
499 if isinstance(bands, slice):
500 # Convert a slice of band names into a slice of array indices
501 # to select the appropriate slice.
502 if bands.start is None:
503 start = None
504 else:
505 start = self.bands.index(bands.start)
506 if bands.stop is None:
507 stop = None
508 else:
509 stop = self.bands.index(bands.stop) + 1
510 return slice(start, stop, bands.step)
512 if isinstance(bands, str):
513 return (self.bands.index(bands),)
515 band_indices = tuple(self.bands.index(band) for band in bands if band in self.bands)
516 return band_indices
518 def matched_spectral_indices(
519 self,
520 other: Image,
521 ) -> tuple[tuple[int, ...] | slice, tuple[int, ...] | slice]:
522 """Match bands between two images
524 Parameters
525 ----------
526 other:
527 The other image to match spectral indices to.
529 Returns
530 -------
531 result:
532 A tuple with a tuple of indices/slices for each dimension,
533 including the spectral dimension.
534 """
535 if self.bands == other.bands and self.n_bands != 0:
536 # The bands match
537 return slice(None), slice(None)
538 if self.n_bands == 0 and other.n_bands == 0:
539 # The images are 2D, so no spectral slicing is necessary
540 return (), ()
541 if self.n_bands == 0 and other.n_bands > 1:
542 err = "Attempted to insert a monochromatic image into a mutli-band image"
543 raise ValueError(err)
544 if other.n_bands == 0:
545 err = "Attempted to insert a multi-band image into a monochromatic image"
546 raise ValueError(err)
548 self_indices = cast(tuple[int, ...], self.spectral_indices(other.bands))
549 matched_bands = tuple(self.bands[bidx] for bidx in self_indices)
550 other_indices = cast(tuple[int, ...], other.spectral_indices(matched_bands))
551 return other_indices, self_indices
553 def matched_slices(self, bbox: Box) -> tuple[tuple[slice, ...], tuple[slice, ...]]:
554 """Get the slices to match this image to a given bounding box
556 Parameters
557 ----------
558 bbox:
559 The bounding box to match this image to.
561 Returns
562 -------
563 result:
564 Tuple of indices/slices to match this image to the given bbox.
565 """
566 if self.bbox == bbox:
567 # No need to slice, since the boxes match
568 _slice = (slice(None),) * bbox.ndim
569 return _slice, _slice
571 slices = self.bbox.overlapped_slices(bbox)
572 return slices
574 def project(
575 self,
576 bands: object | tuple[object] | None = None,
577 bbox: Box | None = None,
578 ) -> Image:
579 """Project this image into a different set of bands
581 Parameters
582 ----------
583 bands:
584 Spectral bands to project this image into.
585 Not all bands have to be contained in the image, and not all
586 bands contained in the image have to be used in the projection.
587 bbox:
588 A bounding box to project the image into.
590 Results
591 -------
592 image:
593 A new image creating by projecting this image into
594 `bbox` and `bands`.
595 """
596 if bands is None:
597 bands = self.bands
598 if not isinstance(bands, tuple):
599 bands = (bands,)
600 if self.is_multiband:
601 indices = self.spectral_indices(bands)
602 data = self.data[indices, :]
603 else:
604 data = self.data
606 if bbox is None:
607 return Image(data, bands=bands, yx0=self.yx0)
609 if self.is_multiband:
610 image = np.zeros((len(bands),) + bbox.shape, dtype=data.dtype)
611 slices = bbox.overlapped_slices(self.bbox)
612 # Insert a slice for the spectral dimension
613 image[(slice(None),) + slices[0]] = data[(slice(None),) + slices[1]]
614 return Image(image, bands=bands, yx0=cast(tuple[int, int], bbox.origin))
616 image = np.zeros(bbox.shape, dtype=data.dtype)
617 slices = bbox.overlapped_slices(self.bbox)
618 image[slices[0]] = data[slices[1]]
619 return Image(image, bands=bands, yx0=cast(tuple[int, int], bbox.origin))
621 @property
622 def multiband_slices(self) -> tuple[tuple[int, ...] | slice, slice, slice]:
623 """Return the slices required to slice a multiband image"""
624 return (self.spectral_indices(self.bands),) + self.bbox.slices # type: ignore
626 def insert_into(
627 self,
628 image: Image,
629 op: Callable = operator.add,
630 ) -> Image:
631 """Insert this image into another image in place.
633 Parameters
634 ----------
635 image:
636 The image to insert this image into.
637 op:
638 The operator to use when combining the images.
640 Returns
641 -------
642 result:
643 `image` updated by inserting this instance.
644 """
645 return insert_image(image, self, op)
647 def insert(self, image: Image, op: Callable = operator.add) -> Image:
648 """Insert another image into this image in place.
650 Parameters
651 ----------
652 image:
653 The image to insert this image into.
654 op:
655 The operator to use when combining the images.
657 Returns
658 -------
659 result:
660 This instance with `image` inserted.
661 """
662 return insert_image(self, image, op)
664 def repeat(self, bands: tuple) -> Image:
665 """Project a 2D image into the spectral dimension
667 Parameters
668 ----------
669 bands:
670 The bands in the projected image.
672 Returns
673 -------
674 result: Image
675 The 2D image repeated in each band in the spectral dimension.
676 """
677 if self.is_multiband:
678 raise ValueError("Image.repeat only works with 2D images")
679 return self.copy_with(
680 np.repeat(self.data[None, :, :], len(bands), axis=0),
681 bands=bands,
682 yx0=self.yx0,
683 )
685 def copy(self, order=None) -> Image:
686 """Make a copy of this image.
688 Parameters
689 ----------
690 order:
691 The ordering to use for storing the bytes.
692 This is unlikely to be needed, and just defaults to
693 the numpy behavior (C) ordering.
695 Returns
696 -------
697 image: Image
698 The copy of this image.
699 """
700 return self.copy_with(order=order)
702 def copy_with(
703 self,
704 data: np.ndarray | None = None,
705 order: str | None = None,
706 bands: tuple[str, ...] | None = None,
707 yx0: tuple[int, int] | None = None,
708 ):
709 """Copy of this image with some parameters updated.
711 Any parameters not specified by the user will be copied from the
712 current image.
714 Parameters
715 ----------
716 data:
717 An update for the data in the image.
718 order:
719 The ordering for stored bytes, from numpy.copy.
720 bands:
721 The bands that the resulting image will have.
722 The number of bands must be the same as the first dimension
723 in the data array.
724 yx0:
725 The lower-left of the image bounding box.
727 Returns
728 -------
729 image: Image
730 The copied image.
731 """
732 if order is None:
733 order = "C"
734 if data is None:
735 data = self.data.copy(order) # type: ignore
736 if bands is None:
737 bands = self.bands
738 if yx0 is None:
739 yx0 = self.yx0
740 return Image(data, bands, yx0)
742 def trimmed(self, threshold: float = 0) -> Image:
743 """Return a copy of the image trimmed to a threshold.
745 This is essentially the smallest image that contains all of the
746 pixels above the threshold.
748 Parameters
749 ----------
750 threshold:
751 The threshold to use for trimming the image.
753 Returns
754 -------
755 image:
756 A copy of the image trimmed to the threshold.
757 """
758 data = self.data.copy()
759 bbox = Box.from_data(data, threshold=threshold)
760 data = data[bbox.slices]
761 y0, x0 = bbox.origin
763 return Image(data, yx0=(y0 + self.y0, x0 + self.x0))
765 def at(self, y: int, x: int) -> ScalarLike | np.ndarray:
766 """The value of the image at a given location.
768 Image does not implment single index access because the
769 result is a scalar, while indexing an image returns another image.
771 Parameters
772 ----------
773 y:
774 The y-coordinate of the location.
775 x:
776 The x-coordinate of the location.
778 Returns
779 -------
780 value:
781 The value of the image at the given location.
782 """
783 _y = y - self.y0
784 _x = x - self.x0
785 if self.ndim == 2:
786 return self.data[_y, _x]
787 return self.data[:, _y, _x]
789 def _i_update(self, op: Callable, other: Image | ScalarLike) -> Image:
790 """Update the data array in place.
792 This is typically implemented by `__i<op>__` methods,
793 like `__iadd__`, to apply an operator and update this image
794 with the data in place.
796 Parameters
797 ----------
798 op:
799 Operator used to combine this image with the `other` image.
800 other:
801 The other image that is combined with this one using the operator
802 `op`.
804 Returns
805 -------
806 image: Image
807 This image, after being updated by the operator
808 """
809 dtype = get_combined_dtype(self.data, other)
810 if self.dtype != dtype:
811 if hasattr(other, "dtype"):
812 _dtype = cast(np.ndarray, other).dtype
813 else:
814 _dtype = type(other) # type: ignore
815 msg = f"Cannot update an array with type {self.dtype} with {_dtype}"
816 raise ValueError(msg)
817 result = op(other)
818 self._data[:] = result.data
819 self._bands = result.bands
820 self._yx0 = result.yx0
821 return self
823 def _check_equality(self, other: Image | ScalarLike, op: Callable) -> Image:
824 """Compare this array to another.
826 This performs an element by element equality check.
828 Parameters
829 ----------
830 other:
831 The image to compare this image to.
832 op:
833 The operator used for the comparision (==, !=, >=, <=).
835 Returns
836 -------
837 image: Image
838 An image made by checking all of the elements in this array with
839 another.
841 Raises
842 ------
843 TypeError:
844 If `other` is not an `Image`.
845 MismatchedBandsError:
846 If `other` has different bands.
847 MismatchedBoxError:
848 if `other` exists in a different bounding box.
849 """
850 if isinstance(other, Image) and other.bands == self.bands and other.bbox == self.bbox:
851 return self.copy_with(data=op(self.data, other.data))
853 if not isinstance(other, Image):
854 if type(other) in ScalarTypes:
855 return self.copy_with(data=op(self.data, other))
856 raise TypeError(f"Cannot compare images to {type(other)}")
858 if other.bands != self.bands:
859 msg = f"Cannot compare images with mismatched bands: {self.bands} vs {other.bands}"
860 raise MismatchedBandsError(msg)
862 raise MismatchedBoxError(
863 f"Cannot compare images with different bounds boxes: {self.bbox} vs. {other.bbox}"
864 )
866 def __eq__(self, other: object) -> Image: # type: ignore
867 """Check if this image is equal to another."""
868 if not isinstance(other, Image) and not isinstance(other, ScalarTypes):
869 raise TypeError(f"Cannot compare an Image to {type(other)}.")
870 return self._check_equality(other, operator.eq) # type: ignore
872 def __ne__(self, other: object) -> Image: # type: ignore
873 """Check if this image is not equal to another."""
874 return ~self.__eq__(other)
876 def __ge__(self, other: Image | ScalarLike) -> Image:
877 """Check if this image is greater than or equal to another."""
878 if type(other) in ScalarTypes:
879 return self.copy_with(data=self.data >= other)
880 return self._check_equality(other, operator.ge)
882 def __le__(self, other: Image | ScalarLike) -> Image:
883 """Check if this image is less than or equal to another."""
884 if type(other) in ScalarTypes:
885 return self.copy_with(data=self.data <= other)
886 return self._check_equality(other, operator.le)
888 def __gt__(self, other: Image | ScalarLike) -> Image:
889 """Check if this image is greater than or equal to another."""
890 if type(other) in ScalarTypes:
891 return self.copy_with(data=self.data > other)
892 return self._check_equality(other, operator.ge)
894 def __lt__(self, other: Image | ScalarLike) -> Image:
895 """Check if this image is less than or equal to another."""
896 if type(other) in ScalarTypes:
897 return self.copy_with(data=self.data < other)
898 return self._check_equality(other, operator.le)
900 def __neg__(self):
901 """Take the negative of the image."""
902 return self.copy_with(data=-self._data)
904 def __pos__(self):
905 """Make a copy using of the image."""
906 return self.copy()
908 def __invert__(self):
909 """Take the inverse (~) of the image."""
910 return self.copy_with(data=~self._data)
912 def __add__(self, other: Image | ScalarLike) -> Image:
913 """Combine this image and another image using addition."""
914 return _operate_on_images(self, other, operator.add)
916 def __iadd__(self, other: Image | ScalarLike) -> Image:
917 """Combine this image and another image using addition and update
918 in place.
919 """
920 return self._i_update(self.__add__, other)
922 def __radd__(self, other: Image | ScalarLike) -> Image:
923 """Combine this image and another image using addition,
924 with this image on the right.
925 """
926 if type(other) in ScalarTypes:
927 return self.copy_with(data=other + self.data)
928 return cast(Image, other).__add__(self)
930 def __sub__(self, other: Image | ScalarLike) -> Image:
931 """Combine this image and another image using subtraction."""
932 return _operate_on_images(self, other, operator.sub)
934 def __isub__(self, other: Image | ScalarLike) -> Image:
935 """Combine this image and another image using subtraction,
936 with this image on the right.
937 """
938 return self._i_update(self.__sub__, other)
940 def __rsub__(self, other: Image | ScalarLike) -> Image:
941 """Combine this image and another image using subtraction,
942 with this image on the right.
943 """
944 if type(other) in ScalarTypes:
945 return self.copy_with(data=other - self.data)
946 return cast(Image, other).__sub__(self)
948 def __mul__(self, other: Image | ScalarLike) -> Image:
949 """Combine this image and another image using multiplication."""
950 return _operate_on_images(self, other, operator.mul)
952 def __imul__(self, other: Image | ScalarLike) -> Image:
953 """Combine this image and another image using multiplication,
954 with this image on the right.
955 """
956 return self._i_update(self.__mul__, other)
958 def __rmul__(self, other: Image | ScalarLike) -> Image:
959 """Combine this image and another image using multiplication,
960 with this image on the right.
961 """
962 if type(other) in ScalarTypes:
963 return self.copy_with(data=other * self.data)
964 return cast(Image, other).__mul__(self)
966 def __truediv__(self, other: Image | ScalarLike) -> Image:
967 """Divide this image by `other`."""
968 return _operate_on_images(self, other, operator.truediv)
970 def __itruediv__(self, other: Image | ScalarLike) -> Image:
971 """Divide this image by `other` in place."""
972 return self._i_update(self.__truediv__, other)
974 def __rtruediv__(self, other: Image | ScalarLike) -> Image:
975 """Divide this image by `other` with this on the right."""
976 if type(other) in ScalarTypes:
977 return self.copy_with(data=other / self.data)
978 return cast(Image, other).__truediv__(self)
980 def __floordiv__(self, other: Image | ScalarLike) -> Image:
981 """Floor divide this image by `other` in place."""
982 return _operate_on_images(self, other, operator.floordiv)
984 def __ifloordiv__(self, other: Image | ScalarLike) -> Image:
985 """Floor divide this image by `other` in place."""
986 return self._i_update(self.__floordiv__, other)
988 def __rfloordiv__(self, other: Image | ScalarLike) -> Image:
989 """Floor divide this image by `other` with this on the right."""
990 if type(other) in ScalarTypes:
991 return self.copy_with(data=other // self.data)
992 return cast(Image, other).__floordiv__(self)
994 def __pow__(self, other: Image | ScalarLike) -> Image:
995 """Raise this image to the `other` power."""
996 return _operate_on_images(self, other, operator.pow)
998 def __ipow__(self, other: Image | ScalarLike) -> Image:
999 """Raise this image to the `other` power in place."""
1000 return self._i_update(self.__pow__, other)
1002 def __rpow__(self, other: Image | ScalarLike) -> Image:
1003 """Raise this other to the power of this image."""
1004 if type(other) in ScalarTypes:
1005 return self.copy_with(data=other**self.data)
1006 return cast(Image, other).__pow__(self)
1008 def __mod__(self, other: Image | ScalarLike) -> Image:
1009 """Take the modulus of this % other."""
1010 return _operate_on_images(self, other, operator.mod)
1012 def __imod__(self, other: Image | ScalarLike) -> Image:
1013 """Take the modulus of this % other in place."""
1014 return self._i_update(self.__mod__, other)
1016 def __rmod__(self, other: Image | ScalarLike) -> Image:
1017 """Take the modulus of other % this."""
1018 if type(other) in ScalarTypes:
1019 return self.copy_with(data=other % self.data)
1020 return cast(Image, other).__mod__(self)
1022 def __and__(self, other: Image | ScalarLike) -> Image:
1023 """Take the bitwise and of this and other."""
1024 return _operate_on_images(self, other, operator.and_)
1026 def __iand__(self, other: Image | ScalarLike) -> Image:
1027 """Take the bitwise and of this and other in place."""
1028 return self._i_update(self.__and__, other)
1030 def __rand__(self, other: Image | ScalarLike) -> Image:
1031 """Take the bitwise and of other and this."""
1032 if type(other) in ScalarTypes:
1033 return self.copy_with(data=other & self.data)
1034 return cast(Image, other).__and__(self)
1036 def __or__(self, other: Image | ScalarLike) -> Image:
1037 """Take the binary or of this or other."""
1038 return _operate_on_images(self, other, operator.or_)
1040 def __ior__(self, other: Image | ScalarLike) -> Image:
1041 """Take the binary or of this or other in place."""
1042 return self._i_update(self.__or__, other)
1044 def __ror__(self, other: Image | ScalarLike) -> Image:
1045 """Take the binary or of other or this."""
1046 if type(other) in ScalarTypes:
1047 return self.copy_with(data=other | self.data)
1048 return cast(Image, other).__or__(self)
1050 def __xor__(self, other: Image | ScalarLike) -> Image:
1051 """Take the binary xor of this xor other."""
1052 return _operate_on_images(self, other, operator.xor)
1054 def __ixor__(self, other: Image | ScalarLike) -> Image:
1055 """Take the binary xor of this xor other in place."""
1056 return self._i_update(self.__xor__, other)
1058 def __rxor__(self, other: Image | ScalarLike) -> Image:
1059 """Take the binary xor of other xor this."""
1060 if type(other) in ScalarTypes:
1061 return self.copy_with(data=other ^ self.data)
1062 return cast(Image, other).__xor__(self)
1064 def __lshift__(self, other: ScalarLike) -> Image:
1065 """Shift this image to the left by other bits."""
1066 if not issubclass(np.dtype(type(other)).type, np.integer):
1067 raise TypeError("Bit shifting an image can only be done with integers")
1068 return self.copy_with(data=self.data << other)
1070 def __ilshift__(self, other: ScalarLike) -> Image:
1071 """Shift this image to the left by other bits in place."""
1072 self[:] = self.__lshift__(other)
1073 return self
1075 def __rlshift__(self, other: ScalarLike) -> Image:
1076 """Shift other to the left by this image bits."""
1077 return self.copy_with(data=other << self.data)
1079 def __rshift__(self, other: ScalarLike) -> Image:
1080 """Shift this image to the right by other bits."""
1081 if not issubclass(np.dtype(type(other)).type, np.integer):
1082 raise TypeError("Bit shifting an image can only be done with integers")
1083 return self.copy_with(data=self.data >> other)
1085 def __irshift__(self, other: ScalarLike) -> Image:
1086 """Shift this image to the right by other bits in place."""
1087 self[:] = self.__rshift__(other)
1088 return self
1090 def __rrshift__(self, other: ScalarLike) -> Image:
1091 """Shift other to the right by this image bits."""
1092 return self.copy_with(data=other >> self.data)
1094 def __str__(self):
1095 """Display the image array, bands, and bounding box."""
1096 return f"Image:\n {str(self.data)}\n bands={self.bands}\n bbox={self.bbox}"
1098 def _is_spectral_index(self, index: Any) -> bool:
1099 """Check to see if an index is a spectral index.
1101 Parameters
1102 ----------
1103 index:
1104 Either a slice, a tuple, or an element in `Image.bands`.
1106 Returns
1107 -------
1108 result:
1109 ``True`` if `index` is band or tuple of bands.
1110 """
1111 bands = self.bands
1112 if isinstance(index, slice):
1113 if index.start in bands or index.stop in bands or (index.start is None and index.stop is None):
1114 return True
1115 return False
1116 if index in self.bands:
1117 return True
1118 if isinstance(index, tuple) and index[0] in self.bands:
1119 return True
1120 return False
1122 def _get_box_slices(self, bbox: Box) -> tuple[slice, slice]:
1123 """Get the slices of the image to insert it into the overlapping
1124 region with `bbox`."""
1125 overlap = self.bbox & bbox
1126 if overlap != bbox:
1127 raise IndexError("Bounding box is outside of the image")
1128 origin = bbox.origin
1129 shape = bbox.shape
1130 y_start = origin[0] - self.yx0[0]
1131 y_stop = origin[0] + shape[0] - self.yx0[0]
1132 x_start = origin[1] - self.yx0[1]
1133 x_stop = origin[1] + shape[1] - self.yx0[1]
1134 y_index = slice(y_start, y_stop)
1135 x_index = slice(x_start, x_stop)
1136 return y_index, x_index
1138 def _get_sliced(self, indices: Any, value: Image | None = None) -> Image:
1139 """Select a subset of an image
1141 Parameters
1142 ----------
1143 indices:
1144 The indices to select a subsection of the image.
1145 The spectral index can either be a tuple of indices,
1146 a slice of indices, or a single index used to select a
1147 single-band 2D image.
1148 The spatial index (if present) is a `Box`.
1150 value:
1151 The value used to set this slice of the image.
1152 This allows the single `_get_sliced` method to be used for
1153 both getting a slice of an image and setting it.
1155 Returns
1156 -------
1157 result: Image | np.ndarray
1158 The resulting image obtained by selecting subsets of the iamge
1159 based on the `indices`.
1160 """
1161 if not isinstance(indices, tuple):
1162 indices = (indices,)
1164 if self.is_multiband:
1165 if self._is_spectral_index(indices[0]):
1166 if len(indices) > 1 and indices[1] in self.bands:
1167 # The indices are all band names,
1168 # so use them all as a spectral indices
1169 bands = indices
1170 spectral_index = self.spectral_indices(bands)
1171 y_index = x_index = slice(None)
1172 elif self._is_spectral_index(indices[0]):
1173 # The first index is a spectral index
1174 spectral_index = self.spectral_indices(indices[0])
1175 if isinstance(spectral_index, slice):
1176 bands = self.bands[spectral_index]
1177 elif len(spectral_index) == 1:
1178 bands = ()
1179 spectral_index = spectral_index[0] # type: ignore
1180 else:
1181 bands = tuple(self.bands[idx] for idx in spectral_index)
1182 indices = indices[1:]
1183 if len(indices) == 1:
1184 # The spatial index must be a bounding box
1185 if not isinstance(indices[0], Box):
1186 raise IndexError(f"Expected a Box for the spatial index but got {indices[1]}")
1187 y_index, x_index = self._get_box_slices(indices[0])
1188 elif len(indices) == 0:
1189 y_index = x_index = slice(None)
1190 else:
1191 raise IndexError(f"Too many spatial indices, expeected a Box bot got {indices}")
1192 full_index = (spectral_index, y_index, x_index)
1193 elif isinstance(indices[0], Box):
1194 bands = self.bands
1195 y_index, x_index = self._get_box_slices(indices[0])
1196 full_index = (slice(None), y_index, x_index)
1197 else:
1198 error = f"3D images can only be indexed by spectral indices or bounding boxes, got {indices}"
1199 raise IndexError(error)
1200 else:
1201 if len(indices) != 1 or not isinstance(indices[0], Box):
1202 raise IndexError(f"2D images can only be sliced by bounding box, got {indices}")
1203 bands = ()
1204 y_index, x_index = self._get_box_slices(indices[0])
1205 full_index = (y_index, x_index) # type: ignore
1207 y0 = y_index.start
1208 if y0 is None:
1209 y0 = 0
1211 x0 = x_index.start
1212 if x0 is None:
1213 x0 = 0
1215 if value is None:
1216 # This is a getter,
1217 # so return an image with the data sliced properly
1218 yx0 = (y0 + self.yx0[0], x0 + self.yx0[1])
1220 data = self.data[full_index]
1222 if data.ndim == 2:
1223 # Only a single band was selected, so return that band
1224 return Image(data, yx0=yx0)
1225 return Image(data, bands=bands, yx0=yx0)
1227 # Set the data
1228 self._data[full_index] = value.data
1229 return self
1231 def overlapped_slices(self, bbox: Box) -> tuple[tuple[slice, ...], tuple[slice, ...]]:
1232 """Get the slices needed to insert this image into a bounding box.
1234 Parameters
1235 ----------
1236 bbox:
1237 The region to insert this image into.
1239 Returns
1240 -------
1241 overlap:
1242 The slice of this image and the slice of the `bbox` required to
1243 insert the overlapping portion of this image.
1245 """
1246 overlap = self.bbox.overlapped_slices(bbox)
1247 if self.is_multiband:
1248 overlap = (slice(None),) + overlap[0], (slice(None),) + overlap[1]
1249 return overlap
1251 def __getitem__(self, indices: Any) -> Image:
1252 """Get the subset of an image
1254 Parameters
1255 ----------
1256 indices:
1257 The indices to select a subsection of the image.
1259 Returns
1260 -------
1261 result:
1262 The resulting image obtained by selecting subsets of the iamge
1263 based on the `indices`.
1264 """
1265 return self._get_sliced(indices)
1267 def __setitem__(self, indices, value: Image) -> Image:
1268 """Set a subset of an image to a given value
1270 Parameters
1271 ----------
1272 indices:
1273 The indices to select a subsection of the image.
1274 value:
1275 The value to use for the subset of the image.
1277 Returns
1278 -------
1279 result:
1280 The resulting image obtained by selecting subsets of the image
1281 based on the `indices`.
1282 """
1283 return self._get_sliced(indices, value)
1286def _operate_on_images(image1: Image, image2: Image | ScalarLike, op: Callable) -> Image:
1287 """Perform an operation on two images, that may or may not be spectrally
1288 and spatially aligned.
1290 Parameters
1291 ----------
1292 image1:
1293 The image on the LHS of the operation
1294 image2:
1295 The image on the RHS of the operation
1296 op:
1297 The operation used to combine the images.
1299 Returns
1300 -------
1301 image:
1302 The resulting combined image.
1303 """
1304 if type(image2) in ScalarTypes:
1305 return image1.copy_with(data=op(image1.data, image2))
1306 image2 = cast(Image, image2)
1307 if image1.bands == image2.bands and image1.bbox == image2.bbox:
1308 # The images perfectly overlap, so just combine their results
1309 with np.errstate(divide="ignore", invalid="ignore"):
1310 result = op(image1.data, image2.data)
1311 return Image(result, bands=image1.bands, yx0=image1.yx0)
1313 if op != operator.add and op != operator.sub and image1.bands != image2.bands:
1314 msg = "Images with different bands can only be combined using addition and subtraction, "
1315 msg += f"got {op}, with bands {image1.bands}, {image2.bands}"
1316 raise ValueError(msg)
1318 # Use all of the bands in the first image
1319 bands = image1.bands
1320 # Add on any bands from the second image not contained in the first image
1321 bands = bands + tuple(band for band in image2.bands if band not in bands)
1322 # Create a box that contains both images
1323 bbox = image1.bbox | image2.bbox
1324 # Create an image that will contain both images
1325 if len(bands) > 0:
1326 shape = (len(bands),) + bbox.shape
1327 else:
1328 shape = bbox.shape
1330 if op == operator.add or op == operator.sub:
1331 dtype = get_combined_dtype(image1, image2)
1332 result = Image(np.zeros(shape, dtype=dtype), bands=bands, yx0=cast(tuple[int, int], bbox.origin))
1333 # Add the first image in place
1334 image1.insert_into(result, operator.add)
1335 # Use the operator to insert the second image
1336 image2.insert_into(result, op)
1337 else:
1338 # Project both images into the full bbox
1339 image1 = image1.project(bbox=bbox)
1340 image2 = image2.project(bbox=bbox)
1341 result = op(image1, image2)
1342 return result
1345def insert_image(
1346 main_image: Image,
1347 sub_image: Image,
1348 op: Callable = operator.add,
1349) -> Image:
1350 """Insert one image into another image
1352 Parameters
1353 ----------
1354 main_image:
1355 The image that will have `sub_image` insertd.
1356 sub_image:
1357 The image that is inserted into `main_image`.
1358 op:
1359 The operator to use for insertion
1360 (addition, subtraction, multiplication, etc.).
1362 Returns
1363 -------
1364 main_image: Image
1365 The `main_image`, with the `sub_image` inserted in place.
1366 """
1367 if len(main_image.bands) == 0 and len(sub_image.bands) == 0:
1368 slices = sub_image.matched_slices(main_image.bbox)
1369 image_slices = slices[1]
1370 self_slices = slices[0]
1371 else:
1372 band_indices = sub_image.matched_spectral_indices(main_image)
1373 slices = sub_image.matched_slices(main_image.bbox)
1374 image_slices = (band_indices[0],) + slices[1] # type: ignore
1375 self_slices = (band_indices[1],) + slices[0] # type: ignore
1377 main_image._data[image_slices] = op(main_image.data[image_slices], sub_image.data[self_slices])
1378 return main_image