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