Parametric maps#

This notebook demonstrates how parametric maps can be made. In such parametric images, pixel intensity corresponds to measurements of the objects, for example area.

import numpy as np
import pyclesperanto_prototype as cle
from skimage.io import imread, imsave
from skimage.measure import regionprops, regionprops_table
from skimage.util import map_array
from napari_simpleitk_image_processing import label_statistics

Starting point for drawing parametric maps is always a label image.

binary = cle.artificial_objects_2d()
labels = cle.voronoi_otsu_labeling(binary, spot_sigma=5)

cle.imshow(labels, labels=True)
../_images/parametric_maps_3_0.png

Parametric maps using scikit-image#

You can also compute your own measurement for each object and then visualize it in a parametric map image. Consider using scikit-image’s regionprops() for computing the measurements.

statistics_table = regionprops_table(cle.pull(labels), properties=('label', 'area',))

Area map#

remapped = map_array(
        cle.pull(labels),
        statistics_table['label'],
        statistics_table['area'],
        )

cle.imshow(remapped, colorbar=True, color_map="jet")
../_images/parametric_maps_7_0.png

Parametric maps in pyclesperanto#

Alternatively, we can use the function cle.replace_intensities for making these parametric maps. It expects a label image and a list of measurements ordered by label, starting with a measurement corresponding to background (index 0). In the list of measurements, no label should be missing. Consider to relabel the label image sequentially if necessary.

statistics = regionprops(cle.pull(labels))
# we're appending a [0] for the "measurement" of the background and attach all values behind
area = [0] + [s.area for s in statistics]

area_map2 = cle.replace_intensities(labels, area)

cle.imshow(area_map2, colorbar=True, color_map="jet")
../_images/parametric_maps_11_0.png

Pixel count map#

pyclesperanto comes with some maps built-in. For example the pixel count map derived from a label image expresses area or volume of objects in colour.

pixel_count_map = cle.label_pixel_count_map(labels)

