Reconstruction

Surface reconstruction & intensity projections module.

Methods and algorithms to reconstruct the outer shape of an object in the image. Can preform projections for intensity and labelled images.

timagetk.algorithms.reconstruction.altimap_from_segmentation(seg_image, axis='z', background_id=1, orientation=None, **kwargs)[source]

Get an altitude map from a segmentation, using the background position.

Height values follows the slices indexing, so 0 for the first slice (hence the importance of orientation). They range in [0, axis_shape - 1] with axis_shape being the size of the Z projected axis. They are thus in voxel unit and not in real units.

Parameters:

  • seg_image (timagetk.LabelledImage or timagetk.TissueImage3D) – The segmented image to use to detect the background position.

  • axis ({'x', 'y', 'z'}, optional) – The axis to use for projection, default to z.

  • background_id (int, optional) – The label associated to the background, defaults to 1. If a TissueImage3D instance is given as seg_image, it will use its background attribute.

  • orientation ({None, -1, 1}, optional) – Defines the orientation of the projection. By default (None), try to guess it with guess_image_orientation. If -1 we consider the z-axis starts at the bottom of the object and goes upward. If 1 we consider the z-axis starts at the top of the object and goes downward.

  • iso_upsampling (bool) – If True, up-sample the image to the minimum voxelsize. This is usefull if you have an axis with a lower resolution, typically the z-axis.

  • interpolation ({'nearest', 'cellbased'}) – Interpolation method to use during iso_upsampling, if any. Defaults to ‘cellbased’.

  • cell_based_sigma (float) – Required when using “cellbased” interpolation method, sigma for cell-based interpolation. In voxel units!

  • altimap_smoothing (bool) – If True, default is False, smooth the altitude map obtained from the mask.

  • gaussian_sigma (float, optional) – Sigma value, in real units, of the Gaussian filter applied to the altutide image. This value is divided by the image voxelsize.

  • connectivity ({1, 2, 3}) – The connectivity of the binary structuring element. Defaults to 1.

  • iterations (int) – The number of iterations of binary dilation by the structuring element. Defaults to 1.

Returns:

timagetk.SpatialImage – YX 2D altitude map, height values starting from the first slice defined by the image orientation.

See also

timagetk.algorithms.reconstruction.contour_from_segmentation

Examples

>>> import numpy as np
>>> from timagetk.algorithms.reconstruction import altimap_from_segmentation
>>> from timagetk import LabelledImage
>>> from timagetk.io import imread
>>> from timagetk.io.util import shared_dataset
>>> from timagetk.visu.mplt import grayscale_imshow
>>> # - Get 'p58-t0' shared segmented image:
>>> seg_image = imread(shared_dataset("p58", "segmented")[0], LabelledImage)
>>> alti = altimap_from_segmentation(seg_image, altimap_smoothing=True)
>>> grayscale_imshow(alti, title="Altitude map / Z-axis", cmap='viridis', val_range='auto')
>>> alti = altimap_from_segmentation(seg_image, axis='x', altimap_smoothing=False)
>>> grayscale_imshow(alti, title="Altitude map / X-axis", cmap='viridis', val_range='auto')
>>> alti = altimap_from_segmentation(seg_image, axis='y', altimap_smoothing=False)
>>> grayscale_imshow(alti, title="Altitude map / Y-axis", cmap='viridis', val_range='auto')
timagetk.algorithms.reconstruction.altimap_surface_projection(image, altimap, axis='z', offset=1, **kwargs)[source]

Return the surface from the altitude map with image values.

Parameters:

  • image (timagetk.SpatialImage or timagetk.LabelledImage) – The image to extract values from (at given altitude).

  • altimap (timagetk.SpatialImage) – The altitude map at witch image values should be taken from.

  • axis ({'x', 'y', 'z'}, optional) – The axis to use for projection, default to z. Should be the same as for the computation of the altimap.

  • offset (int, optional) – The offset to apply to the altimap, defaults to 1.

Returns:

numpy.ndarray – The 2D image at given altitude.

Examples

