Th3 nub3rs 0f a br41n
Mathematics behind popular algorithms in today's Neuromaging
Davide Poggiali
PhD Student, Neurosciences
Roma, Istituto per le Applicazioni del Calcolo M. Picone
April 14, 2016
The most popular commercial, free and FOSS software:¶
- Pmod
(PET, ~2700 euro) - Xinapse-Jim7
(MRI/gen. imaging, ~3000 euro) - freesurfer
(MRI, free) - SPM
( (f-)MRI, FOSS but MATLAB-based) - fsl
(DTI/(f-)MRI, free) - mipav
(MRI, free) - ANTs
(MRI, FOSS)
What is a (medical) image?¶
It is a function, say $$I: \Omega \subset \mathbb{R}^3 \longrightarrow \mathbb{R}$$ sampled over an equispaced grid of points $X=\left\{ x_{1},\dots x_{n} \right\}$ s.t. $x_i \in \Omega$
A segmentation or VOI is an arbitrary subset $\Gamma \subseteq \Omega$ and its mask a function
$$\gamma(x) = \left\{ \begin{array}{c c} 1 & x\in \Gamma \\ 0 & x\notin \Gamma \end{array} \right.$$# let's have a look at the histogram!
I=nib.load("t1w.nii.gz").get_data()
figure(1)
a=hist(I[(I>10)&(I<2500)],100)
f, (ax1, ax2, ax3) = plt.subplots(1, 3,figsize=(18,18))
J1=np.zeros(I.shape); J2=np.zeros(I.shape)
thresh1=(I>1200)&(I<2000)
J1[thresh1]=I[thresh1]
thresh2=(I>500)&(I<1200)
J2[thresh2]=I[thresh2]
ax1.imshow(rot90(I[80,:,:]))
ax2.imshow(rot90(J1[80,:,:]))
ax3.imshow(rot90(J2[80,:,:]))
Or, if more than one sequence is available, multithresholding¶
References¶
- F Kurugollua, B Sankurc, A.E Harmancidl, Color image segmentation using histogram multithresholding and fusion, Image and Vision Computing, 2001
- Z Zeng, R Zwiggelaar, Joint Histogram Modelling for Segmentation Multiple Sclerosis Lesions, Springer Berlin Heidelberg, 2011
- W Burger, MJ Burge, Principles of Digital Image Processing: Fundamental Techniques, Undergraduate Topics in Computer Science, 2009
- D Poggiali, A Favaretto, M Puthenparampil, F Causin, P Gallo, Segmentation of lesions in Multiple Sclerosis: a multithresholding approach, CAEconf 2014
...and after!¶
Why do we need this:¶
- align two images
- to compare them
- to make measurement on the same space
- for motion correction
- to compare images and find changes
- align $n$ images to a reference/template
- to see them in the same orientation
- to make measurement in standard VOIs
- to find longitudinal changes
- ....
There are many ways to align an input volume to a reference¶
in FSL¶
/usr/local/fsl/bin/flirt -in input.nii.gz -ref ref.nii.gz -out output.nii.gz -omat output.mat -bins 256 -cost corratio -searchrx -90 90 -searchry -90 90 -searchrz -90 90 -dof 12 -interp trilinear
in ANTs¶
bash $ANTSPATH/antsRegistrationSyN.sh -d 3 -f reference.nii.gz -m input.nii.gz -t a -o output
This saves two files: the ouput.nii.gz volume and the transformation matrix output.mat¶
You can also apply the same transformation to a second volume, say input2.nii.gz¶
in FSL¶
/usr/local/fsl/bin/flirt -applyxfm -in input2.nii.gz -ref ref.nii.gz -out output2.nii.gz -init output.mat -interp trilinear
in ANTs¶
$ANTSPATH/antsApplyTransforms -d 3 -i input2.nii.gz -r ref.nii.gz -o output_lesions.nii.gz -n NearestNeighbor -t output.mat
What are we doing?¶
Say that the image I is represented by a function $$I: \Omega \subset \mathbb{R}^3 \longrightarrow \mathbb{R}$$ and the reference image is $R$ we look for a function $$f: \mathbb{R}^3 \longrightarrow \mathbb{R}^3 $$
Such that $$\arg\min_{f} \mathcal{D} \left[ I(f(x)), R(x) \right]$$ where $\mathcal{D}\left[ \cdot , \cdot \right]$ is an appropriate distance (or pseudo-distance).
Different kind of transforms¶
Let us consider linear tranformation $$ f(x) = A\cdot x +b$$
DOF:
- 3 we look just for the translation $b$
- 6 we look for $b$ and a 3D rotation $R$
- 7 we look for $b$ and a 3D rotation multiplied by a constant $A=a R$
- 9 we look for $b$ and a similitude $$A=R\cdot\left( \begin{array}{c c c} a & 0 & 0\\ 0 & b & 0 \\ 0 & 0 & c \end{array} \right)$$
- 12 we look for $b$ and any $A$ such that $det(A)\neq 0$
The role of the distance function¶
The distance function $\mathcal{D}\left[ \cdot , \cdot \right]$ plays an important role in registration procedure
Some examples:
SSD, sum squared differences $$SSD \left[ A , B \right]=1/2 \int_{\Omega}\left(A(x)-B(x)\right)^2 \; dx$$
NCC, normalized cross-correlation $$NCC \left[ A , B \right]=1- \frac{\left< A,B\right>^2}{||A||_2 ||B||_2}$$ with the scalar product $$\left< A,B\right>=\int_{\Omega} A(x)B(x) \; dx$$
NGF, normalized gradient-field
$$NGF \left[ A , B \right]= \int_{\Omega} 1 - \left< g_{\eta}(A),g_{\eta}(B)\right>^2 \; dx$$were the gradient-field is defined as $$g_{\eta}(A(x)) = {\bigtriangledown A(x) \over \sqrt{||\bigtriangledown A(x)||^2 + \eta}}$$
KLD, Kullback_Leibner Divergence
$$KLD \left[ A , B \right]= \int_{\Omega} \ln \left( \frac{A(x)}{B(x)} \right)A(x) \; dx$$
References¶
J Modersitzki, Numerical Methods for Image Registration (2004) Oxford University Press
Jan Modersitzki, FAIR: Flexible Algorithms for Image Registration (2009) SIAM
G Dougherty, Digital Image Processing for Medical Applications (2009) Cambridge University Press
JV Hajnal, D Hill and DJ Hawkes, Medical Image Registration (2001) CRC Press
A Passarini, Medical Image Registration for Motion Detection and Correction, thesis, 2015
Let's use DiReCT!¶
Diffeomorphic Registration based Cortical Thickness¶
A brain image $$P(x): \Omega \subseteq \mathbb{R}^3 \longrightarrow \mathbb{R}$$ is segmented in three binary masks $$P_g+P_w+P_c=1$$ representing GM, WM, CSF respectively.
Let us define the gray-white matter interface $$GWI=P_g \cup dil_1(P_w)$$ and the brain tissue mask as $$P_{wg}=P_w \cup P_g$$
We look for the diffeomorphism $$\phi : \Omega \times [0,1] \longrightarrow \Omega $$ such that:
- $\phi(x,0)=x$
- the thickness is contrained, $\left|\phi(x,1)-\phi(x,0)\right|<\tau$ for each $x\in GWI$
- there exists a velocity field $\Omega \times [0,1]\longrightarrow \mathbb{R}^3$ such that $$\frac{d\phi}{dt}=v(\phi(x,t),t)$$ or equivalently $$\phi(x,t)=x+\int_{0}^{t} v(\phi(x,t),t) \; dt$$
So we look for a diffeomorphism $\phi$ which minimizes $$E(\phi(x,1)) = \int_{0}^{1} ||v(x,t)||^2 \; dt \;\; + \;\; \left|\left|P_w(\phi(x,1)) - P_{wg} (x) \right| \right|^2$$ subject to: $$\phi(x,t)=x+\int_{0}^{t} v(\phi(x,t),t) \; dt$$ $$T(x)=\left|\phi(x,1)-\phi(x,0)\right|<\tau \;\; \forall x\in GWI$$
Example¶
References¶
- SR Das, BB Avants, M Grossman, JC Gee, Registration based cortical thickness measurement Neuroimage. 2009 Apr 15;45(3):867-79. doi: 10.1016/j.neuroimage.2008.12.016
Tracer kinetic in brain PET¶
A simple example¶
Let us suppose to have a time-series PET data obtained using $^{18}$F-dg tracer $$P(x) : \Omega \times [0,T_{e}] \longrightarrow \mathbb{R}$$ Sampled over a discerete set of points $X$ and times $T$. If we want to estimate the absolute rate of absorption MRGlu, we need a model of the tracer kinetic.
The Sokoloff model¶
Were C(t) is the mean of the activity inside a given VOI.
This is a ODE compartment model $$\left\lbrace \begin{array}{c c c} \dot{C_p} & = &-k_1 C_p + k_2 C_f + u& \\ \dot{C_f} & = &-k_2 C_f -k_3 C_f& \\ \dot{C_b} & = &+k_3 C_f&\\ \end{array}\right. $$ Once fitted, we can compute $K=\frac{k_1 k_3}{k_2 + k_3}$ and find the glucose absorption rate as $$MRGlu = \frac{PG}{LC} K$$
References¶
- DH Anderson, Compartmental Modeling and Tracer Kinetics, Springer Berlin Heidelberg, 1983.
- RE Carson, P Herscovitch, ME Daube-Witherspoon, Quantitative Func- tional Brain Imaging with Positron Emission Tomography, Academic Press, 1998
- C Cobelli, E Carson, Introduction to modeling in physiology and medicine, AP, 2008
- C Cobelli, D Foster, G Toffolo, Tracer kinetics in biomedical research, Springer, 2001.
- CS Patlak, RG Blasberg, JD Fenstermacher, Graphical evaluation of blood-to-brain transfer constants from multiple-time uptake data. J Cereb Blood Flow Metab 1983, 3(1):1-7