Skip to content
Snippets Groups Projects
Forked from FSL / fslpy
1493 commits behind the upstream repository.
test_featdesign.py 14.34 KiB
#!/usr/bin/env python
#
# test_featdesign.py -
#
# Author: Paul McCarthy <pauldmccarthy@gmail.com>
#
"""Test data sets (in testdata/test_feat) were generated with FSL 5.0.9, and
then 'cleaned' to remove unnecessary files and reduce size, with the following
commands:

First level analyses:

    find . -name "*png" -delete
    find . -name "*ppm" -delete
    find . -name "*html" -delete
    rm -r logs .files tsplot .ramp.gif
    for f in `find . -name "*nii.gz"`; do
        echo $f > $f
    done

Second level analyses:

    find . -name "*png" -delete
    find . -name "*ppm" -delete
    find . -name "*html" -delete
    rm -r logs .files  .ramp.gif
    for f in `find . -name "*nii.gz"`; do
        echo $f > $f
    done

Second level cope?.feats:
    rm -r logs .files tsplot .ramp.gif


Two data sets - `1stlevel_realdata.feat` and `2ndlevel_realdata.gfeat` -
still have real (cut-down) NIFTI images, and have only been cleaned
with the following commands:

    find . -name "*png" -delete
    find . -name "*ppm" -delete
    find . -name "*html" -delete
    rm -r logs .files tsplot .ramp.gif


`1stlevel_realdata.feat`
 - 45 time points
 - Shape (4, 4, 5)
 - 4 evs in total:
   - 2 stimulus EVs
   - 2 temporal derivative EVs
 - 2 contrasts - one on each stimulus EV

`1stlevel_1.feat`
 - 45 time points
 - Shape (64, 64, 5)
 - 10 EVs in total:
   - 2 stimulus EVs
   - 2 temporal derivative EVs
   - 6 Standard motion parameters
 - 2 contrasts - one on each stimulus EV

`1stlevel_2.feat`

 - 45 time points
 - Shape (64, 64, 5)
 - 11 EVs in total:
   - 2 stimulus EVs
   - 2 temporal derivative EVs
   - 1 voxelwise EV
   - 6 Standard motion parameters
 - 2 contrasts - one on each stimulus EV

`1stlevel_3.feat`
 - 45 time points
 - Shape (64, 64, 5)
 - 32 EVs in total:
   - 2 stimulus EVs
   - 1 Temporal derivative EV on first
   - 2 gamma basis functions on second
   - 24 Standard+extended motion parameters
   - 2 Confound EVs
   - 1 Voxelwise confound EV
 - 2 contrasts - one on each stimulus EV

`2ndlevel_1.gfeat`
 - Shape (91, 109, 91)
 - Three subjects
 - Two copes
 - One main EV - group average

`2ndlevel_2.gfeat`
 - Shape (91, 109, 91)
 - Three subjects
 - Two copes
 - One main EV - group average
 - One voxelwise EV

`2ndlevel_realdata.gfeat`
 - Shape (10, 10, 10)
 - Three subjects
 - Two copes
 - One main EV - group average
"""


import os.path as op
import numpy   as np

import pytest

import tests
from   fsl.utils.tempdir     import tempdir
import fsl.data.featdesign   as     featdesign
import fsl.data.featanalysis as     featanalysis


datadir = op.join(op.dirname(__file__), 'testdata', 'test_feat')

featdirs = ['1stlevel_1.feat', '1stlevel_2.feat', '1stlevel_3.feat',
            '1stlevel_realdata.feat',
            '2ndlevel_1.gfeat', '2ndlevel_2.gfeat', '2ndlevel_realdata.gfeat']
nevs     = [10, 11, 32, 4, 1, 2, 1]
shapes   = [(64,  64,  5),
            (64,  64,  5),
            (64,  64,  5),
            ( 4,   4,  5),
            (91, 109, 91),
            (91, 109, 91),
            (10,  10, 10)]
designShapes = [(45, 10), (45, 11), (45, 32), (45, 4), (3, 1), (3, 2), (3, 1)]

TEST_ANALYSES = {}
for i, featdir in enumerate(featdirs):
    TEST_ANALYSES[featdir] = {
        'nevs'        : nevs[i],
        'shape'       : shapes[i],
        'designShape' : designShapes[i]
    }


