Some new FEAT data sets (including real, but small, NIFTI files). Unit tests

for featimage module finished.

tests/__init__.py

@@ -104,6 +104,15 @@ def cleardir(dir):
         elif op.isdir(f):  shutil.rmtree(f)
 
 
+def random_voxels(shape, nvoxels=1):
+    randVoxels = np.vstack([np.random.randint(0, s, nvoxels) for s in shape[:3]]).T
+
+    if nvoxels == 1:
+        return randVoxels[0]
+    else:
+        return randVoxels
+
+
 def make_random_image(filename, dims=(10, 10, 10), xform=None):
     """Creates a NIFTI1 image with random data, saves and
     returns it.

tests/test_featdesign.py

@@ -30,12 +30,29 @@ Second level analyses:
 
 Second level cope?.feats:
     rm -r logs .files tsplot .ramp.gif
-    
 
 
-`1stlevel_1.feat`
+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
@@ -45,6 +62,7 @@ Second level cope?.feats:
 `1stlevel_2.feat`
 
  - 45 time points
+ - Shape (64, 64, 5)
  - 11 EVs in total:
    - 2 stimulus EVs
    - 2 temporal derivative EVs
@@ -54,6 +72,7 @@ Second level cope?.feats:
 
 `1stlevel_3.feat`
  - 45 time points
+ - Shape (64, 64, 5)
  - 32 EVs in total:
    - 2 stimulus EVs
    - 1 Temporal derivative EV on first
@@ -63,16 +82,24 @@ Second level cope?.feats:
    - 1 Voxelwise confound EV
  - 2 contrasts - one on each stimulus EV
 
-`2ndlevel_1.feat`
- - Three inputs
+`2ndlevel_1.gfeat`
+ - Shape (91, 109, 91)
+ - Three subjects
  - Two copes
  - One main EV - group average
 
-`2ndlevel_1.feat`
- - Three inputs
+`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 
 """
 
 
@@ -93,16 +120,35 @@ 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():
 
-    featdirs = ['1stlevel_1.feat', '1stlevel_2.feat', '1stlevel_3.feat',
-                '2ndlevel_1.gfeat', '2ndlevel_2.gfeat']
-    nevs     = [10, 11, 32, 1, 2]
-    shapes   = [(45, 10), (45, 11), (45, 32), (3, 1), (3, 2)]    
-
-    for featdir, nev, shape in zip(featdirs, nevs, shapes):
+    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)
 
@@ -120,9 +166,11 @@ def test_FEATFSFDesign():
                                  settings=settings,
                                  loadVoxelwiseEVs=False)
 
+        rvox = tests.random_voxels(shape)
+
         assert len(des.getEVs()) == nev
-        assert des.getDesign().shape == shape
-        assert des.getDesign((10, 10, 3)).shape == shape
+        assert des.getDesign().shape == desshape
+        assert des.getDesign(rvox).shape == desshape
 
 
 def test_FEATFSFDesign_firstLevelVoxelwiseEV(seed):
@@ -146,7 +194,7 @@ def test_FEATFSFDesign_firstLevelVoxelwiseEV(seed):
             copes=False,
             zstats=False,
             residuals=False,
-            clusterMasks=False)
+            clustMasks=False)
 
         # Now load the design, and make sure that
         # the voxel EVs are filled correctly
@@ -272,6 +320,28 @@ def test_getFirstLevelEVs_3():
         for k, v in atts.items():
             assert getattr(evs[i], k) == v
 
+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):
+
+        print(i, evs[i])
+
+        assert isinstance(evs[i], evtype)
+        for k, v in atts.items():
+            assert getattr(evs[i], k) == v
+ 
+
 
 def test_getHigherLevelEVs_1():
 
@@ -307,12 +377,9 @@ def test_getHigherLevelEVs_2():
 
 def test_loadDesignMat():
 
