From dc9a8906d4e08c7603279f336153dcd0f732834f Mon Sep 17 00:00:00 2001 From: nicholst Date: Thu, 18 Jun 2020 09:53:05 +0100 Subject: [PATCH 1/2] Clean up white spaces, trailing newlines --- .travis.yml | 28 +- @swe_DataType/display.m | 6 +- @swe_DataType/eq.m | 4 +- @swe_DataType/get.m | 10 +- @swe_DataType/swe_DataType.m | 8 +- @swe_cifti/Contents.m | 32 +- @swe_cifti/create.m | 16 +- @swe_cifti/subsasgn.m | 2 +- @swe_cifti/subsref.m | 6 +- @swe_cifti/swe_cifti.m | 8 +- swe.m | 124 +++---- swe_Box_test.m | 2 +- swe_DesRep.m | 34 +- swe_VOI.m | 8 +- swe_batch.m | 2 +- swe_boxCoxTransform.m | 4 +- swe_cfg_batch.m | 4 +- swe_checkCompat.m | 24 +- swe_cifti_clusters.m | 2 +- swe_cifti_convertVolLabels.m | 4 +- swe_cifti_max.m | 10 +- swe_compareVersions.m | 16 +- swe_config_results.m | 2 +- swe_config_smodel.m | 14 +- swe_conman.m | 42 +-- swe_contrasts.m | 116 +++---- swe_contrasts_WB.m | 38 +- swe_cp.m | 176 +++++----- swe_cp_WB.m | 650 +++++++++++++++++------------------ swe_create_vol.m | 6 +- swe_data_hdr_read.m | 12 +- swe_data_hdr_write.m | 16 +- swe_data_read.m | 2 +- swe_data_write.m | 6 +- swe_defaults.m | 8 +- swe_duplication_matrix.m | 2 +- swe_estimateBoxCoxLambda.m | 10 +- swe_getClusterStatistics.m | 12 +- swe_getSPM.m | 252 +++++++------- swe_get_defaults.m | 4 +- swe_get_file_extension.m | 8 +- swe_invNcdf.m | 6 +- swe_jobman.m | 77 ++--- swe_list.m | 216 ++++++------ swe_max.m | 2 +- swe_mesh_clusters.m | 6 +- swe_mesh_max.m | 4 +- swe_progress_bar.m | 8 +- swe_read_cifti_info.m | 10 +- swe_regions.m | 58 ++-- swe_results.m | 10 +- swe_results_ui.m | 274 +++++++-------- swe_rmodel.m | 10 +- swe_run_results.m | 4 +- swe_run_rmodel.m | 2 +- swe_run_smodel.m | 62 ++-- swe_seed.m | 4 +- swe_select.m | 37 +- swe_smodel.m | 10 +- swe_splitCovariate.m | 10 +- swe_tfce_transform.m | 6 +- swe_uc_FDR.m | 25 +- swe_ui_main.m | 24 +- swe_update.m | 14 +- swe_vechCovVechV.m | 28 +- test/runTest.m | 84 ++--- 66 files changed, 1359 insertions(+), 1362 deletions(-) diff --git a/.travis.yml b/.travis.yml index d34542b..84416a5 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,9 +1,9 @@ language: python # use container-based infrastructure sudo: required -services: +services: - docker -cache: +cache: timeout: 1000 directories: - test/data @@ -11,16 +11,16 @@ python: - "3.5" bundler_args: --retry 9 # command to install dependencies, e.g. pip install -r requirements.txt --use-mirrors -install: - # Install git-lfs +install: + # Install git-lfs - cd .. - curl -sLo - https://github.com/github/git-lfs/releases/download/v1.4.1/git-lfs-linux-amd64-1.4.1.tar.gz | tar xzvf - - export PATH=$PATH:`pwd`/git-lfs-1.4.1 - - git lfs install + - git lfs install # Download data. - cd SwE-toolbox/test/data - if ! [ -d .git ]; then git clone https://github.com/NISOx-BDI/SwE-toolbox-testdata.git .; fi - # Delete any previous changes (retry because lfs might download files) + # Delete any previous changes (retry because lfs might download files) - git checkout master - travis_retry git reset --hard origin/master # A second time to allow for 3 more retries as "--retry 9" does not seem to be @@ -31,7 +31,7 @@ jobs: include: # Each '-' with the same indentation as this comment is a job that is run in parallel in the # results generation stage. - - stage: Results generation + - stage: Results generation env: - testname='p_vox_t_nii' script: @@ -41,7 +41,7 @@ jobs: - testname='p_vox_f_nii' script: - did=$(sudo docker run -ti -d --rm -v `pwd`:/swe cmaumet/octave-spm) - - sudo docker exec -t -i $did octave --no-window-system --eval "addpath('/swe/test');runTest('p', 'vox', 'f', 'nii')" + - sudo docker exec -t -i $did octave --no-window-system --eval "addpath('/swe/test');runTest('p', 'vox', 'f', 'nii')" - env: - testname='wb_vox_t_nii' script: @@ -66,7 +66,7 @@ jobs: - testname='wb_tfce_t_nii' script: - did=$(sudo docker run -ti -d --rm -v `pwd`:/swe cmaumet/octave-spm) - - sudo docker exec -t -i $did octave --no-window-system --eval "addpath('/swe/test');runTest('wb', 'tfce', 't', 'nii')" + - sudo docker exec -t -i $did octave --no-window-system --eval "addpath('/swe/test');runTest('wb', 'tfce', 't', 'nii')" - env: - testname='wb_tfce_f_nii' script: @@ -81,7 +81,7 @@ jobs: - testname='p_dat_f_mat' script: - did=$(sudo docker run -ti -d --rm -v `pwd`:/swe cmaumet/octave-spm) - - sudo docker exec -t -i $did octave --no-window-system --eval "addpath('/swe/test');runTest('p', 'dat', 'f', 'mat')" + - sudo docker exec -t -i $did octave --no-window-system --eval "addpath('/swe/test');runTest('p', 'dat', 'f', 'mat')" - env: - testname='wb_dat_t_mat' script: @@ -111,7 +111,7 @@ jobs: - testname='p_dpx_f_gii' script: - did=$(sudo docker run -ti -d --rm -v `pwd`:/swe cmaumet/octave-spm) - - sudo docker exec -t -i $did octave --no-window-system --eval "addpath('/swe/test');runTest('p', 'dpx', 'f', 'gii')" + - sudo docker exec -t -i $did octave --no-window-system --eval "addpath('/swe/test');runTest('p', 'dpx', 'f', 'gii')" - env: - testname='wb_dpx_t_gii' script: @@ -141,7 +141,7 @@ jobs: - testname='p_dpx_f_cii' script: - did=$(sudo docker run -ti -d --rm -v `pwd`:/swe cmaumet/octave-spm) - - sudo docker exec -t -i $did octave --no-window-system --eval "addpath('/swe/test');runTest('p', 'dpx', 'f', 'cii')" + - sudo docker exec -t -i $did octave --no-window-system --eval "addpath('/swe/test');runTest('p', 'dpx', 'f', 'cii')" - env: - testname='wb_dpx_t_cii' script: @@ -166,6 +166,6 @@ before_install: # Update docker version - sudo apt-get update - sudo apt-get -y -o Dpkg::Options::="--force-confnew" install docker-ce - - git config --global user.name "TravisCI" + - git config --global user.name "TravisCI" - git config --global user.email "travis@dummy.com" - - docker pull cmaumet/octave-spm + - docker pull cmaumet/octave-spm diff --git a/@swe_DataType/display.m b/@swe_DataType/display.m index 51feab5..b28c2db 100644 --- a/@swe_DataType/display.m +++ b/@swe_DataType/display.m @@ -3,7 +3,7 @@ function display(obj) % ========================================================================= % Bryan Guillaume % Version Info: $Format:%ci$ $Format:%h$ - + display(get(obj)); - -end \ No newline at end of file + +end diff --git a/@swe_DataType/eq.m b/@swe_DataType/eq.m index f701b21..8be828f 100644 --- a/@swe_DataType/eq.m +++ b/@swe_DataType/eq.m @@ -3,7 +3,7 @@ % ========================================================================= % Bryan Guillaume % Version Info: $Format:%ci$ $Format:%h$ - + out = strcmp(a.value, b.value); -end \ No newline at end of file +end diff --git a/@swe_DataType/get.m b/@swe_DataType/get.m index fa91057..0620201 100644 --- a/@swe_DataType/get.m +++ b/@swe_DataType/get.m @@ -1,9 +1,9 @@ function out = get(obj) - % Get a swe_DataType object - % ========================================================================= + % Get a swe_DataType object + % ========================================================================= % Bryan Guillaume % Version Info: $Format:%ci$ $Format:%h$ - - out = obj.value; -end \ No newline at end of file + out = obj.value; + +end diff --git a/@swe_DataType/swe_DataType.m b/@swe_DataType/swe_DataType.m index 3fec294..2c00c5b 100644 --- a/@swe_DataType/swe_DataType.m +++ b/@swe_DataType/swe_DataType.m @@ -5,9 +5,9 @@ % ========================================================================= % Bryan Guillaume % Version Info: $Format:%ci$ $Format:%h$ - + availableDataTypes = {'Nifti', 'Gifti', 'Cifti', 'Mat', 'VolumeMat', 'SurfaceMat'}; - + if nargin == 0 error( "swe_DataType cannot be initialized empty; choose from 'Nifti', 'Gifti', 'Cifti', 'Mat', 'VolumeMat' or 'SurfaceMat'."); end @@ -17,7 +17,7 @@ end obj.value = dataType; - + obj = class(obj, 'swe_DataType'); -end \ No newline at end of file +end diff --git a/@swe_cifti/Contents.m b/@swe_cifti/Contents.m index c665fe0..323cc6c 100644 --- a/@swe_cifti/Contents.m +++ b/@swe_cifti/Contents.m @@ -18,17 +18,17 @@ % dat.dim = [64 64 32]; % dat.dtype = 'FLOAT64-BE'; % dat.offset = ceil(348/8)*8; -% +% % % alternatively: % % dat = file_array( 'junk.nii',dim,dtype,off,scale,inter) -% +% % disp(dat) -% +% % % Create an empty NIFTI structure % N = nifti; -% +% % fieldnames(N) % Dump fieldnames -% +% % % Creating all the NIFTI header stuff % N.dat = dat; % N.mat = [2 0 0 -110 ; 0 2 0 -110; 0 0 -2 92; 0 0 0 1]; @@ -36,41 +36,41 @@ % N.mat_intent = 'Scanner'; % N.mat0 = N.mat; % N.mat0_intent = 'Aligned'; -% +% % N.diminfo.slice = 3; % N.diminfo.phase = 2; % N.diminfo.frequency = 2; -% N.diminfo.slice_time.code='xxx'; % dump possibilities +% N.diminfo.slice_time.code='xxx'; % dump possibilities % N.diminfo.slice_time.code = 'sequential_increasing'; % N.diminfo.slice_time.start = 1; % N.diminfo.slice_time.end = 32; % N.diminfo.slice_time.duration = 3/32; -% +% % N.intent.code='xxx' ; % dump possibilities % N.intent.code='FTEST'; % or N.intent.code=4; % N.intent.param = [4 8]; -% +% % N.timing.toffset = 28800; % N.timing.tspace=3; % N.descrip = 'This is a NIFTI-1 file'; % N.aux_file='aux-file-name.txt'; % N.cal = [0 1]; -% +% % create(N); % Writes hdr info -% +% % dat(:,:,:)=0; % Write out the data as all zeros -% +% % [i,j,k] = ndgrid(1:64,1:64,1:32); % dat(find((i-32).^2+(j-32).^2+(k*2-32).^2 < 30^2))=1; % Write some ones in the file % dat(find((i-32).^2+(j-32).^2+(k*2-32).^2 < 15^2))=2; -% -% +% +% % % displaying a slice % imagesc(dat(:,:,12));colorbar -% +% % % get a handle to 'junk.nii'; % M=nifti('junk.nii'); -% +% % imagesc(M.dat(:,:,12)); % % _______________________________________________________________________ diff --git a/@swe_cifti/create.m b/@swe_cifti/create.m index f9c87bd..e91e49e 100644 --- a/@swe_cifti/create.m +++ b/@swe_cifti/create.m @@ -13,21 +13,21 @@ function create(obj,wrt) for i=1:numel(obj) - + o = obj(i); - + if ~isa(o.dat,'file_array'), error('Data must be a file_array.'); end - + fname = o.dat.fname; if isempty(fname), error('No filename to write to.'); end - + %-Write NIFTI header sts = write_hdr_raw(fname, o.hdr, o.dat.dtype(end-1)=='B'); if ~sts, error('Unable to write header for "%s".',fname); end - + %-Write extra information write_extras(fname,o.extras); - + %-Create an empty image file if necessary if nargin>1 && any(wrt==1) [pth,nam] = fileparts(fname); @@ -37,8 +37,8 @@ function create(obj,wrt) ext = '.img'; end o.dat.fname = fullfile(pth,[nam ext]); - + initialise(o.dat); end - + end diff --git a/@swe_cifti/subsasgn.m b/@swe_cifti/subsasgn.m index 43fb29d..fd514dc 100644 --- a/@swe_cifti/subsasgn.m +++ b/@swe_cifti/subsasgn.m @@ -27,7 +27,7 @@ if length(subs)>1 t = subsref(objs,subs(1)); % A lot of this stuff is a little flakey, and may cause Matlab to bomb. - % + % %if numel(t) ~= nargin-2, % error('The number of outputs should match the number of inputs.'); %end; diff --git a/@swe_cifti/subsref.m b/@swe_cifti/subsref.m index e4fee69..9b0d658 100644 --- a/@swe_cifti/subsref.m +++ b/@swe_cifti/subsref.m @@ -86,7 +86,7 @@ if s d = findindict(s,'units'); if ~isempty(d) - t = diag([d.rescale*[1 1 1] 1])*t; + t = diag([d.rescale*[1 1 1] 1])*t; end end @@ -151,12 +151,12 @@ se = double(h.slice_end)+1; ss = max(ss,1); se = min(se,double(h.dim(t.slice+1))); - + sd = double(h.slice_duration); s = double(bitand(h.xyzt_units,24)); d = findindict(s,'units'); if d.rescale, sd = sd*d.rescale; end - + ns = (se-ss+1); d = findindict(sc,'sliceorder'); if isempty(d) diff --git a/@swe_cifti/swe_cifti.m b/@swe_cifti/swe_cifti.m index 91b1c8d..1f09865 100644 --- a/@swe_cifti/swe_cifti.m +++ b/@swe_cifti/swe_cifti.m @@ -20,7 +20,7 @@ case {1, 2} if isa(varargin{1},'swe_cifti') h = varargin{1}; - + elseif ischar(varargin{1}) if size(varargin{1},1)>1 h = swe_cifti(cellstr(varargin{1})); @@ -75,7 +75,7 @@ h = class(h,'swe_cifti'); elseif isstruct(varargin{1}) - + h = struct('hdr', varargin{1}.hdr,... 'dat', varargin{1}.dat,... 'extras', varargin{1}.extras); @@ -92,8 +92,8 @@ else error('Don''t know what to do yet.'); end - + otherwise error('Don''t know what to do yet.'); end -end \ No newline at end of file +end diff --git a/swe.m b/swe.m index b128ed1..797a81e 100644 --- a/swe.m +++ b/swe.m @@ -27,8 +27,8 @@ switch lower(Action) %================================================================== - case 'startup' % Startup the toolbox - %================================================================== + case 'startup' % Startup the toolbox + %================================================================== % check installation of toolbox and of SPM8/SPM12 ok = check_installation; if ~ok @@ -36,26 +36,26 @@ fprintf('INSTALLATION PROBLEM!'); return end - + % Welcome message swe('ASCIIwelcome'); - + % Add pathes for SPM functions addpath(fullfile(spm('Dir'),'matlabbatch')) if strncmpi(spm('ver'),'spm12',5) addpath(fullfile(spm('Dir'),'compat')) end - + % Add path to SwE toolbox. addpath(fileparts(mfilename('fullpath'))); - + % launch the main GUI swe_ui_main; - + % print present working directory fprintf('SwE present working directory:\n\t%s\n',pwd) - - + + %================================================================== case 'asciiwelcome' %-ASCII swe banner welcome %================================================================== @@ -64,20 +64,20 @@ a = generateAscii(['SwE v' tmp]); fprintf('%s \n', a{1}, a{2}, a{3}, a{4}); fprintf('swe v%s \n', versionNo); - + case 'colour' - + varargout = {[0.8 0.8 1.0], 'Diluted Blackcurrent Purple'}; - + case 'ver' - + varargout{1}=versionNo; - + %================================================================== otherwise %-Unknown action string %================================================================== error('Unknown action string'); - + end return @@ -116,125 +116,125 @@ end % The following functions are for generating the ascii welcome message for -% SwE. +% SwE. % ------------------------------------------------------------------------- % This method converts character 'char' to it's equivalent 4-line ascii % art. New characters can be added by creating new cases in the below % switch. function aChar = char2ascii(char) - + switch char - + case '0' - + aChar = {' ___ ',... '| |',... '| | |',... '|___|'}; - + case '1' - + aChar = {' _ ',... '/_|',... ' ||',... ' ||'}; - + case '2' - + aChar = {' ___ ',... '(__ \',... '/ __/',... '\___)'}; - + case '3' - + aChar = {' ___ ',... '(__ )',... ' (_ \',... - '(___/'}; + '(___/'}; case '4' - + aChar = {' __ ',... ' /. |',... '(_ _)',... - ' (_) '}; - case '5' - + ' (_) '}; + case '5' + aChar = {' ___ ',... '| __)',... '|__ \',... '(___/'}; - + case '6' - + aChar = {' _ ',... ' / ) ',... '/ _ \',... - '\___/'}; + '\___/'}; case '7' - + aChar = {' ___ ',... '(__ )',... ' / / ',... - '(_/ '}; + '(_/ '}; case '8' - + aChar = {' ___ ',... '| _ |',... '| _ |',... - '|___| '}; + '|___| '}; case '9' - + aChar = {' ___ ',... '/ _ \',... '\_ /',... - ' (_/ '}; - + ' (_/ '}; + case '.' - + aChar = {' ',... ' ',... ' _ ',... - '|_|'}; - + '|_|'}; + case 'S' - + aChar = {' ___ ',... '/ __)',... '\__ \',... - '(___/'}; - + '(___/'}; + case 'w' - + aChar = {' ',... '_ _',... '\\/\//',... - ' \/\/ '}; - + ' \/\/ '}; + case 'E' - + aChar = {' ___ ',... '| __)',... '| __)',... - '|___)'}; - + '|___)'}; + case ' ' - + aChar = {' ',... ' ',... ' ',... ' '}; - + case 'v' - + aChar = {' ',... '_ _ ',... '\\// ',... - ' \/ '}; - - + ' \/ '}; + + end - + end % This function takes two cell arrays containing ascii art in the form @@ -247,11 +247,11 @@ % This function takes in a string as input and recursively generates the % ASCII art representation of said string. function ascii = generateAscii(str) - + if length(str)~=1 ascii = concatAscii(generateAscii(str(1:(end-1))), generateAscii(str(end))); else ascii = char2ascii(str); end - -end \ No newline at end of file + +end diff --git a/swe_Box_test.m b/swe_Box_test.m index f6d923c..f5442b4 100644 --- a/swe_Box_test.m +++ b/swe_Box_test.m @@ -114,4 +114,4 @@ tmp(Q) = p; swe_data_write(CS_test_p,tmp); -end \ No newline at end of file +end diff --git a/swe_DesRep.m b/swe_DesRep.m index ac52d30..fe6be77 100644 --- a/swe_DesRep.m +++ b/swe_DesRep.m @@ -64,7 +64,7 @@ % If (and only if) both vectors have zero mean, i.e. % sum(a)==sum(b)==0, then the cosine of the angle between the % vectors is the same as the correlation between the two variates. -% +% % The design orthogonality matrix is "surfable": Clicking (and % holding or dragging) the cursor around the design orthogonality % image reports the orthogonality of the correponding pair of @@ -83,7 +83,7 @@ % If not an fMRI design, then the Explore sub-menu has two options: % "Files and factors" & "Covariates". % -% * Explore: Files and factors - Multi-page listing of filenames, +% * Explore: Files and factors - Multi-page listing of filenames, % factor indicies and covariates. % % The covariates printed are the raw covariates as entered into @@ -91,8 +91,8 @@ % after any grand mean scaling. % % * Explore: Covariates - Plots of the covariates, showing how they are -% included into the model. -% +% included into the model. +% % Covariates are plotted, one per page, overlaid on the design % matrix. The description strings in the xC covariate structure % array are displayed. The corresponding design matrix column(s) @@ -149,7 +149,7 @@ % .swd - SPM working directory - directory where configuration file resides % [defaults to empty] % .SPMid - (recommended) ID string of creator program. -% h - handle of menu created ('Tag'ged as 'DesRepUI') +% h - handle of menu created ('Tag'ged as 'DesRepUI') % % % FORMAT spm_DesRep('Files&Factors',P,I,xC,sF,xs) @@ -256,13 +256,13 @@ % Version Info: $Format:%ci$ $Format:%h$ -SVNid = '$Rev: 6351 $'; +SVNid = '$Rev: 6351 $'; %-Format arguments %-------------------------------------------------------------------------- if ~nargin SPMid = spm('FnBanner',mfilename,SVNid); - hC = spm_DesRep('DesRepUI'); + hC = spm_DesRep('DesRepUI'); SPM = get(hC,'UserData'); if ~isempty(SPM) cb_menu([],[],'DesMtx',SPM); @@ -630,7 +630,7 @@ hPEstAx = []; end spm_figure('NewPage',[hCovMtx;get(hCovMtx,'Children');hPEstAx;get(hPEstAx,'Children');hCovMtxSc]) - + %-Show components of covariance matrix %-------------------------------------------------------------------------- if isfield(varargin{2},'Vi') @@ -760,7 +760,7 @@ %-Filenames % ( Show at most 32, showing every 2nd/3rd/4th/... as necessary to pair ) -% ( down to <32 items. Always show last item so #images is indicated. ) +% ( down to <32 items. Always show last item so #images is indicated. ) if desmtx && ~isempty(fnames) axes('Position',[.68 .4 .3 .4],'Visible','off',... 'DefaultTextFontSize',FS(8),... @@ -808,7 +808,7 @@ hDesO = axes('Position',[.07 .18 .6 .2]); tmp = 1-abs(O); tmp(logical(tril(ones(nPar),-1))) = 1; hDesOIm = image(tmp*64); - + set(hDesO,'Box','off','TickDir','out',... 'XaxisLocation','top','XTick',PTick,'XTickLabel','',... 'YaxisLocation','right','YTick',PTick,'YTickLabel','',... @@ -859,14 +859,14 @@ set(hAx,'Units','points'); AxPos = get(hAx,'Position'); set(hAx,'YLim',[0,AxPos(4)]) - + dy = FS(9); y0 = floor(AxPos(4)) -dy; y = y0; text(0.3,y,str,... 'HorizontalAlignment','Center',... 'FontWeight','Bold','FontSize',FS(11)) y=y-2*dy; - + for sf = fieldnames(xs)' text(0.3,y,[strrep(sf{1},'_',' '),' :'],... 'HorizontalAlignment','Right','FontWeight','Bold',... @@ -1056,11 +1056,11 @@ %-Descriptions %---------------------------------------------------------------------- hDAx = axes('Position',[0.03,0.1,0.94,0.30],'Visible','off'); - + set(hDAx,'Units','points'); tmp = get(hDAx,'Position'); set(hDAx,'YLim',[0,tmp(4)]) - + dy = FS(9); y0 = floor(tmp(4)) -dy; y = y0; %-Description strings from xC(i).descrip @@ -1145,7 +1145,7 @@ %========================================================================== % spm_DesRep('ScanTick',nScan,lim) % ( Show at most 32, showing every 2nd/3rd/4th/... as necessary to pair ) -% ( down to <32 items. Always show last item so #images is indicated. ) +% ( down to <32 items. Always show last item so #images is indicated. ) if nargin<3, lim=32; else lim=varargin{3}; end if nargin<2, error('insufficient arguments'), end nScan = varargin{2}; @@ -1567,7 +1567,7 @@ function cb_menu(obj,evt,action,SPM,varargin) spm_DesRep('DesMtx',SPM.xX,... filenames,... SPM.xsDes); - + case 'DesOrth' spm_DesRep('DesOrth',SPM.xX); @@ -1584,4 +1584,4 @@ function cb_menu(obj,evt,action,SPM,varargin) case 'xVi' spm_DesRep('xVi', SPM.xVi); -end \ No newline at end of file +end diff --git a/swe_VOI.m b/swe_VOI.m index c252f54..8d90d14 100644 --- a/swe_VOI.m +++ b/swe_VOI.m @@ -4,7 +4,7 @@ % FORMAT [TabDat,xSVC] = swe_VOI(SwE,xSwE,hReg,[xY]) % ------------------------------------------------------------------------- % Inputs: -% +% % SwE - Structure containing analysis details (see spm_spm) % % xSwE - Structure containing SwE, distribution & filtering details @@ -52,7 +52,7 @@ % voxels with values greater than 0. % % See also: spm_list -% Adapted version of `spm_VOI.m`. +% Adapted version of `spm_VOI.m`. % Author of Adaptation: Tom Maullin (07/09/2018) % Version Info: $Format:%ci$ $Format:%h$ %__________________________________________________________________________ @@ -90,7 +90,7 @@ catch xyzmm = spm_results_ui('GetCoords'); end - + %-Specify search volume %-------------------------------------------------------------------------- if isfield(xY,'def') @@ -176,7 +176,7 @@ str = spm_file(D.fname,'short30'); str = regexprep(str, {'\\' '\^' '_' '{' '}'}, ... {'\\\\' '\\^' '\\_' '\\{' '\\}'}); % Escape TeX special characters - str = sprintf('image mask: %s',str); + str = sprintf('image mask: %s',str); VOX = sqrt(sum(D.mat(1:3,1:3).^2)); XYZ = D.mat \ [xSwE.XYZmm; ones(1, size(xSwE.XYZmm, 2))]; j = find(spm_sample_vol(D, XYZ(1,:), XYZ(2,:), XYZ(3,:),0) > 0); diff --git a/swe_batch.m b/swe_batch.m index a3771b4..80f4304 100644 --- a/swe_batch.m +++ b/swe_batch.m @@ -1,7 +1,7 @@ function swe_batch % This function prepares and launches the batch system. % ========================================================================= -% This builds the whole tree for the various tools and their GUI at +% This builds the whole tree for the various tools and their GUI at % the first call to this script. % ========================================================================= % FORMAT bch = swe_batch diff --git a/swe_boxCoxTransform.m b/swe_boxCoxTransform.m index 9ed3bbf..f543009 100644 --- a/swe_boxCoxTransform.m +++ b/swe_boxCoxTransform.m @@ -11,7 +11,7 @@ % ========================================================================= % Bryan Guillaume % Version Info: $Format:%ci$ $Format:%h$ - + if min(data) < 0 error('The data must be positive.') end @@ -22,4 +22,4 @@ boxCoxTransformedData = (data .^ boxCoxLambda - 1) / boxCoxLambda; end - boxCoxTransformedData(data == 0) = -Inf; \ No newline at end of file + boxCoxTransformedData(data == 0) = -Inf; diff --git a/swe_cfg_batch.m b/swe_cfg_batch.m index 0bde3fb..b6a7c81 100644 --- a/swe_cfg_batch.m +++ b/swe_cfg_batch.m @@ -5,7 +5,7 @@ % ========================================================================= % Written by Bryan Guillaume % Version Info: $Format:%ci$ $Format:%h$ - + toolboxDir = spm_str_manip(mfilename('fullpath'), 'h'); addpath(toolboxDir); addpath(fullfile(toolboxDir, 'test')); @@ -18,7 +18,7 @@ 'Method for Neuroimaging Longitudinal Data Analysis.'] }'; swe.values = {swe_config_smodel swe_config_rmodel swe_config_results}; - + %---------------------------------------------------------------------- end diff --git a/swe_checkCompat.m b/swe_checkCompat.m index 4f7a4fc..cf25030 100644 --- a/swe_checkCompat.m +++ b/swe_checkCompat.m @@ -4,7 +4,7 @@ function swe_checkCompat(matVer, tbVer) % ========================================================================= % FORMAT: swe_checkCompat(matVer, tbVer) % ------------------------------------------------------------------------- -% Inputs: +% Inputs: % - matVer: version of SwE recorded in the `.mat` file. % - tbVer: version of SwE toolbox currently being run. % ========================================================================= @@ -12,16 +12,16 @@ function swe_checkCompat(matVer, tbVer) % Version Info: $Format:%ci$ $Format:%h$ if isequal(matVer,tbVer) - return; + return; end - + % The below hashmap records the earliest compatible version for each % release of the SwE toolbox. I.e. when making a new release, say you % are releasing version "y.y.y", please set earliestCompatVer("y.y.y") % equal to "x.x.x" where "x.x.x" is the oldest version of the toolbox % which "y.y.y" can accept `SwE.mat` files from. earliestCompatVer = containers.Map(); - + % These versions did not record version numbers in the SwE.mat file and % therefore cannot be checked in the same way. However, none of these % should be compatible with version 2.0.1 (the version in which this @@ -40,7 +40,7 @@ function swe_checkCompat(matVer, tbVer) earliestCompatVer('1.2.9') = 'NaN.NaN.NaN'; earliestCompatVer('1.2.10') = 'NaN.NaN.NaN'; earliestCompatVer('1.2.11') = 'NaN.NaN.NaN'; - + % Record earliest compatible versions. earliestCompatVer('2.0.0') = '2.0.0'; earliestCompatVer('2.1.0') = '2.0.0'; @@ -54,7 +54,7 @@ function swe_checkCompat(matVer, tbVer) % earliest compatible versions. This code is now redundant but may be % useful in future so has been left in place. Tom Maullin (09/11/2018) %latestCompatVer = latComVer(earliestCompatVer); - + % Check if the `.mat` version is compatible with this version. if ~strcmp(earliestCompatVer(matVer), earliestCompatVer(tbVer)) || ... strcmp(earliestCompatVer(matVer), 'NaN.NaN.NaN') @@ -63,7 +63,7 @@ function swe_checkCompat(matVer, tbVer) 'sion ', tbVer, '). Please re-enter the job specification ',... 'in the batch window.']); end - + end %-------------------------------------------------------------------------- @@ -79,7 +79,7 @@ function swe_checkCompat(matVer, tbVer) vers = ecv.keys; for i = 1:length(vers) ver = vers{i}; - + % If it's an old version we have no recorded version history for % earliest or latest compatible versions. if isequal(ecv(ver), 'NaN.NaN.NaN') @@ -96,14 +96,14 @@ function swe_checkCompat(matVer, tbVer) end end end - + % Second pass for keys that aren't values in ecv. for i = 1:length(vers) ver = vers{i}; - + if ~isKey(lcv, ver) lcv(ver) = lcv(ecv(ver)); end end - -end \ No newline at end of file + +end diff --git a/swe_cifti_clusters.m b/swe_cifti_clusters.m index 89dd4e8..6aedb5e 100644 --- a/swe_cifti_clusters.m +++ b/swe_cifti_clusters.m @@ -68,4 +68,4 @@ clusterSizesInVolume = [clusterSizesInVolume, histc(tmp,1:nClusters)]; end end -end \ No newline at end of file +end diff --git a/swe_cifti_convertVolLabels.m b/swe_cifti_convertVolLabels.m index 8dc53a9..ff199e5 100644 --- a/swe_cifti_convertVolLabels.m +++ b/swe_cifti_convertVolLabels.m @@ -1,7 +1,7 @@ function convertedLabels = swe_cifti_convertVolLabels(originalLabels, abbreviate) % Convert CIFTI struture labels from long labels to abbreviations or convert abbreviations to long labels % FORMAT convertedLabels = swe_cifti_convertVolLabels(originalLabels, abbreviate) - % originalLabels - a cell array of cifti volume labels or a single cifti volume label + % originalLabels - a cell array of cifti volume labels or a single cifti volume label % abbreviate - a boolean indicating if the labels needs to be abbreviated (true) or if the full labels need to be given (false) % convertedLabels - a cell array of converted cifti volume labels or a single cifti converted volume label ciftiVolumeLabels = {... @@ -55,4 +55,4 @@ if strcmp(originalClass, 'char') convertedLabels = char(convertedLabels(1)); -end \ No newline at end of file +end diff --git a/swe_cifti_max.m b/swe_cifti_max.m index f83091d..bf09058 100644 --- a/swe_cifti_max.m +++ b/swe_cifti_max.m @@ -19,9 +19,9 @@ % Bryan Guillaume % Version Info: $Format:%ci$ $Format:%h$ N = []; N_area = []; N_boxcox = []; - Z = []; M = []; A = []; XYZ = []; Mmm = []; + Z = []; M = []; A = []; XYZ = []; Mmm = []; brainStructureShortLabels = []; brainStructureLongLabels = []; - + % for retro-compatibility if isfield(ciftiInformation, 'isClusConstrainedInVolROI') isClusConstrainedInVolROI = ciftiInformation.isClusConstrainedInVolROI; @@ -127,7 +127,7 @@ Z = [Z, Z_tmp]; end % need to convert the volume coordinates into CIfTI coordinates - isMax = ismember(ciftiInformation.volume.XYZ', M_tmp', 'rows')'; + isMax = ismember(ciftiInformation.volume.XYZ', M_tmp', 'rows')'; M = [M, [ciftiInformation.volume.indices(isMax); ones(2,sum(isMax))]]; A = [A; A_tmp + maxA]; % need to convert the volume coordinates into CIfTI coordinates @@ -163,7 +163,7 @@ Z = [Z, Z_tmp]; end % need to convert the volume coordinates into CIfTI coordinates - isMax = ismember(ciftiInformation.volume.XYZ', M_tmp', 'rows')'; + isMax = ismember(ciftiInformation.volume.XYZ', M_tmp', 'rows')'; M = [M, [ciftiInformation.volume.indices(isMax); ones(2,sum(isMax))]]; A = [A; A_tmp + maxA]; % need to convert the volume coordinates into CIfTI coordinates @@ -191,4 +191,4 @@ maxA = max(A); end end -end \ No newline at end of file +end diff --git a/swe_compareVersions.m b/swe_compareVersions.m index 48f75a4..589f50f 100644 --- a/swe_compareVersions.m +++ b/swe_compareVersions.m @@ -3,36 +3,36 @@ % ========================================================================= % FORMAT: result = swe_compareVersions(ver1, ver2, comparisonOperator) % ------------------------------------------------------------------------- -% Inputs: +% Inputs: % - ver1: first version number of the type 'x.x.x' % - ver2: second version number of the type 'x.x.x' % - comparisonOperator: comparison operator ('<', '<=', '==', '>=' or '>') % ========================================================================= % Author: Bryan Guillaume & Tom Maullin (05/07/2019) % Version Info: $Format:%ci$ $Format:%h$ - + % Split version number into parts. ver1 = strsplit(ver1, '.'); ver2 = strsplit(ver2, '.'); - + if numel(ver1) > 3 ver1 = ver1(1:3); end - + if numel(ver2) > 3 ver2 = ver2(1:3); end % Work out base for comparison. base = max(cellfun(@(a) str2num(a), {ver1{:} ver2{:}})) + 1; - + % Convert to num. ver1 = cellfun(@(a) str2num(a), {ver1{:}}); ver2 = cellfun(@(a) str2num(a), {ver2{:}}); - + % Make into number in base `base` ver1 = ver1(1)*base^2 + ver1(2)*base^1 + ver1(3)*base^0; ver2 = ver2(1)*base^2 + ver2(2)*base^1 + ver2(3)*base^0; - + if strcmp(comparisonOperator, '<=') result = (ver1 <= ver2); elseif strcmp(comparisonOperator, '<') @@ -46,4 +46,4 @@ else error('unknown comparison operator') end -end \ No newline at end of file +end diff --git a/swe_config_results.m b/swe_config_results.m index 97131b2..201ac78 100644 --- a/swe_config_results.m +++ b/swe_config_results.m @@ -27,4 +27,4 @@ results.val = {dis}; results.help = {' ' 'Module of the SwE toolbox allowing the display of a previously computed model.'}; -results.prog = @swe_run_results; \ No newline at end of file +results.prog = @swe_run_results; diff --git a/swe_config_smodel.m b/swe_config_smodel.m index bb02879..16ec8fe 100644 --- a/swe_config_smodel.m +++ b/swe_config_smodel.m @@ -104,7 +104,7 @@ ciftiGeomFiles.num = [0 Inf]; % --------------------------------------------------------------------- -% volRoiConstraint +% volRoiConstraint % --------------------------------------------------------------------- volRoiConstraint = cfg_menu; volRoiConstraint.tag = 'volRoiConstraint'; @@ -170,7 +170,7 @@ 'Enter the groups in the ordering of the scans.' 'The groups correspond to subjects sharing a common covariance matrix.' 'The common covariance matrices are allowed to be different across groups.' - + ' '}'; groups.strtype = 'e'; groups.num = [Inf 1]; @@ -473,7 +473,7 @@ ' This choice tends to overestimate the degrees of freedom in some designs, but reduce the quantity of images saved and the computation time.' 'approx I: degrees of freedom estimation with the estimate proposed in Guillaume et al. (2014).' ' This estimate assumes no missing data and does not correct for the presence of a small sample bias and a missing data bias.' - ' Simulations seems to show that the estimate approx II and approx III are better choices (see below).' + ' Simulations seems to show that the estimate approx II and approx III are better choices (see below).' 'approx II: degrees of freedom estimation with an alternative estimate proposed in Guillaume (2015).' ' The estimate accounts partially for the presence of missing data and for a small-sample bias, but does not account for a missing data bias.' ' Simulations seems to indicate that it performs better than approx I, but should be used only under no missing data' @@ -482,7 +482,7 @@ ' Simulations seems to indicate that it systematically performs better than approx II under missing data, but seems slightly less well (slightly conservative) than approx II under no missing data.' ' That is the recommended choice by default. Nevertheless, if there is no missing data, approx II could be selected instead.' ' ' - }'; + }'; % --------------------------------------------------------------------- % modified Modified % --------------------------------------------------------------------- @@ -869,11 +869,11 @@ '' 'Please enter a mask/''.mat'' file containing mesh data.'}; -% --------------------------------------------------------------------- +% --------------------------------------------------------------------- % WB_voxelwise WB clusterwise volumetric .img input % --------------------------------------------------------------------- WB_volumetric = cfg_const; -WB_volumetric.tag = 'WB_volumetric'; +WB_volumetric.tag = 'WB_volumetric'; WB_volumetric.name = 'Volumetric input'; WB_volumetric.val = {0}; WB_volumetric.help = {'' @@ -893,7 +893,7 @@ '' 'Please select whether the ''.mat'' file contains volumetric or surface data.'}; -% --------------------------------------------------------------------- +% --------------------------------------------------------------------- % WB_voxelwise WB clusterwise .img input % --------------------------------------------------------------------- WB_img = cfg_const; diff --git a/swe_conman.m b/swe_conman.m index 4b13909..e5f4d3a 100644 --- a/swe_conman.m +++ b/swe_conman.m @@ -19,8 +19,8 @@ %========================================================================== % % The contrast manager graphicsl user interface (GUI) is a dialog box for -% the specification and selection of contrasts. The contrast selection -% interface is presented initially, pressing the "Define new contrast..." +% the specification and selection of contrasts. The contrast selection +% interface is presented initially, pressing the "Define new contrast..." % button pops up the contrast definition interface: % %__________________________________________________________________________ @@ -28,7 +28,7 @@ % % The contrast selection interface consists of: % -% * "Show" statistic type radio-buttons: "t-contrasts", "F-contrasts", +% * "Show" statistic type radio-buttons: "t-contrasts", "F-contrasts", % "All": % Selects the types of contrast that are presented for selection in the % contrast selection list box. Depending on the use of the contrast @@ -46,7 +46,7 @@ % Selected contrasts are highlit in black. % % * Image of the design matrix: -% A grey-scale representation of the design matrix, to aid +% A grey-scale representation of the design matrix, to aid % interpretation of the contrasts. % % The design matrix is "surfable": Clicking (and holding or dragging) @@ -149,7 +149,7 @@ % spm_input.m for more tips on using evaluated input widgets.) % % * For F-contrasts, a "columns for reduced design" edit widget appears: -% - Here you can specify SwE{F} maps by specifying the reduced design +% - Here you can specify SwE{F} maps by specifying the reduced design % as a sub-partition of the current design matrix. % - Enter the indicies of the design matrix columns you wish to retain % in the reduced design. @@ -247,7 +247,7 @@ % Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging % Andrew Holmes -% Based on: swe_conman.m 4617 2012-01-11 15:46:16Z will +% Based on: swe_conman.m 4617 2012-01-11 15:46:16Z will %========================================================================== @@ -458,11 +458,11 @@ if nargin<3, n=1; else n=varargin{3}; end if nargin<2, STATmode='T|F'; else STATmode=varargin{2}; end if nargin<1, error('no SwE struct specified'); else SwE = varargin{1};end - + %-Change to results directory %---------------------------------------------------------------------- try, cd(SwE.swd); end - + %---------------------------------------------------------------------- try SwE.xCon; @@ -481,7 +481,7 @@ %-Copy SwE structure for later change detecton, eg renamed contrasts tmpSPM = SwE; - + %-Wait until filenames have been selected hDone = findobj(F,'Tag','Done'); waitfor(hDone,'UserData') @@ -496,7 +496,7 @@ I = Q(get(hConList,'Value')); SwE = get(F,'UserData'); xCon = SwE.xCon; - + %-Save SwE.mat only if SwE structure as changed if ~isequal(tmpSPM,SwE) if spm_check_version('matlab','7') >=0 @@ -548,7 +548,7 @@ close(F) varargout = {[],cF}; return - + %------------------------------------------------------------------ case {'off','reset'} %------------------------------------------------------------------ @@ -559,7 +559,7 @@ set(findobj(F,'Tag','Done'),'UserData',-1) end return - + %------------------------------------------------------------------ case 'on' %------------------------------------------------------------------ @@ -1274,7 +1274,7 @@ if fcn==1 %-Parse string from ConMtx widget %------------------------------------------------------------------ set(hD_X1cols,'String','') - + [c,I,emsg,imsg] = swe_conman('ParseCon',str,xX.X,STAT); if all(I) DxCon = spm_FcUtil('Set','',STAT,'c',c,xX.X); @@ -1368,12 +1368,12 @@ DxCon.Vspm2 = []; DxCon.edf = []; DxCon.VspmUncP = []; - + if ~strcmp(DxCon.STAT,STAT), error('STAT & DxCon.STAT mismatch!'), end SwE = get(F,'UserData'); xCon = SwE.xCon; - - % Due to an earlier bug in SPM, a spurious field + + % Due to an earlier bug in SPM, a spurious field % PSTAT may have been created. Remove this if it exists. if isfield(xCon,'PSTAT') xCon=rmfield(xCon,'PSTAT'); @@ -1386,7 +1386,7 @@ SwE.xCon = xCon; set(F,'UserData',SwE); - + %-Redisplay the new list of contrasts, with the new one selected %------------------------------------------------------------------ hConList = findobj(F,'Tag','ConList'); @@ -1406,7 +1406,7 @@ set(get(findobj(gcbf,'Tag','DefineNew'),'UserData'),'Visible','off') set(gcbf,'UIContextMenu',findobj(gcbf,'Tag','ConMan_ConMen')) - + %====================================================================== case 'd_status' %====================================================================== @@ -1511,7 +1511,7 @@ 'HorizontalAlignment','Left',... 'BackgroundColor','w',... 'Position',[042 320 256 016].*WS); - + hh = uicontextmenu('Tag','ConMan_ConMenu'); uimenu(hh,'Label','Rename selected contrast...',... 'Tag','CM_Ren',... @@ -1519,13 +1519,13 @@ 'Interruptible','off','BusyAction','Cancel'); set(hh,'CallBack','swe_conman(''FConMenu_CB'')',... 'Interruptible','off','BusyAction','Cancel'); - + hConList = uicontrol(F,'Style','ListBox','Tag','ConList',... 'ToolTipString',['Select contrast(s) - drag/shift-click',... '/ctrl-click to select multiple contrasts'],... 'String',{'list','not','set','yet'},... 'Max',2,... - 'CallBack','swe_conman(''ConList_CB'')',... + 'CallBack','swe_conman(''ConList_CB'')',... 'UIContextMenu',hh,... 'Interruptible','off','BusyAction','Queue',... 'BackgroundColor','w',... diff --git a/swe_contrasts.m b/swe_contrasts.m index e8462ac..a456610 100644 --- a/swe_contrasts.m +++ b/swe_contrasts.m @@ -1,10 +1,10 @@ function [SwE] = swe_contrasts(SwE,Ic) % This function writes statistic & p value images for parametric analyses. % ========================================================================= -% Fills in SwE.xCon and writes the following images for a parametric +% Fills in SwE.xCon and writes the following images for a parametric % analysis: % - swe_{vox|dpx|dat}_{T|F}stat_c{c#}.{nii|gii|mat}: T/F tatistic image -% - swe_{vox|dpx|dat}_{zT|xF}stat_c{c#}.{nii|gii|mat}: Z/Chi square equivalent +% - swe_{vox|dpx|dat}_{zT|xF}stat_c{c#}.{nii|gii|mat}: Z/Chi square equivalent % statistic image % - swe_{vox|dpx|dat}_{T|F}stat_lp_c{c#}.{nii|gii|mat}: -log10 P image. % - swe_{vox|dpx|dat}_edf_c{c#}.{nii|gii|mat}: error degrees of freedom image. @@ -32,7 +32,7 @@ cd(SwE.swd); end -% For Wild Bootstrap we already made contrast images, we just need to +% For Wild Bootstrap we already made contrast images, we just need to % record them, if we haven't already. %-------------------------------------------------------------------------- if isfield(SwE, 'WB') @@ -64,7 +64,7 @@ isOctave = exist('OCTAVE_VERSION','builtin'); if isCifti - metadata = {'ciftiTemplate', SwE.xY.P{1}}; + metadata = {'ciftiTemplate', SwE.xY.P{1}}; file_data_type = 'dpx'; end @@ -97,7 +97,7 @@ %-Map parameter files %-------------------------------------------------------------------------- - + %-OLS estimators and covariance estimates %-------------------------------------------------------------------------- Vbeta = SwE.Vbeta; @@ -116,7 +116,7 @@ XYZ = SwE.xVol.XYZ; S=size(XYZ,2); for i = 1:length(Ic) - + %-Canonicalise contrast structure with required fields %---------------------------------------------------------------------- ic = Ic(i); @@ -138,13 +138,13 @@ ind = find(any(Co'~=0)); end nCov_beta = (nBeta+1)*nBeta/2; - + % if ".mat" format, load data now if isMat beta = importdata(Vbeta); S = size(beta,2); end - + %-Compute contrast %------------------------------------------------------ fprintf('\t%-32s: %30s',sprintf('contrast image %2d',ic),... @@ -159,14 +159,14 @@ cBeta = zeros(1,S); for j=1:numel(ind) if isMat - cBeta = cBeta + Co(ind(j)) * beta(ind(j),:); + cBeta = cBeta + Co(ind(j)) * beta(ind(j),:); else cBeta = cBeta + Co(ind(j)) * swe_data_read(V(j),'xyz',XYZ); end swe_progress_bar('Set',100*(j/numel(ind))); end swe_progress_bar('Clear') - + if isMat %-save contrasted beta %------------------------------------------------------ @@ -184,7 +184,7 @@ tmp = zeros(SwE.xVol.DIM'); tmp(Q) = cBeta; xCon(ic).Vcon = swe_data_write(xCon(ic).Vcon, tmp); - + clear tmp end if ~isMat @@ -192,7 +192,7 @@ '...written %s', xCon(ic).Vcon.fname)) else fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),sprintf(... - '...written %s', xCon(ic).Vcon)) + '...written %s', xCon(ic).Vcon)) end else if ~isMat @@ -201,7 +201,7 @@ cBeta = zeros(nSizeCon,S); for j=1:numel(ind) if isMat - cBeta = cBeta + Co(ind(j),:)' * beta(ind(j),:); + cBeta = cBeta + Co(ind(j),:)' * beta(ind(j),:); else cBeta = cBeta + Co(ind(j),:)' * swe_data_read(V(j),'xyz',XYZ); end @@ -213,17 +213,17 @@ '...computed')) else fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),sprintf(... - '...computed')) + '...computed')) end - end - + end + % clear beta for memory if isMat clear beta end %-Write inference SwE %====================================================================== - + %-compute the contrasted beta covariances and edof for the contrast switch(xCon(ic).STAT) case 'T' @@ -234,7 +234,7 @@ fprintf('\t%-32s: %30s',sprintf('spm{%c} image %2d',eSTAT,ic),... '...computing'); %-# str = 'contrasted beta covariance computation'; - swe_progress_bar('Init',100,str,''); + swe_progress_bar('Init',100,str,''); it = 0; it2 = 0; @@ -255,7 +255,7 @@ end xCon(ic).edf = sum(SwE.dof.edof_Subj(indSubjInvolved)); end - + % load .mat file(s) if this is the format if isMat cov_beta = importdata(Vcov_beta); @@ -263,19 +263,19 @@ cov_beta_g = importdata(Vcov_beta_g); end end - + for j = 1:nBeta for jj = j:nBeta it = it + 1; if any(j == ind) && any(jj == ind) it2 = it2+1; weight = Co(j,:)'*Co(jj,:); - if (j~=jj) %was wrong (BG - 13/09/13) + if (j~=jj) %was wrong (BG - 13/09/13) weight = weight + weight'; end weight = weight(tril(ones(nSizeCon))==1); if isMat - cCovBc = cCovBc + weight * cov_beta(it,:); + cCovBc = cCovBc + weight * cov_beta(it,:); else cCovBc = cCovBc + weight * swe_data_read(Vcov_beta(it),'xyz',XYZ); end @@ -307,7 +307,7 @@ %---------------------------------------------------------- score = cBeta ./ sqrt(cCovBc); swe_progress_bar('Set',100*(0.1)); - switch dof_type + switch dof_type case 1 tmp = 0; for g = 1:SwE.Gr.nGr @@ -320,7 +320,7 @@ % transform into Z-scores image if any(score>0) % avoid to run the following line when all Z are < 0 (BG - 22/08/2016) uncP(score>0) = spm_Tcdf(-score(score>0),edf(score>0)); - equivalentScore(score>0) = -swe_invNcdf(uncP(score>0)); + equivalentScore(score>0) = -swe_invNcdf(uncP(score>0)); end if any(score<0) % avoid to run the following line when all Z are > 0(BG - 22/08/2016) uncP(score<0) = spm_Tcdf(score(score<0),edf(score<0)); @@ -359,11 +359,11 @@ swe_progress_bar('Set',100*(0.1) + g*80/SwE.Gr.nGr); end clear Wg cov_vis - edf = 2 * cCovBc.^2 ./ CovcCovBc - 2; + edf = 2 * cCovBc.^2 ./ CovcCovBc - 2; clear CovcCovBc if any(score>0) % avoid to run the following line when all Z are < 0 (BG - 22/08/2016) uncP(score>0) = spm_Tcdf(-score(score>0),edf(score>0)); - equivalentScore(score>0) = -swe_invNcdf(uncP(score>0)); + equivalentScore(score>0) = -swe_invNcdf(uncP(score>0)); end if any(score<0) % avoid to run the following line when all Z are > 0(BG - 22/08/2016) uncP(score<0) = spm_Tcdf(score(score<0),edf(score<0)); @@ -373,7 +373,7 @@ %Z = -log10(1-spm_Tcdf(Z,edf)); %transfo into -log10(p) swe_progress_bar('Set',100); case 3 - CovcCovBc = 0; + CovcCovBc = 0; if isMat cov_vis = importdata(Vcov_vis); end @@ -393,7 +393,7 @@ % transform into Z-scores image if any(score>0) % avoid to run the following line when all Z are < 0 (BG - 22/08/2016) uncP(score>0) = spm_Tcdf(-score(score>0),edf(score>0)); - equivalentScore(score>0) = -swe_invNcdf(uncP(score>0)); + equivalentScore(score>0) = -swe_invNcdf(uncP(score>0)); end if any(score<0) % avoid to run the following line when all Z are > 0(BG - 22/08/2016) uncP(score<0) = spm_Tcdf(score(score<0),edf(score<0)); @@ -402,8 +402,8 @@ end %Z = -log10(1-spm_Tcdf(Z,edf)); %transfo into -log10(p) swe_progress_bar('Set',100); - end - + end + case 'F' %-Compute spm{F} image %--------------------------------------------------------- if nSizeCon==1 @@ -451,7 +451,7 @@ swe_progress_bar('Set',100*(g/SwE.Gr.nGr/10+0.1)); end clear Wg cov_vis - edf = 2 * cCovBc.^2 ./ CovcCovBc - 2; + edf = 2 * cCovBc.^2 ./ CovcCovBc - 2; clear CovcCovBc swe_progress_bar('Set',100*(3/4)); % transform into X-scores image @@ -470,12 +470,12 @@ Wg = kron(Co,Co)' * swe_duplication_matrix(nBeta) * SwE.Vis.weight(:,SwE.Vis.iGr_Cov_vis_g==g); Wg = kron(Wg,Wg) * swe_duplication_matrix(SwE.Vis.nCov_vis_g(g)); if isMat - CovcCovBc = CovcCovBc + Wg * swe_vechCovVechV(cov_vis(SwE.Vis.iGr_Cov_vis_g==g,:),SwE.dof.dofMat{g},2); + CovcCovBc = CovcCovBc + Wg * swe_vechCovVechV(cov_vis(SwE.Vis.iGr_Cov_vis_g==g,:),SwE.dof.dofMat{g},2); else - CovcCovBc = CovcCovBc + Wg * swe_vechCovVechV(swe_data_read(Vcov_vis(SwE.Vis.iGr_Cov_vis_g==g),'xyz',XYZ),SwE.dof.dofMat{g},2); + CovcCovBc = CovcCovBc + Wg * swe_vechCovVechV(swe_data_read(Vcov_vis(SwE.Vis.iGr_Cov_vis_g==g),'xyz',XYZ),SwE.dof.dofMat{g},2); end - swe_progress_bar('Set',100*(g/SwE.Gr.nGr/10+0.1)); - end + swe_progress_bar('Set',100*(g/SwE.Gr.nGr/10+0.1)); + end clear Wg cov_vis edf = 2 * cCovBc.^2 ./ CovcCovBc; % transform into X-scores image @@ -489,7 +489,7 @@ end % need to transform in F-score, not in absolute t-score % corrected on 12/05/15 by BG - score = score.^2; + score = score.^2; else score = zeros(1,S); if dof_type ~= 0 @@ -545,19 +545,19 @@ cCovBc_vox = zeros(nSizeCon); cCovBc_vox(tril(ones(nSizeCon))==1) = cCovBc(:,iVox); cCovBc_vox = cCovBc_vox + cCovBc_vox' - diag(diag(cCovBc_vox)); - score(iVox) = cBeta(:,iVox)' / cCovBc_vox * cBeta(:,iVox); - if (dof_type == 1) + score(iVox) = cBeta(:,iVox)' / cCovBc_vox * cBeta(:,iVox); + if (dof_type == 1) tmp = 0; for g = 1:SwE.Gr.nGr cCovBc_g_vox = zeros(nSizeCon); cCovBc_g_vox(tril(ones(nSizeCon))==1) = cCovBc_g(:,iVox,g); cCovBc_g_vox = cCovBc_g_vox + cCovBc_g_vox' - diag(diag(cCovBc_g_vox)); tmp = tmp + (trace(cCovBc_g_vox^2) + (trace(cCovBc_g_vox))^2)/... - SwE.dof.edof_Gr(g); + SwE.dof.edof_Gr(g); end - edf(iVox)=(trace(cCovBc_vox^2) + (trace(cCovBc_vox))^2) / tmp; + edf(iVox)=(trace(cCovBc_vox^2) + (trace(cCovBc_vox))^2) / tmp; end - % update progress_bar only approx 80 times + % update progress_bar only approx 80 times if (mod(iVox,updateEvery) == 0) swe_progress_bar('Set',10 + 80 * (iVox/S)); end @@ -566,7 +566,7 @@ if dof_type ~= 0 clear cCovBc_g score = score .*(edf-xCon(ic).eidf+1)./edf/xCon(ic).eidf; - score(score < 0) = 0; % force negatif F-score to 0 (can happen for very low edf) + score(score < 0) = 0; % force negatif F-score to 0 (can happen for very low edf) % transform into X-scores image % Z2 = chi2inv(spm_Fcdf(Z,xCon(ic).eidf,edf-xCon(ic).eidf+1),1); try % check if the user do not have the fcdf function or one with 'upper' option @@ -574,9 +574,9 @@ equivalentScore(score>1) = swe_invNcdf(uncP(score>1)/2).^2; uncP(score<=1 & score > 0) = 1 - fcdf(score(score<=1 & score > 0),xCon(ic).eidf,edf(score<=1 & score > 0)-xCon(ic).eidf+1); equivalentScore(score<=1 & score > 0) = swe_invNcdf(uncP(score<=1 & score > 0)/2).^2; - catch + catch uncP(score>0) = betainc((edf(score>0) - xCon(ic).eidf + 1)./(edf(score>0) - xCon(ic).eidf + 1 + xCon(ic).eidf * score(score>0)),(edf(score>0)-xCon(ic).eidf+1)/2, xCon(ic).eidf/2); % more accurate to look this way for high scores - equivalentScore(score>0) = swe_invNcdf(uncP(score>0)/2).^2; + equivalentScore(score>0) = swe_invNcdf(uncP(score>0)/2).^2; % Z2 = swe_invNcdf(0.5 - spm_Fcdf(Z,xCon(ic).eidf, edf-xCon(ic).eidf+1)/2).^2; end uncP(score == 0) = 0; @@ -585,7 +585,7 @@ %Z = -log10(1-spm_Fcdf(Z,xCon(ic).eidf,edf)); else score = score *(xCon(ic).edf -xCon(ic).eidf+1)/xCon(ic).edf/xCon(ic).eidf; - score(score < 0) = 0; % force negatif F-score to 0 (can happen for very low edf) + score(score < 0) = 0; % force negatif F-score to 0 (can happen for very low edf) % transform into X-scores image %Z2 = chi2inv(spm_Fcdf(Z,xCon(ic).eidf,xCon(ic).edf-xCon(ic).eidf+1),1); try % check if the user do not have the fcdf function or one with 'upper' options @@ -595,7 +595,7 @@ equivalentScore(score<=1 & score > 0) = swe_invNcdf(uncP(score<=1 & score > 0)/2).^2; catch uncP(score>0) = betainc((xCon(ic).edf - xCon(ic).eidf + 1)./(xCon(ic).edf - xCon(ic).eidf + 1 + xCon(ic).eidf * score(score>0)),(xCon(ic).edf-xCon(ic).eidf+1)/2, xCon(ic).eidf/2); - equivalentScore(score>0) = swe_invNcdf(uncP(score>0)/2).^2; + equivalentScore(score>0) = swe_invNcdf(uncP(score>0)/2).^2; % Z2(Z>0) = swe_invNcdf(0.5 - spm_Fcdf(Z,xCon(ic).eidf,xCon(ic).edf-xCon(ic).eidf+1)/2).^2; end uncP(score == 0) = 0; @@ -609,15 +609,15 @@ luncP = -log10(uncP); swe_progress_bar('Clear') clear cCovBc cB tmp - - + + %-Write SwE - statistic images & edf image if needed %------------------------------------------------------------------ fprintf('%s%30s',repmat(sprintf('\b'),1,30),'...writing'); %-# - + equivalentScore (equivalentScore > realmax('single')) = realmax('single'); equivalentScore (equivalentScore < -realmax('single')) = -realmax('single'); - + if isMat xCon(ic).Vspm = sprintf('swe_%s_%c%cstat_c%02d%s',file_data_type,lower(eSTAT),xCon(ic).STAT,ic, file_ext); save(xCon(ic).Vspm, 'equivalentScore'); @@ -639,9 +639,9 @@ fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),sprintf(... '...written %s',spm_str_manip(xCon(ic).Vspm.fname,'t'))); end - + fprintf('\t%-32s: %30s',sprintf('spm{%c} image %2d',xCon(ic).STAT,ic),... - '...writing'); + '...writing'); if isMat xCon(ic).Vspm2 = sprintf('swe_%s_%cstat_c%02d%s',file_data_type,xCon(ic).STAT,ic,file_ext); @@ -671,12 +671,12 @@ '...written %s',spm_str_manip(xCon(ic).Vspm2,'t'))); %-# else fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),sprintf(... - '...written %s',spm_str_manip(xCon(ic).Vspm2.fname,'t'))); %-# + '...written %s',spm_str_manip(xCon(ic).Vspm2.fname,'t'))); %-# end - + % save raw uncorrected p-values (new on 05/11/2017) fprintf('\t%-32s: %30s',sprintf('-log10(uncP) image %2d',ic),... - '...writing'); + '...writing'); if isMat xCon(ic).VspmUncP = sprintf('swe_%s_%cstat_lp_c%02d%s',file_data_type,xCon(ic).STAT,ic,file_ext); save(xCon(ic).VspmUncP, 'luncP'); @@ -707,7 +707,7 @@ fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),sprintf(... '...written %s',spm_str_manip(xCon(ic).VspmUncP.fname,'t'))); end - + fprintf('\t%-32s: %30s',sprintf('edf image %2d',ic),... '...writing'); @@ -743,7 +743,7 @@ '...written %s',spm_str_manip(xCon(ic).Vedf.fname,'t')))%-# end end - + end % if isempty(xCon(ic).Vspm) end % (for i = 1:length(Ic)) @@ -761,4 +761,4 @@ save('SwE','SwE','-V6'); else save('SwE','SwE'); -end \ No newline at end of file +end diff --git a/swe_contrasts_WB.m b/swe_contrasts_WB.m index ca18c16..e8c79a3 100644 --- a/swe_contrasts_WB.m +++ b/swe_contrasts_WB.m @@ -1,7 +1,7 @@ function [SwE] = swe_contrasts_WB(SwE) % This function creates the missing 'xCon' for wild bootstrap SwE objects. % ========================================================================= -% Note: Whilst contrasts are created in 'swe_contrasts.m', in +% Note: Whilst contrasts are created in 'swe_contrasts.m', in % 'swe_WB_contrasts.m' they are only recorded. They have already been % created by 'swe_cp_WB.m'. % ========================================================================= @@ -18,7 +18,7 @@ try cd(SwE.swd); end - + %-Get file extension %---------------------------------------------------------------------- file_ext = swe_get_file_extension(SwE.xY.P{1}); @@ -42,28 +42,28 @@ file_data_type = 'vox'; end end - + %-Retrieve contrast details. %---------------------------------------------------------------------- STAT = SwE.WB.stat; c = SwE.WB.con; - + %-Retrieve equivalent statistic details. %---------------------------------------------------------------------- if strcmpi(STAT, 't') eSTAT = 'z'; else eSTAT = 'x'; - end - + end + %-Retrieve design matrix. %---------------------------------------------------------------------- X = SwE.xX.X; - + %-Create the xCon object. %---------------------------------------------------------------------- DxCon = spm_FcUtil('Set', 'Contrast 1: Activation', STAT, 'c', c', X); - + % Work out if we are in clusterwise bootstrap or not. %---------------------------------------------------------------------- if isfield(SwE,'WB') @@ -71,7 +71,7 @@ else wbstring = ''; end - + % Add the SwE volumes. %---------------------------------------------------------------------- DxCon.Vspm = swe_data_hdr_read(sprintf('swe_%s_%c%cstat_c%.2d%s', file_data_type, eSTAT, STAT, 1, file_ext)); @@ -83,9 +83,9 @@ % In the most recent versions of the toolbox, VspmFWEP_clusnorm should be used, but, for models estimated % with older versions, it does not exist and VspmFWEP_clus should be used instead. try - DxCon.VspmFWEP_clusnorm = swe_data_hdr_read(sprintf('swe_clusternorm_%cstat_lpFWE%s_c%.2d%s', STAT, wbstring, 1, file_ext)); + DxCon.VspmFWEP_clusnorm = swe_data_hdr_read(sprintf('swe_clusternorm_%cstat_lpFWE%s_c%.2d%s', STAT, wbstring, 1, file_ext)); catch - DxCon.VspmFWEP_clus = swe_data_hdr_read(sprintf('swe_clustere_%cstat_lpFWE%s_c%.2d%s', STAT, wbstring, 1, file_ext)); + DxCon.VspmFWEP_clus = swe_data_hdr_read(sprintf('swe_clustere_%cstat_lpFWE%s_c%.2d%s', STAT, wbstring, 1, file_ext)); end end if isfield(SwE.WB, 'TFCE') @@ -95,20 +95,20 @@ end end DxCon.Vedf = swe_data_hdr_read(sprintf('swe_%s_edf_c%.2d%s', file_data_type, 1, file_ext)); - + % Add the positive contrast. %---------------------------------------------------------------------- if ~isfield(SwE, 'xCon') SwE.xCon = DxCon; else - SwE.xCon = [SwE.xCon DxCon]; + SwE.xCon = [SwE.xCon DxCon]; end - + % We need to add the negative contrast if we are doing a T test. if strcmpi(STAT, 't') - + DxCon = spm_FcUtil('Set', 'Contrast 2: Deactivation', STAT, 'c', -c', X); - + DxCon.Vspm = swe_data_hdr_read(sprintf('swe_%s_%c%cstat_c%.2d%s', file_data_type, eSTAT, STAT, 2, file_ext)); DxCon.VspmUncP = swe_data_hdr_read(sprintf('swe_%s_%cstat_lp%s_c%.2d%s', file_data_type, STAT, wbstring, 2, file_ext)); DxCon.VspmFDRP = swe_data_hdr_read(sprintf('swe_%s_%cstat_lpFDR%s_c%.2d%s', file_data_type, STAT, wbstring, 2, file_ext)); @@ -130,8 +130,8 @@ end end DxCon.Vedf = swe_data_hdr_read(sprintf('swe_%s_edf_c%.2d%s', file_data_type, 1, file_ext)); - - SwE.xCon = [SwE.xCon DxCon]; + + SwE.xCon = [SwE.xCon DxCon]; end -end \ No newline at end of file +end diff --git a/swe_cp.m b/swe_cp.m index d85018e..28d20cb 100644 --- a/swe_cp.m +++ b/swe_cp.m @@ -1,7 +1,7 @@ function swe_cp(SwE) % This function computes covariance and beta maps for parametric analyses. % ========================================================================= -% For a parametric SwE analysis with nifti input, this function computes +% For a parametric SwE analysis with nifti input, this function computes % the following maps: % % - swe_vox_mask: The mask image for the analysis. @@ -13,7 +13,7 @@ function swe_cp(SwE) % - swe_vox_cov_g{g#}_v{v1#}_v{v2#}: The covariance map between betas % {v1#} and {v2#} for group {g#}. % -% For a parametric SwE analysis with GIfTI or CIfTI inputs, this function computes +% For a parametric SwE analysis with GIfTI or CIfTI inputs, this function computes % the following maps: % % - swe_dpx_mask: The mask image for the analysis. @@ -60,7 +60,7 @@ function swe_cp(SwE) swd = fileparts(P{i}); load(fullfile(swd,'SwE.mat')); SwE.swd = swd; - + % detect if this is a WB analysis or a "standard analysis" if isfield(SwE, 'WB') swe_cp_WB(SwE); @@ -139,7 +139,7 @@ function swe_cp(SwE) %-Delete files from previous analyses %-------------------------------------------------------------------------- if exist(fullfile(SwE.swd,sprintf('swe_%s_mask%s',file_data_type,file_ext)),'file') == 2 - + str = {'Current directory contains SwE estimation files:',... 'pwd = ',SwE.swd,... 'Existing results will be overwritten!'}; @@ -150,7 +150,7 @@ function swe_cp(SwE) warning('Overwriting old results\n\t (pwd = %s) ',SwE.swd); %#ok end end - + files = {'^swe_.{3}_mask(\.dtseries)?(\.dscalar)?\..{3}$',... '^swe_.{3}_cov_b\d{2}_b\d{2}(\.dtseries)?(\.dscalar)?\..{3}$',... '^swe_.{3}_cov_vv(\.dtseries)?(\.dscalar)?\..{3}$',... @@ -172,18 +172,18 @@ function swe_cp(SwE) '^swe_clusternorm\d{0,1}_\w{1,2}stat_lp\w{0,3}-WB_c\d{2}(\.dtseries)?(\.dscalar)?\..{3}$',... '^swe_.{3}_resid_y\d{2,4}(\.dtseries)?(\.dscalar)?\..{3}$',... '^swe_.{3}_fit_y\d{2,4}(\.dtseries)?(\.dscalar)?\..{3}$'}; - + for i = 1:length(files) j = spm_select('List',SwE.swd,files{i}); for k = 1:size(j,1) spm_unlink(deblank(j(k,:))); end end - + %========================================================================== % - A N A L Y S I S P R E L I M I N A R I E S %========================================================================== - + %-Initialise %========================================================================== fprintf('%-40s: %30s','Initialising parameters','...computing'); %-# @@ -202,7 +202,7 @@ function swe_cp(SwE) case 0 corr = ones(nScan,1); case 1 - corr = sqrt(nScan/(nScan-nBeta)); % residual correction (type 1) + corr = sqrt(nScan/(nScan-nBeta)); % residual correction (type 1) case 2 corr = (1-diag(Hat)).^(-0.5); % residual correction (type 2) case 3 @@ -214,7 +214,7 @@ function swe_cp(SwE) tmp = I_Hat(iSubj==uSubj(i), iSubj==uSubj(i)); tmp = (tmp + tmp')/2; [tmpV, tmpE] = eig(tmp); - corr{i} = tmpV * diag(1./sqrt(diag(tmpE))) * tmpV'; + corr{i} = tmpV * diag(1./sqrt(diag(tmpE))) * tmpV'; end clear I_Hat tmp case 5 @@ -223,14 +223,14 @@ function swe_cp(SwE) for i = 1:nSubj tmp = I_Hat(iSubj==uSubj(i), iSubj==uSubj(i)); tmp = (tmp + tmp')/2; - corr{i} = inv(tmp); + corr{i} = inv(tmp); end clear I_Hat tmp end %-detect if the design matrix is separable (a little bit messy, but seems to do the job) % -iGr_dof = zeros(1,nScan); +iGr_dof = zeros(1,nScan); iBeta_dof = zeros(1,nBeta); it = 0; while ~all(iGr_dof) @@ -260,18 +260,18 @@ function swe_cp(SwE) ind1 = find(sum(tmp,1)>1,1); % detect the first column in common ind2 = find(tmp(:,ind1)==1); % detect the groups to be fused for ii = ind2' - iGr_dof(iGr_dof==uGr_dof(ii)) = ind2(1); % fuse the groups + iGr_dof(iGr_dof==uGr_dof(ii)) = ind2(1); % fuse the groups end end end -nSubj_dof = zeros(1,nGr_dof); +nSubj_dof = zeros(1,nGr_dof); for i = 1:nGr_dof % renumber to avoid gaps in the numbering iGr_dof(iGr_dof==uGr_dof(i)) = i; iBeta_dof(tmp(i,:)==1) = i; nSubj_dof(i) = length(unique(iSubj(iGr_dof==uGr_dof(i)))); end -pB_dof = zeros(1,nGr_dof); -for i=1:nBeta +pB_dof = zeros(1,nGr_dof); +for i=1:nBeta tmp=1; for ii=1:nSubj if length(unique(xX.X(iSubj==uSubj(ii)&iGr_dof'==iBeta_dof(i),i)))>1 @@ -295,38 +295,38 @@ function swe_cp(SwE) if isfield(SwE.type,'modified') dof_type = SwE.type.modified.dof_mo; else - dof_type = SwE.type.classic.dof_cl; + dof_type = SwE.type.classic.dof_cl; end if dof_type == 0 % so naive estimation is used dof_cov = zeros(1,nBeta); for i = 1:nBeta dof_cov(i) = nSubj_dof(iBeta_dof(i)) - ... - pB_dof(iBeta_dof(i)); + pB_dof(iBeta_dof(i)); end -end +end %-preprocessing for the modified SwE if isfield(SwE.type,'modified') iVis = SwE.Vis.iVis; iGr = SwE.Gr.iGr; - uGr = unique(iGr); + uGr = unique(iGr); nGr = length(uGr); - + % info specific for each group uVis_g = cell(1,nGr); % unique visits for each group nVis_g = zeros(1,nGr); % number of visits for each group uSubj_g = cell(1,nGr); % unique visits for each group nSubj_g = zeros(1,nGr); % number of visits for each group for g = 1:nGr - uVis_g{g} = unique(iVis(iGr==uGr(g))); + uVis_g{g} = unique(iVis(iGr==uGr(g))); nVis_g(g) = length(uVis_g{g}); iSubj_g = iSubj(iGr==uGr(g)); % Subject number for each subject in group for each visit uSubj_g{g} = unique(iSubj_g); % Unique subject numbers of subjects in group nSubj_g(g) = length(uSubj_g{g}); uSubj_g_tmp = uSubj_g{g}; - - for k = 1:nSubj_g(g) + + for k = 1:nSubj_g(g) % The number of visits for subject uSubj_g(k) vis_g_subj(k) = sum(iSubj_g==uSubj_g_tmp(k)); @@ -335,17 +335,17 @@ function swe_cp(SwE) max_nVis_g(g) = max(vis_g_subj); min_nVis_g(g) = min(vis_g_subj); - + clear vis_g_subj - + end - + nCov_vis_g = nVis_g.*(nVis_g+1)/2; % number of covariance elements to be estimated for each group - nCov_vis = sum(nCov_vis_g); % total number of covariance elements to be estimated - - % Flags matrices indicating which residuals have to be used for each covariance element + nCov_vis = sum(nCov_vis_g); % total number of covariance elements to be estimated + + % Flags matrices indicating which residuals have to be used for each covariance element Flagk = false(nCov_vis,nScan); % Flag indicating scans corresponding to visit k for each covariance element - Flagkk = false(nCov_vis,nScan); % Flag indicating scans corresponding to visit kk for each covariance element + Flagkk = false(nCov_vis,nScan); % Flag indicating scans corresponding to visit kk for each covariance element Ind_Cov_vis_diag = nan(1,sum(nVis_g)); % index of the diagonal elements Ind_Cov_vis_off_diag = nan(1,nCov_vis - sum(nVis_g)); % index of the off-diagonal elements Ind_corr_diag=nan(nCov_vis,2); % index of the 2 corresponding diagonal elements @@ -354,19 +354,19 @@ function swe_cp(SwE) for g = 1:nGr for k = 1:nVis_g(g) for kk = k:nVis_g(g) - it = it + 1; + it = it + 1; id = intersect(iSubj(iGr==uGr(g) & iVis==uVis_g{g}(k)),... - iSubj(iGr==uGr(g) & iVis==uVis_g{g}(kk))); % identifiaction of the subjects with both visits k & kk + iSubj(iGr==uGr(g) & iVis==uVis_g{g}(kk))); % identifiaction of the subjects with both visits k & kk Flagk(it,:) = ismember(iSubj,id) & iVis==uVis_g{g}(k); - Flagkk(it,:) = ismember(iSubj,id) & iVis==uVis_g{g}(kk); - if k==kk + Flagkk(it,:) = ismember(iSubj,id) & iVis==uVis_g{g}(kk); + if k==kk it2 = it2+1; it4 = it2; Ind_Cov_vis_diag(it2) = it; else it3 = it3 + 1; it4 = it4 + 1; - Ind_Cov_vis_off_diag(it3) = it; + Ind_Cov_vis_off_diag(it3) = it; end Ind_corr_diag(it,:) = [it2 it4]; iGr_Cov_vis_g(it) = g; @@ -379,12 +379,12 @@ function swe_cp(SwE) for j = 1:nBeta for jj = j:nBeta it=it+1; - for jjj = Ind_Cov_vis_diag + for jjj = Ind_Cov_vis_diag weight(it,jjj) = pX(j,Flagk(jjj,:))*pX(jj,Flagk(jjj,:))'; end - for jjj = Ind_Cov_vis_off_diag + for jjj = Ind_Cov_vis_off_diag weight(it,jjj) = pX(j,Flagk(jjj,:))*pX(jj,Flagkk(jjj,:))' + ... - pX(j,Flagkk(jjj,:))*pX(jj,Flagk(jjj,:))'; + pX(j,Flagkk(jjj,:))*pX(jj,Flagk(jjj,:))'; end end @@ -401,7 +401,7 @@ function swe_cp(SwE) tmp = tmp + 1/edof_Subj(uSubj == uSubj_g{g}(j)); end edof_Gr(g) = nSubj_g(g)^2/tmp; - end + end case {2,3} % compute a matrix containing the variables linked to the degrees of freedom (for test II and III) dofMat = cell(nGr,1); for g = 1:nGr @@ -446,12 +446,12 @@ function swe_cp(SwE) for i = 1:nSubj nVis_i(i) = sum(uSubj(i)==iSubj); end - nCov_vis = sum(nVis_i.*(nVis_i+1)/2); % total number of covariance elements to be estimated + nCov_vis = sum(nVis_i.*(nVis_i+1)/2); % total number of covariance elements to be estimated weight = NaN(nCov_beta,nCov_vis); Ind_Cov_vis_classic = NaN(1,nCov_vis); Indexk = NaN(1,nCov_vis); Indexkk = NaN(1,nCov_vis); - it = 0; + it = 0; for j = 1:nBeta for jj = j:nBeta it = it + 1; @@ -460,7 +460,7 @@ function swe_cp(SwE) ind_i=find(iSubj == uSubj(i)); for ii = 1:nVis_i(i) it2 = it2 + 1; - weight(it,it2) = pX(j,ind_i(ii))*pX(jj,ind_i(ii)); + weight(it,it2) = pX(j,ind_i(ii))*pX(jj,ind_i(ii)); Ind_Cov_vis_classic(it2) = i; Indexk(it2) = ind_i(ii); Indexkk(it2) = ind_i(ii); @@ -474,7 +474,7 @@ function swe_cp(SwE) end end end - end + end %-compute the effective dof from each homogeneous group (here, subject) if dof_type == 1 edof_Gr = edof_Subj; @@ -599,12 +599,12 @@ function swe_cp(SwE) % DIM, M, sprintf('spm_spm:ResI (%02d)', i),... % isMeshData); % end - fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...initialised'); %-# - % + fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...initialised'); %-# + % %========================================================================== % - F I T M O D E L & W R I T E P A R A M E T E R I M A G E S %========================================================================== - + %-Get explicit mask(s) %========================================================================== mask = true(DIM); @@ -650,10 +650,10 @@ function swe_cp(SwE) for iChunk=1:nbchunks chunk = chunks(iChunk):chunks(iChunk+1)-1; - + %-Report progress %====================================================================== - if iChunk > 1, fprintf(repmat(sprintf('\b'),1,72)); end %-# + if iChunk > 1, fprintf(repmat(sprintf('\b'),1,72)); end %-# fprintf('%-40s: %30s', sprintf('Chunk %3d/%-3d',iChunk,nbchunks),... '...processing'); @@ -661,16 +661,16 @@ function swe_cp(SwE) %------------------------------------------------------------------ Y = zeros(nScan, numel(chunk)); cmask = mask(chunk); - + if size(cmask, 2) == 1 cmask = cmask'; end - + for iScan=1:nScan if ~any(cmask), break, end %-Break if empty mask - + Y(iScan, cmask) = swe_data_read(VY(iScan), chunk(cmask));%-Read chunk of data - + cmask(cmask) = Y(iScan, cmask) > xM.TH(iScan); %-Threshold (& NaN) mask if xM.I && ~YNaNrep && xM.TH(iScan) < 0 %-Use implicit mask cmask(cmask) = abs(Y(iScan, cmask)) > eps; @@ -682,12 +682,12 @@ function swe_cp(SwE) % matrix design either in a visit category or within-subject (BG - 27/05/2016) %------------------------------------------------------------------ for g = 1:nGr_dof % first look data for each separable matrix design - if sum(iGr_dof'==g) > 1 % do not look for cases where the separable matrix design is only one row (BG - 05/08/2016) + if sum(iGr_dof'==g) > 1 % do not look for cases where the separable matrix design is only one row (BG - 05/08/2016) cmask(cmask) = any(abs(diff(Y(iGr_dof'==g, cmask),1)) > eps, 1); % mask constant data within separable matrix design g (added by BG on 29/08/16) if isfield(SwE.type,'modified') % added by BG on 29/08/16 for g2 = 1:nGr % then look data for each "homogeneous" group % check if the data is contant over subject for each visit category - for k = 1:nVis_g(g2) + for k = 1:nVis_g(g2) if sum(iGr_dof'==g & iGr == uGr(g2) & iVis == uVis_g{g2}(k)) > 1 % do not look for cases when the data is only one row (BG - 05/08/2016) cmask(cmask) = any(abs(diff(Y(iGr_dof'==g & iGr == uGr(g2) & iVis == uVis_g{g2}(k), cmask),1)) > eps, 1); for kk = k:nVis_g(g2) @@ -731,7 +731,7 @@ function swe_cp(SwE) end clear Y %-Clear to save memory - %-Estimation of the data variance-covariance components (modified SwE) + %-Estimation of the data variance-covariance components (modified SwE) %-SwE estimation (classic version) %-------------------------------------------------------------- c = zeros(numel(chunk),1); @@ -741,7 +741,7 @@ function swe_cp(SwE) for i = Ind_Cov_vis_diag Cov_vis(i,:) = mean(res(Flagk(i,:),:).^2, 1); end - % Check if some voxels have variance < eps and mask them + % Check if some voxels have variance < eps and mask them tmp = ~any(Cov_vis(Ind_Cov_vis_diag,:) < eps); % modified by BG on 29/08/16 if any(~tmp) beta = beta(:,tmp); @@ -778,9 +778,9 @@ function swe_cp(SwE) end end end - + % compute the beta covariance matrice(s) - switch dof_type + switch dof_type case 1 Cov_beta = zeros(nCov_beta, CrS); it = 0; @@ -823,19 +823,19 @@ function swe_cp(SwE) %---------------------------------------------------------------------- mask(chunk) = cmask; VM = swe_data_write(VM, cmask', chunk); - + %-Write beta files %---------------------------------------------------------------------- for iBeta=1:nBeta if CrS c(cmask) = beta(iBeta,:); end - Vbeta(iBeta) = swe_data_write(Vbeta(iBeta), c, chunk); + Vbeta(iBeta) = swe_data_write(Vbeta(iBeta), c, chunk); end %-Write CovVis files if needed %---------------------------------------------------------------------- - if isfield(SwE.type,'modified') + if isfield(SwE.type,'modified') for iCov_vis=1:nCov_vis if CrS c(cmask) = Cov_vis(iCov_vis,:); @@ -852,29 +852,29 @@ function swe_cp(SwE) end Vcov_beta(iCov_beta) = swe_data_write(Vcov_beta(iCov_beta), c, chunk); end - + %-Report progress %====================================================================== fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...done'); %-# swe_progress_bar('Set',i); end % iChunk=1:nbchunks - swe_progress_bar('Clear'); - + swe_progress_bar('Clear'); + %========================================================================== % - P O S T E S T I M A T I O N C L E A N U P %========================================================================== - + S = nnz(mask); if S == 0 error('Please check your data: There are no inmask voxels.'); end - + %-Compute coordinates of voxels within mask %-------------------------------------------------------------------------- [x,y,z] = ind2sub(DIM,find(mask)); XYZ = [x y z]'; - + else % matrix input % check how the data image treat 0 (as NaN or not) YNaNrep = 0; @@ -891,9 +891,9 @@ function swe_cp(SwE) elseif size(Y, 1) ~= SwE.nscan error('The input data does not have %i rows and thus is not compatible with the other specified variables. Please revised the model specification.', SwE.nscan) end - + nVox = size(Y, 2); - + %-Produce the mask cmask = true(1, nVox); %-Use the explicit mask if specified @@ -912,12 +912,12 @@ function swe_cp(SwE) % matrix design either in a visit category or within-subject (BG - 27/05/2016) %------------------------------------------------------------------ for g = 1:nGr_dof % first look data for each separable matrix design - if sum(iGr_dof'==g) > 1 % do not look for cases where the separable matrix design is only one row (BG - 05/08/2016) + if sum(iGr_dof'==g) > 1 % do not look for cases where the separable matrix design is only one row (BG - 05/08/2016) cmask(cmask) = any(abs(diff(Y(iGr_dof'==g,cmask),1)) > eps, 1); % mask constant data within separable matrix design g (added by BG on 29/08/16) if isfield(SwE.type,'modified') % added by BG on 29/08/16 for g2 = 1:nGr % then look data for each "homogeneous" group % check if the data is contant over subject for each visit category - for k = 1:nVis_g(g2) + for k = 1:nVis_g(g2) if sum(iGr_dof'==g & iGr == uGr(g2) & iVis == uVis_g{g2}(k)) > 1 % do not look for cases when the data is only one row (BG - 05/08/2016) cmask(cmask) = any(abs(diff(Y(iGr_dof'==g & iGr == uGr(g2) & iVis == uVis_g{g2}(k) ,cmask),1)) > eps, 1); for kk = k:nVis_g(g2) @@ -1032,12 +1032,12 @@ function swe_cp(SwE) if isfield(SwE.type,'modified') fprintf('\n'); %-# disp('Working on the SwE computation...'); - switch dof_type + switch dof_type case 1 crCov_beta = zeros(nCov_beta,CrS); % initialize SwE for the plane crCov_beta_i = zeros(nGr,nCov_beta,CrS); for g = 1:nGr - crCov_beta_i(g,:,:) = weight(:,iGr_Cov_vis_g==g) * crCov_vis(iGr_Cov_vis_g==g,:); + crCov_beta_i(g,:,:) = weight(:,iGr_Cov_vis_g==g) * crCov_vis(iGr_Cov_vis_g==g,:); Cov_beta = Cov_beta + crCov_beta_i(g,:,:); end case {0 2 3} @@ -1046,12 +1046,12 @@ function swe_cp(SwE) end fprintf('\n'); %-# end - - %-Save betas etc. + + %-Save betas etc. %---------------------------------------------------------- fprintf('%s%30s',repmat(sprintf('\b'),1,30),'...saving results'); %-# - mask = cmask; + mask = cmask; save(sprintf('swe_%s_mask%s',file_data_type,file_ext), 'mask'); clear mask @@ -1082,10 +1082,10 @@ function swe_cp(SwE) save(sprintf('swe_%s_cov_g_bb%s',file_data_type,file_ext), 'cov_beta_g'); clear cov_beta_g crCov_beta_i end - fprintf('%s%30s',repmat(sprintf('\b'),1,30),'...done'); - + fprintf('%s%30s',repmat(sprintf('\b'),1,30),'...done'); + fprintf('%s%30s',repmat(sprintf('\b'),1,30),'...some clean up'); %-# - + XYZ = []; M = []; DIM = []; @@ -1112,7 +1112,7 @@ function swe_cp(SwE) SwE.Vbeta = Vbeta; %-Filehandle - Beta SwE.Vcov_beta = Vcov_beta; %-Filehandle - Beta covariance if isfield(SwE.type,'modified') - SwE.Vcov_vis = Vcov_vis; %-Filehandle - Visit covariance + SwE.Vcov_vis = Vcov_vis; %-Filehandle - Visit covariance end if dof_type == 1 SwE.Vcov_beta_g = Vcov_beta_g; %-Filehandle - Beta covariance contributions @@ -1134,7 +1134,7 @@ function swe_cp(SwE) SwE.Subj.nSubj = nSubj; if isfield(SwE.type,'modified') - + if dof_type >1 SwE.Vis.weight = weight; SwE.Vis.iGr_Cov_vis_g = iGr_Cov_vis_g; @@ -1146,7 +1146,7 @@ function swe_cp(SwE) SwE.Vis.nCov_vis = nCov_vis; SwE.Vis.max_nVis_g = max_nVis_g; SwE.Vis.min_nVis_g = min_nVis_g; - + SwE.Gr.uGr = uGr; SwE.Gr.nGr = nGr; SwE.Gr.nSubj_g = nSubj_g; @@ -1157,20 +1157,20 @@ function swe_cp(SwE) end end -SwE.dof.uGr_dof = uGr_dof; +SwE.dof.uGr_dof = uGr_dof; SwE.dof.nGr_dof = nGr_dof; -SwE.dof.iGr_dof = iGr_dof; +SwE.dof.iGr_dof = iGr_dof; SwE.dof.iBeta_dof = iBeta_dof; SwE.dof.pB_dof = pB_dof; SwE.dof.nSubj_dof = nSubj_dof; SwE.dof.edof_Subj = edof_Subj; SwE.dof.dof_type = dof_type; -if dof_type == 1 +if dof_type == 1 SwE.dof.edof_Gr = edof_Gr; elseif dof_type == 0 SwE.dof.dof_cov = dof_cov; else - SwE.dof.dofMat = dofMat; + SwE.dof.dofMat = dofMat; end % save the version number of the toolbox @@ -1193,4 +1193,4 @@ function swe_cp(SwE) %========================================================================== %spm('FigName','Stats: done',Finter); spm('Pointer','Arrow') fprintf('%-40s: %30s\n','Completed',spm('time')) %-# -fprintf('...use the results section for assessment\n\n') +fprintf('...use the results section for assessment\n\n') diff --git a/swe_cp_WB.m b/swe_cp_WB.m index 6f99393..ddf50d7 100644 --- a/swe_cp_WB.m +++ b/swe_cp_WB.m @@ -4,26 +4,26 @@ function swe_cp_WB(SwE) % For a non-parametric SwE analysis with either NIfTI, GIfTI, CIfTI or '.mat' input, the % following maps are computed: % -% - swe_{unit}_mask: +% - swe_{unit}_mask: % The mask image for the analysis. % -% - swe_{unit}_{T|F}stat_c{c#}: +% - swe_{unit}_{T|F}stat_c{c#}: % Voxelwise parametric statistic map (T or F) for contrast {c#}. % -% - swe_{unit}_{zT|xF}stat_c{c#}: -% Voxelwise parametric equivalent statistic map (Z or Chi Squared) +% - swe_{unit}_{zT|xF}stat_c{c#}: +% Voxelwise parametric equivalent statistic map (Z or Chi Squared) % for contrast {c#}. % -% - swe_{unit}_{zT|xF}stat-WB_c{c#}: -% Voxelwise non-parametric equivalent statistic map (Z or Chi +% - swe_{unit}_{zT|xF}stat-WB_c{c#}: +% Voxelwise non-parametric equivalent statistic map (Z or Chi % Squared) for contrast {c#}. % % - swe_{unit}_{T|F}stat_lp-WB_c{c#}: -% Log10 map of the voxelwise uncorrected P values for contrast +% Log10 map of the voxelwise uncorrected P values for contrast % {c#}. % % - swe_{unit}_{T|F}stat_lpFWE-WB_c{c#}: -% Log10 map of the voxelwise bootstrap-calculated FWE P values for +% Log10 map of the voxelwise bootstrap-calculated FWE P values for % contrast {c#}. % % - swe_{unit}_{T|F}stat_lpFDR-WB_c{c#}: @@ -31,7 +31,7 @@ function swe_cp_WB(SwE) % contrast {c#}. % % - swe_clustere_{T|F}stat_lpFWE-WB_c{c#}: -% Log10 map of the clusterwise bootstrap-calculated FWE P values +% Log10 map of the clusterwise bootstrap-calculated FWE P values % for contrast {c#}. % % - swe_tfce_c{c#}: @@ -40,18 +40,18 @@ function swe_cp_WB(SwE) % - swe_tfce_lp-WB_c{c#}: % Log10 map of the TFCE bootstrap-calculated P values for contrast % {c#}. -% +% % - swe_tfce_lpFWE-WB_c{c#}: % Log10 map of the TFCE bootstrap-calculated FWE P values for % contrast {c#}. % -% The field {unit} used above represents the unit in space in which the statistic is calculated. +% The field {unit} used above represents the unit in space in which the statistic is calculated. % It can be the following strings: % - 'vox' for NifTI files, % - 'dpx' for GIfTI and CIfTI files, % - 'dat' for .mat files. -% Currently (30/08/2018), the only contrasts computed are activation -% (contrast #1) and deactivation (contrast #2) for the contrast vector the +% Currently (30/08/2018), the only contrasts computed are activation +% (contrast #1) and deactivation (contrast #2) for the contrast vector the % user input during the batch entry. % ========================================================================= % FORMAT swe_cp_WB(SwE) @@ -97,7 +97,7 @@ function swe_cp_WB(SwE) isOctave = exist('OCTAVE_VERSION','builtin'); if isCifti - metadata = {'ciftiTemplate', SwE.xY.P{1}}; + metadata = {'ciftiTemplate', SwE.xY.P{1}}; file_data_type = 'dpx'; dataType = swe_DataType('Cifti'); dataTypeSpecificInformation = SwE.cifti; @@ -189,13 +189,13 @@ function swe_cp_WB(SwE) %-Prevent unnecessary octave warning %-------------------------------------------------------------------------- if isOctave - warning ('off', 'histc: empty EDGES specified\n'); + warning ('off', 'histc: empty EDGES specified\n'); end %-Delete files from previous analyses %-------------------------------------------------------------------------- if exist(fullfile(SwE.swd,sprintf('swe_v%s_mask%s',file_data_type,file_ext)),'file') == 2 - + str = {'Current directory contains SwE estimation files:',... 'pwd = ',SwE.swd,... 'Existing results will be overwritten!'}; @@ -237,7 +237,7 @@ function swe_cp_WB(SwE) end % Tolerance for comparing real numbers -tol = 1e-8; +tol = 1e-8; %========================================================================== % - A N A L Y S I S P R E L I M I N A R I E S @@ -341,8 +341,8 @@ function swe_cp_WB(SwE) nSubj_dof(i) = length(unique(iSubj(iGr_dof==uGr_dof(i)))); end -pB_dof = zeros(1,nGr_dof); -for i=1:nBeta +pB_dof = zeros(1,nGr_dof); +for i=1:nBeta tmp=1; for ii=1:nSubj if length(unique(xX.X(iSubj==uSubj(ii)&iGr_dof'==iBeta_dof(i),i)))>1 @@ -366,7 +366,7 @@ function swe_cp_WB(SwE) if isfield(SwE.type,'modified') dof_type = SwE.type.modified.dof_mo; else - dof_type = SwE.type.classic.dof_cl; + dof_type = SwE.type.classic.dof_cl; end if dof_type == 0 % so naive estimation is used @@ -391,9 +391,9 @@ function swe_cp_WB(SwE) dof_cov = zeros(1,nBeta); for i = 1:nBeta dof_cov(i) = nSubj_dof(iBeta_dof(i)) - ... - pB_dof(iBeta_dof(i)); + pB_dof(iBeta_dof(i)); end - + % This variable should be left empty for Niave estimation. dofMat = []; xX.erdf_niave = edf; @@ -409,46 +409,46 @@ function swe_cp_WB(SwE) end iVis = SwE.Vis.iVis; iGr = SwE.Gr.iGr; - uGr = unique(iGr); + uGr = unique(iGr); nGr = length(uGr); SwE.Gr.uGr = uGr; SwE.Gr.nGr = nGr; - + % info specific for each group uVis_g = cell(1,nGr); % unique visits for each group nVis_g = zeros(1,nGr); % number of visits for each group uSubj_g = cell(1,nGr); % unique visits for each group nSubj_g = zeros(1,nGr); % number of visits for each group for g = 1:nGr - uVis_g{g} = unique(iVis(iGr==uGr(g))); + uVis_g{g} = unique(iVis(iGr==uGr(g))); nVis_g(g) = length(uVis_g{g}); iSubj_g = iSubj(iGr==uGr(g)); % Subject number for each subject in group for each visit uSubj_g{g} = unique(iSubj_g); % Unique subject numbers of subjects in group nSubj_g(g) = length(uSubj_g{g}); uSubj_g_tmp = uSubj_g{g}; - + for k = 1:nSubj_g(g) - + % The number of visits for subject uSubj_g(k) vis_g_subj(k) = sum(iSubj_g==uSubj_g_tmp(k)); - + end max_nVis_g(g) = max(vis_g_subj); min_nVis_g(g) = min(vis_g_subj); - + clear vis_g_subj - + end nCov_vis_g = nVis_g.*(nVis_g+1)/2; % number of covariance elements to be estimated for each group nCov_vis = sum(nCov_vis_g); % total number of covariance elements to be estimated - + % Save Vis variables. SwE.Vis.uVis_g = uVis_g; SwE.Vis.nVis_g = nVis_g; SwE.Vis.max_nVis_g = max_nVis_g; SwE.Vis.min_nVis_g = min_nVis_g; - + % Flags matrices indicating which residuals have to be used for each covariance element Flagk = false(nCov_vis,nScan); % Flag indicating scans corresponding to visit k for each covariance element Flagkk = false(nCov_vis,nScan); % Flag indicating scans corresponding to visit kk for each covariance element @@ -479,10 +479,10 @@ function swe_cp_WB(SwE) end end end - + % Record igr_Cov_vis_g. SwE.WB.iGr_Cov_vis_g = iGr_Cov_vis_g; - + % weights for the vectorised SwE weight = NaN(nCov_beta,nCov_vis); @@ -500,22 +500,22 @@ function swe_cp_WB(SwE) end end % Weight giving only the contrasted SwE (WB) - weightR = pinv(swe_duplication_matrix(nSizeCon)) * kron(conWB,conWB) * swe_duplication_matrix(nBeta) * weight; % used to compute the R SwE R' + weightR = pinv(swe_duplication_matrix(nSizeCon)) * kron(conWB,conWB) * swe_duplication_matrix(nBeta) * weight; % used to compute the R SwE R' Wg = cell(nGr,1); Wg_testII = cell(nGr,1); Wg_testIII = cell(nGr,1); - tmp = eye(nSizeCon^2); + tmp = eye(nSizeCon^2); for g = 1:nGr Wg{g} = kron(weightR(:,iGr_Cov_vis_g==g),weightR(:,iGr_Cov_vis_g==g)) * swe_duplication_matrix(nCov_vis_g(g)); Wg_testII{g} = sum(kron(swe_duplication_matrix(nSizeCon),swe_duplication_matrix(nSizeCon)), 1) * Wg{g}; Wg_testIII{g} = tmp(:)' * (kron(swe_duplication_matrix(nSizeCon),swe_duplication_matrix(nSizeCon))) * Wg{g}; end - + SwE.WB.Wg{1} = Wg; SwE.WB.Wg{2} = Wg_testII; SwE.WB.Wg{3} = Wg_testIII; - + %-compute the effective dof from each homogeneous group if dof_type switch dof_type case 1 @@ -579,7 +579,7 @@ function swe_cp_WB(SwE) Ind_Cov_vis_classic = NaN(1,nCov_vis); Indexk = NaN(1,nCov_vis); Indexkk = NaN(1,nCov_vis); - + it=0; for j = 1:nBeta for jj = j:nBeta @@ -608,7 +608,7 @@ function swe_cp_WB(SwE) if dof_type == 1 edof_Gr = edof_Subj; end - weightR = pinv(swe_duplication_matrix(nSizeCon)) * kron(conWB,conWB) * swe_duplication_matrix(nBeta) * weight; % used to compute the R SwE R' + weightR = pinv(swe_duplication_matrix(nSizeCon)) * kron(conWB,conWB) * swe_duplication_matrix(nBeta) * weight; % used to compute the R SwE R' end %-If xM is not a structure then assume it's a vector of thresholds @@ -633,7 +633,7 @@ function swe_cp_WB(SwE) %========================================================================== VY = SwE.xY.VY; spm_check_orientations(VY); - + % check files exists and try pwd %-------------------------------------------------------------------------- for i = 1:numel(VY) @@ -642,19 +642,19 @@ function swe_cp_WB(SwE) VY(i).fname = [n,e]; end end - + M = VY(1).mat; DIM = VY(1).dim; YNaNrep = VY(1).dt(2); - + fprintf('%-40s: %30s','Output images','...initialising'); %-# - + %-Initialise new mask name: current mask & conditions on voxels %---------------------------------------------------------------------- VM = swe_data_hdr_write(sprintf('swe_%s_mask%s', file_data_type, file_ext), DIM, M,... 'swe_cp:resultant analysis mask', metadata, 'uint8'); - + %-Initialise beta image files %---------------------------------------------------------------------- @@ -664,7 +664,7 @@ function swe_cp_WB(SwE) sprintf('swe_cp:beta (%02d) - %s',i,xX.name{i}),... metadata); end - + %-Initialise original parametric score image, T or F %---------------------------------------------------------------------- if WB.stat=='T' @@ -672,63 +672,63 @@ function swe_cp_WB(SwE) else % F stat eSTAT='x'; end - + Vscore = swe_data_hdr_write(sprintf('swe_%s_%cstat_c%02d%s', file_data_type, WB.stat, 1, file_ext), DIM, M,... sprintf('Original parametric %c statistic data.', WB.stat), metadata); - + %-Initialise parametric TFCE score image, if TFCE has been selected. - %---------------------------------------------------------------------- + %---------------------------------------------------------------------- if TFCE Vscore_tfce = swe_data_hdr_write(sprintf('swe_tfce_c%02d%s', 1, file_ext), DIM, M,... - 'Original parametric TFCE statistic data.', metadata); + 'Original parametric TFCE statistic data.', metadata); if WB.stat=='T' Vscore_tfce_neg = swe_data_hdr_write(sprintf('swe_tfce_c%02d%s', 2, file_ext), DIM, M,... - 'Original parametric TFCE statistic data for a negative contrast.', metadata); + 'Original parametric TFCE statistic data for a negative contrast.', metadata); end end - + %-Initialise original parametric edf image %---------------------------------------------------------------------- - + Vedf = swe_data_hdr_write(sprintf('swe_%s_edf_c%02d%s', file_data_type, 1, file_ext), DIM, M,... sprintf('Original parametric %c edf data.', WB.stat), metadata); - + %-Initialise parametric P-Value image %---------------------------------------------------------------------- - + VlP = swe_data_hdr_write(sprintf('swe_%s_%cstat_lp_c%02d%s', file_data_type, WB.stat, 1, file_ext), DIM, M,... 'Original parametric -log10(P) value data (positive).', metadata); - + if WB.stat=='T' VlP_Neg = swe_data_hdr_write(sprintf('swe_%s_%cstat_lp_c%02d%s', file_data_type, WB.stat, 2, file_ext), DIM, M,... 'Original parametric -log10(P) value data (negative).', metadata); end - + %-Initialise converted parametric score image %---------------------------------------------------------------------- VcScore = swe_data_hdr_write(sprintf('swe_%s_%c%cstat_c%02d%s', file_data_type, eSTAT, WB.stat, 1, file_ext), DIM, M,... - sprintf('Parametric %c statistic data derived from %c-Statistic data.', eSTAT, WB.stat), metadata); + sprintf('Parametric %c statistic data derived from %c-Statistic data.', eSTAT, WB.stat), metadata); if WB.stat=='T' VcScore_neg = swe_data_hdr_write(sprintf('swe_%s_%c%cstat_c%02d%s', file_data_type, eSTAT, WB.stat, 2, file_ext), DIM, M,... sprintf('Parametric %c statistic data derived from %c-Statistic data.', eSTAT, WB.stat), metadata); end - + %-Initialise residual images for the resampling %---------------------------------------------------------------------- - + for i = 1:nScan descrip = sprintf('adjusted restricted residuals (%04d)', i); VResWB(i) = swe_data_hdr_write(sprintf('swe_%s_resid_y%04d%s', file_data_type, i, file_ext), DIM, M, descrip, metadata); end - + %-Initialise fitted data images for the resampling %---------------------------------------------------------------------- - + for i = 1:nScan descrip = sprintf('restricted fitted data (%04d)', i); VYWB(i) = swe_data_hdr_write(sprintf('swe_%s_fit_y%04d%s',file_data_type,i,file_ext), DIM, M, descrip, metadata); end - + %-Initialise result images %---------------------------------------------------------------------- VlP_wb_pos = swe_data_hdr_write(sprintf('swe_%s_%cstat_lp-WB_c%02d%s', file_data_type, WB.stat, 1, file_ext), DIM, M,... @@ -736,60 +736,60 @@ function swe_cp_WB(SwE) VlP_wb_FWE_pos = swe_data_hdr_write(sprintf('swe_%s_%cstat_lpFWE-WB_c%02d%s', file_data_type, WB.stat, 1, file_ext), DIM, M,... 'Non-parametric voxelwise FWE -log10(P) value data (positive).', metadata); - + VlP_wb_FDR_pos = swe_data_hdr_write(sprintf('swe_%s_%cstat_lpFDR-WB_c%02d%s', file_data_type, WB.stat, 1, file_ext), DIM, M,... 'Non-parametric voxelwise FDR -log10(P) value data (positive).', metadata); if WB.stat=='T' VlP_wb_neg = swe_data_hdr_write(sprintf('swe_%s_%cstat_lp-WB_c%02d%s', file_data_type, WB.stat, 2, file_ext), DIM, M,... 'Non-parametric voxelwise -log10(P) value data (negative).', metadata); - + VlP_wb_FWE_neg = swe_data_hdr_write(sprintf('swe_%s_%cstat_lpFWE-WB_c%02d%s', file_data_type, WB.stat, 2, file_ext), DIM, M,... 'Non-parametric voxelwise FWE -log10(P) value data (negative).', metadata); - + VlP_wb_FDR_neg = swe_data_hdr_write(sprintf('swe_%s_%cstat_lpFDR-WB_c%02d%s', file_data_type, WB.stat, 2, file_ext), DIM, M,... 'Non-parametric voxelwise FDR -log10(P) value data (negative).', metadata); end - + %-Initialise parametric TFCE results images, if TFCE has been selected. - %---------------------------------------------------------------------- + %---------------------------------------------------------------------- if TFCE VlP_tfce_pos = swe_data_hdr_write(sprintf('swe_tfce_lp-WB_c%02d%s', 1, file_ext), DIM, M,... - 'Non-parametric TFCE -log10(P) value data (positive).', metadata); + 'Non-parametric TFCE -log10(P) value data (positive).', metadata); VlP_tfce_FWE_pos = swe_data_hdr_write(sprintf('swe_tfce_lpFWE-WB_c%02d%s', 1, file_ext), DIM, M,... - 'Non-parametric TFCE FWE -log10(P) value data (positive).', metadata); + 'Non-parametric TFCE FWE -log10(P) value data (positive).', metadata); if WB.stat=='T' VlP_tfce_neg = swe_data_hdr_write(sprintf('swe_tfce_lp-WB_c%02d%s', 2, file_ext), DIM, M,... - 'Non-parametric TFCE -log10(P) value data (negative).', metadata); + 'Non-parametric TFCE -log10(P) value data (negative).', metadata); VlP_tfce_FWE_neg = swe_data_hdr_write(sprintf('swe_tfce_lpFWE-WB_c%02d%s', 2, file_ext), DIM, M,... 'Non-parametric TFCE FWE -log10(P) value data (negative).', metadata); end end - + % Converted score for WB. VcScore_wb_pos = swe_data_hdr_write(sprintf('swe_%s_%c%cstat-WB_c%02d%s', file_data_type, eSTAT, WB.stat, 1, file_ext), DIM, M,... sprintf('Non-parametric %c statistic data derived from %c-Statistic data.', eSTAT, WB.stat), metadata); - + if WB.clusterWise == 1 - + % We also need cluster p value maps here. V_clustere_pos = swe_data_hdr_write(sprintf('swe_clustere_%cstat_c%02d%s', WB.stat, 1, file_ext), DIM, M,... sprintf('Cluster extent (positive, CFT %g).',... SwE.WB.clusterInfo.primaryThreshold), metadata); - + if WB.stat=='T' V_clustere_neg = swe_data_hdr_write(sprintf('swe_clustere_%cstat_c%02d%s', WB.stat, 2, file_ext), DIM, M,... sprintf('Cluster extent (negative, CFT %g).',... - SwE.WB.clusterInfo.primaryThreshold), metadata); + SwE.WB.clusterInfo.primaryThreshold), metadata); end end - + fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...initialised'); %-# %========================================================================== % - F I T M O D E L & W R I T E P A R A M E T E R I M A G E S %========================================================================== - + %-Get explicit mask(s) %========================================================================== mask = true(DIM); @@ -829,7 +829,7 @@ function swe_cp_WB(SwE) chunksize = floor(spm_get_defaults('stats.maxmem') / 8 / nScan); nbchunks = ceil(prod(DIM) / chunksize); chunks = min(cumsum([1 repmat(chunksize,1,nbchunks)]),prod(DIM)+1); - + % activated voxels for cluster-wise inference if (WB.clusterWise == 1) activatedVoxels = false(0); @@ -847,7 +847,7 @@ function swe_cp_WB(SwE) if (WB.stat == 'T') minScore = nan(1, WB.nB + 1); end - + %-Cycle over bunches blocks within planes to avoid memory problems %========================================================================== swe_progress_bar('Init',nbchunks,'Parameter estimation','Chunks'); @@ -857,51 +857,51 @@ function swe_cp_WB(SwE) %-Report progress %====================================================================== - if iChunk > 1, fprintf(repmat(sprintf('\b'),1,72)); end %-# + if iChunk > 1, fprintf(repmat(sprintf('\b'),1,72)); end %-# fprintf('%-40s: %30s', sprintf('Original statistics: Chunk %3d/%-3d',iChunk,nbchunks),... '...processing'); - + %-Get the data in mask, compute threshold & implicit masks %------------------------------------------------------------------ Y = zeros(nScan, numel(chunk)); cmask = mask(chunk); - + if size(cmask, 2) == 1 cmask = cmask'; end for iScan=1:nScan if ~any(cmask), break, end %-Break if empty mask - + Y(iScan, cmask) = swe_data_read(VY(iScan), chunk(cmask));%-Read chunk of data - + cmask(cmask) = Y(iScan, cmask) > xM.TH(iScan); %-Threshold (& NaN) mask if xM.I && ~YNaNrep && xM.TH(iScan) < 0 %-Use implicit mask cmask(cmask) = abs(Y(iScan, cmask)) > eps; end end - cmask(cmask) = any(diff(Y(:,cmask),1)); - + cmask(cmask) = any(diff(Y(:,cmask),1)); + %-Mask out voxels where data is constant in at least one separable % matrix design either in a visit category or within-subject (BG - 27/05/2016) %------------------------------------------------------------------ [cmask, Y, CrS] = swe_mask_seperable(SwE, cmask, Y, iGr_dof); - + %================================================================== %-Proceed with General Linear Model (if there are voxels) %================================================================== if CrS - + %-General linear model: Ordinary least squares estimation - %-------------------------------------------------------------- + %-------------------------------------------------------------- beta = pX*Y; %-Parameter estimates - + % restricted fitted data [resWB, YWB] = swe_fit(SwE, Y, tmpR2, corrWB, beta, SwE.WB.SS); if WB.RSwE == 1 res = swe_fit(SwE, Y, tmpR2, corr, beta, SwE.SS); - else + else res = swe_fit(SwE, Y, xX.X, corr, beta, SwE.SS); end @@ -914,7 +914,7 @@ function swe_cp_WB(SwE) for i = Ind_Cov_vis_diag Cov_vis(i,:) = mean(res(Flagk(i,:),:).^2, 1); end - + % Check if some voxels have variance < eps and mask them tmp = ~any(Cov_vis(Ind_Cov_vis_diag,:) < eps); % modified by BG on 29/08/16 if any(~tmp) @@ -943,7 +943,7 @@ function swe_cp_WB(SwE) %need to check if the eigenvalues of Cov_vis matrices are >=0 if dof_type == 1 tmpSum = zeros(1,CrS); - end + end for g = 1:nGr for iVox = 1:CrS tmp = zeros(nVis_g(g)); @@ -977,13 +977,13 @@ function swe_cp_WB(SwE) cCovBc = 0; if dof_type == 1 tmpSum = zeros(1,CrS); - end + end for i = 1:nSubj Cov_beta_i_tmp = weightR(:,Ind_Cov_vis_classic==i) *... (res(Indexk(Ind_Cov_vis_classic==i),:) .* res(Indexkk(Ind_Cov_vis_classic==i),:)); cCovBc = cCovBc + Cov_beta_i_tmp; if dof_type == 1 - if nSizeCon == 1 + if nSizeCon == 1 tmpSum = tmpSum + Cov_beta_i_tmp.^2/edof_Gr(i); else for iVox = 1:CrS @@ -1019,9 +1019,9 @@ function swe_cp_WB(SwE) % compute the score if (SwE.WB.stat == 'T') - + score = (conWB * beta) ./ sqrt(cCovBc); - + % hypothesis test, using clusterwise threshold if available. if (SwE.WB.clusterWise == 1) if dof_type == 1 @@ -1048,7 +1048,7 @@ function swe_cp_WB(SwE) minScore(1) = min(minScore(1), min(hyptest.positive.conScore)); - + else % need to loop at every voxel cBeta = conWB * beta; @@ -1060,7 +1060,7 @@ function swe_cp_WB(SwE) score(iVox) = cBeta(:,iVox)' / cCovBc_vox * cBeta(:,iVox); end score = score / rankCon; - + % hypothesis test, using clusterwise threshold if available. if (SwE.WB.clusterWise == 1) if dof_type == 1 @@ -1081,9 +1081,9 @@ function swe_cp_WB(SwE) edf = hyptest.positive.edf; end end - + maxScore(1) = max(maxScore(1), max(hyptest.positive.conScore)); - + end % (CrS) %-Write output files @@ -1094,14 +1094,14 @@ function swe_cp_WB(SwE) %---------------------------------------------------------------------- mask(chunk) = cmask; VM = swe_data_write(VM, cmask', chunk); - + %-Write beta files %---------------------------------------------------------------------- for iBeta = 1:nBeta if CrS c(cmask) = beta(iBeta,:); end - Vbeta(iBeta) = swe_data_write(Vbeta(iBeta), c, chunk); + Vbeta(iBeta) = swe_data_write(Vbeta(iBeta), c, chunk); end %-Write WB fitted data images @@ -1112,7 +1112,7 @@ function swe_cp_WB(SwE) end VYWB(iScan) = swe_data_write(VYWB(iScan), c, chunk); end - + %-Write WB residuals %------------------------------------------------------------------ for iScan = 1:nScan @@ -1128,28 +1128,28 @@ function swe_cp_WB(SwE) c(cmask) = score; end Vscore = swe_data_write(Vscore, c, chunk); - + %-Write parametric edf image of the original data %------------------------------------------------------------------ if CrS c(cmask) = hyptest.positive.edf; end Vedf = swe_data_write(Vedf, c, chunk); - + %-Write parametric p-value image %------------------------------------------------------------------ if CrS c(cmask) = -log10(hyptest.positive.p); end VlP = swe_data_write(VlP, c, chunk); - + if WB.stat=='T' if CrS c(cmask) = -log10(hyptest.negative.p); end VlP_Neg = swe_data_write(VlP_Neg, c, chunk); end - + %-Write converted parametric score image of the original data %------------------------------------------------------------------ if CrS @@ -1164,50 +1164,50 @@ function swe_cp_WB(SwE) %-Report progress %====================================================================== swe_progress_bar('Set',iChunk); - end % iChunk=1:nbchunks + end % iChunk=1:nbchunks %========================================================================== % - P O S T E S T I M A T I O N C L E A N U P %========================================================================== - + S = nnz(mask); if S == 0 error('Please check your data: There are no inmask voxels.'); end - + %-Compute coordinates of voxels within mask %-------------------------------------------------------------------------- [x,y,z] = ind2sub(DIM,find(mask)); XYZ = [x y z]'; - + if TFCE % Create parametric TFCE statistic images. if strcmp(WB.stat, 'T') - + % Read in T statistics to get negative and positive TFCE scores. par_tfce = swe_tfce_transform(swe_data_read(VcScore), H, E, C, dh); par_tfce_neg = swe_tfce_transform(-swe_data_read(VcScore), H, E, C, dh); else - + % Convert F statistics to Z scores. scorevol=-swe_invNcdf(10.^(-swe_data_read(VlP))); scorevol(isnan(scorevol))=0; - + % Convert to TFCE. par_tfce = swe_tfce_transform(scorevol, H, E, C, dh); - + clear scorevol end - + % Save parametric TFCE statistic images. swe_data_write(Vscore_tfce, par_tfce); if strcmp(WB.stat, 'T') swe_data_write(Vscore_tfce_neg, par_tfce_neg); end end - + clear beta res Cov_vis c %-Clear to save memory - + % compute the max cluster size if needed (so many ways this can be % done... Not sure this solution is the best) if (SwE.WB.clusterWise == 1) @@ -1220,7 +1220,7 @@ function swe_cp_WB(SwE) end originalClusterStatistics = swe_getClusterStatistics(dataType, LocActivatedVoxels, dataTypeSpecificInformation, giftiAreaFile); - + if originalClusterStatistics.nCluster == 0 warning('no clusters survived the cluster-forming thresholding of the original data for positive effects!') end @@ -1242,7 +1242,7 @@ function swe_cp_WB(SwE) end if (SwE.WB.stat == 'T') - + if dataType == swe_DataType('Gifti') LocActivatedVoxelsNeg = false(nDataElements, 1); LocActivatedVoxelsNeg(mask) = activatedVoxelsNeg; @@ -1280,21 +1280,21 @@ function swe_cp_WB(SwE) else % ".mat" format % check how the data image treat 0 (as NaN or not) YNaNrep = 0; - + fprintf('%-40s: %30s','Output images','...initialising'); %-# - + fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...initialised'); %-# %========================================================================== % - F I T M O D E L & W R I T E P A R A M E T E R I M A G E S %========================================================================== - - + + %-Cycle over bunches blocks within planes to avoid memory problems %========================================================================== str = 'parameter estimation'; swe_progress_bar('Init',100,str,''); - + % activated voxels for cluster-wise inference if (WB.clusterWise == 1) activatedVoxels = false(0); @@ -1308,11 +1308,11 @@ function swe_cp_WB(SwE) if (WB.stat == 'T') minScore = nan(1, WB.nB + 1); end - + %-Get data & construct analysis mask %================================================================= fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...read & mask data') - + %-Get the data in mask, compute threshold & implicit masks %------------------------------------------------------------------ Y = importdata(SwE.xY.P{1}); @@ -1322,9 +1322,9 @@ function swe_cp_WB(SwE) elseif size(Y, 1) ~= SwE.nscan error('The input data does not have %i rows and thus is not compatible with the other specified variables. Please revised the model specification.', SwE.nscan) end - + nDataElements = size(Y, 2); - + %-Produce the mask cmask = true(1, nDataElements); %-Use the explicit mask if specified @@ -1343,7 +1343,7 @@ function swe_cp_WB(SwE) % matrix design either in a visit category or within-subject (BG - 27/05/2016) %------------------------------------------------------------------ [cmask,Y,CrS] = swe_mask_seperable(SwE, cmask, Y, iGr_dof); - + if WB.clusterWise == 1 if isVolumeMat XYZ = XYZ(:,cmask); @@ -1353,21 +1353,21 @@ function swe_cp_WB(SwE) %-Proceed with General Linear Model (if there are voxels) %================================================================== if CrS - + %-General linear model: Ordinary least squares estimation %-------------------------------------------------------------- fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...estimation');%-# - + beta = pX*Y; %-Parameter estimates - + [resWB, YWB]=swe_fit(SwE, Y, tmpR2, corrWB, beta, SwE.WB.SS); - + if WB.RSwE == 1 res=swe_fit(SwE, Y, tmpR2, corr, beta, SwE.SS); - else + else res=swe_fit(SwE, Y, xX.X, corr, beta, SwE.SS); end - + clear Y %-Clear to save memory %-Estimation of the data variance-covariance components (modified SwE) %-SwE estimation (classic version) @@ -1380,7 +1380,7 @@ function swe_cp_WB(SwE) for i = Ind_Cov_vis_diag Cov_vis(i,:) = mean(res(Flagk(i,:),:).^2, 1); end - + % Check if some voxels have variance < eps and mask them tmp = ~any(Cov_vis(Ind_Cov_vis_diag,:) < eps); % modified by BG on 29/08/16 if any(~tmp) @@ -1442,7 +1442,7 @@ function swe_cp_WB(SwE) (res(Indexk(Ind_Cov_vis_classic==i),:) .* res(Indexkk(Ind_Cov_vis_classic==i),:)); cCovBc = cCovBc + Cov_beta_i_tmp; if dof_type == 1 - if nSizeCon == 1 + if nSizeCon == 1 tmpSum = tmpSum + Cov_beta_i_tmp.^2/edof_Gr(i); else for iVox = 1:CrS @@ -1477,16 +1477,16 @@ function swe_cp_WB(SwE) clear tmpSum % compute the score if (SwE.WB.stat == 'T') - + score = (conWB * beta) ./ sqrt(cCovBc); - + if (SwE.WB.clusterWise == 1) if dof_type == 1 hyptest=swe_hyptest(SwE, score, CrS, cCovBc, Cov_vis, edf, activatedVoxels, activatedVoxelsNeg); else hyptest=swe_hyptest(SwE, score, CrS, cCovBc, Cov_vis, dofMat, activatedVoxels, activatedVoxelsNeg); end - p = hyptest.positive.p; + p = hyptest.positive.p; negp = hyptest.negative.p; edf = hyptest.positive.edf; activatedVoxels = hyptest.positive.activatedVoxels; @@ -1502,13 +1502,13 @@ function swe_cp_WB(SwE) edf = hyptest.positive.edf; clear CovcCovBc cCovBc end - + minScore(1) = min(hyptest.positive.conScore); - + if TFCE %%%% TODO: T (MAKE SCORE) end - + else % need to loop at every voxel cBeta = conWB * beta; @@ -1539,30 +1539,30 @@ function swe_cp_WB(SwE) p = hyptest.positive.p; edf = hyptest.positive.edf; end - + if TFCE %%%% TODO: F (MAKE SCORE) end - + end maxScore(1) = max(hyptest.positive.conScore); - + end % (CrS) M = []; DIM = []; - S = CrS; + S = CrS; VM = sprintf('swe_%s_mask%s', file_data_type, file_ext); Vscore = sprintf('swe_%s_%cstat_c%02d%s', file_data_type, WB.stat, 1, file_ext); Vedf = sprintf('swe_%s_edf_c%02d%s', file_data_type, 1, file_ext); - mask = cmask; + mask = cmask; save(sprintf('swe_%s_mask%s', file_data_type, file_ext), 'mask'); - + edf(:,cmask) = edf; save(Vedf, 'edf'); clear Vedf - + if (SwE.WB.stat == 'T') T = zeros(1, nDataElements); T(:,cmask) = score; @@ -1574,52 +1574,52 @@ function swe_cp_WB(SwE) save(Vscore, 'F'); clear F end - + lP = zeros(1, nDataElements); lP(:,cmask) = -log10(hyptest.positive.p); save(sprintf('swe_%s_%cstat_lp_c%02d%s', file_data_type, WB.stat, 1, file_ext), 'lP'); clear lP - + if (SwE.WB.stat == 'T') - + lP_neg = zeros(1, nDataElements); lP_neg(:,cmask) = -log10(hyptest.negative.p); save(sprintf('swe_%s_%cstat_lp_c%02d%s', file_data_type, WB.stat, 2, file_ext), 'lP_neg'); clear lP_neg - + Z = zeros(1, nDataElements); Z(:,cmask) = hyptest.positive.conScore; save(sprintf('swe_%s_z%cstat_c%02d%s', file_data_type, WB.stat, 1, file_ext), 'Z'); clear Z - + else - + X = zeros(1, nDataElements); X(:,cmask) = hyptest.positive.conScore; save(sprintf('swe_%s_x%cstat_c%02d%s', file_data_type, WB.stat, 1, file_ext), 'X'); clear X - + end - + swe_progress_bar('Clear') clear res Cov_vis jj%-Clear to save memory - + % compute the max cluster size if needed (so many ways this can be % done... Not sure this solution is the best) if (SwE.WB.clusterWise == 1) - + if dataType == swe_DataType('VolumeMat') LocActivatedVoxels = XYZ(:,activatedVoxels); else %surface data LocActivatedVoxels = false(nDataElements, 1); LocActivatedVoxels(mask) = activatedVoxels; end - + originalActivatedVoxels = activatedVoxels; originalClusterStatistics = swe_getClusterStatistics(dataType, LocActivatedVoxels, dataTypeSpecificInformation, giftiAreaFile); - + if originalClusterStatistics.nCluster == 0 warning('no clusters survived the cluster-forming thresholding of the original data for positive effects!'); end @@ -1646,8 +1646,8 @@ function swe_cp_WB(SwE) maxClusterSizeNeg(1) = originalClusterStatisticsNeg.maxClusterSize; end - end -end + end +end %========================================================================== % - P O S T E S T I M A T I O N C L E A N U P %========================================================================== @@ -1751,13 +1751,13 @@ function swe_cp_WB(SwE) SwE.swd = pwd; if isfield(SwE.type,'modified') - + SwE.Vis.nCov_vis_g = nCov_vis_g; SwE.Vis.nCov_vis = nCov_vis; SwE.Gr.nSubj_g = nSubj_g; SwE.Gr.uSubj_g = uSubj_g; - + end %-Save analysis parameters in SwE.mat file @@ -1773,7 +1773,7 @@ function swe_cp_WB(SwE) fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...done'); %-# %========================================================================== -%- Produce bootstraps and maximum stats +%- Produce bootstraps and maximum stats %========================================================================== % check whether a resampling Matrix exists in SwE. If yes, use it. If not, produce it. @@ -1800,7 +1800,7 @@ function swe_cp_WB(SwE) originalScore = swe_data_read(VcScore, 'xyz', XYZ); % # blocks nbchunks = ceil(S / chunksize); - chunks = min(cumsum([1 repmat(chunksize, 1, nbchunks)]), S+1); + chunks = min(cumsum([1 repmat(chunksize, 1, nbchunks)]), S+1); end % variables for results uncP = ones(1, S); % one because of the original score @@ -1816,7 +1816,7 @@ function swe_cp_WB(SwE) if SwE.WB.stat == 'T' tfce_uncP_neg = zeros(DIM(1), DIM(2), DIM(3)); end - + % We also need to record the TFCE maximas for TFCE FWE (including the % original parametric max). if TFCE @@ -1827,7 +1827,7 @@ function swe_cp_WB(SwE) maxTFCEScore_neg(1) = max(par_tfce_neg(:)); end end - + end if WB.clusterWise == 1 @@ -1843,7 +1843,7 @@ function swe_cp_WB(SwE) tic str = sprintf('Parameter estimation\nBootstrap # %i', b); swe_progress_bar('Set','xlabel', str) - + % activated voxels for cluster-wise inference if (SwE.WB.clusterWise == 1) activatedVoxels = false(1,S); @@ -1852,35 +1852,35 @@ function swe_cp_WB(SwE) end end if ~isMat - + if TFCE - + % Instantiate volume for TFCE conversion. scorevol = zeros(DIM(1), DIM(2), DIM(3)); - + end - + for iChunk=1:nbchunks chunk = chunks(iChunk):chunks(iChunk+1)-1; sizeChunk = length(chunk); %-Print progress information in command window %------------------------------------------------------------------ - if iChunk > 1, fprintf(repmat(sprintf('\b'),1,72)); end %-# + if iChunk > 1, fprintf(repmat(sprintf('\b'),1,72)); end %-# fprintf('%-40s: %30s', sprintf('Bootstrap # %i: Chunk %3d/%-3d', b, iChunk, nbchunks),... '...processing'); - + Y_b = swe_data_read(VYWB, 'xyz', XYZ(:,chunk)) + ... swe_data_read(VResWB, 'xyz', XYZ(:,chunk)) .* repmat(resamplingMatrix(:,b),1,sizeChunk); - + beta = pX * Y_b; %-Parameter estimates if WB.RSwE == 0 res=swe_fit(SwE, Y_b, xX.X, corr, beta, SwE.SS); - else + else res=swe_fit(SwE, Y_b, tmpR2, corr, beta, SwE.SS); end - + clear Y_b - + %-Estimation of the data variance-covariance components (modified SwE) %-SwE estimation (classic version) %-------------------------------------------------------------- @@ -1939,7 +1939,7 @@ function swe_cp_WB(SwE) (res(Indexk(Ind_Cov_vis_classic==i),:) .* res(Indexkk(Ind_Cov_vis_classic==i),:)); cCovBc = cCovBc + Cov_beta_i_tmp; if dof_type == 1 - if nSizeCon == 1 + if nSizeCon == 1 tmpSum = tmpSum + Cov_beta_i_tmp.^2/edof_Gr(i); else for iVox = 1:sizeChunk @@ -1974,12 +1974,12 @@ function swe_cp_WB(SwE) clear tmpSum % compute the score if (SwE.WB.stat == 'T') - + score = (conWB * beta) ./ sqrt(cCovBc); clear beta - + else - + cBeta = conWB * beta; clear beta score = zeros(1, sizeChunk); @@ -1990,9 +1990,9 @@ function swe_cp_WB(SwE) score(iVox) = cBeta(:,iVox)' / cCovBc_vox * cBeta(:,iVox); end score = score / rankCon; - + end - + if dof_type == 1 hyptest = swe_hyptest(SwE, score, sizeChunk, cCovBc, Cov_vis, edf); else @@ -2012,13 +2012,13 @@ function swe_cp_WB(SwE) if (SwE.WB.stat == 'T') minScore(b+1) = min(minScore(b+1), min(hyptest.positive.conScore)); end - + % Calculate TFCE uncorrected p image. if TFCE % Current XYZ indices currXYZ = XYZ(1:3, chunk); - + % T test already converted to Z if strcmp(WB.stat, 'T') scorevol(sub2ind(DIM',currXYZ(1,:),currXYZ(2,:),currXYZ(3,:))) = hyptest.positive.conScore; @@ -2031,61 +2031,61 @@ function swe_cp_WB(SwE) % Save as scorevol scorevol(sub2ind(DIM',currXYZ(1,:),currXYZ(2,:),currXYZ(3,:))) = sv; end - + end - + end % (iChunk) - + if TFCE - if SwE.WB.stat == 'T' - + if SwE.WB.stat == 'T' + % Bootstrapped tfce vol. tfce = swe_tfce_transform(scorevol,H,E,C,dh); - tfce_neg = swe_tfce_transform(-scorevol,H,E,C,dh); - + tfce_neg = swe_tfce_transform(-scorevol,H,E,C,dh); + else - + % Bootstrapped tfce vol. tfce = swe_tfce_transform(scorevol,H,E,C,dh); - + end - + % Sum how many voxels are lower than the original parametric tfce. tfce_uncP = tfce_uncP + (par_tfce - tol < tfce); if SwE.WB.stat == 'T' tfce_uncP_neg = tfce_uncP_neg + (par_tfce_neg - tol < tfce_neg); end - + % Record maxima for TFCE FWE p values. maxTFCEScore(b+1) = max(tfce(:)); if SwE.WB.stat == 'T' maxTFCEScore_neg(b+1) = max(tfce_neg(:)); end - + clear tfce tfce_neg - + end - + else %isMat - + %-Print progress information in command window %------------------------------------------------------------------ str = sprintf('Bootstrap # %i', b); - + fprintf('%-40s: %1s',str,' '); - + Y_b = YWB + resWB .* repmat(resamplingMatrix(:,b),1,S); - + beta = pX * Y_b; %-Parameter estimates if WB.RSwE == 0 res=swe_fit(SwE, Y_b, xX.X, corr, beta, SwE.SS); - else + else res=swe_fit(SwE, Y_b, tmpR2, corr, beta, SwE.SS); end clear Y_b - + %-Estimation of the data variance-covariance components (modified SwE) %-SwE estimation (classic version) %-------------------------------------------------------------- @@ -2145,7 +2145,7 @@ function swe_cp_WB(SwE) (res(Indexk(Ind_Cov_vis_classic==i),:) .* res(Indexkk(Ind_Cov_vis_classic==i),:)); cCovBc = cCovBc + Cov_beta_i_tmp; if dof_type == 1 - if nSizeCon == 1 + if nSizeCon == 1 tmpSum = tmpSum + Cov_beta_i_tmp.^2/edof_Gr(i); else for iVox = 1:S @@ -2157,13 +2157,13 @@ function swe_cp_WB(SwE) end end end - + % These variables are left empty for classic SwE. Cov_vis = []; dofMat = []; clear Cov_beta_i_tmp end - + % if dof_type == 1, compute the edf now if dof_type == 1 if nSizeCon == 1 @@ -2185,9 +2185,9 @@ function swe_cp_WB(SwE) score = (conWB * beta) ./ sqrt(cCovBc); clear beta - + else - + cBeta = conWB * beta; clear beta score = zeros(1, S); @@ -2197,13 +2197,13 @@ function swe_cp_WB(SwE) cCovBc_vox = cCovBc_vox + cCovBc_vox' - diag(diag(cCovBc_vox)); score(iVox) = cBeta(:,iVox)' / cCovBc_vox * cBeta(:,iVox); end - score = score / rankCon; + score = score / rankCon; end - + if dof_type == 1 - hyptest = swe_hyptest(SwE, score, S, cCovBc, Cov_vis, edf); + hyptest = swe_hyptest(SwE, score, S, cCovBc, Cov_vis, edf); else - hyptest = swe_hyptest(SwE, score, S, cCovBc, Cov_vis, dofMat); + hyptest = swe_hyptest(SwE, score, S, cCovBc, Cov_vis, dofMat); end if (WB.clusterWise == 1) @@ -2218,10 +2218,10 @@ function swe_cp_WB(SwE) maxScore(b+1) = max(hyptest.positive.conScore); if (SwE.WB.stat == 'T') minScore(b+1) = min(hyptest.positive.conScore); - end - + end + end - + % compute the max cluster size if needed (so many ways this can be % done... Not sure this solution is the best) if (WB.clusterWise == 1) @@ -2245,7 +2245,7 @@ function swe_cp_WB(SwE) clusterSizesInSurfacesUnderH0 = [clusterSizesInSurfacesUnderH0, bootstrapedClusterStatistics.clusterSize]; maxClusterSizeInSurfaces(b+1) = bootstrapedClusterStatistics.maxClusterSize; end - + if isfield(bootstrapedClusterStatistics, 'clusterSizesInVolume') clusterSizesInVolumeUnderH0 = [clusterSizesInVolumeUnderH0, bootstrapedClusterStatistics.clusterSizesInVolume]; maxClusterSizeInVolume(b+1) = bootstrapedClusterStatistics.maxClusterSizeInVolume; @@ -2266,11 +2266,11 @@ function swe_cp_WB(SwE) end bootstrapedClusterStatisticsNeg = swe_getClusterStatistics(dataType, LocActivatedVoxelsNeg, dataTypeSpecificInformation, giftiAreaFile); - + if isfield(bootstrapedClusterStatisticsNeg, 'maxClusterSizeInVolume') maxClusterSizeInVolumeNeg(b+1) = bootstrapedClusterStatisticsNeg.maxClusterSizeInVolume; end - + if isfield(bootstrapedClusterStatisticsNeg, 'clusterSizesInSurfaces') clusterSizesInSurfacesNegUnderH0 = [clusterSizesInSurfacesNegUnderH0, bootstrapedClusterStatisticsNeg.clusterSizesInSurfaces]; maxClusterSizeInSurfacesNeg(b+1) = bootstrapedClusterStatisticsNeg.maxClusterSizeInSurfaces; @@ -2281,7 +2281,7 @@ function swe_cp_WB(SwE) clusterSizesInSurfacesNegUnderH0 = [clusterSizesInSurfacesNegUnderH0, bootstrapedClusterStatisticsNeg.clusterSize]; maxClusterSizeInSurfacesNeg(b+1) = bootstrapedClusterStatisticsNeg.maxClusterSize; end - + if isfield(bootstrapedClusterStatisticsNeg, 'clusterSizesInVolume') clusterSizesInVolumeNegUnderH0 = [clusterSizesInVolumeNegUnderH0, bootstrapedClusterStatisticsNeg.clusterSizesInVolume]; maxClusterSizeInVolumeNeg(b+1) = bootstrapedClusterStatisticsNeg.maxClusterSizeInVolume; @@ -2289,7 +2289,7 @@ function swe_cp_WB(SwE) clusterSizesInVolumeNegUnderH0 = [clusterSizesInVolumeNegUnderH0, bootstrapedClusterStatisticsNeg.clusterSize]; maxClusterSizeInVolumeNeg(b+1) = bootstrapedClusterStatisticsNeg.maxClusterSize; end - + maxClusterSizeNeg(b+1) = bootstrapedClusterStatisticsNeg.maxClusterSize; end @@ -2320,7 +2320,7 @@ function swe_cp_WB(SwE) end end -%-For cluster-wise analyses, normalise the cluster sizes +%-For cluster-wise analyses, normalise the cluster sizes %-------------------------------------------------------------------------- if WB.clusterWise == 1 scalingFactorNorm = swe_invNcdf(0.75); @@ -2333,7 +2333,7 @@ function swe_cp_WB(SwE) clusterSizesInSurfacesUnderH0_boxCox_median = median(clusterSizesInSurfacesUnderH0); clusterSizesInSurfacesUnderH0_boxCox_upperHalfIqr = swe_upperHalfIqr(clusterSizesInSurfacesUnderH0); clusterSizesInSurfaces_boxCox = swe_boxCoxTransform(SwE.WB.clusterInfo.clusterSizesInSurfaces, lambdaSurfacesUnderH0); - + SwE.WB.clusterInfo.maxClusterSizeInSurfaces = maxClusterSizeInSurfaces; SwE.WB.clusterInfo.clusterSizesInSurfacesUnderH0_boxCox_lambda = lambdaSurfacesUnderH0; SwE.WB.clusterInfo.clusterSizesInSurfacesUnderH0_boxCox_mean = clusterSizesInSurfacesUnderH0_boxCox_mean; @@ -2341,7 +2341,7 @@ function swe_cp_WB(SwE) SwE.WB.clusterInfo.clusterSizesInSurfacesUnderH0_boxCox_median = clusterSizesInSurfacesUnderH0_boxCox_median; SwE.WB.clusterInfo.clusterSizesInSurfacesUnderH0_boxCox_upperHalfIqr = clusterSizesInSurfacesUnderH0_boxCox_upperHalfIqr; maxClusterSizeInSurfaces_boxCox = swe_boxCoxTransform(maxClusterSizeInSurfaces, lambdaSurfacesUnderH0); - + if (clusterSizesInSurfacesUnderH0_boxCox_upperHalfIqr > 0) SwE.WB.clusterInfo.clusterSizesInSurfaces_norm = scalingFactorNorm * ... (clusterSizesInSurfaces_boxCox - clusterSizesInSurfacesUnderH0_boxCox_median) ./ clusterSizesInSurfacesUnderH0_boxCox_upperHalfIqr; @@ -2360,7 +2360,7 @@ function swe_cp_WB(SwE) end SwE.WB.clusterInfo.maxClusterSizeInSurfaces_norm = nan(1, WB.nB + 1); end - + if numel(clusterSizesInVolumeUnderH0) > 0 SwE.WB.clusterInfo.clusterSizesInVolumeUnderH0 = clusterSizesInVolumeUnderH0; lambdaVolumeUnderH0 = swe_estimateBoxCoxLambda(clusterSizesInVolumeUnderH0); @@ -2370,7 +2370,7 @@ function swe_cp_WB(SwE) clusterSizesInVolumeUnderH0_boxCox_median = median(clusterSizesInVolumeUnderH0); clusterSizesInVolumeUnderH0_boxCox_upperHalfIqr = swe_upperHalfIqr(clusterSizesInVolumeUnderH0); clusterSizesInVolume_boxCox = swe_boxCoxTransform(SwE.WB.clusterInfo.clusterSizesInVolume, lambdaVolumeUnderH0); - + SwE.WB.clusterInfo.maxClusterSizeInVolume = maxClusterSizeInVolume; SwE.WB.clusterInfo.clusterSizesInVolumeUnderH0_boxCox_lambda = lambdaVolumeUnderH0; SwE.WB.clusterInfo.clusterSizesInVolumeUnderH0_boxCox_mean = clusterSizesInVolumeUnderH0_boxCox_mean; @@ -2378,7 +2378,7 @@ function swe_cp_WB(SwE) SwE.WB.clusterInfo.clusterSizesInVolumeUnderH0_boxCox_median = clusterSizesInVolumeUnderH0_boxCox_median; SwE.WB.clusterInfo.clusterSizesInVolumeUnderH0_boxCox_upperHalfIqr = clusterSizesInVolumeUnderH0_boxCox_upperHalfIqr; maxClusterSizeInVolume_boxCox = swe_boxCoxTransform(maxClusterSizeInVolume, lambdaVolumeUnderH0); - + if (clusterSizesInVolumeUnderH0_boxCox_upperHalfIqr > 0) SwE.WB.clusterInfo.clusterSizesInVolume_norm = scalingFactorNorm * ... (clusterSizesInVolume_boxCox - clusterSizesInVolumeUnderH0_boxCox_median) ./ clusterSizesInVolumeUnderH0_boxCox_upperHalfIqr; @@ -2395,9 +2395,9 @@ function swe_cp_WB(SwE) else SwE.WB.clusterInfo.clusterSizesInVolume_norm = []; end - SwE.WB.clusterInfo.maxClusterSizeInVolume_norm = nan(1, WB.nB + 1); + SwE.WB.clusterInfo.maxClusterSizeInVolume_norm = nan(1, WB.nB + 1); end - + SwE.WB.clusterInfo.clusterSize_norm = [SwE.WB.clusterInfo.clusterSizesInSurfaces_norm, SwE.WB.clusterInfo.clusterSizesInVolume_norm]; SwE.WB.clusterInfo.maxClusterSize_norm = max(SwE.WB.clusterInfo.maxClusterSizeInSurfaces_norm, SwE.WB.clusterInfo.maxClusterSizeInVolume_norm); @@ -2413,7 +2413,7 @@ function swe_cp_WB(SwE) clusterSizesInSurfacesNegUnderH0_boxCox_median = median(clusterSizesInSurfacesNegUnderH0); clusterSizesInSurfacesNegUnderH0_boxCox_upperHalfIqr = swe_upperHalfIqr(clusterSizesInSurfacesNegUnderH0); clusterSizesInSurfacesNeg_boxCox = swe_boxCoxTransform(SwE.WB.clusterInfo.clusterSizesInSurfacesNeg, lambdaSurfacesNegUnderH0); - + SwE.WB.clusterInfo.maxClusterSizeInSurfacesNeg = maxClusterSizeInSurfacesNeg; SwE.WB.clusterInfo.clusterSizesInSurfacesNegUnderH0_boxCox_lambda = lambdaSurfacesNegUnderH0; SwE.WB.clusterInfo.clusterSizesInSurfacesNegUnderH0_boxCox_mean = clusterSizesInSurfacesNegUnderH0_boxCox_mean; @@ -2421,7 +2421,7 @@ function swe_cp_WB(SwE) SwE.WB.clusterInfo.clusterSizesInSurfacesNegUnderH0_boxCox_median = clusterSizesInSurfacesNegUnderH0_boxCox_median; SwE.WB.clusterInfo.clusterSizesInSurfacesNegUnderH0_boxCox_upperHalfIqr = clusterSizesInSurfacesNegUnderH0_boxCox_upperHalfIqr; maxClusterSizeInSurfacesNeg_boxCox = swe_boxCoxTransform(maxClusterSizeInSurfacesNeg, lambdaSurfacesNegUnderH0); - + if (clusterSizesInSurfacesNegUnderH0_boxCox_upperHalfIqr > 0) SwE.WB.clusterInfo.clusterSizesInSurfacesNeg_norm = scalingFactorNorm * ... (clusterSizesInSurfacesNeg_boxCox - clusterSizesInSurfacesNegUnderH0_boxCox_median) ./ clusterSizesInSurfacesNegUnderH0_boxCox_upperHalfIqr; @@ -2450,7 +2450,7 @@ function swe_cp_WB(SwE) clusterSizesInVolumeNegUnderH0_boxCox_median = median(clusterSizesInVolumeNegUnderH0); clusterSizesInVolumeNegUnderH0_boxCox_upperHalfIqr = swe_upperHalfIqr(clusterSizesInVolumeNegUnderH0); clusterSizesInVolumeNeg_boxCox = swe_boxCoxTransform(SwE.WB.clusterInfo.clusterSizesInVolumeNeg, lambdaVolumeNegUnderH0); - + SwE.WB.clusterInfo.maxClusterSizeInVolumeNeg = maxClusterSizeInVolumeNeg; SwE.WB.clusterInfo.clusterSizesInVolumeNegUnderH0_boxCox_lambda = lambdaVolumeNegUnderH0; SwE.WB.clusterInfo.clusterSizesInVolumeNegUnderH0_boxCox_mean = clusterSizesInVolumeNegUnderH0_boxCox_mean; @@ -2458,7 +2458,7 @@ function swe_cp_WB(SwE) SwE.WB.clusterInfo.clusterSizesInVolumeNegUnderH0_boxCox_median = clusterSizesInVolumeNegUnderH0_boxCox_median; SwE.WB.clusterInfo.clusterSizesInVolumeNegUnderH0_boxCox_upperHalfIqr = clusterSizesInVolumeNegUnderH0_boxCox_upperHalfIqr; maxClusterSizeInVolumeNeg_boxCox = swe_boxCoxTransform(maxClusterSizeInVolumeNeg, lambdaVolumeNegUnderH0); - + if (clusterSizesInVolumeNegUnderH0_boxCox_upperHalfIqr > 0) SwE.WB.clusterInfo.clusterSizesInVolumeNeg_norm = scalingFactorNorm * ... (clusterSizesInVolumeNeg_boxCox - clusterSizesInVolumeNegUnderH0_boxCox_median) ./ clusterSizesInVolumeNegUnderH0_boxCox_upperHalfIqr; @@ -2495,34 +2495,34 @@ function swe_cp_WB(SwE) lP_wb_pos(:,cmask) = -log10(uncP); save(sprintf('swe_%s_%cstat_lp-WB_c%02d%s', file_data_type, WB.stat, 1, file_ext), 'lP_wb_pos'); clear lP_wb_pos - + if WB.stat == 'T' uncP_neg = 1 + 1/(WB.nB + 1) - uncP; lP_wb_neg = zeros(1, nDataElements); lP_wb_neg(:,cmask) = -log10(uncP_neg); save(sprintf('swe_%s_%cstat_lp-WB_c%02d%s', file_data_type, WB.stat, 2, file_ext), 'lP_wb_neg'); clear lP_wb_neg - + Z_wb = zeros(1, nDataElements); Z_wb(:,cmask) = swe_invNcdf(1 - uncP); save(sprintf('swe_%s_z%cstat-WB_c%02d%s', file_data_type, WB.stat, 1, file_ext), 'Z_wb'); clear Z_wb - + else - + X_wb = zeros(1, nDataElements); X_wb(:,cmask) = spm_invXcdf(1 - uncP,1); save(sprintf('swe_%s_x%cstat-WB_c%02d%s', file_data_type, WB.stat, 1, file_ext), 'X_wb'); clear X_wb - + end - + % % - write out lP_FWE+ and lP_FWE- ; % - + FWERP = ones(1, S); % 1 because the original maxScore is always > original Score - + for b = 1:WB.nB %-FWER-corrected p is proportion of randomisation greater or % equal to statistic. @@ -2535,11 +2535,11 @@ function swe_cp_WB(SwE) lP_wb_FWE_pos(:,cmask) = -log10(FWERP); save(sprintf('swe_%s_%cstat_lpFWE-WB_c%02d%s',file_data_type,WB.stat,1,file_ext), 'lP_wb_FWE_pos'); clear lP_wb_FWE_pos fwerP_pos FWERP - - + + if WB.stat == 'T' FWERPNeg = ones(1, S); % 1 because the original maxScore is always > original Score - + for b = 1:WB.nB %-FWER-corrected p is proportion of randomisation greater or % equal to statistic. @@ -2553,7 +2553,7 @@ function swe_cp_WB(SwE) save(sprintf('swe_%s_%cstat_lpFWE-WB_c%02d%s',file_data_type,WB.stat,2,file_ext), 'lP_wb_FWE_neg'); clear lP_wb_FWE_neg fwerP_neg end - + % % - write out lP_FDR+ and lP_FDR- images; % @@ -2566,7 +2566,7 @@ function swe_cp_WB(SwE) lP_wb_FDR_pos(:,cmask) = -log10(fdrP); save(sprintf('swe_%s_%cstat_lpFDR-WB_c%02d%s',file_data_type,WB.stat,1,file_ext), 'lP_wb_FDR_pos'); clear lP_wb_FDR_pos fdrP_pos fdrP - + if WB.stat =='T' try fdrP = spm_P_FDR(1 + 1/(WB.nB + 1) - uncP); @@ -2578,14 +2578,14 @@ function swe_cp_WB(SwE) save(sprintf('swe_%s_%cstat_lpFDR-WB_c%02d%s',file_data_type,WB.stat,2,file_ext), 'lP_wb_FDR_neg'); clear lP_wb_FDR_neg fdrP_neg fdrP end - + if WB.clusterWise == 1 % Not sure what to output. So might be changed later. % For now, -log(p_{cluster-wise FWER}) image with nan for non-surviving % voxels after the thresholding of the original data - + if canUseBoxCoxNormalisation - + normClusterFwerP_pos_perCluster = ones(1, SwE.WB.clusterInfo.nCluster); % 1 because the original maxScore is always > original Score if (~isempty(SwE.WB.clusterInfo.clusterSize)) for b = 1:WB.nB @@ -2593,7 +2593,7 @@ function swe_cp_WB(SwE) end normClusterFwerP_pos_perCluster = normClusterFwerP_pos_perCluster / (WB.nB + 1); end - + lP_wb_normClusterFWE_pos = zeros(1, nDataElements); tmp = find(cmask); tmp3 = zeros(1, size(SwE.WB.clusterInfo.LocActivatedVoxels,2)); @@ -2602,9 +2602,9 @@ function swe_cp_WB(SwE) end lP_wb_normClusterFWE_pos(tmp(originalActivatedVoxels)) = -log10(tmp3); save(sprintf('swe_clusternorm_%cstat_lpFWE-WB_c%02d%s',WB.stat,1,file_ext), 'lP_wb_normClusterFWE_pos'); - + else - + clusterFwerP_pos_perCluster = ones(1, SwE.WB.clusterInfo.nCluster); % 1 because the original maxScore is always > original Score if (~isempty(SwE.WB.clusterInfo.clusterSize)) for b = 1:WB.nB @@ -2612,7 +2612,7 @@ function swe_cp_WB(SwE) end clusterFwerP_pos_perCluster = clusterFwerP_pos_perCluster / (WB.nB + 1); end - + lP_wb_clusterFWE_pos = zeros(1, nDataElements); tmp = find(cmask); tmp3 = zeros(1, size(SwE.WB.clusterInfo.LocActivatedVoxels,2)); @@ -2620,12 +2620,12 @@ function swe_cp_WB(SwE) tmp3(SwE.WB.clusterInfo.clusterAssignment == iC) = clusterFwerP_pos_perCluster(iC); end lP_wb_clusterFWE_pos(tmp(originalActivatedVoxels)) = -log10(tmp3); - save(sprintf('swe_clustere_%cstat_lpFWE-WB_c%02d%s',WB.stat,1,file_ext), 'lP_wb_clusterFWE_pos'); - + save(sprintf('swe_clustere_%cstat_lpFWE-WB_c%02d%s',WB.stat,1,file_ext), 'lP_wb_clusterFWE_pos'); + end if WB.stat =='T' - + if canUseBoxCoxNormalisationNeg normClusterFwerP_neg_perCluster = ones(1, SwE.WB.clusterInfo.nClusterNeg); % 1 because the original maxScore is always > original Score @@ -2635,7 +2635,7 @@ function swe_cp_WB(SwE) end normClusterFwerP_neg_perCluster = normClusterFwerP_neg_perCluster / (WB.nB + 1); end - + lP_wb_normClusterFWE_neg = zeros(1, nDataElements); tmp = find(cmask); tmp3 = zeros(1, size(SwE.WB.clusterInfo.LocActivatedVoxelsNeg,2)); @@ -2644,9 +2644,9 @@ function swe_cp_WB(SwE) end lP_wb_normClusterFWE_neg(tmp(originalActivatedVoxelsNeg)) = -log10(tmp3); save(sprintf('swe_clusternorm_%cstat_lpFWE-WB_c%02d%s',WB.stat,2,file_ext), 'lP_wb_normClusterFWE_neg'); - + else - + clusterFwerP_neg_perCluster = ones(1, SwE.WB.clusterInfo.nClusterNeg); % 1 because the original maxScore is always > original Score if (~isempty(SwE.WB.clusterInfo.clusterSizeNeg)) for b = 1:WB.nB @@ -2654,7 +2654,7 @@ function swe_cp_WB(SwE) end clusterFwerP_neg_perCluster = clusterFwerP_neg_perCluster / (WB.nB + 1); end - + lP_wb_clusterFWE_neg = zeros(1, nDataElements); tmp = find(cmask); tmp3 = zeros(1, size(SwE.WB.clusterInfo.LocActivatedVoxelsNeg,2)); @@ -2667,7 +2667,7 @@ function swe_cp_WB(SwE) end end else - + Q = cumprod([1,SwE.xVol.DIM(1:2)'])*XYZ - ... sum(cumprod(SwE.xVol.DIM(1:2)')); % @@ -2677,17 +2677,17 @@ function swe_cp_WB(SwE) tmp = zeros(SwE.xVol.DIM'); tmp(Q) = -log10(uncP); swe_data_write(VlP_wb_pos, tmp); - + % If it's F, write out an X map. stat = zeros(SwE.xVol.DIM'); if WB.stat == 'F' stat(Q) = spm_invXcdf(1 - uncP,1); swe_data_write(VcScore_wb_pos, stat); end - + % If it's T, write out a Z map. if WB.stat == 'T' - + % Positive map. stat(Q) = swe_invNcdf(1 - uncP); swe_data_write(VcScore_wb_pos, stat); @@ -2696,13 +2696,13 @@ function swe_cp_WB(SwE) tmp(Q) = -log10(1 + 1/(WB.nB + 1) - uncP); swe_data_write(VlP_wb_neg, tmp); end - + % % - write out lP_FWE+ and lP_FWE- images; % - + FWERP = ones(1, S); % 1 because the original maxScore is always > original Score - + for b = 1:WB.nB %-FWER-corrected p is proportion of randomisation greater or % equal to statistic. @@ -2713,27 +2713,27 @@ function swe_cp_WB(SwE) FWERP = FWERP / (WB.nB + 1); tmp(Q) = -log10(FWERP); swe_data_write(VlP_wb_FWE_pos, tmp); - + % FWE correction for TFCE images. if TFCE - + % Make new tfce fwe p volume. (Initiate to one to account for % original analysis). tfcefwevol = ones(DIM(1), DIM(2), DIM(3)); - + % Calculate FWE p values. for b = 1:WB.nB tfcefwevol = tfcefwevol + (maxTFCEScore(b+1) > par_tfce - tol); end tfcefwevol = tfcefwevol / (WB.nB + 1); - + % Write out volume. swe_data_write(VlP_tfce_FWE_pos, -log10(tfcefwevol)); - + % Same again for negative contrast, if we are using a T statistic. if WB.stat == 'T' - - % Make new negative tfce fwe p volume. (Initiate to one to + + % Make new negative tfce fwe p volume. (Initiate to one to % account for original analysis). tfcefwevol_neg = ones(DIM(1), DIM(2), DIM(3)); @@ -2747,10 +2747,10 @@ function swe_cp_WB(SwE) swe_data_write(VlP_tfce_FWE_neg, -log10(tfcefwevol_neg)); end end - + if WB.stat == 'T' FWERPNeg = ones(1, S); % 1 because the original maxScore is always > original Score - + for b = 1:WB.nB %-FWER-corrected p is proportion of randomisation greater or % equal to statistic. @@ -2762,7 +2762,7 @@ function swe_cp_WB(SwE) tmp(Q) = -log10(FWERPNeg); swe_data_write(VlP_wb_FWE_neg, tmp); end - + % % - write out lP_FDR+ and lP_FDR- images; % @@ -2772,7 +2772,7 @@ function swe_cp_WB(SwE) tmp(Q) = -log10(spm_P_FDR(uncP,[],'P',[],sort(uncP)')); end swe_data_write(VlP_wb_FDR_pos, tmp); - + if WB.stat =='T' try tmp(Q) = -log10(spm_P_FDR(1 + 1/(WB.nB + 1) - uncP)); @@ -2781,7 +2781,7 @@ function swe_cp_WB(SwE) end swe_data_write(VlP_wb_FDR_neg, tmp); end - + if WB.clusterWise == 1 % Not sure what to output. So might be changed later. % For now, -log(p_{cluster-wise FWER}) image with nan for non-surviving @@ -2789,7 +2789,7 @@ function swe_cp_WB(SwE) Q = cumprod([1,SwE.xVol.DIM(1:2)']) * SwE.WB.clusterInfo.LocActivatedVoxels - ... sum(cumprod(SwE.xVol.DIM(1:2)')); tmp= zeros(SwE.xVol.DIM'); - + tmp4 = zeros(1, size(SwE.WB.clusterInfo.LocActivatedVoxels,2)); for iC = 1:SwE.WB.clusterInfo.nCluster tmp4(SwE.WB.clusterInfo.clusterAssignment == iC) = SwE.WB.clusterInfo.clusterSize(iC); @@ -2798,7 +2798,7 @@ function swe_cp_WB(SwE) swe_data_write(V_clustere_pos, tmp); if canUseBoxCoxNormalisation - + VlP_wb_normClusterFWE_pos = swe_data_hdr_write(sprintf('swe_clusternorm_%cstat_lpFWE-WB_c%02d%s', WB.stat, 1, file_ext), DIM, M,... sprintf('Non-parametric normalised clusterwise FWE -log10(P) value data (positive, CFT %g).',... SwE.WB.clusterInfo.primaryThreshold), metadata); @@ -2814,7 +2814,7 @@ function swe_cp_WB(SwE) end normClusterFwerP_pos_perCluster = normClusterFwerP_pos_perCluster / (WB.nB + 1); end - tmp2 = -log10(normClusterFwerP_pos_perCluster); + tmp2 = -log10(normClusterFwerP_pos_perCluster); tmp3 = zeros(1, size(SwE.WB.clusterInfo.LocActivatedVoxels,2)); tmp4 = zeros(1, size(SwE.WB.clusterInfo.LocActivatedVoxels,2)); for iC = 1:SwE.WB.clusterInfo.nCluster @@ -2825,7 +2825,7 @@ function swe_cp_WB(SwE) swe_data_write(VlP_wb_normClusterFWE_pos, tmp); tmp(Q) = tmp4; swe_data_write(V_normCluster_pos, tmp); - + else VlP_wb_clusterFWE_pos = swe_data_hdr_write(sprintf('swe_clustere_%cstat_lpFWE-WB_c%02d%s', WB.stat, 1, file_ext), DIM, M,... @@ -2839,7 +2839,7 @@ function swe_cp_WB(SwE) end clusterFwerP_pos_perCluster = clusterFwerP_pos_perCluster / (WB.nB + 1); end - tmp2 = -log10(clusterFwerP_pos_perCluster); + tmp2 = -log10(clusterFwerP_pos_perCluster); tmp3 = zeros(1, size(SwE.WB.clusterInfo.LocActivatedVoxels,2)); for iC = 1:SwE.WB.clusterInfo.nCluster tmp3(SwE.WB.clusterInfo.clusterAssignment == iC) = tmp2(iC); @@ -2848,12 +2848,12 @@ function swe_cp_WB(SwE) swe_data_write(VlP_wb_clusterFWE_pos, tmp); end - + if WB.stat =='T' Q = cumprod([1,SwE.xVol.DIM(1:2)']) * SwE.WB.clusterInfo.LocActivatedVoxelsNeg - ... sum(cumprod(SwE.xVol.DIM(1:2)')); tmp= zeros(SwE.xVol.DIM'); - + tmp4 = zeros(1, size(SwE.WB.clusterInfo.LocActivatedVoxelsNeg, 2)); for iC = 1:SwE.WB.clusterInfo.nClusterNeg tmp4(SwE.WB.clusterInfo.clusterAssignmentNeg == iC) = SwE.WB.clusterInfo.clusterSizeNeg(iC); @@ -2912,10 +2912,10 @@ function swe_cp_WB(SwE) swe_data_write(VlP_wb_clusterFWE_neg, tmp); end - + end end - + if TFCE tfce_luncP = -log10((tfce_uncP+1)./(WB.nB+1)); swe_data_write(VlP_tfce_pos, tfce_luncP); @@ -2924,7 +2924,7 @@ function swe_cp_WB(SwE) swe_data_write(VlP_tfce_neg, tfce_luncP_neg); end end - + end % save the version number of the toolbox @@ -2966,7 +2966,7 @@ function swe_cp_WB(SwE) %-Mask out voxels where data is constant in at least one separable % matrix design either in a visit category or within-subject (BG - 27/05/2016) function [cmask,Y,CrS]=swe_mask_seperable(SwE, cmask, Y, iGr_dof) - + % Setup nGr_dof = length(unique(iGr_dof)); if isfield(SwE.type,'modified') @@ -3011,7 +3011,7 @@ function swe_cp_WB(SwE) end end end - + Y = Y(:,cmask); %-Data within mask CrS = sum(cmask); %-# current voxels @@ -3030,7 +3030,7 @@ function swe_cp_WB(SwE) if isfield(SwE.type,'modified') dof_type = SwE.type.modified.dof_mo; nGr = length(unique(SwE.Gr.iGr)); - + if nSizeCon == 1 Wg_2 = SwE.WB.Wg{1}; Wg_3 = SwE.WB.Wg{1}; @@ -3038,19 +3038,19 @@ function swe_cp_WB(SwE) Wg_2 = SwE.WB.Wg{2}; Wg_3 = SwE.WB.Wg{3}; end - + else dof_type = SwE.type.classic.dof_cl; - nGr = SwE.Subj.nSubj; + nGr = SwE.Subj.nSubj; end % Convert P values. switch dof_type - + case 0 - + edf = SwE.xX.erdf_niave; - + case 1 % for dof_type = 1, saved in dofMat. Should be changed later edf = dofMat; @@ -3078,7 +3078,7 @@ function swe_cp_WB(SwE) (tmp(:)' * swe_duplication_matrix(nSizeCon) * cCovBc).^2) ./ CovcCovBc; end end - + % P values and activated voxels (if clusterwise). if (SwE.WB.stat == 'T') % We calculate both p and negative p as both will experience overflow @@ -3086,10 +3086,10 @@ function swe_cp_WB(SwE) % results. p = spm_Tcdf(-score, edf); negp = spm_Tcdf(score, edf); - + if SwE.WB.clusterWise~=0 if nargin <=6 - % We may wish to just record the activated voxels. + % We may wish to just record the activated voxels. activatedVoxels = p < (SwE.WB.clusterInfo.primaryThreshold); activatedVoxelsNeg = negp < (SwE.WB.clusterInfo.primaryThreshold); else @@ -3098,7 +3098,7 @@ function swe_cp_WB(SwE) activatedVoxelsNeg = [varargin{2}, negp < (SwE.WB.clusterInfo.primaryThreshold)]; end end - + else scoreTmp = (edf-rankCon+1) ./ edf .* score; scoreTmp(scoreTmp < 0 ) = 0; @@ -3116,31 +3116,31 @@ function swe_cp_WB(SwE) activatedVoxels = [varargin{1}, p < (SwE.WB.clusterInfo.primaryThreshold)]; end end - + negp = 1-p; activatedVoxelsNeg = NaN; - + end % Calculate converted score in a way which minimizes under/overflow if (SwE.WB.stat == 'T') - + conScore = zeros(size(score)); - + switch dof_type % In case 0 edf is the same for everything case 0 - + if any(score>0) conScore(score>0) = -spm_invNcdf(spm_Tcdf(-score(score>0),edf)); end if any(score<=0) conScore(score<=0) = spm_invNcdf(spm_Tcdf(score(score<=0),edf)); end - + % In all other cases edf varies over voxels otherwise - + if any(score>0) conScore(score>0) = -spm_invNcdf(spm_Tcdf(-score(score>0),edf(score>0))); end @@ -3149,11 +3149,11 @@ function swe_cp_WB(SwE) end end - + else conScore = swe_invNcdf(0.5 * p).^2; end - + % Save results hyptest = struct('positive', struct('p', p,... 'edf', edf,... @@ -3163,7 +3163,7 @@ function swe_cp_WB(SwE) if SwE.WB.clusterWise~=0 hyptest.positive.activatedVoxels = activatedVoxels; end - + % If it's a T contrast save negative results as well. if (SwE.WB.stat == 'T') hyptest.negative = struct('p', negp,... @@ -3176,10 +3176,10 @@ function swe_cp_WB(SwE) end else - + % For an F, negative p-values are still useful hyptest.negative = struct('p', negp); - + end end @@ -3193,7 +3193,7 @@ function swe_cp_WB(SwE) conWB = SwE.WB.con; iSubj = SwE.Subj.iSubj; nSubj = SwE.Subj.nSubj; -uSubj = SwE.Subj.uSubj; +uSubj = SwE.Subj.uSubj; rankCon = rank(SwE.WB.con); % This is to prevent tmpR2 being overwritten. @@ -3219,7 +3219,7 @@ function swe_cp_WB(SwE) if restric == 1 corr = repmat(sqrt(nScan/(nScan - nBeta + rankCon)),nScan,1); % residual correction (type 1) else - corr = repmat(sqrt(nScan/(nScan-nBeta)),nScan,1); + corr = repmat(sqrt(nScan/(nScan-nBeta)),nScan,1); end case 2 corr = (1-diag(Hat)).^(-0.5); % residual correction (type 2) @@ -3247,7 +3247,7 @@ function swe_cp_WB(SwE) % This function obtains Y estimates and residuals from fitting data. function [res, Y_est]=swe_fit(SwE, Y, crctX, corr, beta, ss) - + Y_est = crctX * beta; if ss >= 4 % SC2 or SC3 res = zeros(size(Y)); @@ -3264,4 +3264,4 @@ function swe_cp_WB(SwE) function value = swe_upperHalfIqr(X) value = diff( quantile(X, [.5 .75]) ); -end \ No newline at end of file +end diff --git a/swe_create_vol.m b/swe_create_vol.m index 5b4e7bd..64c1f62 100644 --- a/swe_create_vol.m +++ b/swe_create_vol.m @@ -3,7 +3,7 @@ % ========================================================================= % FORMAT vol = swe_create_vol(fname, DIM, M [, descrip]) % ------------------------------------------------------------------------- -% Inputs: +% Inputs: % - fname: Filename of new image % - DIM: Row vector giving image dimensions % - M: 4x4 homogeneous transformation, from V.mat @@ -12,7 +12,7 @@ % - metadata: metadata from GIfTI file (SPM set metadata = {} for NIfTI) % ========================================================================= % Version Info: $Format:%ci$ $Format:%h$ - + if nargin > 3 descrip = varargin{1}; else @@ -45,5 +45,5 @@ else vol = spm_create_vol(vol); end - + end diff --git a/swe_data_hdr_read.m b/swe_data_hdr_read.m index f3eb1f1..8c99fc0 100644 --- a/swe_data_hdr_read.m +++ b/swe_data_hdr_read.m @@ -4,7 +4,7 @@ % P - a char or cell array of filenames % V - a structure array containing data information % - % This function behaves like spm_data_hdr_read but can also read the + % This function behaves like spm_data_hdr_read but can also read the % headers of CIfTI files. % ========================================================================= % Bryan Guillaume @@ -26,14 +26,14 @@ if ~(strcmpi(file_ext,'.dtseries.nii') || strcmpi(file_ext,'.dscalar.nii')) error('Not all files are CIfTI files'); end - + % if some slices of data are requested, remove this from the filename if ~isempty(sliceInd) P2{i} = P2{i}(1:(end-numel(sliceInd))); end - + ciftiObject = swe_cifti(P2{i}, readExt); - + if isempty(sliceInd) sliceInd = 1:ciftiObject.dat.dim(5); else @@ -52,7 +52,7 @@ it = it + 1; end end - else + else V = spm_data_hdr_read(P); end end @@ -67,4 +67,4 @@ 'n', [1 1],... 'descrip', '',... 'private', []); -end \ No newline at end of file +end diff --git a/swe_data_hdr_write.m b/swe_data_hdr_write.m index a502158..f599812 100644 --- a/swe_data_hdr_write.m +++ b/swe_data_hdr_write.m @@ -3,7 +3,7 @@ % ========================================================================= % FORMAT V = swe_data_hdr_write(fname, DIM, M, descrip, metadata[, dataType]) % ------------------------------------------------------------------------- - % Inputs: + % Inputs: % - fname: Filename of new image % - DIM: Row vector giving image dimensions % - M: 4x4 homogeneous transformation, from V.mat @@ -26,7 +26,7 @@ 'pinfo', [1 0 0]',... 'descrip', descrip,... metadata{:}); - + if isfield(V, 'ciftiTemplate') [~, sliceInd] = swe_get_file_extension(V.ciftiTemplate); if isempty(sliceInd) @@ -34,7 +34,7 @@ else sourceName = V.ciftiTemplate(1:( end - numel(sliceInd) )); end - + % make sure we select only one slice V = swe_data_hdr_read(sprintf('%s,1', sourceName), true); @@ -47,7 +47,7 @@ 1,... 0); V.n = [1 1]; - % modify the xml if the number of point is > 1 + % modify the xml if the number of point is > 1 xml = char(V.private.hdr.ext.edata(:)'); ind = strfind(xml, 'NumberOfSeriesPoints="'); if ~isempty(ind) @@ -73,13 +73,13 @@ end end V.private.hdr.ext.edata = uint8(xml)'; - - create(V.private) + + create(V.private) V.private.dat(:) = 0; V = swe_data_hdr_read(fname, false); V.descrip = descrip; else V = spm_data_hdr_write(V); end - -end \ No newline at end of file + +end diff --git a/swe_data_read.m b/swe_data_read.m index 56af48a..ae7f464 100644 --- a/swe_data_read.m +++ b/swe_data_read.m @@ -25,7 +25,7 @@ if ~isstruct(V) V = swe_data_hdr_read(V); end - + if isa(V(1).private, 'swe_cifti') if isempty(varargin) Y = zeros(prod(V(1).dim), numel(V)); % to be coherent with spm_read_vols and gifti format diff --git a/swe_data_write.m b/swe_data_write.m index 5559c47..4ab4545 100644 --- a/swe_data_write.m +++ b/swe_data_write.m @@ -11,7 +11,7 @@ % ========================================================================= % Bryan Guillaume % Version Info: $Format:%ci$ $Format:%h$ - + if isfield(V,'private') cl = class(V.private); elseif isfield(V,'dat') @@ -30,7 +30,7 @@ V.private.dat(varargin{1}) = reshape(Y,size(varargin{1}))'; end end - elseif strcmpi(cl, 'gifti') && ~isempty(varargin) && exist('OCTAVE_VERSION','builtin') + elseif strcmpi(cl, 'gifti') && ~isempty(varargin) && exist('OCTAVE_VERSION','builtin') D = V.private.private.data{1}; tmp = reshape(Y,size(varargin{1})); try @@ -42,4 +42,4 @@ else V = spm_data_write(V, Y, varargin{:}); end -end \ No newline at end of file +end diff --git a/swe_defaults.m b/swe_defaults.m index 99b06ad..d0f868c 100644 --- a/swe_defaults.m +++ b/swe_defaults.m @@ -15,10 +15,10 @@ global SwEdefs -% If true, shuffles the seed of the random number generator to get -% different results every time. Use false, if you want to specify your own -% seed, for instance to insure that results can be replicated or when using +% If true, shuffles the seed of the random number generator to get +% different results every time. Use false, if you want to specify your own +% seed, for instance to insure that results can be replicated or when using % a high performance cluster. %------------------------------------------------------------------------ -SwEdefs.shuffle_seed = true; +SwEdefs.shuffle_seed = true; diff --git a/swe_duplication_matrix.m b/swe_duplication_matrix.m index ed1128c..150c7ce 100644 --- a/swe_duplication_matrix.m +++ b/swe_duplication_matrix.m @@ -22,4 +22,4 @@ end it = it + n - j; end -end \ No newline at end of file +end diff --git a/swe_estimateBoxCoxLambda.m b/swe_estimateBoxCoxLambda.m index 83e7552..4b0e9db 100644 --- a/swe_estimateBoxCoxLambda.m +++ b/swe_estimateBoxCoxLambda.m @@ -10,7 +10,7 @@ % ========================================================================= % Bryan Guillaume % Version Info: $Format:%ci$ $Format:%h$ - + if min(data) <= 0 error('The data must be positive.') end @@ -18,13 +18,13 @@ startingLambdaValue = 0; boxCoxLambda = fminsearch( @(lambda) getMinusBoxCoxLogLikelihood(data, lambda), startingLambdaValue); - function minusBoxCoxLogLikelihood = getMinusBoxCoxLogLikelihood(data, lambda) + function minusBoxCoxLogLikelihood = getMinusBoxCoxLogLikelihood(data, lambda) nData = numel(data); - + transformedData = swe_boxCoxTransform(data, lambda); minusBoxCoxLogLikelihood = 0.5 * nData * log(var(transformedData, 0)) + (1 - lambda) * sum(log(data)); - + end -end \ No newline at end of file +end diff --git a/swe_getClusterStatistics.m b/swe_getClusterStatistics.m index fa32609..fe1afa0 100644 --- a/swe_getClusterStatistics.m +++ b/swe_getClusterStatistics.m @@ -24,14 +24,14 @@ elseif dataType == swe_DataType('Gifti') || dataType == swe_DataType('SurfaceMat') [clusterStatistics.clusterAssignment, ~, clusterAreas] = swe_mesh_clusters(dataTypeSpecificInformation, locationActivatedElements, giftiAreaFile); - + if ~isnan(clusterAreas) clusterStatistics.clusterAreas = clusterAreas; clusterStatistics.maxClusterArea = max(clusterStatistics.clusterAreas); end elseif dataType == swe_DataType('Cifti') - + [clusterStatistics.clusterAssignment, ... clusterStatistics.clusterSizesInSurfaces, ... clusterStatistics.clusterSizesInVolume] = ... @@ -39,13 +39,13 @@ if isempty(clusterStatistics.clusterSizesInSurfaces) clusterStatistics.maxClusterSizeInSurfaces = 0; - else + else clusterStatistics.maxClusterSizeInSurfaces = max(clusterStatistics.clusterSizesInSurfaces); end if isempty(clusterStatistics.clusterSizesInVolume) clusterStatistics.maxClusterSizeInVolume = 0; - else + else clusterStatistics.maxClusterSizeInVolume = max(clusterStatistics.clusterSizesInVolume); end @@ -58,7 +58,7 @@ end clusterStatistics.nCluster = max(clusterStatistics.clusterAssignment); - + if clusterStatistics.nCluster > 0 clusterStatistics.clusterSize = histc(clusterStatistics.clusterAssignment, 1:clusterStatistics.nCluster); clusterStatistics.maxClusterSize = max(clusterStatistics.clusterSize); @@ -66,4 +66,4 @@ clusterStatistics.clusterSize = []; clusterStatistics.maxClusterSize = 0; end -end \ No newline at end of file +end diff --git a/swe_getSPM.m b/swe_getSPM.m index 6fda5b3..dfd761b 100644 --- a/swe_getSPM.m +++ b/swe_getSPM.m @@ -7,13 +7,13 @@ % FORMAT [SwE,xSwE] = swe_getSPM(xSwE); % ------------------------------------------------------------------------- % -% Query SwE in batch mode. See below for a description of fields that -% may be present in xSwE input. Values for missing fields will be +% Query SwE in batch mode. See below for a description of fields that +% may be present in xSwE input. Values for missing fields will be % queried interactively. % -% xSwE - structure containing spm, distribution & filtering +% xSwE - structure containing spm, distribution & filtering % details -% .swd - SwE working directory - directory containing current +% .swd - SwE working directory - directory containing current % SwE.mat % .title - title for comparison (string) % .Z - minimum of Statistics {filtered on u and k} @@ -85,7 +85,7 @@ % .thresDesc - description of height threshold (string) % % In addition, the xCon structure is updated. For newly evaluated -% contrasts, SwE images (swe_{vox|dpx|dat}_{T|F}stat_c{c#}) are written, along +% contrasts, SwE images (swe_{vox|dpx|dat}_{T|F}stat_c{c#}) are written, along % with contrast (swe_{vox|dpx|dat}_beta_c{c#}) images. % % For a parametric analysis the following is added to the xCon @@ -125,18 +125,18 @@ % parameters (with contrast weights c) is estimated by c'b, where b are % the parameter estimates given by b=pinv(X)*Y. % -% For a paramertic analysis, either single contrasts can be examined -% or conjunctions of different contrasts. Contrasts are estimable -% linear combinations of the parameters, and are specified using -% the SwE contrast manager interface [swe_conman.m]. For a +% For a paramertic analysis, either single contrasts can be examined +% or conjunctions of different contrasts. Contrasts are estimable +% linear combinations of the parameters, and are specified using +% the SwE contrast manager interface [swe_conman.m]. For a % non-parametric analysis, two contrasts are recorded; activation % and deactivation for the contrast vector specified in the batch -% window. These are recorded a priori in a seperate function +% window. These are recorded a priori in a seperate function % with certain thresholds applied here [swe_contrasts_WB]. % -% SwE parametric maps are generated for the null hypotheses that -% the contrast is zero (or zero vector in the case of F-contrasts). -% See the help for the contrast manager [swe_conman.m] for a further +% SwE parametric maps are generated for the null hypotheses that +% the contrast is zero (or zero vector in the case of F-contrasts). +% See the help for the contrast manager [swe_conman.m] for a further % details on contrasts and contrast specification. % % A conjunction assesses the conjoint expression of multiple effects. The @@ -248,11 +248,11 @@ try SwE.xVol.XYZ; catch - + %-Check the model has been estimated %---------------------------------------------------------------------- error( 'This model has not been estimated.'); - + end % check format of data @@ -292,7 +292,7 @@ % Tolerance for comparing real numbers for WB analyses % Use a value < to the smallest WB p-value as it will be used to include WB p-values equal to alpha -if isfield(SwE, 'WB') +if isfield(SwE, 'WB') tol = 0.1 / (SwE.WB.nB + 1); end @@ -326,7 +326,7 @@ else save('SwE.mat', 'SwE'); end - xCon = SwE.xCon; + xCon = SwE.xCon; end catch if isfield(SwE, 'WB') && ~exist('OCTAVE_VERSION','builtin') @@ -337,7 +337,7 @@ else save('SwE.mat', 'SwE'); end - xCon = SwE.xCon; + xCon = SwE.xCon; else xCon = {}; end @@ -346,7 +346,7 @@ try Ic = xSwE.Ic; catch - + % If we're not doing wild bootstrap and not in octave, ask for a contrast. if ~isfield(SwE, 'WB') && ~exist('OCTAVE_VERSION','builtin') [Ic,xCon] = swe_conman(SwE,'T&F',Inf,... @@ -408,7 +408,7 @@ end end else - n = 1; + n = 1; end % not sure we want to do that with the SwE (commented for now) @@ -416,9 +416,9 @@ % (Orthogonality within subspace spanned by contrasts) %-------------------------------------------------------------------------- % if nc > 1 && n > 1 && ~spm_FcUtil('|_?',xCon(Ic), xX.xKXs) -% +% % OrthWarn = 0; -% +% % %-Successively orthogonalise % %-NB: This loop is peculiarly controlled to account for the % % possibility that Ic may shrink if some contrasts disappear @@ -426,32 +426,32 @@ % %---------------------------------------------------------------------- % i = 1; % while(i < nc), i = i + 1; -% +% % %-Orthogonalise (subspace spanned by) contrast i w.r.t. previous % %------------------------------------------------------------------ % oxCon = spm_FcUtil('|_',xCon(Ic(i)), xX.xKXs, xCon(Ic(1:i-1))); -% +% % %-See if this orthogonalised contrast has already been entered % % or is colinear with a previous one. Define a new contrast if % % neither is the case. % %------------------------------------------------------------------ % d = spm_FcUtil('In',oxCon,xX.xKXs,xCon); -% +% % if spm_FcUtil('0|[]',oxCon,xX.xKXs) -% +% % %-Contrast was colinear with a previous one - drop it % %-------------------------------------------------------------- % Ic(i) = []; % i = i - 1; -% +% % elseif any(d) -% +% % %-Contrast unchanged or already defined - note index % %-------------------------------------------------------------- % Ic(i) = min(d); -% +% % else -% +% % %-Define orthogonalised contrast as new contrast % %-------------------------------------------------------------- % OrthWarn = OrthWarn + 1; @@ -461,13 +461,13 @@ % xCon = [xCon, oxCon]; % Ic(i) = length(xCon); % end -% +% % end % while... -% +% % if OrthWarn % warning('SwE:ConChange','%d contrasts orthogonalized',OrthWarn) % end -% +% % SwE.xCon = xCon; % end % if nc>1... SwE.xCon = xCon; @@ -604,7 +604,7 @@ % Ask whether to do additional voxelwise or clusterwise inference. try infType = xSwE.infType; - catch + catch if isfield(SwE, 'WB') if isfield(SwE.WB, 'TFCE') infType = spm_input('inference type','+1','b','voxelwise|clusterwise|TFCE',[0,1,2],3); @@ -615,20 +615,20 @@ infType = spm_input('inference type','+1','b','voxelwise|clusterwise',[0,1],1); end end - + if isfield(SwE, 'WB') % Work out the original form of inference performed. This will tell us % which maps have already been generated. Most importantly, whether we % can do FWE p value clusterwise inference. if SwE.WB.voxelWise orig_infType = 'vox'; - elseif SwE.WB.clusterWise + elseif SwE.WB.clusterWise orig_infType = 'clus'; else orig_infType = 'tfce'; end end - + end %-Compute & store contrast parameters, contrast/ESS images, & SwE images @@ -645,8 +645,8 @@ STAT = xCon(Ic(1)).STAT; VspmSv = cat(1,xCon(Ic).Vspm); - - + + %-Check conjunctions - Must be same STAT w/ same df %-------------------------------------------------------------------------- if (nc > 1) && (any(diff(double(cat(1,xCon(Ic).STAT)))) || ... @@ -671,7 +671,7 @@ % We display the equivalent statistics. switch STAT - case 'T' + case 'T' STATstr = sprintf('%c','Z',str); case 'F' STATstr = sprintf('%c','X',str); @@ -692,7 +692,7 @@ else Z = min(Z, swe_data_read(xCon(i).Vspm, 'xyz', XYZ)); end - + end @@ -762,23 +762,23 @@ %-Height threshold - classical inference %-------------------------------------------------------------------------- if STAT ~= 'P' - + % Get the equivalent statistic switch STAT - + case 'T' - + eSTAT = 'Z'; - + case 'F' - + eSTAT = 'X'; - + end - + % If we are doing voxelwise inference on a parametric. if ~isfield(SwE, 'WB') && infType == 0 - + %-Get height threshold %---------------------------------------------------------------- fprintf('%s%30s',repmat(sprintf('\b'),1,30),'...height threshold') %-# @@ -790,7 +790,7 @@ thresDesc = spm_input('p value adjustment to control','+1','b',str,[],1); end - + switch thresDesc case 'FDR' % False discovery rate @@ -803,9 +803,9 @@ thresDesc = ['p<' num2str(u) ' (' thresDesc ')']; switch STAT case 'T' - u = swe_uc_FDR(u,Inf,'Z',n,VspmSv,0); + u = swe_uc_FDR(u,Inf,'Z',n,VspmSv,0); case 'F' - u = swe_uc_FDR(u,[1 1],'X',n,VspmSv,0); + u = swe_uc_FDR(u,[1 1],'X',n,VspmSv,0); end case 'none' % No adjustment: p for conjunctions is p of the conjunction SwE @@ -831,35 +831,35 @@ %-------------------------------------------------------------- fprintf('\n'); %-# error('Unknown control method "%s".',thresDesc); - + end % switch thresDesc - + %-Compute p-values for topological and voxel-wise FDR (all search voxels) %---------------------------------------------------------------- fprintf('%s%30s',repmat(sprintf('\b'),1,30),'...for voxelFDR') %-# switch STAT case 'T' - Ps = (1-spm_Ncdf(Zum)).^n; + Ps = (1-spm_Ncdf(Zum)).^n; case 'F' Ps = (1-spm_Xcdf(Zum,1)).^n; end - + up = NaN; Pp = NaN; uc = NaN; ue = NaN; Pc = []; uu = []; - + Q = find(Z > u); - + % If we are doing clusterwise inference on a parametric. elseif ~isfield(SwE, 'WB') && infType == 1 - + % Record what type of clusterwise inference we are doing. clustWise = 'Uncorr'; - + % No adjustment: p for conjunctions is p of the conjunction SwE %-------------------------------------------------------------- try @@ -878,29 +878,29 @@ else thresDesc = ''; end - + %-Compute p-values for topological and voxel-wise FDR (all search voxels) %---------------------------------------------------------------- fprintf('%s%30s',repmat(sprintf('\b'),1,30),'...for voxelFDR') %-# switch STAT case 'T' - Ps = (1-spm_Ncdf(Zum)).^n; + Ps = (1-spm_Ncdf(Zum)).^n; case 'F' Ps = (1-spm_Xcdf(Zum,1)).^n; end - + up = NaN; Pp = NaN; uc = NaN; ue = NaN; Pc = []; uu = []; - + Q = find(Z > u); % If we are doing voxelwise inference on a WB. elseif isfield(SwE, 'WB') && infType == 0 - + %-Get height threshold %---------------------------------------------------------------------- fprintf('%s%30s',repmat(sprintf('\b'),1,30),'...height threshold') %-# @@ -911,7 +911,7 @@ str = 'FWE|FDR|none'; thresDesc = spm_input('p value adjustment to control','+1','b',str,[],1); end - + switch thresDesc case 'FWE' % Family-wise false positive rate @@ -923,14 +923,14 @@ pu = spm_input('p value (FWE)','+0','r',0.05,1,[0,1]); end thresDesc = ['p<=' num2str(pu) ' (' thresDesc ')']; - + FWE_ps = 10.^-swe_data_read(xCon(Ic).VspmFWEP,'xyz', XYZ); - + % When thresholding on WB FWER p-values, we should include those = to pu % Here, we are using a - tol < b instead of a <= b due to numerical errors % tol was set to 0.1/(nB+1) in order to make sure it is smaller than the smallest WB p-value Q = find(FWE_ps - tol < pu); - + % Obtain the exclusive statistic threshold. This will be the (1-pu)th % percentile of the max. statistic distribution if Ic == 1 @@ -951,14 +951,14 @@ pu = spm_input('p value (FDR)','+0','r',0.05,1,[0,1]); end thresDesc = ['p<=' num2str(pu) ' (' thresDesc ')']; - + % select the WB FDR p-values within the mask FDR_ps = 10.^-swe_data_read(xCon(Ic).VspmFDRP, 'xyz', XYZ); % Here, a parametric score threshold u would differ from voxel to voxel % Thus, setting it to NaN u = NaN; - + % inclusive thresholding for WB Q = find(FDR_ps - tol < pu); @@ -977,7 +977,7 @@ % Here, a parametric score threshold u would differ from voxel to voxel % Thus, setting it to NaN u = NaN - + % inclusive thresholding for WB Q = find(unc_ps - tol < pu); @@ -986,7 +986,7 @@ fprintf('\n'); %-# error('Unknown control method "%s".',thresDesc); end - + up = NaN; Pp = NaN; uc = NaN; @@ -996,7 +996,7 @@ % If we are doing clusterwise WB. elseif isfield(SwE, 'WB') && infType == 1 - + fprintf('%s%30s',repmat(sprintf('\b'),1,30),'...height threshold') %-# try thresDesc = xSwE.thresDesc; @@ -1012,13 +1012,13 @@ thresDesc = 'none'; end end - + switch thresDesc - + case 'FWE' % Family-wise false positive rate % This is performed on the FWE P value map %-------------------------------------------------------- - + try clusterSizeType = xSwE.clusterSizeType; catch @@ -1033,9 +1033,9 @@ % Record what type of clusterwise inference we are doing. clustWise = 'FWE'; - - % For WB we have performed our thresholding and worked out of - % regions of activation. + + % For WB we have performed our thresholding and worked out of + % regions of activation. if Ic == 1 locActVox = SwE.WB.clusterInfo.LocActivatedVoxels; elseif Ic == 2 @@ -1045,7 +1045,7 @@ end [~, index]=ismember(XYZ',locActVox','rows'); Q=find(index~=0); - + % We also should record the cluster forming threshold that was % used. pu = SwE.WB.clusterInfo.primaryThreshold; @@ -1055,7 +1055,7 @@ else u = chi2inv(1-pu, 1); end - + % We should display the cluster forming threshold that % was used. try @@ -1064,13 +1064,13 @@ spm_input(['threshold {p value}'],... '+1','b',['(pre-set: P=' num2str(pu) ')'],[0],0) end - + case 'none' % No adjustment: p for conjunctions is p of the conjunction SwE % This should be performed on the uncorrected WB p-values %-------------------------------------------------------- % Record what type of clusterwise inference we are doing. clustWise = 'Uncorr'; - + % Cluster-forming threshold. try pu = xSwE.u; @@ -1080,39 +1080,39 @@ thresDesc = ['p<=' num2str(pu) ' (unc.)']; % select the WB unc. p-values within the mask unc_ps = 10.^-swe_data_read(xCon(Ic).VspmUncP, 'xyz', XYZ); - + % Here, a parametric score threshold u would differ from voxel to voxel % Thus, setting it to NaN u = NaN - + % inclusive thresholding for WB Q = find(unc_ps - tol < pu); - + up = NaN; Pp = NaN; uc = NaN; ue = NaN; Pc = []; uu = []; - + end - - % If we are doing TFCE. + + % If we are doing TFCE. else - + % Remember we are not doing clusterwise. clustWise = 'None'; - + % Ask user for TFCE FWE alpha. pt = spm_input('p value (TFCE FWE)','+0','r',0.05,1,[0,1]); thresDesc = ['p<=' num2str(pt) ' (FWE)']; - + % Get Tfce Fwe P-values. % In older version of the toolbox, the max TFCE scores were not saved. - % Thus, to avoid retro-compatibility issues, we cannot threshold using + % Thus, to avoid retro-compatibility issues, we cannot threshold using % the (1-pt)th percentile of the max distribution, but only using the FWER p-values tfp = 10.^-swe_data_read(xCon(Ic).VspmTFCEFWEP, 'xyz', XYZ); - + up = NaN; Pp = NaN; uc = NaN; @@ -1124,7 +1124,7 @@ % Here, we are using a - tol < b instead of a <= b due to numerical errors % tol was set to 0.1/(nB+1) in order to make sure it is smaller than the smallest WB p-value Q = find(tfp - tol < pt); - + end end @@ -1142,23 +1142,23 @@ strDataType = 'voxels/vertices'; elseif spm_mesh_detect(xCon(Ic(1)).Vspm) strDataType = 'vertices'; - else - strDataType = 'voxels'; + else + strDataType = 'voxels'; end end if isempty(Q) fprintf('\n'); warning('SwE:NoVoxels','No %s survive thresholding', strDataType); end - + % If we are doing clusterwise ask for threshold. if ~strcmp(clustWise, 'None') %-Extent threshold (disallowed for conjunctions) %-------------------------------------------------------------------------- - if ~isempty(XYZ) && nc == 1 + if ~isempty(XYZ) && nc == 1 fprintf('%s%30s',repmat(sprintf('\b'),1,30),'...extent threshold') %-# - + % Uncorrected threshold. if strcmp(clustWise, 'Uncorr') %-Get extent threshold [default = 0] @@ -1172,7 +1172,7 @@ %-Calculate extent threshold filtering %---------------------------------------------------------------------- if isCifti - A = swe_cifti_clusters(SwE.cifti, XYZ(1,:)); + A = swe_cifti_clusters(SwE.cifti, XYZ(1,:)); elseif ~spm_mesh_detect(xCon(Ic(1)).Vspm) A = spm_clusters(XYZ); else @@ -1187,7 +1187,7 @@ if length(j) >= k; Q = [Q j]; end end end - + % FWE WB threshold. if strcmp(clustWise, 'FWE') %-Get extent threshold [default = 0] @@ -1202,7 +1202,7 @@ %---------------------------------------------------------------------- % recompute the clusters as they may have been reduced due to post-hoc masking if isCifti - [A, clusterSizesInSurfaces, clusterSizesInVolume] = swe_cifti_clusters(SwE.cifti, XYZ(1,:)); + [A, clusterSizesInSurfaces, clusterSizesInVolume] = swe_cifti_clusters(SwE.cifti, XYZ(1,:)); elseif ~spm_mesh_detect(xCon(Ic(1)).Vspm) A = spm_clusters(XYZ); else @@ -1212,7 +1212,7 @@ A = A(XYZ(1,:)); end clusIndices = unique(A); - + % recompute the p-values as they might have increased due to post-hoc masking if Ic == 1 if isCifti && strcmpi(clusterSizeType, 'Box-Cox norm. k_{Z}') @@ -1269,17 +1269,17 @@ else error('Unknown contrast'); end - + ps_fwe = nan(1,numel(A)); for i = 1:length(clusIndices) ind = ( A==clusIndices(i) ); ps_fwe(ind) = sum( maxClusterSize > clusterSizes(i) - 1e-8) / (SwE.WB.nB + 1); end - % select only the voxels surviving the FWER threshold + % select only the voxels surviving the FWER threshold % Here, we use an inclusive p-value threshold to be consistent with an exclusive cluster threshold Q = find(ps_fwe - tol < fwep_c); - + % The exclusive threshold k should be the (1-fwep_c)th percentile of the max cluster size distribution maxClusterSize = sort(maxClusterSize); @@ -1295,20 +1295,20 @@ fprintf('\n'); %-# warning('SwE:NoVoxels','No %s survive cluster extent threshoding', strDataType); end - + else k = 0; - end + end end -end +end % Doftype if isfield(SwE.type,'modified') dof_type = SwE.type.modified.dof_mo; else - dof_type = SwE.type.classic.dof_cl; + dof_type = SwE.type.classic.dof_cl; end %========================================================================== @@ -1360,7 +1360,7 @@ 'nPredict', size(SwE.xX.X, 2),... 'df_Con', rank(xCon(Ic).c),... 'SS', SwE.SS); - + if isfield(SwE.type, 'modified') xSwE.nSubj_g = SwE.Gr.nSubj_g; xSwE.max_nVis_g = SwE.Vis.max_nVis_g; @@ -1370,27 +1370,27 @@ % For WB analyses we have already computed uncorrected, FDR, FWE and % cluster-FWE P values at this point. if isfield(SwE, 'WB') - + % Bootstrap variables. xSwE.WB = 1; xSwE.nB = SwE.WB.nB; - + % Volumes xSwE.VspmUncP = cat(1,xCon(Ic).VspmUncP); xSwE.VspmFDRP = cat(1,xCon(Ic).VspmFDRP); xSwE.VspmFWEP = cat(1,xCon(Ic).VspmFWEP); if SwE.WB.clusterWise - + try xSwE.VspmFWEP_clusnorm = cat(1,xCon(Ic).VspmFWEP_clusnorm); if strcmp(clustWise, 'FWE') xSwE.clusterSizeType = clusterSizeType; end - + hasVolumeData = (SwE.xY.dataType == swe_DataType('Cifti') && numel(SwE.cifti.surfaces) > 0) || SwE.xY.dataType == swe_DataType('Gifti') || SwE.xY.dataType == swe_DataType('SurfaceMat'); hasSurfaceData = (SwE.xY.dataType == swe_DataType('Cifti') && numel(SwE.cifti.volume) > 0) || SwE.xY.dataType == swe_DataType('Nifti') || SwE.xY.dataType == swe_DataType('VolumeMat'); - + if Ic == 1 if hasVolumeData && numel(SwE.WB.clusterInfo.clusterSizesInSurfacesUnderH0) > 0 xSwE.boxcoxInfo.surfaces.lambda = SwE.WB.clusterInfo.clusterSizesInSurfacesUnderH0_boxCox_lambda; @@ -1426,7 +1426,7 @@ catch xSwE.VspmFWEP_clus = cat(1,xCon(Ic).VspmFWEP_clus); - + end end if isfield(SwE.WB, 'TFCE') @@ -1443,13 +1443,13 @@ else xSwE.TFCEthresh = 0; end - + % Uncorrected P values. Ps_vol = xSwE.VspmUncP; Ps = swe_data_read(Ps_vol, 'xyz', SwE.xVol.XYZ); Ps = 10.^(-Ps(~isnan(Ps))); xSwE.Ps = Ps; - + % 95% percentiles if Ic == 1 maxScore = sort(SwE.WB.maxScore); @@ -1458,7 +1458,7 @@ else error('Unknown contrast'); end - xSwE.Pfv = maxScore(ceil(0.95*(xSwE.nB+1))); % Voxelwise FWE P + xSwE.Pfv = maxScore(ceil(0.95*(xSwE.nB+1))); % Voxelwise FWE P if SwE.WB.clusterWise if Ic == 1 maxClusterSize = sort(SwE.WB.clusterInfo.maxClusterSize); @@ -1467,9 +1467,9 @@ else error('Unknown contrast'); end - xSwE.Pfc = maxClusterSize(ceil(0.95*(xSwE.nB+1))); % Clusterwise FWE P + xSwE.Pfc = maxClusterSize(ceil(0.95*(xSwE.nB+1))); % Clusterwise FWE P end - + % edf xSwE.Vedf = cat(1,xCon(Ic).Vedf); else @@ -1477,7 +1477,7 @@ xSwE.edf = xCon(Ic).edf; else xSwE.Vedf = cat(1,xCon(Ic).Vedf); - end + end end % Record clusterwise FWE P value if there is one. @@ -1510,7 +1510,7 @@ xSwE.S_vol = 0; end end - + %-RESELS per voxel (density) if it exists %-------------------------------------------------------------------------- try, xSwE.VRpv = SwE.xVol.VRpv; end diff --git a/swe_get_defaults.m b/swe_get_defaults.m index faec0c7..71bad96 100644 --- a/swe_get_defaults.m +++ b/swe_get_defaults.m @@ -1,8 +1,8 @@ function varargout = swe_get_defaults(defstr, varargin) % Get/set the defaults values associated with an identifier % FORMAT defval = swe_get_defaults(defstr) -% Return the defaults value associated with identifier "defstr". -% Currently, this is a '.' subscript reference into the global +% Return the defaults value associated with identifier "defstr". +% Currently, this is a '.' subscript reference into the global % "defaults" variable defined in swe_defaults.m. % % FORMAT swe_get_defaults(defstr, defval) diff --git a/swe_get_file_extension.m b/swe_get_file_extension.m index 98696df..fc8220c 100644 --- a/swe_get_file_extension.m +++ b/swe_get_file_extension.m @@ -21,7 +21,7 @@ % ========================================================================= % Bryan Guillaume % Version Info: $Format:%ci$ $Format:%h$ - + if nargin == 1 lastDot = false; elseif nargin == 2 @@ -29,10 +29,10 @@ else error('too many inputs'); end - + indexDots = find(filename == '.'); nDots = numel(indexDots); - + if nDots == 0 error('No extension found!'); else @@ -54,4 +54,4 @@ fileExtension = tmp; end end -end \ No newline at end of file +end diff --git a/swe_invNcdf.m b/swe_invNcdf.m index 449c0e4..7d4c0b7 100644 --- a/swe_invNcdf.m +++ b/swe_invNcdf.m @@ -17,8 +17,8 @@ var = 1; end if nargin<2 - mu = 0; + mu = 0; end - z = -sqrt(2).* (sqrt(var) .* erfcinv(2*p)) + mu; - + z = -sqrt(2).* (sqrt(var) .* erfcinv(2*p)) + mu; + end diff --git a/swe_jobman.m b/swe_jobman.m index ba0e5ab..3a2cc2c 100644 --- a/swe_jobman.m +++ b/swe_jobman.m @@ -1,6 +1,6 @@ function varargout = swe_jobman(varargin) % Main interface for SPM Batch System - % Copy of spm_jobman using swe_select instead of spm_select + % Copy of spm_jobman using swe_select instead of spm_select % This function provides a compatibility layer between SPM and matlabbatch. % % FORMAT swe_jobman('initcfg') @@ -17,7 +17,7 @@ % input1,... - optional list of input arguments. These are filled into % open inputs ('X->' marked in the GUI) before a job is % run. When using an "{:}" subscript on a cell array, - % MATLAB expands this cell array into a comma separated + % MATLAB expands this cell array into a comma separated % list of arguments. Therefore, one can collect input % arguments in the right order into a cell array named e.g. % input_array and call swe_jobman('run',job,input_array{:}) @@ -57,8 +57,8 @@ % Copyright (C) 2005-2016 Wellcome Trust Centre for Neuroimaging % Volkmar Glauche % Modified by Bryan Guillaume - % Version Info: $Format:%ci$ $Format:%h$ - + % Version Info: $Format:%ci$ $Format:%h$ + %__________________________________________________________________________ % % Programmers help: @@ -84,21 +84,21 @@ % FORMAT [tag, job] = swe_jobman('harvest', job_id|job|cfg_item|cfg_struct) % Take the job with id job_id in cfg_util and extract what is % needed to save it as a batch job (for experts only). If a (partial) job - % is given instead, the output job is augmented with default settings. + % is given instead, the output job is augmented with default settings. % If the argument is a cfg_item or cfg_struct tree, it will be harvested - % outside cfg_util. + % outside cfg_util. % tag - tag of the root node of the current job/cfg_item tree % job - harvested data from the current job/cfg_item tree %__________________________________________________________________________ - % - % This code is based on earlier versions by John Ashburner, Philippe + % + % This code is based on earlier versions by John Ashburner, Philippe % Ciuciu and Guillaume Flandin. % It now relies on matlabbatch % http://sourceforge.net/projects/matlabbatch/ % Copyright (C) 2008 Freiburg Brain Imaging %__________________________________________________________________________ - - + + %-Force jobs configuration initialisation if needed %-------------------------------------------------------------------------- persistent isInitCfg; @@ -109,7 +109,7 @@ fprintf('done.\n'); end isInitCfg = true; - + %-Open GUI when called without input arguments %-------------------------------------------------------------------------- if ~nargin @@ -117,7 +117,7 @@ if nargout > 0, varargout = {findobj(0,'tag','cfg_ui')}; end return; end - + %-Warn about deprecated syntax %-------------------------------------------------------------------------- action = lower(varargin{1}); @@ -131,7 +131,7 @@ 'Callback ''%s'' is deprecated. Use ''convert'' instead.',action); action = 'convert'; end - + %-Load and convert batch jobs %-------------------------------------------------------------------------- if ismember(action, {'interactive','run','serial'}) @@ -160,7 +160,7 @@ end end end - + %-Perform action %-------------------------------------------------------------------------- switch action @@ -193,7 +193,7 @@ % f2 = uimenu(f0,'Label','xxx', 'Callback',@xxx, ... % 'HandleVisibility','off', 'tag','xxx'); %end - + case {'interactive'} if exist('mljob', 'var') cjob = cfg_util('initjob', mljob); @@ -217,7 +217,7 @@ if nargout > 0 varargout{1} = cjob; end - + case {'serial'} if exist('mljob', 'var') cjob = cfg_util('initjob', mljob); @@ -238,7 +238,7 @@ end cfg_util('deljob', cjob); end - + case {'run'} if ~exist('mljob', 'var') error('Not enough input arguments.'); @@ -254,10 +254,10 @@ end cfg_util('deljob', cjob); end - + case {'convert'} varargout{1} = convert_jobs(varargin{2:end}); - + case {'harvest'} if nargin == 1 error('spm:swe_jobman:CantHarvest', ... @@ -281,7 +281,7 @@ end varargout{1} = tag; varargout{2} = job; - + case {'help'} if (nargin < 2) || isempty(varargin{2}) node = 'spm'; @@ -294,12 +294,12 @@ width = varargin{3}; end varargout{1} = cfg_util('showdocwidth', width, node); - + otherwise error('spm:swe_jobman:unknownOption','Unknown option "%s".',varargin{1}); end - - + + %========================================================================== % function newjobs = load_jobs(job) %========================================================================== @@ -355,8 +355,8 @@ newjobs = [newjobs(:); {[]}]; end end - - + + %========================================================================== % function varargout = convert_jobs(varargin) %========================================================================== @@ -386,8 +386,8 @@ end end varargout = {outnames}; - - + + %========================================================================== % function [mljob, comp] = canonicalise_jobs(job) %========================================================================== @@ -417,8 +417,8 @@ mljob{i} = job{i}; end end - - + + %========================================================================== % function njobs = sub_canonicalise_job(jobs) %========================================================================== @@ -453,8 +453,8 @@ njobs{end+1} = jobs{i0}; end end - - + + %========================================================================== % function sts = fill_run_job(action, cjob, varargin) %========================================================================== @@ -471,8 +471,8 @@ cfg_util('deljob', cjob); error('spm:swe_jobman:jobNotFilled', 'No executable modules, but still unresolved dependencies or incomplete module inputs.'); end - - + + %========================================================================== % function [val sts] = serial_ui(item) %========================================================================== @@ -507,7 +507,7 @@ % enter more (up to varargin{3}(2) values labels = {labels{:} 'Done'}; % values is a cell list of natural numbers, use -1 for Done - values = {values{:} -1}; + values = {values{:} -1}; while numel(val) < item.num(2) val1 = spm_input(sprintf('%s(%d)', item.name, numel(val)+1), 1, ... 'm', labels, values); @@ -530,8 +530,8 @@ error('File selector was closed.'); end end - - + + %========================================================================== % function [code cont] = genscript_run %========================================================================== @@ -542,10 +542,9 @@ code{1} = sprintf('spm(''defaults'', ''%s'');', modality); code{2} = 'swe_jobman(''run'', jobs, inputs{:});'; cont = false; - - + + %-Compatibility layer for SPM5 function varargout = interactive(varargin) function varargout = defaults_edit(varargin) function varargout = run_serial(varargin) - \ No newline at end of file diff --git a/swe_list.m b/swe_list.m index c5981ef..6de75fc 100644 --- a/swe_list.m +++ b/swe_list.m @@ -8,7 +8,7 @@ % Num - number of maxima per cluster [3] % Dis - distance among clusters {mm} [8] % Str - header string -% +% % FORMAT TabDat = swe_list('ListCluster',xSwE,hReg,[Num,Dis,Str]) % List of local maxima for a single suprathreshold cluster % xSwE - Structure containing data (format as below) @@ -45,7 +45,7 @@ % .M - voxels - > mm matrix % .VOX - voxel dimensions {mm} % .DIM - image dimensions {voxels} -% .units - space units +% .units - space units % .VRpv - filehandle - Resels per voxel % .Ps - uncorrected P values in searched volume (for voxel FDR) % .Pp - uncorrected P values of peaks (for peak FDR) @@ -91,9 +91,9 @@ % excursion sets (a collection of face, edge and vertex connected subsets % or clusters). The corrected significance of the results are based on % set, cluster and voxel-level inferences using distributional -% approximations from the Theory of Gaussian Fields and/or common bootstrap -% methods. These distributions assume that the SwE map is a reasonable -% lattice approximation of a continuous random field with known component +% approximations from the Theory of Gaussian Fields and/or common bootstrap +% methods. These distributions assume that the SwE map is a reasonable +% lattice approximation of a continuous random field with known component % field smoothness. % % The p values are based on the probability of obtaining c, or more, @@ -155,40 +155,40 @@ if nargin < 2, error('Not enough input arguments.'); end if nargin < 3, hReg = []; else hReg = varargin{3}; end xSwE = varargin{2}; - + %-Extract results table and display it %---------------------------------------------------------------------- spm('Pointer','Watch') - + TabDat = swe_list('Table',xSwE,varargin{4:end}); - + swe_list('Display',TabDat,hReg); - + spm('Pointer','Arrow') - + %-Return TabDat structure %---------------------------------------------------------------------- varargout = { TabDat }; - - + + %========================================================================== case 'table' %-Table %========================================================================== % FORMAT TabDat = swe_list('table',xSwE,[Num,Dis,Str]) - + %-Parse arguments %---------------------------------------------------------------------- if nargin < 2, error('Not enough input arguments.'); end xSwE = varargin{2}; - + %-Get number of maxima per cluster to be reported %---------------------------------------------------------------------- if length(varargin) > 2, Num = varargin{3}; else Num = spm_get_defaults('stats.results.volume.nbmax'); end - + %-Get minimum distance among clusters (mm) to be reported %---------------------------------------------------------------------- if length(varargin) > 3, Dis = varargin{4}; else Dis = spm_get_defaults('stats.results.volume.distmin'); end - + %-Get header string %---------------------------------------------------------------------- if length(varargin) > 4 && ~isempty(varargin{5}) @@ -200,7 +200,7 @@ Title = 'Posterior Probabilities'; end end - + %-Extract data from xSwE %---------------------------------------------------------------------- isCifti = xSwE.isCifti; @@ -225,7 +225,7 @@ try, QPs = xSwE.Ps; end try, QPp = xSwE.Pp; end try, QPc = xSwE.Pc; end - + % For WB analyses we have already calculated the information for the % table and footer. We just need to read it in. if xSwE.WB @@ -245,7 +245,7 @@ VspmFWEP_clusnorm = []; end end - + % if STAT~='P' % R = full(xSwE.R); % FWHM = full(xSwE.FWHM); @@ -257,7 +257,7 @@ end units{1} = [units{1} ' ']; units{2} = [units{2} ' ']; - + if ~spm_mesh_detect(xSwE.Vspm) DIM = DIM > 1; % non-empty dimensions strDataType = 'voxels'; @@ -267,7 +267,7 @@ end D = sum(DIM); % highest dimension VOX = VOX(DIM); % scaling - + % if STAT ~= 'P' % FWHM = FWHM(DIM); % Full width at half max % FWmm = FWHM.*VOX; % FWHM {units} @@ -281,15 +281,15 @@ % Choose between voxel-wise and topological FDR %---------------------------------------------------------------------- topoFDR = false; %to be checked - + %-Tolerance for p-value underflow, when computing equivalent Z's %---------------------------------------------------------------------- tol = eps*10; - + %-Table Headers %---------------------------------------------------------------------- TabDat.tit = Title; - + % If we are doing a clusterwise/voxelwise analysis the header is the % normal SPM header. if isCifti @@ -299,7 +299,7 @@ additionalField = {}; nColTable = 11; end - + if ~isfield(xSwE, 'TFCEanaly') || ~xSwE.TFCEanaly TabDat.hdr = {... 'set', 'p', '\itp';... @@ -331,9 +331,9 @@ 'peak', STATe, sprintf('\\it%c',STATe);... 'peak', 'p(unc)', '\itp\rm_{uncorr}';... ' ', 'x,y,z {mm}', [units{:}];... - additionalField{:} }'; + additionalField{:} }'; end - + %-Coordinate Precisions %---------------------------------------------------------------------- if isempty(xSwE.XYZmm) || isCifti % empty results or cifti @@ -365,7 +365,7 @@ '%0.3f', '%0.0f','%0.2f','%0.3f',... %-Cluster '%0.3f', '%0.3f', '%6.2f', '%0.3f',... %-Peak xyzfmt, '%s'}; %-XYZ - + %-Table filtering note %---------------------------------------------------------------------- if isinf(Num) @@ -373,8 +373,8 @@ else TabDat.str = sprintf(['table shows %d local maxima ',... 'more than %.1fmm apart'],Num,Dis); - end - + end + %-Footnote with SPM parameters %---------------------------------------------------------------------- if strcmp(STAT, 'T') @@ -384,10 +384,10 @@ Pz = 1-spm_Xcdf(u, 1); eSTAT = 'X'; end - + % Create footer for display. TabDat.ftr = cell(6,2); - + % Number of `extra` lines inserted that don't have to be present in % every display exlns = 0; @@ -446,11 +446,11 @@ Ts(isnan(Ts)) = []; Ts = 10.^-Ts; Ts = sort(Ts(:)); - + % Obtain the FDR p 0.05 value. FDRp_05 = swe_uc_FDR(0.05,Inf,'P',n,Ts); clear Ts - + % Record FWE/FDR/clus FWE p values. (No clus FWE for voxelwise and % TFCE analyses) if xSwE.infType == 1 && strcmp(xSwE.clustWise, 'FWE') @@ -466,14 +466,14 @@ % Record height thresholds. TabDat.ftr{1,1} = ... ['Threshold: Height ' eSTAT ' = %0.2f, p = %0.3f; Extent k = %0.0f ' strDataType '.']; - TabDat.ftr{1,2} = [u,Pz,k]; + TabDat.ftr{1,2} = [u,Pz,k]; % Record FDR p value. TabDat.ftr{2,1} = ... 'vox P(5%% FDR): %0.3f'; TabDat.ftr{2,2} = swe_uc_FDR(0.05,Inf,'P',n,sort(xSwE.Ps)'); - + end - + if xSwE.infType == 1 && strcmp(xSwE.clustWise, 'FWE') && isfield(xSwE, 'boxcoxInfo') TabDat.ftr{(3+exlns),1} = 'k_{Z} = 0.6745 (k_{\\lambda} - Q2(k_{\\lambda}^{H0})) / (Q3(k_{\\lambda}^{H0})-Q2(k_{\\lambda}^{H0}))'; % TabDat.ftr{(3+exlns),1} = 'Null cluster sizes in surfaces: \\lambda_S=%0.2f , \\lambda_V =%0.2f'; @@ -495,7 +495,7 @@ else error('Unknown Box-Cox Info!') end - exlns = exlns + 2; + exlns = exlns + 2; end % If we have groups display group details in ftr. @@ -523,10 +523,10 @@ end TabDat.ftr{3+exlns,1} = [nSubjString nVisitsString]; TabDat.ftr{3+exlns,2} = [xSwE.nSubj_g nVisitsNumbers]; - + exlns = exlns + 1; end - + % Retrieve edf data if isfield(xSwE, 'Vedf') edf = swe_data_read(xSwE.Vedf, 'xyz', xSwE.XYZ_inMask); @@ -534,21 +534,21 @@ edf = xSwE.edf; end edf(isnan(edf)) = []; - + edf_max = max(edf); edf_min = min(edf); edf_med = median(edf); - + % Work out range of dof values diff = abs(edf_max - edf_min); - + % Work out dofType switch xSwE.dofType - - case 0 - + + case 0 + dofTypeStr = 'Naive'; - + case 1 dofTypeStr = 'Approx I'; @@ -566,7 +566,7 @@ error('Unknown degrees of freedom.') end - + % Recording effective Degrees of freedom if xSwE.dofType~=0 && diff > 10^-10 TabDat.ftr{(3+exlns),1}=['Error DF: (' dofTypeStr '): (min) %0.1f, (median) %0.1f, (max) %0.1f']; @@ -575,7 +575,7 @@ TabDat.ftr{(3+exlns),1}=['Error DF: (' dofTypeStr '): %0.1f']; TabDat.ftr{(3+exlns),2}=edf_med; end - + % Record small sample adjustments. TabDat.ftr{(4+exlns),1}='Resid. Adj.: %s'; switch xSwE.SS @@ -589,7 +589,7 @@ % Record contrast degrees of freedom. TabDat.ftr{(5+exlns),1} = 'Contrast DF: %0.0f; Number of predictors: %0.0f'; TabDat.ftr{(5+exlns),2} = [xSwE.df_Con xSwE.nPredict]; - + % Record volume. if isCifti TabDat.ftr{(6+exlns),1} = ... @@ -612,27 +612,27 @@ TabDat.ftr{(7+exlns),2} = sqrt(diag(xSwE.cifti.volume.M(1:3,1:3)'*xSwE.cifti.volume.M(1:3,1:3)))'; elseif isCifti && numel(xSwE.cifti.volume) == 0 exlns = exlns - 1; - elseif ~spm_mesh_detect(xSwE.Vspm) + elseif ~spm_mesh_detect(xSwE.Vspm) TabDat.ftr{(7+exlns),1} = ... ['Voxel size: ' voxfmt units{:}]; TabDat.ftr{(7+exlns),2} = VOX; else exlns = exlns - 1; end - + if isfield(xSwE, 'TFCEanaly') && xSwE.TFCEanaly TabDat.ftr{(8+exlns),1} = 'TFCE: E=%0.1f, H=%0.1f'; TabDat.ftr{(8+exlns),2} = [xSwE.TFCE.E, xSwE.TFCE.H]; exlns = exlns + 1; end if xSwE.WB - + % Recording number of bootstraps. TabDat.ftr{(8+exlns),1}='Bootstrap samples = %0.0f'; TabDat.ftr{(8+exlns),2}= xSwE.nB; - + exlns = exlns + 1; - + end %-Characterize excursion set in terms of maxima % (sorted on Z values and grouped by regions) @@ -669,7 +669,7 @@ %---------------------------------------------------------------------- c = max(A); %-Number of clusters NONSTAT = spm_get_defaults('stats.rft.nonstat'); - + %-Convert maxima locations from voxels to mm %---------------------------------------------------------------------- if isCifti @@ -690,11 +690,11 @@ TabDat.dat = {Pc,c}; TabLin = 1; - + %-Cluster and local maxima p-values & statistics %---------------------------------------------------------------------- while numel(find(isfinite(Z))) - + %-Find largest remaining local maximum %------------------------------------------------------------------ [U,i] = max(Z); %-largest maxima @@ -711,7 +711,7 @@ case 'X' try Pz = 1-chi2cdf(U,1); - catch + catch Pz = 1-spm_Xcdf(U,1); end end @@ -719,7 +719,7 @@ Pz = 10.^-VspmUncP(XYZ(1,i),XYZ(2,i),XYZ(3,i)); end - % If we are not running a wild bootstrap or we are doing a + % If we are not running a wild bootstrap or we are doing a % small volume correction we need to calculate the FDR P value % and leave the other values blank. if ~xSwE.WB || isfield(xSwE,'svc') @@ -741,7 +741,7 @@ Qp = []; ws = warning('off','SPM:outOfRangeNormal'); warning(ws); - + if xSwE.infType == 1 && strcmp(xSwE.clustWise, 'FWE') % only for FWER clusterwise WB if strcmpi(xSwE.clusterSizeType, 'Box-Cox norm. k_{Z}') Pk = 10.^-VspmFWEP_clusnorm(XYZ(1,i),XYZ(2,i),XYZ(3,i)); @@ -756,30 +756,30 @@ % These regions will have NaN for the cluster FWE P-value % when they should have one. So the below is necessary: Pk(isnan(Pk)) = 1; - + elseif xSwE.TFCEanaly - + % Get coordinates of all voxels in the current cluster. currentClus = find(A == A(i)); XYZ_clus = XYZ(:, currentClus); - + % Read in all TFCE FWE P values in this cluser tfp = 10.^-swe_data_read(xSwE.VspmTFCEFWEP, 'xyz', XYZ_clus); - + % Record the minimum TFCE FWE P value in said cluster. Pk = min(tfp); else Pk = []; end - + end - + if i > numel(N_area) % means that this is for volume or there is no area info N_area_tmp = []; else N_area_tmp = N_area(i); end - if i > numel(N_boxcox) + if i > numel(N_boxcox) N_boxcox_tmp = []; else N_boxcox_tmp = N_boxcox(i); @@ -790,7 +790,7 @@ [TabDat.dat{TabLin, 12}] = char(brainStructureShortLabels(i)); end TabLin = TabLin + 1; - + %-Print Num secondary maxima (> Dis mm apart) %------------------------------------------------------------------ [l,q] = sort(-Z(mj)); % sort on Z value @@ -814,11 +814,11 @@ % if Pz < tol % Ze = Inf; % else -% Ze = swe_invNcdf(1 - Pz); +% Ze = swe_invNcdf(1 - Pz); % end % else - % If we are not running a wild bootstrap or if we are + % If we are not running a wild bootstrap or if we are % doing a small volume correction we need to calculate % the FDR P value and leave the other values blank. if ~xSwE.WB || isfield(xSwE,'svc') @@ -835,16 +835,16 @@ % If we are running a wild bootstrap we only need to read in % results we calculated earlier. else - + Pz = 10.^-VspmUncP(XYZ(1,d),XYZ(2,d),XYZ(3,d)); Pu = 10.^-VspmFWEP(XYZ(1,d),XYZ(2,d),XYZ(3,d)); Qu = 10.^-VspmFDRP(XYZ(1,d),XYZ(2,d),XYZ(3,d)); ws = warning('off','SPM:outOfRangeNormal'); Ze = swe_invNcdf(Z(d)); warning(ws); - + end - + D = [D d]; if topoFDR [TabDat.dat{TabLin,7:11}] = ... @@ -859,20 +859,20 @@ end Z(mj) = NaN; % Set local maxima to NaN end - + varargout = {TabDat}; - + %====================================================================== case 'display' %-Display table in Graphics window %====================================================================== % FORMAT swe_list('display',TabDat,hReg) - + %-Parse arguments %---------------------------------------------------------------------- - if nargin < 2, error('Not enough input arguments.'); + if nargin < 2, error('Not enough input arguments.'); else TabDat = varargin{2}; end if nargin < 3, hReg = []; else hReg = varargin{3}; end - + isCifti = (size(TabDat.hdr,2) == 12); if isCifti scalingFactor = 0.9; @@ -885,7 +885,7 @@ %-Get current location (to highlight selected voxel in table) %---------------------------------------------------------------------- xyzmm = swe_results_ui('GetCoords'); - + %-Setup Graphics panel %---------------------------------------------------------------------- Fgraph = spm_figure('FindWin','Satellite'); @@ -897,10 +897,10 @@ ht = 0.4; bot = 0.1; end swe_results_ui('Clear',Fgraph) - + FS = spm('FontSizes'); %-Scaled font sizes PF = spm_platform('fonts'); %-Font names (for this platform) - + %-Table axes & Title %---------------------------------------------------------------------- hAx = axes('Parent',Fgraph,... @@ -916,7 +916,7 @@ hRotate3d = rotate3d(Fgraph); setAllowAxesRotate(hRotate3d, hAx, false); end - + AxPos = get(hAx,'Position'); set(hAx,'YLim',[0,AxPos(4)]) dy = FS(9); y = floor(AxPos(4)) - dy; @@ -934,23 +934,23 @@ h = line([0,0.11]*scalingFactor,[1,1]*(y-dy/4),'LineWidth',0.5,'Color','r'); h = text(0.02*scalingFactor,y-9*dy/8, TabDat.hdr{3,1}); Hs = [Hs,h]; h = text(0.08*scalingFactor,y-9*dy/8, TabDat.hdr{3,2}); Hs = [Hs,h]; - + text(0.22*scalingFactor,y, [TabDat.hdr{1,3} '-level'],'FontSize',FS(9)); line([0.14,0.44]*scalingFactor,[1,1]*(y-dy/4),'LineWidth',0.5,'Color','r'); h = text(0.15*scalingFactor,y-9*dy/8, TabDat.hdr{3,3}); Hc = [Hc,h]; h = text(0.24*scalingFactor,y-9*dy/8, TabDat.hdr{3,4}); Hc = [Hc,h]; h = text(0.31*scalingFactor,y-9*dy/8, TabDat.hdr{3,5}); Hc = [Hc,h]; h = text(0.39*scalingFactor,y-9*dy/8, TabDat.hdr{3,6}); Hc = [Hc,h]; - + text(0.59*scalingFactor,y, [TabDat.hdr{1,8} '-level'],'FontSize',FS(9)); line([0.48,0.80]*scalingFactor,[1,1]*(y-dy/4),'LineWidth',0.5,'Color','r'); h = text(0.49*scalingFactor,y-9*dy/8, TabDat.hdr{3,7}); Hp = [Hp,h]; h = text(0.58*scalingFactor,y-9*dy/8, TabDat.hdr{3,8}); Hp = [Hp,h]; h = text(0.67*scalingFactor,y-9*dy/8, TabDat.hdr{3,9}); Hp = [Hp,h]; h = text(0.74*scalingFactor,y-9*dy/8, TabDat.hdr{3,10}); Hp = [Hp,h]; - + text(0.845*scalingFactor,y - dy/2,TabDat.hdr{3,11},'Fontsize',FS(8)); - + if isCifti text(0.88,y - dy/2,TabDat.hdr{1,12},'Fontsize',FS(8)); end @@ -973,7 +973,7 @@ if ~isempty(TabDat.ftr) set(gca,'DefaultTextFontName',PF.helvetica,... 'DefaultTextInterpreter','Tex','DefaultTextFontSize',FS(8)) - + fx = repmat([0 0.645],ceil(size(TabDat.ftr,1)/2),1); fy = repmat((1:ceil(size(TabDat.ftr,1)/2))',1,2); for i=1:size(TabDat.ftr,1) @@ -982,7 +982,7 @@ 'ButtonDownFcn','get(gcbo,''UserData'')'); end end - + %-Characterize excursion set in terms of maxima % (sorted on Z values and grouped by regions) %====================================================================== @@ -993,7 +993,7 @@ 'FontSize',FS(16),'Color',[1,1,1]*.5); return end - + %-Table proper %====================================================================== @@ -1004,7 +1004,7 @@ 0.49 0.58 0.65 0.74 ... %-Peak 0.84 ... %-XYZ 0.88/scalingFactor ] * scalingFactor; %-Brain structure labels - + %-Pagination variables %---------------------------------------------------------------------- hPage = []; @@ -1015,7 +1015,7 @@ if isempty(TabDat.dat{1,1}) && isempty(TabDat.dat{1,2}) % Pc set(Hs,'Visible','off'); end - + if TabDat.dat{1,2} >= 1 % c h = text(tCol(1),y,sprintf(TabDat.fmt{1},TabDat.dat{1,1}),... 'FontWeight','Bold', 'UserData',TabDat.dat{1,1},... @@ -1028,12 +1028,12 @@ else set(Hs,'Visible','off'); end - + %-Cluster and local maxima p-values & statistics %---------------------------------------------------------------------- HlistXYZ = []; for i=1:size(TabDat.dat,1) - + %-Paginate if necessary %------------------------------------------------------------------ if y < dy @@ -1045,11 +1045,11 @@ hPage = []; y = y0; end - + %-Print cluster and maximum peak-level p values %------------------------------------------------------------------ if ~isempty(TabDat.dat{i,4}), fw = 'Bold'; else fw = 'Normal'; end - + for k=3:nColTable if k < 11 h = text(tCol(k),y,sprintf(TabDat.fmt{k},TabDat.dat{i,k}),... @@ -1074,7 +1074,7 @@ hPage = [hPage, h]; end end - + % Specifically changed so it properly finds hMIPax %------------------------------------------------------------------ tXYZmm = TabDat.dat{i,11}; @@ -1098,7 +1098,7 @@ y = y - dy; end - + %-Number and register last page (if paginated) %---------------------------------------------------------------------- if spm_figure('#page',Fgraph)>1 @@ -1106,7 +1106,7 @@ % 'FontName',PF.helvetica,'FontSize',FS(8),'FontAngle','Italic'); spm_figure('NewPage',hPage) end - + %-End: Store TabDat in UserData of context menu %====================================================================== h = uicontextmenu('Tag','TabDat','UserData',TabDat); @@ -1137,7 +1137,7 @@ spm_XYZreg('Add2Reg',hReg,hAx,'swe_list'); varargout = {}; - + %====================================================================== case 'listcluster' %-List for current cluster only %====================================================================== @@ -1148,7 +1148,7 @@ if nargin < 2, error('Not enough input arguments.'); end if nargin < 3, hReg = []; else hReg = varargin{3}; end xSwE = varargin{2}; - + if isfield(xSwE,'G') warning('"current cluster" option not implemented for meshes.'); varargout = { evalin('base','TabDat') }; @@ -1158,7 +1158,7 @@ %-Get number of maxima per cluster to be reported %------------------------------------------------------------------ if nargin < 4, Num = spm_get_defaults('stats.results.svc.nbmax'); else Num = varargin{4}; end - + %-Get minimum distance among clusters (mm) to be reported %------------------------------------------------------------------ if nargin < 5, Dis = spm_get_defaults('stats.results.svc.distmin'); else Dis = varargin{5}; end @@ -1166,7 +1166,7 @@ %-Get header string %------------------------------------------------------------------ if nargin < 6, Str = ''; else Str = varargin{6}; end - + %-If there are suprathreshold voxels, filter out all but current cluster %------------------------------------------------------------------ if ~isempty(xSwE.Z) @@ -1231,7 +1231,7 @@ end fprintf('%c',repmat('=',1,80)), fprintf('\n\n') - + %====================================================================== case 'xlslist' %-Export table to Excel %====================================================================== @@ -1239,7 +1239,7 @@ if nargin<2, error('Not enough input arguments.'); end TabDat = varargin{2}; - + d = [TabDat.hdr(1:2,:);TabDat.dat]; xyz = d(3:end,end); xyz = num2cell([xyz{:}]'); @@ -1249,7 +1249,7 @@ tmpfile = [tempname '.xls']; xlswrite(tmpfile, d); winopen(tmpfile); - + %====================================================================== case 'csvlist' %-Export table to comma-separated values file %====================================================================== @@ -1257,7 +1257,7 @@ if nargin<2, error('Not enough input arguments.'); end TabDat = varargin{2}; - + tmpfile = [tempname '.csv']; fid = fopen(tmpfile,'wt'); fprintf(fid,[repmat('%s,',1,11) '%d,,\n'],TabDat.hdr{1,:}); @@ -1270,7 +1270,7 @@ end fclose(fid); open(tmpfile); - + %====================================================================== case 'setcoords' %-Coordinate change %====================================================================== diff --git a/swe_max.m b/swe_max.m index 70a124a..f7b3241 100644 --- a/swe_max.m +++ b/swe_max.m @@ -16,7 +16,7 @@ % ========================================================================= % Bryan Guillaume % Version Info: $Format:%ci$ $Format:%h$ - + [N, Z, M, A, XYZ] = spm_max(X, locationsInVoxels); try diff --git a/swe_mesh_clusters.m b/swe_mesh_clusters.m index 088defe..6b9b1bd 100644 --- a/swe_mesh_clusters.m +++ b/swe_mesh_clusters.m @@ -5,7 +5,7 @@ % activatedVertices - a [nVertices x 1] data vector (using NaNs or logicals) % areaFile - an optional file containing the surface area of each vertex % - % clusterAssignment - a [1 x nActivatedVertices] vector of cluster indices + % clusterAssignment - a [1 x nActivatedVertices] vector of cluster indices % clusterNumbersOfVertices - a [1 x nClusters] size of clusters {in vertices} % clusterAreas - a [1 x nClusters] vecot of cluster areas % ========================================================================= @@ -15,7 +15,7 @@ [clusterAssignment, clusterNumbersOfVertices] = spm_mesh_clusters(M, activatedVertices); clusterAssignment = clusterAssignment(activatedVertices)'; clusterNumbersOfVertices = clusterNumbersOfVertices'; - + if nargin > 2 && ~isempty(areaFile) area = swe_data_read(areaFile, find(activatedVertices)); nClusters = numel(clusterExtents); @@ -27,4 +27,4 @@ clusterAreas = NaN; end -end \ No newline at end of file +end diff --git a/swe_mesh_max.m b/swe_mesh_max.m index c409e80..7455480 100644 --- a/swe_mesh_max.m +++ b/swe_mesh_max.m @@ -20,7 +20,7 @@ % Version Info: $Format:%ci$ $Format:%h$ [N, Z, M, A, XYZ] = spm_mesh_max(X, locationsInVertices, G); - + canComputeArea = nargin > 4 && ~isempty(areaFile); if canComputeArea N_area = zeros(numel(N),1); @@ -43,4 +43,4 @@ N_boxcox = []; end -end \ No newline at end of file +end diff --git a/swe_progress_bar.m b/swe_progress_bar.m index aaf8251..2d698f8 100644 --- a/swe_progress_bar.m +++ b/swe_progress_bar.m @@ -1,9 +1,9 @@ function swe_progress_bar(varargin) % Opens the progress bar, if not in octave. % ========================================================================= -% This function was added for testing purposes and prevents the -% spm_progress_bar displaying in octave. The progress bar can be used in -% octave but it creates ascii art which spams the test logs. +% This function was added for testing purposes and prevents the +% spm_progress_bar displaying in octave. The progress bar can be used in +% octave but it creates ascii art which spams the test logs. % ========================================================================= % FORMAT swe_progress_bar('Set',value) % Set the height of the bar itself. @@ -26,4 +26,4 @@ function swe_progress_bar(varargin) spm_progress_bar(varargin{:}); end -end \ No newline at end of file +end diff --git a/swe_read_cifti_info.m b/swe_read_cifti_info.m index d8edd70..21f974d 100644 --- a/swe_read_cifti_info.m +++ b/swe_read_cifti_info.m @@ -11,7 +11,7 @@ % - volumes: information from each volumetric brain structure % Bryan Guillaume % Version Info: $Format:%ci$ $Format:%h$ - + % fetch the cifti information ciftiObject = swe_cifti(filename); xml = char(ciftiObject.hdr.ext.edata(:)'); @@ -19,7 +19,7 @@ root_uid = root(tree); hdr = convert(tree); hdr.attributes = getAttributes(tree, root_uid); - + % Check if input is a CIFTI-2 file switch hdr.attributes.Version case {'2'} @@ -31,7 +31,7 @@ warning('Unknown CIFTI version: %s.', hdr.attributes.Version); return; end - + % extract surfaces and a single volume (by concatenating all brain volumes) surfaces = {}; volumes = {}; @@ -91,10 +91,10 @@ error('Data have to be dense scalars or series.'); end end -end +end %========================================================================== function attrb = getAttributes(tree, uid) attrb = attributes(tree, 'get', uid); if ~isstruct(attrb), attrb = [attrb{:}]; end attrb = cell2struct({attrb.val},{attrb.key},2); -end \ No newline at end of file +end diff --git a/swe_regions.m b/swe_regions.m index 071b5e1..2ea6114 100644 --- a/swe_regions.m +++ b/swe_regions.m @@ -35,7 +35,7 @@ % terms of the first eigenvariate of the filtered and adjusted response in % all suprathreshold voxels within a specified VOI centred on the current % MIP cursor location. Responses are adjusted by removing variance that -% can be predicted by the null space of the F contrast specified (usually +% can be predicted by the null space of the F contrast specified (usually % an F-contrast testing for all effects of interest). % % If temporal filtering has been specified, then the data will be filtered. @@ -44,7 +44,7 @@ % % For a VOI of radius 0, the [adjusted] voxel time-series is returned, and % scaled to have a 2-norm of 1. The actual [adjusted] voxel time series can -% be extracted from xY.y, and will be the same as the [adjusted] data +% be extracted from xY.y, and will be the same as the [adjusted] data % returned by the plotting routine (spm_graph.m) for the same contrast. %__________________________________________________________________________ % Copyright (C) 1999-2016 Wellcome Trust Centre for Neuroimaging @@ -80,7 +80,7 @@ header = get(Finter,'Name'); set(Finter,'Name','VOI time-series extraction'); if ~noGraph, spm_figure('GetWin','Graphics'); end - + %-Find nearest voxel [Euclidean distance] in point list %-------------------------------------------------------------------------- if isempty(xSwE.XYZmm) @@ -95,7 +95,7 @@ spm_XYZreg('GetCoords',hReg),xSwE.XYZmm); xY.xyz = xyz; end - + % and update GUI location %-------------------------------------------------------------------------- if spm_mesh_detect(SwE.xY.VY) @@ -103,8 +103,8 @@ else spm_XYZreg('SetCoords',xyz,hReg); end - - + + %-Get adjustment options and VOI name %-------------------------------------------------------------------------- if ~noGraph @@ -115,11 +115,11 @@ end spm_input(posstr,1,'d','VOI time-series extraction'); end - + if ~isfield(xY,'name') xY.name = spm_input('name of region','!+1','s','VOI'); end - + if ~isfield(xY,'Ic') q(1) = 0; Con = {''}; @@ -137,7 +137,7 @@ i = spm_input('adjust data for (select contrast)','!+1','m',Con); xY.Ic = q(i); end - + %-If fMRI data then ask user to select session %-------------------------------------------------------------------------- if isfield(SwE,'Sess') && ~isfield(xY,'Sess') @@ -154,7 +154,7 @@ [xY, xY.XYZmm, Q] = spm_ROI(xY, xSwE.XYZmm); try, xY = rmfield(xY,'M'); end try, xY = rmfield(xY,'rej'); end - + if isempty(xY.XYZmm) warning('Empty region.'); [Y, xY.y, xY.u, xY.v, xY.s, xY.X0] = deal([]); @@ -175,28 +175,28 @@ return; end end - + %-Extract required data from results files %========================================================================== spm('Pointer','Watch') - -%-Get raw data, whiten and filter + +%-Get raw data, whiten and filter %-------------------------------------------------------------------------- y = swe_data_read(SwE.xY.VY,'xyz',xSwE.XYZ(:,Q)); %y = spm_filter(SwE.xX.K,SwE.xX.W*y); - - + + %-Computation %========================================================================== - + %-Remove null space of contrast %-------------------------------------------------------------------------- if xY.Ic ~= 0 - + %-Parameter estimates: beta = xX.pKX*xX.K*y %---------------------------------------------------------------------- beta = swe_data_read(SwE.Vbeta,'xyz',xSwE.XYZ(:,Q)); - + %-subtract Y0 = XO*beta, Y = Yc + Y0 + e %---------------------------------------------------------------------- if ~isnan(xY.Ic) @@ -204,13 +204,13 @@ else y = y - SwE.xX.xKXs.X * beta; end - + end - + %-Confounds %-------------------------------------------------------------------------- xY.X0 = SwE.xX.xKXs.X(:,[SwE.xX.iB SwE.xX.iG]); - + %-Extract session-specific rows from data and confounds %-------------------------------------------------------------------------- try @@ -218,7 +218,7 @@ y = y(i,:); xY.X0 = xY.X0(i,:); end - + % and add session-specific filter confounds %-------------------------------------------------------------------------- try @@ -227,12 +227,12 @@ try xY.X0 = [xY.X0 SwE.xX.K(xY.Sess).KH]; % Compatibility check end - + %-Remove null space of X0 %-------------------------------------------------------------------------- xY.X0 = xY.X0(:,any(xY.X0)); - - + + %-Compute regional response in terms of first eigenvariate %-------------------------------------------------------------------------- [m,n] = size(y); @@ -251,14 +251,14 @@ u = u*d; v = v*d; Y = u*sqrt(s(1)/n); - + %-Set in structure %-------------------------------------------------------------------------- xY.y = y; xY.u = Y; xY.v = v; xY.s = s; - + %-Display VOI weighting and eigenvariate %========================================================================== if ~noGraph @@ -268,7 +268,7 @@ display_VOI(xY); end end - + %-Save %========================================================================== str = ['VOI_' xY.name '.mat']; @@ -279,7 +279,7 @@ cmd = 'swe_regions(''display'',''%s'')'; fprintf(' VOI saved as %s\n',spm_file(fullfile(SwE.swd,str),'link',cmd)); - + %-Reset title %-------------------------------------------------------------------------- set(Finter,'Name',header); diff --git a/swe_results.m b/swe_results.m index f43bcbd..d0e973d 100644 --- a/swe_results.m +++ b/swe_results.m @@ -1,7 +1,7 @@ function swe_results % Launches batch window containing SwE batch module. % ========================================================================= -% This function prepares (if needed) and launches the batch system with a +% This function prepares (if needed) and launches the batch system with a % job containing the batch module for the specification of % data and design. % ========================================================================= @@ -10,17 +10,17 @@ % Written by Tom Maullin (05/09/2018) % Version Info: $Format:%ci$ $Format:%h$ - % Initiate a job + % Initiate a job if isempty(spm_figure('FindWin','Graphics')) % SPM not running swe_jobman('initcfg') end - + % Launch the batch system with the SwE tab swe_batch % Add the specification module to it swe_jobman('interactive','','swe.results'); return - -end \ No newline at end of file + +end diff --git a/swe_results_ui.m b/swe_results_ui.m index 62bb4ab..46fb021 100644 --- a/swe_results_ui.m +++ b/swe_results_ui.m @@ -6,7 +6,7 @@ % Inputs/Outputs: % - hReg: handle of MIP XYZ registry object % (see spm_XYZreg.m for details) -% - xSwE: structure containing specific SwE data, distribution & filtering +% - xSwE: structure containing specific SwE data, distribution & filtering % details (see spm_getSPM.m for contents) % - SwE: SwE structure containing generic parameters % (see spm_spm.m for contents) @@ -20,11 +20,11 @@ % characterisation of the results of a statistical analysis. % % The user is prompted to select a SwE{T} or SwE{F} image, that is -% thresholded at user specified levels. For a parametric analysis, the -% specification of the contrasts to use and the height and size thresholds -% are described in spm_getSPM.m. For a non-parametric analysis the height +% thresholded at user specified levels. For a parametric analysis, the +% specification of the contrasts to use and the height and size thresholds +% are described in spm_getSPM.m. For a non-parametric analysis the height % threshold may have been set earlier in the matlab batch specification. -% The resulting SwE data is then displayed in the graphics window as a +% The resulting SwE data is then displayed in the graphics window as a % maximum intensity projection, alongside the design matrix and contrasts % employed. % @@ -108,12 +108,12 @@ % body of the code. %__________________________________________________________________________ % Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging - + % Original SPM code by Karl Friston & Andrew Holmes % Adapted by Bryan Guillaume % Based on: swe_results_ui.m 3928 2010-06-16 12:09:22Z guillaume - - + + %========================================================================== % - FORMAT specifications for embedded CallBack functions %========================================================================== @@ -127,7 +127,7 @@ %__________________________________________________________________________ % % FORMAT [hreg,xSwE,SwE] = swe_results_ui('Setup') -% Query SwE and setup GUI. +% Query SwE and setup GUI. % % FORMAT [hreg,xSwE,SwE] = swe_results_ui('Setup',xSwE) % Query SwE and setup GUI using a xSwE input structure. This allows to run @@ -222,7 +222,7 @@ % deletes HandleGraphics objects, but only if they're valid, thus avoiding % warning statements from MATLAB. %__________________________________________________________________________ -% Modified version of spm_results_ui by Bryan Guillaume +% Modified version of spm_results_ui by Bryan Guillaume % Version Info: $Format:%ci$ $Format:%h$ %-Condition arguments @@ -233,7 +233,7 @@ %========================================================================== switch lower(Action), case 'setup' %-Set up results %========================================================================== - + %-Initialise %---------------------------------------------------------------------- [Finter,Fgraph,CmdLine] = spm('FnUIsetup','Stats: Results'); @@ -247,7 +247,7 @@ else [SwE,xSwE] = swe_getSPM; end - + % check format of data file_ext = swe_get_file_extension(SwE.xY.P{1}); isMat = strcmpi(file_ext,'.mat'); @@ -256,12 +256,12 @@ while(isMat) [SwE,xSwE] = swe_getSPM(xSwE); end - - if isempty(xSwE) + + if isempty(xSwE) varargout = {[],[],[]}; return; end - + switch xSwE.STAT case 'T' eSTAT = 'Z'; @@ -279,7 +279,7 @@ %-Ensure pwd = swd so that relative filenames are valid %---------------------------------------------------------------------- cd(SwE.swd) - + %-Get space information %====================================================================== M = SwE.xVol.M; @@ -315,7 +315,7 @@ catch datatype = 'Volumetric (2D/3D)'; end - + switch datatype case 'Volumetric (2D/3D)' units = {'mm' 'mm' 'mm'}; @@ -347,22 +347,22 @@ end xSwE.units = units; SwE.xVol.units = units; - - + + %-Setup Results User Interface; Display MIP, design matrix & parameters %====================================================================== - + %-Setup results GUI %---------------------------------------------------------------------- spm_clf(Finter); spm('FigName',['SwE{',xSwE.STAT,'}: Results'],Finter,CmdLine); hReg = swe_results_ui('SetupGUI',M,DIM,xSwE,Finter); - + %-Setup design interrogation menu %---------------------------------------------------------------------- hDesRepUI = spm_DesRep('DesRepUI',SwE); figure(Finter) - + %-Setup contrast menu %---------------------------------------------------------------------- hC = uimenu(Finter,'Label','Contrasts', 'Tag','ContrastsUI'); @@ -398,10 +398,10 @@ xSwEtmp = xSwE; xSwEtmp.thresDesc = 'p<0.001 (unc.)'; xSwEtmp.infType = 0; uimenu(hC1,'Label','Set to 0.001 (unc.)','UserData',struct('Ic',xSwE.Ic),... 'Callback',{@mychgcon,xSwEtmp}); - + uimenu(hC1,'Label',[xSwE.thresDesc ', k=' num2str(xSwE.k)],... 'Enable','off','Separator','on'); - + %-Setup Maximum intensity projection (MIP) & register %---------------------------------------------------------------------- % create a graphic window for each brain structure @@ -421,7 +421,7 @@ hRotate3d.Enable = 'on'; end for i = 1:numel(SwE.cifti.surfaces) - + if mod(it, gridSize) == 1 offset_y = offset_y + 0.36/gridSize; offset_x = 0; @@ -456,14 +456,14 @@ end end if numel(SwE.cifti.volume) > 0 - + if mod(it, gridSize) == 1 offset_y = offset_y + 0.36/gridSize; offset_x = 0; else offset_x = offset_x + 0.55/gridSize; end - + hMIPax{it} = axes('Parent', Fgraph, 'Position', [0.05 + offset_x, 0.60 + offset_y, 0.55/gridSize, 0.36/gridSize], 'Visible','off'); [isSurviving, indInCifti] = ismember(SwE.cifti.volume.indices, xSwE.XYZ(1,:)); indInCifti = indInCifti(isSurviving); @@ -476,13 +476,13 @@ 'Interpreter','none',... 'FontSize',FS(sizeFont),'Fontweight','Bold',... 'Parent',hMIPax{it}); - + if exist('hRotate3d', 'var') setAllowAxesRotate(hRotate3d, hMIPax{it}, false); end end else - hMIPax = axes('Parent',Fgraph,'Position',[0.05 0.60 0.55 0.36],'Visible','off'); + hMIPax = axes('Parent',Fgraph,'Position',[0.05 0.60 0.55 0.36],'Visible','off'); nBrainStructure = 1; if xSwE.STAT == 'P' str = xSwE.STATstr; @@ -543,8 +543,8 @@ 'HorizontalAlignment','center',... 'VerticalAlignment','baseline',... 'FontWeight','Bold','FontSize',FS(14)) - - + + %-Print SwEresults: Results directory & thresholding info %---------------------------------------------------------------------- hResAx = axes('Parent',Fgraph,... @@ -574,7 +574,7 @@ elseif isfield(xSwE, 'WB') && xSwE.WB && xSwE.infType == 1 && strcmp(xSwE.clustWise, 'Uncorr') text(0,12,sprintf('Wild Bootstrap p-value threshold %s',thresDesc),'Parent',hResAx) else - text(0,12,sprintf('Height threshold %c > %0.6f {%s}',eSTAT,xSwE.u,thresDesc),'Parent',hResAx) + text(0,12,sprintf('Height threshold %c > %0.6f {%s}',eSTAT,xSwE.u,thresDesc),'Parent',hResAx) end else text(0,12,sprintf('Height threshold %c > %0.6f',eSTAT,xSwE.u),'Parent',hResAx) @@ -582,7 +582,7 @@ catch text(0,12,sprintf('Height threshold %c > %0.6f',eSTAT,xSwE.u),'Parent',hResAx) end - if strcmp(xSwE.clustWise, 'FWE') + if strcmp(xSwE.clustWise, 'FWE') if ~isfield(xSwE, 'clusterSizeType') || strcmp(xSwE.clusterSizeType, 'classic k_E') text(0,00,sprintf('Wild Bootstrap extent threshold k > %0.0f %s {p<=%0.3f (FWE)}', xSwE.k, strDataType, xSwE.fwep_c), 'Parent', hResAx) elseif strcmp(xSwE.clusterSizeType, 'Box-Cox norm. k_{Z}') @@ -610,7 +610,7 @@ error('Unknown TFCE threshold.') end end - + %-Plot design matrix %---------------------------------------------------------------------- hDesMtx = axes('Parent',Fgraph,'Position',[0.65 0.55 0.25 0.25]); @@ -621,7 +621,7 @@ 'X', SwE.xX.X,... 'fnames', {reshape({SwE.xY.VY.fname},size(SwE.xY.VY))},... 'Xnames', {SwE.xX.name})) - + if exist('hRotate3d', 'var') setAllowAxesRotate(hRotate3d, hDesMtx, false); end @@ -637,12 +637,12 @@ 'Tag','ConGrphAx','Visible','off'); if exist('hRotate3d', 'var') setAllowAxesRotate(hRotate3d, hConAx, false); - end + end title('contrast(s)') htxt = get(hConAx,'title'); set(htxt,'Visible','on','HandleVisibility','on') end - + for ii = nCon:-1:1 hTmp = axes('Position',[0.65 (0.80 + dy*(nCon - ii +.1)) 0.25 dy*.9]); if xCon(xSwE.Ic(ii)).STAT == 'T' && size(xCon(xSwE.Ic(ii)).c,2) == 1 @@ -657,9 +657,9 @@ 'YTick',[-1,0,+1],'YTickLabel','',... 'YLim',[min(xCon(xSwE.Ic(ii)).c),max(xCon(xSwE.Ic(ii)).c)] +... [-1 +1] * max(abs(xCon(xSwE.Ic(ii)).c))/10 ) - + else - + %-F-contrast - image %-------------------------------------------------------------- h = image((xCon(xSwE.Ic(ii)).c'/max(abs(xCon(xSwE.Ic(ii)).c(:)))+1)*32); @@ -670,7 +670,7 @@ 'YTick',[0:size(SwE.xCon(xSwE.Ic(ii)).c,2)]+0.5,.... 'YTickLabel','',... 'YLim', [0,size(xCon(xSwE.Ic(ii)).c,2)]+0.5 ) - + end ylabel(num2str(xSwE.Ic(ii))) set(h,'ButtonDownFcn','spm_DesRep(''SurfCon_CB'')',... @@ -680,10 +680,10 @@ if exist('hRotate3d', 'var') setAllowAxesRotate(hRotate3d, hTmp, false); - end + end end - - + + %-Store handles of results section Graphics window objects %---------------------------------------------------------------------- H = get(Fgraph,'Children'); @@ -691,11 +691,11 @@ H = findobj(H); Hv = get(H,'Visible'); set(hResAx,'Tag','PermRes','UserData',struct('H',H,'Hv',{Hv})) - + %-Load whole brain table. %------------------------------------------------------------------ swe_list('List',xSwE,hReg); - + %-Finished results setup %---------------------------------------------------------------------- varargout = {hReg,xSwE,SwE}; @@ -703,8 +703,8 @@ assignin ('base','xSwE',xSwE) assignin ('base','SwE',SwE) spm('Pointer','Arrow') - - + + %====================================================================== case 'setupgui' %-Set up results section GUI %====================================================================== @@ -717,7 +717,7 @@ Finter = spm_figure('GetWin',Finter); WS = spm('WinScale'); FS = spm('FontSizes'); - + %-Create frame for Results GUI objects %------------------------------------------------------------------ hPan = uipanel('Parent',Finter,'Title','','Units','Pixels',... @@ -728,11 +728,11 @@ 'BorderType','Etchedin', ... 'Position',[005 005 390 180].*WS,... 'BackgroundColor',[179 179 179]/255); - + %-Initialise registry in hReg frame object %------------------------------------------------------------------ [hReg,xyz] = spm_XYZreg('InitReg',hReg,M,DIM,[0;0;0]); - + %-Setup editable XYZ widgets & cross register with registry %------------------------------------------------------------------ hFxyz = swe_results_ui('DrawXYZgui',M,DIM,varargin{4},xyz,hReg); @@ -741,16 +741,16 @@ %-Set up buttons for results functions %------------------------------------------------------------------ swe_results_ui('DrawButts',hReg,DIM,Finter,WS,FS); - + if spm_check_version('matlab','7.11') ~= 0 drawnow; % required to force "ratio locking" set(findobj(hPan),'Units','Normalized','FontUnits','Normalized'); end - + varargout = {hReg}; - - - + + + %====================================================================== case 'drawbutts' %-Draw results section buttons in Interactive window %====================================================================== @@ -763,7 +763,7 @@ else Finter = varargin{4}; end if nargin < 5, WS = spm('WinScale'); else WS = varargin{5}; end if nargin < 6, FS = spm('FontSizes'); else FS = varargin{6}; end - + %-p-values %------------------------------------------------------------------ hPan = uipanel('Parent',hReg,'Title','p-values','Units','Pixels',... @@ -798,7 +798,7 @@ 'Callback','TabDat = swe_VOI(SwE,xSwE,hReg);',... 'Interruptible','on','Enable','on',... 'Position',[005 005 100 020].*WS); - + %-SwE area - used for Volume of Interest analyses %------------------------------------------------------------------ hPan = uipanel('Parent',hReg,'Title','Multivariate','Units','Pixels',... @@ -827,7 +827,7 @@ 'FontSize',FS(10),... 'ForegroundColor',[1 1 1],... 'BackgroundColor',[179 179 179]/255); - + uicontrol('Parent',hPan,'Style','PushButton','String','plot',... 'FontSize',FS(10),... 'ToolTipString','plot data & contrasts at current voxel',... @@ -835,11 +835,11 @@ 'Interruptible','on','Enable','on',... 'Position',[005 055 100 020].*WS,... 'Tag','plotButton'); - + str = { 'overlays...','slices','sections','montage','render','previous sections','previous render'}; tstr = { 'overlay filtered SPM on another image: ',... '3 slices / ','slice overlay /','ortho sections / ','render /','previous ortho sections /','previous surface rendering'}; - + tmp = { 'spm_transverse(''set'',xSwE,hReg)',... 'spm_sections(xSwE,hReg)',... {@myslover},... @@ -849,12 +849,12 @@ '''dim'', xSwE.DIM))'],... ['global prevsect;','spm_sections(xSwE,hReg,prevsect)'],... ['global prevrend;','if ~isstruct(prevrend)',... - 'prevrend = struct(''rendfile'','''',''brt'',[],''col'',[]); end;',... + 'prevrend = struct(''rendfile'','''',''brt'',[],''col'',[]); end;',... 'spm_render( struct( ''XYZ'', xSwE.XYZ,',... '''t'', xSwE.Z'',',... '''mat'', xSwE.M,',... '''dim'', xSwE.DIM),prevrend.brt,prevrend.rendfile)']}; - + uicontrol('Parent',hPan,'Style','popupmenu','String',str,... 'FontSize',FS(10),... 'ToolTipString',cat(2,tstr{:}),... @@ -862,7 +862,7 @@ 'UserData',tmp,... 'Interruptible','on','Enable','on',... 'Position',[005 030 100 020].*WS); - + str = {'save...',... 'thresholded SPM',... 'all clusters (binary)',... @@ -872,7 +872,7 @@ {@mysavespm, 'binary' },... {@mysavespm, 'n-ary' },... {@mysavespm, 'current'}}; - + uicontrol('Parent',hPan,'Style','popupmenu','String',str,... 'FontSize',FS(10),... 'ToolTipString','Save as image',... @@ -880,7 +880,7 @@ 'UserData',tmp,... 'Interruptible','on','Enable','on',... 'Position',[005 005 100 020].*WS); - + %-ResultsUI controls %------------------------------------------------------------------ uicontrol('Parent',hReg,'Style','PushButton','String','clear',... @@ -892,15 +892,15 @@ 'Interruptible','on','Enable','on',... 'DeleteFcn','spm_clf(''Graphics'')',... 'Position',[280 050 048 020].*WS); - + uicontrol('Parent',hReg,'Style','PushButton','String','exit',... 'ToolTipString','Exit the results section',... 'FontSize',FS(9),'ForegroundColor','r',... 'Callback','swe_results_ui(''close'')',... 'Interruptible','on','Enable','on',... 'Position',[332 050 048 020].*WS); - - + + %====================================================================== case 'drawxyzgui' %-Draw XYZ GUI area %====================================================================== @@ -912,13 +912,13 @@ DIM = varargin{3}; M = varargin{2}; xyz = spm_XYZreg('RoundCoords',xyz,M,DIM); - + %-Font details %------------------------------------------------------------------ WS = spm('WinScale'); FS = spm('FontSizes'); PF = spm_platform('fonts'); - + %-Create XYZ control objects %------------------------------------------------------------------ hFxyz = uicontrol(Finter,'Style','Pushbutton',... @@ -929,7 +929,7 @@ 'FontSize',FS(10),... 'HorizontalAlignment','Left',... 'ForegroundColor','w') - + uicontrol(Finter,'Style','Text','String','x =',... 'Position',[020 015 024 018].*WS,... 'FontName',PF.times,'FontSize',FS(10),'FontAngle','Italic',... @@ -941,7 +941,7 @@ 'HorizontalAlignment','Right',... 'Tag','hX',... 'Callback','swe_results_ui(''EdWidCB'')'); - + uicontrol(Finter,'Style','Text','String','y =',... 'Position',[105 015 024 018].*WS,... 'FontName',PF.times,'FontSize',FS(10),'FontAngle','Italic',... @@ -953,7 +953,7 @@ 'HorizontalAlignment','Right',... 'Tag','hY',... 'Callback','swe_results_ui(''EdWidCB'')'); - + if DIM(3) ~= 1 uicontrol(Finter,'Style','Text','String','z =',... 'Position',[190 015 024 018].*WS,... @@ -969,7 +969,7 @@ else hZ = []; end - + %-Statistic value reporting pane %------------------------------------------------------------------ uicontrol(Finter,'Style','Text','String','statistic',... @@ -982,8 +982,8 @@ 'Position',[285 012 100 020].*WS,... 'FontSize',FS(10),... 'HorizontalAlignment','Center'); - - + + %-Store data %------------------------------------------------------------------ set(hFxyz,'Tag','hFxyz','UserData',struct(... @@ -997,23 +997,23 @@ 'hZ', hZ,... 'hSPM', hSPM,... 'xyz', xyz )); - + set([hX,hY,hZ],'UserData',hFxyz) varargout = {hFxyz}; - - + + %====================================================================== case 'edwidcb' %-Callback for editable widgets %====================================================================== % swe_results_ui('EdWidCB') - + hC = gcbo; d = find(strcmp(get(hC,'Tag'),{'hX','hY','hZ'})); hFxyz = get(hC,'UserData'); UD = get(hFxyz,'UserData'); xyz = UD.xyz; nxyz = xyz; - + o = evalin('base',['[',get(hC,'String'),']'],'sprintf(''error'')'); if ischar(o) || length(o)>1 warning(sprintf('%s: Error evaluating ordinate:\n\t%s',... @@ -1022,15 +1022,15 @@ nxyz(d) = o; nxyz = spm_XYZreg('RoundCoords',nxyz,UD.M,UD.DIM); end - + if abs(xyz(d)-nxyz(d))>0 UD.xyz = nxyz; set(hFxyz,'UserData',UD) if ~isempty(UD.hReg), spm_XYZreg('SetCoords',nxyz,UD.hReg,hFxyz); end set(hC,'String',sprintf('%.3f',nxyz(d))) swe_results_ui('UpdateSPMval',UD) end - - + + %====================================================================== case 'updatespmval' %-Update SwE value in GUI %====================================================================== @@ -1041,8 +1041,8 @@ i = spm_XYZreg('FindXYZ',UD.xyz,UD.XYZ); if isempty(i), str = ''; else str = sprintf('%6.2f',UD.Z(i)); end set(UD.hSPM,'String',str); - - + + %====================================================================== case 'getcoords' % Get current co-ordinates from XYZ widget %====================================================================== @@ -1050,8 +1050,8 @@ if nargin<2, hFxyz='Interactive'; else hFxyz=varargin{2}; end hFxyz = swe_results_ui('FindXYZframe',hFxyz); varargout = {getfield(get(hFxyz,'UserData'),'xyz')}; - - + + %====================================================================== case 'setcoords' % Set co-ordinates to XYZ widget %====================================================================== @@ -1059,12 +1059,12 @@ if nargin<4, hC=0; else hC=varargin{4}; end if nargin<3, hFxyz=swe_results_ui('FindXYZframe'); else hFxyz=varargin{3}; end if nargin<2, error('Set co-ords to what!'); else xyz=varargin{2}; end - + %-If this is an internal call, then don't do anything if hFxyz==hC, return, end - + UD = get(hFxyz,'UserData'); - + %-Check validity of coords only when called without a caller handle %------------------------------------------------------------------ if hC <= 0 @@ -1076,7 +1076,7 @@ else d = []; end - + %-Update xyz information & widget strings %------------------------------------------------------------------ UD.xyz = xyz; set(hFxyz,'UserData',UD) @@ -1084,17 +1084,17 @@ set(UD.hY,'String',sprintf('%.2f',xyz(2))) set(UD.hZ,'String',sprintf('%.2f',xyz(3))) swe_results_ui('UpdateSPMval',UD) - + %-Tell the registry, if we've not been called by the registry... %------------------------------------------------------------------ if (~isempty(UD.hReg) && UD.hReg~=hC) spm_XYZreg('SetCoords',xyz,UD.hReg,hFxyz); end - + %-Return arguments %------------------------------------------------------------------ varargout = {xyz,d}; - + %====================================================================== case 'findxyzframe' % Find hFxyz frame @@ -1115,12 +1115,12 @@ %====================================================================== % swe_results_ui('PlotUi',hAx) if nargin<2, hAx=gca; else hAx=varargin{2}; end - + WS = spm('WinScale'); FS = spm('FontSizes'); Finter=spm_figure('FindWin','Interactive'); figure(Finter) - + %-Check there aren't already controls! %------------------------------------------------------------------ hGraphUI = findobj(Finter,'Tag','hGraphUI'); @@ -1132,7 +1132,7 @@ delete(findobj(Finter,'Tag','hGraphUIbg')) end end - + %-Frames & text %------------------------------------------------------------------ hGraphUIbg = uicontrol(Finter,'Style','Frame','Tag','hGraphUIbg',... @@ -1148,7 +1148,7 @@ 'FontAngle','Italic','FontSize',FS(10),... 'HorizontalAlignment','Left',... 'ForegroundColor','w'); - + %-Controls %------------------------------------------------------------------ h1 = uicontrol(Finter,'Style','CheckBox','String','hold',... @@ -1165,7 +1165,7 @@ 'Tag','holdButton',... 'Position',[015 210 070 020].*WS); set(findobj('Tag','plotButton'),'UserData',h1); - + h2 = uicontrol(Finter,'Style','CheckBox','String','grid',... 'ToolTipString','toggle axes grid',... 'FontSize',FS(10),... @@ -1206,12 +1206,12 @@ 'Callback','swe_results_ui(''PlotUiCB'')',... 'Interruptible','off','Enable','on',... 'Position',[315 210 070 020].*WS); - + %-Handle storage for linking, and DeleteFcns for linked deletion %------------------------------------------------------------------ set(hGraphUI,'UserData',[h1,h2,h3,h4,h5]) set([h1,h2,h3,h4,h5],'UserData',hAx) - + set(hGraphUIbg,'UserData',... [hGraphUI,hGraphUIButtsF,hText,h1,h2,h3,h4,h5],... 'DeleteFcn','swe_results_ui(''Delete'',get(gcbo,''UserData''))') @@ -1228,7 +1228,7 @@ if v==1, return, end str = cellstr(get(hPM,'String')); str = str{v}; - + hAx = get(hPM,'UserData'); switch str case 'Title' @@ -1254,31 +1254,31 @@ otherwise warning(['Unknown action: ',str]) end - + set(hPM,'Value',1) - - + + %====================================================================== case 'clear' %-Clear results subpane %====================================================================== % Fgraph = swe_results_ui('Clear',F,mode) % mode 1 [default] usual, mode 0 - clear & hide Res stuff, 2 - RNP - + if nargin<3, mode=1; else, mode=varargin{3}; end if nargin<2, F='Graphics'; else, F=varargin{2}; end F = spm_figure('FindWin',F); - + %-Clear input objects from 'Interactive' window %------------------------------------------------------------------ %spm_input('!DeleteInputObj') - - + + %-Get handles of objects in Graphics window & note permanent results objects %------------------------------------------------------------------ H = get(F,'Children'); %-Get contents of window H = findobj(H,'flat','HandleVisibility','on'); %-Drop GUI components h = findobj(H,'flat','Tag','PermRes'); %-Look for 'PermRes' object - + if ~isempty(h) %-Found 'PermRes' object % This has handles of permanent results objects in it's UserData @@ -1291,23 +1291,23 @@ HRv = {}; end H = setdiff(H,HR); %-Drop permanent results obj - - + + %-Delete stuff as appropriate %------------------------------------------------------------------ if mode==2 %-Don't delete axes with NextPlot 'add' H = setdiff(H,findobj(H,'flat','Type','axes','NextPlot','add')); end - + delete(H) - + if mode==0 %-Hide the permanent results section stuff set(HR,'Visible','off') else set(HR,{'Visible'},HRv) end - - + + %====================================================================== case 'launchmp' %-Launch multiplanar toolbox %====================================================================== @@ -1316,7 +1316,7 @@ hReg = varargin{4}; DIM = varargin{3}; M = varargin{2}; - + %-Check for existing MultiPlanar toolbox hMP = get(hBmp,'UserData'); if ishandle(hMP) @@ -1324,31 +1324,31 @@ varargout = {hMP}; return end - + %-Initialise and cross-register MultiPlanar toolbox hMP = spm_XYZreg_Ex2('Create',M,DIM); spm_XYZreg('Xreg',hReg,hMP,'spm_XYZreg_Ex2'); - + %-Setup automatic deletion of MultiPlanar on deletion of results controls set(hBmp,'Enable','on','UserData',hMP) set(hBmp,'DeleteFcn','swe_results_ui(''delete'',get(gcbo,''UserData''))') - + varargout = {hMP}; - - + + %====================================================================== case 'delete' %-Delete HandleGraphics objects %====================================================================== % swe_results_ui('Delete',h) h = varargin{2}; delete(h(ishandle(h))); - - + + %====================================================================== otherwise %====================================================================== error('Unknown action string') - + end %========================================================================== @@ -1399,10 +1399,10 @@ function mysavespm(action) switch lower(action) case 'thresh' Z = xSwE.Z; - + case 'binary' Z = ones(size(xSwE.Z)); - + case 'n-ary' if ~isfield(xSwE,'G') Z = spm_clusters(XYZ); @@ -1417,12 +1417,12 @@ function mysavespm(action) C = spm_mesh_clusters(xSwE.G,C); Z = C(xSwE.XYZ(1,:)); end - + case 'current' [xyzmm,i] = spm_XYZreg('NearestXYZ',... spm_results_ui('GetCoords'),xSwE.XYZmm); spm_results_ui('SetCoords',xSwE.XYZmm(:,i)); - + if ~isfield(xSwE,'G') A = spm_clusters(XYZ); j = find(A == A(i)); @@ -1435,7 +1435,7 @@ function mysavespm(action) C = C==C(xSwE.XYZ(1,i)); Z = C(xSwE.XYZ(1,:)); end - + otherwise error('Unknown action.'); end @@ -1466,7 +1466,7 @@ function mysavespm(action) %========================================================================== spm_input('!DeleteInputObj'); xSwE = evalin('base','xSwE;'); - + so = slover; [img,sts] = spm_select(1,'image','Select image for rendering on'); if ~sts, return; end @@ -1476,23 +1476,23 @@ function mysavespm(action) %[mx,mn] = slover('volmaxmin', obj.img.vol); %obj.img.range = [mn mx]; so.img.prop = 1; - + so = add_spm(so,xSwE); - + so.transform = deblank(spm_input('Image orientation', '+1', ... 'Axial|Coronal|Sagittal', char('axial','coronal','sagittal'), 1)); so = fill_defaults(so); slices = so.slices; so.slices = spm_input('Slices to display (mm)', '+1', 'e', ... sprintf('%0.0f:%0.0f:%0.0f',slices(1),mean(diff(slices)),slices(end))); - + so.figure = spm_figure('GetWin', 'SliceOverlay'); so = paint(so); assignin('base','so',so); - + function flag = isNotAnSPMMeshRender(obj, event_obj) if strcmpi(obj.Tag, 'SPMMeshRender') flag = false; else flag = true; - end \ No newline at end of file + end diff --git a/swe_rmodel.m b/swe_rmodel.m index d9aa35e..41f167a 100644 --- a/swe_rmodel.m +++ b/swe_rmodel.m @@ -1,7 +1,7 @@ function swe_rmodel % Launches batch window containing SwE batch module for running a design. % ========================================================================= -% This function prepares (if needed) and launches the batch system with a +% This function prepares (if needed) and launches the batch system with a % job containing the batch module for the computation of a prespecified % design. % ========================================================================= @@ -10,17 +10,17 @@ % Written by Tom Maullin (05/09/2018) % Version Info: $Format:%ci$ $Format:%h$ - % Initiate a job + % Initiate a job if isempty(spm_figure('FindWin','Graphics')) % SPM not running swe_jobman('initcfg') end - + % Launch the batch system with the SwE tab swe_batch % Add the specification module to it swe_jobman('interactive','','swe.rmodel'); return - -end \ No newline at end of file + +end diff --git a/swe_run_results.m b/swe_run_results.m index 9ed9173..89e59a0 100644 --- a/swe_run_results.m +++ b/swe_run_results.m @@ -19,6 +19,6 @@ function swe_run_results(varargin) % Display results. % ------------------------------------------------------------------------- -swe_results_ui; +swe_results_ui; -end \ No newline at end of file +end diff --git a/swe_run_rmodel.m b/swe_run_rmodel.m index 5add3c2..9d01f97 100644 --- a/swe_run_rmodel.m +++ b/swe_run_rmodel.m @@ -33,4 +33,4 @@ function swe_run_rmodel(varargin) cd(SwE.swd); swe_cp(SwE); -end \ No newline at end of file +end diff --git a/swe_run_smodel.m b/swe_run_smodel.m index 0492892..ffe01f3 100644 --- a/swe_run_smodel.m +++ b/swe_run_smodel.m @@ -35,7 +35,7 @@ end % If we've gotten to this point we're committed to overwriting files. -% Delete them so we don't get stuck +% Delete them so we don't get stuck %-------------------------------------------------------------------------- files = {'^mask\..{3}$','^ResMS\..{3}$','^RPV\..{3}$',... '^beta_.{4}\..{3}$','^con_.{4}\..{3}$','^ResI_.{4}\..{3}$',... @@ -157,7 +157,7 @@ switch char(fieldnames(job.type)) case 'modified' - + Vis.iVis = job.type.modified.visits; Vis.nVis = length(unique(Vis.iVis)); Gr.iGr = job.type.modified.groups; @@ -198,7 +198,7 @@ % If it's a manually created structure with field 'R' mandatory % containing matrix and field 'names' optionally containing a % cell array of names. - if isfield(tmp,'R') + if isfield(tmp,'R') R = tmp.R; if isfield(tmp,'names') names = tmp.names; @@ -304,7 +304,7 @@ % dimensions and orientation / voxel size %========================================================================== -fprintf('%-40s: ','Mapping files') +fprintf('%-40s: ','Mapping files') P = job.scans; file_ext = swe_get_file_extension(P{1}); isMat = strcmpi(file_ext,'.mat'); @@ -346,13 +346,13 @@ for ii = 1:nSurfaceBrainStructures if strcmpi(job.ciftiAdditionalInfo.ciftiGeomFile(ii).brainStructureLabel, SwE.cifti.surfaces{i}.brainStructure) SwE.cifti.surfaces{i}.geomFile = char(job.ciftiAdditionalInfo.ciftiGeomFile(ii).geomFile); - if isfield(job.ciftiAdditionalInfo.ciftiGeomFile(ii), 'areaFile') && ~isempty(job.ciftiAdditionalInfo.ciftiGeomFile(ii).areaFile) && ~strcmpi(job.ciftiAdditionalInfo.ciftiGeomFile(ii).areaFile, '') + if isfield(job.ciftiAdditionalInfo.ciftiGeomFile(ii), 'areaFile') && ~isempty(job.ciftiAdditionalInfo.ciftiGeomFile(ii).areaFile) && ~strcmpi(job.ciftiAdditionalInfo.ciftiGeomFile(ii).areaFile, '') SwE.cifti.surfaces{i}.areaFile = char(job.ciftiAdditionalInfo.ciftiGeomFile(ii).areaFile); end break; end if (ii == nSurfaceBrainStructures) - error('At least one of the surface brain structure label in the CIfTI files cannot be found in those specified by the user. Please revise your specification.'); + error('At least one of the surface brain structure label in the CIfTI files cannot be found in those specified by the user. Please revise your specification.'); end end end @@ -369,7 +369,7 @@ end end -fprintf('%30s\n','...done') +fprintf('%30s\n','...done') %-Global values, scaling and global normalisation @@ -623,7 +623,7 @@ 'type', 3,... 'cols',2 + size([H C B G],2),... 'descrip', {str} ); - + G = [G,fB,fW]; Gnames = [Gnames; {gnamesB}; {gnamesW}]; if isempty(xC), xC = [tmpB,tmpW]; else xC = [xC,tmpB,tmpW]; end @@ -738,7 +738,7 @@ if any(abs(ones(N,1)-P_x*ones(N,1))>sqrt(eps)) error(['Input model does not include an intercept. You must '... 'include an intercept in this model.']); - + end %-Design description (an nx2 cellstr) - for saving and display @@ -763,58 +763,58 @@ % - WB configuration - Only if needed %========================================================================== if isfield(job.WB, 'WB_yes') - + WB.SS = job.WB.WB_yes.WB_ss; WB.nB = job.WB.WB_yes.WB_nB; WB.RSwE = job.WB.WB_yes.WB_SwE; WB.voxelWiseInfo = []; switch char(fieldnames(job.WB.WB_yes.WB_infType)) - + case 'WB_voxelwise' WB.clusterWise = 0; WB.voxelWise = 1; - + case 'WB_clusterwise' - + WB.clusterWise = 1; WB.voxelWise = 0; - + % Cluster forming threshold. WB.clusterInfo.primaryThreshold = job.WB.WB_yes.WB_infType.WB_clusterwise.WB_clusThresh; if WB.clusterInfo.primaryThreshold > 1 || WB.clusterInfo.primaryThreshold < 0 error('cluster-forming threshold should be between 0 an 1 (this is a probability)'); end - + % Work out which type of file we are looking at. inputType = job.WB.WB_yes.WB_infType.WB_clusterwise.WB_inputType; - + % If we are looking at '.mat' we need more information. if isfield(inputType, 'WB_mat') - + % Check whether we are looking at surface data. if isfield(inputType.WB_mat, 'WB_surface') WB.clusterInfo.Vfaces = inputType.WB_mat.WB_surface.WB_surfacemask; else WB.clusterInfo.Vxyz = inputType.WB_mat.WB_volumetric.WB_volumetricmask; end - + end - + case 'WB_TFCE' - + % We have no clusterwise results for TFCE - WB.clusterWise = 0; + WB.clusterWise = 0; WB.voxelWise = 0; - + % Create TFCE structure for TFCE analysis. WB.TFCE.H = job.WB.WB_yes.WB_infType.WB_TFCE.WB_TFCE_H; WB.TFCE.E = job.WB.WB_yes.WB_infType.WB_TFCE.WB_TFCE_E; - + % This is by default set to 0.1 as recommended per Smith & % Nichols (2007). If a user wishes to change this value, change % it on the below line: WB.TFCE.dh = 0.1; - + % Error if '.mat' input. if isMat error('TFCE is not currently available for ''.mat'' input.') @@ -828,11 +828,11 @@ if isCifti error('TFCE is not currently available for CIfTI data input.') end - + end - + switch char(fieldnames(job.WB.WB_yes.WB_stat)) - + case 'WB_T' WB.stat = 'T'; WB.con = job.WB.WB_yes.WB_stat.WB_T.WB_T_con; @@ -849,7 +849,7 @@ % Pad with zeros WB.con = [WB.con zeros(1,size(X,2)-size(WB.con,2))]; end - + case 'WB_F' WB.stat = 'F'; WB.con = job.WB.WB_yes.WB_stat.WB_F.WB_F_con; @@ -860,11 +860,11 @@ % Pad with zeros WB.con = [WB.con zeros(size(WB.con,1),size(X,2)-size(WB.con,2))]; end - + otherwise error('unexpected statistic type'); end - + end @@ -906,7 +906,7 @@ %-Display Design report %========================================================================== if ~spm('CmdLine') - fprintf('%-40s: ','Design reporting') + fprintf('%-40s: ','Design reporting') if strcmpi(file_ext,'.mat') fname = cellstr(repmat(' ', nScan, 1)); else diff --git a/swe_seed.m b/swe_seed.m index 3e42d76..ff245e1 100644 --- a/swe_seed.m +++ b/swe_seed.m @@ -12,8 +12,8 @@ function swe_seed(varargin) % ------------------------------------------------------------------------- % Inputs: % - Seed: Value for random number generator seed. -% Deterministically seed the random number generator with value Seed. -% Useful for creating reproducible random numbers for testing; then make +% Deterministically seed the random number generator with value Seed. +% Useful for creating reproducible random numbers for testing; then make % sure that "shuffle_seed" default is false to prevent re-seeding. % ========================================================================= % T. Nichols diff --git a/swe_select.m b/swe_select.m index 8239c8f..bfeefbc 100644 --- a/swe_select.m +++ b/swe_select.m @@ -1,6 +1,6 @@ function varargout = swe_select(varargin) % File selector - % This is a copy of swe_select but able to handle CIfTI files as well. + % This is a copy of swe_select but able to handle CIfTI files as well. % FORMAT [t,sts] = swe_select(n,typ,mesg,sel,wd,filt,frames) % n - number of files [Default: Inf] % A single value or a range. e.g. @@ -43,7 +43,7 @@ % As above, but for selecting frames of 4D NIfTI files % frames - vector of frames to select (defaults to 1, if not specified). % If the frame number is Inf, all frames for the matching images - % are listed. + % are listed. % % FORMAT [files,dirs] = swe_select('FPList',direc,filt) % FORMAT [files,dirs] = swe_select('ExtFPList',direc,filt,frames) @@ -63,8 +63,8 @@ % John Ashburner % Modified by Bryan Guillaume % Version Info: $Format:%ci$ $Format:%h$ - - + + % For developers: %-------------------------------------------------------------------------- % FORMAT cpath = swe_select('CPath',path,cwd) @@ -91,17 +91,17 @@ % % FORMAT files = swe_select('Expand',files) % Return a list of image filenames with appended frame numbers. - + persistent isInitSelect; if isempty(isInitSelect) isInitSelect = true; swe_select('init'); if nargin == 1 && strcmpi(varargin{1},'init'), return; end end - + %-Commands that are not passed to cfg_getfile local_cmds = {'regfilter', 'init', 'expand'}; - + if nargin && ischar(varargin{1}) && any(strcmpi(varargin{1},local_cmds)) switch lower(varargin{1}) case 'init' @@ -110,7 +110,7 @@ end swe_select('regfilter'); swe_select('prevdirs',spm('Dir')); - + case 'regfilter' % Regexp based filters without special handlers cfg_getfile('regfilter', 'mesh', {'\.gii$'}); @@ -126,7 +126,7 @@ cfg_getfile('regfilter', 'image',... {'.*\.nii(,\d+){0,2}$','.*\.img(,\d+){0,2}$'},... false, @swe_select_image, {frames}); - + case 'expand' varargout{1} = swe_select_expand(varargin{2}); if ischar(varargin{2}), varargout{1} = char(varargout{1}); end @@ -151,9 +151,9 @@ elseif nargin > 6 && isnumeric(varargin{1}) varargin{7} = struct('frames', varargin{7}); end - + [t, sts] = cfg_getfile(varargin{:}); - + % cfg_getfile returns cellstr arrays, convert to char arrays if nargin > 0 && ischar(varargin{1}) switch lower(varargin{1}) @@ -168,14 +168,14 @@ else t = char(t); end - + varargout{1} = t; if nargout > 1 varargout{2} = sts; end end - - + + %========================================================================== % FUNCTION varargout = swe_select_image(cmd, varargin) %========================================================================== @@ -211,11 +211,11 @@ elseif numel(frames)==1 && frames(1)==1 [ii{:}] = deal(1); end - + % if domsg % msg(ob,['Listing ' num2str(numel(f)) ' files...']); % end - + % Combine filename and frame number(s) nii = cellfun(@numel, ii); cfiles = cell(sum(nii),1); @@ -259,8 +259,8 @@ ofiles = [ofiles; vfiles]; end end - - + + %========================================================================== % FUNCTION n = swe_select_get_nbframes(file) %========================================================================== @@ -275,4 +275,3 @@ dim = [N.dat.dim 1 1 1 1 1]; n = dim(4); end - \ No newline at end of file diff --git a/swe_smodel.m b/swe_smodel.m index b25e7a8..9bdab81 100644 --- a/swe_smodel.m +++ b/swe_smodel.m @@ -1,7 +1,7 @@ function swe_smodel % Launches batch window containing SwE batch module. % ========================================================================= -% This function prepares (if needed) and launches the batch system with a +% This function prepares (if needed) and launches the batch system with a % job containing the batch module for the specification of % data and design. % ========================================================================= @@ -10,17 +10,17 @@ % Written by Bryan Guillaume % Version Info: $Format:%ci$ $Format:%h$ - % Initiate a job + % Initiate a job if isempty(spm_figure('FindWin','Graphics')) % SPM not running swe_jobman('initcfg') end - + % Launch the batch system with the SwE tab swe_batch % Add the specification module to it swe_jobman('interactive','','swe.smodel'); return - -end \ No newline at end of file + +end diff --git a/swe_splitCovariate.m b/swe_splitCovariate.m index a379dd8..bc64f38 100644 --- a/swe_splitCovariate.m +++ b/swe_splitCovariate.m @@ -1,18 +1,18 @@ function [crossCov, longCov, avgCov] = swe_splitCovariate(cov, subject) % Splits a covariate into pure cross-sectional and longitudinal components. % ========================================================================= -% The split is done first by extracting the within-subect average to build -% the cross-sectional covariate, which is also centered. Second, the -% longitudinal covariate is simply the difference between the original +% The split is done first by extracting the within-subect average to build +% the cross-sectional covariate, which is also centered. Second, the +% longitudinal covariate is simply the difference between the original % covariate and the non-centered cross-sectional covariate. % ========================================================================= % FORMAT: [crossCov, longCov, avgCov] = swe_splitCovariate(cov, subject) % ------------------------------------------------------------------------- % Inputs: % - cov: the covariate to split -% - subject: vector contening the subject numbering +% - subject: vector contening the subject numbering % ------------------------------------------------------------------------- -% Outputs: +% Outputs: % - crossCov: the centered cross-sectional covariate % - longCov: the longitudinal covariate % - avgCrossCov: the average value (scalar) of the time-varying covariate diff --git a/swe_tfce_transform.m b/swe_tfce_transform.m index 14a85a9..58d39da 100644 --- a/swe_tfce_transform.m +++ b/swe_tfce_transform.m @@ -1,7 +1,7 @@ function [tfced] = swe_tfce_transform(img,H,E,C,dh) % This function calculates the TFCE statistical from a Z statistic image. % ========================================================================= -% FORMAT [tfced] = matlab_tfce_transform(img,H,E,C,ndh) +% FORMAT [tfced] = matlab_tfce_transform(img,H,E,C,ndh) % ------------------------------------------------------------------------- % Inputs: % - img the 3D image to be transformed @@ -15,7 +15,7 @@ % % This function was modified from 'matlab_tfce_transform' from MatlabTFCE. % ------------------------------------------------------------------------- -% Additional Octave compatability using SPM functions was added by +% Additional Octave compatability using SPM functions was added by % Tom Maullin (19/09/2018) % Version Info: $Format:%ci$ $Format:%h$ % @@ -56,7 +56,7 @@ % installed. The same computation can be performed with SPM functions but % in much slower time. else - + % find connected components vals = zeros(size(img)); for h = 1:ndh diff --git a/swe_uc_FDR.m b/swe_uc_FDR.m index 6f7a596..e53790d 100644 --- a/swe_uc_FDR.m +++ b/swe_uc_FDR.m @@ -81,15 +81,15 @@ % Copyright (C) 2002-2015 Wellcome Trust Centre for Neuroimaging % Thomas Nichols % Version Info: $Format:%ci$ $Format:%h$ - - + + if (nargin<6), Vm = []; end - + % Set Benjamini & Yeuketeli cV for independence/PosRegDep case %-------------------------------------------------------------------------- cV = 1; - - + + % Load, mask & sort statistic image (if needed) %-------------------------------------------------------------------------- if isstruct(Vs) @@ -110,8 +110,8 @@ Ts = sort(Ts(:)); if STAT ~= 'P', Ts = flipud(Ts); end end - - + + % Calculate p values of image (if needed) %-------------------------------------------------------------------------- if isstruct(Vs) @@ -119,15 +119,15 @@ else Ps = Vs; end - + S = length(Ps); - - + + % Calculate FDR inequality RHS %-------------------------------------------------------------------------- Fi = (1:S)'/S*q/cV; - - + + % Find threshold %-------------------------------------------------------------------------- I = find(Ps<=Fi, 1, 'last' ); @@ -156,4 +156,3 @@ end end end - \ No newline at end of file diff --git a/swe_ui_main.m b/swe_ui_main.m index 6f34aaf..724fbb9 100644 --- a/swe_ui_main.m +++ b/swe_ui_main.m @@ -18,7 +18,7 @@ % ========================================================================= % Author: Tom Maullin (05/09/2018) % Version Info: $Format:%ci$ $Format:%h$ - + %======================================================================= close(findobj(get(0,'Children'),'Tag','SwE Menu')) @@ -48,25 +48,25 @@ 'Rotation',90,... 'VerticalAlignment','bottom','HorizontalAlignment','center',... 'Color',[1 1 1]*.6); - + text(0.5,0.96,'Sandwich Estimator Toolbox',... 'FontName',PF.times,'FontSize',FS(14),... %'FontAngle','Italic',... 'VerticalAlignment','middle','HorizontalAlignment','center',... 'FontWeight','Bold',... 'Color',[1 1 1]*.6); - + text(0.98,0.9,['Version ' swe('ver')],... 'FontName',PF.times,'FontSize',FS(11),'FontAngle','Italic',... 'VerticalAlignment','middle','HorizontalAlignment','right',... 'FontWeight','Bold',... 'Color',[1 1 1]*.6); - + %-Inner Frames %---------------------------------------------------------------------- uicontrol(F,'Style','Frame','Position',[095 005 200 240],... 'BackgroundColor',swe('Colour')); uicontrol(F,'Style','Frame','Position',[105 015 180 220]); - + %-Buttons to launch SwE functions %---------------------------------------------------------------------- uicontrol(F,'String','Specify Model',... @@ -74,7 +74,7 @@ 'CallBack','swe_smodel',... 'Interruptible','on',... 'ForegroundColor','k'); - + uicontrol(F,'String','Run Model',... 'Position',[115 190-(2-1)*55 130 035],... 'CallBack','swe_rmodel',... @@ -85,7 +85,7 @@ 'CallBack','spm_help(''swe_cp.m'')',... 'Interruptible','on',... 'ForegroundColor','b'); - + uicontrol(F,'String','Results',... 'Position',[115 190-(3-1)*55 130 035],... 'CallBack','swe_results',... @@ -96,18 +96,18 @@ 'CallBack','spm_help(''swe_results_ui.m'')',... 'Interruptible','on',... 'ForegroundColor','b'); - + uicontrol(F,'String','Help (Online)',... %'Documentation
  (Online)',... 'Position',[115 190-(4-1)*55 130 035],... 'CallBack','web(''http://www.nisox.org/Software/SwE'')',... 'Interruptible','on',... 'ForegroundColor','k'); - + %-Make visible %---------------------------------------------------------------------- set(F,'Pointer','Arrow','Visible','on'); - - -end \ No newline at end of file + + +end diff --git a/swe_update.m b/swe_update.m index 6b15973..9381b1d 100644 --- a/swe_update.m +++ b/swe_update.m @@ -14,25 +14,25 @@ function swe_update() % Obtain swe version number. vswe=swe('ver'); - + % Look to github for a newer version. url='https://github.com/NISOx-BDI/SwE-toolbox/releases/latest'; [s,stat]=urlread(url); - + % Error if we couldn't contact github. - if ~stat - try + if ~stat + try s = webread(url); catch error('Can''t contact GitHub'); end end - + % Look for latest swe version number. [tok,x]=regexp(s,'Version [0-9.]*','match','tokens'); tok=tok{1}; tok=strrep(tok,'Version ',''); - + % Tell the user whether their version is the newest. if strcmp(tok,vswe) msg = 'Your version of swe is current.'; @@ -42,4 +42,4 @@ function swe_update() if ~nargout, fprintf([blanks(9) msg '\n']); end -end \ No newline at end of file +end diff --git a/swe_vechCovVechV.m b/swe_vechCovVechV.m index cbc1ae0..183fafb 100644 --- a/swe_vechCovVechV.m +++ b/swe_vechCovVechV.m @@ -7,18 +7,18 @@ % % - type: % 1: based on the estimate "II" proposed in Guillaume (2015). -% It accounts partially for missing data (not for a missing data +% It accounts partially for missing data (not for a missing data % bias) and for a small sample bias % 2: based on the estimate "III" proposed in Guillaume (2015). -% It accounts for missing data (included the missing data bias), +% It accounts for missing data (included the missing data bias), % but not for the small sample bias. % 3: based on the estimate "II" proposed in Guillaume (2015), but -% requires no missing data (so accounts for the small sample +% requires no missing data (so accounts for the small sample % bias). -% - covVis: matrix (nCovVis X nVox) containing vech(\hat V) for several -% voxels -% - dofMat: for type 1 and 2, matrix (nCovVis X nCovVis) containing -% degrees of freedom information (precalculated in swe_cp.m). +% - covVis: matrix (nCovVis X nVox) containing vech(\hat V) for several +% voxels +% - dofMat: for type 1 and 2, matrix (nCovVis X nCovVis) containing +% degrees of freedom information (precalculated in swe_cp.m). % for type 3, scalar. % See Guillaume (2015)for more information. % ------------------------------------------------------------------------- @@ -27,7 +27,7 @@ % - vechCovVechV: vech(Cov(vech(\hat V))) as given in Guillaume (2015). % ========================================================================= % -% Note: there are probably ways to compute this quicker, but this will do +% Note: there are probably ways to compute this quicker, but this will do % for now. % By Bryan Guillaume % Version Info: $Format:%ci$ $Format:%h$ @@ -42,28 +42,28 @@ factor = dofMat / (dofMat - 1) / (2* dofMat + 1); end -it = 0; +it = 0; for i = 1:nCovVis for ii = i:nCovVis it = it + 1; - + ind11 = rIndex == rIndex(i) & cIndex == rIndex(ii); ind12 = rIndex == rIndex(i) & cIndex == cIndex(ii); ind22 = (rIndex == cIndex(i) & cIndex == cIndex(ii)) | (cIndex == cIndex(i) & rIndex == cIndex(ii)); ind21 = (rIndex == cIndex(i) & cIndex == rIndex(ii)) | (cIndex == cIndex(i) & rIndex == rIndex(ii)); indi1 = rIndex == rIndex(i) & cIndex == rIndex(i); indi2 = rIndex == cIndex(i) & cIndex == cIndex(i); - indii1 = rIndex == rIndex(ii) & cIndex == rIndex(ii); + indii1 = rIndex == rIndex(ii) & cIndex == rIndex(ii); indii2 = rIndex == cIndex(ii) & cIndex == cIndex(ii); - + switch type - + case 1 tmp = 1 + 2 * dofMat(i,ii) * dofMat(ind11,ind22) * dofMat(ind12,ind21) - ... dofMat(i,ii) * dofMat(ind11,ind22) - ... dofMat(ind11,ind22) * dofMat(ind12,ind21) - ... dofMat(i,ii) * dofMat(ind12,ind21); - + vechCovVechV(it,:) = dofMat(i,ii) / tmp * (... (2 * dofMat(ind11,ind22) * dofMat(ind12,ind21) - dofMat(ind11,ind22) - dofMat(ind12,ind21)) * covVis(i,:) .* covVis(ii,:) + ... (1 - dofMat(ind12,ind21)) * covVis(ind11,:) .* covVis(ind22,:) + ... diff --git a/test/runTest.m b/test/runTest.m index 50b3da8..64a43d3 100644 --- a/test/runTest.m +++ b/test/runTest.m @@ -60,7 +60,7 @@ error(['Test ' testname ' has failed.']) else disp('All tests pass!!') - end + end disp(sprintf('==============================================================\n')) % Teardown method @@ -75,8 +75,8 @@ function testSetup(pOrWb, inferenceType, tOrF, matNiiGiiOrCii) ls; addpath('/swe'); addpath('/swe/test'); - - % Run teardown method just in case some of the files from the previous + + % Run teardown method just in case some of the files from the previous % run managed to get cached. testTearDown(pOrWb, inferenceType, tOrF, matNiiGiiOrCii); @@ -110,30 +110,30 @@ function generateData(pOrWb, inferenceType, tOrF, matNiiGiiOrCii) % Define a contrast. load('SwE.mat'); - + load('/swe/test/data/xCon.mat'); - + if strcmp(tOrF, 't') xCon.name = 'T contrast 1'; xCon.STAT = 'T'; - xCon.c = [1 0 0]'; + xCon.c = [1 0 0]'; else xCon.name = 'F contrast 1'; xCon.STAT = 'F'; - xCon.c = eye(3); + xCon.c = eye(3); end SwE.xCon = xCon; - + save('SwE.mat', 'SwE'); - + % Run swe_getSPM(). For the nii % case we use the xSwE object to % avoid user input. if strcmp(matNiiGiiOrCii, 'nii') || strcmp(matNiiGiiOrCii, 'gii') || strcmp(matNiiGiiOrCii, 'cii') - + load('/swe/test/data/xSwE.mat') - + xSwE.swd = SwE.swd; xSwE.STAT = upper(tOrF); xSwE.S = SwE.xVol.S; @@ -143,13 +143,13 @@ function generateData(pOrWb, inferenceType, tOrF, matNiiGiiOrCii) M = SwE.xVol.M(1:3,1:3); VOX = sqrt(diag(M'*M))'; xSwE.VOX = VOX; - + swe_getSPM(xSwE); - + else - + swe_getSPM(); - + end end @@ -157,12 +157,12 @@ function generateData(pOrWb, inferenceType, tOrF, matNiiGiiOrCii) end function design = createDesignFromTemplate(pOrWb, inferenceType, tOrF, matNiiGiiOrCii) - + % Load template design load('/swe/test/data/design_wb_clus_f_mat.mat'); - + design.dir = {['/swe/test/data/test_' pOrWb '_' inferenceType '_' tOrF '_' matNiiGiiOrCii]}; - + if strcmp(matNiiGiiOrCii, 'nii') fileExtension = 'img'; elseif strcmp(matNiiGiiOrCii, 'cii') @@ -182,11 +182,11 @@ function generateData(pOrWb, inferenceType, tOrF, matNiiGiiOrCii) if strcmp(matNiiGiiOrCii, 'cii') design.ciftiAdditionalInfo = struct; - + design.ciftiAdditionalInfo.ciftiGeomFile(1) = struct; design.ciftiAdditionalInfo.ciftiGeomFile(1).brainStructureLabel = 'CIFTI_STRUCTURE_CORTEX_LEFT'; design.ciftiAdditionalInfo.ciftiGeomFile(1).geomFile = '/swe/test/data/gii_input/L.sphere.4k_fs_LR.surf.gii'; - + design.ciftiAdditionalInfo.volRoiConstraint = 1; end @@ -194,9 +194,9 @@ function generateData(pOrWb, inferenceType, tOrF, matNiiGiiOrCii) design.WB = rmfield(design.WB, 'WB_yes'); design.WB.WB_no = 0; else - + if strcmp(inferenceType, 'clus') - + if ~strcmp(matNiiGiiOrCii, 'mat') design.WB.WB_yes.WB_infType.WB_clusterwise.WB_inputType = ... rmfield(design.WB.WB_yes.WB_infType.WB_clusterwise.WB_inputType, 'WB_mat'); @@ -204,7 +204,7 @@ function generateData(pOrWb, inferenceType, tOrF, matNiiGiiOrCii) end elseif strcmp(inferenceType, 'tfce') - + design.WB.WB_yes.WB_infType = rmfield(design.WB.WB_yes.WB_infType, 'WB_clusterwise'); design.WB.WB_yes.WB_infType.WB_TFCE = struct; design.WB.WB_yes.WB_infType.WB_TFCE.WB_TFCE_E = 0.5; @@ -228,7 +228,7 @@ function generateData(pOrWb, inferenceType, tOrF, matNiiGiiOrCii) end function mapsEqual = verifyMapsUnchanged(pOrWb, inferenceType, tOrF, matNiiGiiOrCii) - + % List all files for testing if strcmp(matNiiGiiOrCii, 'nii') files = ls("*.nii"); @@ -238,7 +238,7 @@ function generateData(pOrWb, inferenceType, tOrF, matNiiGiiOrCii) filetype = 'gii'; elseif strcmp(matNiiGiiOrCii, 'cii') files = ls("*.nii"); - filetype = 'cii'; + filetype = 'cii'; else files = ls("swe_*.mat"); filetype = 'mat'; @@ -281,19 +281,19 @@ function generateData(pOrWb, inferenceType, tOrF, matNiiGiiOrCii) % Remove NaN and zero values file = file(~isnan(file) & file ~= 0); gt_file = gt_file(~isnan(gt_file) & gt_file ~= 0); - + % Check whether the remaining values are equal. % % Footnote: There is some form of machine tolerance % error with the mat cases. Even with all seeds % reset, differences in voxel values of around ~e16 - % ~e16 can occur occasionally on runs. These errors - % can propogate and grow as large as ~e8. To + % ~e16 can occur occasionally on runs. These errors + % can propogate and grow as large as ~e8. To % counter this, in these cases we just look to see % if we are within e-10 of the ground truth. - % Suspected cause: The `beta` and `betainc` - % functions which are only run in these cases and - % are built on old fortran numeric approximations + % Suspected cause: The `beta` and `betainc` + % functions which are only run in these cases and + % are built on old fortran numeric approximations % in octave. (Tom Maullin 07/06/2018) if strcmp(matNiiGiiOrCii, 'mat') result = ~any(abs(file-gt_file) > 10^(-7)); @@ -302,35 +302,35 @@ function generateData(pOrWb, inferenceType, tOrF, matNiiGiiOrCii) result = ~any(file~=gt_file) indexWrong = file~=gt_file; end - + % Useful for debugging. if ~result disp('Length file: ') disp(size(file)) - + disp('Length gt_file: ') disp(size(gt_file)) - + disp('size disagreement values: ') disp(sum(indexWrong)) - + disp('sum of disagreement (file): ') disp(sum(file(indexWrong))) - + disp('sum of disagreement (gt_file): ') disp(sum(gt_file(indexWrong))) - + disp('disagreement values (file)') d1 = file(indexWrong); disp(sprintf('%.9f', d1(1))) - + disp('disagreement values (gt_file)') d2 = gt_file(indexWrong); disp(sprintf('%.9f', d2(1))) - + disp('diff') disp(d2(1)-d1(1)) - + end if result @@ -351,7 +351,7 @@ function generateData(pOrWb, inferenceType, tOrF, matNiiGiiOrCii) end function testTearDown(pOrWb, inferenceType, tOrF, matNiiGiiOrCii) - + % Delete all files from this run. delete('swe_*'); delete('SwE*'); @@ -361,4 +361,4 @@ function testTearDown(pOrWb, inferenceType, tOrF, matNiiGiiOrCii) rename('xSwE_orig.mat', 'xSwE.mat'); end -end \ No newline at end of file +end From 697ec8bc09b5da6a3ffdfc56365a2b572992b1f0 Mon Sep 17 00:00:00 2001 From: nicholst Date: Thu, 18 Jun 2020 09:55:44 +0100 Subject: [PATCH 2/2] Clean up @swe_DataType indenting --- @swe_DataType/display.m | 6 +++--- @swe_DataType/eq.m | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/@swe_DataType/display.m b/@swe_DataType/display.m index b28c2db..b0a4674 100644 --- a/@swe_DataType/display.m +++ b/@swe_DataType/display.m @@ -1,9 +1,9 @@ function display(obj) - % Display a swe_DataType object - % ========================================================================= + % Display a swe_DataType object + % ========================================================================= % Bryan Guillaume % Version Info: $Format:%ci$ $Format:%h$ - display(get(obj)); + display(get(obj)); end diff --git a/@swe_DataType/eq.m b/@swe_DataType/eq.m index 8be828f..d2131d1 100644 --- a/@swe_DataType/eq.m +++ b/@swe_DataType/eq.m @@ -1,9 +1,9 @@ function out = eq(a,b) - % Overload the '==' operator for we_DataType objects + % Overload the '==' operator for we_DataType objects % ========================================================================= % Bryan Guillaume - % Version Info: $Format:%ci$ $Format:%h$ + % Version Info: $Format:%ci$ $Format:%h$ - out = strcmp(a.value, b.value); + out = strcmp(a.value, b.value); end