diff --git a/Addpaths.m b/Addpaths.m
new file mode 100644
index 0000000..575ff28
--- /dev/null
+++ b/Addpaths.m
@@ -0,0 +1,15 @@
+function param = Addpaths
+ % Path to the spm12
+ addpath(genpath('/Users/anjalitarun/Documents/SPM/spm12'));
+
+ addpath(genpath(pwd));
+
+
+ % Code base path
+ param.DSIstudio = '/Users/anjalitarun/Documents/Softwares/DSI_studio/Contents/MacOS';
+
+ % Database path
+ param.HCPDatapath = '/Users/anjalitarun/Downloads/HCP';
+
+
+end
\ No newline at end of file
diff --git a/BrainGraph/main_braingraph.m b/BrainGraph/main_braingraph.m
new file mode 100644
index 0000000..408ebd0
--- /dev/null
+++ b/BrainGraph/main_braingraph.m
@@ -0,0 +1,44 @@
+%% Main code for constructing the voxel-level brain grid
+
+% Created by: Anjali Tarun
+% Date created: 21 February 2018
+% Based from the original code of Hamid Behjat, Martin Larsson and
+% David Abramian
+
+
+% Changes made in the original ODF design of Itturia Medina and HB's code:
+% -- to use whole brain volume instead of WM only
+% -- incorporated structural boost using QA anisotropy
+% -- removes outlier brain voxels to aid in convergence of the
+% eigendecomposition of the Laplacian
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function [WB, G] = main_braingraph(param)
+
+ hcp_root = param.HCPDatapath; % folder in which your HCP subjects are stored in
+
+ dsi_root = param.DSIstudio; % root to dsi_studio's executable
+
+ subjID = param.subject;
+
+ % Settings.
+ odfPow = param.odfPow; % ODF power.
+ neighborhood = param.neighborhood; % 3 or 5 [This is for white matter; gray matter always 3]
+ N_fibers = param.N_fibers; % dsi studio's and M&L's default: 5
+
+ % Preprocess.
+ G = mdh_preprocess(subjID,hcp_root);
+
+ % Whole brain graph.
+ [A_wb,indices_wb,G,O] = mdh_adjacencymatrix_wb(G,odfPow,neighborhood,...
+ dsi_root,param,'N_fibers',N_fibers);
+
+% save(fullfile(G.f.graph,'G_wb.mat'),'G');
+
+ WB.indices = indices_wb;
+ WB.A = A_wb;
+ WB.O = O;
+
+
+end
diff --git a/BrainGraph/mdh/da_corregister_reslice.m b/BrainGraph/mdh/da_corregister_reslice.m
new file mode 100755
index 0000000..f73c5fc
--- /dev/null
+++ b/BrainGraph/mdh/da_corregister_reslice.m
@@ -0,0 +1,30 @@
+function [fileOutput] = da_corregister_reslice( ref, source, interp, prefix )
+%DA_RESIZE Reslice a volume using SPM
+% ref: path of reference volume
+% source: path of volume to reslice
+% interp: 0 for nearest neighbor, 1 for trilinear interpolation
+% prefix: prefix added to the name of the resliced file
+
+if nargin < 3
+ interp = 0;
+end
+if nargin < 4
+ prefix = 'r';
+end
+
+if interp ~=0 && interp ~=1
+ error('Interpolation must be either 0 or 1')
+end
+
+job.ref = {strcat(ref,',1')};
+job.source = {strcat(source,',1')};
+job.roptions.interp = interp;
+job.roptions.wrap = [0 0 0];
+job.roptions.mask = 0;
+job.roptions.prefix = prefix;
+
+out = spm_run_coreg(job);
+fileOutput = strsplit(out.rfiles{1},',');
+fileOutput = fileOutput{1};
+end
+
diff --git a/BrainGraph/mdh/da_resize.m b/BrainGraph/mdh/da_resize.m
new file mode 100755
index 0000000..5909e5b
--- /dev/null
+++ b/BrainGraph/mdh/da_resize.m
@@ -0,0 +1,56 @@
+function [ fileOutput ] = da_resize( fileInput, fileRef )
+%DA_RESIZE_AND_CLEAN Resize a nifti file to match the size of a reference
+%file
+% This function resizes an input nifti file to match the size of a
+% reference file.
+%
+% Inputs:
+% fileInput: full address of nifti file to resize.
+% fileRef: full address of nifti file used as reference for size and
+% coregistration.
+%
+% Outputs:
+% outputFile: full address of created file, with naming convention
+% inputFileName_referenceVoxelSize.nii, and created in the same
+% folder as the input file.
+%
+% David Abramian, May 2017
+
+% Check input for errors
+if ~exist(fileInput,'file')
+ error('Input file does not exist')
+end
+
+if ~exist(fileRef,'file')
+ error('Reference file does not exist')
+end
+
+[dirInput, name, ext] = fileparts(fileInput);
+if ~strcmp(ext, '.nii')
+ error('Input must be a nifti file')
+end
+
+[~, ~, ext] = fileparts(fileRef);
+if ~strcmp(ext, '.nii')
+ error('Reference must be a nifti file')
+end
+
+
+% Determine output voxel size
+hRef = spm_vol(fileRef);
+outputVoxelSize = abs(hRef(1).mat(1,1));
+
+if rem(outputVoxelSize,1) == 0
+ strOutputVoxelSize = num2str(outputVoxelSize,'%.1f');
+else
+ strOutputVoxelSize = num2str(outputVoxelSize);
+end
+
+% Name of output file
+fileOutput = fullfile(dirInput, [name,'_', strOutputVoxelSize, ext]);
+
+% Resize file
+da_corregister_reslice(fileRef, fileInput, 0, 'rnn_');
+movefile(fullfile(dirInput, ['rnn_', name, ext]), fileOutput);
+
+
diff --git a/BrainGraph/mdh/external/argselectAssign.m b/BrainGraph/mdh/external/argselectAssign.m
new file mode 100644
index 0000000..9ffeecd
--- /dev/null
+++ b/BrainGraph/mdh/external/argselectAssign.m
@@ -0,0 +1,47 @@
+% argselectAssign : Assign variables in calling workspace
+%
+% function argselectAssign(variable_value_pairs)
+%
+% Inputs :
+% variable_value_pairs is a cell list of form
+% 'variable1',value1,'variable2',value2,...
+% This function assigns variable1=value1 ... etc in the *callers* workspace
+%
+% This is used at beginning of function to simulate keyword argument
+% passing. Typical usage is
+%
+% argselectAssign(control_params);
+% argselectCheck(control_params,varargin);
+% argselectAssign(varargin);
+%
+% where control_params is a cell list of variable,value pairs containing
+% the default parameter values.
+%
+% See also argselectCheck
+%
+% Author : David K. Hammond, EPFL LTS2
+% Date : December, 2007
+% Project : common utilities
+
+% This file is part of the SGWT toolbox (Spectral Graph Wavelet Transform toolbox)
+% Copyright (C) 2010, David K. Hammond.
+%
+% The SGWT toolbox is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% The SGWT toolbox is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with the SGWT toolbox. If not, see .
+
+function argselectAssign(variable_value_pairs)
+for j =1:2:length(variable_value_pairs)
+ pname=variable_value_pairs{j};
+ pval=variable_value_pairs{j+1};
+ assignin('caller',pname,pval);
+end
diff --git a/BrainGraph/mdh/external/argselectCheck.m b/BrainGraph/mdh/external/argselectCheck.m
new file mode 100644
index 0000000..c79c5fc
--- /dev/null
+++ b/BrainGraph/mdh/external/argselectCheck.m
@@ -0,0 +1,54 @@
+% argselectCheck : Check if control parameters are valid
+%
+% function argselectCheck(control_params,varargin_in)
+%
+% Inputs:
+% control_params and varargin_in are both cell arrays
+% that are lists of pairs 'name1',value1,'name2',value2,...
+%
+% This function checks that every name in varargin_in is one of the names
+% in control_params. It will generate an error if not, otherwise it will
+% do nothing
+%
+% This is used at beginning of function to simulate keyword argument
+% passing. Typical usage is
+%
+% argselectAssign(control_params);
+% argselectCheck(control_params,varargin);
+% argselectAssign(varargin);
+%
+% where control_params is a cell list of variable,value pairs containing
+% the default parameter values.
+%
+% See also argselectAssign
+%
+% Author : David K. Hammond, EPFL LTS2
+% Date : December, 2007
+% Project : common utilities
+
+% This file is part of the SGWT toolbox (Spectral Graph Wavelet Transform toolbox)
+% Copyright (C) 2010, David K. Hammond.
+%
+% The SGWT toolbox is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% The SGWT toolbox is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with the SGWT toolbox. If not, see .
+
+function argselectCheck(control_params,varargin_in)
+[param_names{1:(length(control_params)/2)}]=control_params{1:2:end};
+% check that passed in control parameters are valid
+for j=1:2:length(varargin_in)
+ if(isempty(strmatch(varargin_in{j},param_names)))
+ error(['Invalid control parameter : ',varargin_in{j},...
+ ' Valid options are',sprintf('\n'),...
+ sprintf('%s \n',param_names{:})]);
+ end
+end
\ No newline at end of file
diff --git a/BrainGraph/mdh/hb_adjacencymatrix_gm.m b/BrainGraph/mdh/hb_adjacencymatrix_gm.m
new file mode 100644
index 0000000..99fb750
--- /dev/null
+++ b/BrainGraph/mdh/hb_adjacencymatrix_gm.m
@@ -0,0 +1,109 @@
+function [A,indices,G,A_a,A_ac,A_acp] = hb_adjacencymatrix_gm(mask,conn,connComp,varargin)
+% HB_ADJACENCYMATRIX_GM Constructs gray matter graph A matrix.
+%
+% Inputs:
+% mask: 3D volume, or its full address (/Users/.../xxx.nii).
+% conn: 3D neighbourhood: 6,18 or 26.
+% connComp: Number of connected components in structure.
+%
+% Outputs:
+% A: adjacency matrix
+% indices: corresponding indices in 'mask'
+% Aa:
+% Aac:
+% Aacp:
+%
+% Contributers:
+% Hamid Behjat
+% Martin Larsson
+% March 2017
+
+cntr = {
+ 'weight','no',...
+ 'pialFiles',[],...
+ 'parallelize',0,...
+ 'sourceFile',[],...
+ 'G',[]};
+argselectAssign(cntr);
+argselectCheck(cntr,varargin);
+argselectAssign(varargin);
+
+fprintf('\n.Constructing GM graph adjacency matrix..\n')
+if ~isempty(pialFiles)
+ if ischar(mask)
+ sourceFile_h = spm_vol(mask);
+ elseif ~isempty(sourceFile)
+ sourceFile_h = spm_vol(sourceFile);
+ else
+ error('source file needs to be specified to extract header.')
+ end
+end
+
+% compute adjacency mtrix (A)
+fprintf('..Determining adjacencies..\n')
+[A,indices,mask] = hb_compute_adjacency(mask,conn,...
+ 'weight',weight,'sourceFile',sourceFile,'sS',0);
+fprintf('..done.\n')
+
+A_a = A; % A adjacency
+indices_a = indices;
+
+% Clean A.
+[A,d,mask] = ml_clean_adjacency(A,connComp,mask);
+%fprintf([num2str(numel(indices)-numel(d)),' nodes removed.\n'])
+
+A_ac = A; % A adjacency-cleaned
+
+indices = indices(d);
+
+% sanity check
+if setdiff(indices,find(mask))
+ error('something is fishy..')
+end
+
+if ~isempty(pialFiles)
+ pials = ml_batch_gifti(pialFiles,sourceFile_h.mat);
+ A_tmp = A;
+
+ % Prune A.
+ fprintf('..Pruning adjacency matrix..\n')
+ if conn==6
+ A = ml_prune_adjacency(A,mask,pials,...
+ 'parallelize',parallelize,'conn6',1);
+ else
+ A = ml_prune_adjacency(A,mask,pials,...
+ 'parallelize',parallelize);
+ end
+ fprintf([' ',num2str((numel(find(A_tmp))-numel(find(A)))/2),...
+ ' edges removed.\n'])
+ fprintf('..done.\n')
+
+ A_acp = A; % A adjacency-cleaned-pruned
+
+ % Clean A.
+ fprintf('..Cleaning adjacency matrix..\n')
+ [A,d,mask] = ml_clean_adjacency(A,connComp,mask);
+ fprintf([' ',num2str(numel(indices)-numel(d)),' nodes removed.\n'])
+ fprintf('..done.\n')
+
+ indices = indices(d);
+
+ % sanity check
+ if setdiff(indices,find(mask))
+ error('something is fishy..')
+ end
+end
+
+if ~isempty(G)
+ G.N_gm = size(A,1);
+ G.A_gm = A;
+ G.indices_gm = indices;
+ G.A_gm_nonpruned = A_a;
+ G.indices_gm_nonpruned = indices_a;
+end
+
+A_gm = A; %#ok
+indices_gm = indices; %#ok
+save(fullfile(G.f.graph,'A_gm.mat'),'A_gm');
+save(fullfile(G.f.graph,'indices_gm.mat'),'indices_gm');
+fprintf('done.\n')
diff --git a/BrainGraph/mdh/hb_adjacencymatrix_gmwm.m b/BrainGraph/mdh/hb_adjacencymatrix_gmwm.m
new file mode 100644
index 0000000..ed29a0a
--- /dev/null
+++ b/BrainGraph/mdh/hb_adjacencymatrix_gmwm.m
@@ -0,0 +1,87 @@
+function [A,indices,G] = hb_adjacencymatrix_gmwm(G)
+%HB_ADJACENCYMATRIX_GMWM Constructs GM+WM graph A matrix.
+%
+%
+% Author:
+% Hamid Behjat
+% Feb 2018.
+
+fprintf('\n.Constructing GM+WM graph adjacency matrix..\n')
+
+[A,indices] = hb_compute_adjacency(G.f.gmwm_125,26,'sS',0);
+
+[~,I_gm] = intersect(indices,G.indices_gm);
+[~,I_wm] = intersect(indices,G.indices_wm);
+
+d1 = setdiff(indices,G.indices_gm);
+d2 = setdiff(d1,G.indices_wm);
+[~,I_extra] = intersect(indices,d2);
+
+% sanity checks
+if length(I_gm)~= G.N_gm, error('fishy..'); end
+if length(I_wm)~= G.N_wm, error('fishy..'); end
+if (G.N_gm+G.N_wm+length(I_extra))~=length(indices), error('fishy..'); end
+
+G.N_gmwm = G.N_gm+G.N_wm;
+
+d1 = nnz(A); %#ok
+A(I_extra,:) = 0;
+A(:,I_extra) = 0;
+
+% Find GM-WM boundary vertices.
+d = A;
+d(I_gm,I_gm) = 0;
+d(I_wm,I_wm) = 0;
+[d,~] = ind2sub(size(d),find(d));
+d = unique(d);
+d = indices(d);
+[~,I_gm_boundary] = intersect(d,G.indices_gm);
+[~,I_wm_boundary] = intersect(d,G.indices_wm);
+if ~isempty(intersect(I_gm_boundary,I_wm_boundary)), error('fishy...'); end
+G.indices_gm_boundary = d(I_gm_boundary);
+G.indices_wm_boundary = d(I_wm_boundary);
+
+% Restore GM edges.
+d2 = nnz(A); %#ok
+A(I_gm,I_gm) = G.A_gm;
+d3 = nnz(A); %#ok
+
+% Restore WM edges.
+A(I_wm,I_wm) = G.A_wm;
+d4 = nnz(A); %#ok
+
+% Clean-up A.
+c = find(~sum(A,1));
+A(:,c) = [];
+A(c,:) = [];
+indices(c)=[];
+
+% Store.
+G.A_gmwm = A;
+G.indices_gmwm = indices;
+fprintf('.done.\n')
+
+% Stuff.
+d5 = nnz(A);
+d6 = nnz(G.A_gm)/d5;
+d7 = nnz(G.A_wm)/d5;
+d8 = 1-d6-d7;
+
+% Summary.
+disp(' ');
+disp(' ');
+disp('*** GM+WM graph facts ***')
+disp(['percentage of nodes being GM: ',num2str(100*G.N_gm/G.N_gmwm),'%']);
+disp(['percentage of nodes being WM: ',num2str(100*G.N_wm/G.N_gmwm),'%']);
+disp(' ');
+disp(['percentage of GM nodes at GM-WM boundary: ',...
+ num2str(100*length(G.indices_gm_boundary)/G.N_gm),'%']);
+disp(['percentage of WM nodes at GM-WM boundary: ',...
+ num2str(100*length(G.indices_wm_boundary)/G.N_wm),'%']);
+disp(' ');
+disp(['percentage of edges being intra GM : ',num2str(100*d6),'%']);
+disp(['percentage of edges being intra WM : ',num2str(100*d7),'%']);
+disp(['percentage of edges being inter GM-WM: ',num2str(100*d8),'%']);
+disp(' ');
+disp(['percentage of edges pruned: ',num2str(round(1e4*((d2-d3)/d2))*1e-4*100),'%']);
+
diff --git a/BrainGraph/mdh/hb_compute_adjacency.m b/BrainGraph/mdh/hb_compute_adjacency.m
new file mode 100644
index 0000000..872b208
--- /dev/null
+++ b/BrainGraph/mdh/hb_compute_adjacency.m
@@ -0,0 +1,257 @@
+function [A,indices,mask] = hb_compute_adjacency(mask,conn,varargin)
+% Computes the adjacency matrix for a given 3D (i.e. volume) mask.
+%
+% mask: input mask to compute its adjacency matrix (A).
+% A will be a square, symmetric matrix. Its dimension is equal
+% to the noumber of non-zero elements in the input mask.
+%
+% The input mask should be a volume. For an image, append
+% the image by two zero plane images to convert it to a volume.
+%
+% * 'mask' can also be the full adress to a .nii file, which will
+% then be loaded.
+%
+% conn: neighbourhood connectivity to be considered when compuing
+% the edges
+%
+% - Optional inputs -
+% 'weight': weighted or non-weighted A. Default: 'no'.
+% For an example weigted option us: 'yes_power_additive'
+%
+% 'pow' : the power used in defining weigting for 'yes_power_additive'
+% weigting option.
+%
+% Output:
+% A : Adjacency matrix in sparse format; i.e. only non zeros elements are saved.
+% To view the full matrix, use the 'full' command.
+%
+% ind_rmd : Removed indices after pruning. if no pruning done, [] returned.
+%
+%
+% Examples:
+%
+% Ex.1 Computes non-weighted A, with neighbourhood connectivity 26:
+%
+% A=compute_adjacency(mask,26);
+%
+% Ex.2 Computes non-weighted A, with neighbourhood connectivity 6:
+%
+% A=compute_adjacency(mask,6);
+%
+% Ex.3 The same as in #1:
+%
+% A=compute_adjacency(mask,26,'weight','no');
+%
+% Ex.4 The same as in #1, but with weighted edges, using a power law:
+%
+% A=compute_adjacency(mask,26,'weight','yes_power_additive','pow',3);
+%
+% Ex.5 Define a mask and compute its A matrix:
+%
+% mask = zeros(3,3,3);
+% mask(:,:,1) = [1 1 1; 0 1 0; 0 0 0];
+% mask(:,:,2) = [0 0 0; 0 1 0; 0 0 0];
+% mask(:,:,3) = [0 0 0; 0 1 0; 1 0.5 0.5];
+% A_v1 = compute_adjacency(mask,26);
+% A_v2 = compute_adjacency(mask,6);
+% A_v3 = compute_adjacency(mask,26,'weight','yes_power_additive');
+% full(A_v1)
+% full(A_v2)
+% full(A_v3)
+%
+% Ex.6 Convert a 2D mask to 3D and compute its A matrix:
+%
+% Im = [1 1 1 1 1; 0 0 0 1 0; 0 0 1 0 0; 0 1 0 0 0; 1 1 0.5 0.5 0.5]
+% mask = zeros(5,5,3);
+% mask(:,:,1) = zeros(5,5);
+% mask(:,:,2) = Im;
+% mask(:,:,3) = zeros(5,5);
+% A_v1 = compute_adjacency(mask,26);
+% A_v2 = compute_adjacency(mask,6);
+% A_v3 = compute_adjacency(mask,26,'weight','yes_power_additive','pow',5);
+% full(A_v1)
+% full(A_v2)
+% full(A_v3)
+%--------------------------------------------
+% Author: Hamid Behjat
+%
+% Biomedical Signal Processing Group,
+% Dept. of Biomedical Engineering,
+% Lund University, Sweden
+%
+% Oct 2016 / Feb 2018
+%
+% Acknowledgment:
+% The initial version of part of this function was developed by
+% Elena Najdenovska and Nora Leonardi, 2011.
+% ************************************************
+
+control_params = {
+ 'weight','no',...
+ 'pow',5,...
+ 's',[],...
+ 'distMetric','Eucl',...
+ 'sw',1,...
+ 'indGM',[],...
+ 'sourceFile',[],... % file associated to the input mask
+ 'sS',1};
+argselectAssign(control_params);
+argselectCheck(control_params,varargin);
+argselectAssign(varargin);
+
+if ischar(mask)
+ if sS, fprintf(['- source file: ',mask,'\n']); end
+ mask = spm_read_vols(spm_vol(mask));
+elseif ~isempty(sourceFile)
+ if sS, fprintf(['- source file: ',sourceFile,'\n']); end
+end
+
+dim = size(mask);
+
+N = numel(mask);
+
+indices = find(mask);
+
+[alli,allj,allk] = ind2sub(dim,indices);
+
+switch conn
+ case 6
+ nN = 3;
+ case 18
+ nN = 9;
+ case 26
+ nN = 13;
+ otherwise
+ error('undefined 3D connectivity neighbourhood!')
+end
+
+ci = repmat(alli,nN,1);
+cj = repmat(allj,nN,1);
+ck = repmat(allk,nN,1);
+
+switch conn
+ case 6
+ ni = [alli ; alli+1; alli ];
+ nj = [allj+1; allj ; allj ];
+ nk = [allk ; allk ; allk+1];
+ case 18
+ ni = [alli ;alli+1;alli+1;alli+1;alli ;alli ;alli ;alli+1;alli+1];
+ nj = [allj+1;allj ;allj-1;allj+1;allj ;allj+1;allj+1;allj ;allj ];
+ nk = [allk ;allk ;allk ;allk ;allk+1;allk-1;allk+1;allk-1;allk+1];
+ case 26
+ ni = [alli;alli+1;alli+1;alli+1;alli;alli;alli;alli+1;alli+1;alli+1;alli+1;alli+1;alli+1];
+ nj = [allj+1;allj;allj-1;allj+1;allj;allj+1;allj+1;allj;allj;allj-1;allj+1;allj+1;allj-1];
+ nk = [allk;allk;allk;allk;allk+1;allk-1;allk+1;allk-1;allk+1;allk-1;allk-1;allk+1;allk+1];
+end
+
+maskZ = cat(2,mask,zeros(dim(1),1,dim(3)));
+maskZ = cat(1,maskZ,zeros(1,dim(2)+1,dim(3)));
+maskZ = cat(3,maskZ,zeros(dim(1)+1,dim(2)+1,1));
+
+valid = (ni>=1 & ni<=dim(1) & nj>=1 & nj<=dim(2) & nk>=1 & nk<=dim(3));
+ni = ni(valid);
+nj = nj(valid);
+nk = nk(valid);
+ci = ci(valid);
+cj = cj(valid);
+ck = ck(valid);
+
+tt = sub2ind(size(maskZ),ni,nj,nk);
+ee = maskZ(tt);
+valid = logical(ee);
+ni = ni(valid);
+nj = nj(valid);
+nk = nk(valid);
+ci = ci(valid);
+cj = cj(valid);
+ck = ck(valid);
+
+cInd = sub2ind(dim,ci,cj,ck);
+nInd = sub2ind(dim,ni,nj,nk);
+
+switch weight
+
+ case 'no'
+ % non-weighted adjacency matrix, as used in Behjat et al., 2015 & 2016.
+ if sS, fprintf('- weigting: none.\n'); end
+
+ A = sparse([cInd,nInd],[nInd,cInd],ones(1,2*numel(ni)),N,N);
+
+ case 'yes_power_additive'
+ % Weighted adjacency matrix, as used in Behjat et al., 2013 & 2014.
+ % The GM probability of nodes connected through an edge are used
+ % to defined a weight for the edge.
+
+ pro = 1;
+ p = (mask(cInd).*mask(nInd)*pro).^pow;
+ A = sparse([cInd,nInd],[nInd,cInd],[p,p],N,N);
+
+ case 'yes_power_subtractive'
+
+ pro = 1;
+ p = (mask(cInd).*mask(nInd)*pro).^(-pow);
+ A = sparse([cInd,nInd],[nInd,cInd],[p,p],N,N);
+
+ case 'yes_Nora'
+
+ d = mask(cInd)-mask(nInd);
+ if isempty(s)
+ s = mean(abs(d))*sw;
+ end
+
+ if strcmp(distMetric,'Gaus')
+ sim = exp(-d.^2/(2*s^2));
+ else
+ sim = 1./(1+d.^2/(2*s^2));
+ end
+ A = sparse([cInd,nInd],[nInd,cInd],[sim,sim],N,N);
+
+
+ case 'yes_distance'
+ if sS, fprintf('- weigting: inverse Euclidean ditance.\n'); end
+
+ [aha1, aha2, aha3] = ind2sub(dim, cInd);
+ [nababa1, nababa2, nababa3] = ind2sub(dim, nInd);
+ d = sqrt(sum(([aha1-nababa1,aha2-nababa2,aha3-nababa3]).^2,2));
+
+ if isempty(s)
+ s = mean(abs(d))*sw;
+ end
+
+ if strcmp(distMetric,'Gaus')
+ sim = exp(-d.^2/(2*s^2)); % gaussian similarity
+ elseif strcmp(distMetric,'adjEucl')
+ sim = 1./(1+d.^2/(2*s^2)); % adjusted Euclidean similarity
+ elseif strcmp(distMetric,'Eucl')
+ sim = 1./d; % Euclidean similarity
+ end
+ A = sparse([cInd,nInd],[nInd,cInd],[sim,sim],N,N);
+
+ case 'yesIndexed' % Oct 2016
+
+ d1 = mask(cInd)-mask(nInd);
+ d2 = find(d1);
+ p = zeros(size(d1));
+ p(d2) = -1; %#ok
+ A=sparse([cInd,nInd],[nInd,cInd],[p,p],N,N);
+end
+A = remove_empty_rows_cols(A,indices,'maintain_inds');
+end
+
+function [X,inds,inds_rmd,c] = remove_empty_rows_cols(X,inds,removeType)
+
+switch removeType
+ case 'maintain_inds'
+ c=find(~sum(X,1));
+ c=c(~ismember(c,inds));
+ X(:,c)=[];
+ X(c,:)=[];
+ case 'update_inds'
+ c=find(~sum(X,1));
+ X(:,c)=[];
+ X(c,:)=[];
+ inds_rmd = inds(c);
+ inds(c)=[];
+end
+end
+
diff --git a/BrainGraph/mdh/hb_get_dir.m b/BrainGraph/mdh/hb_get_dir.m
new file mode 100644
index 0000000..8ae45a2
--- /dev/null
+++ b/BrainGraph/mdh/hb_get_dir.m
@@ -0,0 +1,112 @@
+function dirr = hb_get_dir(dirCatgory,varargin)
+%
+%
+%
+% Requires:
+% hb_get_tag.m
+% hb_get_name.m
+%
+% Hamid Behjat
+% March-July 2017
+
+cntr = {
+ 'baseDir',[],...
+ 'subjectID',[],...
+ 'subjectDir',[],... * for 'local_graphs'
+ 'opts',[],...
+ 'G_num',[],...
+ 'iG',[],...
+ 'N_regions',[],...
+ 'ribbonDir',[]};
+argselectAssign(cntr);
+argselectCheck(cntr,varargin);
+argselectAssign(varargin);
+
+switch dirCatgory
+ case 'hcp_t1w_mri'
+ dirr = [baseDir,filesep,num2str(subjectID),filesep,'T1w'];%[12Feb2018] ,filesep,num2str(subjectID),filesep,'mri'];
+
+ % tempDir = [dataDirectory,filesep,num2str(tempIDs(iSubject)),filesep,'T1w',filesep,num2str(tempIDs(iSubject)),filesep,'mri'];
+ case 'hcp_t1w_surf'
+ dirr = [baseDir,filesep,num2str(subjectID),filesep,'T1w',filesep,num2str(subjectID),filesep,'surf'];
+
+ case 'hcp_mni'
+ dirr = [baseDir,filesep,num2str(subjectID),filesep,'MNINonLinear'];
+
+ case 'bjt'
+ dirr = [ribbonDir,filesep,'bjt'];
+
+ case {'local_graphs','local_graphs_lr','local_graphs_nol'}
+ % nol: non-overlapping
+
+ dirr = [...
+ subjectDir,...
+ filesep,...
+ hb_get_name('vol_clusters_long',opts,...
+ 'N_regions',N_regions,'justName',1)];
+
+ if any(strcmp(dirCatgory,{'local_graphs','local_graphs_nol'}))
+ dirr = [dirr,filesep,hb_get_tag('G_size',opts)];
+ elseif strcmp(dirCatgory,'local_graphs_lr')
+ dirr = [dirr,filesep,hb_get_tag('G_size_lr',opts)];
+ end
+
+ if any(strcmp(dirCatgory,{'local_graphs','local_graphs_lr'}))
+ switch opts.localGraphConstructionScheme
+ case 'nearestNeighbsAroundCenter'
+ dirr = [dirr,'nearestNeighbs_'];
+
+ case 'diffuseOutFromCenter'
+ dirr = [dirr,'diffuseOut_'];
+
+ otherwise
+ error('ooopps..')
+ end
+ elseif strcmp(dirCatgory,'local_graphs_nol')
+ dirr = [dirr,'nonOverLapping_'];
+ end
+
+ dummy = hb_get_tag('A_graph_short',opts);
+ dirr = [dirr,dummy(1:end-1)];
+
+ case 'signals'
+ [pathstr,dummy] = fileparts(hb_get_name('mat_local_graph_mini',opts,'subjectDir',subjectDir,'N_regions',N_regions,'iG',iG));
+ dirr = [pathstr,filesep,dummy];
+
+ case 'hcp_mni_fmri'
+ dirr = [baseDir,filesep,num2str(subjectID),filesep,'MNINonLinear',filesep,'Results'];
+
+ otherwise
+ error('ooopps..')
+
+end
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/BrainGraph/mdh/hb_make_single_connected.m b/BrainGraph/mdh/hb_make_single_connected.m
new file mode 100644
index 0000000..4974bf0
--- /dev/null
+++ b/BrainGraph/mdh/hb_make_single_connected.m
@@ -0,0 +1,76 @@
+function [outputFileName,sts,Nr] = hb_make_single_connected(inputFileName,conn,varargin)
+% Checks whether there is only one connected structure in the volume.
+% If the structure in the input volume is not single connected, a new
+% volume will be save in the same location as the original volume which
+% is connected.
+%
+%
+% Inputs:
+% inputFileName: Name of the source .nii file to be checked.
+% The name should be full, i.e. should include full
+% diectory address, file name and .nii at the end.
+%
+% conn: Desired connectivity. conn may have the following
+% scalar values:
+%
+% 4 two-dimensional four-connected neighborhood
+% 8 two-dimensional eight-connected neighborhood
+% 6 three-dimensional six-connected neighborhood
+% 18 three-dimensional 18-connected neighborhood
+% 26 three-dimensional 26-connected neighborhood
+%
+% varargin{1}: Name of the .nii file to be saved.
+% The name should be full, i.e. should include full
+% diectory address, file name and .nii at the end.
+%
+% Outputs:
+% outputFileName: Name of the saved .nii file. The name is full, i.e.
+% includes full diectory address, file name and .nii
+% at the end.
+%
+% sts: 1 if there were voxels removed and the volume saved
+% is diferent from the original volume, otherwise 0.
+%
+% Nr: Number of voxels removed.
+%
+% Hamid Behjat
+% Jan. 2017 / Feb 2018
+
+[p,n,e] = fileparts(inputFileName);
+control_params = {
+ 'outputFileName',fullfile(p,[num2str(conn),'connected_',n,e]),...
+ 'alwaysWriteConnected',0,'sS',1};
+
+argselectAssign(control_params);
+argselectCheck(control_params,varargin);
+argselectAssign(varargin);
+
+if sS, fprintf('\n* Making the structure a single-connected one..\n'); end
+
+vol = spm_read_vols(spm_vol(inputFileName));
+N_v0 = numel(find(vol(:)));
+
+% remove possible small isolated objects
+objSize = floor(numel(find(vol(:)))/3); % a large enough object size
+vol = bwareaopen(vol,objSize,conn);
+N_v1 = numel(find(vol(:)));
+Nr = N_v0-N_v1;
+
+if sS
+ fprintf([' # of voxels removed to make struture single-connected: ',...
+ num2str(Nr),'\n']);
+end
+
+if Nr~=0 || alwaysWriteConnected==1
+ sts=1;
+ V = spm_vol(inputFileName);
+ V.fname = outputFileName;
+ V.private.dat.fname = V.fname;
+ V.descrip = [V.descrip,' -- made connected (conn: ',num2str(conn),')'];
+ if sS, fprintf('- writing volume on disc: \n'); end
+ spm_write_vol(V,vol);
+ if sS, disp(V), end
+else
+ sts = 0;
+end
+if sS, fprintf('done.\n'); end
diff --git a/BrainGraph/mdh/mdh_adjacencymatrix_wb.m b/BrainGraph/mdh/mdh_adjacencymatrix_wb.m
new file mode 100644
index 0000000..4f5dbf9
--- /dev/null
+++ b/BrainGraph/mdh/mdh_adjacencymatrix_wb.m
@@ -0,0 +1,195 @@
+function [A,indices,G,O] = mdh_adjacencymatrix_wb(G,odfPow,neighb,dsi_root,param,varargin)
+% MDH_ADJACENCYMATRIX_WM Constructs white matter graph A matrix.
+%
+%
+%
+%
+% Contributers:
+% Martin Larsson
+% David Abramian
+% Hamid Behjat
+% 2017-2018.
+%
+% Code completely revised by: Anjali Tarun:
+% Changes made:
+% -- to use whole brain volume instead of WM only
+% -- incorporated structural boost using QA anisotropy (Yeh et. al)
+% -- removes outlier brain voxels to aid in convergence of the
+% eigendecomposition of the Laplacian
+
+% February 2018
+
+cntr = {
+ 'N_fibers',5,... % maximum count of the resolving fibers
+ };
+argselectAssign(cntr);
+argselectCheck(cntr,varargin);
+argselectAssign(varargin);
+
+fprintf('\n.Constructing WM graph adjacency matrix..\n')
+
+f = G.f;
+
+[wHeader,wVol] = ml_load_nifti(f.brainmask_fs_125);
+
+% Build ODFs.
+f = mdh_dsistudio_wb(f,dsi_root,'N_fibers',N_fibers);
+
+% Build A.
+[A,odf,index,odf_vertices,odf_faces,f,mask] = mdh_adjacencymatrix_wb_pvt(...
+ f.odfs,wVol,neighb,odfPow,f,param);
+
+wVol = mask;
+indices = find(wVol);
+G.N_wb = size(A,1);
+G.A_wb = A;
+G.A_wb_uniform = double(logical(A)); % non-weighted A
+G.indices_wb = indices;
+G.dim = wHeader.dim;
+[G.xyz(:,1), G.xyz(:,2), G.xyz(:,3)] = ind2sub(G.dim, indices);
+G.f = f;
+
+O = struct;
+O.odf_vertices = odf_vertices';
+O.odf_faces = odf_faces';
+O.index = index(:,indices);
+O.odf = odf;
+
+A_wb = A; %#ok
+indices_wb = indices; %#ok
+save(fullfile(G.f.graph,'A_wb.mat'),'A_wb');
+save(fullfile(G.f.graph,'indices_wb.mat'),'indices_wb');
+fprintf('.done.\n')
+end
+
+function [A,odf,index,odf_vertices,odf_faces,f,mask] = ...
+ mdh_adjacencymatrix_wb_pvt(f_odf,mask,ndim,odf_power,f,param)
+
+dim = size(mask);
+
+% Convert FIB to MAT.
+assert(strcmp(f_odf(end-6:end),'.fib.gz'),'Expected FIB.GZ-file for ODFs.');
+f_odf = f_odf(1:end-length('.fib.gz'));
+f_odf_fib = [f_odf '.fib.gz'];
+f_odf_mat = [f_odf '.mat'];
+
+if ~exist(f_odf_mat,'file')
+ ml_fibgz2mat(f_odf_fib);
+end
+
+% Load ODFs.
+[fa,odf,index,odf_vertices,odf_faces] = ml_dsi_load_odfs_wb(f_odf_mat,true);
+odf = odf.^odf_power;
+odf = [odf; odf];
+
+
+% Quantitative anisotropy, consider only major fiber
+Qa = max(fa,[],1);
+Qa = Qa./max(Qa);
+Qa(isnan(Qa)) = 0;
+
+% Error check, misc.
+if nnz(mask) ~= size(odf,2) || size(odf,1) ~= size(odf_vertices,2)
+ disp('Mismatching sizes.');
+ mask = reshape(Qa, dim);
+ mask = logical(mask);
+end
+if ~isequal(mask(:)',fa(1,:)~=0)
+ error('White matter masks misaligned.');
+end
+
+
+% check if there are voxels with no existing fiber (0 sum odf)
+zerovox = sum(odf,1);
+if(nnz(zerovox)~=size(odf,2))
+ disp(['removing voxels with no odf values.. There are ..',...
+ num2str(length(find(zerovox==0))),'.. of voxels to be removed'])
+ odf(:,zerovox==0)=[];
+ mask2 = zeros(dim);
+ mask2(mask) = zerovox;
+ mask = logical(mask2);
+end
+
+
+% Find neighboring voxels and direction of edges.
+[dirs,ndirs,~,npar] = ml_create_neighborhood(ndim);
+neigh = size(dirs,2);
+
+% Construct graph.
+%omega = 4*pi/neighborhood;
+omega = 4*pi/npar; % We still want 98 for 5x5x5 neighborhood.
+% 2*theta is a the apex angle of a cone.
+costheta = 1-omega/(2*pi);
+cone_set = ndirs'*odf_vertices > costheta;
+cone_set = cone_set./repmat(sum(cone_set,2),1,size(cone_set,2));
+
+% deltaS has no impact due to the normalization that follows.
+% deltaS = 4*pi/size(odf_vertices,2);
+% Pdiff = deltaS*cone_set*odf.*lni;
+
+
+maski = find(mask);
+voxels = length(maski);
+maski_inv = zeros(dim);
+maski_inv(maski) = 1:voxels;
+
+ci = repelem(maski,neigh,1);
+cs = zeros(size(ci,1),3);
+[cs(:,1),cs(:,2),cs(:,3)] = ind2sub(dim,ci);
+ns = cs+repmat(dirs',voxels,1);
+ni = sub2ind(dim,ns(:,1),ns(:,2),ns(:,3));
+ni = reshape(ni,neigh,voxels);
+ci = maski_inv(ci);
+ni = maski_inv(ni);
+lni = logical(ni);
+
+
+%Save copy of QA and updated brainmask
+h = spm_vol(f.white_125);
+hdr = cbiReadNiftiHeader(h.fname);
+niftifile = fullfile(f.extracted, 'QA_max.nii');
+% if ~exist(niftifile,'file')
+ cbiWriteNifti(niftifile,reshape(Qa,dim),hdr,'float32')
+% end
+niftifile2 = fullfile(f.extracted, 'new_brainmask.nii');
+% if ~exist(niftifile2,'file')
+ cbiWriteNifti(niftifile2,mask,hdr,'float32')
+% end
+
+Pdiff = cone_set*odf.*lni;
+Pdiff = 0.5*Pdiff./repmat(max(Pdiff),neigh,1);
+
+Qa = Qa(mask);
+Qa = repmat(Qa, [size(Pdiff,1), size(Qa,1)]);
+
+% Note: Some voxels return nan values upon normalizing the weights locally
+% (Anjali, 4 March 2018), putting them to 0.
+% For some subjects, Pdiff gives 0 weights in some voxels, odf appears 0 (Anjali,5 March 2018)..
+Pdiff(isnan(Pdiff)) = 0;
+
+% Multiply weights with Qa, gamma filtered with alpha, takes only the major fiber
+Pdiff = (Qa.^param.alpha).*Pdiff;
+
+A = sparse(ci(lni),ni(lni),Pdiff(lni),voxels,voxels);
+A = A+A';
+
+% Check if graph is connected
+conn = any(A);
+if length(conn)==size(A,1)
+ disp('Fully connected graph..')
+else
+ disp(['There are ..',num2str(length(conn)),'.. unconnected voxels'])
+end
+
+
+% Check for outlier voxels
+% A2= A(find(A~=0));
+% mean_allvoxel = mean(A2(:));
+% small_num = eps; % order of magnitude
+% thresh = mean_allvoxel*small_num;
+% out = find(A2
+indices_wb = indices; %#ok
+save(fullfile(G.f.graph,'A_wb.mat'),'A_wb');
+save(fullfile(G.f.graph,'indices_wb.mat'),'indices_wb');
+fprintf('.done.\n')
+end
+
+function [A,odf,index,odf_vertices,odf_faces,f,mask] = ...
+ mdh_adjacencymatrix_wb_pvt(f_odf,mask,ndim,odf_power,f,param)
+
+dim = size(mask);
+
+% Convert FIB to MAT.
+assert(strcmp(f_odf(end-6:end),'.fib.gz'),'Expected FIB.GZ-file for ODFs.');
+f_odf = f_odf(1:end-length('.fib.gz'));
+f_odf_fib = [f_odf '.fib.gz'];
+f_odf_mat = [f_odf '.mat'];
+
+if ~exist(f_odf_mat,'file')
+ ml_fibgz2mat(f_odf_fib);
+end
+
+% Load ODFs.
+[fa,odf,index,odf_vertices,odf_faces] = ml_dsi_load_odfs_wb(f_odf_mat,true);
+odf = odf.^odf_power;
+odf = [odf; odf];
+
+
+% Quantitative anisotropy, consider only major fiber
+Qa = max(fa,[],1);
+Qa = Qa./max(Qa);
+Qa(isnan(Qa)) = 0;
+
+% Error check, misc.
+if nnz(mask) ~= size(odf,2) || size(odf,1) ~= size(odf_vertices,2)
+ disp('Mismatching sizes.');
+ mask = reshape(Qa, dim);
+ mask = logical(mask);
+end
+if ~isequal(mask(:)',fa(1,:)~=0)
+ error('White matter masks misaligned.');
+end
+
+
+% check if there are voxels with no existing fiber (0 sum odf)
+zerovox = sum(odf,1);
+if(nnz(zerovox)~=size(odf,2))
+ disp(['removing voxels with no odf values.. There are ..',...
+ num2str(length(find(zerovox==0))),'.. of voxels to be removed'])
+ odf(:,zerovox==0)=[];
+ mask2 = zeros(dim);
+ mask2(mask) = zerovox;
+ mask = logical(mask2);
+end
+
+
+% Find neighboring voxels and direction of edges.
+[dirs,ndirs,~,npar] = ml_create_neighborhood(ndim);
+neigh = size(dirs,2);
+
+% Construct graph.
+%omega = 4*pi/neighborhood;
+omega = 4*pi/npar; % We still want 98 for 5x5x5 neighborhood.
+% 2*theta is a the apex angle of a cone.
+costheta = 1-omega/(2*pi);
+cone_set = ndirs'*odf_vertices > costheta;
+cone_set = cone_set./repmat(sum(cone_set,2),1,size(cone_set,2));
+
+% deltaS has no impact due to the normalization that follows.
+% deltaS = 4*pi/size(odf_vertices,2);
+% Pdiff = deltaS*cone_set*odf.*lni;
+
+
+maski = find(mask);
+voxels = length(maski);
+maski_inv = zeros(dim);
+maski_inv(maski) = 1:voxels;
+
+ci = repelem(maski,neigh,1);
+cs = zeros(size(ci,1),3);
+[cs(:,1),cs(:,2),cs(:,3)] = ind2sub(dim,ci);
+ns = cs+repmat(dirs',voxels,1);
+ni = sub2ind(dim,ns(:,1),ns(:,2),ns(:,3));
+ni = reshape(ni,neigh,voxels);
+ci = maski_inv(ci);
+ni = maski_inv(ni);
+lni = logical(ni);
+
+
+%Save copy of QA and updated brainmask
+h = spm_vol(f.white_125);
+hdr = cbiReadNiftiHeader(h.fname);
+niftifile = fullfile(f.extracted, 'QA_max.nii');
+% if ~exist(niftifile,'file')
+ cbiWriteNifti(niftifile,reshape(Qa,dim),hdr,'float32')
+% end
+niftifile2 = fullfile(f.extracted, 'new_brainmask.nii');
+% if ~exist(niftifile2,'file')
+ cbiWriteNifti(niftifile2,mask,hdr,'float32')
+% end
+
+Pdiff = cone_set*odf.*lni;
+Pdiff = 0.5*Pdiff./repmat(max(Pdiff),neigh,1);
+
+Qa = Qa(mask);
+Qa = repmat(Qa, [size(Pdiff,1), size(Qa,1)]);
+
+% Note: Some voxels return nan values upon normalizing the weights locally
+% (Anjali, 4 March 2018), putting them to 0.
+% For some subjects, Pdiff gives 0 weights in some voxels, odf appears 0 (Anjali,5 March 2018)..
+Pdiff(isnan(Pdiff)) = 0;
+
+% Multiply weights with Qa, gamma filtered with alpha, takes only the major fiber
+Pdiff = (Qa.^param.alpha).*Pdiff;
+
+A = sparse(ci(lni),ni(lni),Pdiff(lni),voxels,voxels);
+A = A+A';
+
+% Check if graph is connected
+conn = any(A);
+if length(conn)==size(A,1)
+ disp('Fully connected graph..')
+else
+ disp(['There are ..',num2str(length(conn)),'.. unconnected voxels'])
+end
+
+
+% Check for outlier voxels
+% A2= A(find(A~=0));
+% mean_allvoxel = mean(A2(:));
+% small_num = eps; % order of magnitude
+% thresh = mean_allvoxel*small_num;
+% out = find(A2
+indices_wm = indices; %#ok
+save(fullfile(G.f.graph,'A_wm.mat'),'A_wm');
+save(fullfile(G.f.graph,'indices_wm.mat'),'indices_wm');
+fprintf('.done.\n')
+end
+
+function [A,odf,index,odf_vertices,odf_faces,f] = ...
+ mdh_adjacencymatrix_wm_pvt(f_odf,mask,ndim,odf_power,f)
+
+dim = size(mask);
+
+% Convert FIB to MAT.
+assert(strcmp(f_odf(end-6:end),'.fib.gz'),'Expected FIB.GZ-file for ODFs.');
+f_odf = f_odf(1:end-length('.fib.gz'));
+f_odf_fib = [f_odf '.fib.gz'];
+f_odf_mat = [f_odf '.mat'];
+
+if ~exist(f_odf_mat,'file')
+ ml_fibgz2mat(f_odf_fib);
+end
+
+% Load ODFs.
+[fa,odf,index,odf_vertices,odf_faces] = ml_dsi_load_odfs(f_odf_mat,true);
+odf = odf.^odf_power;
+odf = [odf; odf];
+
+% Error check, misc.
+if nnz(mask) ~= size(odf,2) || size(odf,1) ~= size(odf_vertices,2)
+ error('Mismatching sizes.');
+end
+if ~isequal(mask(:)',fa(1,:)~=0)
+ error('White matter masks misaligned.');
+end
+
+maski = find(mask);
+voxels = length(maski);
+maski_inv = zeros(dim);
+maski_inv(maski) = 1:voxels;
+
+% Find neighboring voxels and direction of edges.
+[dirs,ndirs,~,npar] = ml_create_neighborhood(ndim);
+neigh = size(dirs,2);
+
+ci = repelem(maski,neigh,1);
+cs = zeros(size(ci,1),3);
+[cs(:,1),cs(:,2),cs(:,3)] = ind2sub(dim,ci);
+ns = cs+repmat(dirs',voxels,1);
+ni = sub2ind(dim,ns(:,1),ns(:,2),ns(:,3));
+ni = reshape(ni,neigh,voxels);
+ci = maski_inv(ci);
+ni = maski_inv(ni);
+lni = logical(ni);
+
+% Construct graph.
+%omega = 4*pi/neighborhood;
+omega = 4*pi/npar; % We still want 98 for 5x5x5 neighborhood.
+% 2*theta is a the apex angle of a cone.
+costheta = 1-omega/(2*pi);
+cone_set = ndirs'*odf_vertices > costheta;
+cone_set = cone_set./repmat(sum(cone_set,2),1,size(cone_set,2));
+
+% deltaS has no impact due to the normalization that follows.
+% deltaS = 4*pi/size(odf_vertices,2);
+% Pdiff = deltaS*cone_set*odf.*lni;
+
+Pdiff = cone_set*odf.*lni;
+Pdiff = 0.5*Pdiff./repmat(max(Pdiff),neigh,1);
+
+A = sparse(ci(lni),ni(lni),Pdiff(lni),voxels,voxels);
+A = A+A';
+
+f.odf_fib = f_odf_fib;
+f.odf_mat = f_odf_mat;
+end
diff --git a/BrainGraph/mdh/mdh_dsistudio.m b/BrainGraph/mdh/mdh_dsistudio.m
new file mode 100755
index 0000000..096601a
--- /dev/null
+++ b/BrainGraph/mdh/mdh_dsistudio.m
@@ -0,0 +1,66 @@
+function [f,logfile1,logfile2] = mdh_dsistudio(f,dsi_root,varargin)
+% MDH_DSISTUDIO Builds ODFs by running DSI STUDIO.
+
+% HB [Feb 2018]
+% - specifying dsi studio root
+% - option to vary N_fibers
+% - input/output file struture; more robust file namings.
+%
+% Contributers:
+% Martin Larsson
+% David Abramian
+% Hamid Behjat
+% 2017-2018.
+
+cntr = {
+ 'N_fibers',5,... % maximum count of the resolving fibers
+ };
+argselectAssign(cntr);
+argselectCheck(cntr,varargin);
+argselectAssign(varargin);
+
+f.source = fullfile(f.diffusion,'data.nii.gz');
+f.out_src = fullfile(f.diffusion,'odfs.src.gz');
+f.mask = fullfile(f.extracted,'white_1.25.nii');
+f.odfs = fullfile(f.diffusion,'odfs.fib.gz');
+
+d = fullfile(fileparts(f.out_src),['odfs.src.gz.odf8.f',...
+ num2str(N_fibers),'rec.bal.012fy.rdi.gqi.1.25.fib.gz']);
+
+if ~exist(f.out_src,'file')
+ fprintf('..Running DSI Studio: obtaining .src file \n')
+ command = [dsi_root,'/dsi_studio --action=src --source=',...
+ f.source,' --output=',f.out_src];
+ [sts,logfile1] = system(command);
+ if ~sts
+ error('error running dsi studio [constructing odfs]');
+ end
+ fprintf('..done.\n');
+
+ fprintf('..Running DSI Studio: obtaining .fib file \n')
+ command = [dsi_root,filesep,'dsi_studio --action=rec --source=',...
+ f.out_src,' --thread=2 --mask=',f.mask,...
+ ' --method=4 --param0=1.25 --odf_order=8 --num_fiber=',...
+ num2str(N_fibers),' --check_btable=1 --record_odf=1'];
+ [sts,logfile2] = system(command);
+%
+ if sts
+% out_string = 'output odfs';
+% I = strfind(logfile2,out_string) + length(out_string) + 1;
+% f.out_fib = logfile2(I:end-1);
+% copyfile(f.out_fib,f.odfs);
+% outfile = fullfile(f.diffusion,'odfs.src.gz.odf*');
+% copyfile(fullfile(f.diffusion,outfile(1).name),f.odfs);
+
+ outfile = dir(fullfile(f.diffusion,'odf.src.gz.odf*'));
+ copyfile(fullfile(f.diffusion,outfile(1).name),f.odfs);
+
+
+ else
+ error('error running dsi studio [constructing odfs]');
+ end
+% fprintf('..done.\n');
+else
+ fprintf('..ODF file already exists.\n')
+end
+
diff --git a/BrainGraph/mdh/mdh_dsistudio_wb.m b/BrainGraph/mdh/mdh_dsistudio_wb.m
new file mode 100644
index 0000000..518ab62
--- /dev/null
+++ b/BrainGraph/mdh/mdh_dsistudio_wb.m
@@ -0,0 +1,67 @@
+function [f,logfile1,logfile2] = mdh_dsistudio_wb(f,dsi_root,varargin)
+% MDH_DSISTUDIO Builds ODFs by running DSI STUDIO.
+
+% HB [Feb 2018]
+% - specifying dsi studio root
+% - option to vary N_fibers
+% - input/output file struture; more robust file namings.
+%
+% Contributers:
+% Martin Larsson
+% David Abramian
+% Hamid Behjat
+% 2017-2018.
+
+% Code revised by: Anjali Tarun:
+% Changes made:
+% -- to use whole brain volume instead of WM only
+% -- obtains whole brain brain mask
+% February 2018
+
+
+cntr = {
+ 'N_fibers',5,... % maximum count of the resolving fibers
+ };
+argselectAssign(cntr);
+argselectCheck(cntr,varargin);
+argselectAssign(varargin);
+
+f.source = fullfile(f.diffusion,'data.nii.gz');
+f.out_src = fullfile(f.diffusion,'odfs.src.gz');
+f.mask = fullfile(f.T1w,'brainmask_fs_125.nii');
+f.odfs = fullfile(f.diffusion,'odfs.fib.gz');
+
+d = fullfile(fileparts(f.out_src),['odfs.src.gz.odf8.f',...
+ num2str(N_fibers),'rec.bal.012fy.rdi.gqi.1.25.fib.gz']);
+
+if ~exist(f.out_src,'file')
+ fprintf('..Running DSI Studio: obtaining .src file \n')
+ command = [dsi_root,'/dsi_studio --action=src --source=',...
+ f.source,' --output=',f.out_src];
+ [sts,logfile1] = system(command);
+ if ~sts
+ error('error running dsi studio [constructing odfs]');
+ end
+ fprintf('..done.\n');
+
+ fprintf('..Running DSI Studio: obtaining .fib file \n')
+ command = [dsi_root,filesep,'dsi_studio --action=rec --source=',...
+ f.out_src,' --thread=2 --mask=',f.mask,...
+ ' --method=4 --param0=1.25 --odf_order=8 --num_fiber=',...
+ num2str(N_fibers),' --check_btable=1 --record_odf=1'];
+ [sts,logfile2] = system(command);
+ if sts
+% out_string = 'output odfs';
+% I = strfind(logfile2,out_string) + length(out_string) + 1;
+% f.out_fib = logfile2(I:end-1);
+% copyfile(f.out_fib,f.odfs);
+ outfile = dir(fullfile(f.diffusion,'odfs.src.gz.odf*'));
+ copyfile(fullfile(f.diffusion,outfile(1).name),f.odfs);
+ else
+ error('error running dsi studio [constructing odfs]');
+ end
+ fprintf('..done.\n');
+else
+ fprintf('..ODF file already exists.\n')
+end
+
diff --git a/BrainGraph/mdh/mdh_preprocess.m b/BrainGraph/mdh/mdh_preprocess.m
new file mode 100644
index 0000000..bac7488
--- /dev/null
+++ b/BrainGraph/mdh/mdh_preprocess.m
@@ -0,0 +1,106 @@
+
+
+
+
+
+function G = mdh_preprocess(subjID,hcp_root) % Based on sgwt_preprocess.m
+% MDH_PREPROCESS Creates necessary files for building GM, WM,GM+WM graphs.
+%
+%
+%
+% Contributers:
+% Martin Larsson
+% David Abramian
+% Hamid Behjat
+% 2017-2018.
+
+
+% Code revised by: Anjali Tarun:
+% Changes made:
+% -- to use whole brain volume instead of WM only
+% -- obtains whole brain brain mask
+% February 2018
+
+f = struct;
+
+% Folders.
+f.T1w = fullfile(hcp_root,subjID,'T1w');
+f.diffusion = fullfile(f.T1w,'Diffusion');
+f.extracted = fullfile(f.T1w,'extracted');
+f.graph = fullfile(f.T1w,'graph');
+if ~exist(f.extracted,'dir')
+ mkdir(f.extracted);
+end
+if ~exist(f.graph,'dir')
+ mkdir(f.graph);
+end
+
+% Files.
+f.ref_125 = fullfile(f.diffusion,'nodif_brain_mask.nii');
+f.aparc_aseg_07 = fullfile(f.T1w,'aparc+aseg.nii');
+f.aparc_aseg_125 = fullfile(f.T1w,'aparc+aseg_1.25.nii');
+f.ribbon_125 = fullfile(f.extracted,'ribbon_1.25.nii');
+f.white_125 = fullfile(f.extracted,'white_1.25.nii');
+f.brainmask_fs = fullfile(f.T1w, 'brainmask_fs.nii');
+f.brainmask_fs_125 = fullfile(f.T1w, 'brainmask_fs_125.nii');
+f.corpus_callosum_125 = fullfile(f.extracted,'corpus_callosum_1.25.nii');
+f.ribbon_125_conn = fullfile(f.extracted,'ribbon_1.25_connected.nii');
+f.gmwm_125 = fullfile(f.extracted,'gmwm_125.nii');
+f.T1_125 = fullfile(f.T1w, 'T1w_acpc_dc_restore_1.25.nii');
+
+% extracts
+% Downsample aparc+aseg.
+ml_extract(f.aparc_aseg_07);
+ml_extract(f.ref_125);
+if ~exist(f.aparc_aseg_125,'file')
+ f.aparc_aseg_125 = da_resize(f.aparc_aseg_07,f.ref_125);
+end
+
+if ~exist(f.aparc_aseg_125,'file')
+ f.aparc_aseg_125 = da_resize(f.aparc_aseg_07,f.ref_125);
+end
+
+if ~exist(f.brainmask_fs_125,'file')
+ ml_extract(f.brainmask_fs);
+
+ f.brainmask_fs_125 = da_resize(f.brainmask_fs,f.ref_125);
+ % Clean mask.
+ [wh,wv] = ml_load_nifti(f.brainmask_fs_125);
+ cc = bwconncomp(wv,26);
+ [~,ind] = max(cellfun(@length,cc.PixelIdxList));
+ wv = zeros(size(wv));
+ wv(cc.PixelIdxList{ind}) = 1;
+ spm_write_vol(wh,wv);
+end
+
+
+% Extract volumes.
+ids_ribbon = [3 42 1000:2999];
+ids_corpus_callosum = [86 251:255];
+ids_white = [2 41 ids_corpus_callosum];
+
+if ~exist(f.corpus_callosum_125,'file')
+ ml_extract_volumes(f.aparc_aseg_125,f.corpus_callosum_125,...
+ ids_corpus_callosum);
+end
+
+if ~exist(f.ribbon_125,'file')
+ ml_extract_volumes(f.aparc_aseg_125,f.ribbon_125,ids_ribbon);
+end
+
+if ~exist(f.white_125,'file')
+ ml_extract_volumes(f.aparc_aseg_125,f.white_125,ids_white);
+
+ % Clean white matter.
+ [wh,wv] = ml_load_nifti(f.white_125);
+ cc = bwconncomp(wv,26);
+ [~,ind] = max(cellfun(@length,cc.PixelIdxList));
+ wv = zeros(size(wv));
+ wv(cc.PixelIdxList{ind}) = 1;
+ spm_write_vol(wh,wv);
+end
+
+
+
+G = struct;
+G.f = f;
diff --git a/BrainGraph/mdh/ml_catvars.m b/BrainGraph/mdh/ml_catvars.m
new file mode 100755
index 0000000..78779d4
--- /dev/null
+++ b/BrainGraph/mdh/ml_catvars.m
@@ -0,0 +1,25 @@
+function ml_catvars(var,dim)
+% ML_CATVARS Concatinates variables in the calling workspace.
+% ML_CATVARS(var,dim)
+% Will concatinate the variables var0,var1,...,varn along the
+% dimension dim in the calling workspace. var0,var1,...,varn will be
+% cleared from the workspace and replaced with a single variable var.
+%
+% See also ML_SPLITVAR.
+%
+% Author:
+% Martin Larsson
+% March 2017
+
+ i = 0;
+ value = [];
+ name = [var '0'];
+ while evalin('caller',['exist(''' name ''',''var'')'])
+ value = cat(dim,value,evalin('caller',name));
+ evalin('caller',['clear(''' name ''')']);
+ i = i+1;
+ name = [var num2str(i)];
+ end
+ assignin('caller',var,value);
+end
+
diff --git a/BrainGraph/mdh/ml_create_neighborhood.m b/BrainGraph/mdh/ml_create_neighborhood.m
new file mode 100755
index 0000000..05ea095
--- /dev/null
+++ b/BrainGraph/mdh/ml_create_neighborhood.m
@@ -0,0 +1,43 @@
+function [dirs,ndirs,norms,npar] = ml_create_neighborhood(dim,rempar)
+% ML_CREATE_NEIGHBORHOOD Creates direction vectors for a 3D neighborhood.
+% [dirs,ndirs,norms] = ML_CREATE_NEIGHBORHOOD(dim)
+% dim specifies the size of the neighborhood cube, eg. dim=3 would
+% define a 26 neighborhood.
+% [dirs,ndirs,norms] = ML_CREATE_NEIGHBORHOOD(dim,rempar)
+% rempar specifies whether or not to remove parallel vectors, eg.
+% ML_CREATE_NEIGHBORHOOD(5,true) would return 98 and not 124 vectors.
+%
+% Author:
+% Martin Larsson
+% June 2017
+
+ r = (dim-1)/2;
+ [y,x,z] = meshgrid(-r:r,-r:r,-r:r);
+ not_origin = true(1,numel(x));
+ not_origin((numel(x)+1)/2) = false;
+ dirs = [x(not_origin); y(not_origin); z(not_origin)];
+
+ if nargin > 1 && rempar || nargout > 1
+ % Normalized vectors.
+ norms = sqrt(sum(dirs.^2));
+ ndirs = dirs./repmat(norms,size(dirs,1),1);
+ end
+
+ if nargout > 3 || nargin > 1 && rempar
+ % Find parallel vectors.
+ same = ndirs'*ndirs > 1-1e-6;
+ include = false(size(dirs,2),1);
+ for i=1:size(dirs,2)
+ include(i) = norms(i) <= min(norms(same(i,:)));
+ end
+ npar = sum(include);
+ end
+
+ if nargin > 1 && rempar
+ % Remove parallel vectors.
+ dirs = dirs(:,include);
+ ndirs = ndirs(:,include);
+ norms = norms(include);
+ end
+end
+
diff --git a/BrainGraph/mdh/ml_dsi_load_odfs.m b/BrainGraph/mdh/ml_dsi_load_odfs.m
new file mode 100755
index 0000000..273fc54
--- /dev/null
+++ b/BrainGraph/mdh/ml_dsi_load_odfs.m
@@ -0,0 +1,51 @@
+function [fa,odf,index,odf_vertices,odf_faces] = ml_dsi_load_odfs(file_name,flipY)
+% ML_DSI_LOAD_ODFS Load ODF data outputted from DSI Studio.
+% [fa,odf,index,odf_vertices,odf_faces] = ML_DSI_LOAD_ODFS(file_name)
+% file_name specifies the MAT-file that contains the ODF data. DSI
+% Studio outputs a FIB.GZ-file that can be converted to a MAT-file
+% using ml_fibgz2mat.
+% [fa,odf,index,odf_vertices,odf_faces] = ML_DSI_LOAD_ODFS(file_name,flipY,dim)
+% flipY specifies whether or not to flip the data along the Y-axis
+% (second dimention). This is required to match the voxel space of
+% DSI Studio (-X,-Y,Z) with the HCP dataset (-X,Y,Z).
+%
+% See also ML_FIBGZ2MAT.
+%
+% Author:
+% Martin Larsson
+% May 2017
+
+ if ~exist('file_name', 'var') || isempty(file_name)
+ [file_name, path_name] = uigetfile({'*.mat','MATLAB files'});
+ if file_name == 0
+ return;
+ end
+ file_name = fullfile(path_name, file_name);
+ end
+
+ load(file_name);
+
+ ml_catvars('fa',1);
+ ml_catvars('odf',2);
+ ml_catvars('index',1);
+ index = index+1;
+
+ if nargin > 1 && flipY
+ fai = find(fa(1,:)~=0);
+ fai_inv = zeros(dimension);
+ fai_inv(fai) = 1:length(fai);
+ fai_inv = flip(fai_inv,2);
+ perm = fai_inv(fai_inv~=0);
+ odf = odf(:,perm);
+ odf_vertices(2,:) = -odf_vertices(2,:);
+
+ for i=1:size(fa,1)
+ fa(i,:) = reshape(flip(reshape(fa(i,:),dimension),2),1,size(fa,2));
+ index(i,:) = reshape(flip(reshape(index(i,:),dimension),2),1,size(index,2));
+ end
+ end
+
+
+
+end
+
diff --git a/BrainGraph/mdh/ml_dsi_load_odfs_wb.m b/BrainGraph/mdh/ml_dsi_load_odfs_wb.m
new file mode 100644
index 0000000..7c622bb
--- /dev/null
+++ b/BrainGraph/mdh/ml_dsi_load_odfs_wb.m
@@ -0,0 +1,50 @@
+function [fa,odf,index,odf_vertices,odf_faces] = ml_dsi_load_odfs_wb(file_name,flipY)
+% ML_DSI_LOAD_ODFS Load ODF data outputted from DSI Studio.
+% [fa,odf,index,odf_vertices,odf_faces] = ML_DSI_LOAD_ODFS(file_name)
+% file_name specifies the MAT-file that contains the ODF data. DSI
+% Studio outputs a FIB.GZ-file that can be converted to a MAT-file
+% using ml_fibgz2mat.
+% [fa,odf,index,odf_vertices,odf_faces] = ML_DSI_LOAD_ODFS(file_name,flipY,dim)
+% flipY specifies whether or not to flip the data along the Y-axis
+% (second dimention). This is required to match the voxel space of
+% DSI Studio (-X,-Y,Z) with the HCP dataset (-X,Y,Z).
+%
+% See also ML_FIBGZ2MAT.
+%
+% Author:
+% Martin Larsson
+% May 2017
+
+ if ~exist('file_name', 'var') || isempty(file_name)
+ [file_name, path_name] = uigetfile({'*.mat','MATLAB files'});
+ if file_name == 0
+ return;
+ end
+ file_name = fullfile(path_name, file_name);
+ end
+
+ load(file_name);
+
+ ml_catvars('fa',1);
+ ml_catvars('odf',2);
+ ml_catvars('index',1);
+ index = index+1;
+
+ if nargin > 1 && flipY
+ fai = find(fa(1,:)~=0);
+ fai_inv = zeros(dimension);
+ fai_inv(fai) = 1:length(fai);
+ fai_inv = flip(fai_inv,2);
+ perm = fai_inv(fai_inv~=0);
+ odf = odf(:,perm);
+ odf_vertices(2,:) = -odf_vertices(2,:);
+
+ for i=1:size(fa,1)
+ fa(i,:) = reshape(flip(reshape(fa(i,:),dimension),2),1,size(fa,2));
+ index(i,:) = reshape(flip(reshape(index(i,:),dimension),2),1,size(index,2));
+ end
+ end
+
+
+end
+
diff --git a/BrainGraph/mdh/ml_extract.m b/BrainGraph/mdh/ml_extract.m
new file mode 100755
index 0000000..8162a32
--- /dev/null
+++ b/BrainGraph/mdh/ml_extract.m
@@ -0,0 +1,19 @@
+function ml_extract(fname)
+% ML_EXTRACT Extracts a file.
+% ML_EXTRACT(fname)
+% If the given files does not exists this function will look for
+% [fname '.gz'] and extract it.
+%
+% Author:
+% Martin Larsson
+% June 2017
+
+ if ~exist(fname,'file')
+ comp_fname = [fname '.gz'];
+ if ~exist(comp_fname,'file')
+ error(['File does not exist: ' fname]);
+ end
+ gunzip(comp_fname);
+ end
+end
+
diff --git a/BrainGraph/mdh/ml_extract_volumes.m b/BrainGraph/mdh/ml_extract_volumes.m
new file mode 100755
index 0000000..1ab32ff
--- /dev/null
+++ b/BrainGraph/mdh/ml_extract_volumes.m
@@ -0,0 +1,31 @@
+function ml_extract_volumes(seg,varargin)
+% ML_EXTRACT_VOLUMES Extracts NIFTI-volumes from a segmentation.
+% ML_EXTRACT_VOLUMES(seg,...)
+% seg is the file name of a segmentation NIFTI-file. Additional
+% arguments specify pairs of file names and ID-vectors. A NIFTI-file
+% will be saved to the specified location containing the voxels in
+% seg matching the IDs provided in the ID-vector.
+%
+% Examples:
+% ML_EXTRACT_VOLUMES('aparc+aseg.nii','ribbon.nii',1000:2999);
+% ML_EXTRACT_VOLUMES('aparc+aseg.nii','thalamus.nii',[9 10 48 49]);
+%
+% Author:
+% Martin Larsson
+% June 2017
+
+ if nargin <= 1
+ return;
+ end
+
+ [segh,segv] = ml_load_nifti(seg);
+
+ h = segh;
+ h.dt(1) = 2; % Data type: unsigned char (the smallest) 12/2017 D
+ for i=1:2:length(varargin)
+ h.fname = varargin{i};
+ v = ismember(segv,varargin{i+1});
+ spm_write_vol(h,v);
+ end
+end
+
diff --git a/BrainGraph/mdh/ml_fibgz2mat.m b/BrainGraph/mdh/ml_fibgz2mat.m
new file mode 100755
index 0000000..939f113
--- /dev/null
+++ b/BrainGraph/mdh/ml_fibgz2mat.m
@@ -0,0 +1,27 @@
+function [mat_file_name] = ml_fibgz2mat(file_name)
+% ML_FIBGZ2MAT Converts FIB.GZ-files into MAT-files.
+% ML_FIBGZ2MAT(file_name)
+% file_name specifies the FIB.GZ-file to convert. If no file name is
+% provided a file dialog will be opened. The name of the new file is
+% returned.
+%
+% Author:
+% Martin Larsson
+% May 2017
+
+ if ~exist('file_name', 'var') || isempty(file_name)
+ [file_name, path_name] = uigetfile({'*.fib.gz', 'Compressed FIB-files'});
+ if file_name == 0
+ return;
+ end
+ file_name = fullfile(path_name, file_name);
+ end
+
+ [pathstr,name,~] = fileparts(file_name);
+ fib_file_name = fullfile(pathstr, name);
+ [pathstr,name,~] = fileparts(fib_file_name);
+ mat_file_name = fullfile(pathstr, strcat(name, '.mat'));
+ gunzip(file_name);
+ movefile(fib_file_name, mat_file_name);
+end
+
diff --git a/BrainGraph/mdh/ml_load_nifti.m b/BrainGraph/mdh/ml_load_nifti.m
new file mode 100755
index 0000000..b070db8
--- /dev/null
+++ b/BrainGraph/mdh/ml_load_nifti.m
@@ -0,0 +1,42 @@
+function [h,v] = ml_load_nifti(fname,volumes)
+% ML_LOAD_NIFTI Extracts and loads a NIFTI-file.
+% ML_LOAD_NIFTI(fname)
+% Makes sure that the file fname exists. If not the function looks
+% for [fname '.gz'] and extracts it.
+% h = ML_LOAD_NIFTI(fname)
+% Loads the header of the NIFTI-file specified by fname.
+% [h,v] = ML_LOAD_NIFTI(fname)
+% Loads the header and volume data of the NIFTI-file specified by
+% fname.
+% ... = ML_LOAD_NIFTI(fname,volumes)
+% If the NIFTI-file is a 4D volume only the volumes specified in
+% volumes are loaded.
+%
+% Author:
+% Martin Larsson
+% June 2017
+
+% Improvements in efficiency 12/2017 D
+
+% ml_extract(fname);
+
+ if nargout > 0
+
+ if nargin > 1
+ if ~iscolumn(volumes)
+ volumes = volumes';
+ end
+ volumes = num2str(volumes);
+ files = strcat(fname,',',volumes);
+ else
+ files = fname;
+ end
+
+ h = spm_vol(files);
+% h = cell2mat(h);
+
+ if nargout > 1
+ v = spm_read_vols(h);
+ end
+ end
+end
diff --git a/BrainGraph/slepEigsLaplacian.m b/BrainGraph/slepEigsLaplacian.m
new file mode 100644
index 0000000..74106e6
--- /dev/null
+++ b/BrainGraph/slepEigsLaplacian.m
@@ -0,0 +1,14 @@
+function [Utr,S1]=slepEigsLaplacian(A,CONST_W,opts)
+
+D = speye(size(A,1));
+
+[Utr,S1]=eigs(@Binline,size(A,1),max(CONST_W),'sa',opts);
+% [Utr,S1]=eigs(@Binline,size(A,1),max(CONST_W),'sr',opts);
+
+
+
+function y=Binline(x)
+ y=D*x-A*x;
+end
+
+end
diff --git a/BrainGraph/slepNormalize.m b/BrainGraph/slepNormalize.m
new file mode 100644
index 0000000..30ecef9
--- /dev/null
+++ b/BrainGraph/slepNormalize.m
@@ -0,0 +1,37 @@
+% Construct normalized (sparse) adjacency matrix
+%
+% returns
+% A - normalized adjacency matrix
+% D - degree vector
+%
+% Dimitri Van De Ville, Sep 2014
+
+% Modified by Anjali Tarun, include random walk normalization, April 3 2016
+
+function [A,D]=slepNormalize(A,CONST_NORMALIZE,CONST_NORMALIZE_type)
+
+if ~issparse(A)
+ warning('Adjacency matrix A is not sparse.');
+end
+
+msize=size(A,1);
+
+% normalize adjacency matrix
+if CONST_NORMALIZE
+ D=sum(A);
+
+ switch CONST_NORMALIZE_type
+ case 1
+ Dn=1./sqrt(D);
+ Dn=spdiags(Dn.',0,spalloc(msize,msize,msize));
+ A=Dn*A*Dn;
+ case 2
+ Dn=1./D;
+ Dn=spdiags(Dn.',0,spalloc(msize,msize,msize));
+ A=Dn*A;
+ end
+end
+
+% construct D
+D=sum(A);
+D=spdiags(D.',0,spalloc(msize,msize,msize));
diff --git a/GraphSignalRecovery/CheckParameterLambda.m b/GraphSignalRecovery/CheckParameterLambda.m
new file mode 100644
index 0000000..2ffd724
--- /dev/null
+++ b/GraphSignalRecovery/CheckParameterLambda.m
@@ -0,0 +1,43 @@
+
+function lambda = CheckParameterLambda(param,BB,L,BC)
+
+ lambdacheck = fullfile(param.EigPath, ['RegularizationParameters_LCurve_',param.subject,'_',num2str(param.smoothR),'.mat']);
+
+ if ~exist(lambdacheck, 'file')
+
+ % Load preliminary fMRI volume
+ Inputs_fMRI
+
+ i = 20; % frame to evaluate
+ vol = V(i,:);
+
+ clear V
+
+ disp('Calculating correct parameter lambda..')
+ lambdas = 1:60;
+% diffs = zeros(length(lambdas),1);
+ x1norms = zeros(length(lambdas),1);
+ parfor c = 1:length(lambdas)
+ lambda = lambdas(c);
+ L_mat = BB+lambda*L;
+ L1 = ichol(L_mat);
+ xsignal = zeros(size(L1,1),1);
+ xsignal(BC) = vol(BC);
+ [x1,~,~,~,~] = pcg(L_mat,xsignal,1e-8,100,L1,L1');
+ x1norms(c) = norm(x1'*(L*x1));
+% diff = B*x1-xsignal;
+% diffs(c) = norm(diff)^2;
+ end
+
+ lambdaStars.lambdas = lambdas;
+ lambdaStars.x1norms = x1norms;
+ save(fullfile(param.EigPath, ['RegularizationParameters_LCurve_',param.subject,'.mat']),'lambdaStars')
+
+ else
+ load(lambdacheck)
+ x1norms = lambdaStars.x1norms;
+ lambdas = lambdaStars.lambdas;
+ end
+
+ lambda = knee_pt(x1norms,lambdas);
+end
\ No newline at end of file
diff --git a/Inputs/Inputs_BrainGrid.m b/Inputs/Inputs_BrainGrid.m
new file mode 100644
index 0000000..08a22b2
--- /dev/null
+++ b/Inputs/Inputs_BrainGrid.m
@@ -0,0 +1,35 @@
+
+% Contains the input parameters for the construction of the brain grid
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+% Specifies the parameters of the braingraph and the eigendecomposition
+param.odfPow = 40; % ODF power.
+param.neighborhood = 3; % 3 or 5 [This is for white matter; gray matter always 3]
+param.N_fibers = 5; % dsi studio's and M&L's default: 5
+% Tune FA/QA
+param.alpha = 2;
+
+% Choose the parameters for the eigendecomposition of the Laplacian to the
+% obtain eigenmodes..
+% whether to have a percent bandwidth or constant
+param.percent = 0;
+param.bandwidth = 1e-3;
+param.c_bandwidth = 5; % just to visualize first 5 eigenmodes
+param.opts.issym=1;
+param.opts.isreal=1;
+param.opts.maxit=2500;
+param.opts.disp=1;
+param.normalize = 1;
+param.normalize_type = 1;
+
+
+% Specifies where to save all results of the brain grid
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+param.title_ODF = ['ODF_Neigh_',num2str(param.neighborhood),'_ODFPower_',num2str(param.odfPow)];
+param.SaveDirectory = fullfile(param.HCPDatapath, 'BrainGraph_results', param.subject,param.title_ODF);
+param.structural = fullfile(param.HCPDatapath, param.subject, 'T1w');
+if ~exist(param.SaveDirectory, 'dir')
+ mkdir(param.SaveDirectory)
+end
\ No newline at end of file
diff --git a/Inputs/Inputs_fMRI.m b/Inputs/Inputs_fMRI.m
new file mode 100644
index 0000000..d95012f
--- /dev/null
+++ b/Inputs/Inputs_fMRI.m
@@ -0,0 +1,49 @@
+%%%%%% INPUT FUNCTIONAL VOLUMES PATH
+
+% filepath of preprocessed functional volumes:
+% should already be in matrix form (timepoints x number of voxel)
+
+param.functional = fullfile(param.HCPDatapath, param.subject,'functional');
+
+if ~isfield(param,'session')
+ funcName = 'BCVolumes_s6_100307_session_rfMRI_REST1_LR.mat';
+else
+ funcName = ['BCVolumes_s6_100307_session_',param.session,'.mat'];
+end
+
+% to detrend volumes (put 1 if volumes were not detrended in preprocessing)
+Detrend = 0;
+
+prefix = 's6';
+% Contains all inputs regarding functional MRI volumes which are subject-specific
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+%% Functional volumes
+
+% Loading functional data matrix
+load(fullfile(param.functional, funcName))
+
+directory = dir(fullfile(param.functional, [prefix,'*']));
+fHeader = spm_vol(fullfile(param.functional,directory(1).name));
+hdr=cbiReadNiftiHeader(fHeader.fname);
+
+
+% If the functional volumes were not detrended, we perform detrending
+% signal inpainting
+if Detrend
+ disp('Let us detrend volume first..')
+ CUTNUMBER=10;
+ SegmentLength = ceil(size(V,2) / CUTNUMBER);
+ for iCut=1:CUTNUMBER
+ if iCut~=CUTNUMBER
+ Segment = (iCut-1)*SegmentLength+1 : iCut*SegmentLength;
+ else
+ Segment = (iCut-1)*SegmentLength+1 : size(V,2);
+ end
+ V(:,Segment) = detrend(V(:,Segment));
+ V(:,Segment) = detrend(V(:,Segment),'constant');
+ fprintf('.');
+ end
+end
+NumScans = size(V,1);
\ No newline at end of file
diff --git a/Inputs/Inputs_fMRI_Mask.m b/Inputs/Inputs_fMRI_Mask.m
new file mode 100644
index 0000000..c491637
--- /dev/null
+++ b/Inputs/Inputs_fMRI_Mask.m
@@ -0,0 +1,38 @@
+
+% Paths
+
+param.EigPath = fullfile(param.HCPDatapath, param.subject, 'T1w','graph');
+
+param.SaveInpainting = fullfile(param.HCPDatapath, 'Inpainting_results');
+
+if ~exist(param.SaveInpainting,'dir')
+ mkdir(param.SaveInpainting)
+end
+
+% Reads gray matter mask from segmented T1 image
+Mask = spm_read_vols(spm_vol(fullfile(param.structural,...
+ 'c1T1w_acpc_dc_restore_1.25.nii')));
+
+Mask(Mask<0.3)=0;
+
+Mask = bwareaopen(logical(Mask),30);
+
+Mask = imfill(Mask,'holes');
+
+BCBrainMask = logical(Mask);
+
+param.smoothR = 6;
+
+load(fullfile(param.EigPath, 'indices_wb.mat'))
+
+param.indices = indices_wb;
+
+BC = find(BCBrainMask);
+
+[sim_ind,~,BC] = intersect(BC, param.indices);
+
+if ~isfield(param.WB,'A')
+ load(fullfile(param.EigPath,'A_wb.mat'))
+end
+
+
diff --git a/RunBrainGraph.m b/RunBrainGraph.m
new file mode 100644
index 0000000..309ce8d
--- /dev/null
+++ b/RunBrainGraph.m
@@ -0,0 +1,41 @@
+function param = RunBrainGraph(param)
+
+
+ % Date and time when the routines are called
+ param.date = strrep(strrep(datestr(datetime('now')),' ','_'),':','_');
+
+ param.constW = param.c_bandwidth;
+ % Original parameters
+ param_init = param;
+
+
+
+ disp(['Extracting whole brain graph using ODF-based design for ..',param.subject])
+
+
+ % Calls the mdh routine to construct the graph and return adjacency matrix
+
+ [param.WB, param.G] = main_braingraph(param);
+
+ [param.WB.A,param.WB.D]=slepNormalize(param.WB.A,param.normalize, param.normalize_type);
+
+ % Computing the eigendecomposition of the Laplacian
+
+ if param.percent
+ param.constW = round(param.bandwidth*param.G.N_wb);
+ else
+ param.constW = param.c_bandwidth;
+ end
+
+ disp('Taking the eigenvalues...')
+
+ [param.WB.Utr,param.WB.S]=slepEigsLaplacian(param.WB.A,param.constW,param.opts);
+
+ % Saves the eigenvalues and eigenvectors
+
+ SaveValues(param)
+
+ SaveToNifti(param, param.WB.Utr,'Spectrum_WB')
+
+
+end
\ No newline at end of file
diff --git a/Utilities/SaveToNifti.m b/Utilities/SaveToNifti.m
new file mode 100644
index 0000000..a07f9da
--- /dev/null
+++ b/Utilities/SaveToNifti.m
@@ -0,0 +1,25 @@
+function [] = SaveToNifti(param, Utr, cases)
+
+
+ wheretosave = fullfile(param.SaveDirectory, cases);
+ if ~exist(wheretosave, 'dir')
+ mkdir(wheretosave)
+ end
+
+ fHeader = spm_vol(fullfile(param.HCPDatapath,param.subject, 'T1w','extracted','new_brainmask.nii'));
+ hdr=cbiReadNiftiHeader(fHeader.fname);
+
+
+ indices = param.G.indices_wb;
+
+ for i = 1:size(Utr,2)
+
+ V = zeros(fHeader.dim);
+ V(indices) = Utr(:,i);
+ niftiFile = fullfile(wheretosave,[sprintf('Eigenmode_%.4d',i),'.nii']);
+
+ cbiWriteNifti(niftiFile,V,hdr,'float32');
+ end
+
+
+end
\ No newline at end of file
diff --git a/Utilities/SaveValues.m b/Utilities/SaveValues.m
new file mode 100644
index 0000000..5e31668
--- /dev/null
+++ b/Utilities/SaveValues.m
@@ -0,0 +1,15 @@
+function [] = SaveValues(param)
+
+
+ if ~exist(fullfile(param.SaveDirectory), 'dir')
+ mkdir(fullfile(param.SaveDirectory))
+ end
+
+
+ Utr = param.WB.Utr;
+ S = param.WB.S;
+
+ save(fullfile(param.SaveDirectory,['WB_Utr_',num2str(param.constW),'.mat']), 'Utr','-v7.3')
+ save(fullfile(param.SaveDirectory,['WB_S_',num2str(param.constW),'.mat']), 'S','-v7.3')
+
+end
\ No newline at end of file
diff --git a/Utilities/knee_pt.m b/Utilities/knee_pt.m
new file mode 100644
index 0000000..0603a03
--- /dev/null
+++ b/Utilities/knee_pt.m
@@ -0,0 +1,189 @@
+function [res_x, idx_of_result] = knee_pt(y,x,just_return)
+%function [res_x, idx_of_result] = knee_pt(y,x,just_return)
+%Returns the x-location of a (single) knee of curve y=f(x)
+% (this is useful for e.g. figuring out where the eigenvalues peter out)
+%
+%Also returns the index of the x-coordinate at the knee
+%
+%Parameters:
+% y (required) vector (>=3 elements)
+% x (optional) vector of the same size as y
+% just_return (optional) boolean
+%
+%If just_return is True, the function will not error out and simply return a Nan on
+%detected error conditions
+%
+%Important: The x and y don't need to be sorted, they just have to
+%correspond: knee_pt([1,2,3],[3,4,5]) = knee_pt([3,1,2],[5,3,4])
+%
+%Important: Because of the way the function operates y must be at least 3
+%elements long and the function will never return either the first or the
+%last point as the answer.
+%
+%Defaults:
+%If x is not specified or is empty, it's assumed to be 1:length(y) -- in
+%this case both returned values are the same.
+%If just_return is not specified or is empty, it's assumed to be false (ie the
+%function will error out)
+%
+%
+%The function operates by walking along the curve one bisection point at a time and
+%fitting two lines, one to all the points to left of the bisection point and one
+%to all the points to the right of of the bisection point.
+%The knee is judged to be at a bisection point which minimizes the
+%sum of errors for the two fits.
+%
+%the errors being used are sum(abs(del_y)) or RMS depending on the
+%(very obvious) internal switch. Experiment with it if the point returned
+%is not to your liking -- it gets pretty subjective...
+%
+%
+%Example: drawing the curve for the submission
+% x=.1:.1:3; y = exp(-x)./sqrt(x); [i,ix]=knee_pt(y,x);
+% figure;plot(x,y);
+% rectangle('curvature',[1,1],'position',[x(ix)-.1,y(ix)-.1,.2,.2])
+% axis('square');
+%
+
+%Food for thought: In the best of possible worlds, per-point errors should
+%be corrected with the confidence interval (i.e. a best-line fit to 2
+%points has a zero per-point fit error which is kind-a wrong).
+%Practially, I found that it doesn't make much difference.
+%
+%dk /2012
+
+
+
+%{
+
+% test vectors:
+
+[i,ix]=knee_pt([30:-3:12,10:-2:0]) %should be 7 and 7
+knee_pt([30:-3:12,10:-2:0]') %should be 7
+knee_pt(rand(3,3)) %should error out
+knee_pt(rand(3,3),[],false) %should error out
+knee_pt(rand(3,3),[],true) %should return Nan
+knee_pt([30:-3:12,10:-2:0],[1:13]) %should be 7
+knee_pt([30:-3:12,10:-2:0],[1:13]*20) %should be 140
+knee_pt([30:-3:12,10:-2:0]+rand(1,13)/10,[1:13]*20) %should be 140
+knee_pt([30:-3:12,10:-2:0]+rand(1,13)/10,[1:13]*20+rand(1,13)) %should be close to 140
+x = 0:.01:pi/2; y = sin(x); [i,ix]=knee_pt(y,x) %should be around .9 andaround 90
+[~,reorder]=sort(rand(size(x)));xr = x(reorder); yr=y(reorder);[i,ix]=knee_pt(yr,xr) %i should be the same as above and xr(ix) should be .91
+knee_pt([10:-1:1]) %degenerate condition -- returns location of the first "knee" error minimum: 2
+
+%}
+
+
+%set internal operation flags
+use_absolute_dev_p = true; %ow quadratic
+
+%deal with issuing or not not issuing errors
+issue_errors_p = true;
+if (nargin > 2 && ~isempty(just_return) && just_return)
+ issue_errors_p = false;
+end
+
+%default answers
+res_x = nan;
+idx_of_result = nan;
+
+%check...
+if (isempty(y))
+ if (issue_errors_p)
+ error('knee_pt: y can not be an empty vector');
+ end
+ return;
+end
+
+%another check
+if (sum(size(y)==1)~=1)
+ if (issue_errors_p)
+ error('knee_pt: y must be a vector');
+ end
+
+ return;
+end
+
+%make a vector
+y = y(:);
+
+%make or read x
+if (nargin < 2 || isempty(x))
+ x = (1:length(y))';
+else
+ x = x(:);
+end
+
+%more checking
+if (ndims(x)~= ndims(y) || ~all(size(x) == size(y)))
+ if (issue_errors_p)
+ error('knee_pt: y and x must have the same dimensions');
+ end
+
+ return;
+end
+
+%and more checking
+if (length(y) < 3)
+ if (issue_errors_p)
+ error('knee_pt: y must be at least 3 elements long');
+ end
+ return;
+end
+
+%make sure the x and y are sorted in increasing X-order
+if (nargin > 1 && any(diff(x)<0))
+ [~,idx]=sort(x);
+ y = y(idx);
+ x = x(idx);
+else
+ idx = 1:length(x);
+end
+
+%the code below "unwraps" the repeated regress(y,x) calls. It's
+%significantly faster than the former for longer y's
+%
+%figure out the m and b (in the y=mx+b sense) for the "left-of-knee"
+sigma_xy = cumsum(x.*y);
+sigma_x = cumsum(x);
+sigma_y = cumsum(y);
+sigma_xx = cumsum(x.*x);
+n = (1:length(y))';
+det = n.*sigma_xx-sigma_x.*sigma_x;
+mfwd = (n.*sigma_xy-sigma_x.*sigma_y)./det;
+bfwd = -(sigma_x.*sigma_xy-sigma_xx.*sigma_y) ./det;
+
+%figure out the m and b (in the y=mx+b sense) for the "right-of-knee"
+sigma_xy = cumsum(x(end:-1:1).*y(end:-1:1));
+sigma_x = cumsum(x(end:-1:1));
+sigma_y = cumsum(y(end:-1:1));
+sigma_xx = cumsum(x(end:-1:1).*x(end:-1:1));
+n = (1:length(y))';
+det = n.*sigma_xx-sigma_x.*sigma_x;
+mbck = flipud((n.*sigma_xy-sigma_x.*sigma_y)./det);
+bbck = flipud(-(sigma_x.*sigma_xy-sigma_xx.*sigma_y) ./det);
+
+%figure out the sum of per-point errors for left- and right- of-knee fits
+error_curve = nan(size(y));
+for breakpt = 2:length(y-1)
+ delsfwd = (mfwd(breakpt).*x(1:breakpt)+bfwd(breakpt))-y(1:breakpt);
+ delsbck = (mbck(breakpt).*x(breakpt:end)+bbck(breakpt))-y(breakpt:end);
+ %disp([sum(abs(delsfwd))/length(delsfwd), sum(abs(delsbck))/length(delsbck)])
+ if (use_absolute_dev_p)
+ % error_curve(breakpt) = sum(abs(delsfwd))/sqrt(length(delsfwd)) + sum(abs(delsbck))/sqrt(length(delsbck));
+ error_curve(breakpt) = sum(abs(delsfwd))+ sum(abs(delsbck));
+ else
+ error_curve(breakpt) = sqrt(sum(delsfwd.*delsfwd)) + sqrt(sum(delsbck.*delsbck));
+ end
+end
+
+%find location of the min of the error curve
+[~,loc] = min(error_curve);
+res_x = x(loc);
+idx_of_result = idx(loc);
+end
+
+
+
+
+
diff --git a/main_inpainting.m b/main_inpainting.m
new file mode 100644
index 0000000..7867037
--- /dev/null
+++ b/main_inpainting.m
@@ -0,0 +1,160 @@
+%% Methodological Pipeline for the Manuscript
+
+% Structural mediation of human brain activity revealed by white-matter
+% interpolation of fMRI
+%
+% Created by: Anjali Tarun
+% Date created: 21 February 2018
+%
+% DESCRIPTION:
+%
+% This is the main manuscript for performing white matter
+% inpainting from fMRI data using a voxel-wise brain connectome.
+% The pipeline works as follows:
+%
+% 1. The voxel-wise brain grid is first constructed. This
+% is done by first extracting the ODF from the diffusion
+% data. After which, the brain graph is built and the
+% Laplacian is computed.
+%
+% 2. Graph signal recovery is performed by minimizing a
+% cost function that is described in the manuscript.
+% The cost function balances between (1) minimizing the RSS
+% and retaining the original GM signals, and (2) imposing
+% smoothness with respect to the structure of the graph.
+%
+%
+% The pipeline is adjusted to take HCP folder as input and outputs
+% the brain graph adjacency matrix, Laplacian, and the corresponding
+% interpolated volumes for each subject and session.
+%
+% Input: HCP folder containing subject IDs
+% The HCP folder should contain (1)preprocessed diffusion data (2)
+% processed fMRI volumes (3) Extended structural volumes
+% Output: Voxel-brain graph and interpolated volumes
+%
+% The struct 'param' contains all the parameter information and paths
+%
+
+%% Specifies code paths and obtain data paths
+
+param = Addpaths();
+
+% Takes the set of subjects to be considered
+param.Subjects = {'100307'};
+
+param.Sessions = {'rfMRI_REST1_LR','rfMRI_REST1_RL','rfMRI_REST2_LR',...
+ 'rfMRI_REST2_RL','tfMRI_MOTOR_LR','tfMRI_MOTOR_RL'};
+
+
+for iS = 1:length(param.Subjects)
+
+ param.subject = param.Subjects{iS};
+ disp(['Constructing brain graph for subject..',param.subject])
+
+
+ %% Construct the voxel-level brain grid from diffusion MRI
+
+ % Calls for the parameters required for constructing the brain grid
+ Inputs_BrainGrid
+
+ % Runs the construction of the brain graph
+ param = RunBrainGraph(param);
+
+ %% Runs the graph signal inpainting
+
+
+ % First checks if the brain graph is already constructed
+ if ~isfield(param.WB,'A')
+ disp('Brain grid not available. Run construction of brain graph first!')
+ end
+
+ % Load parameters for functional volume masks
+ Inputs_fMRI_Mask
+
+ % Initialize L2 norm problem
+
+ A_wb = param.WB.A;
+
+ n = size(A_wb,1);
+
+ % Construct the indicator matrix
+ B = sparse(param.BC,param.BC,ones(1,length(param.BC)),n,n);
+
+ BB = B'*B;
+
+ [A_n,~] = slepNormalize(A_wb,1,1);
+
+ L = speye(size(A_n,1)) - A_n;
+
+ % Computes for the parameter lambda to use
+ lambda = CheckParameterLambda(param,BB,L,BC);
+
+ %% Actual run of the signal inpainting after initialization
+
+ for ses = 1:length(param.Sessions)
+
+ param.session = param.Sessions{ses};
+
+ % Load inputs for functional MRI
+ % Pipeline assumes that the fMRI volumes have already been
+ % preprocessed -- with prefix s6 (smooothed, Gaussian FWHM 6 mm)
+
+ Inputs_fMRI
+
+ disp(['Runing signal recovery for subject..',param.subject, 'and session', param.session])
+
+ if ~exist(fullfile(param.SaveInpainting, ['GSR_Lambda',num2str(lambda),'_',param.subject,'_session_',param.session,'_smoothing',num2str(param.smoothR),'.mat']),'file')
+
+ disp(['Running for lambda..',num2str(lambda)])
+
+ L_mat = BB+lambda*L;
+
+ L1 = ichol(L_mat);
+
+ V_new = zeros(NumScans,size(L1,1));
+
+
+ % Main loop for signal inpainting
+ % This may take a while -- alternatively one can run this
+ % in parallel using parfor
+
+ for i = 1:NumScans
+ disp(['Running for iteration..',num2str(i)])
+ xsignal = zeros(size(L1,1),1);
+
+ xsignal(BC) = V(i,BC);
+
+ [V_new(i,:),~,~,~,~] = pcg(L_mat,xsignal,1e-8,100,L1,L1');
+
+ end
+
+
+ V = V_new;
+ clear V_new
+ save(fullfile(param.SaveInpainting, ['GSR_Lambda',num2str(lambda),'_',param.subject,'_session_',param.session,'_smoothing',num2str(param.smoothR),'.mat']),'V','-v7.3')
+% clear V;
+ param = rmfield(param,'session');
+
+ end
+
+ end
+
+
+
+ %% Visualize inpainted signal
+ % save volumes here
+ path = fullfile(param.SaveInpainting, param.subject);
+ if ~exist(path, 'dir')
+ mkdir(path)
+ end
+
+ for iter = 1:10%NumScans
+
+ vol = zeros(fHeader.dim);
+ vol(indices_wb) = V(iter,:);
+ nfile = fullfile(path,['Inpainting_',param.session,sprintf('_%.3d',iter),'.nii']);
+ cbiWriteNifti(nfile,vol,hdr,'float32');
+ end
+
+end
\ No newline at end of file