Initial import

fsl_course/feat3/feedsRun4.py

+#!/usr/bin/env fslpython
+#
+# A FEEDS test which runs FEAT and then checks the output
+
+import sys
+import os
+import os.path as op
+import utils
+
+# parse args
+
+outDir = sys.argv[1]
+dataDir = sys.argv[2]
+fslDir = os.environ['FSLDIR']
+
+# copy FSF to test location and update paths
+origFSF = op.join(dataDir, "fsl_course_data/fmri/bf/design.fsf")
+origFSLDir = "/opt/fmrib/fsl-alpha"
+origDataDir = "/vols/Data/fsldev/dataSets"
+
+newFSF = op.join(outDir, "design.fsf")
+
+cmd = "./cleanFSF.py {0} {1} {2} {3} {4} {5}".format(origFSF, newFSF, origDataDir, dataDir, outDir, origFSLDir)
+utils.run(cmd)
+
+# run FEAT
+utils.run("unset SGE_ROOT; feat {0}".format(op.join(outDir, "design.fsf")))
+
+sys.exit(0)

fsl_course/feat3/feedsRun5.py

+#!/usr/bin/env fslpython
+#
+# A FEEDS test which runs FEAT and then checks the output
+
+import sys
+import os
+import os.path as op
+import utils
+
+# parse args
+
+outDir = sys.argv[1]
+dataDir = sys.argv[2]
+fslDir = os.environ['FSLDIR']
+
+# copy FSF to test location and update paths
+origFeatDir = op.join(dataDir, "fsl_course_data/fmri/bf/filtdata+.feat")
+origFSF = op.join(origFeatDir, "design.fsf")
+origFSLDir = "/opt/fmrib/fsl-alpha"
+origDataDir = "/vols/Data/fsldev/dataSets"
+
+newFSF = op.join(outDir, "design.fsf")
+
+cmd = "./cleanFSF.py {0} {1} {2} {3} {4} {5}".format(origFSF, newFSF, origDataDir, dataDir, outDir, origFSLDir)
+utils.run(cmd)
+
+# run FEAT
+utils.run("unset SGE_ROOT; feat {0}".format(op.join(outDir, "design.fsf")))
+
+sys.exit(0)

fsl_course/feat3/feedsUtils.py

+#!/usr/bin/env fslpython
+"""Utility functions for use by FEEDS tests. """
+
+
+import os
+import os.path as op
+import uuid
+import tempfile
+import subprocess
+import collections
+
+import numpy as np
+
+
+def fslbin(name):
+    return op.join(os.environ["FSLDIR"], "bin", name)
+
+
+def imrm(imageName):
+    
+    cmd = [fslbin("imrm"), imageName]
+    return subprocess.call(cmd)
+
+
+def readfsf(filename, keys):
+
+    if not isinstance(keys, collections.Sequence):
+        keys = [keys]
+    
+    values = {}
+    
+    with open(filename, 'rt') as f:
+        for line in f:
+            
+            if not line.startswith('set '):
+                continue
+
+            line = line.split()
+            k    = line[1]
+            v    = ' '.join(line[2:])
+
+            if k in keys:
+                values[k] = v
+
+    return [values[key] for key in keys]
+
+
+def writefsf(filename, keys, values):
+    
+    if not isinstance(keys,   collections.Sequence): keys   = [keys]
+    if not isinstance(values, collections.Sequence): values = [values]
+
+    if len(keys) != len(values):
+        raise ValueError('len(keys) != len(values)')
+
+    with open(filename, 'at') as f:
+        for k, v in zip(keys, values):
+
+            f.write('\nset {} {}'.format(k, v))
+
+
+def distance(coord1, coord2):
+    """Returns the euclidean distance between the provided coordinates. """
+    coord1 = np.array(coord1)
+    coord2 = np.array(coord2)
+
+    return np.sqrt(np.sum(coord1 * coord2))
+
+
+def normalisedImageError(image1, image2):
+
+    tempFile = op.join(tempfile.gettempdir(), str(uuid.uuid4()))
+
+    cmd = [fslbin("fslmaths"),
+           image1,
+           "-sub", image2,
+           "-div", image1,
+           "-abs",
+           tempFile]
+    
+    subprocess.check_call(cmd, stderr=subprocess.STDOUT)
+
+    cmd    = [fslbin("fslstats"), tempFile, "-m"]
+    result = float(subprocess.check_output(cmd))
+
+    imrm(tempFile)
+    
+    return result
+
+
+def unscaledImageError(image1, image2):
+
+    cmd = [fslbin("fslstats"),
+           image1,
+           "-d", image2,
+           "-a",
+           "-m"]
+    
+    return float(subprocess.check_output(cmd))
+
+
+def testImage(image1, image2, tolerance, normalised=True):
+    
+    if normalised: error = normalisedImageError(image1, image2)
+    else:          error = unscaledImageError(  image1, image2)
+
+    return (error <= tolerance, error)
+
+
+def testNumber(num1, num2, tolerance, normalised=False):
+
+    num1 = np.array(num1, dtype=np.float)
+    num2 = np.array(num2, dtype=np.float)
+
+    if normalised: error = np.abs(num1 - num2) / num1
+    else:          error = np.abs(num1 - num2)
+
+    if error.size > 1:
+        error = error.mean()
+
+    return (error <= tolerance, error)
+       
+
+def worstError(errors, strict):
+
+    results, errors = zip(*errors)
+
+    if   not all(results):                        return 1
+    elif strict and any([e > 0 for e in errors]): return 1
+    else:                                         return 0

