Commit 15d75bb5 authored by Paul McCarthy's avatar Paul McCarthy 🚵
Browse files

Merge branch 'enh/transform' into 'master'

fslpy tweaks

See merge request !29
parents 24882ffa 30e8c661
......@@ -15,21 +15,16 @@
This practical highlights some of the most useful features provided by
`fslpy`. You may find `fslpy` useful if you are writing Python code to
perform analyses and image processing in conjunction with FSL.
> **Note**: `fslpy` is distinct from `fslpython` - `fslpython` is the Python
> environment that is baked into FSL. `fslpy` is a Python library which is
> installed into the `fslpython` environment.
* [The `Image` class, and other data types](#the-image-class-and-other-data-types)
* [Creating images](#creating-images)
* [Working with image data](#working-with-image-data)
* [Loading other file types](#loading-other-file-types)
* [NIfTI coordinate systems](#nifti-coordinate-systems)
* [Image processing](#image-processing)
* [Transformations and resampling](#transformations-and-resampling)
* [FSL wrapper functions](#fsl-wrapper-functions)
* [In-memory images](#in-memory-images)
* [Loading outputs into Python](#loading-outputs-into-python)
* [The `fslmaths` wrapper](#the-fslmaths-wrapper)
* [The `FileTree`](#the-filetree)
......@@ -45,10 +40,15 @@
* [Querying atlases](#querying-atlases)
* [Loading atlas images](#loading-atlas-images)
* [Working with atlases](#working-with-atlases)
> **Note**: `fslpy` is distinct from `fslpython` - `fslpython` is the Python
> environment that is baked into FSL. `fslpy` is a Python library which is
> installed into the `fslpython` environment.
Let's start with some standard imports and environment set-up:
%% Cell type:code id: tags:
```
......@@ -58,10 +58,11 @@
import os.path as op
import nibabel as nib
import numpy as np
import warnings
warnings.filterwarnings("ignore")
np.set_printoptions(suppress=True, precision=4)
```
%% Cell type:markdown id: tags:
And a little function that we can use to generate a simple orthographic plot:
......@@ -81,10 +82,12 @@
:returns: The figure and orthogaxes (which can be passed back in as the
`fig` argument to plot overlays).
"""
voxel = [int(round(v)) for v in voxel]
data = np.asanyarray(data, dtype=np.float)
data[data <= 0] = np.nan
x, y, z = voxel
xslice = np.flipud(data[x, :, :].T)
......@@ -165,10 +168,18 @@
- Access to affine transformations between the voxel, FSL and world coordinate
systems
- Ability to load metadata from BIDS sidecar files
> The `Image` class behaves differently to the `nibabel.Nifti1Image`. For
> example, when you create an `Image` object, the default behaviour is to load
> the image data into memory. This is configurable however; take a look at
> [the
> documentation](https://users.fmrib.ox.ac.uk/~paulmc/fsleyes/fslpy/latest/fsl.data.image.html#fsl.data.image.Image)
> to explore all of the options.
Some simple image processing routines are also provided - these are covered
[below](#image-processing).
<a class="anchor" id="creating-images"></a>
......@@ -185,12 +196,12 @@
stddir = op.expandvars('${FSLDIR}/data/standard/')
# load a FSL image - the file
# suffix is optional, just like
# in real FSL-land!
img = Image(op.join(stddir, 'MNI152_T1_1mm'))
print(img)
std1mm = Image(op.join(stddir, 'MNI152_T1_1mm'))
print(std1mm)
```
%% Cell type:markdown id: tags:
You can create an `Image` from an existing `nibabel` image:
......@@ -199,32 +210,45 @@
```
# load a nibabel image, and
# convert it into an FSL image
nibimg = nib.load(op.join(stddir, 'MNI152_T1_1mm.nii.gz'))
img = Image(nibimg)
std1mm = Image(nibimg)
```
%% Cell type:markdown id: tags:
Or you can create an `Image` from a `numpy` array:
%% Cell type:code id: tags:
```
data = np.zeros((100, 100, 100))
data = np.zeros((182, 218, 182))
img = Image(data, xform=np.eye(4))
```
%% Cell type:markdown id: tags:
If you have generated some data from another `Image` (or from a
`nibabel.Nifti1Image`) you can use the `header` option to set
the header information on the new image:
%% Cell type:code id: tags:
```
img = Image(data, header=std1mm.header)
```
%% Cell type:markdown id: tags:
You can save an image to file via the `save` method:
%% Cell type:code id: tags:
```
img.save('empty.nii.gz')
img.save('empty')
!ls
```
%% Cell type:markdown id: tags:
`Image` objects have all of the attributes you might expect:
......@@ -325,11 +349,11 @@
and neuroimaging data. Most of these are higher-level wrappers around the
corresponding `nibabel` types:
* The
[`fsl.data.bitmap.Bitmap`](https://users.fmrib.ox.ac.uk/~paulmc/fsleyes/fslpy/latest/fsl.data.bitmap.html)
class can be used to load a bitmap image (e.g. `jpg, `png`, etc) and
class can be used to load a bitmap image (e.g. `jpg`, `png`, etc) and
convert it to a NIfTI image.
* The
[`fsl.data.dicom.DicomImage`](https://users.fmrib.ox.ac.uk/~paulmc/fsleyes/fslpy/latest/fsl.data.dicom.html)
class uses `dcm2niix` to load NIfTI images contained within a DICOM
directory<sup>*</sup>.
......@@ -390,15 +414,24 @@
vox2world = std2mm.getAffine('voxel', 'world')
# Apply the world->voxel
# affine to the coordinates
voxcoords = (np.dot(world2vox[:3, :3], mnicoords.T)).T + world2vox[:3, 3]
```
%% Cell type:markdown id: tags:
The code above is a bit fiddly, so instead of figuring it out, you can just
use the
[`affine.transform`](https://users.fmrib.ox.ac.uk/~paulmc/fsleyes/fslpy/latest/fsl.transform.affine.html#fsl.transform.affine.transform)
function:
%% Cell type:code id: tags:
# The code above is a bit fiddly, so
# instead of figuring it out, you can
# just use the transform() function:
```
from fsl.transform.affine import transform
voxcoords = transform(mnicoords, world2vox)
# just to double check, let's transform
# those voxel coordinates back into world
# coordinates
......@@ -504,15 +537,16 @@
# To combine affines together, we
# have to list them in reverse -
# linear algebra is *weird*.
funcvox2mni = concat(fsl2mni, func2std, vox2fsl)
print(funcvox2mni)
```
%% Cell type:markdown id: tags:
> Below we will use the
> In the next section we will use the
> [`fsl.transform.flirt.fromFlirt`](https://users.fmrib.ox.ac.uk/~paulmc/fsleyes/fslpy/latest/fsl.transform.flirt.html#fsl.transform.flirt.fromFlirt)
> function, which does all of the above for us.
So we've now got some voxel coordinates from our functional data, and an
......@@ -528,55 +562,70 @@
print('Peak activation (MNI voxels): ', mnivoxels)
```
%% Cell type:markdown id: tags:
> Note that in the above example we are only applying a linear transformation
> into MNI space - in reality you would also want to apply your non-linear
> structural-to-standard transformation too. But this is left as [an exercise
> for the
> reader](https://users.fmrib.ox.ac.uk/~paulmc/fsleyes/fslpy/latest/fsl.transform.fnirt.html).
Note that in the above example we are only applying a linear transformation
into MNI space - in reality you would also want to apply your non-linear
structural-to-standard transformation too. This is covered in the next
section.
<a class="anchor" id="image-processing"></a>
### Image processing
<a class="anchor" id="transformations-and-resampling"></a>
### Transformations and resampling
Now, it's all well and good to look at t-statistic values and voxel
coordinates and so on and so forth, but let's spice things up a bit and look
at some images. Let's display our peak activation location in MNI space. To do
this, we're going to resample our functional image into MNI space, so we can
overlay it on the MNI template. This can be done using some handy functions
from the
[`fsl.transform.flirt`](https://users.fmrib.ox.ac.uk/~paulmc/fsleyes/fslpy/latest/fsl.transform.flirt.html)
and
[`fsl.utils.image.resample`](https://users.fmrib.ox.ac.uk/~paulmc/fsleyes/fslpy/latest/fsl.utils.image.resample.html)
module:
modules.
Let's make sure we've got our source and reference images loaded:
%% Cell type:code id: tags:
```
from fsl.transform.flirt import fromFlirt
from fsl.utils.image.resample import resampleToReference
featdir = op.join(op.join('08_fslpy', 'fmri.feat'))
tstat1 = Image(op.join(featdir, 'stats', 'tstat1'))
std = Image(op.expandvars(op.join('$FSLDIR', 'data', 'standard', 'MNI152_T1_2mm')))
```
%% Cell type:markdown id: tags:
Now we'll load the `example_func2standard` FLIRT matrix, and adjust it so that
it transforms from functional *world* coordinates into standard *world*
coordinates - this is what is expected by the `resampleToReference` function,
used below:
%% Cell type:code id: tags:
```
from fsl.transform.flirt import fromFlirt
# Load the func2standard FLIRT matrix, and adjust it
# so that it transforms from functional *world*
# coordinates into standard *world* coordinates -
# this is what is expected by the resampleToReference
# function, used below
func2std = np.loadtxt(op.join(featdir, 'reg', 'example_func2standard.mat'))
func2std = fromFlirt(func2std, tstat1, std, 'world', 'world')
```
%% Cell type:markdown id: tags:
Now we can use `resampleToReference` to resample our functional data into
MNI152 space. This function returns a `numpy` array containing the resampled
data, and an adjusted voxel-to-world affine transformation. But in this case,
we know that the data will be aligned to MNI152, so we can ignore the affine:
%% Cell type:code id: tags:
```
from fsl.utils.image.resample import resampleToReference
# All of the functions in the resample module
# return a numpy array containing the resampled
# data, and an adjusted voxel-to-world affine
# transformation. But when using the
# resampleToReference function, the affine will
# be the same as the MNI152 2mm affine, so we
# can ignore it.
std_tstat1 = resampleToReference(tstat1, std, func2std)[0]
std_tstat1 = Image(std_tstat1, header=std.header)
```
%% Cell type:markdown id: tags:
......@@ -597,15 +646,156 @@
fig = ortho(std_tstat1, mnivoxels, cmap=plt.cm.inferno, fig=fig, cursor=True)
```
%% Cell type:markdown id: tags:
In the example above, we resampled some data from functional space into
standard space using a linear transformation. But we all know that this is not
how things work in the real world - linear transformations are for kids. The
real world is full of lions and tigers and bears and warp fields.
The
[`fsl.transform.fnirt`](https://users.fmrib.ox.ac.uk/~paulmc/fsleyes/fslpy/latest/fsl.transform.fnirt.html#fsl.transform.fnirt.fromFnirt)
and
[`fsl.transform.nonlinear`](https://users.fmrib.ox.ac.uk/~paulmc/fsleyes/fslpy/latest/fsl.transform.nonlinear.html)
modules contain classes and functions for working with FNIRT-style warp fields
(modules for working with lions, tigers, and bears are still under
development).
Let's imagine that we have defined an ROI in MNI152 space, and we want to
project it into the space of our functional data. We can do this by combining
the nonlinear structural to standard registration produced by FNIRT with the
linear functional to structural registration generated by FLIRT. First of
all, we'll load images from each of the functional, structural, and standard
spaces:
%% Cell type:code id: tags:
```
featdir = op.join('08_fslpy', 'fmri.feat')
func = Image(op.join(featdir, 'reg', 'example_func'))
struc = Image(op.join(featdir, 'reg', 'highres'))
std = Image(op.expandvars(op.join('$FSLDIR', 'data', 'standard', 'MNI152_T1_2mm')))
```
%% Cell type:markdown id: tags:
Now, let's say we have obtained our seed location in MNI152 coordinates. Let's
convert them to MNI152 voxels just to double check:
%% Cell type:code id: tags:
```
seedmni = [-48, -74, -9]
seedmnivox = transform(seedmni, std.getAffine('world', 'voxel'))
ortho(std.data, seedmnivox, cursor=True)
```
%% Cell type:markdown id: tags:
Now we'll load the FNIRT warp field, which encodes a nonlinear transformation
from structural space to standard space. FNIRT warp fields are often stored as
*coefficient* fields to reduce the file size, but in order to use it, we must
convert the coefficient field into a *deformation* (a.k.a. *displacement*)
field. This takes a few seconds:
%% Cell type:code id: tags:
```
from fsl.transform.fnirt import readFnirt
from fsl.transform.nonlinear import coefficientFieldToDeformationField
struc2std = readFnirt(op.join(featdir, 'reg', 'highres2standard_warp'), struc, std)
struc2std = coefficientFieldToDeformationField(struc2std)
```
%% Cell type:markdown id: tags:
We'll also load our FLIRT functional to structural transformation, adjust it
so that it transforms between voxel coordinate systems instead of the FSL
coordinate system, and invert so it can transform from structural voxels to
functional voxels:
%% Cell type:code id: tags:
```
from fsl.transform.affine import invert
func2struc = np.loadtxt(op.join(featdir, 'reg', 'example_func2highres.mat'))
func2struc = fromFlirt(func2struc, func, struc, 'voxel', 'voxel')
struc2func = invert(func2struc)
```
%% Cell type:markdown id: tags:
Now we can transform our seed coordinates from MNI152 space into functional
space in two stages. First, we'll use our deformation field to transform from
MNI152 space into structural space:
%% Cell type:code id: tags:
```
seedstruc = struc2std.transform(seedmni, 'world', 'voxel')
seedfunc = transform(seedstruc, struc2func)
print('Seed location in MNI coordinates: ', seedmni)
print('Seed location in functional voxels:', seedfunc)
ortho(func.data, seedfunc, cursor=True)
```
%% Cell type:markdown id: tags:
> FNIRT warp fields kind of work backwards - we can use them to transform
> reference coordinates into source coordinates, but would need to invert the
> warp field using `invwarp` if we wanted to transform from source coordinates
> into referemce coordinates.
Of course, we can also use our deformation field to resample an image from
structural space into MNI152 space. The `applyDeformation` function takes an
`Image` and a `DeformationField`, and returns a `numpy` array containing the
resampled data.
%% Cell type:code id: tags:
```
from fsl.transform.nonlinear import applyDeformation
strucmni = applyDeformation(struc, struc2std)
# remove low-valued voxels,
# just for visualisation below
strucmni[strucmni < 1] = 0
fig = ortho(std.data, [45, 54, 45], cmap=plt.cm.gray)
fig = ortho(strucmni, [45, 54, 45], fig=fig)
```
%% Cell type:markdown id: tags:
The `premat` option to `applyDeformation` can be used to specify our linear
functional to structural transformation, and hence resample a functional image
into MNI152 space:
%% Cell type:code id: tags:
```
tstatmni = applyDeformation(tstat1, struc2std, premat=func2struc)
tstatmni[tstatmni < 3] = 0
fig = ortho(std.data, [45, 54, 45], cmap=plt.cm.gray)
fig = ortho(tstatmni, [45, 54, 45], fig=fig)
```
%% Cell type:markdown id: tags:
There are a few other useful functions tucked away in the
[fsl.utils.image](https://users.fmrib.ox.ac.uk/~paulmc/fsleyes/fslpy/latest/fsl.utils.image.html)
package, with more to be added in the future. The [`fsl.transform`]() package
also contains a wealth of functionality for working with linear (FLIRT) and
non-linear (FNIRT) transformations.
[`fsl.utils.image`](https://users.fmrib.ox.ac.uk/~paulmc/fsleyes/fslpy/latest/fsl.utils.image.html)
and
[`fsl.transform`](https://users.fmrib.ox.ac.uk/~paulmc/fsleyes/fslpy/latest/fsl.transform.html)
packages, with more to be added in the future.
<a class="anchor" id="fsl-wrapper-functions"></a>
## FSL wrapper functions
......@@ -625,11 +815,11 @@
corresponding tool via the command-line:
%% Cell type:code id: tags:
```
from fsl.wrappers import bet, robustfov, LOAD
from fsl.wrappers import robustfov
robustfov('08_fslpy/bighead', 'bighead_cropped')
render('08_fslpy/bighead bighead_cropped -cm blue')
```
......@@ -663,10 +853,12 @@
positional arguments, and pass the additional options as keyword arguments:
%% Cell type:code id: tags:
```
from fsl.wrappers import bet
bet('bighead_cropped', 'bighead_cropped_brain', f=0.3, m=True, s=True)
render('bighead_cropped -b 40 '
'bighead_cropped_brain -cm hot '
'bighead_cropped_brain_skull -ot mask -mc 0.4 0.4 1 '
......@@ -724,10 +916,12 @@
files produced by the tool automatically loaded into memory for you:
%% Cell type:code id: tags:
```
from fsl.wrappers import LOAD
cropped = Image('bighead_cropped')
# The loaded result is called "output",
# because that is the name of the
# argument in the bet wrapper function.
......
......@@ -15,17 +15,12 @@ This practical highlights some of the most useful features provided by
perform analyses and image processing in conjunction with FSL.
> **Note**: `fslpy` is distinct from `fslpython` - `fslpython` is the Python
> environment that is baked into FSL. `fslpy` is a Python library which is
> installed into the `fslpython` environment.
* [The `Image` class, and other data types](#the-image-class-and-other-data-types)
* [Creating images](#creating-images)
* [Working with image data](#working-with-image-data)
* [Loading other file types](#loading-other-file-types)
* [NIfTI coordinate systems](#nifti-coordinate-systems)
* [Image processing](#image-processing)
* [Transformations and resampling](#transformations-and-resampling)
* [FSL wrapper functions](#fsl-wrapper-functions)
* [In-memory images](#in-memory-images)
* [Loading outputs into Python](#loading-outputs-into-python)
......@@ -45,6 +40,11 @@ perform analyses and image processing in conjunction with FSL.
* [Working with atlases](#working-with-atlases)
> **Note**: `fslpy` is distinct from `fslpython` - `fslpython` is the Python
> environment that is baked into FSL. `fslpy` is a Python library which is
> installed into the `fslpython` environment.
Let's start with some standard imports and environment set-up:
......@@ -57,6 +57,7 @@ import nibabel as nib
import numpy as np
import warnings
warnings.filterwarnings("ignore")
np.set_printoptions(suppress=True, precision=4)
```
......@@ -78,6 +79,8 @@ def ortho(data, voxel, fig=None, cursor=False, **kwargs):
`fig` argument to plot overlays).
"""
voxel = [int(round(v)) for v in voxel]
data = np.asanyarray(data, dtype=np.float)
data[data <= 0] = np.nan
......@@ -159,6 +162,14 @@ as:
- Ability to load metadata from BIDS sidecar files
> The `Image` class behaves differently to the `nibabel.Nifti1Image`. For
> example, when you create an `Image` object, the default behaviour is to load
> the image data into memory. This is configurable however; take a look at
> [the
> documentation](https://users.fmrib.ox.ac.uk/~paulmc/fsleyes/fslpy/latest/fsl.data.image.html#fsl.data.image.Image)
> to explore all of the options.
Some simple image processing routines are also provided - these are covered
[below](#image-processing).
......@@ -178,8 +189,8 @@ stddir = op.expandvars('${FSLDIR}/data/standard/')
# load a FSL image - the file
# suffix is optional, just like
# in real FSL-land!
img = Image(op.join(stddir, 'MNI152_T1_1mm'))
print(img)
std1mm = Image(op.join(stddir, 'MNI152_T1_1mm'))
print(std1mm)
```
......@@ -190,7 +201,7 @@ You can create an `Image` from an existing `nibabel` image:
# load a nibabel image, and
# convert it into an FSL image
nibimg = nib.load(op.join(stddir, 'MNI152_T1_1mm.nii.gz'))
img = Image(nibimg)
std1mm = Image(nibimg)
```
......@@ -198,15 +209,26 @@ Or you can create an `Image` from a `numpy` array:
```
data = np.zeros((100, 100, 100))
data = np.zeros((182, 218, 182))
img = Image(data, xform=np.eye(4))
```
If you have generated some data from another `Image` (or from a
`nibabel.Nifti1Image`) you can use the `header` option to set
the header information on the new image:
```
img = Image(data, header=std1mm.header)
```
You can save an image to file via the `save` method:
```
img.save('empty.nii.gz')
img.save('empty')
!ls
```
......@@ -300,7 +322,7 @@ corresponding `nibabel` types:
* The
[`fsl.data.bitmap.Bitmap`](https://users.fmrib.ox.ac.uk/~paulmc/fsleyes/fslpy/latest/fsl.data.bitmap.html)
class can be used to load a bitmap image (e.g. `jpg, `png`, etc) and
class can be used to load a bitmap image (e.g. `jpg`, `png`, etc) and
convert it to a NIfTI image.
* The
[`fsl.data.dicom.DicomImage`](https://users.fmrib.ox.ac.uk/~paulmc/fsleyes/fslpy/latest/fsl.data.dicom.html)
......@@ -364,11 +386,18 @@ vox2world = std2mm.getAffine('voxel', 'world')
# Apply the world->voxel
# affine to the coordinates
voxcoords = (np.dot(world2vox[:3, :3], mnicoords.T)).T + world2vox[:3, 3]
```
# The code above is a bit fiddly, so
# instead of figuring it out, you can
# just use the transform() function:
The code above is a bit fiddly, so instead of figuring it out, you can just
use the
[`affine.transform`](https://users.fmrib.ox.ac.uk/~paulmc/fsleyes/fslpy/latest/fsl.transform.affine.html#fsl.transform.affine.transform)
function:
```
from fsl.transform.affine import transform
voxcoords = transform(mnicoords, world2vox)
# just to double check, let's transform
......@@ -468,9 +497,10 @@ from fsl.transform.affine import concat
# have to list them in reverse -
# linear algebra is *weird*.
funcvox2mni = concat(fsl2mni, func2std, vox2fsl)
print(funcvox2mni)
```
> Below we will use the
> In the next section we will use the
> [`fsl.transform.flirt.fromFlirt`](https://users.fmrib.ox.ac.uk/~paulmc/fsleyes/fslpy/latest/fsl.transform.flirt.html#fsl.transform.flirt.fromFlirt)
> function, which does all of the above for us.
......@@ -488,15 +518,14 @@ print('Peak activation (MNI voxels): ', mnivoxels)
```
> Note that in the above example we are only applying a linear transformation
> into MNI space - in reality you would also want to apply your non-linear
> structural-to-standard transformation too. But this is left as [an exercise
> for the
> reader](https://users.fmrib.ox.ac.uk/~paulmc/fsleyes/fslpy/latest/fsl.transform.fnirt.html).
Note that in the above example we are only applyin