#!/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.
# ------------------------------------------------------------------------------
"""MultiChannel images module."""
import copy as cp
import numpy as np
from pint import UnitRegistry
from skimage import img_as_float
from skimage import img_as_ubyte
from skimage import img_as_uint
from timagetk.bin.logger import get_logger
from timagetk.components.metadata import ImageMetadata
from timagetk.components.metadata import Metadata
from timagetk.components.metadata import ProcessMetadata
from timagetk.components.spatial_image import DEFAULT_SIZE_UNIT
from timagetk.util import _to_list
from timagetk.util import compute_extent
from timagetk.util import compute_voxelsize
from timagetk.util import dimensionality_test
log = get_logger(__name__)
DEF_4D_AXES = "ZYXC"
DEFAULT_AXIS_ORDER_4D = {ax: ax_id for ax_id, ax in enumerate(DEF_4D_AXES)}
DEF_CHANNEL_COLORS = ['red', 'green', 'blue', 'yellow', 'cyan', 'magenta']
[docs]
def check_images_channels(image, channel_names):
"""Check that there is the same number of images and channel names.
Parameters
----------
image : list of str or list of timagetk.SpatialImage
A list of image names, file names or instances.
channel_names : list of str
A list of channel names.
Raises
------
ValueError
If the length of the two lists is not the same.
"""
try:
assert len(image) == len(channel_names)
except AssertionError:
msg = f"Please provide the same number of images ({len(image)}) and channel names ({len(channel_names)})!"
raise ValueError(msg)
return
[docs]
def check_image_shape(image):
"""Check that all images have the same shape.
Parameters
----------
image : dict of timagetk.SpatialImage
A dictionary of ``SpatialImage``.
Raises
------
ValueError
If not all images are of same shape.
"""
ref_ch = list(image.keys())[0]
ref_shape = image[ref_ch].shape
try:
assert all(ref_shape == img.shape for _, img in image.items())
except AssertionError:
msg = "All images should have the same shape!\n"
for ch, img in image.items():
msg += f" - Channel '{ch}' has a shape of {img.shape}\n"
raise ValueError(msg)
return
[docs]
def check_image_voxelsize(image, rtol=1e-5):
"""Check that all images have the same voxel size.
Parameters
----------
image : dict of SpatialImage
A dictionary of SpatialImage.
rtol : float, optional
The relative tolerance parameter, *i.e.* the maximum relative difference allowed.
Raises
------
ValueError
If not all images are of same voxel size.
"""
ref_ch = list(image.keys())[0]
ref_vxs = image[ref_ch].voxelsize
try:
assert all(np.allclose(ref_vxs, img.voxelsize, rtol=rtol) for _, img in image.items())
except AssertionError as err:
log.error(err)
msg = "All images should have the same voxelsize!\n"
for ch, img in image.items():
msg += f" - Channel '{ch}' has a voxelsize of {img.voxelsize}\n"
raise ValueError(msg)
return
[docs]
def check_image_axes_order(image):
"""Check that all images have the same ``axes_order`` attribute.
Parameters
----------
image : dict of SpatialImage
A dictionary of SpatialImage.
Raises
------
ValueError
If not all images have the same ``axes_order`` attribute.
"""
ref_ch = list(image.keys())[0]
ref_axes_order = image[ref_ch].axes_order_dict
try:
assert all(ref_axes_order == img.axes_order_dict for _, img in image.items())
except AssertionError as err:
log.error(err)
msg = "All images should have the same axes ordering!\n"
for ch, img in image.items():
msg += f" - Channel '{ch}' has an axes ordering of {img.voxelsize}\n"
raise ValueError(msg)
return
[docs]
def chnames_generator(n_imgs):
"""Generate channel names as `CH_i` with `i` in range(n_imgs).
Parameters
----------
n_imgs : int
Number of images or channels to generate names for.
Examples
--------
>>> from timagetk.components.multi_channel import chnames_generator
>>> chnames_generator(5)
['CH_0', 'CH_1', 'CH_2', 'CH_3', 'CH_4']
"""
return [f"CH_{i}" for i in range(n_imgs)]
from collections import UserDict
[docs]
class MultiChannelImage(UserDict):
"""Multichannel image data structure.
Attributes
----------
name : str
The name of the multichannel image.
channel_names : list of str
The list of channel names.
metadata : Metadata
A metadata instance for the multichannel image.
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].
voxelsize : list of float
The size of the pixels or voxels (2D/3D), sorted as (planes,) rows & columns [(Z,) Y, X].
ndim : int
The number of dimensions for each channel.
filepath : str
The directory where to find the image, if any.
filename : str
The file name of the image.
Notes
-----
TODO: Add a ``metric_unit`` attribute and associated methods to handle this!
Examples
--------
>>> from timagetk import MultiChannelImage
>>> from timagetk.io import imread
>>> im_url = "https://zenodo.org/record/3737795/files/qDII-CLV3-PIN1-PI-E35-LD-SAM4.czi"
>>> ch_names = ['DII-VENUS-N7', 'pPIN1:PIN1-GFP', 'Propidium Iodide', 'pRPS5a:TagBFP-SV40', 'pCLV3:mCherry-N7']
>>> # IO example #1 - `imread` method can read CZI and return a MultiChannelImage:
>>> im = imread(im_url, channel_names=ch_names)
>>> im.channel_names
>>> isinstance(im, MultiChannelImage)
True
>>> from timagetk.visu.stack import orthogonal_view
>>> orthogonal_view(im.to_blend_image())
>>> # IO example #2 - MultiChannelImage can be instantiated from a filename:
>>> im = MultiChannelImage(im_url, channel_names=ch_names, czi_pattern='.C.ZXY')
>>> print(im)
>>> # Access the channels as you would with a dictionary:
>>> print(im["DII-VENUS-N7"])
SpatialImage object with following metadata:
- shape: (65, 1024, 1024)
- ndim: 3
- dtype: uint16
- origin: [0, 0, 0]
- voxelsize: (6.771450737227882e-07, 2.075664593629475e-07, 2.075664593629475e-07)
- extent: [4.4014429791981236e-05, 0.00021254805438765823, 0.00021254805438765823]
>>> # Access the channels using the `get_channel` method:
>>> ahp6_img = im.get_channel("AHP6")
>>> # Get the name of the multichannel image:
>>> print(im.name)
qDII-CLV3-PIN1-PI-E35-LD-SAM4
>>> # Add or change the name of the multichannel image:
>>> im.set_name("my_experiment1_rep2")
>>> print(im.name)
my_experiment1_rep2
>>> from timagetk.visu.stack import channels_stack_browser
>>> channels_stack_browser(im, im.name)
"""
[docs]
def __init__(self, image, channel_names=None, metadata=None, **kwargs):
"""Constructor.
Parameters
----------
image : str or list of str or list of SpatialImage or dict
If an `str`, should be a multichannel image file name.
If a `list of str`, should be a list of single channel file names.
Lists of SpatialImage are also accepted, they should be of same shape.
If a ``dict``, keys are used as ``channel_names`` and values should be ``SpatialImage``.
channel_names : list of str, optional
List of names to use as channel names.
metadata : dict, optional
Dictionary of image metadata. An empty dictionary by default.
Other Parameters
----------------
name : str
Name of the multichannel image.
czi_pattern : str
Pattern used to read the CZI file, ``".C.ZXY."`` by default.
filepath : str
The directory where to find the image, if any.
filename : str
The file name of the image.
"""
from timagetk.components.image import assert_spatial_image
from timagetk.io.image import imread
if isinstance(image, (list, tuple)):
# Generates channel names if not defined:
if channel_names is None:
channel_names = chnames_generator(len(image))
# Check the number of images and channels:
self._check_images_channels(image, channel_names)
# Load the image if a string (filepath or URL) else assume it's a SpatialImage:
image = [imread(img) if isinstance(img, str) else img for img in image]
# Assert we now have a list of SpatialImage...
[assert_spatial_image(im, 'Channel image #{}') for im in image]
# Generate the {channel_name: image} dictionary:
image = dict(zip(channel_names, image))
else:
assert isinstance(image, dict)
channel_names = list(image.keys())
super().__init__(image)
# -- Make sure images are...
# ... of identical shapes...
self._check_image_shape(image)
# ... have the same voxel sizes...
self._check_image_voxelsize(image)
# ... and axes ordering
self._check_image_axes_order(image)
# -- Initialize attributes:
self.channel_names = None
self.n_ch = None
# - Physical image attributes:
self.voxelsize = None
self.shape = None
self.ndim = None
self.extent = None
self.origin = None
self.dtype = None
self.axes_order = None
# - File attributes:
self.filename = None
self.filepath = None
self.name = None
# - Metadata attributes:
self._metadata = Metadata({})
self.metadata_image = None
self.metadata_acquisition = None
self.metadata_processing = None
# -- Set attributes values:
self.update_attributes(image, channel_names, **kwargs)
def update_attributes(self, image, channel_names, **kwargs):
# -- Add attributes:
self.channel_names = channel_names
self.n_ch = len(self.channel_names)
ref_img = list(image.values())[0]
# - Physical image attributes:
self.voxelsize = ref_img.voxelsize
self.shape = ref_img.shape
self.extent = ref_img.extent
self.origin = ref_img.origin
self.dtype = ref_img.dtype
self.ndim = ref_img.ndim + 1
# - Defines axes order:
axes_order_dict = ref_img.axes_order_dict
# Add `channels` as last axes
axes_order_dict.update({'C': max(axes_order_dict.values()) + 1})
if self.ndim == 4 and axes_order_dict != DEFAULT_AXIS_ORDER_4D:
log.warning(f"You have a non-default 3D multichannel image axes ordering: {axes_order_dict}")
self.axes_order = ref_img.axes_order + 'C'
self.axes_order_dict = axes_order_dict
# - File attributes:
self.filename = kwargs.get("filename", '')
self.filepath = kwargs.get("filepath", '')
self.name = kwargs.get("name", self.filename)
# - Metadata attributes:
metadata = kwargs.get('metadata', self._metadata._dict)
self._metadata = Metadata(metadata)
self.metadata_image = {ch: ImageMetadata(im, origin=self.origin, voxelsize=self.voxelsize) for ch, im in
image.items()}
self.metadata_acquisition = {}
self.metadata_processing = {'class': 'MultiChannelImage'}
[docs]
def __str__(self):
msg = f"MultiChannelImage '{self.name}' has {len(self)} channels: "
msg += ", ".join(self.channel_names)
msg += "\n * "
msg += '\n * '.join(
[f'{ch_name}: {img}' for ch_name, img in self.items()])
return msg
@staticmethod
def _check_images_channels(image, channel_names):
"""Check that we have the same number of images and channel names.
Parameters
----------
image : list
A list of image name or object.
channel_names : list
A list of channel names.
Raises
------
ValueError
If the length of the two lists is not the same.
"""
check_images_channels(image, channel_names)
@staticmethod
def _check_image_shape(image):
"""Check that all images have the same shape.
Parameters
----------
image : dict of SpatialImage
A dictionary of SpatialImage.
Raises
------
ValueError
If not all images are of same shape.
"""
check_image_shape(image)
@staticmethod
def _check_image_voxelsize(image):
"""Check that all images have the same voxel-size.
Parameters
----------
image : list of SpatialImage
A list of SpatialImage.
Raises
------
ValueError
If not all images are of same voxel-size.
"""
check_image_voxelsize(image)
@staticmethod
def _check_image_axes_order(image):
"""Check that all images have the same ``axes_order`` attribute.
Parameters
----------
image : list of SpatialImage
A list of SpatialImage.
Raises
------
ValueError
If not all images have the same ``axes_order`` attribute.
"""
check_image_axes_order(image)
# --------------------------------------------------------------------------
#
# SpatialImage properties:
#
# --------------------------------------------------------------------------
@property
def metadata(self):
"""Get ``MultiChannelImage`` metadata dictionary.
Returns
-------
dict
Metadata dictionary.
"""
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 ``MultiChannelImage`` metadata dictionary."""
self._metadata.update(img_md)
# --------------------------------------------------------------------------
#
# Conversion:
#
# --------------------------------------------------------------------------
[docs]
def get_array(self, return_order=None):
"""Return a NumPy array.
Parameters
----------
return_order : str, optional
Order of the dimension for returned array, same as the ``axes_order`` attribute by default.
Returns
-------
numpy.ndarray
The array that correspond to the requested axes ordering.
Examples
--------
>>> # Example #1: Use a dummy multichannel image:
>>> from timagetk.array_util import dummy_multichannel_image
>>> img = dummy_multichannel_image()
>>> print(img.axes_order)
>>> arr = img.get_array()
>>> print(arr.shape) # same as the ``axes_order`` attribute
(5, 11, 11, 2)
>>> arr = img.get_array("CXYZ") # specify the axes return order
>>> print(arr.shape)
(2, 11, 11, 5)
>>> # Example #2: Use a real multichannel image (downloaded from Zenodo.org):
>>> from timagetk.io import imread
>>> im_fname = "https://zenodo.org/record/3737795/files/qDII-CLV3-PIN1-PI-E35-LD-SAM4.czi"
>>> ch_names = ['DII-VENUS-N7', 'pPIN1:PIN1-GFP', 'Propidium Iodide', 'pRPS5a:TagBFP-SV40', 'pCLV3:mCherry-N7']
>>> img = imread(im_fname, channel_names=ch_names)
>>> arr = img.get_array()
>>> print(arr.shape)
(40, 1024, 1024, 4)
>>> arr = img.get_array("CXYZ")
>>> print(arr.shape)
(4, 1024, 1024, 40)
"""
ch_names = self.channel_names
# - Stack the channel along a new axis (in last position):
arr = np.stack([self[ch].get_array() for ch in ch_names], axis=-1)
if return_order is None:
return_order = "".join(list(self.axes_order_dict.keys()))
axis = [self.axes_order_dict[c.upper()] for c in return_order]
log.debug(f"Found axes transposition to apply: {axis}")
arr = arr.transpose(axis)
return arr
# --------------------------------------------------------------------------
#
# Array methods :
#
# --------------------------------------------------------------------------
[docs]
def invert_axis(self, axis):
"""Revert given axis for all channels.
Parameters
----------
axis : str in {'x', 'y', 'z'}
Axis to invert, can be either 'x', 'y' or 'z' (if 3D)
Returns
-------
MultiChannelImage
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
>>> from timagetk import MultiChannelImage
>>> # Initialize a random (uint8) 3D SpatialImage with one z-slice:
>>> img_a = random_spatial_image((3, 5, 5), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> img_b = random_spatial_image((3, 5, 5), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> img = MultiChannelImage([img_a, img_b])
>>> img.get_channel_names()
['CH_0', 'CH_1']
>>> print(img_a.get_array())
>>> img.invert_axis('z')
>>> print(img['CH_0'].get_array())
"""
for ch in self.channel_names:
self[ch] = self[ch].invert_axis(axis)
return
[docs]
def is_isometric(self):
"""Test image isometry.
Tests if the voxelsize values are the same for every axis or not.
Returns
-------
is_iso : bool
True if the image is isometric, else ``False``.
Examples
--------
>>> from timagetk import MultiChannelImage
>>> from timagetk.io import imread
>>> im_fname = "https://zenodo.org/record/3737795/files/qDII-CLV3-PIN1-PI-E35-LD-SAM4.czi"
>>> ch_names = ['DII-VENUS-N7', 'pPIN1:PIN1-GFP', 'Propidium Iodide', 'pRPS5a:TagBFP-SV40', 'pCLV3:mCherry-N7']
>>> img = imread(im_fname, channel_names=ch_names)
>>> img.is_isometric()
False
"""
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 import MultiChannelImage
>>> from timagetk.io import imread
>>> im_fname = "https://zenodo.org/record/3737795/files/qDII-CLV3-PIN1-PI-E35-LD-SAM4.czi"
>>> ch_names = ['DII-VENUS-N7', 'pPIN1:PIN1-GFP', 'Propidium Iodide', 'pRPS5a:TagBFP-SV40', 'pCLV3:mCherry-N7']
>>> img = imread(im_fname, channel_names=ch_names)
>>> img.is2D()
False
>>> # Take the first z-slice:
>>> img2d = img.get_slice(0, 'z') # Equivalent to img[0]
>>> img2d.is2D()
True
"""
return self.get_channel(self.get_channel_names()[0]).ndim == 2
[docs]
def is3D(self):
"""Return True if the SpatialImage is 3D, else ``False``.
Examples
--------
>>> from timagetk import MultiChannelImage
>>> from timagetk.io import imread
>>> im_fname = "https://zenodo.org/record/3737795/files/qDII-CLV3-PIN1-PI-E35-LD-SAM4.czi"
>>> ch_names = ['DII-VENUS-N7', 'pPIN1:PIN1-GFP', 'Propidium Iodide', 'pRPS5a:TagBFP-SV40', 'pCLV3:mCherry-N7']
>>> img = imread(im_fname, channel_names=ch_names)
>>> img.is3D()
True
"""
return self.get_channel(self.get_channel_names()[0]).ndim == 3
[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 import MultiChannelImage
>>> from timagetk.io import imread
>>> im_fname = "https://zenodo.org/record/3737795/files/qDII-CLV3-PIN1-PI-E35-LD-SAM4.czi"
>>> ch_names = ['DII-VENUS-N7', 'pPIN1:PIN1-GFP', 'Propidium Iodide', 'pRPS5a:TagBFP-SV40', 'pCLV3:mCherry-N7']
>>> img = imread(im_fname, channel_names=ch_names)
>>> print(img.shape)
(5, 4, 3)
>>> img_t = img.transpose()
>>> # Transpose update the shape attribute of the image (here reversed):
>>> print(img_t.shape)
(5, 4, 3)
>>> # Transpose update the voxelsize attribute of the image (here reversed):
>>> print(img_t.voxelsize)
[0.53, 0.52, 0.51]
>>> # Transpose update the metadata dictionary of the image:
>>> print(img_t.metadata)
>>> # -- Transpose accept axe names as input:
>>> img_t = img.transpose('xyz')
>>> print(img_t.shape)
(5, 4, 3)
>>> img_t = img.transpose('x', 'y', 'z')
>>> print(img_t.shape)
(5, 4, 3)
"""
ch_imgs = {}
for ch in self.channel_names:
ch_imgs[ch] = self[ch].transpose(*axes)
return MultiChannelImage(ch_imgs, name=self.name,
filename=self.filename,
filepath=self.filepath,
metadata=self.metadata.get_dict())
# --------------------------------------------------------------------------
#
# GETTERs & SETTERs:
#
# --------------------------------------------------------------------------
[docs]
def set_name(self, name):
"""Set the name of the MultiChannelImage."""
self.name = name
[docs]
def set_channel_name(self, old, new):
"""Change the name of a channel."""
self[new] = self.pop(old)
self.channel_names = list(self.keys())
[docs]
def get_channel(self, channel_name):
"""Return the intensity image for given channel.
Parameters
----------
channel_name : str
Name of the channel to returns.
Returns
-------
SpatialImage
The intensity image for given channel.
Examples
--------
>>> from timagetk import MultiChannelImage
>>> from timagetk.io import imread
>>> im_fname = "https://zenodo.org/record/3737795/files/qDII-CLV3-PIN1-PI-E35-LD-SAM4.czi"
>>> ch_names = ['DII-VENUS-N7', 'pPIN1:PIN1-GFP', 'Propidium Iodide', 'pRPS5a:TagBFP-SV40', 'pCLV3:mCherry-N7']
>>> im = imread(im_fname, channel_names=ch_names)
>>> # Access the channels as you would with a dictionary:
>>> print(im["CLV3"])
SpatialImage object with following metadata:
- shape: (1024, 1024, 40)
- ndim: 3
- dtype: uint16
- origin: [0, 0, 0]
- voxelsize: [0.208, 0.208, 1.171]
- extent: [212.992, 212.992, 46.84]
>>> # Access the channels using the `get_channel` method:
>>> ahp6_img = im.get_channel("AHP6")
>>> print(ahp6_img)
SpatialImage object with following metadata:
- shape: (1024, 1024, 40)
- ndim: 3
- dtype: uint16
- origin: [0, 0, 0]
- voxelsize: [0.208, 0.208, 1.171]
- extent: [212.992, 212.992, 46.84]
"""
return self[channel_name]
[docs]
def get_channel_names(self):
"""Return the list of channel names.
Returns
-------
list
The list of channel names.
Examples
--------
>>> from timagetk import MultiChannelImage
>>> from timagetk.io import imread
>>> im_fname = "https://zenodo.org/record/3737795/files/qDII-CLV3-PIN1-PI-E35-LD-SAM4.czi"
>>> ch_names = ["qDII", "AHP6", "tagBFP", "CLV3"]
>>> img = imread(im_fname, channel_names=ch_names)
>>> img.get_channel_names()
['qDII', 'AHP6', 'tagBFP', 'CLV3']
"""
return cp.copy(self.channel_names)
[docs]
def get_extent(self, axis=None):
"""Return the get_extent of the given axis.
Parameters
----------
axis : str, in {'x', 'y', 'z'}
Axis for which the get_extent should be returned.
Returns
-------
list of float or float
The get_extent of the axis.
Notes
-----
It is related to the array shape and image voxelsize.
Examples
--------
>>> from timagetk import MultiChannelImage
>>> from timagetk.io import imread
>>> im_fname = "https://zenodo.org/record/3737795/files/qDII-CLV3-PIN1-PI-E35-LD-SAM4.czi"
>>> ch_names = ["qDII", "AHP6", "tagBFP", "CLV3"]
>>> img = imread(im_fname, channel_names=ch_names)
>>> img.get_extent()
[43.337284718258445, 212.3404879282953, 212.3404879282953]
>>> img.get_extent('x')
212.3404879282953
"""
if axis is None:
return cp.copy(self.extent)
else:
return cp.copy(self.extent[self.axes_order_dict[axis.upper()]])
[docs]
def get_voxelsize(self, axis=None):
"""Return the voxelsize of the image or of an axis.
Parameters
----------
axis : {'x', 'y', 'z'}, optional
If specified, axis to use for voxel-size, else the voxel-size of the image.
Returns
-------
list or float
The voxel-size of the image or of the axis (if specified).
Examples
--------
>>> from timagetk import MultiChannelImage
>>> from timagetk.io import imread
>>> im_fname = "https://zenodo.org/record/3737795/files/qDII-CLV3-PIN1-PI-E35-LD-SAM4.czi"
>>> ch_names = ["qDII", "AHP6", "tagBFP", "CLV3"]
>>> img = imread(im_fname, channel_names=ch_names)
>>> img.get_voxelsize()
[0.6771450737227882, 0.2075664593629475, 0.2075664593629475]
>>> img.get_voxelsize('z')
0.6771450737227882
"""
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 get_shape(self, axis=None):
"""Return the shape of the image or of an axis.
Parameters
----------
axis : str in {'x', 'y', 'z', 'c'}, optional
If specified, axis to use for shape dimension, else the shape of the image
Returns
-------
list or int
The shape of the array or an integer if an axis is specified.
Examples
--------
>>> from timagetk import MultiChannelImage
>>> from timagetk.io import imread
>>> im_fname = "https://zenodo.org/record/3737795/files/qDII-CLV3-PIN1-PI-E35-LD-SAM4.czi"
>>> ch_names = ["qDII", "AHP6", "tagBFP", "CLV3"]
>>> img = imread(im_fname, channel_names=ch_names)
>>> img.get_shape()
[65, 1024, 1024, 4]
>>> img.get_shape('z')
65
>>> img.get_shape('c')
4
"""
if axis is None:
return list(cp.copy(self.shape)) + [self.n_ch]
else:
try:
assert axis.lower() in ['x', 'y', 'z', 'c']
except AssertionError:
raise ValueError(
"Parameter `axis` take values in {'x', 'y', 'z', 'c'}!")
except AttributeError:
raise ValueError("Parameter `axis` should be string!")
if axis.lower() == 'c':
return self.n_ch
else:
return cp.copy(self.shape[self.axes_order_dict[axis.upper()]])
[docs]
def get_slice(self, slice_id, axis='z'):
"""Return a MultiChannelImage 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 is 'z'.
Returns
-------
MultiChannelImage
2D MultiChannelImage 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
--------
>>> # Example #1: Use a dummy multichannel image:
>>> from timagetk.array_util import dummy_multichannel_image
>>> img = dummy_multichannel_image()
>>> img_z = img.get_slice(1, 'z')
>>> print(img_z.shape)
(11, 11)
>>> print(img_z.axes_order_dict)
{'Y': 0, 'X': 1, 'C': 2}
>>> # 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, 'C': 2}
>>> # Example #2: Use a real multichannel image (downloaded from Zenodo.org):
>>> from timagetk.io import imread
>>> im_fname = "https://zenodo.org/record/3737795/files/qDII-CLV3-PIN1-PI-E35-LD-SAM4.czi"
>>> ch_names = ["qDII", "AHP6", "tagBFP", "CLV3"]
>>> img = imread(im_fname, channel_names=ch_names)
>>> # Taking an existing z-slice from a 3D image works fine:
>>> img_z = img.get_slice(1, 'z')
>>> img_z.is2D()
True
>>> img_z.is3D()
False
>>> # 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, 'C': 2}
>>> # Down sampling x-axis of a 3D image:
>>> img_ds_x2 = img.get_slice(range(0, img.get_shape('x'), 2), 'x')
>>> img_ds_x2.is3D()
True
>>> print(img.get_shape('x'))
1024
>>> print(img_ds_x2.get_shape('x'))
512
>>> print(img.voxelsize)
[1.1712524251243286e-06, 2.075664593629475e-07, 2.075664593629475e-07]
>>> print(img_ds_x2.voxelsize)
[1.1712524251243286e-06, 2.075664593629475e-07, 4.15132918725895e-07]
>>> # 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')
"""
ch_names = self.get_channel_names()
md = self.metadata.get_dict()
name = self.name
filename = self.filename
filepath = self.filepath
ch_sl = {ch: self[ch].get_slice(slice_id, axis=axis) for ch in ch_names}
return MultiChannelImage(ch_sl, metadata=md, name=name, channel_names=ch_names,
filepath=filepath, filename=filename)
[docs]
def to_blend_image(self, **kwargs):
"""Return an RGB BlendImage object used for plotting.
Other Parameters
----------------
colors : list of str, optional
List of colormap to apply to each image
rtype : str in {'uint8', 'uint16', 'float64'}
Type of array to return, `'uint8'` by default.
Returns
-------
BlendImage
The RGB image resulting in the blending of channels.
Examples
--------
>>> from timagetk.io.image import imread
>>> # We use a freely available multichannel image from an online database
>>> url = "https://zenodo.org/record/3737795/files/qDII-CLV3-PIN1-PI-E35-LD-SAM4.czi"
>>> ch_names = ['DII-VENUS-N7', 'pPIN1:PIN1-GFP', 'Propidium Iodide', 'pRPS5a:TagBFP-SV40', 'pCLV3:mCherry-N7']
>>> ch_colors = ['yellow', 'green', 'red', 'blue', 'purple']
>>> mch_img = imread(url, channel_names=ch_names)
>>> rgb_img = mch_img.to_blend_image(colors=ch_colors, rtype='uint8')
>>> from timagetk.visu.stack import orthogonal_view
>>> orthogonal_view(rgb_img)
"""
return BlendImage(
*[self.get_channel(ch) for ch in self.get_channel_names()],
**kwargs)
[docs]
def random_blend_image(size, n_spimg=2, **kwargs):
"""Generate a random ``BlendImage`` from ``SpatialImage``.
Parameters
----------
size : iterable
Shape of the image to generate.
ZYX sorted if 3D else YX sorted.
n_spimg : int, optional
Number of random SpatialImage to generate to create
Other Parameters
----------------
origin : list
Coordinates of the origin in the image, default: [0, 0] or [0, 0, 0]
voxelsize : list
Image voxelsize, default: [1.0, 1.0] or [1.0, 1.0, 1.0]
dtype : str or numpy.dtype
If specified, default is ``np.float32``, change the dtype of the returned image.
Returns
-------
BlendImage
A random blend image from a list of random spatial images.
See Also
--------
timagetk.components.spatial_image.random_spatial_image
Examples
--------
>>> from timagetk.components.multi_channel import random_blend_image
>>> # Example 1 - Generates a BlendImage using 3D SpatialImages:
>>> blend_img = random_blend_image([5, 10, 10], voxelsize=[1., 0.5, 0.5])
>>> type(blend_img)
timagetk.components.multi_channel.BlendImage
>>> blend_img.shape # the last dimension is of size 3, for RGB
(5, 10, 10, 3)
>>> blend_img.axes_order # the last dimension is RGB
'ZYXB'
>>> blend_img.is3D()
True
>>> # Example 2 - Generates a BlendImage using 2D SpatialImages:
>>> blend_img = random_blend_image([10, 10], voxelsize=[0.5, 0.5])
>>> blend_img.shape # the last dimension is of size 3, for RGB
(10, 10, 3)
>>> blend_img.is2D()
True
"""
from timagetk.array_util import random_spatial_image as rsi
return BlendImage(*[rsi(size, **kwargs) for n in range(n_spimg)])
[docs]
class BlendImage(np.ndarray):
"""RGB image class used to represent blended images like ``MultiChannelImage``.
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].
_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 orders, *e.g.* {'x': 2, 'y': 1, 'z': 0} by default for 3D image.
"""
[docs]
def __new__(cls, *images, **kwargs):
"""Blend image constructor.
YX sorted if images are 2D, ZYX sorted if images are 3D.
Parameters
----------
images : numpy.ndarray or list of SpatialImage
SpatialImages to blend.
If a ``numpy.ndarray``, last dimension should be of size 3 (RGB dim).
Other Parameters
----------------
axes_order : str
The physical axes ordering to use with this array.
Defaults to ``default_axes_order(array)``.
origin : list of int
The coordinate of the origin for this array.
Defaults to ``default_origin(array)``.
unit : float
The units associated to the voxel-sizes to use with this array.
Defaults to ``DEFAULT_SIZE_UNIT``.
voxelsize : list of float
The voxel-sizes to use with this array.
Defaults to ``default_voxelsize(array)``.
metadata : dict
The dictionary of metadata to use with this array.
Defaults to ``{}``.
colors : list of str
List of color names to apply to each image, they can be one of the 140 standard HTML colors or an
hexadecimal value like `'#ff0000'` for 'red'.
rtype : {'uint8', 'uint16', 'float64'}
Type of RGB array to return, set to the type of the first image on the list by default.
See Also
--------
timagetk.visu.util.DEF_CHANNEL_COLORS
timagetk.visu.util.combine_channels
Examples
--------
>>> from timagetk.components.multi_channel import BlendImage
>>> # Example 1 - Initialize with some random 3D RGB array:
>>> from timagetk.array_util import random_array
>>> rand_image = random_array([7, 15, 15, 3], dtype='uint8')
>>> rand_image.shape
(7, 15, 15, 3)
>>> blend_img = BlendImage(rand_image, axes_order='ZYXB')
>>> blend_img.shape
(7, 15, 15, 3)
>>> blend_img.axes_order
'ZYXB'
>>> blend_img.voxelsize # Unspecified, so set to default value
[1.0, 1.0, 1.0]
>>> # Example 2 - Initialize with some random SpatialImages
>>> from timagetk.array_util import random_spatial_image
>>> imgA = random_spatial_image((5, 10, 10), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> imgB = random_spatial_image((5, 10, 10), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> blend_img = BlendImage(imgA, imgB)
>>> blend_img.shape # the last dimension is of size 3, for RGB
(5, 10, 10, 3)
>>> blend_img.dtype
dtype('uint8')
>>> # Example 3 - Initialize from an RGB numpy.ndarray slice
>>> import numpy as np
>>> rgb_arr2d = np.array(blend_img[3, :, :, :])
>>> type(rgb_arr2d)
numpy.ndarray
>>> blend_img2d = BlendImage(rgb_arr2d, voxelsize=blend_img.voxelsize[1:])
>>> blend_img2d.voxelsize
[0.5, 0.5]
"""
from timagetk.components.spatial_image import default_axes_order
from timagetk.components.spatial_image import default_origin
from timagetk.components.spatial_image import default_voxelsize
from timagetk.components.image import assert_spatial_image
from timagetk.io.image import default_image_attributes
from timagetk.util import check_type
from timagetk.util import clean_type
# If the input tuple is of size 1, reduce it to its first element!
if len(images) == 1:
images = images[0]
log.debug(f"Got an `{clean_type(images)}` instance as `images` input!")
# Get attributes from `BlendImage` if not overriden at instance call:
if isinstance(images, BlendImage):
if kwargs.get('axes_order', None) is None:
kwargs['axes_order'] = images.axes_order
if kwargs.get('origin', None) is None:
kwargs['origin'] = images.origin
if kwargs.get('unit', None) is None:
kwargs['unit'] = images.unit
if kwargs.get('voxelsize', None) is None:
kwargs['voxelsize'] = images.voxelsize
# Get the list of colors to use for blending:
colors = kwargs.get('colors', None)
if colors is None:
colors = DEF_CHANNEL_COLORS
if isinstance(images, np.ndarray):
# Need to check physical axes ordering:
axes_order = kwargs.get('axes_order', None)
if axes_order is None:
axes_order = default_axes_order(images) + 'B'
log.warning(f"Undefined parameter ``axes_order``, using default: {axes_order}")
try:
assert 'B' in axes_order
except AssertionError:
axes_order += 'B'
log.warning("Blend axis ('B') is not defined in parameter ``axes_order``.")
log.info(f"Defaulted to `{axes_order}` physical axes ordering!")
log.debug(f"Got an `{axes_order}` sorted RGB array as input for `blendImage`.")
kwargs['axes_order'] = axes_order
# Check array dimensionality and the size oth the blend axe:
cls._check_input_rgb_array(images, axes_order.index('B'))
# Now "rename" the array & update the missing attributes (if any):
array = images
kwargs |= default_image_attributes(array, **kwargs)
else:
# Check we have enough color to blend each image
try:
assert len(images) <= len(colors)
except AssertionError:
raise ValueError("Not enough colors to blend provided images!")
# Check we have SpatialImage instance:
[assert_spatial_image(img) for img in images]
log.debug("Got a list of `SpatialImage` instance as `images` input!")
# Check images shape & voxelsize:
cls._check_image_shape(images)
cls._check_image_voxelsize(images)
cls._check_image_axes_order(images)
# TODO: we could re-order if not all the same
# TODO: Check units!
# Re-order to given 'axes_order' if any::
axes_order = kwargs.get('axes_order', None)
if axes_order is not None and isinstance(axes_order, str):
axo = axes_order.replace('B', '')
if images[0].is3D():
try:
assert len(axo) == 3
except AssertionError:
raise ValueError("Input `axes_order` do not match for a 3D image!")
else:
try:
assert len(axo) == 2
except AssertionError:
raise ValueError("Input `axes_order` do not match for a 2D image!")
images = [img.transpose(axo) for img in images]
# Combining channel:
rtype = kwargs.get('rtype', images[0].dtype.name)
kwargs['axes_order'] = images[0].axes_order + 'B' # Blend/RGB is last axe
kwargs['origin'] = images[0].origin
kwargs['unit'] = images[0].unit
kwargs['voxelsize'] = images[0].voxelsize
array = _combine_channels(images, colors=colors, rtype=rtype)
# - 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__``!
obj = np.asarray(array).view(cls)
# ----------------------------------------------------------------------
# Physical axes ordering:
# ----------------------------------------------------------------------
axes_order = kwargs.get('axes_order', None)
if axes_order is None:
axes_order = default_axes_order(array) + 'B' # Blend/RGB is last axe
log.warning(f"Undefined parameter ``axes_order``, using default: {axes_order}")
else:
log.debug(f"Defined parameter ``axes_order`` as: {axes_order}")
# ----------------------------------------------------------------------
# Origin:
# ----------------------------------------------------------------------
origin = kwargs.get('origin', None)
if origin is None:
origin = default_origin(array)
log.warning(f"Undefined parameter ``origin``, using default: {origin}")
else:
log.debug(f"Defined parameter ``origin`` as: {origin}")
# ----------------------------------------------------------------------
# Size unit:
# ----------------------------------------------------------------------
unit = kwargs.get('unit', None)
if unit is None:
unit = DEFAULT_SIZE_UNIT
log.warning(f"Undefined parameter ``unit``, using default: {unit}")
else:
log.debug(f"Defined parameter ``unit`` as: {unit}")
# ----------------------------------------------------------------------
# Voxel size:
# ----------------------------------------------------------------------
voxelsize = kwargs.get('voxelsize', None)
if voxelsize is None:
voxelsize = default_voxelsize(array)
log.warning(f"No ``voxelsize`` given, using default: {voxelsize}")
else:
voxelsize = list(map(float, voxelsize))
log.debug(f"Defined parameter ``voxelsize`` as: {voxelsize}")
# ----------------------------------------------------------------------
# Metadata:
# ----------------------------------------------------------------------
metadata = kwargs.get('metadata', None)
# - Initialize (empty) or check type of given metadata dictionary:
if metadata is None:
metadata = {}
log.debug("No metadata dictionary provided!")
else:
check_type(metadata, 'metadata', dict)
# - Set 'voxelsize', 'origin', 'unit' & 'extent' hidden attributes:
obj._axes_order = axes_order
obj._voxelsize = voxelsize
obj._origin = origin
obj._extent = compute_extent(voxelsize, array.shape[:-1])
obj._unit = unit
obj._unit_reg = UnitRegistry() # initialize unit registry
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.components.multi_channel import random_blend_image
>>> # Initialize a random (uint8) 3D BlendImage:
>>> blend_img = random_blend_image((5, 10, 10), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> type(blend_img)
timagetk.components.multi_channel.BlendImage
>>> # Taking the first z-slice automatically return a BlendImage instance:
>>> blend_img2d = blend_img[3, :, :, :][0]
>>> type(blend_img2d)
timagetk.components.multi_channel.BlendImage
"""
self._axes_order = getattr(obj, '_axes_order', None)
self._origin = getattr(obj, '_origin', None)
self._voxelsize = getattr(obj, '_voxelsize', None)
self._unit_reg = getattr(obj, '_unit_reg', None)
self._unit = getattr(obj, '_unit', None)
self._extent = getattr(obj, '_extent', None)
@staticmethod
def _check_input_rgb_array(array, blend_axis_idx):
"""Check the requirements for an RGB array as input."""
# - The last dimension should be the RGB axe of size 3
try:
assert array.shape[blend_axis_idx] == 3
except AssertionError:
msg = f"Input RGB array do NOT have a size 3 axis at the specified index ({blend_axis_idx})!"
msg += f"Got an array of 'shape': {array.shape}"
raise ValueError(msg)
# - Test input array dimensionality, should be of dimension 2 or 3:
try:
assert array.ndim in [3, 4]
except AssertionError:
msg = "Input 'RGB array' must be 2D or 3D (excluding the last axis)!"
for attr in ['ndim', 'shape']:
try:
msg += "Got '{}': {}".format(attr, getattr(array, attr))
except:
pass
raise ValueError(msg)
@staticmethod
def _check_image_axes_order(images):
"""Check that all images have the same physical axes ordering.
Parameters
----------
images : list of SpatialImage
A list of SpatialImage.
Raises
------
ValueError
If not all images have the same physical axes ordering.
"""
check_image_axes_order(dict(zip(chnames_generator(len(images)), images)))
@staticmethod
def _check_image_shape(images):
"""Check that all images have the same shape.
Parameters
----------
images : list of SpatialImage
A list of SpatialImage.
Raises
------
ValueError
If not all images are of same shape.
"""
check_image_shape(dict(zip(chnames_generator(len(images)), images)))
@staticmethod
def _check_image_voxelsize(images):
"""Check that all images have the same voxelsize.
Parameters
----------
images : list of SpatialImage
A list of SpatialImage.
Raises
------
ValueError
If not all images are of same voxelsize.
"""
check_image_voxelsize(dict(zip(chnames_generator(len(images)), images)))
# --------------------------------------------------------------------------
#
# BlendImage properties:
#
# --------------------------------------------------------------------------
@property
def axes_order(self):
"""Get ``BlendImage`` physical axes ordering.
Returns
-------
str
``BlendImage`` physical axes ordering.
Example
-------
>>> from timagetk.components.multi_channel import random_blend_image
>>> # Initialize a random (uint8) 3D BlendImage:
>>> blend_img = random_blend_image((5, 10, 10), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> print(blend_img.axes_order)
ZYX
"""
return cp.copy(self._axes_order)
@property
def axes_order_dict(self):
"""Get ``BlendImage`` physical axes ordering.
Returns
-------
dict
``BlendImage`` physical axes ordering.
Example
-------
>>> from timagetk.components.multi_channel import random_blend_image
>>> # Initialize a random (uint8) 3D BlendImage:
>>> blend_img = random_blend_image((5, 10, 10), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> print(blend_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 ``BlendImage`` physical extent.
Returns
-------
list
``BlendImage`` physical extent.
Notes
-----
It is related to the array shape and image voxelsize.
Example
-------
>>> from timagetk.components.multi_channel import random_blend_image
>>> # Initialize a random (uint8) 3D BlendImage:
>>> blend_img = random_blend_image((5, 10, 10), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> print(blend_img.extent)
[3.0, 2.5, 2.5]
"""
return cp.copy(self._extent)
@extent.setter
def extent(self, img_extent):
"""Set ``BlendImage`` physical extent.
Parameters
----------
img_extent : list
``BlendImage`` 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.components.multi_channel import random_blend_image
>>> # Initialize a random (uint8) 3D BlendImage:
>>> blend_img = random_blend_image((5, 10, 10), voxelsize=[1., 0.5, 0.5],dtype='uint8')
>>> print(blend_img.voxelsize)
[1.0, 0.5, 0.5]
>>> blend_img.extent = [10., 10., 10.]
>>> print(blend_img.extent)
[10., 10., 10.]
>>> print(blend_img.voxelsize)
[2.0, 1.0, 1.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 - 1, img_extent)
# - Update 'extent' hidden attribute:
self._extent = _to_list(img_extent)
# - Update 'voxelsize' hidden attribute:
self._voxelsize = compute_voxelsize(self.get_shape(), self.extent)
log.info("Set extent to '{}'".format(self.extent))
log.info("Changed voxelsize to '{}'".format(self.voxelsize))
return
@property
def origin(self):
"""Get ``BlendImage`` origin.
Returns
-------
list
``BlendImage`` origin coordinates
Example
-------
>>> from timagetk.components.multi_channel import random_blend_image
>>> # Initialize a random (uint8) 3D BlendImage:
>>> blend_img = random_blend_image((5, 10, 10), voxelsize=[1., 0.5, 0.5],dtype='uint8')
>>> print(blend_img.origin)
[0, 0, 0]
"""
return cp.copy(self._origin)
@origin.setter
def origin(self, img_origin):
"""Set ``BlendImage`` origin.
Given list should be of same length than the image dimensionality.
Parameters
----------
img_origin : list
``BlendImage`` origin coordinates,
Example
-------
>>> from timagetk.components.multi_channel import random_blend_image
>>> # Initialize a random (uint8) 3D BlendImage:
>>> blend_img = random_blend_image((5, 10, 10), voxelsize=[1., 0.5, 0.5],dtype='uint8')
>>> print(blend_img.origin)
[0, 0, 0]
>>> blend_img.origin = [2, 2, 2]
>>> print(blend_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 - 1, img_origin)
# - Update hidden attribute 'origin':
self._origin = _to_list(img_origin)
log.info("Set origin to '{}'".format(self.origin))
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.components.multi_channel import random_blend_image
>>> # Initialize a random (uint8) 3D BlendImage:
>>> blend_img = random_blend_image((5, 10, 10), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> print(blend_img.get_unit()) # default unit is micrometer
'micrometer'
>>> blend_img.unit = 1e-3 # change to millimeter
>>> print(blend_img.get_unit())
'millimeter'
"""
self._unit = unit
@property
def voxelsize(self):
"""Get ``BlendImage`` voxelsize, sorted as planes, columns & rows.
Returns
-------
list
``BlendImage`` voxelsize
Example
-------
>>> from timagetk.components.multi_channel import random_blend_image
>>> # Initialize a random (uint8) 3D BlendImage:
>>> blend_img = random_blend_image((5, 10, 10), voxelsize=[1., 0.5, 0.5],dtype='uint8')
>>> print(blend_img.voxelsize)
[1., 0.5, 0.5]
"""
return cp.copy(self._voxelsize)
@voxelsize.setter
def voxelsize(self, img_vxs):
"""Change ``BlendImage`` 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.components.multi_channel import random_blend_image
>>> # Initialize a random (uint8) 3D BlendImage:
>>> blend_img = random_blend_image((5, 10, 10), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> print(blend_img.voxelsize)
[1., 0.5, 0.5]
>>> print(blend_img.extent)
[5.0, 5.0, 5.0]
>>> # - Change blend_img voxelsize:
>>> blend_img.voxelsize = [0.5, 0.5, 0.5]
>>> print(blend_img.voxelsize)
[0.5, 0.5, 0.5]
>>> print(blend_img.extent)
[2.5, 5.0, 5.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 - 1, img_vxs)
# - Update 'voxelsize' hidden attribute:
self._voxelsize = _to_list(img_vxs)
# - Update 'extent' hidden attribute:
self._extent = compute_extent(self.voxelsize, self.get_shape())
log.info("Set voxelsize to '{}'".format(self.voxelsize))
log.info("Changed extent to '{}'".format(self.extent))
return
# --------------------------------------------------------------------------
#
# BlendImage methods :
#
# --------------------------------------------------------------------------
[docs]
def is_isometric(self):
"""Test image isometry.
Tests if the voxelsize values are the same for every axis or not.
Returns
-------
is_iso : bool
True if the image is isometric, else ``False``.
Examples
--------
>>> from timagetk.components.multi_channel import random_blend_image
>>> # Initialize a random (uint8) 3D BlendImage:
>>> blend_img = random_blend_image((5, 10, 10), voxelsize=[1., 0.5, 0.5],dtype='uint8')
>>> blend_img.is_isometric()
False
>>> blend_img.voxelsize = [1., 1., 1.]
>>> blend_img.is_isometric()
True
"""
vxs = [self.voxelsize[0]] * (self.ndim - 1)
return np.allclose(vxs, self.voxelsize, atol=1.e-9)
[docs]
def is2D(self):
"""Return True if the BlendImage is 2D, else ``False``.
Examples
--------
>>> from timagetk.components.multi_channel import random_blend_image
>>> # Initialize a random (uint8) 3D BlendImage:
>>> blend_img = random_blend_image((5, 10, 10), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> blend_img.is2D()
False
>>> # Take the first z-slice:
>>> blend_img2d = blend_img.get_slice(0, 'z') # Equivalent to img[0]
>>> blend_img2d.is2D()
True
"""
return self.ndim - 1 == 2
[docs]
def is3D(self):
"""Return True if the BlendImage is 3D, else ``False``.
Examples
--------
>>> from timagetk.components.multi_channel import random_blend_image
>>> # Initialize a random (uint8) 3D BlendImage:
>>> blend_img = random_blend_image((5, 10, 10), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> blend_img.is3D()
True
"""
return self.ndim - 1 == 3
# --------------------------------------------------------------------------
#
# GETTERs & SETTERs:
#
# --------------------------------------------------------------------------
def _get_unit(self):
return (self._unit * self._unit_reg.meter).to_compact()
[docs]
def get_unit(self):
"""Get the human-readable image unit in International System of Units (SI)."""
return str(self._get_unit().u)
[docs]
def get_dtype(self):
"""Get the human-readable image type.
Examples
--------
>>> from timagetk.components.multi_channel import random_blend_image
>>> # Initialize a random (uint8) 3D BlendImage:
>>> blend_img = random_blend_image((5, 10, 10), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> blend_img.get_dtype()
'uint8'
"""
return self.dtype.name
[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.
"""
return np.asarray(self, dtype=dtype).copy()
[docs]
def get_voxelsize(self, axis=None):
"""Return the voxelsize of the image or of an axis.
Parameters
----------
axis : {'x', 'y', 'z'}, optional
If specified, axis to use for voxel-size, else the voxel-size of the image.
Returns
-------
list or float
The voxel-size of the image or of the axis (if specified).
Examples
--------
>>> from timagetk.components.multi_channel import random_blend_image
>>> # Initialize a random (uint8) 3D BlendImage:
>>> blend_img = random_blend_image((5, 10, 10), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> blend_img.get_voxelsize()
[1. , 0.5, 0.5]
>>> blend_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 get_extent(self, axis=None):
"""Return the get_extent of the given axis.
Parameters
----------
axis : str, in {'x', 'y', 'z'}
Axis for which the get_extent should be returned.
Returns
-------
list of float or float
The get_extent of the axis.
Notes
-----
It is related to the array shape and image voxelsize.
Examples
--------
>>> from timagetk.components.multi_channel import random_blend_image
>>> # Initialize a random (uint8) 3D BlendImage:
>>> blend_img = random_blend_image((5, 10, 10), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> blend_img.get_extent()
[4. , 4.5, 4.5]
>>> blend_img.get_extent('x')
4.5
>>> blend_img.get_extent('y')
4.5
>>> blend_img.get_extent('z')
4.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 an axis.
Parameters
----------
axis : str in {'x', 'y', 'z'}, optional
If specified, axis to use for shape dimension, else the shape of the image
Returns
-------
list or int
The shape of the array or an integer if an axis is specified.
Examples
--------
>>> from timagetk.components.multi_channel import random_blend_image
>>> # Initialize a random (uint8) 3D BlendImage:
>>> blend_img = random_blend_image((5, 10, 10), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> blend_img.get_shape()
[3, 5, 5]
>>> blend_img.get_shape('z')
5
>>> blend_img.get_shape('x')
10
"""
if axis is None:
return list(cp.copy(self.shape[:-1]))
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_slice(self, slice_id, axis='z'):
"""Return a ``BlendImage`` with only one slice for given axis.
Parameters
----------
slice_id : int or Iterable
Index of the slice to return.
axis : int or str in {'x', 'y', 'z'}, optional
Axis to use for slicing, default is 'z'.
Returns
-------
BlendImage
2D BlendImage 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.components.multi_channel import random_blend_image
>>> # Initialize a random (uint8) 3D SpatialImage:
>>> blend_img = random_blend_image((5, 10, 10), voxelsize=[1., 0.5, 0.5], dtype='uint8')
>>> blend_img.axes_order_dict
{'Z': 0, 'Y': 1, 'X': 2, 'B': 3}
>>> blend_img.get_shape()
[5, 10, 10]
>>> # Taking an existing z-slice from a 3D image works fine:
>>> blend_img_z = blend_img.get_slice(1, 'z')
>>> blend_img_z.axes_order_dict
{'Y': 0, 'X': 1, 'B': 2}
>>> blend_img_z.shape
(10, 10, 3)
>>> blend_img_z.get_shape()
[10, 10]
>>> # Taking an existing x-slice from a 3D image works fine:
>>> blend_img_x = blend_img.get_slice(3, 'x')
>>> blend_img_x.get_shape()
[5, 10]
>>> # Down-sampling x-axis of a 3D image:
>>> nx = blend_img.get_shape('x')
>>> img_ds_x2 = blend_img.get_slice(range(0, nx, 2), 'x')
>>> img_ds_x2.get_shape()
[5, 10, 5]
>>> # Taking an NON-existing z-slice from a 3D image raises an error:
>>> blend_img.get_slice(50, 'z')
ValueError: Required z-slice (50) does not exists (max: 5)
>>> # Taking a z-slice from a 2D image raises an error:
>>> blend_img_z.get_slice(5, 'z')
"""
if isinstance(axis, str):
axis_id = self.axes_order_dict[axis.upper()]
else:
axis_id = cp.copy(axis)
axis = {v: k for k, v in self.axes_order_dict.items()}[axis]
if axis.lower() == 'z':
# Make sure we have a 3D image:
try:
assert self.is3D()
except AssertionError:
raise ValueError("Can only extract z-slice from 3D images!")
# If a list with just one integer, ``slice_id`` should be an integer:
if isinstance(slice_id, list) and len(slice_id) == 1:
slice_id = slice_id[0]
# Make sure this slice exists:
max_slice = self.get_shape(axis)
if isinstance(slice_id, int):
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))
arr = np.take(self, slice_id, axis=axis_id)
vxs = self.voxelsize
ori = self.origin
axo = self.axes_order
# If ``slice_id`` is an integer, one axis will "disappear", its info should thus be removed:
if isinstance(slice_id, int):
# Remove value corresponding to returned slice axis:
vxs.pop(axis_id)
ori.pop(axis_id)
axo = axo.replace(axis.upper(), '')
else:
# Assume we have a constant compression ratio along the selected axis:
ratio = self.get_shape(axis) / arr.shape[axis_id]
vxs[axis_id] = vxs[axis_id] * ratio
axo = axo
log.debug(f"Slice shape: {arr.shape}")
log.debug(f"Slice ``origin``: {ori}")
log.debug(f"Slice ``voxelsize``: {vxs}")
log.debug(f"Slice ``axes_order``: {axo}")
return BlendImage(arr, origin=ori, voxelsize=vxs, axes_order=axo, unit=self.unit)
[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
-------
BlendImage
The image with permuted axes.
Examples
--------
>>> from timagetk.components.multi_channel import random_blend_image
>>> # Initialize a random (uint8) 3D BlendImage:
>>> blend_img = random_blend_image((5, 6, 7), voxelsize=[1., 0.51, 0.52], dtype='uint8')
>>> blend_img.axes_order
'ZYXB'
>>> blend_img_t = blend_img.transpose()
>>> blend_img_t.axes_order
'XYZB'
>>> # Transpose update the shape attribute of the image (here reversed):
>>> print(img_t.get_shape())
[7, 6, 5]
>>> # Transpose update the voxelsize attribute of the image (here reversed):
>>> print(img_t.get_voxelsize())
[0.52, 0.51, 1.0]
>>> # -- Transpose accept axe names as input:
>>> blend_img_t = blend_img.transpose('xzy') # transpose ZYXB to XZYB
>>> blend_img_t.axes_order
'XZYB'
>>> print(img_t.get_shape())
[7, 6, 5]
>>> img_t = blend_img.transpose('x', 'y', 'z') # transpose ZYXB to XYZB
>>> print(img_t.shape)
[7, 6, 5]
"""
# Handle order change for X, Y & Z dims
if axes == ():
# If no axes order is specified, reverse the axes.
axes = list(range(self.ndim - 1))[::-1]
else:
if len(axes) == 1 and isinstance(axes[0], str) and len(axes[0]) == self.ndim - 1:
# If axes are given as str, translate them to "axe ids":
axes = [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":
axes = [self.axes_order_dict[axe.upper()] for axe in axes]
else:
raise ValueError(f"Could not make sense of given `axes`: {axes}")
# Uses ``axes`` to reorder the 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])
# Transpose the array:
arr = np.transpose(self.get_array(), axes + [max(axes) + 1]) # add the RGB as last dimension
return BlendImage(arr, origin=ori, voxelsize=vxs, axes_order=axo, unit=self.unit)
[docs]
def label_blending(label, image, colors=None, alpha=0.3, bg_label=1, bg_color=(0, 0, 0)):
"""Return an RGB image where color-coded labels are painted over the image.
This is a shorter version of 'skimage.color' module 'label2rgb' function.
Parameters
----------
label : numpy.ndarray, shape (M, N [, P])
Integer array of labels with the same shape as `image`.
image : numpy.ndarray, shape (M, N [, P])
Image used as underlay for labels.
colors : list, optional
List of colors. If the number of labels exceeds the number of colors, then the colors are cycled.
alpha : float [0, 1], optional
Opacity of colorized labels. Ignored if image is ``None``.
bg_label : int, optional
Label that's treated as the background.
bg_color : str or array, optional
Background color. Must be a name in `color_dict` or RGB float values between [0, 1].
Returns
-------
BlendImage
The result of blending a cycling colormap (`colors`) for each distinct value in `label` with the image, at a certain alpha value.
An array of float of shape (M, N, 3)
Notes
-----
Compared to the `label2rgb` function of ``skimage``, we set the `image_alpha` to `1 - alpha`.
Examples
--------
>>> from timagetk import LabelledImage
>>> from timagetk.components.multi_channel import label_blending
>>> from timagetk.io import imread
>>> from timagetk.io.util import shared_dataset
>>> from timagetk.visu.stack import rgb_stack_browser
>>> from timagetk.visu.stack import stack_panel
>>> int_path = shared_dataset("sphere", 'intensity')[0]
>>> seg_path = shared_dataset("sphere", 'segmented')[0]
>>> int_img = imread(int_path)
>>> seg_img = LabelledImage(imread(seg_path), not_a_label=0)
>>> blend = label_blending(seg_img, int_img)
>>> rgb_stack_browser(blend, "Intensity & segmentation blending", val_range=[None, None])
>>> stack_panel(blend, step=10, suptitle="Intensity & segmentation blending", val_range=[None, None])
"""
from skimage.color import gray2rgb
from skimage.color import label2rgb
from timagetk.components.image import assert_spatial_image
assert_spatial_image(label)
assert_spatial_image(image)
if label.ndim == 3:
label_rgb = label2rgb(label, colors=colors, alpha=1, bg_label=bg_label,
bg_color=bg_color)
image_rgb = gray2rgb(img_as_float(image))
return BlendImage(label_rgb * alpha + image_rgb * (1 - alpha))
else:
return BlendImage(label2rgb(label, image, colors=colors, alpha=alpha,
bg_label=bg_label, bg_color=bg_color,
image_alpha=1 - alpha))
def _combine_channels(images, colors=None, rtype='uint8'):
"""Return an RGB image combining the given channels after RGB coloring.
Parameters
----------
images : list of numpy.ndarray or list of SpatialImage
List of grayscale images to blend
colors : list of str, optional
List of color names to apply to each image, they can be one of the 140
standard HTML colors or an hexadecimal value like '#ff0000' for 'red'.
rtype : str in {'uint8', 'uint16', 'float64'}
Type of array to return, 'uint8' by default.
Returns
-------
numpy.ndarray
An RGB array obtained by combining each image with their assigned ``colors``, of shape (M[, N][, P], 3).
"""
from skimage.color import gray2rgb
from timagetk.visu.util import sum_no_overflow
from timagetk.visu.util import color_code
try:
assert rtype in ['uint8', 'uint16', 'float64']
except AssertionError:
raise ValueError(f"Return type is wrong: {rtype}.")
else:
log.debug(f"Selected return type: {rtype}!")
if isinstance(images, tuple):
images = list(images)
if not isinstance(images, list):
raise TypeError(f"Got a {type(images)} for `images` instead of a list!")
n_img = len(images)
if colors is None:
colors = DEF_CHANNEL_COLORS
n_colors = len(colors)
try:
assert n_colors >= n_img
except AssertionError:
msg = f"Not enough colors ({n_colors}) to blend images ({n_img})!"
raise ValueError(msg)
ctype = images[0].get_dtype() # current dtype
# Initialize an empty RGB array (of same type than the first image of the list):
blend = np.zeros_like(gray2rgb(images[0]))
log.debug(f"Initialize an empty blend array of shape {blend.shape}!")
log.debug(f"Initialize an empty blend array of type {blend.dtype}!")
# Loop trough channel to "blend":
# 1. to RGB array
# 2. apply colormap
# 3. add them to blend RGB array
for n_img, img in enumerate(images):
log.debug(f"Adding image #{n_img} with '{colors[n_img]}' colormap...")
c_arr = color_code(colors[n_img])
if ctype == 'uint8':
type_conv = img_as_ubyte
elif ctype == 'uint16':
type_conv = img_as_uint
else:
type_conv = img_as_float
log.debug(f"Converted '{colors[n_img]}' colormap to {c_arr} color array")
blend = sum_no_overflow(blend, type_conv(np.multiply(c_arr, img_as_float(gray2rgb(img)))))
# Convert returned type:
if rtype == 'uint8' and ctype != 'uint8':
blend = img_as_ubyte(blend)
elif rtype == 'uint16' and ctype != 'uint16':
blend = img_as_uint(blend)
elif rtype == 'float64' and ctype != 'float64':
blend = img_as_float(blend)
return blend
[docs]
def combine_channels(images, colors=None, rtype='uint8', **kwargs):
"""Return an RGB image combining the given channels after RGB coloring.
Example of available ``colors`` are: 'red', 'green', 'blue', 'yellow', 'cyan', 'magenta'.
Parameters
----------
images : list(array), shape (M[, N][, P])
List of grayscale images to blend
colors : list of str, optional
List of color names to apply to each image, they can be one of the 140
standard HTML colors or an hexadecimal value like '#ff0000' for 'red'.
rtype : str in {'uint8', 'uint16', 'float64'}
Type of array to return, 'uint8' by default.
Returns
-------
BlendImage
An RGB array obtained by combining each image with their assigned ``colors``, of shape (M[, N][, P], 3).
See Also
--------
timagetk.visu.util.DEF_CHANNEL_COLORS
timagetk.visu.util.color_code
Notes
-----
Beware that channels with co-localized signals will lead to "new color" by
combination, *e.g.* by combinations of red and green, yellow will appear.
Matplotlib method ``imshow`` accept ``float`` or ``uint8`` data types for
RGB images.
In this method we use a ``float64`` array to combine the images.
Since this method is mostly used for image representation, it is converted
back to ``uint8`` in order to save some memory.
Examples
--------
>>> from timagetk.synthetic_data.wall_image import example_layered_sphere_wall_image
>>> from timagetk.synthetic_data.nuclei_image import nuclei_image_from_point_positions
>>> from timagetk.components.multi_channel import combine_channels
>>> from timagetk.visu.mplt import grayscale_imshow
>>> from timagetk.visu.stack import orthogonal_view
>>> # EXAMPLE - Synthetic wall (green) & nuclei (red) image
>>> kwargs = {"extent": 50., "voxelsize": (0.5, 0.25, 0.25), "dtype": 'uint8'}
>>> wall_img, pts = example_layered_sphere_wall_image(return_points=True, wall_intensity=250, **kwargs)
>>> nuclei_img = nuclei_image_from_point_positions(pts, nuclei_intensity=250., **kwargs)
>>> blend = combine_channels([nuclei_img, wall_img])
>>> blend.get_dtype()
'uint8'
>>> # Get the middle z-slice and plot it:
>>> z_sh = blend.get_shape('z')
>>> z_sl = z_sh // 2
>>> grayscale_imshow(blend.get_slice(50), title=f"z-slice{z_sl}/{z_sh}")
>>> # Show an orthogonal view:
>>> orthogonal_view(blend, suptitle="Synthetic wall (green) & nuclei (red) image")
"""
return BlendImage(images, colors=colors, rtype=rtype, **kwargs)