fsl_course/feat3/utils.py

+from subprocess import check_output
+import nibabel as nb
+import numpy as np
+import hashlib
+import glob
+import os.path as op
+from feedsUtils import normalisedImageError
+
+def calcHash(value):
+    hashObj = hashlib.md5()
+    hashObj.update(value)
+    return hashObj.hexdigest()
+
+def hashCompare(file1, file2):
+    hash1 = calcHash(open(file1, 'rb').read())
+    hash2 = calcHash(open(file2, 'rb').read())
+    return (hash1 == hash2)
+
+def normImageError(im1, im2):
+    im1data = nb.load(im1).get_data()
+    im2data = nb.load(im2).get_data()
+    return np.abs((im1data-im2data)).mean() / np.abs(im1data).mean()
+
+def run(cmd):
+    print("RUNNING: " + cmd)
+    check_output(cmd, shell=True)
+
+# compare design files (excluding design.fsf) between original (benchmark) and new FEAT dirs
+def checkFeatDesigns(origFeatDir,newFeatDir,strict):
+    newDesigns = glob.glob(op.join(newFeatDir, "**", "design.*[!fsf]"),recursive=True)
+    for d in newDesigns:
+        print("TESTING: " + d)
+        isEqual = hashCompare(d, d.replace(newFeatDir, origFeatDir))
+        if not strict or isEqual:
+            continue
+        else:
+            print("FAIL: Design test failed: {0}".format(d))
+            sys.exit(1)
+
+# compare image files between original (benchmark) and new FEAT dirs
+def checkFeatImages(origFeatDir,newFeatDir,strict):
+    newImages = glob.glob(op.join(newFeatDir,"**", "*.nii.gz"), recursive=True)
+    for im in newImages:
+        print("TESTING: " + im)
+        imgError = normalisedImageError(im, im.replace(newFeatDir, origFeatDir))
+        if (imgError < 0.001):
+            continue
+        elif not strict and (imgError < 0.01):
+            continue
+        else:
+            print("FAIL: Image test failed: {0}".format(im))
+            sys.exit(1)

fsl_course/fslnets/feedsInputs

+fsl_course_data/rest/FSLNets
+fsl_course_data/rest/ICA/Group/groupmelodic60_fix.ica
+fsl_course_data/rest/ICA/Group/Designs
+fsl_course_data/rest/octave
\ No newline at end of file

