Mentions légales du service

Skip to content
Snippets Groups Projects

Python examples

Last edited by LEGRAND Jonathan

Blockmatching:

Performs image registration using two intensity images, one as a reference, the other as a floating image (to register) and returns the transformation and the registered image according to the selected registration method and parameters.

from timagetk.util import data_path
from timagetk.io import imread
from timagetk.algorithms.blockmatching import blockmatching
from timagetk.visu.mplt import grayscale_imshow

flo_path = data_path('p58-t0_INT_down_interp_2x.inr.gz')
floating_image = imread(flo_path)
ref_path = data_path('p58-t1_INT_down_interp_2x.inr.gz')
reference_image = imread(ref_path)

# - RIGID registration
param_str_2 = '-trsf-type rigid -py-ll 2'
trsf_rig, res_rig = blockmatching(floating_image, reference_image,
                                  param_str_2=param_str_2)

# - Display rigid registration effect
img2plot = [floating_image, reference_image, res_rig, reference_image]
img_titles = ["t0", "t1", "Registered t0 on t1", "t1"]
grayscale_imshow(img2plot, "Effect of rigid registration", img_titles, vmin=0, vmax=255, max_per_line=2)

# - DEFORMABLE registration
param_str_2 = '-trsf-type vectorfield -py-ll 1'
trsf_def, res_def = blockmatching(floating_image, reference_image,
                                  init_result_trsf=trsf_rig,
                                  param_str_2=param_str_2)

# - Display deformable registration effect
img2plot = [floating_image, reference_image, res_def, reference_image]
img_titles = ["t0", "t1", "Registered t0 on t1", "t1"]
grayscale_imshow(img2plot, "Effect of deformable registration", img_titles, vmin=0, vmax=255, max_per_line=2)

This corresponds to the following line commands:

blockmatching -ref p58-t1_INT_down_interp_2x.inr.gz -flo p58-t0_INT_down_interp_2x.inr.gz \
  -trsf-type rigid -py-ll 2 -res result-rigid.inr -res-trsf result-rigid.trsf

# 'result-rigid.inr' is the file name of the result image, ie the floating image
# (specified by '-flo') resampled in the same frame than the reference image
# (specified by '-ref'). It has then the same dimensions and voxel sizes.
# 'result-rigid.trsf' is the file name of the result transformation 
# (here a linear transformation, hence a matrix, hence a file of small size)
#
# Note that, 'result-rigid.trsf' being known, 'result-rigid.inr' can be re-computed with
# 'applyTrsf p58-t0_INT_down_interp_2x.inr.gz result-rigid.inr \
#    -trsf result-rigid.trsf -ref p58-t1_INT_down_interp_2x.inr.gz'

blockmatching -ref p58-t1_INT_down_interp_2x.inr.gz -flo p58-t0_INT_down_interp_2x.inr.gz \
  -trsf-type vectorfield -py-ll 1 -res result-def.inr -res-def result-def.trsf \
  -init-trsf result-rigid.trsf -composition-with-initial

# 'result-def.inr' is the file name of the result image
# 'result-def.trsf' is the file name of the result transformation 
# (here a non-linear transformation, hence a vectorial image, hence a file of large size)
# Note: using an initial transformation (here 'result-rigid.trsf') can be done by two means
# - either '-init-trsf result-rigid.trsf -composition-with-initial'
# - or '-init-res-trsf  result-rigid.trsf'
# and behaviors may be different. I (GM) am in favor of the first solution, but one should
# look what has been done in the above used wrapping

Linear Filtering:

from timagetk.util import data_path
from timagetk.io import imread
from timagetk.algorithms import linearfilter
img = imread(data_path('p58-t0_INT_down_interp_2x.inr.gz'))

# - Apply linear filtering with default parameters:
# -- Gaussian smoothing, sigma=1.0 in voxel units, param_str_1 = "-smoothing -sigma 1.0"
out_img = linearfilter(img)

grayscale_imshow([img, out_img], "Gaussian smoothing, sigma=1.0 voxel", ["Original", "Smoothed"], vmin=0, vmax=255)

