Manual landmarks definition using intensity projections and altitude maps¶

%matplotlib inline
import matplotlib.pyplot as plt

import numpy as np
from timagetk.io import imread
from timagetk.io.image import _image_from_url
from timagetk.algorithms.trsf import apply_trsf
from timagetk.algorithms.trsf import inv_trsf
from timagetk.algorithms.pointmatching import apply_trsf_to_points
from timagetk.algorithms.pointmatching import pointmatching
from timagetk.algorithms.trsf import create_trsf
from timagetk.components.spatial_image import SpatialImage
from timagetk.gui.mpl_image_surface_landmarks import SurfaceLandmarksPlot

/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

INFO     [timagetk.visu.util] Custom colormap 'glasbey' is already registered, use `force=True` to register again!

---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
Cell In[2], line 10
      8 from timagetk.algorithms.trsf import create_trsf
      9 from timagetk.components.spatial_image import SpatialImage
---> 10 from timagetk.gui.mpl_image_surface_landmarks import SurfaceLandmarksPlot

ImportError: cannot import name 'SurfaceLandmarksPlot' from 'timagetk.gui.mpl_image_surface_landmarks' (/builds/J-gEBwyb/0/mosaic/timagetk/src/timagetk/gui/mpl_image_surface_landmarks.py)

import plotly.graph_objects as go
from plotly.subplots import make_subplots

def coordinates_tables(tables, names=None, title=None):
    n_tables = len(tables)
    fig = make_subplots(rows=1, cols=n_tables, shared_yaxes=True, horizontal_spacing=0.03, specs=[[{"type": "table"}]*n_tables]    )
    if names is not None:
        try:
            assert len(names) == n_tables
        except AssertionError:
            raise ValueError(f"Not the same number of tables ({n_tables}) and names ({len(names)})!")
        
    for i, pts in enumerate(tables):
        table = go.Table(header={'values': ["X", "Y", "Z"], 'font': {'size': 10}, 'align': "center", 'prefix': names[i]}, 
                         cells={'values': np.round(pts, 3).T, 'align': "left"})
        fig.add_trace(table, row=1, col=i+1)

    fig.update_layout(height=650, showlegend=False, title_text=title)
    fig.show()

import plotly.graph_objects as go

def alti2surface(alti_map, proj=None, z_invert=True):
    import plotly.graph_objects as go
    # Get altitude map and create plotly surface object
    if z_invert:
        alti_map = -alti_map  # revert the map to see it from top
    sh_0, sh_1 = alti_map.shape
    if proj is not None:
        ye, xe = proj.get_extent()
        x, y = np.linspace(0, xe, sh_0), np.linspace(0, ye, sh_1)
        surf_fig = go.Surface(z=alti_map, x=x, y=y, surfacecolor=proj,
                              colorbar={'title': 'Intensity'},
                              colorscale='Greys', reversescale=True)
    else:
        x, y = np.linspace(0, sh_0, sh_0), np.linspace(0, sh_1, sh_1)
        surf_fig = go.Surface(z=alti_map, x=x, y=y,
                              colorbar={'title': 'Altitude'},
                              colorscale='Viridis')
    return surf_fig

def pts2scatter(pts, z_invert=True):
    px, py, pz = np.array(pts).T
    if z_invert:
        pz = -pz
    pts_fig = go.Scatter3d(x=px, y=py, z=pz, mode='markers', marker_symbol='diamond',
                           marker={'size':6, "line": {'width': 2, 'color':'DarkSlateGrey'}})
    return pts_fig

Use known transformation from `create_trsf`¶

To test the GUI & the algorithms, we use a known rigid transformation and:

apply it to the reference image to create the floating image
apply the invert of the transformation to the reference points to create floating points as pairs of landmarks

Load the reference image¶

ref_img = imread(_image_from_url("https://zenodo.org/record/7151866/files/090223-p58-flo-top.lsm",
                                 hash_value="5548917d3d1490200d0d56cbb31c0d35", hash_method='md5'))
ref_pts = np.loadtxt('p58_t0_reference_ldmk-01.txt')

Create the rigid transformation, floating image & points¶

trsf = create_trsf('random', template_img=ref_img, trsf_type='rigid', seed=12, angle_range=[0.0, 0.1])
flo_img = apply_trsf(ref_img, trsf, template_img=ref_img)
flo_img.filename = 'registered_'+flo_img.filename
flo_pts = apply_trsf_to_points(ref_pts, inv_trsf(trsf, template_img=ref_img))

# print(ref_pts)
# print(flo_pts)
flo_img.filename

Surface landmarks GUI:¶

gui = SurfaceLandmarksPlot(img_left=ref_img, img_right=flo_img, surface_threshold=np.percentile(ref_img, 60),
                           points={'left': ref_pts, 'right': flo_pts}, orientation={'left': 1, 'right': 1},
                           altimap_smoothing=True)

Access the defined 3D landmarks¶

From the 2D points, in real coordinates, we return the position in 3D using the corresponding altitude map:

pred_ref_pts = gui.points_3d_coordinates()['left']
pred_flo_pts = gui.points_3d_coordinates()['right']

Access the estimated transformations¶

Using these 3D points (except the last one from each list), the GUI compute a landmarks based transformation using the pointmatching algorithm. The green channel is the floating (right) image registered by this transformation.