>>> import numpy as np
>>> from timagetk.algorithms.reconstruction import altimap_surface_projection
>>> from timagetk.algorithms.reconstruction import altimap_from_segmentation
>>> from timagetk import LabelledImage
>>> from timagetk.io import imread
>>> from timagetk.io.util import shared_dataset
>>> from timagetk.visu.mplt import grayscale_imshow
>>> seg_image = imread(shared_dataset("p58", "segmented")[1], LabelledImage)  # get 'p58-t1' shared labelled image
>>> int_image = imread(shared_dataset("p58", "intensity")[1])  # get 'p58-t1' shared intensity image
>>> # Example #1 - Project the labelled image along the Z-axis:
>>> alti = altimap_from_segmentation(seg_image)
>>> proj = altimap_surface_projection(seg_image, alti)
>>> f = grayscale_imshow([alti, proj], suptitle='Z-axis projection', title=["Altitude map", "Contour map"], cmap=['viridis', 'glasbey'], val_range='auto')

>>> # Example #2 - Project the intensity image along the Z-axis:
>>> alti = altimap_from_segmentation(seg_image)
>>> proj = altimap_surface_projection(int_image, alti)
>>> f = grayscale_imshow([alti, proj], suptitle='Z-axis projection', title=["Altitude map", "Intensity contour map"], cmap=['viridis', 'gray'], val_range=['auto', 'dtype'])

>>> # Example #3 - Project along the X-axis:
>>> alti = altimap_from_segmentation(seg_image, axis='x')
>>> proj = altimap_surface_projection(int_image, alti, axis='x')
>>> f = grayscale_imshow([alti, proj], suptitle='X-axis projection', title=["Altitude map", "Intensity contour map"], cmap=['viridis', 'gray'], val_range=['auto', 'dtype'])

>>> # Example #4 - Project along the Y-axis:
>>> alti = altimap_from_segmentation(seg_image, axis='y')
>>> proj = altimap_surface_projection(int_image, alti, axis='y')
>>> f = grayscale_imshow([alti, proj], suptitle='Y-axis projection', title=["Altitude map", "Intensity contour map"], cmap=['viridis', 'gray'], val_range=['auto', 'dtype'])

>>> # Example #5 - Project the labelled image along the Y-axis:
>>> alti = altimap_from_segmentation(seg_image, axis='y')
>>> proj = altimap_surface_projection(seg_image, alti, axis='y')
>>> f = grayscale_imshow([alti, proj], suptitle='Y-axis projection', title=["Altitude map", "Contour map"], cmap=['viridis', 'glasbey'], val_range='auto')
timagetk.algorithms.reconstruction.altitude_map(image, gaussian_sigma=1.0, height=3.0, threshold=None, orientation=None, **kwargs)[source]

Compute a surface projection image and associated altitude map.

Height values follows the slices indexing, so 0 for the first slice (hence the importance of orientation). They range in [0, axis_shape - 1] with axis_shape being the size of the Z projected axis. They are thus in voxel unit and not in real units.

Parameters:

  • image (timagetk.SpatialImage or timagetk.MultiChannelImage) – The intensity image to project.

  • gaussian_sigma (float, optional) – Sigma value, in real units, of the Gaussian filter applied to the image. Skip this step by setting the value to 0 or None. Defaults to 1. This value is divided by the image voxelsize.

  • height (float, optional) – Height parameter.

  • threshold (int or float, optional) – Instensity value threshold to create the mask. By default, try to guess it with guess_intensity_threshold.

  • orientation ({-1, 1}, optional) – Defines the orientation of the projection, try to guess it by default. If -1 we consider the z-axis starts at the bottom of the object and goes upward. If 1 we consider the z-axis starts at the top of the object and goes downward.

  • axis (str) – The axis to use for projection, defaults to “Z”.

  • channel (str) – If a MultiChannelImage is used as input image, select the channel to use with this algorithm.

  • altimap_smoothing (bool) – If True, default is False, smooth the altitude map obtained from the mask.

Returns:

