Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • fsl/pytreat-practicals-2020
  • mchiew/pytreat-practicals-2020
  • ndcn0236/pytreat-practicals-2020
  • nichols/pytreat-practicals-2020
4 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 1522 additions and 381 deletions
getting_started/05_nifti.ipynb
View file @ 15d75bb5
%% Cell type:markdown id: tags:
# NIfTI images and python
The [nibabel](http://nipy.org/nibabel/) module is used to read and write NIfTI images and also some other medical imaging formats (e.g., ANALYZE, GIFTI, MINC, MGH). This module is included within the FSL python environment.
The [`nibabel`](http://nipy.org/nibabel/) module is used to read and write NIfTI
images and also some other medical imaging formats (e.g., ANALYZE, GIFTI,
MINC, MGH). `nibabel` is included within the FSL python environment.
Building upon `nibabel`, the
[`fslpy`](https://users.fmrib.ox.ac.uk/~paulmc/fsleyes/fslpy/latest/) library
contains a number of FSL-specific classes and functions which you may find
useful. But let's start with `nibabel` - `fslpy` is introduced in a different
practical (`advanced_topics/08_fslpy.ipynb`).
## Contents
* [Reading images](#reading-images)
* [Header info](#header-info)
* [Voxel sizes](#voxel-sizes)
* [Coordinate orientations and mappings](#orientation-info)
* [Voxel sizes](#voxel-sizes)
* [Coordinate orientations and mappings](#orientation-info)
* [Writing images](#writing-images)
* [Exercise](#exercise)
---
<a class="anchor" id="reading-images"></a>
## Reading images
It is easy to read an image:
For most neuroimaging dataformats reading an image is as simple as calling `nibabel.load`.
%% Cell type:code id: tags:
```
import numpy as np
import nibabel as nib
import os.path as op
filename = op.expandvars('${FSLDIR}/data/standard/MNI152_T1_1mm.nii.gz')
imobj = nib.load(filename, mmap=False)
# display header object
imhdr = imobj.header
# extract data (as an numpy array)
imdat = imobj.get_data().astype(float)
print(imdat.shape)
print('header', imhdr)
# extract data (as a numpy array)
imdat = imobj.get_fdata()
print('data', imdat.shape)
```
%% Cell type:markdown id: tags:
> Make sure you use the full filename, including the .nii.gz extension.
> Make sure you use the full filename, including the `.nii.gz` extension.
> `fslpy` provides FSL-like automatic file suffix detection though.
> We use the `expandvars()` function above to insert the FSLDIR
> environmental variable into our string. This function is
> discussed more fully in the file management practical.
> We use the expandvars() function above to insert the FSLDIR
>environmental variable into our string. This function is
>discussed more fully in the file management practical.
Reading the data off the disk is not done until `get_data()` is called.
Reading the data off the disk is not done until `get_fdata()` is called.
> Pitfall:
>
> The option `mmap=False`is necessary as turns off memory mapping, which otherwise would be invoked for uncompressed NIfTI files but not for compressed files. Since some functionality behaves differently on memory mapped objects, it is advisable to turn this off.
> The option `mmap=False` disables memory mapping, which would otherwise be
> invoked for uncompressed NIfTI files but not for compressed files. Since
> some functionality behaves differently on memory mapped objects, it is
> advisable to turn this off unless you specifically want it.
Once the data is read into a numpy array then it is easily manipulated.
> We recommend converting it to float at the start to avoid problems with integer arithmetic and overflow, though this is not compulsory.
> The `get_fdata` method will return floating point data, regardless of the
> underlying image data type. If you want the image data in the type that it
> is stored (e.g. integer ROI labels), then use
> `imdat = np.asanyarray(imobj.dataobj)` instead.
---
<a class="anchor" id="header-info"></a>
## Header info
There are many methods available on the header object - for example, look at `dir(imhdr)` or `help(imhdr)` or the [nibabel webpage about NIfTI images](http://nipy.org/nibabel/nifti_images.html)
There are many methods available on the header object - for example, look at
`dir(imhdr)` or `help(imhdr)` or the [nibabel webpage about NIfTI
images](http://nipy.org/nibabel/nifti_images.html)
<a class="anchor" id="voxel-sizes"></a>
### Voxel sizes
Dimensions of the voxels, in mm, can be found from:
%% Cell type:code id: tags:
```
voxsize = imhdr.get_zooms()
print(voxsize)
```
%% Cell type:markdown id: tags:
<a class="anchor" id="orientation-info"></a>
### Coordinate orientations and mappings
Information about the NIfTI qform and sform matrices can be extracted like this:
%% Cell type:code id: tags:
```
sform = imhdr.get_sform()
sformcode = imhdr['sform_code']
qform = imhdr.get_qform()
qformcode = imhdr['qform_code']
print(qformcode)
print(qform)
```
%% Cell type:markdown id: tags:
You can also get both code and matrix together like this:
%% Cell type:code id: tags:
```
affine, code = imhdr.get_qform(coded=True)
print(affine, code)
```
%% Cell type:markdown id: tags:
If you don't want to have to worry about the difference between `qform` and `sform`,
you can just let `nibabel` return what it thinks is the appropriate `affine`:
%% Cell type:code id: tags:
```
print('affine', imobj.affine)
```
%% Cell type:markdown id: tags:
> Note that we access the `affine` attribute from the image object here, not the image header (like above).
> Accessing the affine this way has the advantage that it will also work for data types, where the affine is stored in a different way in the header.
---
<a class="anchor" id="writing-images"></a>
## Writing images
If you have created a modified image by making or modifying a numpy array then you need to put this into a NIfTI image object in order to save it to a file. The easiest way to do this is to copy all the header info from an existing image like this:
If you have created a modified image by making or modifying a numpy array then
you need to put this into a NIfTI image object in order to save it to a file.
The easiest way to do this is to copy all the header info from an existing
image like this:
%% Cell type:code id: tags:
```
newdata = imdat * imdat
newhdr = imhdr.copy()
newobj = nib.nifti1.Nifti1Image(newdata, None, header=newhdr)
nib.save(newobj, "mynewname.nii.gz")
```
%% Cell type:markdown id: tags:
where `newdata` is the numpy array (the above is a random example only) and `imhdr` is the existing image header (as above).
where `newdata` is the numpy array (the above is a random example only) and
`imhdr` is the existing image header (as above).
> It is possible to also just pass in an affine matrix rather than a
> copied header, but we *strongly* recommend against this if you are
> processing an existing image as otherwise you run the risk of
> swapping the left-right orientation. Those that have used
> `save_avw` in matlab may well have been bitten in this way in the
> past. Therefore, copy a header from one of the input images
> whenever possible, and just use the affine matrix option if you are
> creating an entirely separate image, like a simulation.
If the voxel size of the image is different, then extra modifications will be required. For this, or for building an image from scratch, see the [nibabel documentation](http://nipy.org/nibabel/nifti_images.html) on NIfTI images.
If the voxel size of the image is different, then extra modifications will be
required. Take a look at the `fslpy` practical for some extra image
manipulation options, including cropping and resampling
(`advanced_topics/08_fslpy.ipynb`).
---
<a class="anchor" id="exercise"></a>
<a class="anchor" id="exercises"></a>
## Exercise
Write some code to read in a 4D fMRI image (you can find one [here] if
you don't have one handy), calculate the tSNR and then save the 3D result.
Write some code to read in a 4D fMRI image (you can find one
[here](http://www.fmrib.ox.ac.uk/~mark/files/av.nii.gz) if you don't have one
handy), calculate the tSNR and then save the 3D result.
> The tSNR of a time series signal is simply its mean divided by its standard
> deviation.
%% Cell type:code id: tags:
```
# Calculate tSNR
```
......
getting_started/05_nifti.md
View file @ 15d75bb5
# NIfTI images and python
The [nibabel](http://nipy.org/nibabel/) module is used to read and write NIfTI images and also some other medical imaging formats (e.g., ANALYZE, GIFTI, MINC, MGH). This module is included within the FSL python environment.
The [`nibabel`](http://nipy.org/nibabel/) module is used to read and write NIfTI
images and also some other medical imaging formats (e.g., ANALYZE, GIFTI,
MINC, MGH). `nibabel` is included within the FSL python environment.
Building upon `nibabel`, the
[`fslpy`](https://users.fmrib.ox.ac.uk/~paulmc/fsleyes/fslpy/latest/) library
contains a number of FSL-specific classes and functions which you may find
useful. But let's start with `nibabel` - `fslpy` is introduced in a different
practical (`advanced_topics/08_fslpy.ipynb`).
## Contents
* [Reading images](#reading-images)
* [Header info](#header-info)
* [Voxel sizes](#voxel-sizes)
* [Coordinate orientations and mappings](#orientation-info)
* [Voxel sizes](#voxel-sizes)
* [Coordinate orientations and mappings](#orientation-info)
* [Writing images](#writing-images)
* [Exercise](#exercise)
......@@ -16,7 +26,7 @@ The [nibabel](http://nipy.org/nibabel/) module is used to read and write NIfTI i
<a class="anchor" id="reading-images"></a>
## Reading images
It is easy to read an image:
For most neuroimaging dataformats reading an image is as simple as calling `nibabel.load`.
```
import numpy as np
......@@ -24,36 +34,47 @@ import nibabel as nib
import os.path as op
filename = op.expandvars('${FSLDIR}/data/standard/MNI152_T1_1mm.nii.gz')
imobj = nib.load(filename, mmap=False)
# display header object
imhdr = imobj.header
# extract data (as an numpy array)
imdat = imobj.get_data().astype(float)
print(imdat.shape)
print('header', imhdr)
# extract data (as a numpy array)
imdat = imobj.get_fdata()
print('data', imdat.shape)
```
> Make sure you use the full filename, including the .nii.gz extension.
> Make sure you use the full filename, including the `.nii.gz` extension.
> `fslpy` provides FSL-like automatic file suffix detection though.
> We use the `expandvars()` function above to insert the FSLDIR
> environmental variable into our string. This function is
> discussed more fully in the file management practical.
> We use the expandvars() function above to insert the FSLDIR
>environmental variable into our string. This function is
>discussed more fully in the file management practical.
Reading the data off the disk is not done until `get_data()` is called.
Reading the data off the disk is not done until `get_fdata()` is called.
> Pitfall:
>
> The option `mmap=False`is necessary as turns off memory mapping, which otherwise would be invoked for uncompressed NIfTI files but not for compressed files. Since some functionality behaves differently on memory mapped objects, it is advisable to turn this off.
> The option `mmap=False` disables memory mapping, which would otherwise be
> invoked for uncompressed NIfTI files but not for compressed files. Since
> some functionality behaves differently on memory mapped objects, it is
> advisable to turn this off unless you specifically want it.
Once the data is read into a numpy array then it is easily manipulated.
> We recommend converting it to float at the start to avoid problems with integer arithmetic and overflow, though this is not compulsory.
> The `get_fdata` method will return floating point data, regardless of the
> underlying image data type. If you want the image data in the type that it
> is stored (e.g. integer ROI labels), then use
> `imdat = np.asanyarray(imobj.dataobj)` instead.
---
<a class="anchor" id="header-info"></a>
## Header info
There are many methods available on the header object - for example, look at `dir(imhdr)` or `help(imhdr)` or the [nibabel webpage about NIfTI images](http://nipy.org/nibabel/nifti_images.html)
There are many methods available on the header object - for example, look at
`dir(imhdr)` or `help(imhdr)` or the [nibabel webpage about NIfTI
images](http://nipy.org/nibabel/nifti_images.html)
<a class="anchor" id="voxel-sizes"></a>
### Voxel sizes
......@@ -85,13 +106,24 @@ affine, code = imhdr.get_qform(coded=True)
print(affine, code)
```
If you don't want to have to worry about the difference between `qform` and `sform`,
you can just let `nibabel` return what it thinks is the appropriate `affine`:
```
print('affine', imobj.affine)
```
> Note that we access the `affine` attribute from the image object here, not the image header (like above).
> Accessing the affine this way has the advantage that it will also work for data types, where the affine is stored in a different way in the header.
---
<a class="anchor" id="writing-images"></a>
## Writing images
If you have created a modified image by making or modifying a numpy array then you need to put this into a NIfTI image object in order to save it to a file. The easiest way to do this is to copy all the header info from an existing image like this:
If you have created a modified image by making or modifying a numpy array then
you need to put this into a NIfTI image object in order to save it to a file.
The easiest way to do this is to copy all the header info from an existing
image like this:
```
newdata = imdat * imdat
......@@ -99,7 +131,9 @@ newhdr = imhdr.copy()
newobj = nib.nifti1.Nifti1Image(newdata, None, header=newhdr)
nib.save(newobj, "mynewname.nii.gz")
```
where `newdata` is the numpy array (the above is a random example only) and `imhdr` is the existing image header (as above).
where `newdata` is the numpy array (the above is a random example only) and
`imhdr` is the existing image header (as above).
> It is possible to also just pass in an affine matrix rather than a
> copied header, but we *strongly* recommend against this if you are
......@@ -110,17 +144,25 @@ where `newdata` is the numpy array (the above is a random example only) and `imh
> whenever possible, and just use the affine matrix option if you are
> creating an entirely separate image, like a simulation.
If the voxel size of the image is different, then extra modifications will be required. For this, or for building an image from scratch, see the [nibabel documentation](http://nipy.org/nibabel/nifti_images.html) on NIfTI images.
If the voxel size of the image is different, then extra modifications will be
required. Take a look at the `fslpy` practical for some extra image
manipulation options, including cropping and resampling
(`advanced_topics/08_fslpy.ipynb`).
---
<a class="anchor" id="exercise"></a>
<a class="anchor" id="exercises"></a>
## Exercise
Write some code to read in a 4D fMRI image (you can find one [here](http://www.fmrib.ox.ac.uk/~mark/files/av.nii.gz) if
you don't have one handy), calculate the tSNR and then save the 3D result.
Write some code to read in a 4D fMRI image (you can find one
[here](http://www.fmrib.ox.ac.uk/~mark/files/av.nii.gz) if you don't have one
handy), calculate the tSNR and then save the 3D result.
> The tSNR of a time series signal is simply its mean divided by its standard
> deviation.
```
# Calculate tSNR
```
getting_started/06_plotting.ipynb
View file @ 15d75bb5
Source diff could not be displayed: it is too large. Options to address this: view the blob.
getting_started/06_plotting.md
View file @ 15d75bb5
......@@ -109,7 +109,7 @@ there is also an alternative: `scatter()`
```
fig, ax = plt.subplots()
# setup some sizes for each point (arbitrarily example here)
ssize = 100*abs(samp1-samp2) + 10
ssize = 100*abs(samp1-samp2) + 10
ax.scatter(samp1, samp2, s=ssize, alpha=0.5)
# now add the y=x line
allsamps = np.hstack((samp1,samp2))
......@@ -153,11 +153,28 @@ import nibabel as nib
import os.path as op
nim = nib.load(op.expandvars('${FSLDIR}/data/standard/MNI152_T1_1mm.nii.gz'), mmap=False)
imdat = nim.get_data().astype(float)
plt.imshow(imdat[:,:,70], cmap=plt.cm.gray)
imslc = imdat[:,:,70]
plt.imshow(imslc, cmap=plt.cm.gray)
plt.colorbar()
plt.grid('off')
```
Note that matplotlib will use the **voxel data orientation**, and that
configuring the plot orientation is **your responsibility**. To rotate a
slice, simply transpose the data (`.T`). To invert the data along along an
axis, you don't need to modify the data - simply swap the axis limits around:
```
plt.imshow(imslc.T, cmap=plt.cm.gray)
plt.xlim(reversed(plt.xlim()))
plt.ylim(reversed(plt.ylim()))
plt.colorbar()
plt.grid('off')
```
<a class="anchor" id="3D-plots"></a>
### 3D plots
......@@ -208,5 +225,3 @@ example code from the docs).
```
# Make up some data and do the funky plot
```
getting_started/07_jupyter.ipynb
View file @ 15d75bb5
%% Cell type:markdown id: tags:
# Jupyter notebook and IPython
Our main interaction with python so far has been through the [Jupyter notebook](http://jupyter.org/).
These notebooks are extremely popular these days within the python scientific community, however they support many more languages, such as R and octave (and even matlab with the right [plugin](https://github.com/Calysto/matlab_kernel)).
They allow for interactive analysis of your data interspersed by explanatory notes (including LaTeX) with inline plotting.
However, they can not be called as scripts on the command line or be imported from other python code, which makes them rather stand-alone.
This makes them more useful for analysis that needs to be reproducible, but does not need to be replicated on different datasets (e.g., making a plot for a paper).
For more ad-hoc analysis it can be useful to just use the command line (i.e., a REPL).
We strongly recommend to use the IPython (available as `ipython` or `fslipython`) rather than default python REPL (available through `python` or `fslpython`)
Our main interaction with python so far has been through the [Jupyter
notebook](http://jupyter.org/). These notebooks are extremely popular these
days within the python scientific community, however they support many more
languages, such as R and octave (and even matlab with the right
[plugin](https://github.com/Calysto/matlab_kernel)). They allow for
interactive analysis of your data interspersed by explanatory notes (including
LaTeX) with inline plotting. However, they can not be called as scripts on
the command line or be imported from other python code, which makes them
rather stand-alone. This makes them more useful for analysis that needs to be
reproducible, but does not need to be replicated on different datasets (e.g.,
making a plot for a paper).
For more ad-hoc analysis it can be useful to just use the command line (i.e.,
a REPL<sup>*</sup>). We strongly recommend to use the IPython (available as
`fslipython` or `ipython`) rather than default python REPL (available through
`fslpython` or `python`), as IPython is much more user-friendly.
> <sup>*</sup>REPL = **R**ead-**E**val-**P**rint-**L**oop - the geeky term for
> an interactive prompt. You may hear younger generations using the term
> [ESRR](https://www.youtube.com/watch?v=wBoRkg5-Ieg) instead.
Both Ipython and the jupyter notebook offer a whole range of magic commands, which all start with a `%` sign.
* A magic command starting with a single `%` sign will only affect the single line.
* A magic command starting with two '%' signs will affect the whole block of code.
Note that the normal python interpreter will not understand these magic commands, so you will have to take them out when writing a python script or library.
Here we will discuss some of the many features available to you in Ipython and the Jupyter notebook
---
## Getting help
To get the documentation for any object or method simply append a question mark
%% Cell type:code id: tags:
```
import string
string.capwords?
```
%% Cell type:markdown id: tags:
This also works for any of the magic commands discussed below
%% Cell type:code id: tags:
```
%run?
```
%% Cell type:markdown id: tags:
Alternatively you can put two questions marks to get the complete code for the method or object class
%% Cell type:code id: tags:
```
import string
string.capwords??
```
%% Cell type:markdown id: tags:
Both Ipython and Jupyter also come with autocomplete, which is available at any time by pressing the tab key
---
## Running shell commands
Commands starting with a `!` will be sent to the shell rather than the python interpreter.
%% Cell type:code id: tags:
```
!fslstats ${FSLDIR}/data/standard/FMRIB58_FA_1mm.nii.gz -r
```
%% Cell type:markdown id: tags:
You can even capture the output from the shell command in a variable:
%% Cell type:code id: tags:
```
r = !fslstats ${FSLDIR}/data/standard/FMRIB58_FA_1mm.nii.gz -r
r_lower, r_upper = [float(element) for element in r[0].split()]
print('Bounds are ({:.0f}, {:.0f})'.format(r_lower, r_upper))
```
%% Cell type:markdown id: tags:
---
## Running python scripts
We could run a python script as a shell command above. However, it will often be more convenient to use `%run` instead.
> ```
> %run <python script> <arguments...>
> ```
Arguments are provided in exactly the same way as if you called `python` in the shell. The main advantages are:
- Any top-level variables will be made available to you after the script finishes
- All the debugging/timing/profiling tools discussed below will be available to you
A common workflow, when writing a python script is to have an Ipython terminal open next to your text editor and regularly use %run to test the script
---
## Running other programming languages
In the notebook you can include a whole code block using another language by using `%%<language>` (for many languages you will have to install a toolkit first, just google your favorite language besides python)
%% Cell type:code id: tags:
```
%%bash
for filename in `ls *.md` ; do
head -n 1 ${filename}
done
```
%% Cell type:markdown id: tags:
---
## Timing/profiling code
We can time a line of code with `%time` or a whole code block using `%%time`.
To get the time needed to calculate the sine of a million random numbers:
%% Cell type:code id: tags:
```
import numpy as np
numbers = np.random.rand(int(1e6))
%time np.sin(numbers)
```
%% Cell type:markdown id: tags:
For very fast evaluation, you might need to run it multiple times to get an accurate estimate. The `%timeit` (or `%%timeit` for a code block) takes care of this for you.
%% Cell type:code id: tags:
```
import numpy as np
numbers = np.random.rand(10)
%timeit np.sin(numbers) # this will take a few seconds to run
```
%% Cell type:markdown id: tags:
Finally, if you want to figure out what part of the code is actually slowing you down you can use `%prun`, which gives you an overview of how long the interpreter spent in each method:
%% Cell type:code id: tags:
```
import nibabel as nib
import os.path as op
%prun nib.load(op.expandvars('${FSLDIR}/data/standard/FMRIB58_FA_1mm.nii.gz'))
```
%% Cell type:markdown id: tags:
---
## Debugging
Despite your best efforts in many cases some error will crop up
%% Cell type:code id: tags:
```
import numpy as np
def total(a_list):
"""Calculate the total of a list.
This is a very naive (not recommended) and bugged implementation
"""
# create local copy befor changing the input
local_list = list(a_list)
total = 0.
while len(local_list) > 0:
total += local_list.pop(1) # returns element at index=1 and removes it
return total
print(total([2, 3, 4]))
```
%% Cell type:markdown id: tags:
You can always open a debugger at the location of the last error by using the `%debug` magic command. You can find a list of commands available in the debugger [here](http://www.georgejhunt.com/olpc/pydebug/pydebug/ipdb.html)
%% Cell type:code id: tags:
```
%debug
```
%% Cell type:markdown id: tags:
Try to check the value of `a_list` and `local_list` from within the debugger.
> WARNING: you need to quit the debugger before any further commands will run (type `q` into the prompt)!
If you always want to enter the debugger when an error is raised you can call `%pdb on` at any time (call `%pdf off` to reverse this)
---
## Enabling plotting
By far the most popular scientific plotting library is [matplotlib](https://matplotlib.org/).
You can enable plotting in Ipython or the jupyter notebook using `%matplotlib <backend>`, where [backend](https://matplotlib.org/faq/usage_faq.html#what-is-a-backend) is the system that will be used to display the plots.
When failing to provide a backend it will simply use the default (which is usually fine).
* In the jupyter notebook use the `nbagg` backend for interactive plots or the `inline` backend for non-interactive plots
* Otherwise on Mac OSx use the `macosx` backend
%% Cell type:code id: tags:
```
%matplotlib nbagg
```
%% Cell type:markdown id: tags:
> Keep in mind that as soon as you have started plotting you can no longer change your backend without exiting the python interpreter and restarting `python` (note that in the jupyter notebook you can just press `Restart` in the `Kernel` menu).
To do the equivalent in a python script would look like
> ```
> import matplotlib as mpl
> mpl.use(<backend>)
> ```
For interactive use it can be handy to have all the `numpy` numeric functions and `matplotlib` plotting functions directly available without importing them explicitly.
This can be achieved using the `%pylab <backend>` magic command.
%% Cell type:code id: tags:
```
%pylab nbagg
```
%% Cell type:markdown id: tags:
This is equivalent in python code to:
> ```
> import matplotlib as mpl
> mpl.use(<backend>)
> from matplotlib.pylab import *
> ```
> The last line imports everything from the matplotlib.pylab module into the namespace.
I start most of my notebooks or terminals with the `%pylab` command, because afterwards I can just do stuff like:
%% Cell type:code id: tags:
```
x = linspace(0, pi, 301)
y = sin(x)
plot(x, y, 'r-')
```
%% Cell type:markdown id: tags:
The main disadvantage is that it will not be obvious to the naive reader of this code, whether functions like `linspace`, `sin`, or `plot` are originate from numpy, matplotlib, or are built-in.
This is why we dont recommend `from <module> import *` statements in any longer code or code you intend to share.
---
## Exporting code from the Jupyter notebook
If you have a code cell in the jupyter notebook, that you want to convert into a script, you can use the %%writefile
%% Cell type:code id: tags:
```
%%writefile script_from_notebook.py
# a bunch of imports
import numpy as np
from datetime import datetime
```
%% Cell type:markdown id: tags:
Any additional code cells need to contain the `-a` flag to stop jupyter from overwriting the original code
%% Cell type:code id: tags:
```
%%writefile -a script_from_notebook.py
print('today is ', datetime.now())
print('sin(3) is ', np.sin(3))
```
%% Cell type:markdown id: tags:
We can now run this script
%% Cell type:code id: tags:
```
%run script_from_notebook.py
```
%% Cell type:markdown id: tags:
---
## Exporting code from the Ipython terminal
You can access the full history of your session using `%history`.
To save the history to a file use `%history -f <filename>`.
You will probably have to clean a lot of erroneous commands you typed from that file before you are able to run it as a script.
......
getting_started/07_jupyter.md
View file @ 15d75bb5
# Jupyter notebook and IPython
Our main interaction with python so far has been through the [Jupyter notebook](http://jupyter.org/).
These notebooks are extremely popular these days within the python scientific community, however they support many more languages, such as R and octave (and even matlab with the right [plugin](https://github.com/Calysto/matlab_kernel)).
They allow for interactive analysis of your data interspersed by explanatory notes (including LaTeX) with inline plotting.
However, they can not be called as scripts on the command line or be imported from other python code, which makes them rather stand-alone.
This makes them more useful for analysis that needs to be reproducible, but does not need to be replicated on different datasets (e.g., making a plot for a paper).
For more ad-hoc analysis it can be useful to just use the command line (i.e., a REPL).
We strongly recommend to use the IPython (available as `ipython` or `fslipython`) rather than default python REPL (available through `python` or `fslpython`)
Our main interaction with python so far has been through the [Jupyter
notebook](http://jupyter.org/). These notebooks are extremely popular these
days within the python scientific community, however they support many more
languages, such as R and octave (and even matlab with the right
[plugin](https://github.com/Calysto/matlab_kernel)). They allow for
interactive analysis of your data interspersed by explanatory notes (including
LaTeX) with inline plotting. However, they can not be called as scripts on
the command line or be imported from other python code, which makes them
rather stand-alone. This makes them more useful for analysis that needs to be
reproducible, but does not need to be replicated on different datasets (e.g.,
making a plot for a paper).
For more ad-hoc analysis it can be useful to just use the command line (i.e.,
a REPL<sup>*</sup>). We strongly recommend to use the IPython (available as
`fslipython` or `ipython`) rather than default python REPL (available through
`fslpython` or `python`), as IPython is much more user-friendly.
> <sup>*</sup>REPL = **R**ead-**E**val-**P**rint-**L**oop - the geeky term for
> an interactive prompt. You may hear younger generations using the term
> [ESRR](https://www.youtube.com/watch?v=wBoRkg5-Ieg) instead.
Both Ipython and the jupyter notebook offer a whole range of magic commands, which all start with a `%` sign.
* A magic command starting with a single `%` sign will only affect the single line.
......@@ -206,4 +219,3 @@ We can now run this script
You can access the full history of your session using `%history`.
To save the history to a file use `%history -f <filename>`.
You will probably have to clean a lot of erroneous commands you typed from that file before you are able to run it as a script.
getting_started/08_scripts.ipynb
View file @ 15d75bb5
%% Cell type:markdown id: tags:
# Callable scripts in python
In this tutorial we will cover how to write simple stand-alone scripts in python that can be used as alternatives to bash scripts.
In this tutorial we will cover how to write simple stand-alone scripts in
python that can be used as alternatives to bash scripts.
There are some code blocks within this webpage, but for this practical we _**strongly
recommend that you write the code in an IDE or editor**_ instead and then run the scripts from a terminal.
**Important**: Throughout this series of practicals we have been working
entirely within the Jupyter notebook environment. But it's now time to
graduate to writing *real* Python scripts, and running them within a
*real* enviromnent.
So within this practical there are some code blocks, but instead of running
them inside the notebook we **strongly recommend that you write the code in
an IDE or editor**,and then run the scripts from a terminal. [Don't
panic](https://www.youtube.com/watch?v=KojYatpLPSE), we're right here,
ready to help.
## Contents
* [Basic script](#basic-script)
* [Calling other executables](#calling-other-executables)
* [Command line arguments](#command-line-arguments)
* [Example script](#example-script)
* [Exercise](#exercise)
---
<a class="anchor" id="basic-script"></a>
## Basic script
The first line of a python script is usually:
%% Cell type:code id: tags:
``` python
```
#!/usr/bin/env python
```
%% Cell type:markdown id: tags:
which invokes whichever version of python can be found by `/usr/bin/env` since python can be located in many different places.
For FSL scripts we use an alternative, to ensure that we pick up the version of python (and associated packages) that we ship with FSL. To do this we use the line:
%% Cell type:code id: tags:
``` python
```
#!/usr/bin/env fslpython
```
%% Cell type:markdown id: tags:
After this line the rest of the file just uses regular python syntax, as in the other tutorials. Make sure you make the file executable - just like a bash script.
<a class="anchor" id="calling-other-executables"></a>
## Calling other executables
The most essential call that you need to use to replicate the way a bash script calls executables is `subprocess.run()`. A simple call looks like this:
%% Cell type:code id: tags:
``` python
```
import subprocess as sp
import shlex
sp.run(['ls', '-la'])
```
%% Cell type:markdown id: tags:
> Passing the arguments as a list is good practice and improves the safety of
> the call.
To suppress the output do this:
%% Cell type:code id: tags:
``` python
```
spobj = sp.run(['ls'], stdout = sp.PIPE)
```
%% Cell type:markdown id: tags:
To store the output do this:
%% Cell type:code id: tags:
``` python
spobj = sp.run('ls -la'.split(), stdout = sp.PIPE)
```
spobj = sp.run(shlex.split('ls -la'), stdout = sp.PIPE)
sout = spobj.stdout.decode('utf-8')
print(sout)
```
%% Cell type:markdown id: tags:
> Note that the `decode` call in the middle line converts the string from a byte string to a normal string. In Python 3 there is a distinction between strings (sequences of characters, possibly using multiple bytes to store each character) and bytes (sequences of bytes). The world has moved on from ASCII, so in this day and age, this distinction is absolutely necessary, and Python does a fairly good job of it.
> shlex.split and shlex.quote are functions designed to break up and quote
> (respectively) shell command lines and arguments. Quoting of user provided
> arguments helps to prevent unintended consequences from inappropriate inputs.
>
> Note that the `decode` call in the middle line converts the string from a byte
> string to a normal string. In Python 3 there is a distinction between strings
> (sequences of characters, possibly using multiple bytes to store each
> character) and bytes (sequences of bytes). The world has moved on from ASCII,
> so in this day and age, this distinction is absolutely necessary, and Python
> does a fairly good job of it.
If the output is numerical then this can be extracted like this:
%% Cell type:code id: tags:
``` python
```
import os
fsldir = os.getenv('FSLDIR')
spobj = sp.run([fsldir+'/bin/fslstats', fsldir+'/data/standard/MNI152_T1_1mm_brain', '-V'], stdout = sp.PIPE)
sout = spobj.stdout.decode('utf-8')
vol_vox = float(sout.split()[0])
vol_mm = float(sout.split()[1])
results = sout.split()
vol_vox = float(results[0])
vol_mm = float(results[1])
print('Volumes are: ', vol_vox, ' in voxels and ', vol_mm, ' in mm')
```
%% Cell type:markdown id: tags:
An alternative way to run a set of commands would be like this:
%% Cell type:code id: tags:
``` python
```
import shlex
commands = """
{fsldir}/bin/fslmaths {t1} -bin {t1_mask}
{fsldir}/bin/fslmaths {t2} -mas {t1_mask} {t2_masked}
"""
fsldirpath = os.getenv('FSLDIR')
commands = commands.format(t1 = 't1.nii.gz', t1_mask = 't1_mask', t2 = 't2', t2_masked = 't2_masked', fsldir = fsldirpath)
sout=[]
for cmd in commands.split('\n'):
for cmd in commands.splitlines():
if cmd: # avoids empty strings getting passed to sp.run()
print('Running command: ', cmd)
spobj = sp.run(cmd.split(), stdout = sp.PIPE)
spobj = sp.run(shlex.split(cmd), stdout = sp.PIPE)
sout.append(spobj.stdout.decode('utf-8'))
```
%% Cell type:markdown id: tags:
> Don't be tempted to use the shell=True argument to subprocess.run, especially
> if you are dealing with user input - if the user gave
> *myfile; rm -f ~*
> as a file name and you called the command with shell=True **and** you
> passed the command in as a string then bad things happen!
>
> The safe way to use these kinds of inputs is to pass them through shlex.quote()
> before sending.
>
> ```a = shlex.quote('myfile; rm -f ~')
> cmd = "ls {}".format(a)
> sp.run(shlex.split(cmd))```
> If you're calling lots of FSL tools, the `fslpy` library has a number of
> *wrapper* functions, which can be used to call an FSL command directly
> from Python - check out `advanced_topics/08_fslpy.ipynb`.
<a class="anchor" id="command-line-arguments"></a>
## Command line arguments
The simplest way of dealing with command line arguments is use the module `sys`, which gives access to an `argv` list:
%% Cell type:code id: tags:
``` python
```
import sys
print(len(sys.argv))
print(sys.argv[0])
```
%% Cell type:markdown id: tags:
For more sophisticated argument parsing you can use `argparse` - good documentation and examples of this can be found on the web.
> argparse can automatically produce help text for the user, validate input etc., so it is strongly recommended.
---
<a class="anchor" id="example-script"></a>
## Example script
Here is a simple bash script (it masks an image and calculates volumes - just as a random example). DO NOT execute the code blocks here within the notebook/webpage:
%% Cell type:code id: tags:
``` python
```
#!/bin/bash
if [ $# -lt 2 ] ; then
echo "Usage: $0 <input image> <output image>"
exit 1
fi
infile=$1
outfile=$2
# mask input image with MNI
$FSLDIR/bin/fslmaths $infile -mas $FSLDIR/data/standard/MNI152_T1_1mm_brain $outfile
# calculate volumes of masked image
vv=`$FSLDIR/bin/fslstats $outfile -V`
vol_vox=`echo $vv | awk '{ print $1 }'`
vol_mm=`echo $vv | awk '{ print $2 }'`
echo "Volumes are: $vol_vox in voxels and $vol_mm in mm"
```
%% Cell type:markdown id: tags:
And an alternative in python:
%% Cell type:code id: tags:
``` python
```
#!/usr/bin/env fslpython
import os, sys
import subprocess as sp
fsldir=os.getenv('FSLDIR')
if len(sys.argv)<2:
print('Usage: ', sys.argv[0], ' <input image> <output image>')
sys.exit(1)
infile = sys.argv[1]
outfile = sys.argv[2]
# mask input image with MNI
spobj = sp.run([fsldir+'/bin/fslmaths', infile, '-mas', fsldir+'/data/standard/MNI152_T1_1mm_brain', outfile], stdout = sp.PIPE)
# calculate volumes of masked image
spobj = sp.run([fsldir+'/bin/fslstats', outfile, '-V'], stdout = sp.PIPE)
sout = spobj.stdout.decode('utf-8')
vol_vox = float(sout.split()[0])
vol_mm = float(sout.split()[1])
results = sout.split()
vol_vox = float(results[0])
vol_mm = float(results[1])
print('Volumes are: ', vol_vox, ' in voxels and ', vol_mm, ' in mm')
```
%% Cell type:markdown id: tags:
---
<a class="anchor" id="exercise"></a>
## Exercise
Write a simple version of fslstats that is able to calculate either a
mean or a _sum_ (and hence can do something that fslstats cannot!)
%% Cell type:code id: tags:
``` python
```
# Don't write anything here - do it in a standalone script!
```
......
getting_started/08_scripts.md
View file @ 15d75bb5
# Callable scripts in python
In this tutorial we will cover how to write simple stand-alone scripts in python that can be used as alternatives to bash scripts.
In this tutorial we will cover how to write simple stand-alone scripts in
python that can be used as alternatives to bash scripts.
There are some code blocks within this webpage, but for this practical we _**strongly
recommend that you write the code in an IDE or editor**_ instead and then run the scripts from a terminal.
**Important**: Throughout this series of practicals we have been working
entirely within the Jupyter notebook environment. But it's now time to
graduate to writing *real* Python scripts, and running them within a
*real* enviromnent.
So within this practical there are some code blocks, but instead of running
them inside the notebook we **strongly recommend that you write the code in
an IDE or editor**,and then run the scripts from a terminal. [Don't
panic](https://www.youtube.com/watch?v=KojYatpLPSE), we're right here,
ready to help.
## Contents
......@@ -63,10 +72,10 @@ print(sout)
> arguments helps to prevent unintended consequences from inappropriate inputs.
>
> Note that the `decode` call in the middle line converts the string from a byte
> string to a normal string. In Python 3 there is a distinction between strings
> (sequences of characters, possibly using multiple bytes to store each
> character) and bytes (sequences of bytes). The world has moved on from ASCII,
> so in this day and age, this distinction is absolutely necessary, and Python
> string to a normal string. In Python 3 there is a distinction between strings
> (sequences of characters, possibly using multiple bytes to store each
> character) and bytes (sequences of bytes). The world has moved on from ASCII,
> so in this day and age, this distinction is absolutely necessary, and Python
> does a fairly good job of it.
If the output is numerical then this can be extracted like this:
......@@ -100,7 +109,7 @@ for cmd in commands.splitlines():
sout.append(spobj.stdout.decode('utf-8'))
```
> Don't be tempted to use the shell=True argument to subprocess.run, especially
> Don't be tempted to use the shell=True argument to subprocess.run, especially
> if you are dealing with user input - if the user gave
> *myfile; rm -f ~*
> as a file name and you called the command with shell=True **and** you
......@@ -113,6 +122,11 @@ for cmd in commands.splitlines():
> cmd = "ls {}".format(a)
> sp.run(shlex.split(cmd))```
> If you're calling lots of FSL tools, the `fslpy` library has a number of
> *wrapper* functions, which can be used to call an FSL command directly
> from Python - check out `advanced_topics/08_fslpy.ipynb`.
<a class="anchor" id="command-line-arguments"></a>
## Command line arguments
......@@ -144,7 +158,7 @@ infile=$1
outfile=$2
# mask input image with MNI
$FSLDIR/bin/fslmaths $infile -mas $FSLDIR/data/standard/MNI152_T1_1mm_brain $outfile
# calculate volumes of masked image
# calculate volumes of masked image
vv=`$FSLDIR/bin/fslstats $outfile -V`
vol_vox=`echo $vv | awk '{ print $1 }'`
vol_mm=`echo $vv | awk '{ print $2 }'`
......@@ -166,7 +180,7 @@ infile = sys.argv[1]
outfile = sys.argv[2]
# mask input image with MNI
spobj = sp.run([fsldir+'/bin/fslmaths', infile, '-mas', fsldir+'/data/standard/MNI152_T1_1mm_brain', outfile], stdout = sp.PIPE)
# calculate volumes of masked image
# calculate volumes of masked image
spobj = sp.run([fsldir+'/bin/fslstats', outfile, '-V'], stdout = sp.PIPE)
sout = spobj.stdout.decode('utf-8')
results = sout.split()
......@@ -186,4 +200,3 @@ mean or a _sum_ (and hence can do something that fslstats cannot!)
```
# Don't write anything here - do it in a standalone script!
```
getting_started/README.md
View file @ 15d75bb5
Getting started
===============
This directory contains the "Getting started" practical.
This directory contains the _Getting started_ practical, intended for those of
you who are new to Python as a language, or want a refresher on how it all
works.
If you are already comfortable with Python, `numpy`, `nibabel`, and `jupyter`.
then you might want to check out the `advanced_topics` folder.
1. Python basics (MJ)
2. Text input/output (MC)
......@@ -9,5 +17,5 @@ This directory contains the "Getting started" practical.
4. Numpy (PM)
5. Nifti images (MJ)
6. Image manipulation (MC)
7. Jupyter notebook and Ipython(MC)
7. Jupyter notebook and IPython (MC)
8. Writing a callable script (MJ)
talks/introduction/pytreat_intro.ipynb
View file @ 15d75bb5
%% Cell type:markdown id: tags:
# Welcome to the WIN PyTreat 2018!
# Welcome to the WIN PyTreat 2020!
Program: https://docs.google.com/document/d/10CwLEhUi-YiwfC2F40QCVm6eEVwKiaXkfTKz67xWAfM/edit?usp=sharing
Program: https://docs.google.com/document/d/11LQgxC-LZPG_TXS3MP9tYXNWABAQdTI00iavk-tKttU/edit
__We need to do some setting up, so get your laptop ready!__
__Make sure you have FSL 5.0.10 installed and working!__
__Make sure you have FSL 6.0.3 installed and working!__
__End all sentences with an exclamation mark!__
__Open this page in your web browser!__
https://git.fmrib.ox.ac.uk/fsl/pytreat-2018-practicals/tree/master/talks/introduction/pytreat_intro.md
https://git.fmrib.ox.ac.uk/fsl/pytreat-practicals-2020/tree/master/talks/introduction/pytreat_intro.md
## Overview
* [Python in a nutshell](#python-in-a-nutshell)
* [`fslpython`](#fslpython)
* [Running Python scripts](#running-python-scripts)
* [Interactive Python: IPython](#interactive-python-ipython)
* [Python editors](#python-editors)
* [Python in your browser: Jupyter Notebook](#python-in-your-browser-jupyter-notebook)
* [Git](#git)
* [The PyTreat practicals](#the-pytreat-practicals)
<a class="anchor" id="python-in-a-nutshell"></a>
## Python in a nutshell
### Pros
* _Flexible_ Feel free to use functions, classes, objects, modules and
packages. Or don't - it's up to you!
* _Fast_ If you do things right (in other words, if you use `numpy`)
* _Dynamically typed_ No need to declare your variables, or specify their
types.
%% Cell type:code id: tags:
```
a = 'Great, I am _so_ sick of writing "char *a;"!'
print(a)
a = 12345
print('a is now an number!', a)
```
%% Cell type:markdown id: tags:
* _Intuitive syntax_ How do I run some code for each of the elements in my
list?
%% Cell type:code id: tags:
```
a = [1, 2, 3, 4, 5]
for elem in a:
print(elem)
```
%% Cell type:markdown id: tags:
### Cons
* _Dynamically typed_ Easier to make mistakes, harder to catch them
* _No compiler_ See above
* _Slow_ if you don't do things the right way
* _Python 2 is not the same as Python 3_ But there's an easy solution: Forget
that Python 2 exists.
* _Hard to manage different versions of python_ But we have a solution for
you: `fslpython`.
<a class="anchor" id="fslpython"></a>
## `fslpython`
FSL 5.0.10 and newer comes with its own version of Python, bundled with nearly
all of the scientific libraries that you are likely to need.
So if you use `fslpython` for all of your development, you can be sure that it
will work in FSL!
> `fslpython` is based on _conda_ and, in FSL 5.0.10, is Python version
> 3.5.2. You can read more about conda [here](https://conda.io/docs/).
> `fslpython` is based on _conda_ and, in FSL 6.0.3, is Python version
> 3.7.3. You can read more about conda [here](https://conda.io/docs/).
<a class="anchor" id="running-python-scripts"></a>
## Running Python scripts
Here's a basic Python script - a _Hello world_ for neuroimaging:
> ```
> #!/usr/bin/env fslpython
>
> # That first line up there ensures that your
> # script will be executed in the fslpython
> # environment. If you are writing a general
> # Python script, you should use this line
> # instead: #!/usr/bin/env python
>
> # In Python, we need to "import" libraries
> # (called modules) before we can use them.
> import sys
> import nibabel as nib
>
> # We can get to our command
> # line arguments via sys.argv
> fpath = sys.argv[1]
>
> # We can use nibabel to load
> # NIFTI images (and other
> # neuroimaging data formats)
> img = nib.load(fpath)
> data = img.get_data()
>
> # Now we're working with a
> # numpy array.
> nzmean = data[data != 0].mean()
>
> print('mean:', nzmean)
> ```
__Exercise__ Save the above code to a file called `script.py`, then run this
in a terminal (replace `/path/to/some/image/on/your/computer.nii.gz` with a
path to some image on your computer):
> ```
> chmod a+x script.py
> ./script.py /path/to/some/image/on/your/computer.nii.gz
> ```
<a class="anchor" id="interactive-python-ipython"></a>
## Interactive Python: IPython
Python is an [_interpreted
language_](https://en.wikipedia.org/wiki/Interpreted_language), like MATLAB. So
you can either write your code into a file, and then run that file, or you can
type code directly into a Python interpreter.
Python has a standard interpreter built-in - run `fslpython` in a terminal,
and see what happens (use CTRL+D to exit).
__But__ there is another interpreter called [`ipython`](https://ipython.org/)
which is vastly superior to the standard Python interpreter. Use `ipython`
instead! It is already installed in `fslpython`, but we just need to create a
link to it - do this now in a terminal:
> You might need to "sudo" this if your version of FSL needs admin privileges
> to modify.
>
> ln -s $FSLDIR/fslpython/envs/fslpython/bin/ipython $FSLDIR/bin/fslipython
Now if you want to do some interactive work, you can use `fslipython` in a
terminal.
__But__ there is another interpreter called [IPython](https://ipython.org/)
which is vastly superior to the standard Python interpreter. Use IPython
instead! It is already installed in `fslpython`, so if you want to do some
interactive work, you can use `fslipython` in a terminal.
__Exercise__ Do it now! Start `fslipython`, then copy/paste this code into the
prompt!
> ```
> # this line is not python - it is
> # specific to ipython/jupyter notebook
> %matplotlib
>
> import numpy as np
>
> x = np.concatenate(([0.25, 0.75], np.arange(0.1, 1.0, 0.1)))
> y = np.concatenate(([0.75, 0.75], -np.sin(np.linspace(np.pi / 4, 3 * np.pi / 4, 9))))
>
> import matplotlib.pyplot as plt
>
> fig = plt.figure()
> ax = fig.add_subplot(111)
>
> ax.scatter(x, y)
> ```
<a class="anchor" id="python-editors"></a>
## Python editors
> Summary:
> - Make your tab key insert four spaces. Don't use tab characters in Python
> code.
>
> - Use [Spyder](https://pythonhosted.org/spyder/) if you want a MATLAB-like
> envionment (focus on analysis, rather than coding).
>
> - Use [PyCharm](https://www.jetbrains.com/pycharm/) if you want an IDE-like
> environment (focus on coding, rather than analysis).
>
> - Use [Atom](https://atom.io/) or [VS Code](https://code.visualstudio.com/)
> if you like using the latest and greatest.
>
> - If you like your existing editor, use it. But you will be better off if
> you can integrate it with `fslpython`, [pylint](https://www.pylint.org/)
> and [pyflakes](https://github.com/PyCQA/pyflakes).
You can use any text editor that you want to edit Python files. But the one
golden rule that you must follow, no matter what editor you use:
__Configure your tab key to insert four spaces. Don't use tab characters in
Python code!__
This is the [standard
convention](https://www.python.org/dev/peps/pep-0008/#indentation) for Python
code. If you deviate from this convention, and somebody else needs to work
with your code, they will be angry at you!
Now, with that out of the way, there are several good Python editors available
to you. If you are getting started with Python, we recommend
[PyCharm](https://www.jetbrains.com/pycharm/) or
[Spyder](https://pythonhosted.org/spyder/).
If you are used to MATLAB, and you do a lot of experimenting or interactive
work, then you might like Spyder. If you spend most of your time writing code
rather than experimenting, then go with PyCharm.
Importantly, both PyCharm and Spyder will correctly indent your Python code!
> If you are going to stick with Emacs for your Python development, then it
> should correctly indent Python code by default. But if it isn't, add
> the following to your `~/.emacs` file:
>
> (defun my-python-mode-hook ()
> (setq indent-tabs-mode nil)
> (setq python-indent 4) ; for versions prior to 24.3
> (setq python-indent-offset 4)) ; for versions 24.3 or newer
> (add-hook 'python-mode-hook 'my-python-mode-hook)
### Spyder
Spyder is a MATLAB-like environment for Python. It has a code editor and an
interactive IPython prompt. You can inspect variables that are in your
workspace, plot data, and so on and so forth.
Beyond that, Spyder is fairly simple - it does not have much in the way of
project management tools, or integration with version control (i.e. `git`).
Spyder can be installed directly into `fslpython`:
> If your FSL installation requires administrative privileges to modify, you
> will need to prefix these commands with sudo.
>
> Install Spyder:
>
> $FSLDIR/fslpython/bin/conda install -n fslpython -y spyder
> Create a link so you can call it easily:
>
> ln -s $FSLDIR/fslpython/envs/fslpython/bin/spyder $FSLDIR/bin/fslspyder
Now to run Spyder, you can just type:
> ```
> fslspyder &
> ```
Now you need to make sure that Spyder is using the `fslpython` environment to
run your code.
1. Go to _python_ (the menu) > _Preferences_ > Python Interpreter
2. Make sure that _Use the following Python interpreter_ is selected
3. Make sure that the path is `$FSLDIR/fslpython/envs/fslpython/bin/python`
(for your specific value of `$FSLDIR`).
Type the following into the console to test that everything is working:
%% Cell type:code id: tags:
```
import matplotlib.pyplot as plt
plt.plot([1, 2, 3], [4, 5, 6])
```
%% Cell type:markdown id: tags:
### PyCharm
PyCharm is a general-purpose Python development environment. When compared to
Spyder, it is less geared towards interactive analysis, but has better code
editing tools (e.g. autocomplete and refactoring), and better file
management/version control integration.
And it is also easy to install - simply download the Community edition from
the [home page](https://www.jetbrains.com/pycharm/). Then, if you are on a
mac, double-click the `.dmg` file, and drag PyCharm into your `/Applications/`
folder. Then double-click on PyCharm to start it.
Once you have chosen a theme, you will be asked if you would like to create a
_Launcher script_ - do this, because you will then be able to open files from
the terminal by typing `charm [filename]`.
Now you will be presented with a _Welcome to PyCharm_ window. __Before doing
anything else__, click on the _Configure_ button down in the bottom right, and
choose _Preferences_. Then in the _Project Interpreter_ section:
1. Click on the gear button and choose _Add local..._
2. Choose _Existing environment_
3. Click on the _..._ button, and navigate to
`$FSLDIR/fslpython/envs/fslpython/bin/python`.
### Other options
Emacs is capable of being used as a fully-fledged Python IDE, if you have the
time and patience to configure it.
If you are the fashionable sort, try one of these:
- [Visual Studio Code](https://code.visualstudio.com/)
- [Atom](https://atom.io/)
<a class="anchor" id="python-in-your-browser-jupyter-notebook"></a>
## Python in your browser: Jupyter Notebook
It is possible to do your Python-based development and experimentation inside
your web browser, thanks to the [Jupyter project](https://jupyter.org/).
Jupyter works in the following way:
1. You start a Juptyer web server on your computer. This web server provides
the environment in which your Python code is executed.
2. You open https://localhost:8888 (or similar) in a web browser - this opens
a connection to the (locally running) web server.
3. You start a "Notebook" (Jupyter's version of a file), and start typing.
You can put text, LaTeX, and of course Python code into a notebook.
4. All of the code that you write gets sent to the web server and
executed. Then the results get sent back to the web browser, and displayed
in your notebook - magic!
All of the PyTreat practicals are written in a Jupyter notebook. Some of the
talks are too - you're looking at a Jupyter Notebook right now!
So you're going to need to install Jupyter to get the most out of the
practicals that we have prepared for you.
__Exercise__ Install Jupyter into `fslpython` - run these commands in a
terminal:
> Remember to prefix with sudo if your FSL install needs admin to modify.
>
> $FSLDIR/fslpython/bin/conda install -n fslpython -y jupyter
> And add a link so you can call it easily:
>
> ln -s $FSLDIR/fslpython/envs/fslpython/bin/jupyter $FSLDIR/bin/fsljupyter
<a class="anchor" id="git"></a>
## Git
All the cool kids these days use [git](https://git-scm.com/) to
collaboratively work on their Python code. The PyTreat is a great opportunity
to start learning and using it for your own work!
Git is different from CVS and SVN in that it is _distributed_. In CVS and SVN,
there is only one central repository. You check out a copy of the source from
the repository, make some changes, and then commit those changes back to the
central repository.
In git, there are multiple repositories. There is usually one repository which
acts as the central one, but you will _clone_ (or _fork_) that central
repository to create your own full copy of the repository.
Then, you can make changes and commit them in your own repository. And at any
point, you can push your changes back to the central repository.
### WIN's gitlab server
https://git.fmrib.ox.ac.uk/
If you have a FMRIB account, then you also have a Gitlab account. Gitlab is a
git server that you can use to store your "central" git repository. Gitlab is
very similar to https://www.github.com, but it is managed by WIN, and your
code is not publicly visible (although it can be if you want). Gitlab backs
up your code automatically, and has a nice web interface.
You can have up to 10 projects on your gitlab account - talk to the FMRIB IT
people if you need more.
### Using git and gitlab
> We need to go through a couple of intiial configuration steps before
> proceeding. You only need to do this once (for each computer that you use).
>
> First, run these commands to configure git on your system.
>
> ```
> git config --global user.name "Your name"
> git config --global user.email "Your email address"
> ```
>
> Now you need to create a SSH key pair, so your computer can talk to the
> gitlab server without you having to log in. Don't be scared - there are
> detailed instructions on doing this at
> https://git.fmrib.ox.ac.uk/help/ssh/README.md - follow the instructions
> under the section entitled __Generating a new SSH key pair__.
Working with git and gitlab generally involves these steps:
1. Add new files to your local repository, or make changes to existing files.
2. Run `git add` on the new/changed files to _stage_ them.
3. Run `git commit` to commit all staged changes to your local repository.
4. Run `git push` to push those commits to your gitlab repository.
When you start working on a new project (or if you have an existing project
that you want to put into git):
__1. Organise all of your project files into their own folder__
This sounds obvious, but just to be sure.
__2. Create a repository for your project on gitlab__
Log in to gitlab (https://git.fmrib.ox.ac.uk/), then click on the _+_ button
towards the top right, and select _New project_. Give the project a name and
choose its visiblity (note that _Public_ means your project will be visible to
the world).
__3. Turn your project folder into a git repository__
Now, follow the instructions that are listd on your gitlab project home page,
under __Existing folder__ (repeated here):
> cd existing_folder
> git init
> git remote add origin git@git.fmrib.ox.ac.uk:username/project.git
> git add .
> git commit -m "Initial commit"
> git push -u origin master
The `git add .` line will add _all_ of the files in your project directory
into git. If you only want certain files in git, then `git add` them one by
one (or use standard bash file patterns, e.g. `git add *.py`).
You should avoid putting large binary files or data files into git - it works
best with plain text. Talk to the FMRIB IT people if you really need to store
large files in git, as they can help you with this.
__4. Develop your super cool project!__
Now you can get to work! Whenever you make changes to your code that you want
saved, follow the `git add`, `git commit`, `git push` steps described above.
For example, let's say we have added a new file called `cool_module.py`. To
get this into git, we would do the following:
%% Cell type:code id: tags:
```
git add cool_module.py
git commit -m "Added cool module. It's super cool"
git push origin master
```
%% Cell type:markdown id: tags:
<a class="anchor" id="the-pytreat-practicals"></a>
## The PyTreat practicals
All of the practicals for PyTreat are hosted in a git repository at:
https://git.fmrib.ox.ac.uk/fsl/pytreat-2018-practicals
https://git.fmrib.ox.ac.uk/fsl/pytreat-practicals-2020
So let's go get them!
> ```
> git clone git@git.fmrib.ox.ac.uk:fsl/pytreat-2018-practicals.git
> cd pytreat-2018-practicals
> fsljupyter notebook
> git clone https://git.fmrib.ox.ac.uk/fsl/pytreat-practicals-2020.git
> cd pytreat-practicals-2020
> fslpython -m notebook
> ```
......
talks/introduction/pytreat_intro.md
View file @ 15d75bb5
# Welcome to the WIN PyTreat 2018!
# Welcome to the WIN PyTreat 2020!
Program: https://docs.google.com/document/d/10CwLEhUi-YiwfC2F40QCVm6eEVwKiaXkfTKz67xWAfM/edit?usp=sharing
Program: https://docs.google.com/document/d/11LQgxC-LZPG_TXS3MP9tYXNWABAQdTI00iavk-tKttU/edit
__We need to do some setting up, so get your laptop ready!__
__Make sure you have FSL 5.0.10 installed and working!__
__Make sure you have FSL 6.0.3 installed and working!__
__End all sentences with an exclamation mark!__
__Open this page in your web browser!__
https://git.fmrib.ox.ac.uk/fsl/pytreat-2018-practicals/tree/master/talks/introduction/pytreat_intro.md
https://git.fmrib.ox.ac.uk/fsl/pytreat-practicals-2020/tree/master/talks/introduction/pytreat_intro.md
## Overview
......@@ -92,8 +92,8 @@ So if you use `fslpython` for all of your development, you can be sure that it
will work in FSL!
> `fslpython` is based on _conda_ and, in FSL 5.0.10, is Python version
> 3.5.2. You can read more about conda [here](https://conda.io/docs/).
> `fslpython` is based on _conda_ and, in FSL 6.0.3, is Python version
> 3.7.3. You can read more about conda [here](https://conda.io/docs/).
<a class="anchor" id="running-python-scripts"></a>
......@@ -160,20 +160,10 @@ Python has a standard interpreter built-in - run `fslpython` in a terminal,
and see what happens (use CTRL+D to exit).
__But__ there is another interpreter called [`ipython`](https://ipython.org/)
which is vastly superior to the standard Python interpreter. Use `ipython`
instead! It is already installed in `fslpython`, but we just need to create a
link to it - do this now in a terminal:
> You might need to "sudo" this if your version of FSL needs admin privileges
> to modify.
>
> ln -s $FSLDIR/fslpython/envs/fslpython/bin/ipython $FSLDIR/bin/fslipython
Now if you want to do some interactive work, you can use `fslipython` in a
terminal.
__But__ there is another interpreter called [IPython](https://ipython.org/)
which is vastly superior to the standard Python interpreter. Use IPython
instead! It is already installed in `fslpython`, so if you want to do some
interactive work, you can use `fslipython` in a terminal.
__Exercise__ Do it now! Start `fslipython`, then copy/paste this code into the
......@@ -384,24 +374,6 @@ All of the PyTreat practicals are written in a Jupyter notebook. Some of the
talks are too - you're looking at a Jupyter Notebook right now!
So you're going to need to install Jupyter to get the most out of the
practicals that we have prepared for you.
__Exercise__ Install Jupyter into `fslpython` - run these commands in a
terminal:
> Remember to prefix with sudo if your FSL install needs admin to modify.
>
> $FSLDIR/fslpython/bin/conda install -n fslpython -y jupyter
> And add a link so you can call it easily:
>
> ln -s $FSLDIR/fslpython/envs/fslpython/bin/jupyter $FSLDIR/bin/fsljupyter
<a class="anchor" id="git"></a>
## Git
......@@ -538,14 +510,14 @@ git push origin master
All of the practicals for PyTreat are hosted in a git repository at:
https://git.fmrib.ox.ac.uk/fsl/pytreat-2018-practicals
https://git.fmrib.ox.ac.uk/fsl/pytreat-practicals-2020
So let's go get them!
> ```
> git clone git@git.fmrib.ox.ac.uk:fsl/pytreat-2018-practicals.git
> cd pytreat-2018-practicals
> fsljupyter notebook
> git clone https://git.fmrib.ox.ac.uk/fsl/pytreat-practicals-2020.git
> cd pytreat-practicals-2020
> fslpython -m notebook
> ```
talks/matlab_vs_python/bloch/bloch.ipynb 0 → 100644
View file @ 15d75bb5
%% Cell type:markdown id: tags:
Imports
%% Cell type:code id: tags:
``` python
import numpy as np
from scipy.integrate import ode, solve_ivp
import matplotlib.pyplot as plt
```
%% Cell type:markdown id: tags:
# Define the Bloch equation
$$\frac{d\vec{M}}{dt} = \vec{M}\times \vec{B} - \frac{M_x + M_y}{T2} - \frac{M_z - M_0}{T1}$$
In this section:
- define a function
- numpy functions like np.cross
%% Cell type:code id: tags:
``` python
# define bloch equation
def bloch(t, M, T1, T2):
# get effective B-field at time t
B = B_eff(t)
# cross product of M and B, add T1 and T2 relaxation terms
return np.array([M[1]*B[2] - M[2]*B[1] - M[0]/T2,
M[2]*B[0] - M[0]*B[2] - M[1]/T2,
M[0]*B[1] - M[1]*B[0] - (M[2]-1)/T1])
# alternatively
#return np.cross(M,B) - np.array([M[0]/T2, M[1]/T2, (M[2]-1)/T1])
```
%% Cell type:markdown id: tags:
# Define the pulse sequence
We work in the rotating frame, so we only need the amplitude envelope of the RF pulse
Typically, B1 excitation fields point in the x- and/or y-directions
Gradient fields point in the z-direction
In this simple example, a simple sinc-pulse excitation pulse is applied for 1 ms along the x-axis
then a gradient is turned on for 1.5 ms after that.
In this section:
- constants such as np.pi
- functions like np.sinc
%% Cell type:code id: tags:
``` python
# define effective B-field
def B_eff(t):
# Do nothing for 0.25 ms
if t < 0.25:
return np.array([0, 0, 0])
# Sinc RF along x-axis and slice-select gradient on for 1.00 ms
elif t < 1.25:
return np.array([np.pi*np.sinc((t-0.75)*4), 0, np.pi])
# Do nothing for 0.25 ms
elif t < 1.50:
return np.array([0, 0, 0])
# Slice refocusing gradient on for 1.50 ms
# Half the area of the slice-select gradient lobe
elif t < 3.00:
return np.array([0, 0, -(1/3)*np.pi])
# Pulse sequence finished
else:
return np.array([0, 0, 0])
```
%% Cell type:markdown id: tags:
# Plot the pulse sequence
In this section:
- unpacking return values
- unwanted return values
- list comprehension
- basic plotting
%% Cell type:code id: tags:
``` python
# Create 2 vertical subplots
_, ax = plt.subplots(2, 1, figsize=(12,12))
# Get pulse sequence B-fields from 0 - 5 ms
pulseq = [B_eff(t) for t in np.linspace(0,5,1000)]
pulseq = np.array(pulseq)
# Plot RF
ax[0].plot(pulseq[:,0])
ax[0].set_ylabel('B1')
# Plot gradient
ax[1].plot(pulseq[:,2])
ax[1].set_ylabel('Gradient')
```
%% Cell type:markdown id: tags:
# Integrate ODE
This uses a Runge-Kutta variant called the "Dormand-Prince method" by default. Other ode integration methods are available.
In this section:
- ode solvers
- lambdas (anonymous functions)
%% Cell type:code id: tags:
``` python
# Set the initial conditions
# equilibrium magnetization (M) = (0, 0, 1)
M = [0, 0, 1]
# Set time interval for integration
t = [0, 5]
# Set max step size
dt = 0.005
# Set T1 and T2
T1, T2 = 1500, 50
# Integrate ODE
# In Scipy 1.2.0, the first argument to solve_ivp must be a function that takes exactly 2 arguments
sol = solve_ivp(lambda t, M : bloch(t, M, T1, T2), t, M, max_step=dt)
# Grab output
t = sol.t
M = sol.y
```
%% Cell type:markdown id: tags:
# Plot Results
In this section:
- more plotting
%% Cell type:code id: tags:
``` python
# Create single axis
_, ax = plt.subplots(figsize=(12,12))
# Plot x, y and z components of Magnetisation
ax.plot(t, M[0,:], label='Mx')
ax.plot(t, M[1,:], label='My')
ax.plot(t, M[2,:], label='Mz')
# Add legend and grid
ax.legend()
ax.grid()
```
talks/matlab_vs_python/bloch/bloch.m 0 → 100644
View file @ 15d75bb5
% Plot the pulse sequence
% create figure
figure();
% get pulse sequence B-fields from 0-5 ms
pulseq = zeros(1000,3);
for i = 1:1000
pulseq(i,:) = B_eff(i*0.005);
end
% plot RF
subplot(2,1,1);
plot(pulseq(:,1));
ylabel('B1');
% plot gradient
subplot(2,1,2);
plot(pulseq(:,3));
ylabel('Gradient');
%% Integrate ODE
T1 = 1500;
T2 = 50;
t0 = 0;
t1 = 5;
dt = 0.005;
M0 = [0; 0; 1];
[t, M] = ode45(@(t,M)bloch_ode(t, M, T1, T2), linspace(t0, t1, (t1-t0)/dt), M0);
%% Plot Results
% create figure
figure();hold on;
% plot x, y and z components of Magnetisation
plot(t, M(:,1), 'linewidth', 2);
plot(t, M(:,2), 'linewidth', 2);
plot(t, M(:,3), 'linewidth', 2);
% add legend and grid
legend({'Mx','My','Mz'});
grid on;
%% define the bloch equation
function dM = bloch_ode(t, M, T1, T2)
% get effective B-field at time t
B = B_eff(t);
% cross product of M and B, add T1 and T2 relaxation terms
dM = [M(2)*B(3) - M(3)*B(2) - M(1)/T2;
M(3)*B(1) - M(1)*B(3) - M(2)/T2;
M(1)*B(2) - M(2)*B(1) - (M(3)-1)/T1];
end
%% define effective B-field
function b = B_eff(t)
% Do nothing for 0.25 ms
if t < 0.25
b = [0, 0, 0];
% Sinc RF along x-axis and slice-select gradient on for 1.00 ms
elseif t < 1.25
b = [pi*sinc((t-0.75)*4), 0, pi];
% Do nothing for 0.25 ms
elseif t < 1.50
b = [0, 0, 0];
% Slice refocusing gradient on for 1.50 ms
% Half the area of the slice-select gradient lobe
elseif t < 3.00
b = [0, 0, -(1/3)*pi];
% pulse sequence finished
else
b = [0, 0, 0];
end
end
talks/matlab_vs_python/fitting/fit_model.ipynb 0 → 100644
View file @ 15d75bb5
%% Cell type:markdown id: tags:
## Short example of Model fitting
Here we fit a MRI-style model to some simulated data.
The model is: $\textrm{Signal} = M_0\exp\left[-R_2\textrm{TE}\right]\left(1-\exp\left[-R_1\textrm{TI}\right]\right)$
The parameters that we will be fitting are $(M_0,R_1,R_2)$.
Basic imports:
%% Cell type:code id: tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
```
%% Cell type:markdown id: tags:
In this section:
- defining a numpy array
- double list comprehension
%% Cell type:code id: tags:
``` python
TEs = np.array([10,40,50,60,80]) # TE values in ms
TRs = np.array([.8,1,1.5,2]) # TR in seconds (I know this is bad)
# All combinations of TEs/TRs
comb = np.array([(x,y) for x in TEs for y in TRs])
TEs,TRs = comb[:,0],comb[:,1]
```
%% Cell type:markdown id: tags:
Now we define our forward model
In this section:
- inline function definition
- random number generation
%% Cell type:code id: tags:
``` python
# function for our model
def forward(p):
M0,R1,R2 = p
return M0*np.exp(-R2*TEs)*(1-np.exp(-R1*TRs))
# simulate data using model
true_p = [100,1/.8,1/50] # M0, R1=1/T1,R2=1/T2
data = forward(true_p)
snr = 50
noise_std = true_p[0]/snr
noise = np.random.randn(data.size)*noise_std
data = data + noise
# quick plot of the data
plt.figure()
plt.plot(data)
```
%% Output
[<matplotlib.lines.Line2D at 0x121b13978>]
%% Cell type:markdown id: tags:
Now we have the data and our forward model we are almost ready to begin fitting.
We need a cost function to optimise. We will use mean squared error.
In this section:
- '**' operation
- np.mean
%% Cell type:code id: tags:
``` python
# cost function is mean square error divided by 2
def cf(p):
pred = forward(p)
return np.mean((pred-data)**2)/2.0
```
%% Cell type:markdown id: tags:
Now we can set up our optimisation.
In this section:
- scipy minimize
- dictionary
- keyword arguments
%% Cell type:code id: tags:
``` python
# get ready to minimize
p0 = [200,1/1,1/70] # some random initial guess
method = 'powell' # pick a method. scipy has loads!
kw_args = {'x0':p0,'method':method}
result = minimize(cf,**kw_args)
```
%% Cell type:markdown id: tags:
Plot the data with the model prediction.
In this section
- printing
- text formatting
%% Cell type:code id: tags:
``` python
plt.figure()
plt.plot(data,'o')
plt.plot(forward(result.x))
print('fitted = {}'.format(result.x))
print('true = {}'.format(true_p))
```
%% Output
fitted = [1.02562035e+02 1.16194491e+00 1.98071179e-02]
true = [100, 1.25, 0.02]
%% Cell type:markdown id: tags:
## Optional: use gradients and hessians to help with the optimisation
In this example the forward model is simple enough that calculating the derivatives of the cost function is relatively easy to do analytically. Below is an example of how you could define these and use them in the fitting
%% Cell type:code id: tags:
``` python
# gradient of the forward model
def forward_deriv(p):
M0,R1,R2 = p
E1,E2 = np.exp(-R1*TRs),np.exp(-R2*TEs)
dE1 = -TRs*E1
dE2 = -TEs*E2
# f = M0*E2*(1-E1)
dfdM0 = E2*(1-E1)
dfdR1 = M0*E2*(-dE1)
dfdR2 = M0*dE2*(1-E1)
return np.array([dfdM0,dfdR1,dfdR2])
# hessian of the forward model
def forward_deriv2(p):
M0,R1,R2 = p
E1,E2 = np.exp(-R1*TRs),np.exp(-R2*TEs)
dE1 = -TRs*E1
dE2 = -TEs*E2
ddE1 = (TRs**2)*E1
ddE2 = (TEs**2)*E2
dfdM0dM0 = np.zeros(E1.shape)
dfdM0dR1 = E2*(-dE1)
dfdM0dR2 = dE2*(1-E1)
dfdR1dM0 = E2*(-dE1)
dfdR1dR1 = M0*E2*(-ddE1)
dfdR1dR2 = M0*(dE2)*(-dE1)
dfdR2dM0 = dE2*(1-E1)
dfdR2dR1 = M0*dE2*(-dE1)
dfdR2dR2 = M0*ddE2*(1-E1)
return np.array([[dfdM0dM0,dfdM0dR1,dfdM0dR2],
[dfdR1dM0,dfdR1dR1,dfdR1dR2],
[dfdR2dM0,dfdR2dR1,dfdR2dR2]])
# cost function is mean square error divided by 2
def cf(p):
pred = forward(p)
return np.mean((pred-data)**2)/2.0
def cf_grad(p):
pred = forward(p)
deriv = forward_deriv(p)
return np.mean( deriv * (pred-data)[None,:],axis=1)
def cf_hess(p):
pred = forward(p)
deriv = forward_deriv(p)
deriv2 = forward_deriv2(p)
H = np.zeros((len(p),len(p)))
for i in range(H.shape[0]):
for j in range(H.shape[1]):
H[i,j] = np.mean(deriv2[i,j]*(pred-data) + deriv[i]*deriv[j])
return H
```
%% Cell type:code id: tags:
``` python
# get ready to minimize
p0 = [200,1/1,1/70] # some random guess
method = 'trust-ncg'
kw_args = {'x0':p0,'method':method,'jac':cf_grad,'hess':cf_hess}
result = minimize(cf,**kw_args)
plt.figure()
plt.plot(data,'o')
plt.plot(forward(result.x))
print('fitted = {}'.format(result.x))
print('true = {}'.format(true_p))
```
%% Output
fitted = [1.02576306e+02 1.16153921e+00 1.98044006e-02]
true = [100, 1.25, 0.02]
%% Cell type:markdown id: tags:
## The end.
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
talks/matlab_vs_python/fitting/fit_model.m 0 → 100644
View file @ 15d75bb5
%% Play with model fitting
% experimental parameters
TEs = [10 40 50 60 80];
TRs = [.8 1 1.5 2];
[TEs,TRs] = meshgrid(TEs,TRs);
TEs = TEs(:)'; TRs = TRs(:)';
% forward model
forward = @(p)( p(1)*exp(-p(3)*TEs).*(1-exp(-p(2)*TRs)));
% simulate data
true_p = [100,1/.8,1/50];
data = forward(true_p);
snr = 50;
noise_std = 100/snr;
noise = randn(size(data))*noise_std;
data = data+noise;
plot(data)
%%
% cost function is mean squared error (MSE)
cf = @(x)( mean( (forward(x)-data).^2 ) );
% initial guess
p0 = [200,1/1,1/70];
% using fminsearch (Nelder-Mead)
p = fminsearch(@(x) cf(x),p0);
% plot result
figure,hold on
plot(data,'.')
plot(forward(p))
%% The below uses fminunc, which allows morre flexibility
% (like choosing the algorithm or providing gradients and Hessians)
options = optimoptions('fminunc','Display','off','Algorithm','quasi-newton');
[x,fval] = fminunc(cf,p0,options);
figure,hold on
plot(data,'.')
plot(forward(x))
talks/matlab_vs_python/migp/MIGP.ipynb 0 → 100644
View file @ 15d75bb5
%% Cell type:markdown id: tags:
# MIGP
For group ICA, `melodic` uses multi-session temporal concatenation. This will perform a single 2D ICA run on the concatenated data matrix (obtained by stacking all 2D data matrices of every single data set on top of each other).
![temporal concatenation](concat_diag.png)
Resulting in **high dimension** datasets!
Furthermore, with ICA we are typically only interested in a comparitively low dimension decomposition so that we can capture spatially extended networks.
Therefore the first step is to reduce the dimensionality of the data. This can be achieved in a number of ways, but `melodic`, by default, uses `MIGP`.
> MIGP is an incremental approach that aims to provide a very close approximation to full temporal concatenation followed by PCA, but without the large memory requirements *(Smith et al., 2014)*.
Essentially, MIGP stacks the datasets incrementally in the temporal dimension, and whenever the temporal dimension exceeds a specified size, a PCA-based temporal reduction is performed.
> MIGP does not increase at all in memory requirement with increasing numbers of subjects, no large matrices are ever formed, and the computation time scales linearly with the number of subjects. It is easily parallelisable, simply by applying the approach in parallel to subsets of subjects, and then combining across these with the same “concatenate and reduce” approach described above *(Smith et al., 2014)*.
## This notebook
This notebook will download an open fMRI dataset (~50MB) for use in the MIGP demo, regresses confounds from the data, performs spatial smoothing with 10mm FWHM, and then runs group `melodic` with `MIGP`.
* [Fetch the data](#download-the-data)
* [Clean the data](#clean-the-data)
* [Run `melodic`](#run-melodic)
* [Plot group ICs](#plot-group-ics)
Firstly we will import the necessary packages for this notebook:
%% Cell type:code id: tags:
``` python
from nilearn import datasets
from nilearn import image
from nilearn import plotting
import nibabel as nb
import numpy as np
import os.path as op
import os
import glob
import matplotlib.pyplot as plt
```
%% Cell type:markdown id: tags:
<a class="anchor" id="download-the-data"></a>
## Fetch the data
This data is a derivative from the [COBRE](http://fcon_1000.projects.nitrc.org/indi/retro/cobre.html) sample found in the International Neuroimaging Data-sharing Initiative, originally released under Creative Commons - Attribution Non-Commercial.
It comprises 10 preprocessed resting-state fMRI selected from 72 patients diagnosed with schizophrenia and 74 healthy controls (6mm isotropic, TR=2s, 150 volumes).
Create a directory in the users home directory to store the downloaded data:
> **NOTE:** [`expanduser`](https://docs.python.org/3.7/library/os.path.html#os.path.expanduser) will expand the `~` to the be users home directory:
%% Cell type:code id: tags:
``` python
data_dir = op.expanduser('~/nilearn_data')
if not op.exists(data_dir):
os.makedirs(data_dir)
```
%% Cell type:markdown id: tags:
Download the data (if not already downloaded):
> **Note:** We use a method from [`nilearn`](https://nilearn.github.io/index.html) called [`fetch_cobre`](https://nilearn.github.io/modules/generated/nilearn.datasets.fetch_cobre.html) to download the fMRI data
%% Cell type:code id: tags:
``` python
d = datasets.fetch_cobre(data_dir=data_dir)
```
%% Cell type:markdown id: tags:
<a class="anchor" id="clean-the-data"></a>
## Clean the data
Regress confounds from the data and to spatially smooth the data with a gaussian filter of 10mm FWHM.
> **Note:**
> 1. We use [`clean_img`](https://nilearn.github.io/modules/generated/nilearn.image.clean_img.html) from the [`nilearn`](https://nilearn.github.io/index.html) package to regress confounds from the data
> 2. We use [`smooth_img`](https://nilearn.github.io/modules/generated/nilearn.image.smooth_img.html) from the [`nilearn`](https://nilearn.github.io/index.html) package to spatially smooth the data
> 3. [`zip`](https://docs.python.org/3.3/library/functions.html#zip) takes iterables and aggregates them in a tuple. Here it is used to iterate through four lists simultaneously
> 4. We use list comprehension to loop through all the filenames and append suffixes
%% Cell type:code id: tags:
``` python
# Create a list of filenames for cleaned and smoothed data
clean = [f.replace('.nii.gz', '_clean.nii.gz') for f in d.func]
smooth = [f.replace('.nii.gz', '_clean_smooth.nii.gz') for f in d.func]
# loop through each subject, regress confounds and smooth
for img, cleaned, smoothed, conf in zip(d.func, clean, smooth, d.confounds):
print(f'{img}: regress confounds: ', end='')
image.clean_img(img, confounds=conf).to_filename(cleaned)
print(f'smooth.')
image.smooth_img(img, 10).to_filename(smoothed)
```
%% Cell type:markdown id: tags:
To run ```melodic``` we will need a brain mask in MNI152 space at the same resolution as the fMRI.
> **Note:**
> 1. We use [`load_mni152_brain_mask`](https://nilearn.github.io/modules/generated/nilearn.datasets.load_mni152_brain_mask.html) from the [`nilearn`](https://nilearn.github.io/index.html) package to load the MNI152 mask
> 2. We use [`resample_to_img`](https://nilearn.github.io/modules/generated/nilearn.image.resample_to_img.html) from the [`nilearn`](https://nilearn.github.io/index.html) package to resample the mask to the resolution of the fMRI
> 3. We use [`math_img`](https://nilearn.github.io/modules/generated/nilearn.image.math_img.html) from the [`nilearn`](https://nilearn.github.io/index.html) package to binarize the resample mask
> 4. The mask is plotted using [`plot_anat`](https://nilearn.github.io/modules/generated/nilearn.plotting.plot_anat.html) from the [`nilearn`](https://nilearn.github.io/index.html) package
%% Cell type:code id: tags:
``` python
# load a single fMRI dataset (func0)
func0 = nb.load(d.func[0].replace('.nii.gz', '_clean_smooth.nii.gz'))
# load MNI153 brainmask, resample to func0 resolution, binarize, and save to nifti
mask = datasets.load_mni152_brain_mask()
mask = image.resample_to_img(mask, func0)
mask = image.math_img('img > 0.5', img=mask)
mask.to_filename(op.join(data_dir, 'brain_mask.nii.gz'))
# plot brainmask to make sure it looks OK
disp = plotting.plot_anat(image.index_img(func0, 0))
disp.add_contours(mask, threshold=0.5)
```
%% Cell type:markdown id: tags:
<a class="anchor" id="run-melodic"></a>
### Run ```melodic```
Generate a command line string and run group ```melodic``` on the smoothed fMRI with a dimension of 10 components:
> **Note**:
> 1. Here we use python [f-strings](https://www.python.org/dev/peps/pep-0498/), formally known as literal string interpolation, which allow for easy formatting
> 2. [`op.join`](https://docs.python.org/3.7/library/os.path.html#os.path.join) will join path strings using the platform-specific directory separator
> 3. [`','.join(smooth)`](https://docs.python.org/3/library/stdtypes.html#str.join) will create a comma seprated string of all the items in the list `smooth`
%% Cell type:code id: tags:
``` python
# generate melodic command line string
melodic_cmd = f"melodic -i {','.join(smooth)} --mask={op.join(data_dir, 'brain_mask.nii.gz')} -d 10 -v -o cobre.gica "
print(melodic_cmd)
```
%% Cell type:markdown id: tags:
> **Note:**
> 1. Here we use the `!` operator to execute the command in the shell
> 2. The `{}` will expand the contained python variable in the shell
%% Cell type:code id: tags:
``` python
# run melodic
! {melodic_cmd}
```
%% Cell type:markdown id: tags:
<a class="anchor" id="plot-group-ics"></a>
### Plot group ICs
Now we can load and plot the group ICs generated by ```melodic```.
This function will be used to plot ICs:
> **NOTE:**
> 1. Here we use [`plot_stat_map`](https://nilearn.github.io/modules/generated/nilearn.plotting.plot_stat_map.html) from the [`nilearn`](https://nilearn.github.io/index.html) package to plot the orthographic images
> 2. [`subplots`](https://matplotlib.org/api/_as_gen/matplotlib.pyplot.subplots.html) from `matplotlib.pyplot` creates a figure and multiple subplots
> 3. [`find_xyz_cut_coords`](https://nilearn.github.io/modules/generated/nilearn.plotting.find_xyz_cut_coords.html) from the [`nilearn`](https://nilearn.github.io/index.html) package will find the image coordinates of the center of the largest activation connected component
> 4. [`zip`](https://docs.python.org/3.3/library/functions.html#zip) takes iterables and aggregates them in a tuple. Here it is used to iterate through two lists simultaneously
> 5. [`iter_img`](https://nilearn.github.io/modules/generated/nilearn.image.iter_img.html) from the [`nilearn`](https://nilearn.github.io/index.html) package creates an iterator from an image that steps through each volume/time-point of the image
%% Cell type:code id: tags:
``` python
def map_plot(d):
N = d.shape[-1]
fig, ax = plt.subplots(int(np.ceil((N/2))),2, figsize=(12, N))
for img, ax0 in zip(image.iter_img(d), ax.ravel()):
coord = plotting.find_xyz_cut_coords(img, activation_threshold=3.5)
plotting.plot_stat_map(img, cut_coords=coord, vmax=10, axes=ax0)
return fig
```
%% Cell type:markdown id: tags:
Hopefully you can see some familiar looking RSN spatial patterns:
%% Cell type:code id: tags:
``` python
# Load ICs
ics = nb.load('cobre.gica/melodic_IC.nii.gz')
# plot
fig = map_plot(ics)
```
talks/matlab_vs_python/migp/concat_diag.png 0 → 100644
View file @ 15d75bb5
talks/matlab_vs_python/migp/concat_diag.png

27.1 KiB

talks/matlab_vs_python/migp/demo_MIGP_NIFTI.m 0 → 100644
View file @ 15d75bb5
%% MIGP - MELODIC's Incremental Group-PCA
% Steve Smith, FMRIB, 2012-2013
% not yet finally tweaked/evaluated, or published - please do not pass on.
% adapted by Sean Fitz, 2020
%% User Options (equivalent to melodic defaults)
INlist=dir('~/nilearn_data/cobre/fmri_*_smooth.nii.gz'); % list of input 4D NIFTI standard space datasets
INmask='~/nilearn_data/brain_mask.nii.gz'; % 3D NIFTI mask image
GO='matMIGP'; % output filename
dPCAint=299; % internal number of components - typically 2-4 times number of timepoints in each run (if you have enough RAM for that)
dPCAout=299; % number of eigenvectors to output - should be less than dPCAint and more than the final ICA dimensionality
sep_vn=false; % switch on separate variance nomalisation for each input dataset
%% Run MIGP
Nsub=length(INlist); [~,r]=sort(rand(Nsub,1)); % will process subjects in random order
mask=read_avw(INmask);
masksize=size(mask);
mask=reshape(mask,prod(masksize),1);
demean = @(x) x - repmat(mean(x,1),[size(x,1), 1, 1]);
for i = 1:Nsub
% read data
filename=[INlist(r(i)).folder, '/', INlist(r(i)).name];
fprintf('Reading data file %s\n', filename);
grot=double(read_avw(filename));
grot=reshape(grot,prod(masksize),size(grot,4));
% demean
fprintf('\tRemoving mean image\n');
grot=demean(grot(mask~=0,:)');
% var-norm separately for each subject
if sep_vn==true
fprintf('\tNormalising by voxel-wise variance\r');
[uu,ss,vv]=svds(grot,30);
vv(abs(vv)<2.3*std(vv(:)))=0;
stddevs=max(std(grot-uu*ss*vv'),0.001);
grot=grot./repmat(stddevs,size(grot,1),1);
end
if (i==1)
W=demean(grot); clear grot;
else
% concat
W=[W; demean(grot)]; clear grot;
% reduce W to dPCAint eigenvectors
if size(W,1)-10 > dPCAint
fprintf('\tReducing data matrix to a %i dimensional subspace\n', dPCAint);
[uu,dd]=eigs(W*W',dPCAint);
W=uu'*W;
clear uu;
end
end
end
%% Reshape and save
grot=zeros(prod(masksize),dPCAout);
grot(mask~=0,:)=W(1:dPCAout,:)'; grot=reshape(grot,[masksize ,dPCAout]);
fprintf('Save to %s.nii.gz\n', GO)
save_avw(grot,GO,'f',[1 1 1 1]);
system(sprintf('fslcpgeom %s %s -d',filename,GO));
talks/matlab_vs_python/migp/matlab_MIGP.ipynb 0 → 100644
View file @ 15d75bb5
%% Cell type:markdown id: tags:
## Matlab MIGP
This notebook will load the dimension reduced data from Matlab MIGP, run group ICA, and then plot the group ICs.
* [Run `melodic`](#run-matlab-melodic)
* [Plot group ICs](#plot-matlab-group-ics)
Firstly we will import the necessary packages for this notebook:
%% Cell type:code id: tags:
``` python
from nilearn import plotting
from nilearn import image
import nibabel as nb
import matplotlib.pyplot as plt
import numpy as np
import os.path as op
```
%% Cell type:markdown id: tags:
It will be necessary to know the location where the data was stored so that we can load the brainmask:
> **Note**: [`expanduser`](https://docs.python.org/3.7/library/os.path.html#os.path.expanduser) will expand the `~` to the be users home directory
%% Cell type:code id: tags:
``` python
data_dir = op.expanduser('~/nilearn_data')
```
%% Cell type:markdown id: tags:
<a class="anchor" id="run-matlab-melodic"></a>
### Run ```melodic```
Generate a command line string and run group ```melodic``` on the Matlab MIGP dimension reduced data with a dimension of 10 components. We disable MIGP because it was already run separately in Matlab.
> **Note**:
> 1. Here we use python [f-strings](https://www.python.org/dev/peps/pep-0498/), formally known as literal string interpolation, which allow for easy formatting
> 2. [`op.join`](https://docs.python.org/3.7/library/os.path.html#os.path.join) will join path strings using the platform-specific directory separator
%% Cell type:code id: tags:
``` python
# generate melodic command line string
melodic_cmd = f"melodic -i matMIGP.nii.gz --mask={op.join(data_dir, 'brain_mask.nii.gz')} -d 10 -v --nobet --disableMigp -o matmigp.gica"
print(melodic_cmd)
```
%% Cell type:markdown id: tags:
> **Note:**
> 1. Here we use the `!` operator to execute the command in the shell
> 2. The `{}` will expand the contained python variable in the shell
%% Cell type:code id: tags:
``` python
# run melodic
! {melodic_cmd}
```
%% Cell type:markdown id: tags:
<a class="anchor" id="plot-matlab-group-ics"></a>
### Plot group ICs
Now we can load and plot the group ICs generated by ```melodic```.
This function will be used to plot ICs:
> **NOTE:**
> 1. Here we use [`plot_stat_map`](https://nilearn.github.io/modules/generated/nilearn.plotting.plot_stat_map.html) from the [`nilearn`](https://nilearn.github.io/index.html) package to plot the orthographic images
> 2. [`subplots`](https://matplotlib.org/api/_as_gen/matplotlib.pyplot.subplots.html) from `matplotlib.pyplot` creates a figure and multiple subplots
> 3. [`find_xyz_cut_coords`](https://nilearn.github.io/modules/generated/nilearn.plotting.find_xyz_cut_coords.html) from the [`nilearn`](https://nilearn.github.io/index.html) package will find the image coordinates of the center of the largest activation connected component
> 4. [`zip`](https://docs.python.org/3.3/library/functions.html#zip) takes iterables and aggregates them in a tuple. Here it is used to iterate through two lists simultaneously
> 5. [`iter_img`](https://nilearn.github.io/modules/generated/nilearn.image.iter_img.html) from the [`nilearn`](https://nilearn.github.io/index.html) package creates an iterator from an image that steps through each volume/time-point of the image
%% Cell type:code id: tags:
``` python
def map_plot(d):
N = d.shape[-1]
fig, ax = plt.subplots(int(np.ceil((N/2))),2, figsize=(12, N))
for img, ax0 in zip(image.iter_img(d), ax.ravel()):
coord = plotting.find_xyz_cut_coords(img, activation_threshold=3.5)
plotting.plot_stat_map(img, cut_coords=coord, vmax=10, axes=ax0)
return fig
```
%% Cell type:markdown id: tags:
Hopefully you can see some familiar looking RSN spatial patterns:
%% Cell type:code id: tags:
``` python
ics = nb.load('matmigp.gica/melodic_IC.nii.gz')
fig = map_plot(ics)
```
talks/matlab_vs_python/migp/python_MIGP.ipynb 0 → 100644
View file @ 15d75bb5
%% Cell type:markdown id: tags:
## Python MIGP
This notebook will perform *python* MIGP dimension reduction, run group ICA, and then plot the group ICs.
* [Run python `MIGP`](#run-python-migp)
* [Run `melodic`](#run-python-melodic)
* [Plot group ICs](#plot-python-group-ics)
Firstly we will import the necessary packages for this notebook:
%% Cell type:code id: tags:
``` python
import glob
import random
import nibabel as nb
import numpy as np
from scipy.sparse.linalg import svds, eigs
import matplotlib.pyplot as plt
from nilearn import plotting
from nilearn import image
import os.path as op
```
%% Cell type:markdown id: tags:
It will be necessary to know the location where the data was stored so that we can load the brainmask.
> **Note**: [`expanduser`](https://docs.python.org/3.7/library/os.path.html#os.path.expanduser) will expand the `~` to the be users home directory
%% Cell type:code id: tags:
``` python
data_dir = op.expanduser('~/nilearn_data')
```
%% Cell type:markdown id: tags:
<a class="anchor" id="run-python-migp"></a>
### Run python `MIGP`
Firstly we need to set the MIGP parameters:
> **Note:**
> 1. [`glob.glob`](https://docs.python.org/3/library/glob.html) will create a list of filenames that match the glob/wildcard pattern
> 2. [`nb.load`](https://nipy.org/nibabel/gettingstarted.html) from the [`nibabel`](https://nipy.org/nibabel/index.html) package will load the image into [`nibabel.Nifti1Image`](https://nipy.org/nibabel/reference/nibabel.nifti1.html) object. This will not load the actual data though.
> 3. We use a list comprehension to loop through all the filenames and load them with [`nibabel`](https://nipy.org/nibabel/index.html)
%% Cell type:code id: tags:
``` python
# create lists of (nibabel) image objects
in_list = [nb.load(f) for f in glob.glob(f'{data_dir}/cobre/fmri_*_smooth.nii.gz')]
in_mask = nb.load(f'{data_dir}/brain_mask.nii.gz')
# set user parameters (equivalent to melodic defaults)
GO = 'pyMIGP.nii.gz' # output filename
dPCA_int = 299 # internal number of components - typically 2-4 times number of timepoints in each run (if you have enough RAM for that)
dPCA_out = 299 # number of eigenvectors to output - should be less than dPCAint and more than the final ICA dimensionality
sep_vn = False # switch on separate variance nomalisation for each input dataset
```
%% Cell type:markdown id: tags:
> **Note:**
> 1. [`random.shuffle`](https://docs.python.org/3.7/library/random.html#random.shuffle) will shuffle a list, in this instance it shuffles the list of [`nibabel.Nifti1Image`](https://nipy.org/nibabel/reference/nibabel.nifti1.html) objects
> 2. [`ravel`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.ravel.html) will unfold a n-d array into vector. Similar to the `:` operator in Matlab
> 3. [`reshape`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.reshape.html) works similarly to reshape in Matlab, but be careful becase the default order is different from Matlab.
> 4. [`.T`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.T.html) does a transpose in `numpy`
> 5. The final element of an array is indexed with `-1` in `numpy`, as opposed to `end` in Matlab
> 6. [`svds`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.svds.html) and [`eigs`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.eigs.html?highlight=eigs#scipy.sparse.linalg.eigs) come from the [`scipy.sparse.linalg`](https://docs.scipy.org/doc/scipy/reference/sparse.linalg.html) package
> 7. [`svds`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.svds.html) and [`eigs`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.eigs.html?highlight=eigs#scipy.sparse.linalg.eigs) are very similar to their Matlab counterparts, but be careful because Matlab `svds` returns $U$, $S$, and $V$, whereas python `svds` returns $U$, $S$, and $V^T$
> 8. We index into the output of `eigs(W@W.T, dPCA_int)[1]` to only return the 2nd output (index 1)
%% Cell type:code id: tags:
``` python
# randomise the subject order
random.shuffle(in_list)
# load and unravel brainmask
mask = in_mask.get_fdata(caching='unchanged') != 0.0
# function to demean the data
def demean(x):
return x - np.mean(x, axis=0)
# loop through input files/subjects
for i, f in enumerate(in_list):
# read data
print(f'Reading data file {f.get_filename()}')
grot = f.get_fdata(caching='unchanged')[mask].T
# demean
print(f'\tRemoving mean image')
grot = demean(grot)
# var-norm
if sep_vn:
print(f'\tNormalising by voxel-wise variance')
[uu, ss, vt] = svds(grot, k=30)
vt[np.abs(vt) < (2.3 * np.std(vt))] = 0;
stddevs = np.maximum(np.std(grot - (uu @ np.diag(ss) @ vt), axis=0), 0.001)
grot = grot/stddevs
if i == 0:
W = demean(grot)
else:
# concat
W = np.concatenate((W, demean(grot)), axis=0)
# reduce W to dPCA_int eigenvectors
if W.shape[0]-10 > dPCA_int:
print(f'\tReducing data matrix to a {dPCA_int} dimensional subspace')
uu = eigs(np.dot(W, W.T), dPCA_int)[1]
uu = np.real(uu)
W = np.dot(uu.T, W)
# reshape and save
grot = np.zeros([mask.size, dPCA_out])
grot[mask.ravel(), :] = W[:dPCA_out, :].T
grot = np.reshape(grot, in_list[0].shape[:3] + (dPCA_out,))
print(f'Save to {GO}')
nb.Nifti1Image(grot, affine=in_list[0].affine).to_filename(GO)
```
%% Cell type:markdown id: tags:
<a class="anchor" id="run-python-melodic"></a>
### Run ```melodic```
Generate a command line string and run group ```melodic``` on the Matlab MIGP dimension reduced data with a dimension of 10 components. We disable MIGP because it was already run separately in Matlab.
> **Note**:
> 1. Here we use python [f-strings](https://www.python.org/dev/peps/pep-0498/), formally known as literal string interpolation, which allow for easy formatting
> 2. [`op.join`](https://docs.python.org/3.7/library/os.path.html#os.path.join) will join path strings using the platform-specific directory separator
%% Cell type:code id: tags:
``` python
# generate melodic command line string
melodic_cmd = f"melodic -i pyMIGP.nii.gz --mask={op.join(data_dir, 'brain_mask.nii.gz')} -d 10 -v --nobet --disableMigp -o pymigp.gica"
print(melodic_cmd)
```
%% Cell type:markdown id: tags:
> **Note:**
> 1. Here we use the `!` operator to execute the command in the shell
> 2. The `{}` will expand the contained python variable in the shell
%% Cell type:code id: tags:
``` python
# run melodic
! {melodic_cmd}
```
%% Cell type:markdown id: tags:
<a class="anchor" id="plot-python-group-ics"></a>
### Plot group ICs
Now we can load and plot the group ICs generated by ```melodic```.
This function will be used to plot ICs:
> **NOTE:**
> 1. Here we use [`plot_stat_map`](https://nilearn.github.io/modules/generated/nilearn.plotting.plot_stat_map.html) from the [`nilearn`](https://nilearn.github.io/index.html) package to plot the orthographic images
> 2. [`subplots`](https://matplotlib.org/api/_as_gen/matplotlib.pyplot.subplots.html) from `matplotlib.pyplot` creates a figure and multiple subplots
> 3. [`find_xyz_cut_coords`](https://nilearn.github.io/modules/generated/nilearn.plotting.find_xyz_cut_coords.html) from the [`nilearn`](https://nilearn.github.io/index.html) package will find the image coordinates of the center of the largest activation connected component
> 4. [`zip`](https://docs.python.org/3.3/library/functions.html#zip) takes iterables and aggregates them in a tuple. Here it is used to iterate through two lists simultaneously
> 5. [`iter_img`](https://nilearn.github.io/modules/generated/nilearn.image.iter_img.html) from the [`nilearn`](https://nilearn.github.io/index.html) package creates an iterator from an image that steps through each volume/time-point of the image
%% Cell type:code id: tags:
``` python
def map_plot(d):
N = d.shape[-1]
fig, ax = plt.subplots(int(np.ceil((N/2))),2, figsize=(12, N))
for img, ax0 in zip(image.iter_img(d), ax.ravel()):
coord = plotting.find_xyz_cut_coords(img, activation_threshold=3.5)
plotting.plot_stat_map(img, cut_coords=coord, vmax=10, axes=ax0)
return fig
```
%% Cell type:markdown id: tags:
Hopefully you can see some familiar looking RSN spatial patterns:
%% Cell type:code id: tags:
``` python
ics = nb.load('pymigp.gica/melodic_IC.nii.gz')
fig = map_plot(ics)
```
%% Cell type:code id: tags:
``` python
```