Usually, we you place landmarks, it is to initialize a registration process, often starting with a rigid block-matching registration. To mimmic this behaviour and anticipate the outcome of the next step, the GUI also estimate a rigid transformation with the landmarks transformation as initialisation. The blue channel is the floating (right) image registered by the composotion of manual & rigid transformations.

gui._update_trsf()  # Make sure the transformations are "up-to-date"

# The landmark based rigid transformation
np.round(gui.trsf.get_array(), 3)

# The `blockmatching` rigid transformation initialized with the landmark transformation
np.round(gui.rigid_trsf.get_array(), 3)

Access the projection images & altitudes maps¶

ref_alti = gui.altitude_maps['left']
ref_proj = gui.projected_images['left']

flo_alti = gui.altitude_maps['right']
flo_proj = gui.projected_images['right']

Show associated altitude maps¶

import matplotlib.gridspec as gridspec
from matplotlib import colors

fig = plt.figure(figsize=(15, 10), constrained_layout=True)
# Create two panels
gs = gridspec.GridSpec(ncols=2, nrows=1, figure=fig)
ref_ax = fig.add_subplot(gs[0, 0])
flo_ax = fig.add_subplot(gs[0, 1])
# Add the altitude map as images:
images = []
images.append(ref_ax.imshow(ref_alti, cmap='viridis', extent=gui._mpl_panel_extent('left')))
images.append(flo_ax.imshow(flo_alti, cmap='viridis', extent=gui._mpl_panel_extent('right')))
ref_ax.xaxis.tick_top()
flo_ax.xaxis.tick_top()
# Add the landmarks points:
px, py, c = np.array([(p[0], p[1], n) for n, p in enumerate(pred_ref_pts)]).T
ref_ax.scatter(px, py, c=c, cmap='Set3', ec='k', s=100)
px, py, c = np.array([(p[0], p[1], n) for n, p in enumerate(pred_flo_pts)]).T
flo_ax.scatter(px, py, c=c, cmap='Set3', ec='k', s=100)

# Find the min and max of all colors for use in setting the color scale.
vmin = min(image.get_array().min() for image in images)
vmax = max(image.get_array().max() for image in images)
norm = colors.Normalize(vmin=vmin, vmax=vmax)
for im in images:
    im.set_norm(norm)
# Add a colorbar at the bottom
fig.colorbar(images[0], ax=[ref_ax, flo_ax], orientation='horizontal', fraction=.1)

Estimated transformations from landmarks¶

Ground truth transformation¶

Using the generated set of landmarks from a random rigid transformation, we can get the “ground truth” transformation:

# To compute a transformation to apply to the floating image using landmarks:
gt_trsf = pointmatching(flo_pts, ref_pts, template_img=ref_img, method='rigid', real=True)
np.round(gt_trsf.get_array(), 3)

# The initial 'random' rigid trsf:
np.round(trsf.get_array(), 3)

# The GT transfo is the invert of the initial 'random' rigid trsf:
np.round(inv_trsf(trsf).get_array(), 3)

Predicted transformation from landmarks point-matching¶

#pred_ref_pts = gui.points_3d_coordinates()['left']
#pred_flo_pts = gui.points_3d_coordinates()['right']
pred_trsf = pointmatching(pred_flo_pts, pred_ref_pts, template_img=ref_img, method='rigid', real=True)
np.round(pred_trsf.get_array(), 3)

Note that it may be slightly different from the gui.trsf as we use all the 3D points in the previous cell, whereas the GUI do not use the last one (as you may be in the process of editing it).

# Deviation from "ground truth" (GT - PRED)
np.round(gt_trsf.get_array() - pred_trsf.get_array(), 3)

Estimate the deviation from ‘ground truth’ landmarks coordinates:¶

Knowing the true floating point coordiantes, we can estimate the loss of precision induced by the use of projection and altitude map in the GUI. Only the Z coordinates (last column) will differs as XY remains unchanged by the GUI.

pred_pts = gui.points_3d_coordinates()['right']
# print(pred_pts)
diff = pred_pts - flo_pts

coordinates_tables([flo_pts, pred_pts, diff], ['GT_', 'Pred_', 'Diff_'], "Estimation bias on floating points by GUI")

3D scatter plot of ground truth & gui floating image landmarks¶

import plotly.express as px
import pandas as pd

pred_pts = gui.points_3d_coordinates()['right']
group_names = ['Ground Truth', 'Prediction']
pts_list = []
for n, pts in enumerate([flo_pts, pred_pts]):
    pts_list.extend([list(p)+[group_names[n]] for p in pts])

df = pd.DataFrame(pts_list, columns=['x', 'y', 'z', 'group'])

fig = px.scatter_3d(df, x='x', y='y', z='z', color='group')
fig.show()

3D rendering of the altitude map & projection¶

from plotly.subplots import make_subplots

fig = make_subplots(rows=1, cols=2, specs=[[{'type': 'surface'}, {'type': 'surface'}]])
fig.add_trace(alti2surface(ref_alti, ref_proj), col=1, row=1)
fig.add_trace(pts2scatter(ref_pts), col=1, row=1)
fig.add_trace(alti2surface(flo_alti, flo_proj), col=2, row=1)
fig.add_trace(pts2scatter(pred_pts), col=2, row=1)

Table of Contents