timagetk.SpatialImage – YX 2D altitude map, height values starting from the first slice defined by the image orientation.

See also

timagetk.array_util.guess_intensity_threshold, timagetk.array_util.guess_image_orientation

Examples

>>> from timagetk.algorithms.reconstruction import altitude_map
>>> from timagetk.io import imread
>>> from timagetk.io.util import shared_dataset
>>> from timagetk.visu.mplt import grayscale_imshow
>>> img = imread(shared_dataset("p58")[0])
>>> # EXAMPLE #1 - Get z-altitude map:
>>> alti = altitude_map(img, threshold=45, orientation=1)
>>> grayscale_imshow(alti, title="Altitude map", cmap='viridis', val_range='auto')
>>> # EXAMPLE #2 - Smooth the returned z-altitude map with Gaussian filter:
>>> alti = altitude_map(img, threshold=45, orientation=1, altimap_smoothing=True)
>>> grayscale_imshow(alti, title="Altitude map", cmap='viridis', val_range='auto')
>>> # EXAMPLE #3 - Get y-altitude map:
>>> alti = altitude_map(img, threshold=45, orientation=1, axis="y")
>>> grayscale_imshow(alti, title="Altitude map", cmap='viridis', val_range='auto')
>>> # EXAMPLE #4 - Get -y-altitude map:
>>> alti = altitude_map(img, threshold=45, orientation=-1, axis="y")
>>> grayscale_imshow(alti, title="Altitude map", cmap='viridis', val_range='auto')
timagetk.algorithms.reconstruction.contour_from_segmentation(seg_image, background_id=1, orientation=None, **kwargs)[source]

Get the contour image from the labelled image, that is the background outline.

Parameters:

  • seg_image (timagetk.LabelledImage or timagetk.TissueImage3D) – The labelled image to use to detect the background position.

  • background_id (int, optional) – The label associated to the background, defaults to 1. If a TissueImage3D instance is given as seg_image, it will use its background attribute.

  • orientation ({None, -1, 1}, optional) – Defines the orientation of the projection. By default (None), try to guess it with guess_image_orientation. If -1 we consider the z-axis starts at the bottom of the object and goes upward. If 1 we consider the z-axis starts at the top of the object and goes downward.

  • iso_upsampling (bool) – If True, up-sample the image to the minimum voxelsize. This is usefull if you have an axis with a lower resolution, typically the z-axis.

  • interpolation ({'nearest', 'cellbased'}) – Interpolation method to use during iso_upsampling, if any. Defaults to ‘cellbased’.

  • cell_based_sigma (float) – Required when using “cellbased” interpolation method, sigma for cell-based interpolation. In voxel units!

  • connectivity ({1, 2, 3}) – The connectivity of the binary structuring element. Defaults to 1.

  • iterations (int) – The number of iterations of binary dilation by the structuring element. Defaults to 1.

Returns:

timagetk.SpatialImage – Contour image, that is the background outline.

See also

scipy.ndimage.generate_binary_structure, scipy.ndimage.binary_dilation

Examples

>>> import numpy as np
>>> from timagetk.algorithms.reconstruction import contour_from_segmentation
>>> from timagetk import LabelledImage
>>> from timagetk.io import imread
>>> from timagetk.io.util import shared_dataset
>>> from timagetk.visu.stack import stack_browser
>>> # - Get 'p58-t1' shared segmented image:
>>> seg_image = imread(shared_dataset("p58", "segmented")[1], LabelledImage, not_a_label=0)
>>> # - Get the contour:
>>> contour = contour_from_segmentation(seg_image, orientation=1, connectivity=1, iterations=2)
>>> b = stack_browser(contour, title="Contour map", cmap='gray', val_range=[0, 1])
>>> # - Get the labelled contour:
>>> labelled_contour = LabelledImage(seg_image * contour)
>>> b = stack_browser(labelled_contour, title="Labelled contour map", cmap='glasbey', val_range='auto', colorbar=True)
timagetk.algorithms.reconstruction.end_margin(img, width, axis=None, null_value=0)[source]

