dmri_connect_parm_calc

 Compute binarised anatomical connectivity matrix and fiber length matrix

 [Units]
       v: mm/ms (=m/s)
     tau: ms
   Delta: ms
 (f_leng: mm)

 2012/12/03 M.Fukushima
 2012/12/05 M.Fukushima
 2013/02/12 M.Fukushima
 2014/06/02 rhayashi      Nlabel, ofile -> parcel_dir
 2014/12/25 rhayashi      added parcels_file for saving parcellation info
                                (brain_file, vbmeg_area_ix)
 2017/10/25 Y. Takeda     added c_streng for saving connectivity strength info
 2018/10/23 rhayashi      save 'Area_cortex'  

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function [Ind, Delta] = dmri_connect_parm_calc(...
0002   parcels_file, parcel_dir, cmatfile, output_dmri_file, ...
0003   thre, v, tau, uni_flag, save_flag)
0004 % Compute binarised anatomical connectivity matrix and fiber length matrix
0005 %
0006 % [Units]
0007 %       v: mm/ms (=m/s)
0008 %     tau: ms
0009 %   Delta: ms
0010 % (f_leng: mm)
0011 %
0012 % 2012/12/03 M.Fukushima
0013 % 2012/12/05 M.Fukushima
0014 % 2013/02/12 M.Fukushima
0015 % 2014/06/02 rhayashi      Nlabel, ofile -> parcel_dir
0016 % 2014/12/25 rhayashi      added parcels_file for saving parcellation info
0017 %                                (brain_file, vbmeg_area_ix)
0018 % 2017/10/25 Y. Takeda     added c_streng for saving connectivity strength info
0019 % 2018/10/23 rhayashi      save 'Area_cortex'
0020 %
0021 % Copyright (C) 2011, ATR All Rights Reserved.
0022 % License : New BSD License(see VBMEG_LICENSE.txt)
0023 
0024 %
0025 % --- Previous check
0026 %
0027 if nargin < 4
0028     error('Please check input arguments.');
0029 end
0030 
0031 if ~exist('thre', 'var')
0032     thre = 1e-4;
0033 end
0034 if ~exist('v', 'var')
0035     v = 6; % mm/ms = (m/s)
0036 end
0037 if ~exist('tau', 'var')
0038     tau = 27; % mm
0039 end
0040 if ~exist('uni_flag', 'var')
0041     uni_flag = 0;
0042 end
0043 if ~exist('save_flag', 'var')
0044     save_flag = 1;
0045 end
0046 output_dir = fileparts(output_dmri_file);
0047 if exist(output_dir, 'dir') ~= 7
0048     mkdir(output_dir);
0049 end
0050 
0051 %
0052 % --- Main Procedure
0053 %
0054 
0055 % Load parcellation information
0056 warning('off', 'MATLAB:load:variableNotFound');
0057 load(parcels_file, 'brain_file', 'vbmeg_area_ix', 'membershipmat', 'Area_cortex');
0058 warning('on',  'MATLAB:load:variableNotFound');
0059 
0060 % Load matrix
0061 load(cmatfile, 'c_count', 'f_leng', 'Ntrack', 't_count');
0062 
0063 % Load mask
0064 parcel_file = fullfile(parcel_dir, filesep, 'parcel_mm_coord.mat');
0065 
0066 load(parcel_file, 'label', 'Nlabel');
0067 c_label = zeros(Nlabel,1);
0068 for ii = 1:Nlabel
0069   c_label(ii) = sum(label==ii);
0070 end
0071 
0072 fprintf('calculating connectivity matrix...');
0073 
0074 %%% Metrics calculation
0075 % Connectivity strength
0076 c_streng = c_count.*repmat(c_label,[1 Nlabel])' ...
0077   ./(repmat(c_label,[1 Nlabel]).*repmat(t_count,[1 Nlabel])');
0078 c_streng(isnan(c_streng)) = 0; % for missing label voxels
0079 c_streng(isinf(c_streng)) = 0;
0080 
0081 % Symmetrize
0082 [c_streng, f_leng] = make_symmetry_matrix(c_streng, f_leng);
0083 
0084 % Binarised anatomical connectivity matrix
0085 Ind = sparse(c_streng>thre);
0086 Ind = Ind - diag(diag(Ind));
0087 
0088 % Fiber length matrix
0089 Delta = (f_leng.*Ind)/v + tau*Ind;
0090 
0091 % For diagonal entries
0092 Ind = Ind + speye(size(Ind,1));
0093 Delta = Delta + tau*speye(size(Ind,1));
0094 
0095 dmri_parm = [];
0096 dmri_parm.Nlabel = Nlabel;
0097 dmri_parm.ofile = parcel_file;
0098 dmri_parm.cmatfile = cmatfile;
0099 dmri_parm.thre = thre;
0100 dmri_parm.v = v;
0101 dmri_parm.tau = tau;
0102 dmri_parm.uni_flag = uni_flag;
0103 
0104 if exist('uni_flag','var')
0105   if uni_flag
0106     % Uniform Delta
0107     Delta = tau*(Delta>0);
0108   end
0109 end
0110 
0111 if exist('save_flag','var')
0112   if save_flag
0113     connectivity_matrix = Ind;
0114     delay_ms_matrix    = Delta;
0115 
0116     fiber_length_m_matrix = zeros(size(delay_ms_matrix));
0117     ix = find(delay_ms_matrix);
0118     fiber_length_m_matrix(ix) = (delay_ms_matrix(ix) - dmri_parm.tau)/1000 * dmri_parm.v;
0119     fiber_length_m_matrix = sparse(fiber_length_m_matrix); 
0120 
0121     % Save
0122     save(output_dmri_file,'connectivity_matrix', ...
0123                           'delay_ms_matrix', ...
0124                           'fiber_length_m_matrix', ...
0125                           'brain_file', 'vbmeg_area_ix',...
0126                           'membershipmat', 'dmri_parm','c_streng');
0127                       
0128     if exist('Area_cortex', 'var')
0129         vb_save(output_dmri_file, 'Area_cortex');
0130     end
0131   end
0132 end
0133 fprintf('done\n');