How to - Iterative multi-angle images fusion

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from timagetk.algorithms.blockmatching import blockmatching
from timagetk.tasks.fusion import iterative_fusion
from timagetk.algorithms.pointmatching import pointmatching
from timagetk.algorithms.resample import resample
from timagetk.algorithms.trsf import inv_trsf
from timagetk.algorithms.trsf import create_trsf
from timagetk.io import imread
from timagetk.io.image import _image_from_url
from timagetk.visu.stack import stack_browser
from timagetk.visu.stack import orthogonal_view
from timagetk.visu.mplt import grayscale_imshow
/builds/J-gEBwyb/0/mosaic/timagetk/src/timagetk/components/labelled_image.py:31: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from tqdm.autonotebook import tqdm

Get the shared multi-angle intensity images

To access the shared multi-angle intensity images, we use 3 view angle images from the first time-point of the p58 time-serie.

file_urls = [
    "https://zenodo.org/record/7151866/files/090223-p58-flo-top.lsm",
    "https://zenodo.org/record/7151866/files/090223-p58-flo-tilt1.lsm",
    "https://zenodo.org/record/7151866/files/090223-p58-flo-tilt2.lsm"
]
file_hashs = [
    "5548917d3d1490200d0d56cbb31c0d35",
    "fa004a84100fcd5ff1d2f9d7c5cbda30",
    "006bf3bacb667eb72d0fb9a6419edcc1"
]
fnames = [_image_from_url(url ,hash_value=md5_hash, hash_method='md5') for url, md5_hash in zip(file_urls, file_hashs)]
  0%|          | 0.00/14.7M [00:00<?, ?B/s]
  7%|▋         | 1.06M/14.7M [00:00<00:01, 11.1MB/s]
 14%|█▍        | 2.12M/14.7M [00:00<00:01, 10.3MB/s]
 21%|██▏       | 3.12M/14.7M [00:00<00:01, 9.82MB/s]
 29%|██▊       | 4.22M/14.7M [00:00<00:01, 10.3MB/s]
 36%|███▌      | 5.25M/14.7M [00:00<00:00, 10.4MB/s]
 44%|████▎     | 6.41M/14.7M [00:00<00:00, 10.9MB/s]
 51%|█████▏    | 7.56M/14.7M [00:00<00:00, 11.1MB/s]
 59%|█████▉    | 8.66M/14.7M [00:00<00:00, 11.1MB/s]
 67%|██████▋   | 9.84M/14.7M [00:00<00:00, 11.3MB/s]
 75%|███████▍  | 11.0M/14.7M [00:01<00:00, 11.4MB/s]
 83%|████████▎ | 12.2M/14.7M [00:01<00:00, 11.6MB/s]
 91%|█████████ | 13.3M/14.7M [00:01<00:00, 11.5MB/s]
 98%|█████████▊| 14.4M/14.7M [00:01<00:00, 11.5MB/s]
14.7MB [00:01, 11.2MB/s]                            

  0%|          | 0.00/12.5M [00:00<?, ?B/s]
  7%|▋         | 864k/12.5M [00:00<00:01, 8.83MB/s]
 16%|█▌        | 1.94M/12.5M [00:00<00:01, 10.1MB/s]
 25%|██▍       | 3.09M/12.5M [00:00<00:00, 10.9MB/s]
 33%|███▎      | 4.16M/12.5M [00:00<00:00, 10.8MB/s]
 42%|████▏     | 5.19M/12.5M [00:00<00:00, 10.7MB/s]
 50%|████▉     | 6.22M/12.5M [00:00<00:00, 10.7MB/s]
 59%|█████▊    | 7.31M/12.5M [00:00<00:00, 10.8MB/s]
 67%|██████▋   | 8.38M/12.5M [00:00<00:00, 10.7MB/s]
 76%|███████▌  | 9.41M/12.5M [00:00<00:00, 10.3MB/s]
 84%|████████▎ | 10.4M/12.5M [00:01<00:00, 10.3MB/s]
 92%|█████████▏| 11.4M/12.5M [00:01<00:00, 10.3MB/s]
12.5MB [00:01, 10.5MB/s]                            

  0%|          | 0.00/13.5M [00:00<?, ?B/s]
  7%|▋         | 1.00M/13.5M [00:00<00:01, 10.3MB/s]
 16%|█▌        | 2.09M/13.5M [00:00<00:01, 10.8MB/s]
 24%|██▍       | 3.22M/13.5M [00:00<00:00, 10.9MB/s]
 32%|███▏      | 4.28M/13.5M [00:00<00:00, 10.8MB/s]
 40%|████      | 5.41M/13.5M [00:00<00:00, 11.0MB/s]
 48%|████▊     | 6.47M/13.5M [00:00<00:00, 10.9MB/s]
 56%|█████▌    | 7.53M/13.5M [00:00<00:00, 10.8MB/s]
 64%|██████▎   | 8.56M/13.5M [00:00<00:00, 10.8MB/s]
 72%|███████▏  | 9.62M/13.5M [00:00<00:00, 10.7MB/s]
 80%|███████▉  | 10.8M/13.5M [00:01<00:00, 10.8MB/s]
 88%|████████▊ | 11.8M/13.5M [00:01<00:00, 10.8MB/s]
 96%|█████████▌| 12.9M/13.5M [00:01<00:00, 11.0MB/s]