Set selected end margin axis to a null value.

Parameters:

  • img (numpy.ndarray) – A 3D array to modify.

  • width (int) – Size of the margin.

  • axis (int, optional) – Axis along which the margin is edited. Edit them all by default.

  • null_value (int, optional) – The null value to use for the end margins, 0 by default.

Returns:

numpy.ndarray – Input array with the end margin set to a null value.

See also

stroke

Example

>>> from timagetk.algorithms.reconstruction import end_margin
>>> from timagetk.array_util import random_spatial_image
>>> from timagetk import SpatialImage
>>> from timagetk.visu.stack import orthogonal_view
>>> img = random_spatial_image((5, 10, 10), dtype='uint8')
>>> img_m = SpatialImage(end_margin(img, 1), voxelsize=img.get_voxelsize())
>>> orthogonal_view(img, cmap='viridis', suptitle="Original")
>>> orthogonal_view(img_m, cmap='viridis', suptitle="+end_margin")
timagetk.algorithms.reconstruction.im2surface(image, threshold_value=45, front_half_only=False, **kwargs)[source]

Contour altitude map and intensity projection.

This function computes a surfacic view of the object, according to a revisited version of the method described in [Barbier].

Parameters:

  • image (timagetk.SpatialImage or timagetk.MultiChannelImage) – Image to use for altitude map and intensity projection.

  • threshold_value (int or float) – Consider intensities greater than this threshold.

  • front_half_only (bool) – If True, only consider the first half of all slices along the last physical axis. Defaults to False.

  • altimap_smoothing (bool) – If True, default is False, smooth the altitude map obtained from the mask.

Returns:

  • timagetk.SpatialImage – Contour intensity projection.

  • timagetk.SpatialImage – Altitude of maximum intensity projection.

References

[Barbier]

Barbier de Reuille, P. , Bohn‐Courseau, I. , Godin, C. and Traas, J. (2005), A protocol to analyse cellular dynamics during plant development. The Plant Journal, 44: 1045-1053. DOI.

Example

>>> from timagetk.io.util import shared_dataset
>>> from timagetk.io import imread
>>> from timagetk.algorithms.reconstruction import im2surface
>>> from timagetk.visu.stack import stack_panel
>>> from timagetk.visu.mplt import grayscale_imshow
>>> # - Get 'p58' shared intensity image:
>>> image = imread(shared_dataset("p58")[0])
>>> mip_img, alt_img = im2surface(image, 60, altimap_smoothing=True)
>>> grayscale_imshow([mip_img, alt_img], suptitle="im2surface", title=["Contour Projection", "Altitude Map"], val_range='auto', cmap=['gray', 'viridis'])
timagetk.algorithms.reconstruction.image_surface_projection(image, gaussian_sigma=1.0, height=3.0, threshold=None, orientation=None, **kwargs)[source]

Compute a surface projection image and associated altitude map.

Parameters:

  • image (timagetk.SpatialImage or timagetk.MultiChannelImage) – The intensity image to project.

  • gaussian_sigma (float, optional) – Sigma value, in real units, of the Gaussian filter applied to the image. Skip this step by setting the value to 0 or None. Defaults to 1. This value is divided by the image voxelsize.

  • height (float, optional) – Height parameter.

  • threshold (int or float, optional) – Intensity value threshold to create the mask. By default, try to guess it with guess_intensity_threshold.

  • orientation ({None, -1, 1}, optional) – Defines the orientation of the projection. By default (None), try to guess it with guess_image_orientation. If -1 we consider the z-axis starts at the bottom of the object and goes upward. If 1 we consider the z-axis starts at the top of the object and goes downward.

  • axis (str) – The axis to use for projection, defaults to “Z”.

  • iso_upsampling (bool) – If True (default), up-sample the image to the minimum voxelsize. This is usefull if you have an axis with a lower resolution, typically the z-axis.

  • interpolation ({'linear', 'cspline'}) – Interpolation method to use during iso_upsampling, if any. Defaults to ‘linear’.

  • channel (str) – If a MultiChannelImage is used as input image, select the channel to use with this algorithm.

  • altimap_smoothing (bool) – If True, default is False, smooth the altitude map obtained from the mask.

