Add FSL dual regression as an option

Now computes the proper `melodic` & `dual_regression` pipeline, which is a good
sanity check.

IO/loadDualRegression.m

+function [ P, A ] = loadDualRegression( dualregDir, params )
+
+% Subject maps / timecourses
+P = cell(params.S, 1);
+A = cell(params.S, 1);
+sr = 0;
+for s = 1:params.S
+    P{s} = 0.0;
+    A{s} = cell(params.R(s), 1);
+    for r = 1:params.R(s)
+        subj = sprintf('subject%05d', sr);
+        
+        % Collate run-specific maps
+        Psr = read_avw(fullfile( ...
+            dualregDir, ['dr_stage2_' subj '.nii.gz']));
+        P{s} = P{s} + reshape(Psr, params.V, params.iN);
+        
+        % Timecourses
+        A{s}{r} = load(fullfile( ...
+            dualregDir, ['dr_stage1_' subj '.txt']))';
+        
+        sr = sr + 1;
+    end
+
+    % Take mean map
+    P{s} = P{s} ./ params.R(s);
+end
+
+end

IO/saveNIfTIs.m

@@ -31,5 +31,8 @@ fclose(fileID);
 fileID = fopen(fullfile(directory, 'MELODIC_SpecFile.txt'), 'w');
 fprintf(fileID, '%s\n', fileNames{:});
 fclose(fileID);
+fileID = fopen(fullfile(directory, 'DualRegression_SpecFile.txt'), 'w');
+fprintf(fileID, '%s ', fileNames{:});
+fclose(fileID);
 
 end

Methods/DualRegression.sh

+#!/bin/sh
+
+# Runs dual regression on simulated dataset
+# Usage: DualRegression.sh <output_dir> <nifti_dir> <group_maps>
+
+output_dir=$1
+nifti_dir=$2
+group_maps=$3
+
+mkdir "${output_dir}"
+
+# Run all FSL commands locally
+unset SGE_ROOT
+
+# Run dual_regression
+des_norm=1
+design=-1
+n_perm=0
+dual_regression "${group_maps}" \
+    ${des_norm} ${design} ${n_perm} \
+    "${output_dir}" \
+    $(cat "${nifti_dir}/DualRegression_SpecFile.txt") \
+    > "${output_dir}/TerminalOutput.txt"

RunTests.m

@@ -359,16 +359,21 @@ for n = 1:params.nRepeats
         
     end
     
-    %% Run MELODIC %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    %% Run MELODIC / DR %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
     try
         melodicDir = fullfile(repeatDir, 'MELODIC.gica');
         system(sprintf( ...
             'sh Methods/MELODIC.sh %s %s %d %1.2f', ...
             melodicDir, niftiDir, params.iN, params.TR));
+        dualregDir = fullfile(repeatDir, 'MELODIC.dr');
+        system(sprintf( ...
+            'sh Methods/DualRegression.sh %s %s %s', ...
+            dualregDir, niftiDir, fullfile(melodicDir, 'melodic_IC.nii.gz')));
         
         icadr_Pg = loadMELODIC(melodicDir, params);
-        [icadr_P, icadr_A] = runDR(D, icadr_Pg, params);
+        [icadr_P, icadr_A] = loadDualRegression(dualregDir, params);
+        %[icadr_P, icadr_A] = runDR(D, icadr_Pg, params);
         icadr_pcA = calculateNetmats(icadr_A, params);
         
         scores.ICA_DR(n) = calculateDecompositionAccuracy( ...