def test_FEATFSFDesign():

    for featdir in TEST_ANALYSES.keys():

        nev      = TEST_ANALYSES[featdir]['nevs']
        shape    = TEST_ANALYSES[featdir]['shape']
        desshape = TEST_ANALYSES[featdir]['designShape']
        featdir  = op.join(datadir, featdir)
        settings = featanalysis.loadSettings(featdir)

        # We can't load the voxelwise EVs
        # here, because all of the .nii.gz
        # files in the test directory are
        # stubs. Voxelwise EVs get tested
        # in the test_FEATFSFDesign_firstLevelVoxelwiseEV
        # function
        des = featdesign.FEATFSFDesign(featdir,
                                       loadVoxelwiseEVs=False)

        # Can also specify the design.fsf settings
        featdesign.FEATFSFDesign(featdir,
                                 settings=settings,
                                 loadVoxelwiseEVs=False)

        rvox = tests.random_voxels(shape)

        assert len(des.getEVs()) == nev
        assert des.getDesign().shape == desshape
        assert des.getDesign(rvox).shape == desshape

        del des
        des = None


def test_FEATFSFDesign_firstLevelVoxelwiseEV(seed):

    template = op.join(datadir, '1stlevel_2.feat')

    with tests.testdir() as testdir:

        # Set up some test data - we make
        # a copy of the feat directory,
        # and generate some dummy data for
        # the voxelwise EVs.
        featdir    = op.join(testdir, '1stlevel_2.feat')
        shape4D    = (64, 64, 5, 45)
        shape      = shape4D[:3]

        featdir = tests.make_mock_feat_analysis(
            template, testdir, shape4D,
            indata=False,
            pes=False,
            copes=False,
            zstats=False,
            residuals=False,
            clustMasks=False)

        # Now load the design, and make sure that
        # the voxel EVs are filled correctly
        design    = featdesign.FEATFSFDesign(featdir)
        voxevIdxs = [ev.index for ev in design.getEVs()
                     if isinstance(ev, (featdesign.VoxelwiseEV,
                                        featdesign.VoxelwiseConfoundEV))]

        randVoxels = np.vstack([np.random.randint(0, s, 10) for s in shape]).T

        for voxel in randVoxels:

            voxel  = [int(v) for v in voxel]
            offset = np.ravel_multi_index(voxel, shape)
            matrix = design.getDesign(voxel)

            for i, evidx in enumerate(voxevIdxs):
                expect = np.arange(i, i + 45) + offset
                assert np.all(np.isclose(matrix[:, evidx], expect))
        del design
        design = None


def test_getFirstLevelEVs_1():
    featdir  = op.join(datadir, '1stlevel_1.feat')
    settings = featanalysis.loadSettings(featdir)
    matrix   = featdesign.loadDesignMat(op.join(featdir, 'design.mat'))

    expected = [(featdesign.NormalEV,             {'index' : 0, 'origIndex' : 0}),
                (featdesign.TemporalDerivativeEV, {'index' : 1}),
                (featdesign.NormalEV,             {'index' : 2, 'origIndex' : 1}),
                (featdesign.TemporalDerivativeEV, {'index' : 3}),
                (featdesign.MotionParameterEV,    {'index' : 4, 'motionIndex' : 0}),
                (featdesign.MotionParameterEV,    {'index' : 5, 'motionIndex' : 1}),
                (featdesign.MotionParameterEV,    {'index' : 6, 'motionIndex' : 2}),
                (featdesign.MotionParameterEV,    {'index' : 7, 'motionIndex' : 3}),
                (featdesign.MotionParameterEV,    {'index' : 8, 'motionIndex' : 4}),
                (featdesign.MotionParameterEV,    {'index' : 9, 'motionIndex' : 5})]

    evs = featdesign.getFirstLevelEVs(featdir, settings, matrix)
    assert len(evs) == 10

    for i, (evtype, atts) in enumerate(expected):

        assert isinstance(evs[i], evtype)
        for k, v in atts.items():
            assert getattr(evs[i], k) == v
    del evs
    evs = None