See also

timagetk.array_util.guess_intensity_threshold, timagetk.array_util.guess_image_orientation

Returns:

  • timagetk.SpatialImage – YX 2D intensity projection image.

  • timagetk.SpatialImage – YX 2D altitude map, height values starting from the first slice defined by the image orientation.

Examples

>>> import numpy as np
>>> from timagetk.algorithms.reconstruction import image_surface_projection
>>> from timagetk import SpatialImage
>>> from timagetk.io import imread
>>> from timagetk.io.util import shared_dataset
>>> from timagetk.visu.mplt import grayscale_imshow
>>> # - Get 'p58' shared intensity image:
>>> img = imread(shared_dataset("p58")[0])
>>> # -- EXAMPLE #1 - Use automatic thresholding
>>> proj, alti = image_surface_projection(img, orientation=1)
>>> grayscale_imshow([proj, SpatialImage(alti)], title=["Intensity projection map", "Altitude map"], cmap=['gray', 'viridis'], val_range=['type', 'auto'])
>>> # -- EXAMPLE #2 - Use 60-th percentile as threshold
>>> threshold = np.percentile(img, 60)
>>> print(f"60-th percentile intensity threshold: {threshold}")
>>> proj, alti = image_surface_projection(img, threshold=threshold, orientation=1)
>>> grayscale_imshow([proj, SpatialImage(alti)], title=["Intensity projection map", "Altitude map"], cmap=['gray', 'viridis'], val_range=['type', 'auto'])
>>> # -- EXAMPLE #3 - Project along x-axis
>>> proj, alti = image_surface_projection(img, orientation=1, axis='x')
>>> grayscale_imshow([proj, SpatialImage(alti)], suptitle="X-axis projection", title=["Intensity projection map", "Altitude map"], cmap=['gray', 'viridis'], val_range=['type', 'auto'])
>>> # -- EXAMPLE #4 - Project along y-axis
>>> proj, alti = image_surface_projection(img, orientation=1, axis='y')
>>> grayscale_imshow([proj, SpatialImage(alti)], suptitle="Y-axis projection", title=["Intensity projection map", "Altitude map"], cmap=['gray', 'viridis'], val_range=['type', 'auto'])
>>> # -- EXAMPLE #5 - Project along -y-axis
>>> y_shape = img.get_shape('y')
>>> threshold = np.percentile(img[:, y_shape//2: y_shape-1, :], 70)  # brighter on this side, so higher threshold
>>> proj, alti = image_surface_projection(img, threshold=threshold, orientation=-1, axis='y')
>>> grayscale_imshow([proj, SpatialImage(alti)], suptitle="Y-axis projection", title=["Intensity projection map", "Altitude map"], cmap=['gray', 'viridis'], val_range=['type', 'auto'])
timagetk.algorithms.reconstruction.padding(image, z_pad=(0, 0), y_pad=(0, 0), x_pad=(0, 0), pad_value=0)[source]

Set an outline of values to an image.

Parameters:

  • image (timagetk.SpatialImage or timagetk.LabelledImage or timagetk.TissueImage3D or timagetk.MultiChannelImage or numpy.ndarray) – A 3D image to modify by applying padding. In case of a numpy array, consider the axes to be ZYX sorted.

  • z_pad ((int, int), optional) – Padding to apply to the beginning and end of the z-axis. Defaults to (0, 0).

  • y_pad ((int, int), optional) – Padding to apply to the beginning and end of the y-axis. Defaults to (0, 0).

  • x_pad ((int, int), optional) – Padding to apply to the beginning and end of the x-axis. Defaults to (0, 0).

  • pad_value (int, optional) – The padding value to use. Defaults to 0.

Returns:

timagetk.SpatialImage or timagetk.LabelledImage or timagetk.TissueImage3D or timagetk.MultiChannelImage or numpy.ndarray – Input image with the padding.

See also

end_margins

Example

>>> from timagetk.array_util import random_spatial_image
>>> from timagetk.algorithms.reconstruction import padding
>>> from timagetk import SpatialImage
>>> from timagetk.visu.stack import orthogonal_view
>>> img = random_spatial_image((5, 8, 8), dtype='uint8')
>>> # Example: Add a slice of null value on top of stack:
>>> img_s = padding(img, (1, 0))
>>> print(img.shape)
(5, 8, 8)
>>> print(img_s.shape)
(6, 8, 8)
>>> orthogonal_view(img, cmap='viridis', suptitle="Original")
>>> orthogonal_view(img_s, cmap='viridis', suptitle="Padded image")
timagetk.algorithms.reconstruction.project_segmentation(seg_img, axis='z', orientation=1, background=1, return_background=True)[source]

Project a segmented image for given axis and direction.

This method consider the background as transparent and mimic a surfacic projection.

Parameters:

  • seg_img (timagetk.LabelledImage) – The segmented image to project.

  • axis ({'x', 'y', 'z'}, optional) – The axis to use for projection, default to z.

  • orientation ({1, -1}, optional) – Direction in which to move along the selected axis, 1 by default 1 follow increasing values along selected axis. -1 follow decreasing values along selected axis.

  • background (int, optional) – The label value associated to the background, 1 by default

  • return_background (bool, optional) – If True (default), add the background in the returned projection, else leave it to 0

Returns:

timagetk.LabelledImage – The projected segmentation along given axis & direction

Notes

The returned LabelledImage declare 0 as not_a_label. Voxel-sizes, origin and metadata are kept from original segmented image.

Examples

>>> from timagetk.algorithms.reconstruction import project_segmentation
>>> from timagetk import TissueImage3D
>>> from timagetk.io import imread
>>> from timagetk.io.util import shared_dataset
>>> # Get the 'sphere' example labelled image:
>>> seg_img = imread(shared_dataset('sphere', 'segmented')[0], TissueImage3D, not_a_label=0, background=1)
>>> # Create a 'top-view' z-projection:
>>> z_proj = project_segmentation(seg_img)
>>> from timagetk.visu.util import greedy_colormap
>>> from timagetk.visu.mplt import grayscale_imshow
>>> grayscale_imshow(z_proj, cmap=greedy_colormap(z_proj), title="Segmentation projection map")
>>> # Get a realistic labelled image 'p58'
>>> seg_img = imread(shared_dataset('p58', 'segmented')[0], TissueImage3D, not_a_label=0, background=1)
>>> x_projs = [project_segmentation(seg_img,'x',1), project_segmentation(seg_img,'x',-1)]
>>> y_projs = [project_segmentation(seg_img,'y',1), project_segmentation(seg_img,'y',-1)]
>>> grayscale_imshow([z_proj] + x_projs + y_projs, suptitle="Segmentation projection map", title=["z-proj", "x-proj", "-x-proj", "y-proj", "-y-proj"], cmap='glasbey', val_range='auto', max_per_line=5)
timagetk.algorithms.reconstruction.surface2im(points, altitude_map, vz=1)[source]

Convert points from contour projection to the real world.

Parameters:

  • points (list) – List of points from the maximum intensity projection.

  • altitude_map (timagetk.SpatialImage) – Altitude of maximum intensity projection.

Returns:

list – List of points in the real world

Examples

>>> from timagetk.algorithms.reconstruction import surface2im
>>> from timagetk.algorithms.reconstruction import altitude_map
>>> from timagetk.io import imread
>>> from timagetk.io.util import shared_dataset
>>> ref_path = shared_dataset("p58")[0]
>>> ref_img = imread(ref_path)
>>> alti = altitude_map(ref_img, threshold=45, orientation=1)
>>> xyz = surface2im([[50, 50], [200, 200]], alti, ref_img.get_voxelsize('z'))
>>> print(xyz)
[[20.03200054168701, 20.03200054168701, 63.7017617225647], [80.12800216674805, 80.12800216674805, 19.23072052001953]]