#!/usr/bin/env python
# -*- coding: utf-8 -*-
# ------------------------------------------------------------------------------
# Copyright (c) 2022 Univ. Lyon, ENS de Lyon, UCB Lyon 1, CNRS, INRAe, Inria
# All rights reserved.
# This file is part of the TimageTK library, and is released under the "GPLv3"
# license. Please see the LICENSE.md file that should have been included as
# part of this package.
# ------------------------------------------------------------------------------
"""Base image data structure."""
import copy as cp
from collections.abc import Iterable
import numpy as np
import vt
from pint import UnitRegistry
from timagetk.bin.logger import get_logger
from timagetk.components.metadata import IMAGE_MD_TAGS
from timagetk.components.metadata import ImageMetadata
from timagetk.components.metadata import Metadata
from timagetk.components.metadata import ProcessMetadata
from timagetk.third_party.vt_image import vtImage
from timagetk.util import _to_list
from timagetk.util import check_dimensionality
from timagetk.util import check_type
from timagetk.util import compute_extent
from timagetk.util import compute_voxelsize
from timagetk.util import dimensionality_test
#: Default axes order for 2D image.
DEFAULT_AXES_2D = "YX"
#: Default axes order for 3D image.
DEFAULT_AXES_3D = "ZYX"
#: Dictionary of default axes order for 2D image.
DEFAULT_AXIS_ORDER_2D = {ax: ax_id for ax_id, ax in enumerate(DEFAULT_AXES_2D)}
#: Dictionary of default axes order for 3D image.
DEFAULT_AXIS_ORDER_3D = {ax: ax_id for ax_id, ax in enumerate(DEFAULT_AXES_3D)}
#: Default voxels size for 2D image.
DEFAULT_VXS_2D = [1.0, 1.0]
#: Default voxels size for 3D image.
DEFAULT_VXS_3D = [1.0, 1.0, 1.0]
#: Default origin for 2D image.
DEFAULT_ORIG_2D = [0, 0]
#: Default origin for 3D image.
DEFAULT_ORIG_3D = [0, 0, 0]
#: Default image size unit is micrometer (1e-6).
DEFAULT_SIZE_UNIT = 1e-6
#: Valid data type dictionary for images, for details about numpy types, see here https://docs.scipy.org/doc/numpy-1.13.0/user/basics.types.html
DICT_TYPES = {0: np.uint8, 1: np.int8, 2: np.uint16, 3: np.int16, 4: np.uint32,
5: np.int32, 6: np.uint64, 7: np.int64, 8: np.float32,
9: np.float64, 10: np.float_, 11: np.complex64, 12: np.complex128,
13: np.complex_, 'uint8': np.uint8, 'uint16': np.uint16,
'ushort': np.uint16, 'uint32': np.uint32, 'uint64': np.uint64,
'uint': np.uint64, 'ulonglong': np.uint64, 'int8': np.int8,
'int16': np.int16, 'short': np.int16, 'int32': np.int32,
'int64': np.int64, 'int': np.int64, 'longlong': np.int64,
'float16': np.float16, 'float32': np.float32,
'single': np.float32, 'float64': np.float64, 'double': np.float64,
'float': np.float64, 'float128': np.float_,
'longdouble': np.float_, 'longfloat': np.float_,
'complex64': np.complex64, 'singlecomplex': np.complex64,
'complex128': np.complex128, 'cdouble': np.complex128,
'cfloat': np.complex128, 'complex': np.complex128,
'complex256': np.complex_, 'clongdouble': np.complex_,
'clongfloat': np.complex_, 'longcomplex': np.complex_}
# For details about numpy types, see https://docs.scipy.org/doc/numpy-1.13.0/user/basics.types.html
#: List of valid data type names.
AVAIL_TYPES = sorted([k for k in DICT_TYPES if isinstance(k, str)])
#: Default data type.
DEFAULT_DTYPE = DICT_TYPES[0]
#: List of protected image attribute or property.
PROTECT_PPTY = ['shape', 'min', 'max', 'mean']
DEFAULT_VALUE_MSG = "Set '{}' to default value: {}"
log = get_logger(__name__)
[docs]
def default_origin(input_array):
"""Return the default origin depending on the array dimensionality.
Parameters
----------
input_array : numpy.ndarray
2D or 3D array defining an image, *e.g.* intensity or labelled image.
Returns
-------
list
Default origin coordinates.
"""
if input_array.ndim == 2:
orig = DEFAULT_ORIG_2D
else:
orig = DEFAULT_ORIG_3D
log.debug(DEFAULT_VALUE_MSG.format("origin", orig))
return orig
[docs]
def default_voxelsize(input_array):
"""Return the default pixel or voxel size depending on the array dimensionality.
Parameters
----------
input_array : numpy.ndarray
2D or 3D array defining an image, *e.g.* intensity or labelled image.
Returns
-------
list
Default voxel size.
"""
if input_array.ndim == 2:
vxs = DEFAULT_VXS_2D
else:
vxs = DEFAULT_VXS_3D
log.debug(DEFAULT_VALUE_MSG.format("voxelsize", vxs))
return vxs
[docs]
def default_axes_order(input_array):
"""Return the default physical axes order depending on the array dimensionality.
Parameters
----------
input_array : numpy.ndarray
2D or 3D array defining an image, *e.g.* intensity or labelled image.
Returns
-------
str
Default physical axes ordering.
"""
if input_array.ndim == 2:
axo = DEFAULT_AXES_2D
else:
axo = DEFAULT_AXES_3D
log.debug(DEFAULT_VALUE_MSG.format("axes order", axo))
return axo
[docs]
class SpatialImage(np.ndarray):
"""Data structure of 2D and 3D images.
A ``SpatialImage`` gathers a numpy array and some metadata (such as voxel size, physical extent, origin, type, etc.).
Through a ``numpy.ndarray`` inheritance, all usual operations on
`numpy.ndarray <http://docs.scipy.org/doc/numpy-1.10.1/reference/generated/numpy.ndarray.html>`_
objects (sum, product, transposition, etc.) are available.
All image processing operations are performed on this data structure, that
is also used to solve inputs (read) and outputs (write).
We `subclass <https://docs.scipy.org/doc/numpy-1.14.0/user/basics.subclassing.html>`_
ndarray using view casting in the ``__new__`` section.
View casting is the standard ``ndarray`` mechanism by which you take an
``ndarray`` of any subclass, and return a view of the array as another
(specified) subclass, here a ``SpatialImage``.
.. warning::
* Modification of a property, when possible, will change the object, not return a new one!
* Attributes ``shape`` and ``ndim`` are defined by the input array.
* Specifying ``dtype`` at object instantiation will change the array dtype to the given one.
* Images tags can be set or updated by a metadata dictionary.
Attributes
----------
_extent : list of float
The physical extent of the image, sorted as (planes,) rows & columns [(Z,) Y, X].
_origin : list of int
The coordinates of origin, sorted as (planes,) rows & columns [(Z,) Y, X].
_metadata : dict
A dictionary of metadata.
_voxelsize : list of float
The size of the pixels or voxels (2D/3D), sorted as (planes,) rows & columns [(Z,) Y, X].
axes_order : dict
Dictionary of axes order, *e.g.* {'x': 2, 'y': 1, 'z': 0} by default for 3D image.
filepath : str
The directory where to find the image, if any.
filename : str
The file name of the image.
metadata_image : ImageMetadata
A self generated list of attribute, contains basic image information such as: 'shape', 'ndim', 'dtype', 'origin', 'voxelsize' & 'extent'.
See ``ImageMetadata`` class docstring in ``timagetk.components.metadata``.
metadata_acquisition : Metadata
Acquisition metadata is a list of attributes related to the measurement settings, the machine setup, or even
used equipments and parameters.
Use the 'scan_info' key in ``metadata`` dict to pass them to constructor.
metadata_processing : ProcessMetadata
A self generated hierarchical dictionary of the algorithms applied to the image.
It contains process names, the parameters use as well as their values and, if any, the sub-process called with
the same degree of info.
.. todo::
* Add a ``metric_unit`` attribute and associated methods.
* Defines a ``personal_metadata`` attribute to store other metadata ?
"""
[docs]
def __new__(cls, array, origin=None, voxelsize=None, dtype=None, metadata=None, **kwargs):
"""Data structure of 2D and 3D images.
Parameters
----------
array : numpy.ndarray or timagetk.SpatialImage
2D (YX) or 3D (ZYX) array defining an image, *e.g.* intensity or labelled image.
origin : list, optional
(Z)YX coordinates of the origin in the image, default is ``[0, 0]`` if 2D or ``[0, 0, 0]`` if 3D.
voxelsize : list, optional.
(Z)YX image voxelsize, default is ``[1.0, 1.0]`` if 2D or ``[1.0, 1.0, 1.0]`` if 3D.
dtype : str, optional
If given, should be in ``AVAIL_TYPES``, and will modify the input `array` type.
metadata : dict, optional
Dictionary of image metadata, default is an empty dict.
Other Parameters
----------------
unit : float
The size unit, in International System of Units (meters), associated to the image.
axes_order : dict
Dictionary of axes orders, *e.g.* {'x': 2, 'y': 1, 'z': 0} by default for 3D image.
Returns
-------
timagetk.SpatialImage
Image with metadata.
See Also
--------
timagetk.components.spatial_image.AVAIL_TYPES
Notes
-----
By default numpy arrays are C-contiguous with index starting at ``0``.
We follow this convention, meaning that 2D images are column then row sorted.
For 3D images planes are thus the first index, *i.e.* to access plane 3 of a ``SpatialImage`` named ``img``,
do: ``plane = img[2]``.
To use the array dtype (from numpy), leave ``dtype`` to ``None``.
Example
-------
>>> import numpy as np
>>> from timagetk.io import imread
>>> from timagetk.io.util import shared_dataset
>>> from timagetk import SpatialImage
>>> image_path = shared_dataset("p58", "intensity")[0]
>>> spim = imread(image_path)
>>> im = SpatialImage(spim, voxelsize=[1., 1., 1.])
>>> from timagetk.array_util import random_spatial_image
>>> img = random_spatial_image((3, 5, 5), voxelsize=[1., 0.5, 0.5], unit=1e-6)
>>> isinstance(img, SpatialImage) # Test `SpatialImage` type
True
>>> isinstance(img, np.ndarray) # Test `np.ndarray` type => INHERITANCE
True
>>> print(img.voxelsize) # Access the voxel size information
[1., 0.5, 0.5]
>>> print(img.unit) # Access the unit information
micrometer
>>> print(img.shape) # Access the plane information
(3, 5, 5)
>>> img[0] # Access the first plane
"""
from timagetk.third_party.vt_converter import _new_from_vtimage
from timagetk.third_party.vt_converter import attr_from_spatialimage
# - Accept ``SpatialImage`` or ``numpy.ndarray`` instances:
log.debug(f"``SpatialImage`` got a ``{type(array)}`` as input!")
kwargs['origin'] = origin
kwargs['voxelsize'] = voxelsize
kwargs['metadata'] = metadata
if isinstance(array, SpatialImage):
# - Override keyword arguments if defined as attributes in `array` instance:
kwargs |= attr_from_spatialimage(array, **kwargs)
elif isinstance(array, vt.vtImage):
# Do NOT assign to `origin` & `metadata` variables to avoid overriding input!
# Furthermore, vtImage does not (yet?) have those attributes.
array, voxelsize, _, _ = _new_from_vtimage(array)
kwargs['voxelsize'] = voxelsize
else:
# - Test input array type, should be a numpy array :
check_type(array, 'array', np.ndarray)
# - Take care of image attributes:
from timagetk.io.image import default_image_attributes
kwargs |= default_image_attributes(array, log_undef=log.warning, **kwargs)
# ----------------------------------------------------------------------
# ARRAY:
# ----------------------------------------------------------------------
# - Test input array dimensionality, should be of dimension 2 or 3:
try:
assert array.ndim in [2, 3]
except AssertionError:
msg = "Input 'array' must be 2D or 3D!\n"
for attr in ['ndim', 'shape', 'dtype']:
try:
msg += "Got '{}': {}\n".format(attr, getattr(array, attr))
except:
pass
raise ValueError(msg)
# - View casting (see class documentation) of the array as a SpatialImage subclass:
# NOTE: it converts the array to given dtype if not None:
# This call ``__array_finalize__``!
if array.flags.c_contiguous:
obj = np.asarray(array, dtype=dtype).view(cls)
else:
obj = np.ascontiguousarray(array, dtype=dtype).view(cls)
# ----------------------------------------------------------------------
# METADATA:
# ----------------------------------------------------------------------
metadata = kwargs['metadata']
# - Get acquisition metadata using 'scan_info' key:
scan_info = metadata.pop('scan_info', {})
obj.metadata_acquisition = Metadata(scan_info)
if obj.metadata_acquisition is None:
msg = "Could not set 'metadata_acquisition' attribute!"
raise AttributeError(msg)
# - Get ImageMetadata:
image_metadata = ImageMetadata(array, **kwargs)
# - Update the metadata dictionary for image IMAGE_MD_TAGS & extent:
metadata = image_metadata.update_metadata(metadata)
if 'acquisition_date' in metadata:
image_metadata.update({'acquisition_date': metadata['acquisition_date']})
obj.metadata_image = image_metadata
obj._metadata = {}
# - Set 'axes_order', 'voxelsize', 'origin', 'extent' & 'unit' hidden attributes:
obj._axes_order = image_metadata['axes_order']
obj._voxelsize = _to_list(image_metadata['voxelsize'])
obj._origin = _to_list(image_metadata['origin'])
obj._extent = _to_list(image_metadata['extent'])
obj._unit = image_metadata['unit']
obj._unit_reg = UnitRegistry() # initialize unit registry
# - Get processing metadata :
obj.metadata_processing = {}
# -- Set the 'class' metadata if not defined:
try:
assert 'class' in obj.metadata_processing
except AssertionError:
obj.metadata_processing.update({'class': 'SpatialImage'})
else:
# DO NOT try to change returned class, eg. 'LabelledImage', ...
pass
# - Defines some attributes:
# -- File name and path, if known:
obj.filename = kwargs.get('filename', metadata.get('filename', ''))
obj.filepath = kwargs.get('filepath', metadata.get('filepath', ''))
return obj
[docs]
def __array_finalize__(self, obj):
"""Mechanism provided by numpy to allow subclasses instanciation.
Parameters
----------
obj : object
The object returned by the ``__new__`` method.
Examples
--------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((3, 5, 5), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> type(img)
timagetk.components.spatial_image.SpatialImage
>>> # Taking the first z-slice automatically return a SpatialImage instance:
>>> img2d = img[0]
>>> type(img2d)
timagetk.components.spatial_image.SpatialImage
"""
self._axes_order = getattr(obj, '_axes_order', None)
self._origin = getattr(obj, '_origin', None)
self._voxelsize = getattr(obj, '_voxelsize', None)
self._extent = getattr(obj, '_extent', None)
self._unit_reg = getattr(obj, '_unit_reg', None)
self._unit = getattr(obj, '_unit', None)
self.filename = getattr(obj, 'filename', None)
self.filepath = getattr(obj, 'filepath', None)
self._metadata = getattr(obj, '_metadata', None)
self.metadata_image = getattr(obj, 'metadata_image', None)
self.metadata_acquisition = getattr(obj, 'metadata_acquisition', None)
self.metadata_processing = getattr(obj, 'metadata_processing', None)
# Case where we are doing view casting of a SpatialImage and we changed the array:
# md_img = self.metadata_image
# View Casting DEBUG:
# if md_img is not None:
# print("Got image metadata: {}".format(md_img.get_dict()))
# print("With md_img['shape']: {}".format(md_img['shape']))
# print("With obj.shape: {}".format(obj.shape))
# -----
# if md_img is not None and md_img['shape'] != obj.shape:
# # Update image metadata:
# self.metadata_image.get_from_image(obj)
# self.metadata = self.metadata_image.update_metadata(self.metadata)
# # Reset some attributes
# for attr in ['min', 'max', 'mean']:
# setattr(self, '_' + attr, None)
[docs]
def __str__(self):
"""Print information about the object."""
msg = "SpatialImage object with following metadata:\n"
md = self.metadata_image.get_dict()
msg += '\n'.join([' - {}: {}'.format(k, v) for k, v in md.items()])
return msg
[docs]
def __setstate__(self, state, /):
"""Restore the state of the object for deserialization with pickle"""
sp_state, array_state = state
super().__setstate__(array_state)
self.__dict__.update(sp_state[0])
self.metadata_image = ImageMetadata(self)
self.metadata_image.update(sp_state[1])
self.metadata_acquisition = Metadata()
self.metadata_acquisition.update(sp_state[2])
self._unit_reg = UnitRegistry()
[docs]
def __reduce_ex__(self, __protocol):
"""Returns reconstruction info for pickle serialization"""
reconstructor, reconstructor_args, state = super().__reduce_ex__(__protocol)
state_dict = {
key: val for key, val in self.__dict__.items()
if key not in ["metadata_image", "metadata_acquisition", "_unit_reg"]
}
sp_state = (state_dict, self.metadata_image.get_dict(), self.metadata_acquisition.get_dict())
new_state = (sp_state, state)
return reconstructor, reconstructor_args, new_state
# --------------------------------------------------------------------------
#
# Conversion:
#
# --------------------------------------------------------------------------
[docs]
def to_vtimage(self):
"""Return the same object converted as ``vtImage`` instance.
Returns
-------
vtImage
A vtImage object.
Examples
--------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((3, 4, 5), voxelsize=[1., 0.51, 0.5], dtype='uint8')
>>> vt_im = img.to_vtimage()
>>> from vt import vtImage
>>> isinstance(vt_im, vtImage)
True
>>> vt_im.shape()
>>> vt_im.spacing()
>>> # Initialize a random (uint8) 2D SpatialImage:
>>> img = random_spatial_image((4, 5), voxelsize=[0.51, 0.5], dtype='uint8')
>>> vt_im = img.to_vtimage()
>>> from vt import vtImage
>>> isinstance(vt_im, vtImage)
True
>>> vt_im.shape()
>>> vt_im.spacing()
"""
return vtImage(self, self.voxelsize)
[docs]
def get_array(self, dtype=None):
"""Return a copy of the image as a ``numpy.ndarray`` object.
Parameters
----------
dtype : str or numpy.dtype, optional
If specified, default is ``None``, change the dtype of the returned array.
Returns
-------
numpy.ndarray
Numpy array linked to the object (Z)YX sorted.
Notes
-----
Changing the `dtype` of the array force the return of a copy.
Examples
--------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((3, 5, 5), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> from timagetk import SpatialImage
>>> isinstance(img, SpatialImage)
True
>>> # Get the array linked to the SpatialImage:
>>> array = img.get_array()
>>> isinstance(array, SpatialImage)
False
>>> import numpy as np
>>> isinstance(array, np.ndarray)
True
>>> # Return the pointer:
>>> array.fill(0)
>>> print(array)
>>> print(img.get_array())
"""
return np.asarray(self, dtype=dtype).copy()
[docs]
def get_xyz_array(self, dtype=None):
"""Get an XY(Z) sorted ``numpy.ndarray`` from a ``SpatialImage``.
Parameters
----------
dtype : str or numpy.dtype, optional
If specified, default is ``None``, change the dtype of the returned array.
Returns
-------
numpy.ndarray
Numpy array linked to the SpatialImage, XY(Z) sorted.
.. deprecated::
This method should be removed for the next major release
Examples
--------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((3, 5, 5), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> from timagetk import SpatialImage
>>> isinstance(img, SpatialImage)
True
>>> # Get the array linked to the SpatialImage:
>>> array = img.get_xyz_array()
>>> print(array.shape)
"""
import warnings
warnings.warn("This method maintain legacy compatibility with timagetk version anterior to 3!",
DeprecationWarning)
ax = self.axes_order_dict
if self.is2D():
return self.get_array(dtype).transpose((ax['X'], ax['Y']))
else:
return self.get_array(dtype).transpose((ax['X'], ax['Y'], ax['Z']))
# --------------------------------------------------------------------------
#
# SpatialImage properties:
#
# --------------------------------------------------------------------------
@property
def axes_order(self):
"""Get ``SpatialImage`` physical axes ordering.
Returns
-------
str
``SpatialImage`` physical axes ordering.
Example
-------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((3, 5, 5), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> print(img.axes_order)
ZYX
"""
return cp.copy(self._axes_order)
@property
def axes_order_dict(self):
"""Get ``SpatialImage`` physical axes ordering.
Returns
-------
dict
``SpatialImage`` physical axes ordering.
Example
-------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((3, 5, 5), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> print(img.axes_order_dict)
{'Z': 0, 'Y': 1, 'X': 2}
"""
return {ax: axi for axi, ax in enumerate(self._axes_order)}
@property
def extent(self):
"""Get ``SpatialImage`` physical extent.
Returns
-------
list
``SpatialImage`` physical extent.
Notes
-----
It is related to the array shape and image voxelsize.
Example
-------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((3, 5, 5), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> print(img.extent)
[3.0, 2.5, 2.5]
"""
return cp.copy(self._extent)
@extent.setter
def extent(self, img_extent):
"""Set ``SpatialImage`` physical extent.
Parameters
----------
img_extent : list
``SpatialImage`` new physical extent.
Notes
-----
Do NOT use this method except if you know EXACTLY what you are doing!
This will change voxelsize based on array shape.
Metadata are updated according to the new physical extent and voxelsize.
Example
-------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((3, 5, 5),voxelsize=[1., 0.5, 0.5],dtype='uint8')
>>> print(img.voxelsize)
[1.0, 0.5, 0.5]
>>> img.extent = [10., 10., 10.]
>>> print(img.extent)
[10., 10., 10.]
>>> print(img.voxelsize)
[3.3333333333333335, 2.0, 2.0]
"""
log.warning("You are changing a physical property of the image!")
log.warning("This should not happen!")
log.warning("You better know what you are doing!")
dimensionality_test(self.ndim, img_extent)
# - Update 'extent' hidden attribute:
self._extent = _to_list(img_extent)
# - Update 'voxelsize' hidden attribute:
self._voxelsize = compute_voxelsize(self.shape, self.extent)
# - Update 'extent' & 'voxelsize' metadata:
self.metadata_image['extent'] = self.extent
self.metadata_image['voxelsize'] = self.voxelsize
log.info("Set extent to '{}'".format(self.extent))
log.info("Changed voxelsize to '{}'".format(self.voxelsize))
return
@property
def metadata(self):
"""Get ``SpatialImage`` metadata dictionary.
Returns
-------
dict
Metadata dictionary.
Example
-------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((3, 5, 5),voxelsize=[1., 0.5, 0.5],dtype='uint8',metadata={'info': "test OK"})
>>> print(img.metadata['info'])
'test OK'
>>> print(img.metadata)
{'info': 'test OK',
'shape': (3, 5, 5),
'ndim': 3,
'dtype': dtype('uint8'),
'origin': [0, 0, 0],
'voxelsize': [1.0, 0.5, 0.5],
'get_extent': [3.0, 2.5, 2.5]}
"""
md = {}
if isinstance(self.metadata_image, ImageMetadata):
md.update(self.metadata_image.get_dict())
if isinstance(self.metadata_acquisition, Metadata):
acq = self.metadata_acquisition.get_dict()
if acq != {}:
md['acquisition'] = acq
if isinstance(self.metadata_processing, ProcessMetadata):
proc = self.metadata_processing.get_dict()
if proc != {}:
md['processing'] = proc
self._metadata.update(md)
return cp.copy(self._metadata)
@metadata.setter
def metadata(self, img_md):
"""Update ``SpatialImage`` metadata dictionary.
Example
-------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((3, 5, 5),voxelsize=[1., 0.5, 0.5],dtype='uint8')
>>> print(img.metadata)
{'shape': (3, 5, 5),
'ndim': 3,
'dtype': dtype('uint8'),
'origin': [0, 0, 0],
'voxelsize': [1.0, 0.5, 0.5],
'get_extent': [3.0, 2.5, 2.5]}
>>> img.metadata = {'info': "test OK"}
>>> print(img.metadata)
{'shape': (3, 5, 5),
'ndim': 3,
'dtype': dtype('uint8'),
'origin': [0, 0, 0],
'voxelsize': [1.0, 0.5, 0.5],
'get_extent': [3.0, 2.5, 2.5],
'info': 'test OK'}
"""
self._metadata.update(img_md)
@property
def origin(self):
"""Get ``SpatialImage`` origin.
Returns
-------
list
``SpatialImage`` origin coordinates.
Example
-------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((3, 5, 5),voxelsize=[1., 0.5, 0.5],dtype='uint8')
>>> print(img.origin)
[0, 0, 0]
"""
return cp.copy(self._origin)
@origin.setter
def origin(self, img_origin):
"""Set ``SpatialImage`` origin.
Given list should be of same length than the image dimensionality.
Parameters
----------
img_origin : list
``SpatialImage`` origin coordinates.
Example
-------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((3, 5, 5),voxelsize=[1., 0.5, 0.5],dtype='uint8')
>>> print(img.origin)
[0, 0, 0]
>>> img.origin = [2, 2, 2]
>>> print(img.origin)
[2, 2, 2]
"""
log.warning("You are changing a physical property of the image!")
log.warning("This should not happen!")
log.warning("You better know what you are doing!")
dimensionality_test(self.ndim, img_origin)
# - Update hidden attribute 'origin':
self._origin = _to_list(img_origin)
# - Update hidden attribute metadata key 'origin':
self.metadata_image['origin'] = self.origin
log.info("Set origin to '{}'".format(self.origin))
return
@property
def voxelsize(self):
"""Get ``SpatialImage`` voxelsize, sorted as planes, columns & rows.
Returns
-------
list
``SpatialImage`` voxelsize
Example
-------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((3, 5, 5),voxelsize=[1., 0.5, 0.5],dtype='uint8')
>>> print(img.voxelsize)
[1., 0.5, 0.5]
"""
return cp.copy(self._voxelsize)
@voxelsize.setter
def voxelsize(self, img_vxs):
"""Change ``SpatialImage`` voxelsize by preserving 'extent'.
Parameters
----------
img_vxs : list
Image new voxelsize.
Notes
-----
Do NOT use this method except if you know EXACTLY what you are doing!
This will change voxelsize based on array shape.
Metadata are updated according to the new physical extent and voxelsize.
Example
-------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((3, 5, 5), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> print(img.voxelsize)
[1., 0.5, 0.5]
>>> print(img.extent)
[3.0, 2.5, 2.5]
>>> print(img.shape)
(3, 5, 5)
>>> # - Change img voxelsize:
>>> img.voxelsize = [0.5, 0.5, 0.5]
>>> print(img.voxelsize)
[0.5, 0.5, 0.5]
>>> print(img.extent)
[1.5, 2.5, 2.5]
>>> print(img.shape)
(3, 5, 5)
"""
log.warning("You are changing a physical property of the image!")
log.warning("This should not happen!")
log.warning("You better know what you are doing!")
dimensionality_test(self.ndim, img_vxs)
# - Update 'voxelsize' hidden attribute:
self._voxelsize = _to_list(img_vxs)
# - Update 'extent' hidden attribute:
self._extent = compute_extent(self.voxelsize, self.shape)
# - Update 'extent' & 'voxelsize' metadata:
self.metadata_image['voxelsize'] = self.voxelsize
self.metadata_image['extent'] = self.extent
log.info("Set voxelsize to '{}'".format(self.voxelsize))
log.info("Changed extent to '{}'".format(self.extent))
return
@property
def unit(self):
"""Return the unit in International System of Units (SI).
Returns
-------
float
The unit in International System of Units (SI).
"""
return cp.copy(self._unit)
@unit.setter
def unit(self, unit):
"""Set the unit of the image size, in International System of Units (SI).
Parameters
----------
unit : int or float
The unit, in meters, associated to the image.
Examples
--------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((3, 5, 5), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> print(img.get_unit()) # default unit is micrometer
'micrometer'
>>> img.unit = 1e-3 # change to millimeter
>>> print(img.get_unit())
'millimeter'
"""
self._unit = unit
# --------------------------------------------------------------------------
#
# SpatialImage methods :
#
# --------------------------------------------------------------------------
[docs]
def is_isometric(self):
"""Test image isometry.
Tests if the voxelsize values are the same for every axis or not.
Returns
-------
bool
``True`` if the image is isometric, else ``False``.
Examples
--------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((3, 5, 5),voxelsize=[1., 0.5, 0.5],dtype='uint8')
>>> img.is_isometric()
False
>>> from timagetk.algorithms.resample import isometric_resampling
>>> # Transform to small isometric 3D SpatialImage:
>>> img_iso = isometric_resampling(img, method='max')
>>> img_iso.is_isometric()
True
>>> img_iso.voxelsize
[1.0, 1.0, 1.0]
"""
vxs = [self.voxelsize[0]] * self.ndim
return np.allclose(vxs, self.voxelsize, atol=1.e-9)
[docs]
def is2D(self):
"""Return True if the SpatialImage is 2D, else ``False``.
Examples
--------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((3, 5, 5), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> img.is2D()
False
>>> # Take the first z-slice:
>>> img2d = img.get_slice(0, 'z') # Equivalent to img[0]
>>> img2d.is2D()
True
"""
return self.ndim == 2
[docs]
def is3D(self):
"""Return True if the SpatialImage is 3D, else ``False``.
Examples
--------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((3, 5, 5), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> img.is3D()
True
"""
return self.ndim == 3
# --------------------------------------------------------------------------
#
# GETTERs & SETTERs:
#
# --------------------------------------------------------------------------
def _get_unit(self):
return (self._unit * self._unit_reg.meter).to_compact()
[docs]
def get_unit(self, short=False):
"""Get the human-readable image unit in International System of Units (SI).
Examples
--------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((3, 5, 5), voxelsize=[1., 0.5, 0.5])
>>> print(img.get_unit())
micrometer
>>> print(img.get_unit(short=True))
µm
"""
if short:
return f"{self._get_unit().u:P~}"
else:
return f"{self._get_unit().u}"
[docs]
def get_dtype(self):
"""Get the human-readable image type.
Examples
--------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((3, 5, 5), voxelsize=[1., 0.5, 0.5])
>>> img.get_dtype()
'uint8'
"""
return self.dtype.name
[docs]
def get_extent(self, axis=None):
"""Return the extent of the image or of the selected axis, if any.
Parameters
----------
axis : {'x', 'y', 'z'}, optional
If specified, axis for which the extent should be returned.
Returns
-------
list of float or float
The extent of the image or of the selected axis, if any.
The list is sorted like the physical axes order.
Notes
-----
The extent of an image is related to its shape and voxelsize.
See Also
--------
timagetk.util.compute_extent
Examples
--------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((3, 5, 5), voxelsize=[1., 0.5, 0.5])
>>> img.get_extent()
[2., 2., 2.]
>>> img.get_extent('x')
2.0
>>> img.get_extent('y')
2.0
>>> img.get_extent('z')
2.0
"""
if axis is None:
return cp.copy(self.extent)
else:
return cp.copy(self.extent[self.axes_order_dict[axis.upper()]])
[docs]
def get_shape(self, axis=None):
"""Return the shape of the image or of the selected axis, if any.
Parameters
----------
axis : {'x', 'y', 'z'}, optional
If specified, axis for which the shape should be returned.
Returns
-------
list of int or int
The shape of the image or of the selected axis, if any.
The list is sorted like the physical axes order.
Examples
--------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((3, 5, 5), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> img.get_shape()
[3, 5, 5]
>>> img.get_shape('z')
3
"""
if axis is None:
return list(cp.copy(self.shape))
else:
try:
assert axis.lower() in ['x', 'y', 'z']
except AssertionError:
raise ValueError("Parameter `axis` take values in {'x', 'y', 'z'}!")
except AttributeError:
raise ValueError("Parameter `axis` should be string!")
return cp.copy(self.shape[self.axes_order_dict[axis.upper()]])
[docs]
def get_voxelsize(self, axis=None):
"""Return the voxel-size of the image or of the selected axis, if any.
Parameters
----------
axis : {'x', 'y', 'z'}, optional
If specified, axis for which the voxel-size should be returned.
Returns
-------
list of float or float
The voxel-size of the image or of the selected axis, if any.
The list is sorted like the physical axes order.
Examples
--------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((3, 5, 5), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> img.get_voxelsize()
[1., 0.5, 0.5]
>>> img.get_voxelsize('z')
1.0
"""
if axis is None:
return list(cp.copy(self.voxelsize))
else:
try:
assert axis.lower() in ['x', 'y', 'z']
except AssertionError:
raise ValueError(
"Parameter `axis` take values in {'x', 'y', 'z'}!")
except AttributeError:
raise ValueError("Parameter `axis` should be string!")
return cp.copy(self.voxelsize[self.axes_order_dict[axis.upper()]])
[docs]
def set_slice(self, slice_id, array, axis='z'):
"""Set the values of an axis slice to given array.
Parameters
----------
slice_id : int
Slice number to modify.
array : numpy.ndarray or int
Values that replace those at given slice id, if an integer make an array with this value.
axis : str in {'x', 'y', 'z'}, optional
Axis to use for slicing, default to 'z'.
Examples
--------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((3, 5, 5), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> sl = img.get_slice(0, axis='z')
>>> print(sl.get_array())
>>> new_sl = random_spatial_image((5, 5), voxelsize=[0.5, 0.5], dtype='uint8')
>>> img.set_slice(0, new_sl, axis='z')
>>> print(new_sl)
>>> print(img.get_slice(0, 'z'))
"""
sl = self.get_slice(slice_id, axis)
if isinstance(array, int):
array = np.ones_like(sl) * array
try:
assert sl.shape == array.shape
except AssertionError:
msg = "Given array {} has a wrong shape for {}-slice {} {}!"
raise ValueError(msg.format(array.shape, axis, slice_id, sl.shape))
ax_id = self.axes_order_dict[axis.upper()]
if ax_id == 2:
self[:, :, slice_id] = array
elif ax_id == 1:
self[:, slice_id, :] = array
elif ax_id == 0:
self[slice_id, :, :] = array
else:
raise ValueError("Parameter `axis` take values in {'x', 'y', 'z'}!")
def _get_slice(self, slice_id, axis):
"""Take one or more slices along given axis.
Parameters
----------
slice_id : int or Iterable
Slice to return.
axis : int or str in {'x', 'y', 'z'}
Axis to use for slicing, default to 'z'.
Returns
-------
numpy.ndarray
2D/3D array with only the required slice(s).
dict
New axes ordering to use with slice(s).
numpy.ndarray
Origin of the slice(s).
numpy.ndarray
Voxel size of the slice(s).
dict
Dictionary of slice(s) metadata.
Raises
------
ValueError
If the image is not 3D and ``axis='z'``.
If ``slice_id`` does not exists, *i.e.* should satisfy: ``0 < slice_id < max(len(axis))``.
"""
# -- Correspondence between `axis` string and `axis_idx` index:
if isinstance(axis, str):
axis_idx = self.axes_order_dict[axis.upper()]
else:
axis_idx = cp.copy(axis)
axis = {v: k for k, v in self.axes_order_dict.items()}[axis]
# -- Check `axis` is defined for this image:
try:
assert axis.upper() in set(self.axes_order_dict.keys())
except AssertionError:
raise ValueError(f"Could not find required {axis.upper()}-axis in this image!")
# -- If a list with just one integer, ``slice_id`` should be an integer:
if isinstance(slice_id, Iterable) and len(slice_id) == 1:
slice_id = slice_id[0]
# -- Make sure this slice exists:
max_slice = self.shape[axis_idx]
if not isinstance(slice_id, Iterable):
try:
assert slice_id <= max_slice
except AssertionError:
msg = "Required {}-slice ({}) does not exists (max: {})"
raise ValueError(msg.format(axis, slice_id, max_slice))
else:
try:
assert np.max(slice_id) <= max_slice
except AssertionError:
msg = "Required {}-slice range ({}-{}) does not exists (max: {})"
raise ValueError(msg.format(axis, np.min(slice_id), np.max(slice_id), max_slice))
# -- Take the slice(s) in the array:
arr = np.take(self.get_array(), slice_id, axis=axis_idx)
vxs = self.voxelsize
ori = self.origin
md = self.metadata
if not isinstance(slice_id, Iterable):
# If ``slice_id`` is an integer, one axis will "disappear"...
# its info should thus be removed:
new_axes_order = ''.join([ax for ax, _ in self.axes_order_dict.items() if ax.upper() != axis.upper()])
# Remove value corresponding to returned slice axis:
vxs.pop(axis_idx)
ori.pop(axis_idx)
# Clear image metadata tags, they will be recomputed by SpatialImage.__new__
for attr in IMAGE_MD_TAGS + ['extent']:
md.pop(attr)
else:
new_axes_order = self.axes_order.copy()
# Assume we have a constant compression ratio along the selected axis:
ratio = self.get_shape(axis) / arr.shape[axis_idx]
vxs[axis_idx] = vxs[axis_idx] * ratio
for attr in ['voxelsize', 'extent']:
md.pop(attr)
log.debug("Re-slicing image with:")
log.debug(f"New array shape: {arr.shape}")
log.debug(f"New ``axes_order``: {new_axes_order}")
log.debug(f"New ``ori``: {ori}")
log.debug(f"New ``vxs``: {vxs}")
log.debug(f"New ``md``: {md}")
return arr, new_axes_order, ori, vxs, md
[docs]
def get_slice(self, slice_id, axis='z'):
"""Return a SpatialImage with only one slice for given axis.
Parameters
----------
slice_id : int
Slice to return.
axis : int or str in {'x', 'y', 'z'}, optional
Axis to use for slicing, default to 'z'.
Returns
-------
SpatialImage
2D SpatialImage with only the required slice.
Raises
------
ValueError
If the image is not 3D and ``axis='z'``.
If ``slice_id`` does not exist, *i.e.* should satisfy: ``0 < slice_id < max(len(axis))``.
Examples
--------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((3, 5, 5), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> print(img.axes_order_dict)
{'Z': 0, 'Y': 1, 'X': 2}
>>> print(img.get_shape())
[3, 5, 5]
>>> # Taking an existing z-slice from a 3D image works fine:
>>> img_z = img.get_slice(1, 'z')
>>> print(img_z.axes_order_dict)
{'Y': 0, 'X': 1}
>>> print(img_z.get_shape())
[5, 5]
>>> # Taking an existing x-slice from a 3D image works fine:
>>> img_x = img.get_slice(3, 'x')
>>> print(img_x.axes_order_dict)
{'Z': 0, 'Y': 1}
>>> print(img_x.get_shape())
[3, 5]
>>> # Down-sampling x-axis of a 3D image:
>>> nx = img.get_shape('x')
>>> img_ds_x2 = img.get_slice(range(0, nx, 2), 'x')
>>> # Taking an NON-existing z-slice from a 3D image raises an error:
>>> img.get_slice(50, 'z')
>>> # Taking a z-slice from a 2D image raises an error:
>>> img_z.get_slice(5, 'z')
"""
arr, axes_order, ori, vxs, md = self._get_slice(slice_id, axis)
return SpatialImage(arr, origin=ori, voxelsize=vxs, metadata=md, unit=self.unit, axes_order=axes_order)
[docs]
@staticmethod
def get_available_types():
"""Print the available bits type dictionary.
See Also
--------
DICT_TYPES : dict
Dictionary of available bits type.
Returns
-------
dict
Dictionary of available bits type.
"""
return DICT_TYPES
[docs]
def get_pixel(self, idx):
"""Get ``SpatialImage`` pixel value at given array coordinates.
Parameters
----------
idx : list
Indices as list of integers.
Returns
-------
value
Pixel value.
Raises
------
TypeError
If the given `idx` is not a list.
ValueError
If the number of `idx` is wrong, should be the image dimensionality.
If the `idx` coordinates are not within the image boundaries.
Example
-------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((3, 5, 5), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> img.get_pixel([1,1])
"""
check_type(idx, 'idx', list)
# Check we have the right number of `idx`:
ndim = self.ndim
try:
assert len(idx) == ndim
except AssertionError:
msg = "Input `idx` must have a length equal to the image dimensionality!"
raise ValueError(msg)
# Check the `ids` are within the image boundaries:
sh = self.shape
ids = range(0, ndim)
try:
assert all((idx[i] >= 0) & (idx[i] < sh[i] + 1) for i in ids)
except AssertionError:
raise ValueError(f"Input `idx` {idx} are not within the image shape {self.shape}!")
if ndim == 2:
pix_val = self[idx[0], idx[1]]
else:
pix_val = self[idx[0], idx[1], idx[2]]
return pix_val
[docs]
def set_pixel(self, idx, value):
"""Change ``SpatialImage`` pixel|voxel value at given array coordinates.
Parameters
----------
idx : list
Array coordinates as list of integers.
value : array.dtype
New value for the pixel|voxel, should be compatible with the array `dtype`.
Raises
------
TypeError
If the given `idx` is not a list.
ValueError
If the number of `idx` is wrong, should be the image dimensionality.
If the `idx` coordinates are not within the image boundaries.
Example
-------
>>> import numpy as np
>>> from timagetk import SpatialImage
>>> test_array = np.ones((3,3), dtype=np.uint8)
>>> image = SpatialImage(test_array)
>>> image.set_pixel([1, 1], 2)
>>> image.get_array()
array([[1, 1, 1],
[1, 2, 1],
[1, 1, 1]], dtype=uint8)
>>> # Trying to set a value not compatible with the array.dtype:
>>> image.set_pixel([2, 2], 0.5)
>>> image.get_array()
array([[1, 1, 1],
[1, 2, 1],
[1, 1, 0]], dtype=uint8)
>>> # Trying to set a value not compatible with the array.dtype:
>>> image.set_pixel([0, 0], -6)
>>> image.get_array()
array([[250, 1, 1],
[ 1, 2, 1],
[ 1, 1, 0]], dtype=uint8)
>>> image.set_pixel([0, 0], -6)
"""
check_type(idx, 'idx', list)
# Check we have the right number of `idx`:
ndim = self.ndim
try:
assert len(idx) == ndim
except AssertionError:
msg = "Input `idx` must have a length equal to the image dimensionality!"
raise ValueError(msg)
# Check the `ids` are within the image boundaries:
sh = self.shape
ids = range(0, ndim)
try:
assert all((idx[i] >= 0) & (idx[i] < sh[i] + 1) for i in ids)
except AssertionError:
raise ValueError(f"Input `idx` {idx} are not within the image shape {self.shape}!")
if ndim == 2:
self[idx[0], idx[1]] = value
else:
self[idx[0], idx[1], idx[2]] = value
return
[docs]
def get_region(self, region):
"""Extract a region using list of start & stop indexes.
There should be two values per dimension in `region`, *e.g.* ``region=[5, 8, 5, 8]`` for a 2D image.
If the image is 3D and, in one dimension the start and stop indexes only differ by one (one layer of voxels), the returned image will be transformed to 2D!
Parameters
----------
region : list
Indexes as list of integers, *e.g.* `[y-start, y-stop, x-start, x-stop]` for a 2D image.
Returns
-------
SpatialImage
Output image.
Raises
------
TypeError
If the given `region` is not a list.
ValueError
If the number of indexes in `region` is wrong, should be twice the image dimensionality.
If the `region` coordinates are not within the array boundaries.
Example
-------
>>> from timagetk import SpatialImage
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 2D SpatialImage:
>>> img = random_spatial_image((5, 5), voxelsize=[0.5, 0.5], dtype='uint8')
>>> region = [1, 3, 1, -1] # y-start, y-stop, x-start, x-stop
>>> out_img = img.get_region(region)
>>> isinstance(out_img, SpatialImage)
True
>>> out_img.get_array()
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((5, 5, 5), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> region = [1, 2, 1, 3, 1, 3] # z-start, z-stop, y-start, y-stop, x-start, x-stop
>>> out_img = img.get_region(region)
>>> isinstance(out_img, SpatialImage)
True
>>> out_img.is2D()
True
"""
# TODO: use slice instead of 'indices' list?
check_type(region, 'idx', list)
# Check we have the right number of `idx`:
try:
assert len(region) == 2 * self.ndim
except AssertionError:
msg = "Input `region` must have a length equal to twice the number of dimension of the image!"
raise ValueError(msg)
region_start = region[::2]
region_stop = region[1::2]
sh = self.shape
# Check the "stop" values and replace negative ones:
for n, r in enumerate(region_stop):
if r < 0:
region_stop[n] = sh[n] + r
try:
assert 0 < region_stop[n] <= sh[n] - 1
except AssertionError:
raise ValueError(f"Input `region` {region} are not within the image shape {self.shape}!")
# Check the "start" values:
for n, r in enumerate(region_start):
try:
assert 0 <= r < sh[n] - 1
except AssertionError:
raise ValueError(f"Input `region` {region} are not within the image shape {self.shape}!")
bbox = [slice(region_start[i], region_stop[i] + 1) for i in range(self.ndim)]
tmp_arr, tmp_vox = self.get_array(), self.voxelsize
reg_val = tmp_arr[tuple(bbox)]
out_sp_img = SpatialImage(reg_val, voxelsize=tmp_vox)
if 1 in out_sp_img.shape:
out_sp_img = out_sp_img.to_2d()
return out_sp_img
[docs]
def set_region(self, region, value):
"""Replace a region of the image by given value.
Parameters
----------
region : list
Indices as list of integers.
value : array.dtype or ndarray
New value for the selected pixels, type of ``SpatialImage`` array.
Returns
-------
SpatialImage
``SpatialImage`` instance.
Raises
------
TypeError
If the given `region` is not a list.
ValueError
If the number of indexes in `region` is wrong, should be twice the image dimensionality.
If the `region` coordinates are not within the array boundaries.
If the `region` and `value` shape mismatch.
If the given `value` is not of the same dtype as the array or not a ndarray.
Example
-------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((5, 5), voxelsize=[0.5, 0.5], dtype='uint8')
>>> region = [1, 3, 1, 3]
>>> img.set_region(region, 3)
>>> img.get_array()
"""
check_type(region, 'idx', list)
# Check we have the right number of `idx`:
ndim = self.ndim
try:
assert len(region) == 2 * ndim
except AssertionError:
msg = "Input 'region' must have a length equal to twice the number of dimension of the image!"
raise ValueError(msg)
# Check the `region` is within the image boundaries:
sh = self.shape
ids = range(0, 2 * ndim, 2)
try:
assert all((region[i] >= 0) & (region[i + 1] < sh[i // 2] + 1) for i in ids)
except AssertionError:
err = "Input 'idx' {} are not within the image shape {}!"
raise ValueError(err.format(region, self.shape))
if isinstance(value, np.ndarray):
region_shape = [region[i + 1] - region[i] for i in ids]
try:
np.testing.assert_equal(region_shape, value.shape)
except AssertionError:
msg = "Given region shape ({}) and value shape ({}) mismatch!"
raise ValueError(msg.format(region_shape, value.shape))
dtype = self.dtype
if isinstance(value, np.ndarray):
if ndim == 2:
self[region[0]: region[1], region[2]: region[3]] = value[:, :].astype(dtype)
else:
self[region[0]: region[1], region[2]: region[3], region[4]: region[5]] = value[:, :, :].astype(dtype)
elif isinstance(value, dtype):
if ndim == 2:
self[region[0]: region[1], region[2]: region[3]] = value
else:
self[region[0]: region[1], region[2]: region[3], region[4]: region[5]] = value
else:
msg = "The given `value` is not of the same dtype than the array or not a ndarray"
raise ValueError(msg)
return
# --------------------------------------------------------------------------
#
# SpatialImage transformation functions:
#
# --------------------------------------------------------------------------
[docs]
def astype(self, dtype, **kwargs):
"""Copy of the SpatialImage with updated data type.
Parameters
----------
dtype : str
New type of data to apply.
Returns
-------
SpatialImage
Image with the new data type.
Notes
-----
Keyword arguments are passed to ``numpy.astype()`` method.
Examples
--------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((5, 5), voxelsize=[0.5, 0.5], dtype='uint8')
>>> # Check the type:
>>> img.dtype
dtype('uint8')
>>> print(img)
SpatialImage object with following metadata:
- shape: (5, 5)
- ndim: 2
- dtype: uint8
- axes_order: YX
- voxelsize: [0.5, 0.5]
- unit: 1e-06
- origin: [0, 0]
- extent: [2.0, 2.0]
- acquisition_date: None
>>> # Convert it to 16bits unsigned integers:
>>> img16 = img.astype('uint16')
>>> # Check the type:
>>> img16.dtype
dtype('uint16')
>>> print(img16)
SpatialImage object with following metadata:
- shape: (5, 5)
- ndim: 2
- dtype: uint16
- axes_order: YX
- voxelsize: [0.5, 0.5]
- unit: 1e-06
- origin: [0, 0]
- extent: [2.0, 2.0]
- acquisition_date: None
"""
# Get the list of attribute to preserve:
img_attr = ['axes_order', 'origin', 'voxelsize', 'metadata', 'not_a_label', 'background', 'unit', 'filename']
img_kwargs = {attr: getattr(self, attr, None) for attr in img_attr}
# - Convert the numpy array:
array = self.get_array().astype(dtype, **kwargs)
return SpatialImage(array, **img_kwargs)
[docs]
def crop(self, start, stop, axis='z'):
"""Crop the selected axis and preserve data between given start & stop slice id.
Parameters
----------
start : int
The slice id to start cropping.
stop : int
The slice id to stop cropping.
axis : {'x', 'y', 'z'} or int, optional
Name of the axis to crop, or its index. By default, the ``z`` axis.
Returns
-------
SpatialImage
The cropped image
Raises
------
ValueError
If the image is not 3D and ``axis='z'``.
If ``start`` and/or ``stop`` indexes does not exist on this axis.
Examples
--------
>>> from timagetk.io.util import shared_dataset
>>> from timagetk.io import imread
>>> from timagetk.visu.stack import orthogonal_view
>>> img = imread(shared_dataset("p58", "intensity")[0])
>>> # Example - Get the first half of the image on the x-axis
>>> cropped_x = img.crop(0, img.get_shape('x')//2, 'x' )
>>> print(img.get_extent('x'))
91.74656248092651
>>> print(cropped_x.get_extent('x'))
45.67296123504639
>>> print(img)
SpatialImage object with following metadata:
- shape: (160, 230, 230)
- ndim: 3
- dtype: uint8
- axes_order: ZYX
- voxelsize: [0.40064001083374023, 0.40064001083374023, 0.40064001083374023]
- unit: 1e-06
- origin: [0, 0, 0]
- extent: [63.7017617225647, 91.74656248092651, 91.74656248092651]
- acquisition_date: None
>>> print(cropped_x)
SpatialImage object with following metadata:
- shape: (160, 230, 115)
- ndim: 3
- dtype: uint8
- axes_order: ZYX
- voxelsize: [0.40064001083374023, 0.40064001083374023, 0.40064001083374023]
- unit: 1e-06
- origin: [0, 0, 0]
- extent: [63.7017617225647, 91.74656248092651, 45.67296123504639]
- acquisition_date: None
>>> orthogonal_view(img)
>>> orthogonal_view(cropped_x)
"""
axis_shape = self.get_shape(axis=axis)
try:
assert 0 <= start <= stop - 1 < stop <= axis_shape - 1
except AssertionError:
raise ValueError("Check given copping limits!")
if stop - start == 1:
return self.get_slice(start, axis=axis)
if isinstance(axis, str):
axis_id = self.axes_order_dict[axis.upper()]
else:
axis_id = cp.copy(axis)
if axis_id > self.ndim:
raise ValueError(f"Can not crop axis {axis_id} on a {self.ndim}D image!")
arr = np.take(self.get_array(), range(start, stop, 1), axis=axis_id)
# Get the list of attribute to preserve:
img_attr = ['axes_order', 'origin', 'voxelsize', 'metadata', 'not_a_label', 'background', 'unit', 'filename']
img_kwargs = {attr: getattr(self, attr, None) for attr in img_attr}
log.debug(f"Array shape: {arr.shape}")
log.debug(f"Used image attributes: {img_kwargs}")
return SpatialImage(arr, **img_kwargs)
[docs]
def to_2d(self):
"""Convert a 3D SpatialImage to a 2D SpatialImage.
Conversion is possible only if there is a "flat" axis (*i.e.* with only one slice along this axis).
Returns
-------
SpatialImage
The 2D SpatialImage.
Examples
--------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage with one z-slice:
>>> img = random_spatial_image((1, 5, 5), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> img.to_2d()
>>> # Initialize a random (uint8) 3D SpatialImage with one y-slice:
>>> img = random_spatial_image((5, 1, 5), voxelsize=[0.5, 1., 0.5], dtype='uint8')
>>> img.to_2d()
>>> # Initialize a random (uint8) 3D SpatialImage with one x-slice:
>>> img = random_spatial_image((5, 5, 1), voxelsize=[0.5, 0.5, 1.], dtype='uint8')
>>> img.to_2d()
"""
from timagetk.io.image import default_image_attributes
if self.is3D() and 1 in self.shape:
# Get the array and its shape:
shape, array = self.get_shape(), self.get_array()
# Get attribute 'voxelsize':
voxelsize = self.voxelsize
# Get attributes 'origin' & 'metadata':
ori, md = self.origin, self.metadata
# Search axe to squeeze:
if self.get_shape('x') == 1:
ax2squeeze = self.axes_order_dict['X']
elif self.get_shape('y') == 1:
ax2squeeze = self.axes_order_dict['Y']
else:
ax2squeeze = self.axes_order_dict['Z']
# Squeeze found axe:
new_arr = np.squeeze(array, axis=(ax2squeeze,))
# Edit attributes 'voxelsize' & 'origin':
voxelsize.pop(ax2squeeze)
ori.pop(ax2squeeze)
return SpatialImage(new_arr, **default_image_attributes(new_arr, voxelsize=voxelsize, origin=ori, metadata=md))
elif self.is2D():
log.error("This image is already a 2D image!")
return None
else:
log.error('This 3D SpatialImage can not be reshaped to 2D.')
return None
[docs]
def to_3d(self, flat_axis='z'):
"""Convert a 2D SpatialImage to a 3D SpatialImage.
Parameters
----------
flat_axis : {'x', 'y', 'z'}, optional
Flat axis to add to the 2D image. By default, `'z'` axis.
Returns
-------
SpatialImage
The 3D image with a flat axis.
Notes
-----
Obtained 3D SpatialImage has a "flat" axis (i.e. with only one slice along this axis).
The voxelsize of this axis is set to ``1.``.
The origin of this axis is set to ``0``.
Examples
--------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 2D SpatialImage:
>>> img = random_spatial_image((5, 5), voxelsize=[0.5, 0.5], dtype='uint8')
>>> z_flat = img.to_3d(flat_axis='z')
>>> z_flat.get_shape()
>>> y_flat = img.to_3d(flat_axis='y')
>>> y_flat.get_shape()
>>> x_flat = img.to_3d(flat_axis='x')
>>> x_flat.get_shape()
>>> x_flat.origin
>>> x_flat.voxelsize
"""
from timagetk.io.image import default_image_attributes
if not self.is2D():
log.error('This SpatialImage is not 2D.')
return None
# Get the 'id' of the axis to reverse:
axis_id = DEFAULT_AXIS_ORDER_3D[flat_axis.upper()]
# Get the array and its shape:
shape, array = self.get_shape(), self.get_array()
# Get attribute 'voxelsize':
voxelsize = self.voxelsize
# Get attributes 'origin' & 'metadata':
ori, md = self.origin, self.metadata
# Update attributes:
voxelsize.insert(axis_id, 1.)
ori.insert(axis_id, 0)
shape.insert(axis_id, 1)
# Add new flat "axe":
new_arr = np.reshape(array, shape)
return SpatialImage(new_arr, **default_image_attributes(new_arr, voxelsize=voxelsize, origin=ori, metadata=md))
def _new_order(self, *axes) -> list:
"""Return the list of integers indexing axes to re-order."""
if len(axes) == 1 and isinstance(axes[0], str) and len(axes[0]) == self.ndim:
# If axes are given as str, translate them to "axe ids":
return [self.axes_order_dict[axe.upper()] for i, axe in enumerate(axes[0])]
elif isinstance(axes, (list, tuple)) and all(isinstance(axe, str) for axe in axes):
# If axes are given as list of str, translate them to "axe ids":
return [self.axes_order_dict[axe.upper()] for axe in axes]
else:
raise ValueError(f"Could not make sense of given `axes`: {axes}")
[docs]
def transpose(self, *axes):
"""Permute image axes to given order, reverse by default.
Parameters
----------
axes : list of int or list of str, optional
By default, reverse the dimensions, otherwise permute the axes according to the values given.
Returns
-------
SpatialImage
The image with permuted axes.
Examples
--------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((3, 4, 5), voxelsize=[0.51, 0.52, 0.53], dtype='uint8')
>>> img.metadata
{'shape': (3, 4, 5),
'ndim': 3,
'dtype': dtype('uint8'),
'axes_order': 'ZYX',
'voxelsize': [0.51, 0.52, 0.53],
'unit': 1e-06,
'origin': (0.0, 0.0, 0.0),
'extent': [1.02, 1.56, 2.12],
'acquisition_date': None}
>>> img_t = img.transpose()
>>> img_t.metadata
{'shape': (5, 4, 3),
'ndim': 3,
'dtype': dtype('uint8'),
'axes_order': 'XYZ',
'voxelsize': [0.53, 0.52, 0.51],
'unit': 1e-06,
'origin': [0.0, 0.0, 0.0],
'extent': [2.12, 1.56, 1.02],
'acquisition_date': None}
>>> print(img_t.shape) # `shape` attribute has been updated
(5, 4, 3)
>>> print(img_t.voxelsize) # `voxelsize` attribute has been updated
[0.53, 0.52, 0.51]
>>> # Method `transpose` accept axes names as input
>>> img_t = img.transpose('xzy')
>>> print(img_t.shape)
(5, 3, 4)
>>> # Method `transpose` accept a tuple of axes names as input
>>> img_t = img.transpose('x', 'z', 'y')
>>> print(img_t.shape)
(5, 3, 4)
>>> # Initialize a random (uint8) 2D SpatialImage:
>>> img = random_spatial_image((4, 5), voxelsize=[0.52, 0.53], dtype='uint8')
>>> img_t = img.transpose()
>>> # Transpose update the shape attribute of the image (here reversed):
>>> print(img_t.shape)
(5, 4)
"""
if axes == ():
# If no axes order is specified, reverse the axes.
axes = list(range(self.ndim))[::-1]
elif isinstance(axes, (list, tuple)) and all(isinstance(ax, int) for ax in axes):
pass
else:
axes = self._new_order(*axes)
# Uses ``axes`` to reorder the `axes_order`, `origin` and `voxelsize` attributes:
ori = [self.origin[axe] for axe in axes]
vxs = [self.voxelsize[axe] for axe in axes]
axo = ''.join([self.axes_order[axe] for axe in axes])
# Get the list of attribute to preserve:
img_attr = ['metadata', 'not_a_label', 'background', 'unit', 'filename']
img_kwargs = {attr: getattr(self, attr, None) for attr in img_attr}
# Transpose the array:
arr = np.transpose(self.get_array(), axes)
return SpatialImage(arr, origin=ori, voxelsize=vxs, axes_order=axo, **img_kwargs)
[docs]
def invert_axis(self, axis):
"""Revert given axis.
Parameters
----------
axis : {'x', 'y', 'z'}
Axis to invert, can be either 'x', 'y' or 'z' (if 3D).
Returns
-------
SpatialImage
Image with reverted array for selected axis.
Raises
------
ValueError
If given ``axis`` is not in {'x', 'y', 'z'} for 3D images or not in {'x', 'y'} for 2D images.
Example
-------
>>> from timagetk.array_util import random_spatial_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> img = random_spatial_image((3, 5, 5), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> print(img.get_slice(0, "z").get_array())
>>> inv_img = img.invert_axis(axis='z')
>>> print(inv_img.get_slice(0, "z").get_array())
"""
# Get the list of attribute to preserve:
img_attr = ['axes_order', 'origin', 'voxelsize', 'metadata', 'not_a_label', 'background', 'unit', 'filename']
img_kwargs = {attr: getattr(self, attr, None) for attr in img_attr}
if self.is2D():
arr = self._invert_2d(axis)
else:
arr = self._invert_3d(axis)
return SpatialImage(arr, **img_kwargs)
def _invert_3d(self, axis):
"""Revert x, y, or z axis of the given 3D array.
Parameters
----------
axis : str in {'x', 'y', 'z'}
Array axis to revert.
Returns
-------
numpy.ndarray
Reverted array for selected axis.
Raises
------
ValueError
If given ``axis`` is not in {'x', 'y', 'z'}
"""
if self.axes_order_dict[axis.upper()] == 2:
return self.get_array()[:, :, ::-1]
elif self.axes_order_dict[axis.upper()] == 1:
return self.get_array()[:, ::-1, :]
elif self.axes_order_dict[axis.upper()] == 0:
return self.get_array()[::-1, :, :]
else:
raise ValueError("Unknown axis '{}' for a 3D array.".format(axis))
def _invert_2d(self, axis):
"""Revert x-axis or y-axis of the given 2D array.
Parameters
----------
axis : str in {'x', 'y'}
Array axis to revert.
Returns
-------
numpy.ndarray
Reverted array for selected axis.
Raises
------
ValueError
If given ``axis`` is not in {'x', 'y'}
"""
if self.axes_order_dict[axis.upper()] == 1:
return self.get_array()[:, ::-1]
elif self.axes_order_dict[axis.upper()] == 0:
return self.get_array()[::-1, :]
else:
raise ValueError("Unknown axis '{}' for a 2D array.".format(axis))
