Using TIRL to map between MRI and microscopy data in BigMac

The BigMac MRI and microscopy data is co-registered using the Tensor Image Registration Library (TIRL) in FSL (Huszar2019). This tutorial provides a basic overview of how to map between MRI and microscopy data in BigMac.

More information on TIRL can be found at https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/TIRL and https://git.fmrib.ox.ac.uk/ihuszar/tirl. The git also includes a tutorial. The relevant paper is https://www.biorxiv.org/content/10.1101/849570v1.full.pdf.

In BigMac, TIRL has been run on every microscopy slide in turn, to create a warpfield between the microscopy and structural MRI. This warpfield is stored as a TIRL image or .timg file.

Here the structural MRI is used as reference due to its high spatial resolution and good image contrast (clear gm/wm delineation).

MRI-MRI warpfields generated using FSLs FLIRT/FNIRT can be used to map between the microscopy and other MRI contrasts such as dMRI.

Note, TIRL is all python-based. It is included in recent versions of FSL ann can be run using fslpython. Alternatively tirl can be installed via github above. If doing the latter, please use version 2.1.3b1 or above

The TIRL transformations in each timg describe going from the native/voxel space to the tirl/physical space where the data are co-aligned.

Scripts for TIRL mapping in BigMac can be found at: https://git.fmrib.ox.ac.uk/amyh/bigmacanalysis/-/tree/main/TIRL

# Libraries
import tirl
import numpy as np
from tirl import fsl
from tirl.timage import TImage
import os
import sys
from PIL import Image
import nibabel as nib

Load data and warpfields

# BigMac directory
ddir = "/vols/Data/km/ahoward/BigMac-clean/"

# Folder containing microscopy data and tirl warpfields
base = ddir + "/Microscopy/PLI/Anterior/P089x/"

# timg containing information to transform microscopy to physical space
# (i.e. tirl intemediary space) where the data are coregistered
micro_timg = base + "/Reg2MRI/Full/slice.timg"

# timg containing information to transfrom mri to physical space
mri_timg = base + "/Reg2MRI/Full/volume.timg"

# High resolution PLI images
pli_inplane = base+"/In_plane.tif" # The data we want to warp
pli_inclination = base+"/Inclination.tif" # The image used to drive the registration

# FNIRT transformations, dwi -> struct (opposite to what you think you need)
# b10k 1mm data
dwi = ddir + "/MRI/Postmortem/dwi/b10k/1.0mm/data/S0.nii.gz"
field = ddir + "/MRI/Postmortem/dwi/b10k/1.0mm/reg/dwi_1.0mm_2_struct_fnirt/dwi_1.0mm_2_struct_fnirt_field.nii.gz"

# Load microscopy timg
micro = tirl.load(micro_timg)

Inbuilt functions

TIRL Timages contain lots of useful inbuilt functions. Use autocomplete to see full list.

Timage structure

• micro.domain.chain = transformation chain native => physical space where data is coregistered

• micro.data = access / manipulate data

Useful functions

• micro.evaluate(domain) # => interpolate data at coords of some other domain

• micro.snapshot(fout) # => save snapshot image of data to fout

• micro.domain.chain # => access transformation chain. Output from tirl defined as from native/voxel space -> physical/tirl space

• micro.domain.chain.inverse() # => invert chain

• micro.domain.copy_transformations(domain) # => copy transformations from one timg to another

• micro.domain.get_physical_coordinates() # => get microscopy coordinates in physical space where data is coregistered

• micro.domain.get_voxel_coordinates() # => get microscopy pixels / mri voxels

• micro.domain.chain.map_vectors(vectors,coordinates) # => map vectors (at voxel coordinates) through transformation chain

• micro.domain.chain.map(coordinates) # => map coordinates through transformation chain

• micro.domain.chain.map_vector_and_coordinates(vectors,coordinates) # => map vectors and voxel coordinates

• micro.domain.map_physical_coordinates() # => map coordinates in physical space to voxel/native space

• micro.data.flatten() # => returns flattened data

• micro.shape # => shape of data (after offset transformations)

• micro.header # => header information

• micro.interpolator(coords) # => interpolate data at coordinate positions

# Note we can load any data as a TImage to access these functions
# from tirl.timage import TImage
fa = TImage('/vols/Data/BigMac/MRI/Postmortem/dwi/b10k/1.0mm/data.dtifit/dti_FA.nii.gz')
inplane = TImage(pli_inplane)
inplane.domain.copy_transformations(micro.domain)

# accessing help
help(micro.domain)

Combining warpfields from tirl and flirt/fnirt

example: create warpfield microscopy to dwi

# 1. Load microscopy slide that we have registered to structural MRI space using TIRL
# (microscopy -> physical)
micro = tirl.load(micro_timg)
chain = micro.domain.chain.copy()

# 2. Create warpfield from microscopy to structural
# (microscopy -> physical -> mri)
struct = tirl.load(mri_timg)
# Invert warpfield in mri_timg and add to warpfield from micro_timg
tirl2struct = struct.domain.chain.inverse()
chain += tirl2struct

