Skip to content

Commit

Permalink
BLD Minor polish
Browse files Browse the repository at this point in the history
  • Loading branch information
@spinicist
spinicist committed Jun 24, 2020
1 parent 78fc967 commit a4263ad
Show file tree
Hide file tree
Showing 6 changed files with 246 additions and 669 deletions.
237 changes: 237 additions & 0 deletions Python/T1_mapping.ipynb
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,237 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# T1-Mapping with DESPOT1 / VFA\n",
"\n",
"This notebook deomstrates how T1 mapping with a basic `nipype` pipeline with `QUIT` programs.\n",
"\n",
"This pipeline (http://github.com/spinicist/nanslice) downloads the BrainWeb brain & B1 phantoms from http://brainweb.bic.mni.mcgill.ca. It then replaces the tissue classification labels with values of Proton Density and T1, simulates an SPGR/FLASH image with some added noise, and finally uses that simulated data to fit for T1 and PD again.\n",
"\n",
"In order to display the images, and to do some "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Imports"
]
},
{
"cell_type": "code",
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"from QUIT.interfaces.relax import DESPOT1, DESPOT1Sim, DESPOT2\n",
"from nanslice import Layer\n",
"import nanslice.jupyter as ns\n",
"import nibabel as nib\n",
"import numpy as np\n",
"import requests\n",
"import gzip\n",
"import os.path"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"## Create Phantom Data\n",
"\n",
"This section downloads the Brainweb phantoms. These are stored MINC format so we load the images using `nibabel` to get the raw data later."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Data Fetching\n",
"\n",
"We now download the Brainweb brain & T1 phantoms, if they are not already present in the current directory. These are stored in MINC format, so we use "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Display the phantom data as a check"
]
},
{
"cell_type": "code",
"metadata": {},
"outputs": [],
"source": [
"if not isfile('classes.mnc'):\n",
" print('Downloading classes')\n",
" params = {'download_for_real':'[Start Download!]',\n",
" 'do_download_alias':'phantom_1.0mm_normal_crisp',\n",
" 'format_value':'minc',\n",
" 'who_name': 'Tobias Wood',\n",
" 'who_institution': 'KCL',\n",
" 'who_email': '[email protected]'}\n",
" response = requests.get(url='http://brainweb.bic.mni.mcgill.ca/cgi/brainweb1', params=params)\n",
" minc_file = open('classes.mnc', 'wb')\n",
" minc_file.write(response.content)\n",
"if not isfile('rf20_C.mnc'):\n",
" print('Downloading B1')\n",
" params = {'download_for_real':'[Start Download!]',\n",
" 'do_download_alias':'rf20_C.mnc.gz',\n",
" 'format_value':'minc',\n",
" 'zip_value':'none',\n",
" 'who_name': 'Tobias Wood',\n",
" 'who_institution': 'KCL',\n",
" 'who_email': '[email protected]'}\n",
" response = requests.get(url='https://brainweb.bic.mni.mcgill.ca/cgi/brainweb1', params=params)\n",
" b1_file = open('rf20_C.mnc', 'wb')\n",
" b1_file.write(response.content)\n",
"classes = nanslice.Layer('classes.mnc')\n",
"b0_minc = nanslice.Layer('rf20_B.mnc')\n",
"b1_minc = nanslice.Layer('rf20_C.mnc')"
]
},
{
"cell_type": "code",
"metadata": {},
"outputs": [],
"source": [
"ns.three_plane(classes)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we define a function that allows us to convert the tissue classes to T1 and PD/M0 values, and use it to create our reference T1, PD and B1 images. These will be saved into the current directory. This function also allows us to subsample the images, which is useful when investigating methods that take longer to fit.\n",
"\n",
"The BrainWeb phantom uses the following mapping:\n",
"\n",
"0. Background\n",
"1. CSF\n",
"2. Grey Matter\n",
"3. White Matter\n",
"4. Fat\n",
"5. Muscle/Skin\n",
"6. Skin\n",
"7. Skull\n",
"8. Glial Matter\n",
"9. Connective\n",
"\n",
"To keep things simple, we only specify values for CSF, GM & WM, and set the other tissue types to zero. The B1 data is used as is from the phantom."
]
},
{
"cell_type": "code",
"metadata": {},
"outputs": [],
"source": [
"def classes_to_params(M0Vals, T1Vals, subsamp=1):\n",
" class_data = classes.image.get_data().astype('int32')\n",
" M0data = np.choose(class_data[::subsamp,::subsamp,::subsamp], M0vals).astype('float32')\n",
" T1data = np.choose(class_data[::subsamp,::subsamp,::subsamp], T1vals).astype('float32')\n",
" T2data = np.choose(class_data[::subsamp,::subsamp,::subsamp], T2vals).astype('float32')\n",
" B1data = b1_minc.image.get_data().astype('float32')[::subsamp,::subsamp,::subsamp]\n",
" # PDdata = np.array(list(map(PDFunc, classes.image.get_data())))\n",
" M0image = nib.nifti1.Nifti1Image(M0data, affine=classes.image.affine)\n",
" T1image = nib.nifti1.Nifti1Image(T1data, affine=classes.image.affine)\n",
" B1image = nib.nifti1.Nifti1Image(B1data, affine=classes.image.affine)\n",
" nib.save(M0image, 'M0.nii.gz')\n",
" nib.save(T1image, 'T1.nii.gz')\n",
" nib.save(B1image, 'B1.nii.gz')\n",
"\n",
"M0vals = np.array([0, 1, 0.8, 0.7, 0, 0, 0, 0, 0, 0])\n",
"T1vals = np.array([0, 3.0, 1.3, 0.9, 0, 0, 0, 0, 0, 0])\n",
"classes_to_params(M0vals, T1vals, 1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Simulate Image\n",
"\n",
"As well as fitting image data to find parameters, QUIT allows you to simulate the images from the paramters. Here we use the phantom images from the step above to create some simple gradient-echo (also called FLASH, SPGR or FFE depending on vendor) images.\n",
"\n",
"First we define a dictionary that sets the sequence parameters for the type of scan we are simulating. In this case we only need the TR and the flip-angle (FA). The TR is specified in seconds, not milliseconds, and there are multiple values for the flip-angle as we need to simulate multiple images.\n",
"\n",
"We then define and run the simulator. This will return a results object that contains the paths to the output files, which we use to display both volumes in the output image. The noise is defined as an SNR level relative assuming a nominal PD of 1 (this is why the PD of CSF was set to 1 above)."
]
},
{
"cell_type": "code",
"metadata": {},
"outputs": [],
"source": [
"sequence_parameters = {'SPGR': {'TR': 10e-3, 'FA': [3, 18]}}\n",
"simulator_results = DESPOT1Sim(sequence=d1seq, in_file='sim_spgr.nii.gz', noise=0.001, PD='PD.nii.gz', T1='T1.nii.gz').run()\n",
"\n",
"simulated_spgr = Layer(simulator_results.outputs.out_file, volume=0)\n",
"display(ns.three_plane(spgr))\n",
"spgr.volume = 1\n",
"display(ns.three_plane(spgr))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Fit Data and Compare to Reference\n",
"\n",
"Now we have simulated data we can run the DESPOT1/VFA fitting code and compare the results to our phantom reference data."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## B1 Mapping\n",
"\n",
"SPGR images are affected by B1 inhomogeneity. In the above simulation, we did not specify a B1 map and assumed that B1 was flat."
]
},
{
"cell_type": "code",
"metadata": {},
"outputs": [],
"source": [
"fitting_results = DESPOT1(sequence=d1seq, in_file='sim_spgr.nii.gz').run()\n",
"\n",
"display(ns.compare('T1.nii.gz', fitting_results.outputs.t1_map, title='T1 Comparison'))\n",
"display(ns.compare('PD.nii.gz', fitting_results.outputs.pd_map, title='PD Comparison'))"
]
}
],
"metadata": {
"kernel_info": {
"name": "python3"
},
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.4"
},
"nteract": {
"version": "0.8.4"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Loading

0 comments on commit a4263ad

Please sign in to comment.