# - Apply linear filtering with other parameters:
# -- Laplacian operator with sigma=0.6 voxels
param_str_2 = '-laplacian -sigma 0.6'
out_img = linearfilter(img, param_str_2=param_str_2)

grayscale_imshow([img.get_z_slice(40), out_img.get_z_slice(40)], "Laplacian operator, sigma=0.6 voxel", ["Original", "Laplacian"], vmin=0, vmax=255)

This corresponds to the following line commands:

linearFilter p58-t0_INT_down_interp_2x.inr.gz result-image.inr -smoothing -sigma 1.0

# it seems that the python wrapping has already a default behavior with 
# '-smoothing -sigma 1.0'. It is the reason of the two strings in all the functions
# named 'API_something()' in the C code: the first is dedicated to some standard behavior
# while the second is for the user options (that somehow overwrite the standard behavior)

linearFilter p58-t0_INT_down_interp_2x.inr.gz result-image.inr -laplacian -sigma 0.6

Morphological operators on grayscale image:

from timagetk.util import data_path
from timagetk.io import imread
from timagetk.algorithms import morpho
from timagetk.visu.mplt import grayscale_imshow
input_image = imread(data_path('p58-t0_INT_down_interp_2x.inr.gz'))

# - Apply morphological dilation on grayscale image by a sphere of radius 1
out_image = morpho(input_image, param_str_1="-dilation -sphere -radius 1 -iterations 1")

# - Display
grayscale_imshow([input_image.get_z_slice(50), out_image.get_z_slice(50)], "Morphological dilation by a sphere of radius 1", ["Original slice", "Dilated slice"])

This corresponds to the following line commands:

morpho p58-t0_INT_down_interp_2x.inr.gz result-image.inr -dilation -sphere -radius 1 \
  -iterations 1

# '-sphere' means that a true sphere will be used as structuring element
# since its radius is 1, it makes no sense. Using either the 6- or the 18-connectivity will 
# yield the same result at a lower computational cost

Apply Trsf:

from timagetk.util import data_path
from timagetk.io import imread
from timagetk.algorithms.blockmatching import blockmatching
from timagetk.algorithms.trsf import apply_trsf
from timagetk.visu.mplt import grayscale_imshow

flo_path = data_path('p58-t0_INT_down_interp_2x.inr.gz')
floating_image = imread(flo_path)
ref_path = data_path('p58-t1_INT_down_interp_2x.inr.gz')
reference_image = imread(ref_path)

# - Compute RIGID registration, get only the transformation:
param_str_2 = '-trsf-type rigid -py-ll 2'
trsf_rig, _ = blockmatching(floating_image, reference_image, param_str_2=param_str_2)

# - Apply the computed transformation to the original image, using the reference as template:
res_rig = apply_trsf(floating_image, trsf_rig, template_img=reference_image)

# - Display rigid registration effect
img2plot = [floating_image, reference_image, res_rig, reference_image]
img_titles = ["t0", "t1", "Registered t0 on t1", "t1"]
grayscale_imshow(img2plot, "Effect of rigid registration", img_titles, vmin=0, vmax=255, max_per_line=2)

This corresponds to the following line commands:

blockmatching -ref p58-t1_INT_down_interp_2x.inr.gz -flo p58-t0_INT_down_interp_2x.inr.gz \
  -trsf-type rigid -py-ll 2 -res-trsf result-rigid.trsf

applyTrsf p58-t0_INT_down_interp_2x.inr.gz result-rigid.inr \
  -trsf result-rigid.trsf -ref p58-t1_INT_down_interp_2x.inr.gz

Mean image:

from timagetk.util import data_path
from timagetk.io import imread
from timagetk.algorithms.blockmatching import blockmatching
from timagetk.algorithms.mean_images import mean_images
from timagetk.visu.mplt import grayscale_imshow

flo_path = data_path('p58-t0_INT_down_interp_2x.inr.gz')
floating_image = imread(flo_path)
ref_path = data_path('p58-t1_INT_down_interp_2x.inr.gz')
reference_image = imread(ref_path)

# - Compute RIGID registration:
param_str_2 = '-trsf-type rigid -py-ll 2'
trsf_rig, res_rig = blockmatching(floating_image, reference_image, param_str_2=param_str_2)