# 3. Import FNIRT transformation & add to transformation chain
# (microscopy -> physical -> structural -> dwi)
# Note, opposite transform to what you might expect due to how fsl defines warpfields
warp = fsl.load_warp(field, moving=dwi)
chain += warp

Mapping coordinates

a) map pixel coordinates to structural MRI voxel locations

# Map coordinates of high resolution microscopy image through tirl chain to dwi space
micro_hr = TImage(pli_inclination)
micro_hr.domain.copy_transformations(micro.domain)

# Find hr pixel coordinates in physical (tirl) space
pixel_physical_coords = micro_hr.domain.get_physical_coordinates()

# Map to structural space
pixel_coords_in_struct = tirl2struct.map(pixel_physical_coords[0:10,:])
# Here only mapping first 10 pixels as an example
print(pixel_coords_in_struct)

[[229.49050031 107.73321454 218.01824585]
 [229.47817781 107.73373977 218.01531657]
 [229.46585674 107.73426504 218.01238556]
 [229.45353709 107.73479037 218.00945283]
 [229.44121889 107.73531575 218.00651839]
 [229.42890211 107.73584118 218.00358222]
 [229.41658677 107.73636666 218.00064433]
 [229.40427287 107.73689219 217.99770471]
 [229.3919604  107.73741778 217.99476338]
 [229.37964936 107.73794342 217.99182032]]

# Note, this is the same as:
micro = tirl.load(micro_timg)
micro2struct = micro.domain.chain + tirl2struct
pixel_native_coords = micro_hr.domain.get_voxel_coordinates()  # coordinates in native space
pixel_coords_in_struct = micro2struct.map(pixel_native_coords[0:10,:])
print(pixel_coords_in_struct)
# This version is less memory intensive as we are only ever mapping 10 pixels. In comparison,
# the function get_pixel_coordinates() maps all microscopy pixels to real space.

[[229.49050031 107.73321454 218.01824585]
 [229.47817781 107.73373977 218.01531657]
 [229.46585674 107.73426504 218.01238556]
 [229.45353709 107.73479037 218.00945283]
 [229.44121889 107.73531575 218.00651839]
 [229.42890211 107.73584118 218.00358222]
 [229.41658677 107.73636666 218.00064433]
 [229.40427287 107.73689219 217.99770471]
 [229.3919604  107.73741778 217.99476338]
 [229.37964936 107.73794342 217.99182032]]

b) map pixel coordinates to dwi voxels

# Define transformation chain from tirl/physical space to dwi space
tirl2dwi = tirl2struct + warp

# Map to dwi
pixel_coords_in_dwi = tirl2dwi.map(pixel_physical_coords[0:10,:])
print(pixel_coords_in_dwi)

# Again this is equivalent to (where the following is less memory intensive)
# pixel_coords_in_dwi = chain.map(pixel_native_coords[0:10,:])
# where chain maps micro->dwi (see chain construction above)
# print(pixel_coords_in_dwi)

[[68.44825822 31.06628612 66.21822589]
 [68.44448256 31.06629764 66.21740675]
 [68.44070676 31.06630796 66.21658779]
 [68.43693085 31.0663171  66.21576896]
 [68.43315482 31.06632511 66.21495027]
 [68.42937871 31.06633201 66.21413169]
 [68.42560253 31.06633783 66.2133132 ]
 [68.42182628 31.06634262 66.21249479]
 [68.41804999 31.0663464  66.21167644]
 [68.41427366 31.0663492  66.21085813]]

Mapping vectors

TIRL can also map vectors and account for the rotations in the transformation chain.

Mapping the normal from the microscopy plane is slightly more complex as the TIRL transformation chain has some 2D and some 3D transformations.  In general, we must take care to make sure we have a good microscopy -> mri mapping.

a) mapping the orientation angle from PLI through to dwi space

# Map 2D PLI vector (orientation angle) through tirl chain into diffusion space

# Load in-plane angle
inplane = TImage(pli_inplane)
# Change scale 0-pi
inplane.normalise()
inplane = inplane*np.pi

# Convert inplane angle to x/y

x = -np.sin(inplane[:]).ravel()
y = np.cos(inplane[:]).ravel()
xy = np.c_[x, y]

# map pli vector through img.domain.chain -> dwi space
pixel_coords = inplane.domain.get_voxel_coordinates()
inplane_in_dwi = chain.map_vector(xy[0:10,:], pixel_coords[0:10,:], rule="fs")

# flip vectors LR so they are stored in a similar manner to FSL vectors (e.g. the output from dtifit)
inplane_in_dwi = inplane_in_dwi*[-1,1,1]
print(inplane_in_dwi)