def test_getFirstLevelEVs_2():
    featdir  = op.join(datadir, '1stlevel_2.feat')
    settings = featanalysis.loadSettings(featdir)
    matrix   = featdesign.loadDesignMat(op.join(featdir, 'design.mat'))

    expected = [(featdesign.NormalEV,             {'index' : 0, 'origIndex' : 0}),
                (featdesign.TemporalDerivativeEV, {'index' : 1}),
                (featdesign.NormalEV,             {'index' : 2, 'origIndex' : 1}),
                (featdesign.TemporalDerivativeEV, {'index' : 3}),
                (featdesign.VoxelwiseEV,          {'index' : 4, 'origIndex' : 2}),
                (featdesign.MotionParameterEV,    {'index' : 5,  'motionIndex' : 0}),
                (featdesign.MotionParameterEV,    {'index' : 6,  'motionIndex' : 1}),
                (featdesign.MotionParameterEV,    {'index' : 7,  'motionIndex' : 2}),
                (featdesign.MotionParameterEV,    {'index' : 8,  'motionIndex' : 3}),
                (featdesign.MotionParameterEV,    {'index' : 9,  'motionIndex' : 4}),
                (featdesign.MotionParameterEV,    {'index' : 10, 'motionIndex' : 5})]

    evs = featdesign.getFirstLevelEVs(featdir, settings, matrix)
    assert len(evs) == 11

    for i, (evtype, atts) in enumerate(expected):

        assert isinstance(evs[i], evtype)
        for k, v in atts.items():
            assert getattr(evs[i], k) == v
    del evs
    evs = None


def test_getFirstLevelEVs_3():
    featdir  = op.join(datadir, '1stlevel_3.feat')
    settings = featanalysis.loadSettings(featdir)
    matrix   = featdesign.loadDesignMat(op.join(featdir, 'design.mat'))

    expected = [(featdesign.NormalEV,             {'index' : 0, 'origIndex' : 0}),
                (featdesign.TemporalDerivativeEV, {'index' : 1}),
                (featdesign.NormalEV,             {'index' : 2, 'origIndex' : 1}),
                (featdesign.BasisFunctionEV,      {'index' : 3}),
                (featdesign.BasisFunctionEV,      {'index' : 4}),
                (featdesign.VoxelwiseConfoundEV,  {'index' : 5, 'voxIndex'     : 0}),
                (featdesign.MotionParameterEV,    {'index' : 6,  'motionIndex' : 0}),
                (featdesign.MotionParameterEV,    {'index' : 7,  'motionIndex' : 1}),
                (featdesign.MotionParameterEV,    {'index' : 8,  'motionIndex' : 2}),
                (featdesign.MotionParameterEV,    {'index' : 9,  'motionIndex' : 3}),
                (featdesign.MotionParameterEV,    {'index' : 10, 'motionIndex' : 4}),
                (featdesign.MotionParameterEV,    {'index' : 11, 'motionIndex' : 5}),
                (featdesign.MotionParameterEV,    {'index' : 12, 'motionIndex' : 6}),
                (featdesign.MotionParameterEV,    {'index' : 13, 'motionIndex' : 7}),
                (featdesign.MotionParameterEV,    {'index' : 14, 'motionIndex' : 8}),
                (featdesign.MotionParameterEV,    {'index' : 15, 'motionIndex' : 9}),
                (featdesign.MotionParameterEV,    {'index' : 16, 'motionIndex' : 10}),
                (featdesign.MotionParameterEV,    {'index' : 17, 'motionIndex' : 11}),
                (featdesign.MotionParameterEV,    {'index' : 18, 'motionIndex' : 12}),
                (featdesign.MotionParameterEV,    {'index' : 19, 'motionIndex' : 13}),
                (featdesign.MotionParameterEV,    {'index' : 20, 'motionIndex' : 14}),
                (featdesign.MotionParameterEV,    {'index' : 21, 'motionIndex' : 15}),
                (featdesign.MotionParameterEV,    {'index' : 22, 'motionIndex' : 16}),
                (featdesign.MotionParameterEV,    {'index' : 23, 'motionIndex' : 17}),
                (featdesign.MotionParameterEV,    {'index' : 24, 'motionIndex' : 18}),
                (featdesign.MotionParameterEV,    {'index' : 25, 'motionIndex' : 19}),
                (featdesign.MotionParameterEV,    {'index' : 26, 'motionIndex' : 20}),
                (featdesign.MotionParameterEV,    {'index' : 27, 'motionIndex' : 21}),
                (featdesign.MotionParameterEV,    {'index' : 28, 'motionIndex' : 22}),
                (featdesign.MotionParameterEV,    {'index' : 29, 'motionIndex' : 23}),
                (featdesign.ConfoundEV,           {'index' : 30, 'confIndex'   : 0}),
                (featdesign.ConfoundEV,           {'index' : 31, 'confIndex'   : 1})]

    evs = featdesign.getFirstLevelEVs(featdir, settings, matrix)

    assert len(evs) == 32

    for i, (evtype, atts) in enumerate(expected):

        assert isinstance(evs[i], evtype)
        for k, v in atts.items():
            assert getattr(evs[i], k) == v

    del evs
    evs = None