-    analyses = ['1stlevel_1.feat', '1stlevel_2.feat', '1stlevel_3.feat',
-                '2ndlevel_1.gfeat', '2ndlevel_2.gfeat']
-    shapes   = [(45, 10), (45, 11), (45, 32), (3, 1), (3, 2)]
-
-    for analysis, shape in zip(analyses, shapes):
-        featdir = op.join(datadir, analysis)
+    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)
 

tests/test_featimage.py

@@ -24,47 +24,73 @@ import fsl.data.featanalysis as featanalysis
 
 datadir    = op.join(op.dirname(__file__), 'testdata', 'test_feat')
 featdirs   = ['1stlevel_1.feat', '1stlevel_2.feat', '1stlevel_2.feat',
+              '1stlevel_realdata.feat',
               '2ndlevel_1.gfeat/cope1.feat', '2ndlevel_1.gfeat/cope2.feat',
-              '2ndlevel_2.gfeat/cope1.feat', '2ndlevel_2.gfeat/cope2.feat']
+              '2ndlevel_2.gfeat/cope1.feat', '2ndlevel_2.gfeat/cope2.feat',
+              '2ndlevel_realdata.gfeat/cope1.feat',
+              '2ndlevel_realdata.gfeat/cope2.feat']
 shapes = [(64,  64,  5, 45),
           (64,  64,  5, 45),
           (64,  64,  5, 45),
+          ( 4,   4,  5, 45),
           (91, 109, 91,  3),
           (91, 109, 91,  3),
           (91, 109, 91,  3),
-          (91, 109, 91,  3)]
+          (91, 109, 91,  3),
+          (10,  10, 10,  3),
+          (10,  10, 10,  3)]
 xforms = [[[-4, 0, 0,    0],
            [ 0, 4, 0,    0],
            [ 0, 0, 6,    0],
            [ 0, 0, 0,    1]]] * 3 + \
+         [[[-4, 0, 0, -120],
+           [ 0, 4, 0,  120],
+           [ 0, 0, 6,    0],
+           [ 0, 0, 0,    1]]] * 1 + \
          [[[-2, 0, 0,   90],
            [ 0, 2, 0, -126],
            [ 0, 0, 2,  -72],
-           [ 0, 0, 0,    1]]] * 4
+           [ 0, 0, 0,    1]]] * 4 + \
+         [[[-2, 0, 0,   20],
+           [ 0, 2, 0, -102],
+           [ 0, 0, 2,   -6],
+           [ 0, 0, 0,    1]]] * 2
 xforms = [np.array(xf) for xf in xforms]
-              
-TEST_ANALYSES = list(zip(
-    featdirs,
-    shapes,
-    xforms))
+
+TEST_ANALYSES = {}
+for i, featdir in enumerate(featdirs):
+    TEST_ANALYSES[featdir] = {
+        'shape' : shapes[i],
+        'xform' : xforms[i]
+    }
 
 
 def test_FEATImage_attributes():
 
-    for i, featdir in enumerate(featdirs):
-        with tests.testdir() as testdir:
+    # TEst bad input
+    with pytest.raises(Exception):
+        featimage.FEATImage('baddir')
+
+    for featdir in TEST_ANALYSES.keys():
 
-            featdir = tests.make_mock_feat_analysis(
-                op.join(datadir, featdir),
-                testdir,
-                shapes[i],
-                xforms[i],
-                pes=False,
-                copes=False,
-                zstats=False,
-                residuals=False,
-                clustMasks=False)
+        shape = TEST_ANALYSES[featdir]['shape']
+        xform = TEST_ANALYSES[featdir]['xform']
+        with tests.testdir() as testdir:
 
+            if 'realdata' not in featdir:
+                featdir = tests.make_mock_feat_analysis(
+                    op.join(datadir, featdir),
+                    testdir,
+                    shape,
+                    xform,
+                    pes=False,
+                    copes=False,
+                    zstats=False,
+                    residuals=False,
+                    clustMasks=False)
+            else:
+                featdir = op.join(datadir, featdir)
+                
             # Now create a FEATImage. We validate its
             # attributes against the values returned by
             # the functions in featdesign/featanalysis.