[[-0.98966752  0.00262689  0.14335721]
 [-0.98966752  0.00262689  0.14335721]
 [-0.98966752  0.00262689  0.14335721]
 [-0.98966752  0.00262689  0.14335721]
 [-0.98966752  0.00262689  0.14335721]
 [-0.98966752  0.00262689  0.14335721]
 [-0.98966752  0.00262689  0.14335721]
 [-0.98966046  0.00263049  0.14340588]
 [-0.98966046  0.00263049  0.14340588]
 [-0.98966046  0.00263049  0.14340588]]

b) mapping microscopy normal to dwi space

# Need to separate 2D and 3D transformations in tirl
chain2D = chain[:7]
chain3D = chain[7:]

# Define normal after 2D transformations
norm = np.zeros(inplane_in_dwi.shape)
norm[:,2] = 1
pixel_coords_2D = chain2D.map(pixel_coords)
norm_in_dwi = chain3D.map_vector(norm[0:10,:], pixel_coords_2D[0:10,:], rule="fs")

# flip vectors LR so they are stored in a similar manner to FSL vectorss (e.g. the output from dtifit)
norm_in_dwi = norm_in_dwi*[-1,1,1]
print(norm_in_dwi)

[[-0.01226941 -0.99771633 -0.06641971]
 [-0.01226941 -0.99771633 -0.06641971]
 [-0.01226941 -0.99771633 -0.06641971]
 [-0.01226941 -0.99771633 -0.06641971]
 [-0.01226941 -0.99771633 -0.06641971]
 [-0.01226941 -0.99771633 -0.06641971]
 [-0.01226941 -0.99771633 -0.06641971]
 [-0.01227313 -0.99771779 -0.06639714]
 [-0.01227313 -0.99771779 -0.06639714]
 [-0.01227313 -0.99771779 -0.06639714]]

Interpolating / resampling data

e.g. resampling dwi data on microscopy slide

# Load FA image as TImage and compile correct transformation chain (dwi -> real)
fa = TImage(ddir + '/MRI/Postmortem/dwi/b10k/1.0mm/data.dtifit/dti_FA.nii.gz')

# Since mapping to HR microscopy space can be time consuming, can check mapping by instead interpolating
# in the downsampled microscopy space, as defined by micro (instead of micro_hr)
# e.g
micro = TImage(micro_timg)
micro.domain.chain = chain
fa_in_micro = fa.evaluate(micro.domain)
fa_in_micro.snapshot('./tmp_FA.tif',overwrite=True)
micro.snapshot('./tmp_pli.tif',overwrite=True)

# OR
micro = TImage(micro_timg)
fa_in_micro = fa.interpolator(tirl2dwi.map(micro.domain.get_physical_coordinates()))
fa_in_micro = TImage(fa_in_micro.reshape(micro.shape))
fa_in_micro.snapshot('./tmp_FA2.tif',overwrite=True)

# OR can change interpolator
from tirl.interpolators.scipyinterpolator import ScipyNearestNeighbours
fa.interpolator = ScipyNearestNeighbours
micro.domain.chain = chain
fa_in_micro_nn = fa.evaluate(micro.domain)
fa_in_micro_nn.snapshot('./tmp_FA3.tif',overwrite=True)

# OR to map to high resolution microscopy (computationally expensive):

# Load FA image as TImage and compile correct transformation chain (dwi -> real)
fa = TImage('/vols/Data/BigMac/MRI/Postmortem/dwi/b10k/1.0mm/data.dtifit/dti_FA.nii.gz')

# Resample FA at microscopy pixel coordinates
micro_hr.domain.chain = chain # chain defined above & transforms from microscopy -> dwi
fa_in_micro = fa.evaluate(micro_hr.domain)
fa_in_micro.snapshot('./tmp.tif',overwrite=True)

# This is equivalent to
fa_in_micro = fa.interpolator(chain.map(micro_hr.domain.get_voxel_coordinates()))
fa_in_micro = TImage(fa_in_micro.reshape(micro_hr.shape))
fa_in_micro.snapshot('./tmp2.tif')

# OR
fa.domain.chain = tirl2dwi.inverse() # This will take a long time to compute. Better to compile from fnirt/flirt & struct2tirl
micro = TImage(micro_timg)
micro_hr.domain.copy_transformations(micro.domain)
fa_in_micro = fa.evaluate(micro_hr.domain)
fa_in_micro.snapshot('./tmp3.tif',overwrite=True)

# OR
fa_in_micro = fa.interpolator(tirl2dwi.map(micro_hr.domain.get_physical_coordinates()))
fa_in_micro = TImage(fa_in_micro.reshape(micro_hr.shape))
fa_in_micro.snapshot('./tmp4.tif')

Manipulating nifti volumes

e.g. to save some output as a nifti

# Copy header information from dwi image and change resolution
dwi_im = nib.load(dwi)
new_hdr = dwi_im.header
new_hdr.set_zooms([0.004,0.004,0.35])

# Save pixel coordinates as a nifti volume
niftiobj = nib.Nifti1Image(pixel_coords_in_dwi, new_hdr.get_best_affine(), new_hdr)
nib.save(niftiobj, './tmp.nii.gz')