cle.imshow(pixel_count_map, color_map='jet', colorbar=True)
/Users/haase/code/pyclesperanto_prototype/pyclesperanto_prototype/_tier9/_statistics_of_labelled_pixels.py:283: RuntimeWarning: invalid value encountered in true_divide
  region_props['mean_max_distance_to_centroid_ratio'] = region_props['max_distance_to_centroid'] / region_props[
/Users/haase/code/pyclesperanto_prototype/pyclesperanto_prototype/_tier9/_statistics_of_labelled_pixels.py:285: RuntimeWarning: invalid value encountered in true_divide
  region_props['mean_max_distance_to_mass_center_ratio'] = region_props['max_distance_to_mass_center'] / region_props[
../_images/parametric_maps_13_1.png

Parametric maps using SimpleITK-based measurements#

Furthermore, also SimpleITK comes with quantitative measurements for labeled images. For convenience reasons, we will use the scriptable napari plugin napari-simpleitk-image-processing for deriving measurements.

statistics_sitk = label_statistics(labels, labels, size=True, shape=True, perimeter=True, intensity=False)

print(statistics_sitk.keys())
dict_keys(['label', 'elongation', 'feret_diameter', 'flatness', 'roundness', 'equivalent_ellipsoid_diameter_0', 'equivalent_ellipsoid_diameter_1', 'equivalent_spherical_perimeter', 'equivalent_spherical_radius', 'number_of_pixels', 'number_of_pixels_on_border', 'perimeter', 'perimeter_on_border', 'perimeter_on_border_ratio'])

Number of pixels map#

# we're appending a [0] for the "measurement" of the background and attach all values behind
number_of_pixels = [0] + statistics_sitk['number_of_pixels']

number_of_pixels_map = cle.replace_intensities(labels, number_of_pixels)

cle.imshow(number_of_pixels_map, colorbar=True, color_map="jet")
../_images/parametric_maps_17_0.png

Extension ratio map#

The extension ratio is a shape descriptor derived from the maximum distance of pixels to their object’s centroid divided by the average distance of pixels to the centroid. It alllows to differentiate elongated from roundish objects.

extension_ratio_map = cle.extension_ratio_map(labels)

cle.imshow(extension_ratio_map, 
           color_map='jet', 
           colorbar=True, 
           min_display_intensity=1.5, 
           max_display_intensity=2.5)
../_images/parametric_maps_19_0.png

Mean / minimum / maximum / standard-deviation intensity map#

If we use additionally an intensity image, we can also produce parametric images showing intensity measurements.

blobs = cle.scale(imread('../../data/blobs.tif'), factor_x=2, factor_y=2, auto_size=True)

cle.imshow(blobs)
../_images/parametric_maps_21_0.png
mean_intensity_map = cle.label_mean_intensity_map(blobs, labels)
cle.imshow(mean_intensity_map, 
           color_map='jet', 
           colorbar=True,
           min_display_intensity=100, 
           max_display_intensity=250)
../_images/parametric_maps_22_0.png
maximum_intensity_map = cle.maximum_intensity_map(blobs, labels)
cle.imshow(maximum_intensity_map, 
           color_map='jet', 
           colorbar=True, 
           min_display_intensity=200, 
           max_display_intensity=270)
../_images/parametric_maps_23_0.png
stddev_intensity_map = cle.standard_deviation_intensity_map(blobs, labels)
cle.imshow(stddev_intensity_map, color_map='jet', colorbar=True)
../_images/parametric_maps_24_0.png

Aspect ratio map#

# we're appending a [0] for the "measurement" of the background and attach all values behind
aspect_ratio = [0]

for s in statistics:
    if s.minor_axis_length:
        aspect_ratio.append(s.major_axis_length / s.minor_axis_length)
    else:
        aspect_ratio.append(0) # note: an aspect ratio of 0 is an "invalid" value
aspect_ratio_map = cle.replace_intensities(labels, aspect_ratio)

cle.imshow(aspect_ratio_map, colorbar=True, color_map="jet")
../_images/parametric_maps_27_0.png

Elongation map#

# we're appending a [0] for the "measurement" of the background and attach all values behind
eccentricity = [0] + [s.eccentricity for s in statistics]

eccentricity_map = cle.replace_intensities(labels, eccentricity)

cle.imshow(eccentricity_map, colorbar=True, color_map="jet")
../_images/parametric_maps_29_0.png

Extent map#

# we're appending a [0] for the "measurement" of the background and attach all values behind
extent = [0] + [s.extent for s in statistics]

extent_map = cle.replace_intensities(labels, extent)

cle.imshow(extent_map, colorbar=True, color_map="jet")
../_images/parametric_maps_31_0.png

Feret diameter map#

# we're appending a [0] for the "measurement" of the background and attach all values behind
feret_diameter = [0] + [s.feret_diameter_max for s in statistics]

feret_diameter_map = cle.replace_intensities(labels, feret_diameter)

cle.imshow(feret_diameter_map, colorbar=True, color_map="jet")
../_images/parametric_maps_33_0.png

Orientation map#

# we're appending a [0] for the "measurement" of the background and attach all values behind
orientation = [0] + [s.orientation for s in statistics]

orientation_map = cle.replace_intensities(labels, orientation)

cle.imshow(orientation_map, colorbar=True, color_map="jet")
../_images/parametric_maps_35_0.png

Perimeter map#

# we're appending a [0] for the "measurement" of the background and attach all values behind
perimeter_ski = [0] + [s.perimeter for s in statistics]

perimeter_ski_map = cle.replace_intensities(labels, perimeter_ski)

cle.imshow(perimeter_ski_map, colorbar=True, color_map="jet")
../_images/parametric_maps_37_0.png

Solidity map#

# we're appending a [0] for the "measurement" of the background and attach all values behind
solidity = [0] + [s.solidity for s in statistics]

solidity_map = cle.replace_intensities(labels, solidity)

cle.imshow(solidity_map, colorbar=True, color_map="jet")
../_images/parametric_maps_39_0.png

Elongation map#

# we're appending a [0] for the "measurement" of the background and attach all values behind
elongation = [0] + statistics_sitk['elongation']

elongation_map = cle.replace_intensities(labels, elongation)

cle.imshow(elongation_map, colorbar=True, color_map="jet")
../_images/parametric_maps_41_0.png

Flatness map#

# we're appending a [0] for the "measurement" of the background and attach all values behind
flatness = [0] + statistics_sitk['flatness']

flatness_map = cle.replace_intensities(labels, flatness)

cle.imshow(flatness_map, colorbar=True, color_map="jet")
../_images/parametric_maps_43_0.png

Roundness map#

# we're appending a [0] for the "measurement" of the background and attach all values behind
roundness = [0] + statistics_sitk['roundness']

roundness_map = cle.replace_intensities(labels, roundness)

cle.imshow(roundness_map, colorbar=True, color_map="jet")
../_images/parametric_maps_45_0.png

Perimeter map#

# we're appending a [0] for the "measurement" of the background and attach all values behind
perimeter_sitk = [0] + statistics_sitk['perimeter']

perimeter_sitk_map = cle.replace_intensities(labels, perimeter_sitk)

cle.imshow(perimeter_sitk_map, colorbar=True, color_map="jet")
../_images/parametric_maps_47_0.png

Quality assurance#

If you generate a parametric image with the “label” column, the parametric image should actually be equal to the label input image.

# we're appending a [0] for the "measurement" of the background and attach all values behind
label = [0] + statistics_sitk['label']

label_map = cle.replace_intensities(labels, label)

cle.imshow(label_map, labels=True)
../_images/parametric_maps_49_0.png
label_difference = labels - label_map

cle.imshow(label_difference)
../_images/parametric_maps_50_0.png
label_difference.min(), label_difference.max()
(0.0, 0.0)

When comparing perimeter measurements from scikit-image and SimpleITK, we could see small differences. Here we can visualize which objects are affected and to what degree.

perimeter_difference = perimeter_sitk_map - perimeter_ski_map

cle.imshow(perimeter_difference, colorbar=True, color_map="jet")
../_images/parametric_maps_53_0.png