@@ -75,6 +101,9 @@ def test_FEATImage_attributes():
             evnames  = [ev.title for ev in design.getEVs()]
             contrastnames, contrasts = featanalysis.loadContrasts(featdir)
 
+            assert np.all(np.isclose(fi.shape,         shape))
+            assert np.all(np.isclose(fi.voxToWorldMat, xform))
+
             assert fi.getFEATDir()             == featdir
             assert fi.getAnalysisName()        == op.splitext(op.basename(featdir))[0]
             assert fi.isFirstLevelAnalysis()   == featanalysis.isFirstLevelAnalysis(settings)
@@ -100,16 +129,24 @@ def test_FEATImage_attributes():
 
 def test_FEATImage_imageAccessors():
 
-    for i, featdir in enumerate(featdirs):
+    for featdir in TEST_ANALYSES.keys():
+        
+        shape = TEST_ANALYSES[featdir]['shape']
+        xform = TEST_ANALYSES[featdir]['xform']
+        
         with tests.testdir() as testdir:
-            featdir = tests.make_mock_feat_analysis(
-                op.join(datadir, featdir),
-                testdir,
-                shapes[i],
-                xforms[i])
 
-            shape4D    = shapes[  i]
-            shape      = shape4D[:3]
+            if 'realdata' not in featdir:
+                featdir = tests.make_mock_feat_analysis(
+                    op.join(datadir, featdir),
+                    testdir,
+                    shape,
+                    xform)
+            else:
+                featdir = op.join(datadir, featdir)
+
+            shape4D = shape
+            shape   = shape4D[:3]
 
             fi    = featimage.FEATImage(featdir)
             nevs  = fi.numEVs()
@@ -127,24 +164,114 @@ def test_FEATImage_imageAccessors():
                     assert fi.getClusterMask(con).shape == shape
             
 
+def test_FEATImage_nostats():
+    
+    featdir = op.join(datadir, '1stlevel_nostats.feat')
+    shape   = (4, 4, 5, 45)
+
+    with tests.testdir() as testdir:
+        
+        featdir = tests.make_mock_feat_analysis(featdir, testdir, shape)
+        fi      = featimage.FEATImage(featdir)
+
+        assert fi.getDesign() is None
+        assert fi.numPoints() == 0
+        assert fi.numEVs()    == 0
+        assert fi.evNames()   == []
+
+        with pytest.raises(Exception):
+            fi.fit([1, 2, 3], (2, 2, 2))
+            
+        with pytest.raises(Exception):
+            fi.partialFit([1, 2, 3], (2, 2, 2)) 
+                    
 
 def test_FEATImage_fit_firstLevel():
-    featdir = ''
-    fi      = featimage.FEATImage(featdir) 
-    fi.fit(0, (0, 0, 0))
+    
+    featdir  = op.join(datadir, '1stlevel_realdata.feat')
+    fi       = featimage.FEATImage(featdir)
+    expect   = np.array([
+        10625.35273455,  10625.35263602,  10625.35248499,  10625.35272602,
+        10625.35286707,  10625.35237145,  10625.35244999,  10625.35270435,
+        10625.35272762,  10629.03397661,  10685.36428581,  10658.64633521,
+        10524.89226543,  10415.61794156,  10373.6671577 ,  10437.10001383,
+        10403.88611746,  10226.98548936,  10080.14323091,  10012.89132265,
+        9936.530395   ,   9957.92598556,  10090.51140821,  10199.6446317 ,
+        10246.67689405,  10261.45133255,  10265.36943466,  10266.22514043,
+        10266.32736048,  10264.92716455,  10243.05011597,  10245.62798475,
+        10287.91883737,  10325.38456267,  10341.92299781,  10347.17916861,
+        10348.58339616,  10348.89634025,  10348.93522057,  10345.25397481,
+        10288.9236822 ,  10315.64160242,  10449.39567496,  10558.66999883,
+        10597.64918744]) 
+    
+    # bad contrast
+    with pytest.raises(Exception):
+        fi.fit([1, 2, 3, 4, 5, 6, 7], (2, 2, 2))
+
+    # bad voxel
+    with pytest.raises(Exception):
+        fi.fit([1, 0, 0, 0], (6, 7, 7))
+
+    result = fi.fit([1, 1, 1, 1], (2, 2, 2))
+    assert np.all(np.isclose(result, expect))
 
 
 def test_FEATImage_fit_higherLevel():
