diff --git a/advanced_topics/08_fslpy.ipynb b/advanced_topics/08_fslpy.ipynb
index bd6881587408aacc837eb5993c71d5c686ad80e7..4895b8427fd4d443cab3c951228c6d9b794de17b 100644
--- a/advanced_topics/08_fslpy.ipynb
+++ b/advanced_topics/08_fslpy.ipynb
@@ -26,13 +26,112 @@
" * [Creating images](#creating-images)\n",
" * [Working with image data](#working-with-image-data)\n",
" * [Loading other file types](#loading-other-file-types)\n",
- "* [FSL atlases](#fsl-atlases)\n",
+ " * [NIfTI coordinate systems](#nifti-coordinate-systems)\n",
+ " * [Image processing](#image-processing)\n",
"* [The `filetree`](#the-filetree)\n",
- "* [Image processing](#image-processing)\n",
"* [FSL wrapper functions](#fsl-wrapper-functions)\n",
- "* [NIfTI coordinate systems](#nifti-coordinate-systems)\n",
+ "* [FSL atlases](#fsl-atlases)\n",
+ " * [Querying atlases](#querying-atlases)\n",
+ " * [Loading atlas images](#loading-atlas-images)\n",
+ " * [Working with atlases](#working-with-atlases)\n",
+ "\n",
+ "\n",
+ "Let's start with some standard imports and environment set-up:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "%matplotlib inline\n",
+ "import matplotlib.pyplot as plt\n",
+ "import os\n",
+ "import os.path as op\n",
+ "import nibabel as nib\n",
+ "import numpy as np\n",
+ "import warnings\n",
+ "warnings.filterwarnings(\"ignore\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "And a little function that we can use to generate a simple orthographic plot:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def ortho(data, voxel, fig=None, **kwargs):\n",
+ " \"\"\"Simple orthographic plot of a 3D array using matplotlib.\n",
+ "\n",
+ " :arg data: 3D numpy array\n",
+ " :arg voxel: XYZ coordinates for each slice\n",
+ " :arg fig: Existing figure and axes for overlay plotting\n",
"\n",
+ " All other arguments are passed through to the `imshow` function.\n",
"\n",
+ " :returns: The figure and axes (which can be passed back in as the\n",
+ " `fig` argument to plot overlays).\n",
+ " \"\"\"\n",
+ " x, y, z = voxel\n",
+ " xslice = np.flipud(data[x, :, :].T)\n",
+ " yslice = np.flipud(data[:, y, :].T)\n",
+ " zslice = np.flipud(data[:, :, z].T)\n",
+ "\n",
+ " if fig is None:\n",
+ " fig = plt.figure()\n",
+ " xax = fig.add_subplot(1, 3, 1)\n",
+ " yax = fig.add_subplot(1, 3, 2)\n",
+ " zax = fig.add_subplot(1, 3, 3)\n",
+ " else:\n",
+ " fig, xax, yax, zax = fig\n",
+ "\n",
+ " xax.imshow(xslice, **kwargs)\n",
+ " yax.imshow(yslice, **kwargs)\n",
+ " zax.imshow(zslice, **kwargs)\n",
+ "\n",
+ " for ax in (xax, yax, zax):\n",
+ " ax.set_xticks([])\n",
+ " ax.set_yticks([])\n",
+ " fig.tight_layout(pad=0)\n",
+ "\n",
+ " return (fig, xax, yax, zax)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "And another function which uses FSLeyes for more complex plots:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import shlex\n",
+ "import IPython.display as display\n",
+ "from fsleyes.render import main\n",
+ "\n",
+ "def render(cmdline):\n",
+ " prefix = '-of screenshot.png -hl -c 2 '\n",
+ " main(shlex.split(prefix + cmdline))\n",
+ " return display.Image('screenshot.png')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
"<a class=\"anchor\" id=\"the-image-class-and-other-data-types\"></a>\n",
"## The `Image` class, and other data types\n",
"\n",
@@ -67,7 +166,6 @@
"metadata": {},
"outputs": [],
"source": [
- "import os.path as op\n",
"from fsl.data.image import Image\n",
"\n",
"stddir = op.expandvars('${FSLDIR}/data/standard/')\n",
@@ -92,8 +190,6 @@
"metadata": {},
"outputs": [],
"source": [
- "import nibabel as nib\n",
- "\n",
"# load a nibabel image, and\n",
"# convert it into an FSL image\n",
"nibimg = nib.load(op.join(stddir, 'MNI152_T1_1mm.nii.gz'))\n",
@@ -113,8 +209,6 @@
"metadata": {},
"outputs": [],
"source": [
- "import numpy as np\n",
- "\n",
"data = np.zeros((100, 100, 100))\n",
"img = Image(data, xform=np.eye(4))"
]
@@ -175,7 +269,7 @@
"outputs": [],
"source": [
"std2mm = Image(op.join(stddir, 'MNI152_T1_2mm'))\n",
- "mask2mm = Image(op.join(stddir, 'MNI152_T1_2mm_mask'))\n",
+ "mask2mm = Image(op.join(stddir, 'MNI152_T1_2mm_brain_mask'))\n",
"\n",
"print(std1mm.sameSpace(std2mm))\n",
"print(std2mm.sameSpace(mask2mm))\n",
@@ -219,7 +313,7 @@
"outputs": [],
"source": [
"data = std2mm.data\n",
- "print(data.min, data.max())"
+ "print(data.min(), data.max())"
]
},
{
@@ -240,7 +334,7 @@
"outputs": [],
"source": [
"slc = std2mm[:, :, 45]\n",
- "std2mm[0:10, :, :] = 0"
+ "std2mm[0:10, :, :] *= 2"
]
},
{
@@ -308,18 +402,73 @@
"> in order to use this class.\n",
"\n",
"\n",
- "<a class=\"anchor\" id=\"fsl-atlases\"></a>\n",
- "## FSL atlases\n",
+ "<a class=\"anchor\" id=\"nifti-coordinate-systems\"></a>\n",
+ "### NIfTI coordinate systems\n",
"\n",
"\n",
- "The\n",
- "[`fsl.data.atlases`](https://users.fmrib.ox.ac.uk/~paulmc/fsleyes/fslpy/latest/fsl.data.atlases.html)\n",
- "module provides access to all of the atlas images that are stored in the\n",
- "`$FSLDIR/data/atlases/` directory of a standard FSL installation. It can be\n",
- "used to load and query probabilistic and label-based atlases.\n",
+ "The `Image.getAffine` method gives you access to affine transformations which\n",
+ "can be used to convert coordinates between the different coordinate systems\n",
+ "associated with a NIfTI image. Have some MNI coordinates you'd like to convert\n",
+ "to voxels? Easy!"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "stddir = op.expandvars('${FSLDIR}/data/standard/')\n",
+ "std2mm = Image(op.join(stddir, 'MNI152_T1_2mm'))\n",
+ "\n",
+ "mnicoords = np.array([[0, 0, 0],\n",
+ " [0, -18, 18]])\n",
"\n",
+ "world2vox = std2mm.getAffine('world', 'voxel')\n",
+ "vox2world = std2mm.getAffine('voxel', 'world')\n",
"\n",
- "The `atlases` module needs to be initialised using the `rescanAtlases` function:"
+ "# Apply the world->voxel\n",
+ "# affine to the coordinates\n",
+ "voxcoords = (np.dot(world2vox[:3, :3], mnicoords.T)).T + world2vox[:3, 3]\n",
+ "\n",
+ "# The code above is a bit fiddly, so\n",
+ "# instead of figuring it out, you can\n",
+ "# just use the transform() function:\n",
+ "from fsl.transform.affine import transform\n",
+ "voxcoords = transform(mnicoords, world2vox)\n",
+ "\n",
+ "# just to double check, let's transform\n",
+ "# those voxel coordinates back into world\n",
+ "# coordinates\n",
+ "backtomni = transform(voxcoords, vox2world)\n",
+ "\n",
+ "for m, v, b in zip(mnicoords, voxcoords, backtomni):\n",
+ " print(m, '->', v, '->', b)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "> The `Image.getAffine` method can give you transformation matrices\n",
+ "> between any of these coordinate systems:\n",
+ ">\n",
+ "> - `'voxel'`: Image data voxel coordinates\n",
+ "> - `'world'`: mm coordinates, defined by the sform/qform of an image\n",
+ "> - `'fsl'`: The FSL coordinate system, used internally by many FSL tools\n",
+ "> (e.g. FLIRT)\n",
+ "\n",
+ "\n",
+ "Oh, that example was too easy I hear you say? Try this one on for size. Let's\n",
+ "say we have run FEAT on some task fMRI data, and want to get the MNI\n",
+ "coordinates of the voxel with peak activation.\n",
+ "\n",
+ "\n",
+ "> This is what people used to use `Featquery` for, back in the un-enlightened\n",
+ "> days.\n",
+ "\n",
+ "\n",
+ "Let's start by identifying the voxel with the biggest t-statistic:"
]
},
{
@@ -328,15 +477,27 @@
"metadata": {},
"outputs": [],
"source": [
- "import fsl.data.atlases as atlases\n",
- "atlases.rescanAtlases()"
+ "featdir = op.join(op.join('08_fslpy', 'fmri.feat'))\n",
+ "\n",
+ "tstat1 = Image(op.join(featdir, 'stats', 'tstat1')).data\n",
+ "\n",
+ "# Recall from the numpy practical that\n",
+ "# argmax gives us a 1D index into a\n",
+ "# flattened view of the array. We can\n",
+ "# use the unravel_index function to\n",
+ "# convert it into a 3D index.\n",
+ "peakvox = np.abs(tstat1).argmax()\n",
+ "peakvox = np.unravel_index(peakvox, tstat1.shape)\n",
+ "print('Peak voxel coordinates for tstat1:', peakvox, tstat1[peakvox])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "You can list all of the available atlases using `listAtlases`:"
+ "Now that we've got the voxel coordinates in functional space, we need to\n",
+ "transform them into MNI space. FEAT provides a transformation which goes\n",
+ "directly from functional to standard space, in the `reg` directory:"
]
},
{
@@ -345,18 +506,19 @@
"metadata": {},
"outputs": [],
"source": [
- "for desc in atlases.listAtlases():\n",
- " print(desc)"
+ "func2std = np.loadtxt(op.join(featdir, 'reg', 'example_func2standard.mat'))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "`listAtlases` returns a list of `AtlasDescription` objects, each of which\n",
- "contains descriptive information about one atlas. You can retrieve the\n",
- "`AtlasDescription` for a specific atlas via the `getAtlasDescription`\n",
- "function:"
+ "But ... wait a minute ... this is a FLIRT matrix! We can't just plug voxel\n",
+ "coordinates into a FLIRT matrix and expect to get sensible results, because\n",
+ "FLIRT works in an internal FSL coordinate system, which is not quite\n",
+ "`'voxel'`, and not quite `'world'`. So we need to do a little more work.\n",
+ "Let's start by loading our functional image, and the MNI152 template (the\n",
+ "source and reference images of our FLIRT matrix):"
]
},
{
@@ -365,20 +527,21 @@
"metadata": {},
"outputs": [],
"source": [
- "desc = atlases.getAtlasDescription('harvardoxford-cortical')\n",
- "print(desc.name)\n",
- "print(desc.atlasID)\n",
- "print(desc.specPath)\n",
- "print(desc.atlasType)"
+ "func = Image(op.join(featdir, 'reg', 'example_func'))\n",
+ "std = Image(op.expandvars(op.join('$FSLDIR', 'data', 'standard', 'MNI152_T1_2mm')))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "Each `AtlasDescription` maintains a list of `AtlasLabel` objects, each of\n",
- "which represents one region that is defined in the atlas. You can access all\n",
- "of the `AtlasLabel` objects via the `labels` attribute:"
+ "Now we can use them to get affines which convert between all of the different\n",
+ "coordinate systems - we're going to combine them into a single uber-affine,\n",
+ "which transforms our functional-space voxels into MNI world coordinates via:\n",
+ "\n",
+ " 1. functional voxels -> FLIRT source space\n",
+ " 2. FLIRT source space -> FLIRT reference space\n",
+ " 3. FLIRT referece space -> MNI world coordinates"
]
},
{
@@ -387,15 +550,16 @@
"metadata": {},
"outputs": [],
"source": [
- "for lbl in desc.labels[:5]:\n",
- " print(lbl)"
+ "vox2fsl = func.getAffine('voxel', 'fsl')\n",
+ "fsl2mni = std .getAffine('fsl', 'world')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "Or you can retrieve a specific label using the `find` method:"
+ "Combining two affines into one is just a simple dot-product. There is a\n",
+ "`concat()` function which does this for us, for any number of affines:"
]
},
{
@@ -404,18 +568,20 @@
"metadata": {},
"outputs": [],
"source": [
- "# search by region name\n",
- "print(desc.find(name='Occipital Pole'))\n",
+ "from fsl.transform.affine import concat\n",
"\n",
- "# or by label value\n",
- "print(desc.find(value=48))"
+ "# To combine affines together, we\n",
+ "# have to list them in reverse -\n",
+ "# linear algebra is *weird*.\n",
+ "funcvox2mni = concat(fsl2mni, func2std, vox2fsl)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "The `loadAtlas` function can be used to load the atlas image:"
+ "So we've now got some voxel coordinates from our functional data, and an\n",
+ "affine to transform into MNI world coordinates. The rest is easy:"
]
},
{
@@ -424,34 +590,36 @@
"metadata": {},
"outputs": [],
"source": [
- "# For probabilistic atlases, you\n",
- "# can ask for the 3D ROI image\n",
- "# by setting loadSummary=True.\n",
- "# You can also request a\n",
- "# resolution - by default the\n",
- "# highest resolution version\n",
- "# will be loaded.\n",
- "lblatlas = atlases.loadAtlas('harvardoxford-cortical',\n",
- " loadSummary=True,\n",
- " resolution=2)\n",
- "\n",
- "# By default you will get the 4D\n",
- "# probabilistic atlas image (for\n",
- "# atlases for which this is\n",
- "# available).\n",
- "probatlas = atlases.loadAtlas('harvardoxford-cortical',\n",
- " resolution=2)\n",
- "\n",
- "print(lblatlas)\n",
- "print(probatlas)"
+ "mnicoords = transform(peakvox, funcvox2mni)\n",
+ "mnivoxels = transform(mnicoords, std.getAffine('world', 'voxel'))\n",
+ "mnivoxels = [int(round(v)) for v in mnivoxels]\n",
+ "print('Peak activation (MNI coordinates):', mnicoords)\n",
+ "print('Peak activation (MNI voxels): ', mnivoxels)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "`LabelAtlas` objects have a method called `label`, which can be used to\n",
- "interrogate the atlas at specific locations:"
+ "> Note that in the above example we are only applying a linear transformation\n",
+ "> into MNI space - in reality you would also want to apply your non-linear\n",
+ "> structural-to-standard transformation too. But this is left as [an exercise\n",
+ "> for the\n",
+ "> reader](https://users.fmrib.ox.ac.uk/~paulmc/fsleyes/fslpy/latest/fsl.transform.fnirt.html).\n",
+ "\n",
+ "\n",
+ "<a class=\"anchor\" id=\"image-processing\"></a>\n",
+ "### Image processing\n",
+ "\n",
+ "\n",
+ "Now, it's all well and good to look at t-statistic values and voxel\n",
+ "coordinates and so on and so forth, but let's spice things up a bit and look\n",
+ "at some images. Let's display our peak activation location in MNI space. To do\n",
+ "this, we're going to resample our functional image into MNI space, so we can\n",
+ "overlay it on the MNI template. This can be done using some handy functions\n",
+ "from the\n",
+ "[`fsl.utils.image.resample`](https://users.fmrib.ox.ac.uk/~paulmc/fsleyes/fslpy/latest/fsl.utils.image.resample.html)\n",
+ "module:"
]
},
{
@@ -460,30 +628,38 @@
"metadata": {},
"outputs": [],
"source": [
- "# The label method accepts 3D\n",
- "# voxel or world coordinates\n",
- "val = lblatlas.label((25, 52, 43), voxel=True)\n",
- "lbl = lblatlas.find(value=val)\n",
- "print('Region at voxel [25, 52, 43]: {} [{}]'.format(val, lbl.name))\n",
+ "from fsl.transform.flirt import fromFlirt\n",
+ "from fsl.utils.image.resample import resampleToReference\n",
"\n",
+ "featdir = op.join(op.join('08_fslpy', 'fmri.feat'))\n",
+ "tstat1 = Image(op.join(featdir, 'stats', 'tstat1'))\n",
+ "std = Image(op.expandvars(op.join('$FSLDIR', 'data', 'standard', 'MNI152_T1_2mm')))\n",
"\n",
- "# or a 3D weighted or binary mask\n",
- "mask = np.zeros(lblatlas.shape)\n",
- "mask[30:60, 30:60, 30:60] = 1\n",
- "mask = Image(mask, header=lblatlas.header)\n",
+ "# Load the func2standard FLIRT matrix, and adjust it\n",
+ "# so that it transforms from functional *world*\n",
+ "# coordinates into standard *world* coordinates -\n",
+ "# this is what is expected by the resampleToReference\n",
+ "# function, used below\n",
+ "func2std = np.loadtxt(op.join(featdir, 'reg', 'example_func2standard.mat'))\n",
+ "func2std = fromFlirt(func2std, tstat1, std, 'world', 'world')\n",
"\n",
- "lbls, props = lblatlas.label(mask)\n",
- "print('Labels in mask:')\n",
- "for lbl, prop in zip(lbls, props):\n",
- " lblname = lblatlas.find(value=lbl).name\n",
- " print(' {} [{}]: {:0.2f}%'.format(lbl, lblname, prop))"
+ "# All of the functions in the resample module\n",
+ "# return a numpy array containing the resampled\n",
+ "# data, and an adjusted voxel-to-world affine\n",
+ "# transformation. But when using the\n",
+ "# resampleToReference function, the affine will\n",
+ "# be the same as the MNI152 2mm affine, so we\n",
+ "# can ignore it.\n",
+ "std_tstat1 = resampleToReference(tstat1, std, func2std)[0]\n",
+ "std_tstat1 = Image(std_tstat1, header=std.header)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "`ProbabilisticAtlas` objects have an analogous method called `values`:"
+ "Now that we have our t-statistic image in MNI152 space, we can plot it in\n",
+ "standard space using `matplotlib`:"
]
},
{
@@ -492,48 +668,43 @@
"metadata": {},
"outputs": [],
"source": [
- "vals = probatlas.values((25, 52, 43), voxel=True)\n",
- "print('Regions at voxel [25, 52, 43]:')\n",
- "for idx, val in enumerate(vals):\n",
- " if val > 0:\n",
- " lbl = probatlas.find(index=idx)\n",
- " print(' {} [{}]: {:0.2f}%'.format(lbl.value, lbl.name, val))\n",
+ "stddir = op.expandvars('${FSLDIR}/data/standard/')\n",
+ "std2mm = Image(op.join(stddir, 'MNI152_T1_2mm'))\n",
"\n",
- "print('Average proportions of regions within mask:')\n",
- "vals = probatlas.values(mask)\n",
- "for idx, val in enumerate(vals):\n",
- " if val > 0:\n",
- " lbl = probatlas.find(index=idx)\n",
- " print(' {} [{}]: {:0.2f}%'.format(lbl.value, lbl.name, val))"
+ "std_tstat1 = std_tstat1.data\n",
+ "std_tstat1 = np.ma.masked_where(std_tstat1 < 3, std_tstat1)\n",
+ "\n",
+ "fig = ortho(std2mm, mnivoxels, cmap=plt.cm.gray)\n",
+ "fig = ortho(std_tstat1, mnivoxels, cmap=plt.cm.inferno, fig=fig)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
+ "There are a few other useful functions tucked away in the\n",
+ "[fsl.utils.image](https://users.fmrib.ox.ac.uk/~paulmc/fsleyes/fslpy/latest/fsl.utils.image.html)\n",
+ "package, with more to be added in the future. The [`fsl.transform`]() package\n",
+ "also contains a wealth of functionality for working with linear (FLIRT) and\n",
+ "non-linear (FNIRT) transformations.\n",
+ "\n",
+ "\n",
"<a class=\"anchor\" id=\"the-filetree\"></a>\n",
"## The `filetree`\n",
"\n",
- "<a class=\"anchor\" id=\"nifti-coordinate-systems\"></a>\n",
- "## NIfTI coordinate systems\n",
- "\n",
- "<a class=\"anchor\" id=\"image-processing\"></a>\n",
- "## Image processing\n",
"\n",
"<a class=\"anchor\" id=\"fsl-wrapper-functions\"></a>\n",
"## FSL wrapper functions\n",
"\n",
"\n",
- "\n",
- "<a class=\"anchor\" id=\"nifti-coordinate-systems\"></a>\n",
- "## NIfTI coordinate systems\n",
- "\n",
+ "The\n",
+ "[fsl.wrappers](https://users.fmrib.ox.ac.uk/~paulmc/fsleyes/fslpy/latest/fsl.wrappers.html)\n",
+ "package is the home of \"wrapper\" functions for a range of FSL tools. You can\n",
+ "use them to call an FSL tool from Python code, without having to worry about\n",
+ "constructing a command-line, or saving/loading input/output images.\n",
"\n",
"\n",
- "The `getAffine` method gives you access to affine transformations which can be\n",
- "used to convert coordinates between the different coordinate systems\n",
- "associated with an image. Have some MNI coordinates you'd like to convert to\n",
- "voxels? Easy!"
+ "You can use the FSL wrapper functions with file names:"
]
},
{
@@ -542,54 +713,49 @@
"metadata": {},
"outputs": [],
"source": [
- "mnicoords = np.array([[0, 0, 0],\n",
- " [0, -18, 18]])\n",
- "\n",
- "world2vox = img.getAffine('world', 'voxel')\n",
- "vox2world = img.getAffine('voxel', 'world')\n",
- "\n",
- "# Apply the world->voxel\n",
- "# affine to the coordinates\n",
- "voxcoords = (np.dot(world2vox[:3, :3], mnicoords.T)).T + world2vox[:3, 3]\n",
- "\n",
- "# The code above is a bit fiddly, so\n",
- "# instead of figuring it out, you can\n",
- "# just use the transform() function:\n",
- "from fsl.transform.affine import transform\n",
- "voxcoords = transform(mnicoords, world2vox)\n",
- "\n",
- "# just to double check, let's transform\n",
- "# those voxel coordinates back into world\n",
- "# coordinates\n",
- "backtomni = transform(voxcoords, vox2world)\n",
- "\n",
- "for m, v, b in zip(mnicoords, voxcoords, backtomni):\n",
- " print(m, '->', v, '->', b)"
+ "from fsl.wrappers import bet, robustfov, LOAD\n",
+ "os.chdir('08_fslpy')\n",
+ "robustfov('bighead', 'bighead_cropped')\n",
+ "render('bighead bighead_cropped -cm blue')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "> The `Image.getAffine` method can give you transformation matrices\n",
- "> between any of these coordinate systems:\n",
- ">\n",
- "> - `'voxel'`: Image data voxel coordinates\n",
- "> - `'world'`: mm coordinates, defined by the sform/qform of an image\n",
- "> - `'fsl'`: The FSL coordinate system, used internally by many FSL tools\n",
- "> (e.g. FLIRT)\n",
- "\n",
+ "Or, if you have images in memory, you can pass the `Image` objects in,\n",
+ "and have the results automatically loaded in too:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "cropped = Image('bighead_cropped')\n",
+ "betted = bet(cropped, LOAD)['output']\n",
"\n",
+ "fig = ortho(cropped, (80, 112, 85))\n",
+ "fig = ortho(betted, (80, 112, 85), fig=fig, cmap=plt.cm.inferno)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "<a class=\"anchor\" id=\"fsl-atlases\"></a>\n",
+ "## FSL atlases\n",
"\n",
- "Oh, that example was too easy I hear you say? Try this one on for size. Let's\n",
- "say we have run FEAT on some task fMRI data, and want to get the MNI\n",
- "coordinates of the voxel with peak activation.\n",
"\n",
- "> This is what people used to use `Featquery` for, back in the un-enlightened\n",
- "> days.\n",
+ "The\n",
+ "[`fsl.data.atlases`](https://users.fmrib.ox.ac.uk/~paulmc/fsleyes/fslpy/latest/fsl.data.atlases.html)\n",
+ "module provides access to all of the atlas images that are stored in the\n",
+ "`$FSLDIR/data/atlases/` directory of a standard FSL installation. It can be\n",
+ "used to load and query probabilistic and label-based atlases.\n",
"\n",
"\n",
- "Let's start by identifying the voxel with the biggest t-statistic:"
+ "The `atlases` module needs to be initialised using the `rescanAtlases` function:"
]
},
{
@@ -598,30 +764,39 @@
"metadata": {},
"outputs": [],
"source": [
- "featdir = op.join(op.join('05_nifti', 'fmri.feat'))\n",
+ "import fsl.data.atlases as atlases\n",
+ "atlases.rescanAtlases()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "<a class=\"anchor\" id=\"querying-atlases\"></a>\n",
+ "### Querying atlases\n",
"\n",
- "# The Image.data attribute returns a\n",
- "# numpy array containing, well, the\n",
- "# image data.\n",
- "tstat1 = Image(op.join(featdir, 'stats', 'tstat1')).data\n",
"\n",
- "# Recall from the numpy practical that\n",
- "# argmax gives us a 1D index into a\n",
- "# flattened view of the array. We can\n",
- "# use the unravel_index function to\n",
- "# convert it into a 3D index.\n",
- "peakvox = np.abs(tstat1).argmax()\n",
- "peakvox = np.unravel_index(peakvox, tstat1.shape)\n",
- "print('Peak voxel coordinates for tstat1:', peakvox, tstat1[peakvox])"
+ "You can list all of the available atlases using `listAtlases`:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for desc in atlases.listAtlases():\n",
+ " print(desc)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "Now that we've got the voxel coordinates in functional space, we need to\n",
- "transform them into MNI space. FEAT provides a transformation which goes\n",
- "directly from functional to standard space, in the `reg` directory:"
+ "`listAtlases` returns a list of `AtlasDescription` objects, each of which\n",
+ "contains descriptive information about one atlas. You can retrieve the\n",
+ "`AtlasDescription` for a specific atlas via the `getAtlasDescription`\n",
+ "function:"
]
},
{
@@ -630,19 +805,20 @@
"metadata": {},
"outputs": [],
"source": [
- "func2std = np.loadtxt(op.join(featdir, 'reg', 'example_func2standard.mat'))"
+ "desc = atlases.getAtlasDescription('harvardoxford-cortical')\n",
+ "print(desc.name)\n",
+ "print(desc.atlasID)\n",
+ "print(desc.specPath)\n",
+ "print(desc.atlasType)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "But ... wait a minute ... this is a FLIRT matrix! We can't just plug voxel\n",
- "coordinates into a FLIRT matrix and expect to get sensible results, because\n",
- "FLIRT works in an internal FSL coordinate system, which is not quite\n",
- "`'voxel'`, and not quite `'world'`. So we need to do a little more work.\n",
- "Let's start by loading our functional image, and the MNI152 template (the\n",
- "source and reference images of our FLIRT matrix):"
+ "Each `AtlasDescription` maintains a list of `AtlasLabel` objects, each of\n",
+ "which represents one region that is defined in the atlas. You can access all\n",
+ "of the `AtlasLabel` objects via the `labels` attribute:"
]
},
{
@@ -651,21 +827,15 @@
"metadata": {},
"outputs": [],
"source": [
- "func = Image(op.join(featdir, 'reg', 'example_func'))\n",
- "std = Image(op.expandvars(op.join('$FSLDIR', 'data', 'standard', 'MNI152_T1_2mm')))"
+ "for lbl in desc.labels[:5]:\n",
+ " print(lbl)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "Now we can use them to get affines which convert between all of the different\n",
- "coordinate systems - we're going to combine them into a single uber-affine,\n",
- "which transforms our functional-space voxels into MNI world coordinates via:\n",
- "\n",
- " 1. functional voxels -> FLIRT source space\n",
- " 2. FLIRT source space -> FLIRT reference space\n",
- " 3. FLIRT referece space -> MNI world coordinates"
+ "Or you can retrieve a specific label using the `find` method:"
]
},
{
@@ -674,16 +844,22 @@
"metadata": {},
"outputs": [],
"source": [
- "vox2fsl = func.getAffine('voxel', 'fsl')\n",
- "fsl2mni = std .getAffine('fsl', 'world')"
+ "# search by region name\n",
+ "print(desc.find(name='Occipital Pole'))\n",
+ "\n",
+ "# or by label value\n",
+ "print(desc.find(value=48))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "Combining two affines into one is just a simple dot-product. There is a\n",
- "`concat()` function which does this for us, for any number of affines:"
+ "<a class=\"anchor\" id=\"loading-atlas-images\"></a>\n",
+ "### Loading atlas images\n",
+ "\n",
+ "\n",
+ "The `loadAtlas` function can be used to load the atlas image:"
]
},
{
@@ -692,20 +868,38 @@
"metadata": {},
"outputs": [],
"source": [
- "from fsl.transform.affine import concat\n",
+ "# For probabilistic atlases, you\n",
+ "# can ask for the 3D ROI image\n",
+ "# by setting loadSummary=True.\n",
+ "# You can also request a\n",
+ "# resolution - by default the\n",
+ "# highest resolution version\n",
+ "# will be loaded.\n",
+ "lblatlas = atlases.loadAtlas('harvardoxford-cortical',\n",
+ " loadSummary=True,\n",
+ " resolution=2)\n",
"\n",
- "# To combine affines together, we\n",
- "# have to list them in reverse -\n",
- "# linear algebra is *weird*.\n",
- "funcvox2mni = concat(fsl2mni, func2std, vox2fsl)"
+ "# By default you will get the 4D\n",
+ "# probabilistic atlas image (for\n",
+ "# atlases for which this is\n",
+ "# available).\n",
+ "probatlas = atlases.loadAtlas('harvardoxford-cortical',\n",
+ " resolution=2)\n",
+ "\n",
+ "print(lblatlas)\n",
+ "print(probatlas)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "So we've now got some voxel coordinates from our functional data, and an affine\n",
- "to transform into MNI world coordinates. The rest is easy:"
+ "<a class=\"anchor\" id=\"working-with-atlases\"></a>\n",
+ "### Working with atlases\n",
+ "\n",
+ "\n",
+ "Both `LabelAtlas` and `ProbabilisticAtlas` objects have a method called `get`,\n",
+ "which can be used to extract ROI images for a specific region:"
]
},
{
@@ -714,27 +908,52 @@
"metadata": {},
"outputs": [],
"source": [
- "mnicoords = transform(peakvox, funcvox2mni)\n",
- "print('Peak activation (MNI coordinates):', mnicoords)"
+ "stddir = op.expandvars('${FSLDIR}/data/standard/')\n",
+ "std2mm = Image(op.join(stddir, 'MNI152_T1_2mm'))\n",
+ "\n",
+ "frontal = lblatlas.get(name='Frontal Pole').data\n",
+ "frontal = np.ma.masked_where(frontal < 1, frontal)\n",
+ "\n",
+ "fig = ortho(std2mm, (45, 54, 45), cmap=plt.cm.gray)\n",
+ "fig = ortho(frontal, (45, 54, 45), cmap=plt.cm.winter, fig=fig)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "> Note that in the above example we are only applying a linear transformation\n",
- "> into MNI space - in reality you would also want to apply your non-linear\n",
- "> structural-to-standard transformation too. But this is left as an exercise\n",
- "> for the reader ;).\n",
+ "Calling `get` on a :meth:`ProbabilisticAtlas` will return a probability image:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "stddir = op.expandvars('${FSLDIR}/data/standard/')\n",
+ "std2mm = Image(op.join(stddir, 'MNI152_T1_2mm'))\n",
"\n",
+ "frontal = probatlas.get(name='Frontal Pole')\n",
+ "frontal = np.ma.masked_where(frontal < 1, frontal)\n",
"\n",
- "<a class=\"anchor\" id=\"image-processing\"></a>\n",
- "## Image processing\n",
+ "fig = ortho(std2mm, (45, 54, 45), cmap=plt.cm.gray)\n",
+ "fig = ortho(frontal, (45, 54, 45), cmap=plt.cm.inferno, fig=fig)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The `get` method can be used to retrieve an image for a region by:\n",
+ "- an `AtlasLabel` object\n",
+ "- The region index\n",
+ "- The region value\n",
+ "- The region name\n",
"\n",
- "Now, it's all well and good to look at t-statistric values and voxel\n",
- "coordinates and so on and so forth, but let's spice things up a bit and look\n",
- "at some images. Let's display our peak activation location in MNI space. To\n",
- "do this, we're going to resample our functional image into"
+ "\n",
+ "`LabelAtlas` objects have a method called `label`, which can be used to\n",
+ "interrogate the atlas at specific locations:"
]
},
{
@@ -743,18 +962,51 @@
"metadata": {},
"outputs": [],
"source": [
- "from IPython.display import Image as Screenshot\n",
- "!fsleyes render -of screenshot.png -std"
+ "# The label method accepts 3D\n",
+ "# voxel or world coordinates\n",
+ "val = lblatlas.label((25, 52, 43), voxel=True)\n",
+ "lbl = lblatlas.find(value=val)\n",
+ "print('Region at voxel [25, 52, 43]: {} [{}]'.format(val, lbl.name))\n",
+ "\n",
+ "\n",
+ "# or a 3D weighted or binary mask\n",
+ "mask = np.zeros(lblatlas.shape)\n",
+ "mask[30:60, 30:60, 30:60] = 1\n",
+ "mask = Image(mask, header=lblatlas.header)\n",
+ "\n",
+ "lbls, props = lblatlas.label(mask)\n",
+ "print('Labels in mask:')\n",
+ "for lbl, prop in zip(lbls, props):\n",
+ " lblname = lblatlas.find(value=lbl).name\n",
+ " print(' {} [{}]: {:0.2f}%'.format(lbl, lblname, prop))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "### (Advanced) Transform coordinates with nonlinear warpfields\n",
- "\n",
+ "`ProbabilisticAtlas` objects have an analogous method called `values`:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "vals = probatlas.values((25, 52, 43), voxel=True)\n",
+ "print('Regions at voxel [25, 52, 43]:')\n",
+ "for idx, val in enumerate(vals):\n",
+ " if val > 0:\n",
+ " lbl = probatlas.find(index=idx)\n",
+ " print(' {} [{}]: {:0.2f}%'.format(lbl.value, lbl.name, val))\n",
"\n",
- "have to use your own dataset"
+ "print('Average proportions of regions within mask:')\n",
+ "vals = probatlas.values(mask)\n",
+ "for idx, val in enumerate(vals):\n",
+ " if val > 0:\n",
+ " lbl = probatlas.find(index=idx)\n",
+ " print(' {} [{}]: {:0.2f}%'.format(lbl.value, lbl.name, val))"
]
}
],
diff --git a/advanced_topics/08_fslpy.md b/advanced_topics/08_fslpy.md
index ef8b860f4fa964fb8ff9e5f2bbd3391edaef1181..f79cc08addc91bdaf5421e5c4c43eff5db032ec3 100644
--- a/advanced_topics/08_fslpy.md
+++ b/advanced_topics/08_fslpy.md
@@ -20,11 +20,86 @@ perform analyses and image processing in conjunction with FSL.
* [Creating images](#creating-images)
* [Working with image data](#working-with-image-data)
* [Loading other file types](#loading-other-file-types)
-* [FSL atlases](#fsl-atlases)
+ * [NIfTI coordinate systems](#nifti-coordinate-systems)
+ * [Image processing](#image-processing)
* [The `filetree`](#the-filetree)
-* [Image processing](#image-processing)
* [FSL wrapper functions](#fsl-wrapper-functions)
-* [NIfTI coordinate systems](#nifti-coordinate-systems)
+* [FSL atlases](#fsl-atlases)
+ * [Querying atlases](#querying-atlases)
+ * [Loading atlas images](#loading-atlas-images)
+ * [Working with atlases](#working-with-atlases)
+
+
+Let's start with some standard imports and environment set-up:
+
+
+```
+%matplotlib inline
+import matplotlib.pyplot as plt
+import os
+import os.path as op
+import nibabel as nib
+import numpy as np
+import warnings
+warnings.filterwarnings("ignore")
+```
+
+
+And a little function that we can use to generate a simple orthographic plot:
+
+
+```
+def ortho(data, voxel, fig=None, **kwargs):
+ """Simple orthographic plot of a 3D array using matplotlib.
+
+ :arg data: 3D numpy array
+ :arg voxel: XYZ coordinates for each slice
+ :arg fig: Existing figure and axes for overlay plotting
+
+ All other arguments are passed through to the `imshow` function.
+
+ :returns: The figure and axes (which can be passed back in as the
+ `fig` argument to plot overlays).
+ """
+ x, y, z = voxel
+ xslice = np.flipud(data[x, :, :].T)
+ yslice = np.flipud(data[:, y, :].T)
+ zslice = np.flipud(data[:, :, z].T)
+
+ if fig is None:
+ fig = plt.figure()
+ xax = fig.add_subplot(1, 3, 1)
+ yax = fig.add_subplot(1, 3, 2)
+ zax = fig.add_subplot(1, 3, 3)
+ else:
+ fig, xax, yax, zax = fig
+
+ xax.imshow(xslice, **kwargs)
+ yax.imshow(yslice, **kwargs)
+ zax.imshow(zslice, **kwargs)
+
+ for ax in (xax, yax, zax):
+ ax.set_xticks([])
+ ax.set_yticks([])
+ fig.tight_layout(pad=0)
+
+ return (fig, xax, yax, zax)
+```
+
+
+And another function which uses FSLeyes for more complex plots:
+
+
+```
+import shlex
+import IPython.display as display
+from fsleyes.render import main
+
+def render(cmdline):
+ prefix = '-of screenshot.png -hl -c 2 '
+ main(shlex.split(prefix + cmdline))
+ return display.Image('screenshot.png')
+```
<a class="anchor" id="the-image-class-and-other-data-types"></a>
@@ -56,7 +131,6 @@ It's easy to create an `Image` - you can create one from a file name:
```
-import os.path as op
from fsl.data.image import Image
stddir = op.expandvars('${FSLDIR}/data/standard/')
@@ -73,8 +147,6 @@ You can create an `Image` from an existing `nibabel` image:
```
-import nibabel as nib
-
# load a nibabel image, and
# convert it into an FSL image
nibimg = nib.load(op.join(stddir, 'MNI152_T1_1mm.nii.gz'))
@@ -86,8 +158,6 @@ Or you can create an `Image` from a `numpy` array:
```
-import numpy as np
-
data = np.zeros((100, 100, 100))
img = Image(data, xform=np.eye(4))
```
@@ -123,7 +193,7 @@ and a number of useful methods:
```
std2mm = Image(op.join(stddir, 'MNI152_T1_2mm'))
-mask2mm = Image(op.join(stddir, 'MNI152_T1_2mm_mask'))
+mask2mm = Image(op.join(stddir, 'MNI152_T1_2mm_brain_mask'))
print(std1mm.sameSpace(std2mm))
print(std2mm.sameSpace(mask2mm))
@@ -151,7 +221,7 @@ You can get the image data as a `numpy` array via the `data` attribute:
```
data = std2mm.data
-print(data.min, data.max())
+print(data.min(), data.max())
```
> Note that this will give you the data in its underlying type, unlike the
@@ -163,7 +233,7 @@ You can also read and write data directly via the `Image` object:
```
slc = std2mm[:, :, 45]
-std2mm[0:10, :, :] = 0
+std2mm[0:10, :, :] *= 2
```
@@ -228,6 +298,257 @@ corresponding `nibabel` types:
> in order to use this class.
+<a class="anchor" id="nifti-coordinate-systems"></a>
+### NIfTI coordinate systems
+
+
+The `Image.getAffine` method gives you access to affine transformations which
+can be used to convert coordinates between the different coordinate systems
+associated with a NIfTI image. Have some MNI coordinates you'd like to convert
+to voxels? Easy!
+
+
+```
+stddir = op.expandvars('${FSLDIR}/data/standard/')
+std2mm = Image(op.join(stddir, 'MNI152_T1_2mm'))
+
+mnicoords = np.array([[0, 0, 0],
+ [0, -18, 18]])
+
+world2vox = std2mm.getAffine('world', 'voxel')
+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:
+from fsl.transform.affine import transform
+voxcoords = transform(mnicoords, world2vox)
+
+# just to double check, let's transform
+# those voxel coordinates back into world
+# coordinates
+backtomni = transform(voxcoords, vox2world)
+
+for m, v, b in zip(mnicoords, voxcoords, backtomni):
+ print(m, '->', v, '->', b)
+```
+
+
+> The `Image.getAffine` method can give you transformation matrices
+> between any of these coordinate systems:
+>
+> - `'voxel'`: Image data voxel coordinates
+> - `'world'`: mm coordinates, defined by the sform/qform of an image
+> - `'fsl'`: The FSL coordinate system, used internally by many FSL tools
+> (e.g. FLIRT)
+
+
+Oh, that example was too easy I hear you say? Try this one on for size. Let's
+say we have run FEAT on some task fMRI data, and want to get the MNI
+coordinates of the voxel with peak activation.
+
+
+> This is what people used to use `Featquery` for, back in the un-enlightened
+> days.
+
+
+Let's start by identifying the voxel with the biggest t-statistic:
+
+
+```
+featdir = op.join(op.join('08_fslpy', 'fmri.feat'))
+
+tstat1 = Image(op.join(featdir, 'stats', 'tstat1')).data
+
+# Recall from the numpy practical that
+# argmax gives us a 1D index into a
+# flattened view of the array. We can
+# use the unravel_index function to
+# convert it into a 3D index.
+peakvox = np.abs(tstat1).argmax()
+peakvox = np.unravel_index(peakvox, tstat1.shape)
+print('Peak voxel coordinates for tstat1:', peakvox, tstat1[peakvox])
+```
+
+
+Now that we've got the voxel coordinates in functional space, we need to
+transform them into MNI space. FEAT provides a transformation which goes
+directly from functional to standard space, in the `reg` directory:
+
+
+```
+func2std = np.loadtxt(op.join(featdir, 'reg', 'example_func2standard.mat'))
+```
+
+
+But ... wait a minute ... this is a FLIRT matrix! We can't just plug voxel
+coordinates into a FLIRT matrix and expect to get sensible results, because
+FLIRT works in an internal FSL coordinate system, which is not quite
+`'voxel'`, and not quite `'world'`. So we need to do a little more work.
+Let's start by loading our functional image, and the MNI152 template (the
+source and reference images of our FLIRT matrix):
+
+
+```
+func = Image(op.join(featdir, 'reg', 'example_func'))
+std = Image(op.expandvars(op.join('$FSLDIR', 'data', 'standard', 'MNI152_T1_2mm')))
+```
+
+
+Now we can use them to get affines which convert between all of the different
+coordinate systems - we're going to combine them into a single uber-affine,
+which transforms our functional-space voxels into MNI world coordinates via:
+
+ 1. functional voxels -> FLIRT source space
+ 2. FLIRT source space -> FLIRT reference space
+ 3. FLIRT referece space -> MNI world coordinates
+
+
+```
+vox2fsl = func.getAffine('voxel', 'fsl')
+fsl2mni = std .getAffine('fsl', 'world')
+```
+
+
+Combining two affines into one is just a simple dot-product. There is a
+`concat()` function which does this for us, for any number of affines:
+
+
+```
+from fsl.transform.affine import concat
+
+# To combine affines together, we
+# have to list them in reverse -
+# linear algebra is *weird*.
+funcvox2mni = concat(fsl2mni, func2std, vox2fsl)
+```
+
+
+So we've now got some voxel coordinates from our functional data, and an
+affine to transform into MNI world coordinates. The rest is easy:
+
+
+```
+mnicoords = transform(peakvox, funcvox2mni)
+mnivoxels = transform(mnicoords, std.getAffine('world', 'voxel'))
+mnivoxels = [int(round(v)) for v in mnivoxels]
+print('Peak activation (MNI coordinates):', mnicoords)
+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).
+
+
+<a class="anchor" id="image-processing"></a>
+### Image processing
+
+
+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.utils.image.resample`](https://users.fmrib.ox.ac.uk/~paulmc/fsleyes/fslpy/latest/fsl.utils.image.resample.html)
+module:
+
+
+```
+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')))
+
+# 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')
+
+# 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)
+```
+
+
+Now that we have our t-statistic image in MNI152 space, we can plot it in
+standard space using `matplotlib`:
+
+
+```
+stddir = op.expandvars('${FSLDIR}/data/standard/')
+std2mm = Image(op.join(stddir, 'MNI152_T1_2mm'))
+
+std_tstat1 = std_tstat1.data
+std_tstat1 = np.ma.masked_where(std_tstat1 < 3, std_tstat1)
+
+fig = ortho(std2mm, mnivoxels, cmap=plt.cm.gray)
+fig = ortho(std_tstat1, mnivoxels, cmap=plt.cm.inferno, fig=fig)
+```
+
+
+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.
+
+
+<a class="anchor" id="the-filetree"></a>
+## The `filetree`
+
+
+<a class="anchor" id="fsl-wrapper-functions"></a>
+## FSL wrapper functions
+
+
+The
+[fsl.wrappers](https://users.fmrib.ox.ac.uk/~paulmc/fsleyes/fslpy/latest/fsl.wrappers.html)
+package is the home of "wrapper" functions for a range of FSL tools. You can
+use them to call an FSL tool from Python code, without having to worry about
+constructing a command-line, or saving/loading input/output images.
+
+
+You can use the FSL wrapper functions with file names:
+
+```
+from fsl.wrappers import bet, robustfov, LOAD
+os.chdir('08_fslpy')
+robustfov('bighead', 'bighead_cropped')
+render('bighead bighead_cropped -cm blue')
+```
+
+Or, if you have images in memory, you can pass the `Image` objects in,
+and have the results automatically loaded in too:
+
+```
+cropped = Image('bighead_cropped')
+betted = bet(cropped, LOAD)['output']
+
+fig = ortho(cropped, (80, 112, 85))
+fig = ortho(betted, (80, 112, 85), fig=fig, cmap=plt.cm.inferno)
+```
+
+
<a class="anchor" id="fsl-atlases"></a>
## FSL atlases
@@ -248,6 +569,10 @@ atlases.rescanAtlases()
```
+<a class="anchor" id="querying-atlases"></a>
+### Querying atlases
+
+
You can list all of the available atlases using `listAtlases`:
@@ -295,6 +620,10 @@ print(desc.find(value=48))
```
+<a class="anchor" id="loading-atlas-images"></a>
+### Loading atlas images
+
+
The `loadAtlas` function can be used to load the atlas image:
@@ -322,6 +651,48 @@ print(probatlas)
```
+<a class="anchor" id="working-with-atlases"></a>
+### Working with atlases
+
+
+Both `LabelAtlas` and `ProbabilisticAtlas` objects have a method called `get`,
+which can be used to extract ROI images for a specific region:
+
+
+```
+stddir = op.expandvars('${FSLDIR}/data/standard/')
+std2mm = Image(op.join(stddir, 'MNI152_T1_2mm'))
+
+frontal = lblatlas.get(name='Frontal Pole').data
+frontal = np.ma.masked_where(frontal < 1, frontal)
+
+fig = ortho(std2mm, (45, 54, 45), cmap=plt.cm.gray)
+fig = ortho(frontal, (45, 54, 45), cmap=plt.cm.winter, fig=fig)
+```
+
+
+Calling `get` on a :meth:`ProbabilisticAtlas` will return a probability image:
+
+
+```
+stddir = op.expandvars('${FSLDIR}/data/standard/')
+std2mm = Image(op.join(stddir, 'MNI152_T1_2mm'))
+
+frontal = probatlas.get(name='Frontal Pole')
+frontal = np.ma.masked_where(frontal < 1, frontal)
+
+fig = ortho(std2mm, (45, 54, 45), cmap=plt.cm.gray)
+fig = ortho(frontal, (45, 54, 45), cmap=plt.cm.inferno, fig=fig)
+```
+
+
+The `get` method can be used to retrieve an image for a region by:
+- an `AtlasLabel` object
+- The region index
+- The region value
+- The region name
+
+
`LabelAtlas` objects have a method called `label`, which can be used to
interrogate the atlas at specific locations:
@@ -365,174 +736,3 @@ for idx, val in enumerate(vals):
lbl = probatlas.find(index=idx)
print(' {} [{}]: {:0.2f}%'.format(lbl.value, lbl.name, val))
```
-
-
-<a class="anchor" id="the-filetree"></a>
-## The `filetree`
-
-<a class="anchor" id="nifti-coordinate-systems"></a>
-## NIfTI coordinate systems
-
-<a class="anchor" id="image-processing"></a>
-## Image processing
-
-<a class="anchor" id="fsl-wrapper-functions"></a>
-## FSL wrapper functions
-
-
-
-<a class="anchor" id="nifti-coordinate-systems"></a>
-## NIfTI coordinate systems
-
-
-
-The `getAffine` method gives you access to affine transformations which can be
-used to convert coordinates between the different coordinate systems
-associated with an image. Have some MNI coordinates you'd like to convert to
-voxels? Easy!
-
-```
-mnicoords = np.array([[0, 0, 0],
- [0, -18, 18]])
-
-world2vox = img.getAffine('world', 'voxel')
-vox2world = img.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:
-from fsl.transform.affine import transform
-voxcoords = transform(mnicoords, world2vox)
-
-# just to double check, let's transform
-# those voxel coordinates back into world
-# coordinates
-backtomni = transform(voxcoords, vox2world)
-
-for m, v, b in zip(mnicoords, voxcoords, backtomni):
- print(m, '->', v, '->', b)
-```
-
-> The `Image.getAffine` method can give you transformation matrices
-> between any of these coordinate systems:
->
-> - `'voxel'`: Image data voxel coordinates
-> - `'world'`: mm coordinates, defined by the sform/qform of an image
-> - `'fsl'`: The FSL coordinate system, used internally by many FSL tools
-> (e.g. FLIRT)
-
-
-
-Oh, that example was too easy I hear you say? Try this one on for size. Let's
-say we have run FEAT on some task fMRI data, and want to get the MNI
-coordinates of the voxel with peak activation.
-
-> This is what people used to use `Featquery` for, back in the un-enlightened
-> days.
-
-
-Let's start by identifying the voxel with the biggest t-statistic:
-
-```
-featdir = op.join(op.join('05_nifti', 'fmri.feat'))
-
-# The Image.data attribute returns a
-# numpy array containing, well, the
-# image data.
-tstat1 = Image(op.join(featdir, 'stats', 'tstat1')).data
-
-# Recall from the numpy practical that
-# argmax gives us a 1D index into a
-# flattened view of the array. We can
-# use the unravel_index function to
-# convert it into a 3D index.
-peakvox = np.abs(tstat1).argmax()
-peakvox = np.unravel_index(peakvox, tstat1.shape)
-print('Peak voxel coordinates for tstat1:', peakvox, tstat1[peakvox])
-```
-
-
-Now that we've got the voxel coordinates in functional space, we need to
-transform them into MNI space. FEAT provides a transformation which goes
-directly from functional to standard space, in the `reg` directory:
-
-```
-func2std = np.loadtxt(op.join(featdir, 'reg', 'example_func2standard.mat'))
-```
-
-But ... wait a minute ... this is a FLIRT matrix! We can't just plug voxel
-coordinates into a FLIRT matrix and expect to get sensible results, because
-FLIRT works in an internal FSL coordinate system, which is not quite
-`'voxel'`, and not quite `'world'`. So we need to do a little more work.
-Let's start by loading our functional image, and the MNI152 template (the
-source and reference images of our FLIRT matrix):
-
-```
-func = Image(op.join(featdir, 'reg', 'example_func'))
-std = Image(op.expandvars(op.join('$FSLDIR', 'data', 'standard', 'MNI152_T1_2mm')))
-```
-
-
-Now we can use them to get affines which convert between all of the different
-coordinate systems - we're going to combine them into a single uber-affine,
-which transforms our functional-space voxels into MNI world coordinates via:
-
- 1. functional voxels -> FLIRT source space
- 2. FLIRT source space -> FLIRT reference space
- 3. FLIRT referece space -> MNI world coordinates
-
-
-```
-vox2fsl = func.getAffine('voxel', 'fsl')
-fsl2mni = std .getAffine('fsl', 'world')
-```
-
-Combining two affines into one is just a simple dot-product. There is a
-`concat()` function which does this for us, for any number of affines:
-
-```
-from fsl.transform.affine import concat
-
-# To combine affines together, we
-# have to list them in reverse -
-# linear algebra is *weird*.
-funcvox2mni = concat(fsl2mni, func2std, vox2fsl)
-```
-
-So we've now got some voxel coordinates from our functional data, and an affine
-to transform into MNI world coordinates. The rest is easy:
-
-```
-mnicoords = transform(peakvox, funcvox2mni)
-print('Peak activation (MNI coordinates):', mnicoords)
-```
-
-
-> 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 ;).
-
-
-<a class="anchor" id="image-processing"></a>
-## Image processing
-
-Now, it's all well and good to look at t-statistric 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
-
-```
-from IPython.display import Image as Screenshot
-!fsleyes render -of screenshot.png -std
-```
-
-
-### (Advanced) Transform coordinates with nonlinear warpfields
-
-
-have to use your own dataset
diff --git a/advanced_topics/08_fslpy/bighead.nii.gz b/advanced_topics/08_fslpy/bighead.nii.gz
new file mode 100755
index 0000000000000000000000000000000000000000..6aff1e3c77ded97c81bacb06cd1d84b10220d6ed
Binary files /dev/null and b/advanced_topics/08_fslpy/bighead.nii.gz differ
diff --git a/advanced_topics/08_fslpy/fmri.feat/reg/example_func.nii.gz b/advanced_topics/08_fslpy/fmri.feat/reg/example_func.nii.gz
new file mode 100755
index 0000000000000000000000000000000000000000..68f581786315d538a87884115620bec392ff9d4f
Binary files /dev/null and b/advanced_topics/08_fslpy/fmri.feat/reg/example_func.nii.gz differ
diff --git a/advanced_topics/08_fslpy/fmri.feat/reg/example_func2standard.mat b/advanced_topics/08_fslpy/fmri.feat/reg/example_func2standard.mat
new file mode 100755
index 0000000000000000000000000000000000000000..e0c43af6147ac26c51e16077d5efbedb7196c219
--- /dev/null
+++ b/advanced_topics/08_fslpy/fmri.feat/reg/example_func2standard.mat
@@ -0,0 +1,4 @@
+1.057622308 0.05073972589 0.008769375553 -31.51452272
+-0.0541050194 0.9680196522 0.1445326796 2.872941273
+-0.01020603009 -0.2324151706 1.127557283 21.35031106
+0 0 0 1
diff --git a/advanced_topics/08_fslpy/fmri.feat/stats/tstat1.nii.gz b/advanced_topics/08_fslpy/fmri.feat/stats/tstat1.nii.gz
new file mode 100755
index 0000000000000000000000000000000000000000..3f09fd681deb1e1123ebc4714f31a3665405be47
Binary files /dev/null and b/advanced_topics/08_fslpy/fmri.feat/stats/tstat1.nii.gz differ