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