-    featdir = ''
-    fi      = featimage.FEATImage(featdir) 
-    fi.fit(0, (0, 0, 0)) 
+    
+    featdir = op.join(datadir, '2ndlevel_realdata.gfeat/cope1.feat')
+    fi      = featimage.FEATImage(featdir)
+    expect  = np.array([86.37929535, 86.37929535, 86.37929535])
+    result = fi.fit([1], (5, 5, 5))
+    
+    assert np.all(np.isclose(result, expect))
 
 
 def test_FEATImage_partialFit():
-    featdir = ''
-    fi      = featimage.FEATImage(featdir) 
-    fi.fit(0, (0, 0, 0))
+    
+    featdir = op.join(datadir, '1stlevel_realdata.feat')
+    fi      = featimage.FEATImage(featdir)
+    expect  = np.array([
+        10476.23185443,  10825.13542716,  10781.15297632,  11032.96902851,
+        11019.47536463,  10856.61537194,  10463.38391362,  10664.79227954,
+        10948.11850521,  10712.01568132,  10882.88149773,  10745.65913733,
+        10590.78057109,  10524.89948148,  10744.82941967,  10422.96359453,
+        10156.39446402,  10219.8339014 ,  10115.32145738,  10494.28109315,
+        10121.89555309,  10165.92560409,  10556.70058668,  10198.67478569,
+        10045.04934583,   9906.60233353,  10015.75569565,   9864.38786016,
+        10241.91219554,  10099.08538293,  10444.7826294 ,  10152.72622847,
+        10077.09075815,  10128.63320464,  10087.93203101,  10450.29667654,
+        10207.89144059,  10137.98334586,  10387.09965691,  10194.79436483,
+        10203.21032619,  10136.1942605 ,  10128.23728873,  10416.78984136,
+        10118.51262128])
+    result = fi.partialFit([1, 1, 1, 1], (2, 2, 2))
+    
+    assert np.all(np.isclose(result, expect))
 
 
-def test_FEATImage_nostats():
-    pass 
+def test_modelFit(seed):
+
+    for i in range(500):
+        
+        # 2 evs, 20 timepoints
+        # First EV is a boxcar,
+        # second is a random regressor
+        design       = np.zeros((20, 2))
+        design[:, 0] = np.tile([1, 1, 0, 0], 5)
+        design[:, 1] = np.random.random(20)
+
+        # demean columns of design matrix
+        for ev in range(design.shape[1]):
+            design[:, ev] = design[:, ev] - design[:, ev].mean()
+
+        # Generate some random PEs, and
+        # generate the data that would 
+        # have resulted in them
+        pes      = np.random.random(2)
+        expect   = np.dot(design, pes)
+        contrast = [1] * design.shape[1]
+        
+        result1 = featimage.modelFit(expect, design, contrast, pes, True)
+        result2 = featimage.modelFit(expect, design, contrast, pes, False)
+
+        assert np.all(np.isclose(result1 - expect.mean(), expect))
+        assert np.all(np.isclose(result2,                 expect))

tests/testdata/test_feat/1stlevel_nostats.feat/absbrainthresh.txt

+1000.0

tests/testdata/test_feat/1stlevel_nostats.feat/design.con