fsl_course/fslnets/feedsRun

+#!/bin/bash -e
+#
+# This test calls Octave with the script test_run.m
+# The tests are run in that script directly.
+
+outdir=$1
+datadir=$2
+subdatadir=$datadir/fsl_course_data/rest
+
+# Environment variables required by Matlab/Octave
+export FEEDS_FSLNETS_OUTPUT=$outdir
+export FEEDS_FSLNETS_INPUT=$subdatadir/ICA/Group
+export FEEDS_FSLNETS_ROOT=$subdatadir/FSLNets
+export FEEDS_FSLNETS_OCTAVE_LIBS=$subdatadir/octave
+
+# Run the test
+octave < ./test_run.m
+status=$?
+
+if [ "$status" -eq "1" ]; then
+    exit 1
+fi
+
+for file in $outdir/*.matlab; do
+    echo $file
+    tail -n +2 $file > ${file/.matlab/.mat}
+done
+
+exit 0
\ No newline at end of file

fsl_course/fslnets/test_compare.m

+function success = test_compare()
+
+    % Environment variables
+    path_fslnets = getenv('FEEDS_FSLNETS_ROOT');
+    path_out     = getenv('FEEDS_FSLNETS_OUTPUT');
+    path_ref     = fullfile(path_fslnets, 'test_data_ref' );
+    
+    assert( ~isempty(path_out), 'Environment variables are not set properly.' )
+
+    % Define comparisons
+    compare(1).file   = 'netmats';
+    compare(1).names  = {'netmats1','netmats2','netmats3'};
+    compare(1).thresh = 1e-12;
+
+    compare(2).file   = 'group';
+    compare(2).names  = {'Znet1','Znet2','Znet3'};
+    compare(2).thresh = 1e-12;
+
+    compare(3).file   = 'group';
+    compare(3).names  = {'Mnet1','Mnet2','Mnet3'};
+    compare(3).thresh = 1e-12;
+
+    compare(4).file   = 'glm';
+    compare(4).names  = {'p_uncorrected3','p_corrected3'};
+    compare(4).thresh = 1e-12;
+
+    compare(5).file   = 'lda';
+    compare(5).names  = {'accuracy3_FLD','accuracy3_SVM'};
+    compare(5).thresh = 1e-12;
+
+    % Function used to compare arrays 
+    difference = @(x,y) max(abs(x(:) - y(:)));
+
+    % Iterate on comparisons
+    success = true;
+    ncomp   = numel(compare);
+    for i = 1:ncomp
+
+        % extract data
+        comp   = compare(i);
+        file   = comp.file;
+        names  = comp.names;
+        thresh = comp.thresh;
+
+        fprintf( 'Comparison #%d (file %s.matlab)...\n', i, file );
+ 
+        % load appropriate files
+        reference = load(fullfile( path_ref, [file '.matlab'] ));
+        target    = load(fullfile( path_out, [file '.matlab'] ));
+
+        % compare each variables
+        nvar = numel(names);
+        err  = zeros(1,nvar);
+
+        for j = 1:nvar            
+            err(j) = difference( reference.(names{j}), target.(names{j}) );
+            fprintf( '\t- Diff(%s) : %g\n', names{j}, err(j) );
+        end
+
+        % show status for this comparison
+        status = {'FAIL','PASS'};
+        nfail  = sum( err > thresh );
+        status = status{ 1+(nfail==0) };
+        fprintf( '\t+ Status: %s (%d/%d above %g)\n', status, nfail, nvar, thresh );
+
+        % update overall status
+        success = success && (nfail == 0); 
+
+    end
+
+end
+

fsl_course/fslnets/test_practical.m

+function test_practical( graphics )
+
+    if nargin < 1, graphics = 0; end
+
+    % Environment variables
+    env_input   = getenv('FEEDS_FSLNETS_INPUT');  % path to ICA/Group folder
+    env_output  = getenv('FEEDS_FSLNETS_OUTPUT'); % path to folder in which results should be saved
+
+    assert(~isempty(env_input) && ~isempty(env_output), 'Environment variables are not set properly.' );
+
+    % Setup paths
+    path_here   = fileparts(mfilename('fullpath'));
+    path_input  = env_input;
+    path_output = env_output;
+
+    %-----------------------------------
+
+
+    % Load data from dual regression
+    group_maps = fullfile(path_input,'groupmelodic60_fix.ica/melodic_IC');
+    data_dir   = fullfile(path_input,'groupmelodic60_fix.ica/Age.DR');
+
+    ts = nets_load( data_dir, 2, 1 );
+    
+    % Plot and check the spectra
+    if graphics
+        ts_spectra = nets_spectra(ts);
+    end
+
+    % Select "good" time-courses
+    tsid_remove = [12 14 15 23 25 26 27 29 31 32 34 35 37 38 39 40 41 42 43 44 45 46 47 49 51 52 53 54 55 56 58 59 60];
+    tsid_keep   = setdiff( 1:60, tsid_remove );
+
+    ts.DD = tsid_keep;
+    ts = nets_tsclean(ts,1); 
+
+
+    % Compute netmats
+    do_rtoz  = 0; % disable for testing
+    netmats1 = nets_netmats(ts,do_rtoz,'corr');
+    netmats2 = nets_netmats(ts,do_rtoz,'icov');
+    netmats3 = nets_netmats(ts,do_rtoz,'ridgep',1);   
+    
+    save(fullfile(path_output,'netmats.matlab'),'netmats1','netmats2','netmats3');
+
+    %-----------------------------------
+
+
+    % Group average
+    [Znet1,Mnet1] = nets_groupmean(netmats1,0);
+    [Znet2,Mnet2] = nets_groupmean(netmats2,0);
+    [Znet3,Mnet3] = nets_groupmean(netmats3,0);
+
+    save(fullfile(path_output,'group.matlab'),'Znet1','Znet2','Znet3','Mnet1','Mnet2','Mnet3');
+    # nets_netweb( Znet1, Znet2, ts.DD, group_maps, 'netweb' );
+
+    % Visualise the hierarchy
+    if graphics    
+        
+        % Using ward linkage
+        nets_hierarchy(Znet1,Znet2,ts.DD,group_maps);
+        
+        % Or with netjs
+        system('firefox netjs/index.html &');
+    end
+
+    %-----------------------------------
+
+
+    % Univariate cross-subject modelling
+    design   = fullfile( path_input, 'Designs/Age/design.mat' );
+    contrast = fullfile( path_input, 'Designs/Age/design.con' );
+    [p_uncorrected3,p_corrected3] = nets_glm(netmats3,design,contrast,graphics);
+    
+    save(fullfile(path_output,'glm.matlab'),'p_uncorrected3','p_corrected3');
+
+    % Investigate the significant differences
+    if graphics
+        
+        % Most significant edges
+        nets_edgepics(ts,group_maps,Znet3,reshape(p_corrected3(4,:),ts.Nnodes,ts.Nnodes),6);
+    
+        % Visualise details of the t-test for particular edge
+        nets_boxplots(ts,netmats3,23,7,10);
+    end
+
+    %-----------------------------------
+
+
+    % Run LDA
+    accuracy3_FLD = nets_lda(netmats3,10,2);
+    accuracy3_SVM = nets_lda(netmats3,10,8);
+    
+    save(fullfile(path_output,'lda.matlab'),'accuracy3_FLD','accuracy3_SVM');
+
+end

fsl_course/fslnets/test_run.m

+env_fsl     = getenv('FSLDIR');
+env_fslnets = getenv('FEEDS_FSLNETS_ROOT');
+env_libs = getenv('FEEDS_FSLNETS_OCTAVE_LIBS')
+
+path_fsl     = fullfile( env_fsl,     'etc/matlab' );
+path_fslnets = fullfile( env_fslnets);
+path_stats   = fullfile( env_libs, 'statistics-1.2.4/inst' );
+path_svm     = fullfile( env_libs, 'libsvm/matlab' );
+
+addpath(path_fsl);
+addpath(path_fslnets);
+addpath(path_stats);
+addpath(path_svm);
+
+path
+
+test_practical(0);
+

fsl_course/intro/feedsInputs

+fsl_course_data/intro

fsl_course/intro/feedsRun

+#!/bin/bash -e
+
+./test_fslutils1.sh $@
+./test_fslutils2.sh $@
+./test_fslutils3.sh $@
+./test_fslutils4.sh $@
+./test_bet.sh       $@

fsl_course/intro/test_bet.sh

+#!/bin/bash
+# Test from the FSL course intro practical.
+# Tests a few instantiations of BET, and
+# robustfov.
+
+set -e
+
+outdir=$1
+datadir=$2
+subdatadir=$datadir/fsl_course_data/intro
+
+bet $subdatadir/structural $outdir/structural_brain       -f 0.5 -g 0 -m -s
+bet $subdatadir/structural $outdir/structural_brain_lowf  -f 0.2
+bet $subdatadir/structural $outdir/structural_brain_highf -f 0.8
+
+
+bet $subdatadir/sub3m0 $outdir/sub3m0_brain
+bet $subdatadir/sub3m0 $outdir/sub3m0_brain_g0.5      -g 0.5
+bet $subdatadir/sub3m0 $outdir/sub3m0_brain_g0.2_f0.3 -g 0.2 -f 0.3
+
+bet $subdatadir/bighead   $outdir/bighead_brain
+bet $subdatadir/bighead   $outdir/bighead_brain_R          -R
+bet $subdatadir/bighead   $outdir/bighead_brain_c89_97_161 -c 89 97 161
+
+robustfov -i $subdatadir/bighead -r $outdir/bighead_crop > $outdir/robustfov_bighead.output
+fslhd $outdir/bighead_crop | tail -n+2 > $outdir/fslhd_bighead_crop.output
+bet $outdir/bighead_crop  $outdir/bighead_crop_brain
+
+exit 0

fsl_course/intro/test_fslutils1.sh

+#!/bin/bash
+# Tests from the FSL course intro practical.
+# Does a few basic tests on these commands:
+# 
+#   - fslinfo
+#   - fslhd
+#   - fslstats
+
+set -e
+
+outdir=$1
+datadir=$2
+subdatadir=$datadir/fsl_course_data/intro
+
+# First, we're going to run a bunch 
+# of basic commands on three images.
+testImages="highres thresh_zstat1 egepi"
+
+# We run fslhd through tail,
+# because it outputs the file
+# name (as specified in the
+# call) on the first line.
+for image in $testImages; do
+  fslinfo  $subdatadir/$image              > $outdir/fslinfo_${image}.output
+  fslhd    $subdatadir/$image  | tail -n+2 > $outdir/fslhd_${image}.output
+  fslstats $subdatadir/$image -R -r        > $outdir/fslstats_Rr_${image}.result
+  fslstats $subdatadir/$image -l 1 -R -r   > $outdir/fslstats_l1Rr_${image}.result
+  fslstats $subdatadir/$image -m -M        > $outdir/fslstats_mM_${image}.result
+done
+
+exit 0

fsl_course/intro/test_fslutils2.sh

+#!/bin/bash
+# 
+# Tests from the FSL course intro practical.
+# Does a few basic fslmaths and fslstats tests.
+
+set -e
+
+outdir=$1
+datadir=$2
+subdatadir=$datadir/fsl_course_data/intro
+
+# Run some fslmaths commands
+fslmaths $subdatadir/image0 -sub $subdatadir/image1 $outdir/imdiff
+fslmaths $outdir/imdiff -div $subdatadir/image0 -mul 100 $outdir/imdiffpercent
+fslmaths $subdatadir/LThal_mask_func -thr 0.5 -bin $outdir/LThal_mask_func_bin 
+
+# Test fslstats on the above images.
+# Using md5sum for an exact match
+fslstats $outdir/LThal_mask_func_bin -V                              > $outdir/fslstats_LThal_mask_1.result
+fslstats $subdatadir/thresh_zstat1 -k $outdir/LThal_mask_func_bin -V > $outdir/fslstats_LThal_mask_2.result
+fslstats $subdatadir/thresh_zstat1 -k $outdir/LThal_mask_func_bin -M > $outdir/fslstats_LThal_mask_3.result
+
+exit 0

fsl_course/intro/test_fslutils3.sh

+#!/bin/bash
+# Tests from the FSL course intro practical.
+# Tests the fslsplit and fslmerge cmmands.
+
+set -e
+
+outdir=$1
+datadir=$2
+subdatadir=$datadir/fsl_course_data/intro
+
+# Split a diffusion image (assumed to be 64 volumes)
+fslsplit $subdatadir/diffdata $outdir/vol
+
+# Rename the nodif volumes
+mv $outdir/vol0000.nii.gz $outdir/b0vol0000.nii.gz
+mv $outdir/vol0016.nii.gz $outdir/b0vol0016.nii.gz
+mv $outdir/vol0032.nii.gz $outdir/b0vol0032.nii.gz
+mv $outdir/vol0048.nii.gz $outdir/b0vol0048.nii.gz
+
+# Merge the diffusion and nodif volumes.
+# Be specific on the globbing, because
+# testimage will have created vol*_error
+# and b0vol*_error images.
+fslmerge -t $outdir/dwidata    $outdir/vol????.*
+fslmerge -t $outdir/nodiffdata $outdir/b0vol????.*
+
+exit 0

fsl_course/intro/test_fslutils4.sh

+#!/bin/bash
+# Tests from the FSL course intro practical.
+# Tests the fslroi command.
+
+set -e
+
+outdir=$1
+datadir=$2
+subdatadir=$datadir/fsl_course_data/intro
+
+# Do some fslroi-ing
+fslroi $subdatadir/wrapped $outdir/back  0 256 70 186 0 128
+fslroi $subdatadir/wrapped $outdir/front 0 256 0  70  0 128
+fslmerge -y $outdir/fix_wrapped $outdir/back.nii.gz $outdir/front.nii.gz
+fslroi $subdatadir/highres $outdir/highres_roi `fslstats $subdatadir/highres -w`
+
+# Save the output of fslhd too, to
+# make sure fslroi sets affines
+# correctly.
+# Crop the first line of fslhd,
+# because it prints the file name,
+# which may not match.
+fslhd $outdir/back        | tail -n+2 > $outdir/fslhd_back.output
+fslhd $outdir/front       | tail -n+2 > $outdir/fslhd_front.output
+fslhd $outdir/fix_wrapped | tail -n+2 > $outdir/fslhd_fix_wrapped.output
+fslhd $outdir/highres_roi | tail -n+2 > $outdir/fslhd_highres_roi.output
+
+exit 0
+

fsl_course/melodic/feedsInputs

+fsl_course_data/rest/ICA/Group

fsl_course/melodic/feedsRun

+#!/bin/bash -e
+
+./testGroupMelodic.sh $@

fsl_course/melodic/testGroupMelodic.sh

+#!/bin/bash -e
+#
+# A FEEDS test for melodic
+
+outdir=$1
+datadir=$2
+subdatadir=$datadir/fsl_course_data/rest/ICA/Group
+
+unset SGE_ROOT
+
+#Generate full paths to clean functionals, in standard space, for melodic
+awk -v DATADIR="$subdatadir/" '{print DATADIR$0}' ${subdatadir}/inputlist_clean.txt > $outdir/inputlist_clean.txt
+
+#Run melodic
+melodic -i $outdir/inputlist_clean.txt -o $outdir/groupmelodic_fix.ica -d 25 -a concat --tr=2.000 --nobet --bgthreshold=10 -v --report --mmthresh=0.5 --Ostats --seed=1
+
+echo "Finished group Melodic."
+
+exit 0