13.5MB [00:01, 10.9MB/s]                            

Load the multi-angle images:

img_list = [imread(fname) for fname in fnames]

We use the first image of the list, which is a ‘top’ view, as the initial reference image.

ref_img = img_list[0]
float_imgs = img_list[1:]

Down-sample the images to speed-up the computations

print("# - Voxelsize (ZYX) for the initial image list:")
print("\n".join([f"Image '{im.filename}' (#{n}): {im.voxelsize}" for n, im in enumerate(img_list)]))
print("# - Shape (ZYX) for the initial image list:")
print("\n".join([f"Image '{im.filename}' (#{n}): {im.shape}" for n, im in enumerate(img_list)]))
# - Voxelsize (ZYX) for the initial image list:
Image '090223-p58-flo-top.lsm' (#0): [1.0, 0.20031953747828882, 0.20031953747828882]
Image '090223-p58-flo-tilt1.lsm' (#1): [1.0, 0.20031953747828882, 0.20031953747828882]
Image '090223-p58-flo-tilt2.lsm' (#2): [1.0, 0.20031953747828882, 0.20031953747828882]
# - Shape (ZYX) for the initial image list:
Image '090223-p58-flo-top.lsm' (#0): (59, 460, 460)
Image '090223-p58-flo-tilt1.lsm' (#1): (50, 460, 460)
Image '090223-p58-flo-tilt2.lsm' (#2): (54, 460, 460)
img_list = [resample(im, voxelsize=[1., 0.5, 0.5]) for im in img_list]
INFO     [timagetk.algorithms.resample] Resampling using provided voxel-sizes: [1.0, 0.5, 0.5]
INFO     [timagetk.algorithms.resample] Resampling using provided voxel-sizes: [1.0, 0.5, 0.5]
INFO     [timagetk.algorithms.resample] Resampling using provided voxel-sizes: [1.0, 0.5, 0.5]
print("# - Voxelsize (ZYX) for the resampled image list:")
print("\n".join([f"Image '{im.filename}' (#{n}): {im.voxelsize}" for n, im in enumerate(img_list)]))
print("# - Shape (ZYX) for the resampled image list:")
print("\n".join([f"Image '{im.filename}' (#{n}): {im.shape}" for n, im in enumerate(img_list)]))
# - Voxelsize (ZYX) for the resampled image list:
Image '090223-p58-flo-top.lsm' (#0): [1.0, 0.5, 0.5]
Image '090223-p58-flo-tilt1.lsm' (#1): [1.0, 0.5, 0.5]
Image '090223-p58-flo-tilt2.lsm' (#2): [1.0, 0.5, 0.5]
# - Shape (ZYX) for the resampled image list:
Image '090223-p58-flo-top.lsm' (#0): (59, 184, 184)
Image '090223-p58-flo-tilt1.lsm' (#1): (50, 184, 184)
Image '090223-p58-flo-tilt2.lsm' (#2): (54, 184, 184)

Intensity projections for multi-angle images

f = grayscale_imshow(img_list, suptitle="Multi-angle images - Contour projections", title=[img.filename for img in img_list], val_range=[0, 255], threshold=40)
  0%|          | 0/3 [00:00<?, ?image/s]
 33%|███▎      | 1/3 [00:03<00:07,  3.82s/image]
WARNING  [guess_intensity_threshold] Could not detect intensity threshold with 'bimodal search'...
INFO     [guess_intensity_threshold] Using value obtained from 'percentile search'...
 67%|██████▋   | 2/3 [00:07<00:03,  4.00s/image]
100%|██████████| 3/3 [00:12<00:00,  4.07s/image]
100%|██████████| 3/3 [00:12<00:00,  4.04s/image]

../_images/5544ed89c8c610d2757689e6665f1e90b7634721150bb15c6964e5cc7985548f.png

Load the manual registration landmarks

As the angles between the “views” is important, we start with a manual initialisation using landmarks.

flo_pts_01 = np.loadtxt('p58_t0_floating_ldmk-01.txt')
flo_pts_02 = np.loadtxt('p58_t0_floating_ldmk-02.txt')
ref_pts_01 = np.loadtxt('p58_t0_reference_ldmk-01.txt')
ref_pts_02 = np.loadtxt('p58_t0_reference_ldmk-02.txt')

Create manual initial transformations with pointmatching algorithm

man_trsf_01 = pointmatching(flo_pts_01, ref_pts_01, template_img=ref_img, method='rigid')
man_trsf_02 = pointmatching(flo_pts_02, ref_pts_02, template_img=ref_img, method='rigid')
init_trsfs = [man_trsf_01, man_trsf_02]

Iterative fusion

Fast iterative fusion

In this version of iterative fusion, we performs a vectorfield registration ONLY for the last iteration step to save computation & time.

vf_fused_img = iterative_fusion(img_list, init_trsfs=init_trsfs, vectorfield_at_last=True)
INFO     [timagetk.tasks.fusion] Initial rigid registration on the first image of the list...
INFO     [timagetk.algorithms.blockmatching] RIGID registration with a 'left' initialisation matrix
INFO     [timagetk.algorithms.blockmatching] RIGID registration with a 'left' initialisation matrix
INFO     [timagetk.tasks.fusion] Fusion - Initial registration (iteration 0)...
INFO     [timagetk.tasks.fusion] Reference image shape: [184, 184]
INFO     [timagetk.tasks.fusion] Reference image voxel-sizes: [0.5, 0.5]
INFO     [timagetk.tasks.fusion] Floating image #0 shape: [184, 184]
INFO     [timagetk.tasks.fusion] Floating image #0 voxel-sizes: [0.5, 0.5]
INFO     [timagetk.tasks.fusion] Floating image #1 shape: [184, 184]
INFO     [timagetk.tasks.fusion] Floating image #1 voxel-sizes: [0.5, 0.5]
INFO     [timagetk.tasks.fusion] Template image shape: [92, 92]
INFO     [timagetk.tasks.fusion] Template image voxel-sizes: [1.0, 1.0]
INFO     [timagetk.algorithms.blockmatching] AFFINE registration with a 'left' initialisation matrix
INFO     [timagetk.algorithms.blockmatching] AFFINE registration with a 'left' initialisation matrix
WARNING  [timagetk.third_party.vt_parser] Input `image` is a SpatialImage, option should be in ['gray', 'grey', 'linear', 'cspline'] but got 'nearest'!
WARNING  [timagetk.third_party.vt_parser] Input `image` is a SpatialImage, option should be in ['gray', 'grey', 'linear', 'cspline'] but got 'nearest'!
WARNING  [timagetk.third_party.vt_parser] Input `image` is a SpatialImage, option should be in ['gray', 'grey', 'linear', 'cspline'] but got 'nearest'!
INFO     [timagetk.tasks.fusion] Average image shape: [92, 92]
INFO     [timagetk.tasks.fusion] Average image voxel-sizes: [1.0, 1.0]
INFO     [timagetk.tasks.fusion] Fused image shape: [92, 92]
INFO     [timagetk.tasks.fusion] Fused image voxel-sizes: [1.0, 1.0]
INFO     [timagetk.tasks.fusion] Multi-angle fusion done in 8.773s
INFO     [timagetk.tasks.fusion] Fusion - Registration iteration 1...
INFO     [timagetk.tasks.fusion] Reference image shape: [92, 92]
INFO     [timagetk.tasks.fusion] Reference image voxel-sizes: [1.0, 1.0]
INFO     [timagetk.tasks.fusion] Floating image #0 shape: [184, 184]
INFO     [timagetk.tasks.fusion] Floating image #0 voxel-sizes: [0.5, 0.5]
INFO     [timagetk.tasks.fusion] Floating image #1 shape: [184, 184]
INFO     [timagetk.tasks.fusion] Floating image #1 voxel-sizes: [0.5, 0.5]
INFO     [timagetk.tasks.fusion] Floating image #2 shape: [184, 184]
INFO     [timagetk.tasks.fusion] Floating image #2 voxel-sizes: [0.5, 0.5]
INFO     [timagetk.tasks.fusion] Template image shape: [92, 92]
INFO     [timagetk.tasks.fusion] Template image voxel-sizes: [1.0, 1.0]
INFO     [timagetk.algorithms.blockmatching] AFFINE registration with a 'left' initialisation matrix
INFO     [timagetk.algorithms.blockmatching] AFFINE registration with a 'left' initialisation matrix
INFO     [timagetk.algorithms.blockmatching] AFFINE registration with a 'left' initialisation matrix
WARNING  [timagetk.third_party.vt_parser] Input `image` is a SpatialImage, option should be in ['gray', 'grey', 'linear', 'cspline'] but got 'nearest'!
WARNING  [timagetk.third_party.vt_parser] Input `image` is a SpatialImage, option should be in ['gray', 'grey', 'linear', 'cspline'] but got 'nearest'!
WARNING  [timagetk.third_party.vt_parser] Input `image` is a SpatialImage, option should be in ['gray', 'grey', 'linear', 'cspline'] but got 'nearest'!
INFO     [timagetk.tasks.fusion] Average image shape: [92, 92]
INFO     [timagetk.tasks.fusion] Average image voxel-sizes: [1.0, 1.0]
INFO     [timagetk.tasks.fusion] Fused image shape: [92, 92]
INFO     [timagetk.tasks.fusion] Fused image voxel-sizes: [1.0, 1.0]
INFO     [timagetk.tasks.fusion] Multi-angle fusion done in 8.656s
INFO     [timagetk.tasks.fusion] Fusion - Registration iteration 2...
INFO     [timagetk.tasks.fusion] Reference image shape: [92, 92]
INFO     [timagetk.tasks.fusion] Reference image voxel-sizes: [1.0, 1.0]
INFO     [timagetk.tasks.fusion] Floating image #0 shape: [184, 184]
INFO     [timagetk.tasks.fusion] Floating image #0 voxel-sizes: [0.5, 0.5]
INFO     [timagetk.tasks.fusion] Floating image #1 shape: [184, 184]
INFO     [timagetk.tasks.fusion] Floating image #1 voxel-sizes: [0.5, 0.5]
INFO     [timagetk.tasks.fusion] Floating image #2 shape: [184, 184]
INFO     [timagetk.tasks.fusion] Floating image #2 voxel-sizes: [0.5, 0.5]
INFO     [timagetk.tasks.fusion] Template image shape: [184, 184]
INFO     [timagetk.tasks.fusion] Template image voxel-sizes: [0.5, 0.5]
INFO     [timagetk.algorithms.blockmatching] AFFINE registration with a 'left' initialisation matrix
INFO     [timagetk.algorithms.resample] Resampling using provided voxel-sizes: [0.5, 0.5]
INFO     [timagetk.algorithms.blockmatching] VECTORFIELD registration with a 'left' initialisation matrix
INFO     [timagetk.algorithms.blockmatching] AFFINE registration with a 'left' initialisation matrix
INFO     [timagetk.algorithms.resample] Resampling using provided voxel-sizes: [0.5, 0.5]
INFO     [timagetk.algorithms.blockmatching] VECTORFIELD registration with a 'left' initialisation matrix
INFO     [timagetk.algorithms.blockmatching] AFFINE registration with a 'left' initialisation matrix
INFO     [timagetk.algorithms.resample] Resampling using provided voxel-sizes: [0.5, 0.5]
INFO     [timagetk.algorithms.blockmatching] VECTORFIELD registration with a 'left' initialisation matrix
WARNING  [timagetk.third_party.vt_parser] Input `image` is a SpatialImage, option should be in ['gray', 'grey', 'linear', 'cspline'] but got 'nearest'!
WARNING  [timagetk.third_party.vt_parser] Input `image` is a SpatialImage, option should be in ['gray', 'grey', 'linear', 'cspline'] but got 'nearest'!
WARNING  [timagetk.third_party.vt_parser] Input `image` is a SpatialImage, option should be in ['gray', 'grey', 'linear', 'cspline'] but got 'nearest'!
INFO     [timagetk.tasks.fusion] Average image shape: [184, 184]
INFO     [timagetk.tasks.fusion] Average image voxel-sizes: [0.5, 0.5]
BAL_Inverse2DVectorField: divergence: 2958/33856, non-convergence: 247/33856
INFO     [timagetk.tasks.fusion] Fused image shape: [184, 184]
INFO     [timagetk.tasks.fusion] Fused image voxel-sizes: [0.5, 0.5]
INFO     [timagetk.tasks.fusion] Multi-angle fusion done in 12.313s
print(f"Fused image voxel-size: {vf_fused_img.get_voxelsize()}")
print(f"Fused image shape: {vf_fused_img.get_shape()}")
Fused image voxel-size: [0.5, 0.5]
Fused image shape: [184, 184]
orthogonal_view(vf_fused_img, figsize=(15, 15), suptitle="p58-t0 multi-angle iterative fusion in super-resolution with manual initialization")
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[15], line 1
----> 1 orthogonal_view(vf_fused_img, figsize=(15, 15), suptitle="p58-t0 multi-angle iterative fusion in super-resolution with manual initialization")

File /builds/J-gEBwyb/0/mosaic/timagetk/src/timagetk/visu/stack.py:723, in orthogonal_view(image, x_slice, y_slice, z_slice, suptitle, figname, cmap, **kwargs)
    720 if isinstance(image, MultiChannelImage):
    721     image = BlendImage([image[ch] for ch in image.get_channel_names()])
--> 723 x_sh, y_sh, z_sh = image.get_shape()[::-1]
    724 x_ext, y_ext, z_ext = image.get_extent()[::-1]
    725 x_vxs, y_vxs, z_vxs = image.get_voxelsize()[::-1]

ValueError: not enough values to unpack (expected 3, got 2)