max_im = mean_images([reference_image, res_rig], param_str_2="-max")
min_im = mean_images([reference_image, res_rig], param_str_2="-min")
mean_im = mean_images([reference_image, res_rig], param_str_2="-mean")
median_im = mean_images([reference_image, res_rig], param_str_2="-median")

# - Display image averaging effect
img2plot = [floating_image, reference_image, res_rig, reference_image]
img_titles = ["t0", "t1", "Registered t0 on t1", "t1"]
grayscale_imshow([reference_image, res_rig, max_im, min_im, mean_im, median_im], "Image averaging", ["Image1", "Image2", "max", "min", "mean", "median"], vmin=0, vmax=255, max_per_line=2)

This corresponds to the following line commands:

blockmatching -ref p58-t1_INT_down_interp_2x.inr.gz -flo p58-t0_INT_down_interp_2x.inr.gz \
  -trsf-type rigid -py-ll 2 -res result-rigid.inr -res-trsf result-rigid.trsf

meanImages -images p58-t1_INT_down_interp_2x.inr.gz result-rigid.inr \
  -res image-max.inr -max
meanImages -images p58-t1_INT_down_interp_2x.inr.gz result-rigid.inr \
  -res image-min.inr -min
meanImages -images p58-t1_INT_down_interp_2x.inr.gz result-rigid.inr \
  -res image-mean.inr -mean
meanImages -images p58-t1_INT_down_interp_2x.inr.gz result-rigid.inr \
  -res image-median.inr -median

# here taking the median value out of two values makes no sense, of course.

Create trsf:

from timagetk.algorithms.trsf import create_trsf

# - Generate identity transformation matrix:
id_trsf = create_trsf(param_str_1='-identity')
# - Generate random transformation matrix:
rand_trsf = create_trsf(param_str_1='-random')

# - Print the transformation matrix:
id_trsf.c_display('identity')
rand_trsf.c_display('random')

This corresponds to the following line commands:

createTrsf identity.trsf -identity
createTrsf random.trsf -random

# Note the default type for created transformations is 'affine'
# If one want a reproducible random, the random seed can be specified with '-srandom %ld'

Invert trsf:

from timagetk.algorithms.trsf import inv_trsf

# - Generate random transformation matrix:
rand_trsf = create_trsf(param_str_1='-random')
# - Invert it:
inv_rand_trsf = inv_trsf(rand_trsf)

# - Print the original and inverted transformation matrix:
rand_trsf.c_display('random')
inv_rand_trsf.c_display('inverted random')

This corresponds to the following line commands:

createTrsf random.trsf -random
invTrsf random.trsf inv-random.trsf

Mean trsf:

Command lines:

./createTrsf random1.trsf -random
./createTrsf random2.trsf -random
echo './random1.trsf
./random2.trsf' > trsf_list.txt
./meanTrsfs trsf_list.txt -res mean_trsf.trsf -trsf-type affine3D -mean

Compose trsf:

from timagetk.algorithms.trsf import compose_trsf
from timagetk.algorithms.trsf import inv_trsf

# - Generate random transformation matrix:
rand_trsf = create_trsf(param_str_1='-random')
# - Invert it:
inv_rand_trsf = inv_trsf(rand_trsf)

# - Compose the two (should return identity!)
res_trsf = compose_trsf([rand_trsf, inv_rand_trsf])

# - Print the composed transformation matrix:
res_trsf.c_display('composed')

This corresponds to the following line commands:

createTrsf random.trsf -random
invTrsf random.trsf inv-random.trsf
composeTrsf -res res.trsf -trsfs random.trsf inv-random.trsf

# Note: the order of transformations after '-trsfs' is important

regionalext:

from timagetk.util import data_path
from timagetk.io import imread
from timagetk.algorithms import regionalext
from timagetk.visu.mplt import grayscale_imshow
input_image = imread(data_path('p58-t0_INT_down_interp_2x.inr.gz'))

# - Compute the minimal height map
output_image = regionalext(input_image, param_str_2="-minima -connectivity 26 -h 5")

