#!/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.
# ------------------------------------------------------------------------------
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import gridspec
from skimage.exposure import cumulative_distribution
from skimage.exposure import histogram
from tqdm import tqdm
from timagetk import MultiChannelImage
from timagetk.algorithms.connexe import connected_components
from timagetk.algorithms.regionalext import regional_extrema
from timagetk.algorithms.watershed import watershed
from timagetk.bin.logger import get_logger
from timagetk.components.multi_channel import BlendImage
from timagetk.components.multi_channel import combine_channels
from timagetk.components.multi_channel import label_blending
from timagetk.components.spatial_image import SpatialImage
from timagetk.util import type_to_range
from timagetk.visu.util import _force_aspect_ratio
from timagetk.visu.util import convert_str_range
log = get_logger(__name__)
[docs]
def image_plot(image, axe=None, cmap="gray", **kwargs):
"""Display a 2D image with size unit and a coordinate centered pixel.
Parameters
----------
image : timagetk.components.spatial_image.SpatialImage
2D image to represent.
axe : matplotlib.axes.Axes
Use it to combine plots in the same subfigure.
cmap : Any
The Colormap instance or registered colormap name used to map scalar data to colors.
This parameter is ignored for RGB(A) data.
Other Parameters
----------------
val_range : str or list of int
Define the range of values used by the colormap, by default `'type'`.
norm : matplotlib.colors.BoundaryNorm
The `.Normalize` instance used to scale scalar data to the [0, 1] range before mapping to colors using *cmap*.
By default, a linear scaling mapping the lowest value to 0 and the highest to 1 is used.
This parameter is ignored for RGB(A) data.
extent : list of float
If provided, set the extent of the displayed image.
By default, try to use the real extent of the SpatialImage instance or set to ``image.shape`` (voxel unit).
voxelsize : list of float
If provided, set the voxelsize of the displayed image.
By default, try to use the real voxelsize of the SpatialImage instance or set to ``[1., 1.]`` (voxel unit).
unit : str
If provided, set the unit on the image axes.
By default, try to use the real unit of the SpatialImage instance or set to ``'voxels'`` (voxel unit).
interpoplation : str
The interpolation method used. See notes for more details. Defaults to 'none'.
Returns
-------
matplotlib.axes.Axes
The updated `Axes` object.
matplotlib.axes.Axes
The sub-figure object.
Notes
-----
Supported values for `interpolation` are 'none', 'antialiased', 'nearest', 'bilinear', 'bicubic', 'spline16',
'spline36', 'hanning', 'hamming', 'hermite', 'kaiser', 'quadric', 'catrom', 'gaussian', 'bessel', 'mitchell',
'sinc', 'lanczos', 'blackman'.
See ``matplotlib.axes._axes.Axes.imshow`` documentation for more details.
See Also
--------
matplotlib.axes._axes.Axes.imshow
Examples
--------
>>> import matplotlib.pyplot as plt
>>> from timagetk.io import imread
>>> from timagetk.io.util import shared_data
>>> from timagetk.visu.mplt import image_plot
>>> fname = 'p58-t0-a0.lsm'
>>> image_path = shared_data(fname, "p58")
>>> img = imread(image_path)
>>> # - Display the middle z-slice of the 3D image, with 'viridis' colormap:
>>> fig = plt.figure()
>>> ax_img = plt.subplot()
>>> mid_z = img.get_slice(img.get_shape('z')//2, 'z') # get the middle z-slice
>>> ax_img, fig_img = image_plot(mid_z,ax_img,cmap='viridis')
>>> axe.set_title(fname) # add the file name as title
>>> fig.colorbar(fig_img, ax=axe, label='intensity') # add the intensity colorbar
"""
val_range = kwargs.pop('val_range', 'type')
if isinstance(val_range, str):
val_range = convert_str_range(image, val_range)
else:
assert isinstance(val_range, (list, tuple)) and len(val_range) == 2
mini, maxi = val_range
norm = kwargs.pop('norm', None)
if norm is not None:
mini, maxi = None, None # Using `norm` override the `vmin` & `vmax` parameters in `imshow`
try:
y_lab, x_lab = image.axes_order
except ValueError:
y_lab, x_lab = image.axes_order[:2]
if 'extent' in kwargs and len(kwargs['extent']) == 2:
# Override by keyword argument:
extent = kwargs.pop('extent')
else:
try:
extent = image.get_extent() # This is real world units
except:
extent = kwargs.get('extent', image.shape) # By default, assume we are in voxel units
if 'voxelsize' in kwargs and len(kwargs['voxelsize']) == 2:
# Override by keyword argument:
vxs = kwargs.pop('voxelsize')
else:
try:
vxs = image.get_voxelsize() # This is real world units
except:
vxs = kwargs.get('voxelsize', [1., 1.]) # By default, assume we are in voxel units
if 'unit' in kwargs and isinstance(kwargs['unit'], str):
# Override by keyword argument:
unit = kwargs.pop('unit')
else:
try:
unit = image.get_unit() # This is real world units
except:
unit = kwargs.get('unit', 'voxels') # By default, assume we are in voxel units
x_lab += f" ({unit})" # add the size unit to the axis label
y_lab += f" ({unit})" # add the size unit to the axis label
# The pixel coordinate should be at its center:
x_offset = vxs[1] / 2
y_offset = vxs[0] / 2
mpl_extent = (-x_offset, extent[1] + x_offset, extent[0] + y_offset, -y_offset)
if isinstance(image, SpatialImage):
arr = image.get_array() # get an array instance
elif isinstance(image, MultiChannelImage):
arr = combine_channels([image.get_channel(ch) for ch in image.get_channel_names()])
else:
arr = image
interp = kwargs.get('interpolation', 'none')
fig_img = axe.imshow(arr, cmap=cmap, vmin=mini, vmax=maxi, extent=mpl_extent,
interpolation=interp, origin='upper', norm=norm)
axe.xaxis.tick_top() # move the x-axis to the top
if kwargs.get('title', False):
axe.set_title(kwargs['title'])
if kwargs.get('no_xtick', False):
fig_img.axes.xaxis.set_ticklabels([])
if kwargs.get('no_ytick', False):
fig_img.axes.yaxis.set_ticklabels([])
if kwargs.get('no_xlab', False):
axe.set_xlabel("")
else:
axe.set_xlabel(x_lab) # add the axis label
if kwargs.get('no_ylab', False):
axe.set_ylabel("")
else:
axe.set_ylabel(y_lab) # add the axis label
return axe, fig_img
[docs]
def vignette(image, cmap="gray", **kwargs):
"""Create a vignette of a 2D image.
Parameters
----------
image : timagetk.components.spatial_image.SpatialImage
2D image to represent.
cmap : Any
The Colormap instance or registered colormap name used to map scalar data to colors.
This parameter is ignored for RGB(A) data.
Other Parameters
----------------
figname : str
If a string is given, consider it's a valid file location (path, name & extension)
Do not show the figure, just create, write & close.
interpoplation : str
The interpolation method used. See notes for more details. Defaults to 'none'.
"""
fig = plt.figure(figsize=np.array(image.get_extent()) / 10.)
ax = plt.Axes(fig, [0., 0., 1., 1.])
ax.set_axis_off()
fig.add_axes(ax)
ax, fig_img = image_plot(image, ax, cmap=cmap, no_xtick=True, no_ytick=True, no_xlab=True, no_ylab=True, **kwargs)
figname = kwargs.get('figname', "")
# Save the figure if a 'figname' is given:
if figname != "":
plt.savefig(figname)
return fig_img
else:
plt.show()
return None
[docs]
def imshow_property(image, ppty, slice_id=None, axis="z", val_range="auto", cmap='viridis', **kwargs):
"""Plot a labelled image with property values instead of labels.
Parameters
----------
image : timagetk.LabelledImage
Labelled image to map property values.
ppty : dict
Label based dictionary of values.
slice_id : int, optional
Use all slices by default (None), else the z-slice to get for the histogram
axis : str, optional
Axis to use when specifying a `slice_id`, ``'z'`` by default.
val_range : 2-list of int, optional
Define the range of values used by the colormap.
By default, ``"auto"`` scale it to the range of mapped property.
cmap : str, optional
Name of the colormap to use, defaults to ``"viridis"``.
Examples
--------
>>> from timagetk.io import imread
>>> from timagetk import TissueImage3D
>>> from timagetk.io.util import shared_dataset
>>> from timagetk.visu.mplt import imshow_property
>>> from timagetk.visu.mplt import grayscale_imshow
>>> from timagetk.components.labelled_image import labels_at_stack_margins
>>> from timagetk.features.cells import volume
>>> image_path = shared_dataset("p58", "segmented")[0]
>>> image = imread(image_path, TissueImage3D, background=1, not_a_label=0)
>>> # List cells too close to stack margin, as they should be excluded:
>>> margin_labels = labels_at_stack_margins(image, 5)
>>> # Compute the cell volumes:
>>> vol = volume(image, cell_ids=set(image.cell_ids())-set(margin_labels))
>>> axis = 'y'
>>> mid_sh = image.get_shape(axis) // 2
>>> imshow_property(image, vol, mid_sh, axis, cmap='viridis', ppty_name="Volume")
>>> imshow_property(image, vol, cmap='viridis')
"""
thumb_size = kwargs.pop('thumb_size', 8.)
title = kwargs.pop('title', '')
ppty_name = kwargs.pop('ppty_name', '')
if title == "":
title = image.filename
title_suffix = ""
if image.is3D():
if slice_id is not None:
max_slice = image.get_shape(axis)
image = image.get_slice(slice_id, axis=axis)
title_suffix = '{}-slice ({}/{})'.format(axis.upper(), slice_id, max_slice)
else:
title_suffix = "Z-projection"
from timagetk.algorithms.reconstruction import project_segmentation
image = project_segmentation(image, axis=axis, orientation=kwargs.pop('orientation', 1),
background=image.background)
# labels = np.unique(image) # get the list of labels in the 2D image
# assert all([l in ppty for l in labels]) # make sure we have all of them mapped to a value
if val_range == 'auto':
val_range = [min(ppty.values()), max(ppty.values())]
else:
assert isinstance(val_range, (list, tuple)) and len(val_range) == 2
from timagetk.components.labelled_image import relabel_with_property
img = relabel_with_property(image, ppty)
fig, ax = plt.subplots()
ye, xe = img.extent
fig.set_size_inches(w=thumb_size, h=thumb_size * ye / xe)
if ppty_name != "":
fig.suptitle(ppty_name)
ax, fig_img = image_plot(img, ax, cmap=cmap, val_range=val_range)
fig.colorbar(fig_img, ax=ax)
if title_suffix != '':
title += " - " + title_suffix
ax.set_title(title)
figname = kwargs.get('figname', "")
# Save the figure if a 'figname' is given:
if figname != "":
plt.savefig(figname)
return None
elif kwargs.get('no_show', False):
# Using 'no_show' keyword argument to build _multiple_plots, w
return fig_img
else:
plt.show()
return None
[docs]
def grayscale_imshow(image, slice_id=None, axis="z", val_range='type', cmap='gray', **kwargs):
"""Display a 2D grayscale image.
Uses ``matplotlib.imshow()`` and ``_multiple_plots`` if more than one `image` is given.
Parameters
----------
image : timagetk.components.spatial_image.SpatialImage or list of timagetk.components.spatial_image.SpatialImage
2D image(s) to represent, if 3D image(s), a ``contour_projection`` is first computed.
slice_id : int, optional
Use all slices by default (None), else the z-slice to get for the histogram
axis : str, optional
Axis to use when specifying a `slice_id`, ``'z'`` by default.
val_range : str or list of str, optional
Define the range of values used by the colormap, by default *type*
cmap : matplotlib.colors.ListedColormap or str or list of str, optional
Colormap to use, see the ``stack_browser`` notes for advised colormaps
Other Parameters
----------------
title : str
If provided (default is empty), add this string of characters as title
suptitle : str
A general title placed above the sub-figures titles, usually used when a list of images is given
threshold : int
The intensity threshold to use with contour projection when using 3D images.
Setting `slice_id` override the use of this parameter.
orientation : {-1, 1}
The orinetation of the image, use `-1` with an inverted microscope.
In that case the top of the object is at the bottome of the image.
axe : `.axes.Axes` object
Use it to combine plots in the same subfigure
cbar : bool
If ``False``, default is ``True``, no colorbar will be displayed
extent : list of float
If provided (default, ``None``), set the extent of the displayed image.
By default, use the real unit extent of the SpatialImage istance.
max_per_line : int
Number of figure per line when using more than one images, default is 4
thumb_size : float
Image size in inch (default=5.)
no_show : bool
If ``True``, default is ``False``, do not call the blocking 'show' method
figname : str
If a string is given, consider it's a valid file location (path, name & extension)
Do not show the figure, just create, write & close.
interpoplation : str
The interpolation method used. See notes for more details. Defaults to 'none'.
Returns
-------
AxesSubplot
Use it to combine plots in the same subfigure.
Notes
-----
Supported values for `val_range` are:
- **auto**, get the min and max value of the image;
- **type**, get the maximum range from the `image.dtype`, *e.g.* 'uint8'=[0, 255];
- length-2 list of value, *e.g.* [10, 200];
If a list of `image` is given, a list of `title` should also be given.
If a list of `image` is given and only one `val_range` is given, we use the same range for all images.
Supported values for `interpolation` are 'none', 'antialiased', 'nearest', 'bilinear', 'bicubic', 'spline16',
'spline36', 'hanning', 'hamming', 'hermite', 'kaiser', 'quadric', 'catrom', 'gaussian', 'bessel', 'mitchell',
'sinc', 'lanczos', 'blackman'.
See ``matplotlib.axes._axes.Axes.imshow`` documentation for more details.
See Also
--------
matplotlib.axes._axes.Axes.imshow
Examples
--------
>>> from timagetk.io import imread
>>> from timagetk.io.util import shared_dataset
>>> from timagetk.visu.mplt import grayscale_imshow
>>> # - Get 'p58' shared intensity image:
>>> image = imread(shared_dataset("p58")[0])
>>> # Example 1 - Display the middle z-slice of the 3D image stack, with 'viridis' colormap:
>>> mid_z = image.get_shape('z')//2
>>> v = grayscale_imshow(image, mid_z, 'z', cmap='viridis', colorbar=True)
>>> # Example 2 - Display the middle x-slice of two 3D image stacks:
>>> from timagetk.algorithms.exposure import global_contrast_stretch
>>> stretch_img = global_contrast_stretch(image)
>>> mid_x = image.get_shape('x')//2
>>> v = grayscale_imshow([image, stretch_img], mid_x, 'x', suptitle="Effect of global stretching")
>>> # Example 3 - Display a z-axis projection of the 3D image:
>>> v = grayscale_imshow(image, title=image.filename)
"""
def _to_2d(img, slice_id=None):
if img.is2D():
return img, ""
from timagetk.algorithms.reconstruction import image_surface_projection
if slice_id is not None:
max_slice = img.get_shape(axis)
img = img.get_slice(slice_id, axis=axis)
title_suffix = '{}-slice ({}/{})'.format(axis, slice_id, max_slice)
else:
img, altimap = image_surface_projection(img, threshold=kwargs.pop('threshold', None),
orientation=kwargs.pop('orientation', 1))
title_suffix = '{}-axis projection'.format(axis)
return img, title_suffix
thumb_size = kwargs.pop('thumb_size', 5.)
title = kwargs.pop('title', '')
if isinstance(image, list):
# Handles the title for each graph:
if isinstance(title, str):
titles = [title] * len(image)
else:
titles = title
# Make sure we have a 2D image to plot:
for n, img in enumerate(tqdm(image, unit="image")):
image[n], title_suffix = _to_2d(img, slice_id)
if title_suffix != "":
titles[n] = f"{titles[n]} - {title_suffix}" if titles[n] != "" else title_suffix
# Update kwargs to pass them to ``_multiple_plots``
kwargs.update({'val_range': val_range, 'cmap': cmap, 'title': titles})
return _multiple_plots(grayscale_imshow, image, **kwargs)
else:
image, title_suffix = _to_2d(image, slice_id)
if title_suffix != "":
title = f"{title} - {title_suffix}" if title != "" else title_suffix
axe = kwargs.pop('axe', None)
if axe is None:
# Initialise figure:
fig, ax_img = plt.subplots(figsize=[thumb_size, thumb_size])
else:
fig = None
# Use given plt.axe from plt.subplot():
ax_img = axe
norm, bounds = None, None
if isinstance(image, BlendImage):
cmap = None
val_range = [None, None]
else:
from matplotlib import colors
if isinstance(cmap, str):
if cmap == 'glasbey':
from timagetk.visu.util import get_glasbey
cmap = get_glasbey(np.unique(image), not_a_label=image.not_a_label,
background=getattr(image, 'background', 1))
elif isinstance(cmap, colors.LinearSegmentedColormap):
pass
else:
# Assume its a valid matplotlib ListedColormap
mini, maxi = int(image.min()), int(image.max())
values_spread = maxi - mini + 1
try:
assert values_spread == cmap.N
except AssertionError:
raise ValueError(f"Number of colors {cmap.N} and values {values_spread} does not match!")
else:
val_range = [None, None]
bounds = np.linspace(mini - 0.5, maxi + 0.5, cmap.N + 1)
from matplotlib.colors import BoundaryNorm
norm = BoundaryNorm(bounds, cmap.N)
interp = kwargs.get('interpolation', 'none')
ax_img, fig_img = image_plot(image, ax_img, cmap=cmap, val_range=val_range, norm=norm, interpolation=interp)
if kwargs.get('axis_off', False):
plt.axis('off')
if kwargs.get('no_xtick', False):
fig_img.axes.xaxis.set_ticklabels([])
if kwargs.get('no_ytick', False):
fig_img.axes.yaxis.set_ticklabels([])
if kwargs.get('no_xlab', False):
fig_img.axes.xaxis.set_label_text("")
if kwargs.get('no_ylab', False):
fig_img.axes.yaxis.set_label_text("")
ax_img.set_title(title)
if kwargs.get('colorbar', False):
if bounds is not None:
if maxi < 100:
ticks_step = 1
while len(bounds) / ticks_step > 20:
ticks_step += 1
else:
ticks_step = 10
while len(bounds) / ticks_step > 20:
ticks_step += 10
plt.colorbar(fig_img, ax=ax_img, ticks=bounds[::ticks_step][:-1] + 0.5, boundaries=bounds)
else:
plt.colorbar(fig_img, ax=ax_img, orientation="horizontal")
# plt.tight_layout()
subplots_adjust = kwargs.get('subplots_adjust', None)
if subplots_adjust is not None:
plt.subplots_adjust(**subplots_adjust)
figname = kwargs.get('figname', "")
# Save the figure if a 'figname' is given:
if figname != "":
plt.savefig(figname)
elif kwargs.get('no_show', False):
# Using 'no_show' keyword argument to build _multiple_plots, w
return fig_img
else:
plt.show()
return fig
[docs]
def slice_n_hist(image, title="", img_title="", figname="", vmin=None, vmax=None, **kwargs):
"""Display a 2D image with value histogram and cumulative histogram.
Parameters
----------
image : numpy.ndarray or timagetk.components.spatial_image.SpatialImage
2D image to represent
title : str, optional
If provided (default is empty), add this string of characters as title
img_title : str, optional
If provided (default is empty), add this string of characters as title
figname : str, optional
If provided (default is empty), the image will be saved under this filename.
vmin : int or float, optional
Minimum value to use in colormap
vmax : int or float, optional
Maximum value to use in colormap
Other Parameters
----------------
bins: int
The number of bins in the histogram. 256 by default.
val_range : str or list of str, optional
Define the range of values used by the colormap, by default *type*
cmap : str or list of str, optional
Colormap to use, see the ``stack_browser`` notes for advised colormaps
See Also
--------
skimage.exposure.histogram
skimage.exposure.cumulative_distribution
Examples
--------
>>> from timagetk.io import imread
>>> from timagetk.io.util import shared_data
>>> from timagetk.visu.mplt import slice_n_hist
>>> fname = 'p58-t0-a0.lsm'
>>> image_path = shared_data(fname, "p58")
>>> img = imread(image_path)
>>> # - Display the middle z-slice of the 3D image, with 'viridis' colormap:
>>> mid_z = img.get_slice(img.get_shape('z')//2, 'z')
>>> slice_n_hist(mid_z, title=img.filename, cmap='viridis')
"""
try:
assert image.ndim == 2
except:
raise ValueError("Input ``image`` should be 2D!")
mini, maxi = type_to_range(image)
if vmin is None:
vmin = mini
if vmax is None:
vmax = maxi
# Initialise figure:
fig = plt.figure(constrained_layout=True)
gs = fig.add_gridspec(2, 2, width_ratios=[6, 3], height_ratios=[1, 1], hspace=0.2)
plt.suptitle(title)
# Display 2D image:
ax_img = fig.add_subplot(gs[:, 0])
ax_img, fig_img = image_plot(image, ax_img, val_range=[vmin, vmax], **kwargs)
ax_img.set_title(img_title)
# Plot intensity histogram
ax_hist = fig.add_subplot(gs[0, 1])
# plt.hist(image.flatten(), bins=kwargs.get('bins', 256), range=(vmin, vmax + 1), normed=True)
hist, hist_centers = histogram(image.get_array(), nbins=kwargs.get('bins', 256), normalize=True)
ax_hist.plot(hist_centers, hist, lw=2)
ax_hist.set_title('Intensity histogram')
# Plot intensity cumulative histogram
ax_chist = fig.add_subplot(gs[1, 1])
# plt.hist(image.flatten(), bins=kwargs.get('bins', 256), range=(vmin, vmax + 1), cumulative=True,
# histtype='step', normed=True)
chist, chist_centers = cumulative_distribution(image.get_array(), nbins=kwargs.get('bins', 256))
ax_chist.plot(chist_centers, chist, lw=2)
ax_chist.set_title('Cumulative histogram')
if figname != "":
plt.savefig(figname)
return
[docs]
def plot_img_and_hist(image, cmap='gray', suptitle="", title="", **kwargs):
"""Plot an image along with its histogram and cumulative histogram.
Example from .. [#] https://scikit-image.org/docs/stable/auto_examples/color_exposure/plot_local_equalize.html#sphx-glr-auto-examples-color-exposure-plot-local-equalize-py
Parameters
----------
image : timagetk.components.spatial_image.SpatialImage
2D image to represent.
cmap : str or list of str, optional
Colormap to use, see the ``stack_browser`` notes for advised colormaps.
suptitle : str, optional
If provided (default is empty), add this string of characters as figure title.
title : str, optional
If provided (default is empty), add this string of characters as image title.
Other Parameters
----------------
axe : plt.axes
Two matplotlib axes to use to plot the image and the histogram.
extent : list of float
If provided (default, ``None``), set the extent of the displayed image.
By default, use the real unit extent of the SpatialImage istance.
data : numpy.ndarray
An array of date used to create the histograms.
threshold : None or float
A theshold value to display on the histogram.
Notes
-----
As Numpy and Matplotlib have different axis order convention, we transpose the
first two axis of the given numpy array to display the first axis (rows, 'X')
horizontally per the usual plotting convention.
Examples
--------
>>> from timagetk.io import imread
>>> from timagetk.io.util import shared_dataset
>>> from timagetk.visu.mplt import plot_img_and_hist
>>> fname = 'p58-t0-a0.lsm'
>>> image_path = shared_dataset("p58",'intensity')[0]
>>> img = imread(image_path)
>>> # - Display the middle z-slice of the 3D image, with 'viridis' colormap:
>>> mid_z = img.get_slice(img.get_shape('z')//2, 'z')
>>> ax_img, ax_hist, ax_cdf = plot_img_and_hist(mid_z, cmap='viridis')
>>> # - Create a top z-axis projection to display the 3D image and show the 3D image histograms:
>>> from timagetk.visu.projection import projection
>>> proj = projection(img, method='contour', orientation=1)
>>> ax_img, ax_hist, ax_cdf = plot_img_and_hist(proj, data=img)
"""
import matplotlib.ticker as ticker
assert image.ndim == 2
thumb_size = kwargs.pop('thumb_size', 5.)
axe = kwargs.get('axe', None)
if axe is None:
# Initialise figure:
fig, axes = plt.subplots(ncols=2)
fig.set_size_inches(h=thumb_size, w=thumb_size * 2)
ax_img, ax_hist = axes
else:
ax_img, ax_hist = axe
# Share x-axis for histogram and cumulative histogram
ax_cdf = ax_hist.twinx()
# Get minimal and maximal values available depending on image dtype:
vmin, vmax = type_to_range(image)
pv = "Pixels" if image.is2D() else "Voxels"
# Display image
ax_img, fig_img = image_plot(image, ax_img, cmap=cmap, val_range='type', colorbar=True)
ax_img.set_title(title)
fig.colorbar(fig_img, ax=axe, orientation='horizontal', location='bottom', label=f"{pv} intensity")
fig.suptitle(suptitle)
# Plot intensity histogram:
data = kwargs.get('data', image.get_array())
data = data.flatten()
ax_hist.hist(data, bins=min(vmax + 1, 256), range=(vmin, vmax + 1))
ax_hist.ticklabel_format(axis='y', style='scientific', scilimits=(0, 0))
ax_hist.set_xlabel(f"{pv} intensity")
ax_hist.set_ylabel(f"Number of {pv.lower()}")
if image.dtype == 'uint8':
ax_hist.xaxis.set_major_locator(ticker.MultipleLocator(50))
ax_hist.xaxis.set_minor_locator(ticker.MultipleLocator(10))
elif image.dtype == 'uint16':
ax_hist.xaxis.set_major_locator(ticker.MultipleLocator(50 * 255))
ax_hist.xaxis.set_minor_locator(ticker.MultipleLocator(10 * 255))
# Plot intensity cumulative histogram:
ax_cdf.hist(data, bins=min(vmax + 1, 256), range=(vmin, vmax + 1),
cumulative=True, histtype='step', density=True, color='red')
ax_cdf.set_ylabel("Cumulative distribution")
ax_hist.set_xlim(vmin, vmax)
# Add a grid
ax_hist.grid(which='both', linestyle=':')
# Add the threshold vertical line if any:
threshold = kwargs.get('threshold', None)
if threshold is not None:
ax_cdf.vlines([threshold], 0, 1, linestyles='dashdot', color='green')
# Save the figure if a 'figname' is given:
figname = kwargs.get('figname', "")
if figname != "":
plt.savefig(figname)
else:
plt.show()
return ax_img, ax_hist, ax_cdf
[docs]
def image_n_hist(image, cmap='gray', title="", img_title="", figname="",
**kwargs):
"""Plot an image(s) along with its histogram and cumulative histogram.
Parameters
----------
image : numpy.ndarray or timagetk.components.spatial_image.SpatialImage or list
Image to use to generate the histogram, can be a list.
cmap : str, optional
Colormap to use, see the ``stack_browser`` notes for advised colormaps.
title : str, optional
If provided (default is empty), add this string of characters as title.
img_title : str or list of str, optional
If provided (default is empty), add this string of characters as title.
figname : str, optional
If provided (default is empty), the image will be saved under this filename.
Other Parameters
----------------
thumb_size : float
Image size in inch (default=5.)
bins : integer or array_like
If an integer is given, `bins + 1` bin edges are returned.
Unequally spaced bins are supported if ``bins`` is a sequence.
Defaults to ``256``.
Notes
-----
If a list of ``image`` is given, a list of ``img_title`` should also be given.
See Also
--------
timagetk.visu.stack.stack_browser for list of colormaps and some explanations
Examples
--------
>>> from timagetk.io import imread
>>> from timagetk.io.util import shared_data
>>> from timagetk.visu.mplt import image_n_hist
>>> from timagetk.algorithms.exposure import equalize_adapthist
>>> image_path = shared_dataset("p58")[0]
>>> image = imread(image_path)
>>> ks_range = [1/12., 1/10., 1/8., 1/6.]
>>> sh = image.shape[:2]
>>> out_imgs = [equalize_adapthist(image,kernel_size=[ks*sh[0], ks*sh[1]]) for ks in ks_range]
>>> subtitles = ["Original"] + ["Kernel-size {}".format(ks) for ks in ks_range]
>>> image_n_hist([image.get_slice(50, 'z')] + [img.get_slice(50, 'z') for img in out_imgs], title="adaptive histogram equalization (z-slice=50/{})".format(image.shape[-1]), img_title=subtitles)
"""
thumb_size = kwargs.get('thumb_size', 5.)
aspect_ratio = kwargs.get('xy_ratio', None)
# Get grid shape depending on number of input 'image' to display:
if isinstance(image, list):
n_figs = len(image)
# - Check 'img_title'
if isinstance(img_title, list):
assert len(img_title) == n_figs
else:
img_title = [img_title] * n_figs
n_rows = 2
else:
n_figs = 1
n_rows = 2
# Adjust the size of the figure according to the size of the grid
plt.figure(figsize=[thumb_size * n_figs, thumb_size * n_rows])
plt.suptitle(title)
axes = np.zeros((n_rows, n_figs), dtype=object)
if isinstance(image, list):
# - Create the grid with known number of rows and columns:
gs = gridspec.GridSpec(n_rows, n_figs)
# - Add each grayscale image to the grid:
for n, img in enumerate(image):
axes[0, n] = plt.subplot(gs[0, n]) # image (plt.imshow)
axes[1, n] = plt.subplot(gs[1, n]) # histogram (plt.hist)
ax_img, ax_hist, ax_cdf = plot_img_and_hist(img, cmap=cmap, title=img_title[n], axe=axes[:, n],
xy_ratio=aspect_ratio)
else:
axes[0, 0] = plt.subplot(n_rows, n_figs, 1) # image (plt.imshow)
axes[1, 0] = plt.subplot(n_rows, n_figs, 2) # histogram (plt.hist)
ax_img, ax_hist, ax_cdf = plot_img_and_hist(image, cmap=cmap, title=img_title, axe=axes[:, 0],
xy_ratio=aspect_ratio)
# Give more space at the top of the figure for suptitle:
if title != "":
plt.tight_layout(rect=[0.01, 0.01, .99, 0.93])
else:
plt.tight_layout()
# Save the figure if a 'figname' is given:
if figname != "":
plt.savefig(figname)
else:
plt.show()
return
[docs]
def plot_profiles(image, axe, x=None, y=None, z=None, title="", vmin=None, vmax=None):
"""Plot a profile of given line coordinates.
The profile is the line at the intersection of the two planes.
You thus need to specify two , and only two of the three axis.
They act as numpy array index.
Parameters
----------
axe
title
image : timagetk.components.spatial_image.SpatialImage
Image to use to plot the profile.
x : int, optional
X-axis plane coordinates.
y : int, optional
Y-axis plane coordinates.
z : int, optional
Z-axis plane coordinates.
vmin : int | float, optional
Minimum value to use in colormap.
vmax : int | float, optional
Maximum value to use in colormap.
Examples
--------
>>> import matplotlib.pyplot as plt
>>> from timagetk.io import imread
>>> from timagetk.io.util import shared_data
>>> from timagetk.visu.mplt import plot_profiles
>>> image_path = shared_dataset("p58")[0]
>>> image = imread(image_path)
>>> plt.figure(figsize=[9., 5.])
>>> ax = plt.subplot()
>>> # Plot the original profile and the equalized one:
>>> ax = plot_profiles(image, ax, x=int(image.get_shape('x')/2.), z=40, title="original")
>>> from timagetk.algorithms.exposure import equalize_adapthist
>>> eq_img = equalize_adapthist(image)
>>> ax = plot_profiles(eq_img, ax, x=int(image.get_shape('x')/2.), z=40, title="equalized")
>>> from timagetk.algorithms.linearfilter import gaussian_filter
>>> std = 1.
>>> smooth_img = gaussian_filter(image, sigma=std, real=True)
>>> ax = plot_profiles(smooth_img, ax, x=int(image.get_shape('x')/2.), z=40, title="smoothed({}um)".format(std))
>>> eq_smooth_img = equalize_adapthist(gaussian_filter(image, sigma=std, real=True))
>>> ax = plot_profiles(eq_smooth_img, ax, x=int(image.get_shape('x')/2.), z=40, title="equalized & smoothed({}um)".format(std))
>>> plt.legend()
>>> plt.show()
"""
assert sum([axis is None for axis in [x, y, z]]) == 1
from skimage.util.dtype import dtype_range
xmin, xmax = dtype_range[image.dtype.type]
if vmin is None:
vmin = xmin
if vmax is None:
vmax = xmax
# plt.figure(figsize=[9., 5.])
xs, ys, zs = image.shape
if x is None:
profile = image.get_array()[:, y, z]
y_label = "Y: {}/{}; Z:{}/{}".format(y, ys, z, zs)
elif y is None:
profile = image.get_array()[x, :, z]
y_label = "X: {}/{}; Z:{}/{}".format(x, xs, z, zs)
else:
profile = image.get_array()[x, y, :]
y_label = "X: {}/{}; Y:{}/{}".format(x, xs, y, ys)
ax_line = _profiles(axe, profile, title, y_label, vmin, vmax)
return ax_line
def _profiles(axe, profile, legend_label, y_label, vmin, vmax):
"""Profile sub-function.
Parameters
----------
axe : plt.axes
Matplotlib axes to use.
profile : numpy.ndarray
A profile to plot, see ``plot_profile`` for profile definition.
legend_label : str,
Label of the profile, to use with 'plt.legend()'.
y_label : str,
Label to add to y-axis.
vmin : int or float, optional
Minimum value to use in colormap.
vmax : int or float, optional
Maximum value to use in colormap.
Returns
-------
plt.axes
Matplotlib axes
"""
lp = len(profile)
ax_line = axe
ax_line.plot(range(lp), profile, label=legend_label)
ax_line.set_xlim(0, lp)
ax_line.set_ylim(vmin, vmax)
ax_line.set_xlabel(y_label)
ax_line.set_ylabel('Voxel intensity')
return ax_line
[docs]
def profile_details(image, x=None, y=None, z=None, title="", vmin=None,
vmax=None, plane="z", figname="", **kwargs):
"""Plot a profile of given line coordinates and an image zoom-in around it.
The profile is the line at the intersection of the two planes.
You thus need to specify two , and only two of the three axis.
They act as numpy array index.
The 'plane' define which one will be used for the image display.
Parameters
----------
image : timagetk.components.spatial_image.SpatialImage
Image to use to plot the profile.
x : int, optional
X-axis plane coordinates.
y : int, optional
Y-axis plane coordinates.
z : int, optional
Z-axis plane coordinates.
title : str, optional
Title to give to the figure.
vmin : int | float, optional
Minimum value to use in colormap.
vmax : int | float, optional
Maximum value to use in colormap.
plane : str, optional
Plane to use for image display.
figname : str, optional
If provided (default is empty), the image will be saved under this filename.
Other Parameters
----------------
zone : int
Width of the zone around the line.
cmap : str
Colormap to use, see the notes of ``stack_browser`` for advised colormaps.
Examples
--------
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from matplotlib import gridspec
>>> from timagetk.io import imread
>>> from timagetk.io.util import shared_data
>>> from timagetk.visu.mplt import profile_details
>>> image_path = shared_dataset("p58")[0]
>>> image = imread(image_path)
>>> # Plot the profile details:
>>> profile_details(image, x=int(image.get_shape('x')/2.), z=40, title="original")
"""
plt.figure(figsize=[9., 5.])
plt.suptitle(title)
nrow = 2
gs = gridspec.GridSpec(nrow, 1, height_ratios=[1, 8])
gs.update(hspace=0.01)
axes = [plt.subplot(gs[n, 0]) for n in range(nrow)]
axes = _profile_details(image, axes, x, y, z, "", vmin, vmax, plane,
**kwargs)
# Save the figure if a 'figname' is given:
if figname != "":
plt.savefig(figname)
else:
plt.show()
def _profile_details(image, axes, x=None, y=None, z=None, title="", vmin=None,
vmax=None, plane="z", **kwargs):
"""Plot a profile of given line coordinates and an image zoom-in around it.
The profile is the line at the intersection of the two planes.
You thus need to specify two , and only two of the three axis.
They act as numpy array index.
The 'plane' define which one will be used for the image display.
Parameters
----------
image : timagetk.components.spatial_image.SpatialImage
Image to use to plot the profile.
axes : plt.axes
Matplotlib axes to use.
x : int, optional
X-axis plane coordinates.
y : int, optional
Y-axis plane coordinates.
z : int, optional
Z-axis plane coordinates.
vmin : int | float, optional
Minimum value to use in colormap.
vmax : int | float, optional
Maximum value to use in colormap.
plane : str, optional
Plane to use for image display.
Other Parameters
----------------
zone : int
Width of the zone around the line.
cmap : str
Colormap to use, see the notes of ``stack_browser`` for advised colormaps.
Returns
-------
"""
zone = kwargs.get('zone', 5)
cmap = kwargs.get('camp', 'gray')
assert sum([axis is None for axis in [x, y, z]]) == 1
ax_img, ax_line = axes
from skimage.util.dtype import dtype_range
xmin, xmax = dtype_range[image.dtype.type]
if vmin is None:
vmin = xmin
if vmax is None:
vmax = xmax
# plt.figure(figsize=[9., 5.])
xs, ys, zs = image.shape
if x is None:
if plane == 'y':
arr = image.get_array()[:, y, z - zone:z + zone]
else:
arr = image.get_array()[:, y - zone:y + zone, z]
plane = "z"
profile = image.get_array()[:, y, z]
y_label = "Y: {}/{}; Z:{}/{}".format(y, ys, z, zs)
elif y is None:
if plane == 'x':
arr = image.get_array()[x, :, z - zone:z + zone]
else:
arr = image.get_array()[x - zone:x + zone, :, z]
plane = "z"
profile = image.get_array()[x, :, z]
y_label = "X: {}/{}; Z:{}/{}".format(x, xs, z, zs)
else:
if plane == 'x':
arr = image.get_array()[x, y - zone:y + zone, :]
else:
arr = image.get_array()[x - zone:x + zone, y, :]
plane = "y"
profile = image.get_array()[x, y, :]
y_label = "X: {}/{}; Y:{}/{}".format(x, xs, y, ys)
if plane == 'x':
aspect_ratio = image.extent[1] / image.extent[2]
elif plane == 'y':
aspect_ratio = image.extent[0] / image.extent[2]
else:
aspect_ratio = image.extent[0] / image.extent[1]
ax_line = _profiles(ax_line, profile, title, y_label, vmin, vmax)
ax_img.plot([0, arr.shape[1]], [zone - 1, zone - 1], color='red')
ax_img.imshow(arr, cmap=cmap, vmin=vmin, vmax=vmax, aspect=aspect_ratio,
interpolation='none')
ax_img.set_xticks([])
ax_img.set_yticks([])
ax_img.set_ylabel("{}-plane".format(plane.upper()), rotation=0,
verticalalignment='center', horizontalalignment='right')
return ax_img, ax_line
def _get_zone_array(image, x, y, z, plane, zone):
"""Get a 2D sub-SpatialImage from a 3D SpatialImage from given coordinates and plane.
Parameters ``plane`` and ``zone`` indicate the plane and the sampling around the given coordinate of this plane.
Parameters
----------
image : timagetk.components.spatial_image.SpatialImage
Image from which the sub-array should be taken from.
x : int, None
X-axis plane coordinates.
y : int, None
Y-axis plane coordinates.
z : int, None
Z-axis plane coordinates.
plane : str
Plane to use for image display.
zone : int
Width of the zone around the line.
Returns
-------
timagetk.components.spatial_image.SpatialImage
Selected 2D sub-array, with the second axis of dimension ``zone * 2``.
Examples
--------
>>> import numpy as np
>>> from timagetk.io import imread
>>> from timagetk.io.util import shared_data
>>> from timagetk.visu.mplt import _get_zone_array
>>> fname = 'p58-t0-a0.lsm'
>>> image_path = shared_data(fname, "p58")
>>> image = imread(image_path)
>>> mid_x = int(image.get_shape('x')/2.)
>>> mid_y = int(image.get_shape('y')/2.)
>>> mid_z = int(image.get_shape('z')/2.)
>>> arr = _get_zone_array(image, x=mid_x, y=mid_y, z=None, plane='x', zone=5)
>>> print(arr.shape, arr.voxelsize)
>>> arr = _get_zone_array(image, x=mid_x, y=mid_y, z=None, plane='y', zone=5)
>>> print(arr.shape, arr.voxelsize)
>>> arr = _get_zone_array(image, x=mid_x, y=None, z=mid_z, plane='x', zone=5)
>>> print(arr.shape, arr.voxelsize)
>>> arr = _get_zone_array(image, x=mid_x, y=None, z=mid_z, plane='z', zone=5)
>>> print(arr.shape, arr.voxelsize)
>>> arr = _get_zone_array(image, x=None, y=mid_y, z=mid_z, plane='y', zone=5)
>>> print(arr.shape, arr.voxelsize)
>>> arr = _get_zone_array(image, x=None, y=mid_y, z=mid_z, plane='z', zone=5)
>>> print(arr.shape, arr.voxelsize)
"""
vx, vy, vz = image.voxelsize
if x is None:
if plane == 'y':
arr = image.get_array()[:, y, z - zone:z + zone]
arr = SpatialImage(arr, voxelsize=[vx, vz])
else: # z-plane
arr = image.get_array()[:, y - zone:y + zone, z]
arr = SpatialImage(arr, voxelsize=[vx, vy])
elif y is None:
if plane == 'x':
arr = image.get_array()[x, :, z - zone:z + zone]
arr = SpatialImage(arr, voxelsize=[vy, vz])
# arr = SpatialImage(arr, voxelsize=[vy, vz])
else: # z-plane
arr = image.get_array()[x - zone:x + zone, :, z].transpose((1, 0))
arr = SpatialImage(arr, voxelsize=[vy, vx])
else:
if plane == 'x':
arr = image.get_array()[x, y - zone:y + zone, :].transpose((1, 0))
arr = SpatialImage(arr, voxelsize=[vz, vy])
else: # y-plane
arr = image.get_array()[x - zone:x + zone, y, :].transpose((1, 0))
arr = SpatialImage(arr, voxelsize=[vz, vx])
return arr
[docs]
def profile_hmin(image, x=None, y=None, z=None, title="", intensity_range='auto', plane="z", **kwargs):
"""Plot a profile of given line coordinates and an image zoom-in around it.
The profile is the line at the intersection of the two planes.
You thus need to specify two , and only two of the three axis.
They act as numpy array index.
The 'plane' define which one will be used for the image display.
Parameters
----------
image : timagetk.components.spatial_image.SpatialImage
Image to use to plot the profile
x : int, optional
X-axis plane coordinates.
y : int, optional
Y-axis plane coordinates.
z : int, optional
Z-axis plane coordinates.
title : str, optional
Title to give to the figure, *e.g.* the file name, default is empty.
intensity_range : str | list, optional
Define the range of values used by the colormap, by default 'type'.
plane : str, optional
Plane to use for image display.
Other Parameters
----------------
zone : int
Width of the zone around the line.
cmap : str
Colormap to use, see the notes of ``stack_browser`` for advised colormaps.
no_show : bool
If True, do not call `plt.show()` automatically.
hmin_max : int
Change the slider max value of h-minima.
hmin_init : int
Change the slider initial value of h-minima.
hmin_step : int
Change the slider step value of h-minima.
Examples
--------
>>> import numpy as np
>>> from timagetk.io import imread
>>> from timagetk.io.util import shared_data
>>> from timagetk.visu.mplt import profile_hmin
>>> fname = 'p58-t0-a0.lsm'
>>> image_path = shared_data(fname, "p58")
>>> image = imread(image_path)
>>> zslice = 40
>>> x_coord = int(image.get_shape('x')/2.)
>>> # - Interactive 'h-min' search on 'original' image:
>>> profile_hmin(image, x=x_coord, z=zslice, title=fname)
>>> # - Interactive 'h-min' search on 'smoothed' image:
>>> from timagetk.algorithms.linearfilter import gaussian_filter
>>> smooth_img = gaussian_filter(image, sigma=1., real=True)
>>> profile_hmin(smooth_img, x=x_coord, z=zslice, title=fname+" - smoothed")
>>> # - Interactive 'h-min' search on 'isometric & smoothed' image:
>>> from timagetk.algorithms.linearfilter import gaussian_filter
>>> from timagetk.algorithms.resample import isometric_resampling
>>> iso_smooth_img = gaussian_filter(isometric_resampling(image), sigma=0.3, real=True)
>>> iso_zslice = iso_smooth_img.get_shape('z') / image.get_shape('z') * zslice
>>> profile_hmin(iso_smooth_img, x=x_coord, z=iso_zslice, title=fname+" - isometric & smoothed")
>>> # - Interactive 'h-min' search on 'isometric, smoothed & equalized' image:
>>> from timagetk.algorithms.exposure import equalize_adapthist
>>> iso_smooth_eq_img = equalize_adapthist(iso_smooth_img)
>>> profile_hmin(iso_smooth_eq_img, x=x_coord, z=iso_zslice, title=fname+" - isometric, smoothed & equalized'")
"""
zone = kwargs.get('zone', 10)
cmap = kwargs.get('cmap', 'gray')
if image.dtype == np.uint8:
hmin_init = 5
hmin_max = 45
hmin_step = 1
else:
hmin_init = 500
hmin_max = 4500
hmin_step = 100
hmin_max = kwargs.get('hmin_max', hmin_max)
hmin_init = kwargs.get('hmin_init', hmin_init)
hmin_step = kwargs.get('hmin_step', hmin_step)
xs, ys, zs = image.shape
if x is None:
if plane != 'y':
plane = "z"
profile = image.get_array()[:, y, z]
y_label = "Y: {}/{}; Z:{}/{}".format(y, ys, z, zs)
elif y is None:
if plane != 'x':
plane = "z"
profile = image.get_array()[x, :, z]
y_label = "X: {}/{}; Z:{}/{}".format(x, xs, z, zs)
else:
if plane != 'x':
plane = "y"
profile = image.get_array()[x, y, :]
y_label = "X: {}/{}; Y:{}/{}".format(x, xs, y, ys)
# - Compute segmented sub-image:
subimg = _get_zone_array(image, x, y, z, plane, zone)
aspect_ratio = subimg.extent[0] / subimg.extent[1]
if isinstance(intensity_range, str):
intensity_range = convert_str_range(subimg, intensity_range)
else:
assert isinstance(intensity_range, list) and len(intensity_range) == 2
mini, maxi = intensity_range
# - h-transform:
ext_img = regional_extrema(subimg, height=hmin_init, method='min', connectivity=8)
if ext_img is None:
ext_img = np.zeros_like(subimg.get_array())
# - connexe component labelling:
seed_img = connected_components(ext_img, low_threshold=1, high_threshold=hmin_init, connectivity=8)
n_seeds = len(np.unique(seed_img)) - 1 # '0' is in the list!
log.info("Detected {} seeds!".format(n_seeds))
if seed_img is None:
seed_img = np.zeros_like(subimg.get_array())
# - watershed segmentation:
seg_img = watershed(subimg, seed_img)
if seg_img is None:
seg_img = np.zeros_like(subimg)
# - blend intensity and labelled image (return a numpy array):
blend = label_blending(seg_img, subimg)
# -- Initialise figure
fig = plt.figure(figsize=[12., 6.])
plt.subplots_adjust(bottom=0.1)
plt.suptitle(title)
nrow = 5
gs = gridspec.GridSpec(nrow, 1)
gs.update(hspace=0.01)
ax_img, ax_ext, ax_seg, ax_blend, ax_line = [plt.subplot(gs[n, 0]) for n in range(nrow)]
# - Plot intensity sub-image:
ax_img.plot([0, subimg.get_shape('x')], [zone - 1, zone - 1], color='red')
ax_img.imshow(subimg.get_array().transpose((1, 0)), cmap=cmap, vmin=mini, vmax=maxi, interpolation='none')
ax_img.set_ylabel("{}-plane".format(plane.upper()), rotation=0, verticalalignment='center',
horizontalalignment='right')
# - Plot detected regional-minima in sub-image:
ax_ext.plot([0, subimg.get_shape('x')], [zone - 1, zone - 1], color='black')
ext_fig = ax_ext.imshow(ext_img.get_array().transpose((1, 0)), cmap='viridis', interpolation='none', vmin=0,
vmax=np.max(ext_img))
ax_ext.set_ylabel("h-minima", rotation=0, verticalalignment='center', horizontalalignment='right')
# - Plot segmented sub-image:
ax_seg.plot([0, subimg.get_shape('x')], [zone - 1, zone - 1], color='black')
seg_fig = ax_seg.imshow(seg_img.get_array().transpose((1, 0)), cmap='Set1', interpolation='none', vmin=0,
vmax=np.max(seed_img))
ax_seg.set_ylabel("segmentation", rotation=0, verticalalignment='center', horizontalalignment='right')
# - Plot blending of intensity & segmented sub-image:
ax_blend.plot([0, subimg.get_shape('x')], [zone - 1, zone - 1], color='black')
blend_fig = ax_blend.imshow(blend.transpose((1, 0, 2)), interpolation='none')
# blend_fig = ax_blend.imshow(seg_img.get_array().transpose((1, 0)), cmap='Set1', aspect=aspect_ratio,
# interpolation='none')
ax_blend.set_ylabel("blending", rotation=0, verticalalignment='center', horizontalalignment='right')
# Force correct aspect ratio of image and remove x & y ticks:
for ax in [ax_img, ax_ext, ax_seg, ax_blend]:
_force_aspect_ratio(ax, aspect_ratio)
ax.set_xticks([])
ax.set_yticks([])
# - Plot profile:
ax_line = _profiles(ax_line, profile, "profile", y_label, mini, maxi)
# - Plot the h-extrema:
hext_fig, = ax_line.plot(range(len(profile)), ext_img[:, zone - 1], label="h-extrema")
from timagetk.visu.stack import slider
hs = slider(label='h-minima', mini=1, maxi=hmin_max, init=hmin_init, step=hmin_step, fmt="%1.0f")
def _update(val):
hmin = hs.val
# - h-transform:
ext_img = regional_extrema(subimg, height=hmin, method='min')
if ext_img is None:
ext_img = np.zeros_like(subimg)
# - connexe component labelling:
seed_img = connected_components(ext_img, low_threshold=1, high_threshold=hmin)
n_seeds = len(np.unique(seed_img)) - 1 # '0' is in the list!
log.info("Detected {} seeds!".format(n_seeds))
if seed_img is None:
seed_img = np.zeros_like(subimg)
# - watershed segmentation:
seg_img = watershed(subimg, seed_img)
if seg_img is None:
seg_img = np.zeros_like(subimg)
# - blend intensity and labelled image:
blend = label_blending(seg_img, subimg)
# - Update figures:
ext_fig.set_data(ext_img.get_array().transpose((1, 0)))
hext_fig.set_ydata(ext_img[:, zone - 1])
seg_fig.set_data(seg_img.get_array().transpose((1, 0)))
seg_fig.set_clim([np.min(seg_img.get_array()), np.max(seg_img.get_array())])
blend_fig.set_data(blend.transpose((1, 0, 2)))
fig.canvas.draw_idle()
hs.on_changed(_update)
ax_line.legend(fontsize='small', framealpha=0.2, frameon=True)
if kwargs.get('no_show', False):
return None
else:
plt.show()
return int(hs.val)
# def segmented_profiles(image, axe=None, x=None, y=None, z=None, hmin=20,
# title="", vmin=None,
# vmax=None):
# """Plot a profile of given line coordinates.
#
# The profile is the line at the intersection of the two planes.
# You thus need to specify two , and only two of the three axis.
# They act as numpy array index.
#
# Parameters
# ----------
# image : timagetk.components.spatial_image.SpatialImage
# Image to use to plot the profile
# x : int, optional
# x-axis plane coordinates
# y : int, optional
# x-axis plane coordinates
# z : int, optional
# x-axis plane coordinates
# vmin : int | float, optional
# minimum value to use in colormap
# vmax : int | float, optional
# maximum value to use in colormap
#
# Examples
# --------
# >>> import matplotlib.pyplot as plt
# >>> from timagetk.io import imread
# >>> from timagetk.io.util import shared_data
# >>> from timagetk.visu.mplt import plot_profiles
# >>> image_path = shared_dataset("p58")[0]
# >>> image = imread(image_path)
#
# >>> plt.figure(figsize=[9., 5.])
# >>> # Plot the original profile and the equalized one:
# >>> ax = plot_profiles(image, x=int(image.get_shape('x')/2.), z=40, title="original")
# >>> from timagetk.algorithms.exposure import z_slice_equalize_adapthist
# >>> eq_img = z_slice_equalize_adapthist(image)
# >>> ax = plot_profiles(eq_img, ax, x=int(image.get_shape('x')/2.), z=40, title="equalized")
# >>> from timagetk.plugins import linear_filtering
# >>> std = 1.
# >>> smooth_img = linear_filtering(image, method="gaussian_smoothing", sigma=std, real=True)
# >>> ax = plot_profiles(smooth_img, ax, x=int(image.get_shape('x')/2.), z=40, title="smoothed({}um)".format(std))
# >>> eq_smooth_img = z_slice_equalize_adapthist(linear_filtering(image, method="gaussian_smoothing", sigma=std, real=True))
# >>> ax = plot_profiles(eq_smooth_img, ax, x=int(image.get_shape('x')/2.), z=40, title="equalized & smoothed({}um)".format(std))
# >>> plt.legend()
# >>> plt.show()
#
# Returns
# -------
#
# """
# assert sum([axis is None for axis in [x, y, z]]) == 1
# from skimage.util.dtype import dtype_range
#
# xmin, xmax = dtype_range[image.dtype.type]
# if vmin is None:
# vmin = xmin
# if vmax is None:
# vmax = xmax
#
# # plt.figure(figsize=[9., 5.])
# xs, ys, zs = image.shape
# if x is None:
# profile = image.get_array()[:, y, z]
# y_label = "Y: {}/{}; Z:{}/{}".format(y, ys, z, zs)
# elif y is None:
# profile = image.get_array()[x, :, z]
# y_label = "X: {}/{}; Z:{}/{}".format(x, xs, z, zs)
# else:
# profile = image.get_array()[x, y, :]
# y_label = "X: {}/{}; Y:{}/{}".format(x, xs, y, ys)
#
# if axe is None:
# axe = plt.subplot()
#
# x = range(len(profile))
# y = profile
# seg_profile = seg_img[zone - 1, :]
#
# points = np.array([x, y]).T.reshape(-1, 1, 2)
# segments = np.concatenate([points[:-1], points[1:]], axis=1)
#
# norm = plt.Normalize(seg_profile.min(), seg_profile.max())
# lc = LineCollection(segments, cmap='Set1', norm=norm)
# lc.set_array(seg_profile)
# lc.set_linewidth(2)
# lc.set_label('segmented profile')
# segline = ax_blendline.add_collection(lc)
# print("here {}".format(n))
#
# n += 1
#
# return ax_line
[docs]
def intensity_histogram(image, slice_id=None, axis="z", **kwargs):
"""Intensity distribution as histogram and cumulative histogram.
Uses ``matplotlib.hist()`` and ``_multiple_plots`` if more than one `image` is given.
Parameters
----------
image : numpy.ndarray or timagetk.components.spatial_image.SpatialImage
Image to use to generate the histogram, can also be a list.
slice_id : int, optional
Use all slices by default (None), else the z-slice to get for the histogram.
axis : str, optional
Axis to use when specifying a `slice_id`, ``'z'`` by default.
Other Parameters
----------------
title : list of str
A list of titles to use as sub-figures titles, usually used when a list of images is given.
suptitle : str
A general title placed above the sub-figures titles, usually used when a list of images is given.
axe : plt.axes
Matplotlib axes to use with ``_multiple_plots``, then ``no_show`` should be ``True``.
max_per_line : int
Number of figure per line when using more than one images, default is ``4``.
thumb_size : float
Image size in inch, default is ``5``.
no_show : bool
If ``True``do not call the blocking 'show' method.
Default is ``False``.
no_left_ylab : bool
If ``True``remove the left y-axis label.
Default is ``False``.
no_right_ylab : bool
If ``True``remove the right y-axis label.
Default is ``False``.
Examples
--------
>>> from timagetk.io.util import shared_data
>>> from timagetk.io import imread
>>> from timagetk.visu.mplt import intensity_histogram
>>> import matplotlib.pyplot as plt
>>> img = imread(shared_data("p58-t0-a0.lsm", "p58"))
>>> intensity_histogram(img)
>>> intensity_histogram(img, slice_id=5, axis="z")
>>> from timagetk.algorithms.exposure import global_contrast_stretch
>>> stretch_img = global_contrast_stretch(img)
>>> intensity_histogram([img, stretch_img], suptitle="Effect of global stretching on intensity distribution")
"""
if isinstance(image, list):
return _multiple_plots(intensity_histogram, image, cbar=False, **kwargs)
thumb_size = kwargs.pop('thumb_size', 5.)
axe = kwargs.get('axe', None)
if axe is None:
# Initialise figure:
fig, ax_hist = plt.subplots()
fig.set_size_inches(w=thumb_size, h=thumb_size)
else:
# Use given plt.axe from plt.subplot():
ax_hist = axe
# Share x-axis for histogram and cumulative histogram
ax_cdf = ax_hist.twinx()
# Get minimal and maximal values available depending on image dtype:
vmin, vmax = type_to_range(image)
if slice_id is not None:
max_slice = image.get_shape(axis)
image = image.get_slice(slice_id, axis=axis)
title_suffix = '{}-slice ({}/{})'.format(axis, slice_id, max_slice)
else:
title_suffix = ""
pv = "Pixels" if image.is2D() else "Voxels"
# Plot intensity histogram:
hist_fig = ax_hist.hist(image.flatten(), bins=kwargs.get('bins', 256), range=(vmin, vmax + 1))
ax_hist.ticklabel_format(axis='y', style='scientific', scilimits=(0, 0))
ax_hist.set_xlabel("{} intensity".format(pv))
if not kwargs.get('no_left_ylab', False):
ax_hist.set_ylabel("Number of {}".format(pv.lower()))
# Plot intensity cumulative histogram:
ax_cdf.hist(image.flatten(), bins=kwargs.get('bins', 256), range=(vmin, vmax + 1), cumulative=True, histtype='step',
density=True,
color='red')
if not kwargs.get('no_right_ylab', False):
ax_cdf.set_ylabel("Cumulative distribution")
ax_hist.set_xlim(vmin, vmax)
title = kwargs.get('title', '')
if title == '':
try:
title = image.filename
except AttributeError:
pass
if title == '':
title += title_suffix
elif title_suffix != '':
title += " - " + title_suffix
ax_hist.set_title(title)
if title != "":
plt.subplots_adjust(top=0.95)
else:
plt.subplots_adjust(top=0.99)
figname = kwargs.get('figname', "")
# Save the figure if a 'figname' is given:
if figname != "":
plt.savefig(figname)
return None
elif kwargs.get('no_show', False):
return hist_fig, 1.
else:
plt.show()
return None
def _multiple_plots(plotting_func, image, **kwargs):
"""Iterate `plotting_func` to match number of given `image`.
Parameters
----------
plotting_func : func
Plotting function to iterate.
image : list of timagetk.components.spatial_image.SpatialImage
List of image to which to apply the plotting function.
Other Parameters
----------------
title : str, optional
If provided (default is empty), add this string of characters as title.
suptitle : str
A general title placed above the sub-figures titles, usually used when a list of images is given.
max_per_line : int
Maximum number of figure per line (default=4).
thumb_size : float
Image size in inch (default=5.).
Returns
-------
matplotlib.pyplot.figure.Figure
A matplotlib figure
"""
assert isinstance(image, list)
n_figs = len(image)
suptitle = kwargs.pop('suptitle', '')
figname = kwargs.pop('figname', '')
no_show = kwargs.pop('no_show', False)
thumb_size = kwargs.pop('thumb_size', 5.)
max_per_line = int(kwargs.pop('max_per_line', 4))
title = kwargs.pop('title', None)
if title:
assert len(title) == n_figs
val_range = kwargs.pop('val_range', None)
if isinstance(val_range, list) and not isinstance(val_range[0], int):
assert len(val_range) == n_figs
elif val_range:
val_range = [val_range] * n_figs
# Control 'threshold', used by `grayscale_imshow`,
cmap = kwargs.pop('cmap', None)
if isinstance(cmap, list):
assert len(cmap) == n_figs
elif cmap:
cmap = [cmap] * n_figs
# Control 'threshold', used by `grayscale_imshow` for `contour_projection`
threshold = kwargs.pop('threshold', None)
if isinstance(threshold, list):
assert len(threshold) == n_figs
elif threshold:
threshold = [threshold] * n_figs
# Control 'colorbar', used by `grayscale_imshow`, if defined manually, no auto-control
if 'cbar' not in kwargs:
all_same_range = np.alltrue([vrange == val_range[0] for vrange in val_range])
all_same_cmap = np.alltrue([cm == cmap[0] for cm in cmap])
if all_same_range and all_same_cmap:
# Colorbars are only displayed for the last figure if value range is always the same:
kwargs.update({'cbar': False})
else:
# If auto-scaling or specified list of value range for colorbar, need to always have it
kwargs.update({'cbar': True})
else:
all_same_range = False
all_same_cmap = False
# - Get the number of figure to display per rows, *i.e.* the number of columns:
n_columns = min([n_figs, max_per_line])
# - Get the number of rows:
n_rows = n_figs // max_per_line
n_rows += 1 if n_rows != (n_figs / max_per_line) else 0
ref_extent = image[0].get_extent()
same_extent = all(im.get_extent() == ref_extent for im in image)
try:
ref_height, ref_width = image[0].get_extent()
except ValueError:
ref_height, ref_width, _ = image[0].get_extent()
same_width = all(im.get_extent()[1] == ref_width for im in image)
if not same_width:
max_width = max([im.get_extent()[1] for im in image])
width_ratios = [thumb_size * im.get_extent()[1] / max_width for im in image]
else:
width_ratios = [thumb_size] * len(image)
hw_ratio = [np.divide(*im.get_extent()) for im in image]
max_hw_ratio = max(hw_ratio)
if max_hw_ratio > 1.:
thumb_height = max([hwr * thumb_size for hwr in hw_ratio])
else:
thumb_height = thumb_size
fig, axes = plt.subplots(n_rows, n_columns, gridspec_kw={'width_ratios': width_ratios})
if n_rows == 1 or n_columns == 1:
axes = axes.reshape(n_rows, n_columns)
plt.suptitle(suptitle)
id_col, id_row = 0, 0
sub_figs = []
for n in range(n_figs):
img = image[n]
if title is not None:
kwargs['title'] = title[n]
if cmap is not None:
kwargs['cmap'] = cmap[n]
if threshold is not None:
kwargs['threshold'] = threshold[n]
if val_range:
kwargs['val_range'] = val_range[n]
# Left y-label are displayed only for the first plot of each row:
first_col = id_col == 0
kwargs['no_ylab'] = True if not first_col else False
# Bottom x-label are displayed only for the last row:
first_row = id_row == 0
last_row = id_row == n_rows - 1
kwargs['no_xlab'] = True if not last_row else False
if same_extent:
# If all images have the same extent we can remove the Y-ticks for images after the first one (per row)
kwargs['no_ytick'] = True if not first_col else False
# If all images have the same extent we can remove the X-ticks for images after the first one (per columns)
kwargs['no_xtick'] = True if not first_row else False
# Call plotting function:
s_fig = plotting_func(img, axe=axes[id_row, id_col], no_show=True, **kwargs)
sub_figs.append(s_fig)
if id_col == max_per_line - 1:
id_col = 0
id_row += 1
else:
id_col += 1
h_offset = 0.5
if suptitle != "":
h_offset += 0.5 # save some inches to display suptitle
if all_same_range and all_same_cmap:
plt.colorbar(sub_figs[0], ax=axes, orientation="horizontal", fraction=.1, aspect=40)
h_offset += 0.8 # save some inches to display colorbar
fig.set_size_inches(w=thumb_size * n_columns, h=(thumb_height * n_rows) + h_offset)
# plt.tight_layout()
# Give more space at the top of the figure for suptitle:
if suptitle != "" or title:
plt.subplots_adjust(top=0.95)
else:
plt.subplots_adjust(top=0.99)
if plotting_func.__name__ == 'grayscale_imshow':
if all_same_range and all_same_cmap:
# Case where we put a colorbar at the bottom below all subplots:
plt.subplots_adjust(bottom=0.2)
plt.subplots_adjust(hspace=0.2)
else:
plt.subplots_adjust(bottom=0.1)
plt.subplots_adjust(hspace=0.2)
plt.subplots_adjust(wspace=0.09)
# Save the figure if a 'figname' is given:
if figname != "":
plt.savefig(figname)
plt.close()
elif no_show:
pass
else:
plt.show()
return fig, axes