def test_getFirstLevelEVs_realdata():
    featdir  = op.join(datadir, '1stlevel_realdata.feat')
    settings = featanalysis.loadSettings(featdir)
    matrix   = featdesign.loadDesignMat(op.join(featdir, 'design.mat'))
    expected = [(featdesign.NormalEV,             {'index' : 0, 'origIndex' : 0}),
                (featdesign.TemporalDerivativeEV, {'index' : 1}),
                (featdesign.NormalEV,             {'index' : 2, 'origIndex' : 1}),
                (featdesign.TemporalDerivativeEV, {'index' : 3})]

    evs = featdesign.getFirstLevelEVs(featdir, settings, matrix)

    assert len(evs) == 4

    for i, (evtype, atts) in enumerate(expected):

        assert isinstance(evs[i], evtype)
        for k, v in atts.items():
            assert getattr(evs[i], k) == v
    del evs
    evs = None


def test_getHigherLevelEVs_1():

    featdir  = op.join(datadir, '2ndlevel_1.gfeat')
    settings = featanalysis.loadSettings(featdir)
    matrix   = featdesign.loadDesignMat(op.join(featdir, 'design.mat'))

    evs = featdesign.getHigherLevelEVs(featdir, settings, matrix)

    assert len(evs) == 1
    assert isinstance(evs[0], featdesign.NormalEV)
    assert evs[0].index     == 0
    assert evs[0].origIndex == 0
    del evs
    evs = None


def test_getHigherLevelEVs_2():

    featdir  = op.join(datadir, '2ndlevel_2.gfeat')
    settings = featanalysis.loadSettings(featdir)
    matrix   = featdesign.loadDesignMat(op.join(featdir, 'design.mat'))

    evs = featdesign.getHigherLevelEVs(featdir, settings, matrix)

    assert len(evs) == 2
    assert isinstance(evs[0], featdesign.NormalEV)
    assert evs[0].index     == 0
    assert evs[0].origIndex == 0
    assert isinstance(evs[1], featdesign.VoxelwiseEV)
    del evs
    evs = None


def test_loadDesignMat():

    for featdir in TEST_ANALYSES.keys():
        shape   = TEST_ANALYSES[featdir]['designShape']
        featdir = op.join(datadir, featdir)
        fname   = op.join(featdir, 'design.mat')
        mat     = featdesign.loadDesignMat(fname)

        assert mat.shape == shape

    nonfile = op.join(datadir, 'non-existent-file')
    badfile = op.join(datadir, '1stlevel_1.feat', 'design.fsf')
    with pytest.raises(Exception):
        featdesign.loadDesignMat(nonfile)

    with pytest.raises(Exception):
        featdesign.loadDesignMat(badfile)


def test_VoxelwiseEVs():
    with tempdir():
        img = tests.make_random_image('image.nii.gz', (10, 10, 10, 10))

        ev1 = featdesign.VoxelwiseEV(        0, 0, 'ev1', 'image.nii.gz')
        ev2 = featdesign.VoxelwiseConfoundEV(0, 0, 'ev2', 'image.nii.gz')

        for xyz in tests.random_voxels((10, 10, 10), 10):
            x, y, z = map(int, xyz)
            exp = img.dataobj[x, y, z, :]
            assert np.all(ev1.image[x, y, z, :] == exp)
            assert np.all(ev2.image[x, y, z, :] == exp)