-
Notifications
You must be signed in to change notification settings - Fork 21
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
34 changed files
with
1,070 additions
and
1,321 deletions.
There are no files selected for viewing
Empty file.
This file was deleted.
Oops, something went wrong.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -10,7 +10,7 @@ | |
| "\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 " | ||
| "In order to display the images, please install the `nanslice` package from `pip`." | ||
| ] | ||
| }, | ||
| { | ||
|
|
@@ -22,25 +22,22 @@ | |
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "metadata": {}, | ||
| "metadata": { | ||
| "tags": [] | ||
| }, | ||
| "outputs": [], | ||
| "source": [ | ||
| "%matplotlib inline\n", | ||
| "from QUIT.interfaces.relax import DESPOT1, DESPOT1Sim, DESPOT2\n", | ||
| "from nanslice import Layer\n", | ||
| "from qipype.interfaces.relax import DESPOT1, DESPOT1Sim\n", | ||
| "from qipype.sims import init_brainweb, make_phantom\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" | ||
| ] | ||
| "import nibabel as nib" | ||
| ], | ||
| "execution_count": null | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "## Create Phantom Data\n", | ||
| "\n", | ||
|
|
@@ -53,64 +50,25 @@ | |
| "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')" | ||
| "We now download the BrainWeb brain & T1 phantoms, if they are not already present in the current directory. QUIT provides a utility function for this." | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "metadata": {}, | ||
| "execution_count": null, | ||
| "metadata": { | ||
| "tags": [] | ||
| }, | ||
| "outputs": [], | ||
| "source": [ | ||
| "ns.three_plane(classes)" | ||
| "init_brainweb()" | ||
| ] | ||
| }, | ||
| { | ||
| "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", | ||
| "Now we use another QUIT function to convert the BrainWeb phantoms (which are tissue class maps) into parameter maps, in this case Proton Density and T1. There are 10 tissue classes:\n", | ||
| "\n", | ||
| "0. Background\n", | ||
| "1. CSF\n", | ||
|
|
@@ -123,32 +81,23 @@ | |
| "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." | ||
| "The `make_phantom` function requires a dictionary, where every parameter map we want to make is a key, and the value is a 10 element array with the values for each class. To keep things simple here, we only specify values for CSF, GM & WM, and set the other tissue types to zero. We also specify that we want a B1 map for later use, and sub-sample the images by a factor of 2 to make it faster." | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "metadata": {}, | ||
| "metadata": { | ||
| "tags": [] | ||
| }, | ||
| "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)" | ||
| ] | ||
| "param_dict = {'PD': [0, 1, 0.8, 0.7, 0, 0, 0, 0, 0, 0],\n", | ||
| " 'T1': [0, 3.0, 1.3, 0.9, 0, 0, 0, 0, 0, 0]}\n", | ||
| "make_phantom(param_dict, B1=True, subsamp=2)\n", | ||
| "display(ns.three_plane('T1.nii.gz'))\n", | ||
| "display(ns.three_plane('B1.nii.gz'))" | ||
| ], | ||
| "execution_count": null | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
|
|
@@ -169,13 +118,12 @@ | |
| "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", | ||
| "simulator_results = DESPOT1Sim(sequence=sequence_parameters, out_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))" | ||
| ] | ||
| "display(ns.three_plane(simulator_results.outputs.out_file, volume=0))\n", | ||
| "display(ns.three_plane(simulator_results.outputs.out_file, volume=1))" | ||
| ], | ||
| "execution_count": null | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
|
|
@@ -186,25 +134,85 @@ | |
| "Now we have simulated data we can run the DESPOT1/VFA fitting code and compare the results to our phantom reference data." | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "metadata": { | ||
| "tags": [] | ||
| }, | ||
| "outputs": [], | ||
| "source": [ | ||
| "fit = DESPOT1(sequence=sequence_parameters, in_file=simulator_results.outputs.out_file)\n", | ||
| "print(fit)\n", | ||
| "print(dir(fit))\n", | ||
| "fitting_results = fit.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'))" | ||
| ], | ||
| "execution_count": null | ||
| }, | ||
| { | ||
| "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." | ||
| "SPGR images are affected by B1 inhomogeneity. In the above simulation, we did not specify a B1 map and assumed that B1 was flat, but this is not true in real images. We will now add B1 into both our simulation and fitting.\n", | ||
| "\n", | ||
| "For DESPOT1, and most other methods, B1 is a \"fixed\" parameter, in contrast to T1 & PD which are \"varying\" parameters. This means that DESPOT1 does not fit a value for B1, it must be measured with another scan and is then fixed within each voxel. There are two methods in QUIT for calculating B1 maps (see https://quit.readthedocs.io/en/latest/Docs/B1.html).\n", | ||
| "\n", | ||
| "First we simulate our input images again, supplying the B1 map, but fit without supplying the B1 map. We add the `prefix` parameter to the fit to avoid overwriting the output of the previous fitting run." | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "metadata": { | ||
| "tags": [] | ||
| }, | ||
| "outputs": [], | ||
| "source": [ | ||
| "b1_sim_results = DESPOT1Sim(sequence=sequence_parameters, out_file='sim_b1.nii.gz', b1_map='B1.nii.gz', noise=0.001, PD='PD.nii.gz', T1='T1.nii.gz').run()\n", | ||
| "no_b1_fit_results = DESPOT1(prefix='no_b1', sequence=sequence_parameters, in_file=b1_sim_results.outputs.out_file).run()\n", | ||
| "display(ns.compare('T1.nii.gz', no_b1_fit_results.outputs.t1_map, clim=(0.8,1.5), title='Phantom Vs No-B1 in Fitting'))" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "metadata": { | ||
| "tags": [] | ||
| }, | ||
| "outputs": [], | ||
| "source": [ | ||
| "print(no_b1_fit_results.outputs.t1_map)" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "Finally, we fit again, this time supplying the B1 map, and compare the results" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "metadata": { | ||
| "tags": [] | ||
| }, | ||
| "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'))" | ||
| "b1_fit_results = DESPOT1(prefix='b1', sequence=sequence_parameters, in_file=b1_sim_results.outputs.out_file, b1_map='B1.nii.gz').run()\n", | ||
| "display(ns.compare(b1_fit_results.outputs.t1_map, no_b1_fit_results.outputs.t1_map, title='B1 in fitting versus no B1'))" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [] | ||
| } | ||
| ], | ||
| "metadata": { | ||
|
|
@@ -226,7 +234,7 @@ | |
| "name": "python", | ||
| "nbconvert_exporter": "python", | ||
| "pygments_lexer": "ipython3", | ||
| "version": "3.7.4" | ||
| "version": "3.7.6-final" | ||
| }, | ||
| "nteract": { | ||
| "version": "0.8.4" | ||
|
|
||
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.