# - Display it:
grayscale_imshow(output_image.get_z_slice(50), "Minimal height map", "Height map", vmin=0, vmax=5)

This corresponds to the following line commands:

regionalext p58-t0_INT_down_interp_2x.inr.gz hminima-image.inr \
  -minima -connectivity 26 -h 5

connexe:

from timagetk.util import data_path
from timagetk.io import imread
from timagetk.algorithms import regionalext
from timagetk.algorithms import connexe
from timagetk.visu.mplt import grayscale_imshow
input_image = imread(data_path('p58-t0_INT_down_interp_2x.inr.gz'))

# - Compute the minimal height map
h_image = regionalext(input_image, param_str_2="-minima -connectivity 26 -h 5")

# - Get the connected components
cx_image = connexe(h_image, param_str_2="-low-threshold 1 -high-threshold 5 -labels -connectivity 26")

# - Display it:
grayscale_imshow([input_image.get_z_slice(50), cx_image.get_z_slice(50)], "Local minima detection in 3D", ["Original slice", "Detected mimina"])

This corresponds to the following line commands:

regionalext p58-t0_INT_down_interp_2x.inr.gz hminima-image.inr \
  -minima -connectivity 26 -h 5
connexe hminima-image.inr seeds-image.inr -low-threshold 1 -high-threshold 5 \
  -labels -connectivity 26

# Note that 'connexe' here does a hysteresis thresholding with the specifications
# of the two thresholds '-low-threshold 1 -high-threshold 5' (5 is the h value of
# the h-minima computation)
# Since we are interested in labeled connected components, the '-labels' option
# forces their labeling.

watershed:

Make use of regionalext & connexe to generate seed image. Often also require linearFilter to filter image before seed detection and watershed (not used in the following example).

from timagetk.util import data_path
from timagetk.io import imread
from timagetk.algorithms import regionalext
from timagetk.algorithms import connexe
from timagetk.algorithms import watershed
from timagetk.visu.mplt import grayscale_imshow
input_image = imread(data_path('p58-t0_INT_down_interp_2x.inr.gz'))

# - Compute the minimal height map
h_image = regionalext(input_image, param_str_2="-minima -connectivity 26 -h 5")

# - Get the connected components as seeds:
cx_image = connexe(h_image, param_str_2="-low-threshold 1 -high-threshold 5 -labels -connectivity 26")

# - Use intensity image and seeds to performs watershed:
lab_img = watershed(input_image, cx_image, param_str_1='-labelchoice most')

grayscale_imshow([input_image.get_z_slice(50), lab_img.get_z_slice(50)], "3D watershed", ["Original slice", "Labelled slice"])

This corresponds to the following line commands:

regionalext p58-t0_INT_down_interp_2x.inr.gz hminima-image.inr \
  -minima -connectivity 26 -h 5
connexe hminima-image.inr seeds-image.inr -low-threshold 1 -high-threshold 5 \
  -labels -connectivity 26
watershed -seeds seeds-image.inr -gradient p58-t0_INT_down_interp_2x.inr.gz \
  watershed-image.inr -labelchoice most

Morphological operators on labeled cell image:

from timagetk.util import data_path
from timagetk.io import imread
from timagetk.algorithms import cell_filter
from timagetk.visu.mplt import grayscale_imshow
input_image = imread(data_path('p58-t0_SEG_down_interp_2x.inr.gz'))

# - Get the range of values
import numpy as np
uniq = np.unique(input_image)

# - Apply morphological erosion on grayscale image by a sphere of radius 1
out_image = cell_filter(input_image, param_str_1="-erosion -sphere -radius 1 -iterations 1")

# - Display
grayscale_imshow([input_image.get_z_slice(50), out_image.get_z_slice(50)], "Morphological erosion by a sphere of radius 1", ["Original slice", "Dilated slice"], vmin=uniq[0], vmax=uniq[-1])

This corresponds to the following line commands:

cellfilter p58-t0_SEG_down_interp_2x.inr.gz cellfilter-image.inr \
  -erosion -sphere -radius 1 -iterations 1

Comments

    Please register or sign in to add a comment.