+/ContrastName1	auditory 
+/ContrastName2	visual 
+/NumWaves	4
+/NumContrasts	2
+/PPheights		1.139390e+00	1.211980e+00
+/RequiredEffect		2.611	3.237
+
+/Matrix
+1.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 
+0.000000e+00 0.000000e+00 1.000000e+00 0.000000e+00 

tests/testdata/test_feat/1stlevel_nostats.feat/design.frf

+1
+1
+2
+2

tests/testdata/test_feat/1stlevel_nostats.feat/design.fsf

+
+# FEAT version number
+set fmri(version) 6.00
+
+# Are we in MELODIC?
+set fmri(inmelodic) 0
+
+# Analysis level
+# 1 : First-level analysis
+# 2 : Higher-level analysis
+set fmri(level) 1
+
+# Which stages to run
+# 0 : No first-level analysis (registration and/or group stats only)
+# 7 : Full first-level analysis
+# 1 : Pre-processing
+# 2 : Statistics
+set fmri(analysis) 1
+
+# Use relative filenames
+set fmri(relative_yn) 0
+
+# Balloon help
+set fmri(help_yn) 1
+
+# Run Featwatcher
+set fmri(featwatcher_yn) 1
+
+# Cleanup first-level standard-space images
+set fmri(sscleanup_yn) 0
+
+# Output directory
+set fmri(outputdir) "/Users/paulmc/Projects/fslpy/tests/testdata/test_feat/1stlevel_nostats"
+
+# TR(s)
+set fmri(tr) 3.000000 
+
+# Total volumes
+set fmri(npts) 45
+
+# Delete volumes
+set fmri(ndelete) 0
+
+# Perfusion tag/control order
+set fmri(tagfirst) 1
+
+# Number of first-level analyses
+set fmri(multiple) 1
+
+# Higher-level input type
+# 1 : Inputs are lower-level FEAT directories
+# 2 : Inputs are cope images from FEAT directories
+set fmri(inputtype) 2
+
+# Carry out pre-stats processing?
+set fmri(filtering_yn) 1
+
+# Brain/background threshold, %
+set fmri(brain_thresh) 10
+
+# Critical z for design efficiency calculation
+set fmri(critical_z) 5.3
+
+# Noise level
+set fmri(noise) 0.66
+
+# Noise AR(1)
+set fmri(noisear) 0.34
+
+# Motion correction
+# 0 : None
+# 1 : MCFLIRT
+set fmri(mc) 1
+
+# Spin-history (currently obsolete)
+set fmri(sh_yn) 0
+
+# B0 fieldmap unwarping?
+set fmri(regunwarp_yn) 0
+
+# EPI dwell time (ms)
+set fmri(dwell) 0.7
+
+# EPI TE (ms)
+set fmri(te) 35
+
+# % Signal loss threshold
+set fmri(signallossthresh) 10
+
+# Unwarp direction
+set fmri(unwarp_dir) y-
+
+# Slice timing correction
+# 0 : None
+# 1 : Regular up (0, 1, 2, 3, ...)
+# 2 : Regular down
+# 3 : Use slice order file
+# 4 : Use slice timings file
+# 5 : Interleaved (0, 2, 4 ... 1, 3, 5 ... )
+set fmri(st) 0
+
+# Slice timings file
+set fmri(st_file) ""
+
+# BET brain extraction
+set fmri(bet_yn) 1
+
+# Spatial smoothing FWHM (mm)
+set fmri(smooth) 5
+
+# Intensity normalization
+set fmri(norm_yn) 0
+
+# Perfusion subtraction
+set fmri(perfsub_yn) 0
+
+# Highpass temporal filtering
+set fmri(temphp_yn) 1
+
+# Lowpass temporal filtering
+set fmri(templp_yn) 0
+
+# MELODIC ICA data exploration
+set fmri(melodic_yn) 0
+
+# Carry out main stats?
+set fmri(stats_yn) 0
+
+# Carry out prewhitening?
+set fmri(prewhiten_yn) 1
+
+# Add motion parameters to model