Compute matrix to predict directional diffusion in each voxel from fibers fe = feConnectomeBuildModel(fe) INPUTS: fe - An fe structure, see feCreate.m See also: feFitModel.m, feComputePredictedSignal.m Example: Copyright (2013-2014), Franco Pestilli, Stanford University, pestillifranco@gmail.com. -- The LiFE Model -- LiFE stands for Linear Fascicle Evaluation. The quality of a white-matter connectome is evaluated by using the connectome to model the diffusion signal in the volume of the connectome. The LiFE model is set up in the following way. The diffusion signal in the voxels comprising the volume of the connectome is modelled as a linear equation of the following form: dSig = [M Miso] * w Where dSig is a nVoxels * nBvecs x 1 vector containing the predicted diffusion signal in each voxel and each diffusion direction measured. M is a block matrix that has the following shape: M = [fiber_component, isotropic_component] Where fiber_component is a nVoxels * nBvecs x nFibers matrix containing the contribution of each fiber to the diffusion signal in each voxel and direction. This matrix is sparse, because if a fiber does not pass through a voxel, the contribution in that combination is set to 0. The isotropic_component accounts for the part of the diffusion signal that is not generated by the fiber, but rather by other elements of the tissue in that voxel and freely diffusing water. The fiber_component is the *sum* over contributions from individual unique fibers in the voxel. The isotropic_component is simply the mean diffusion signal over all directions in the voxel (assumed to be spherical in shape...). When a fiber passes through a voxel, the element of M at (voxel * bvec, fiber) will contain the following the predicted diffusion signal as follows: S = S0 exp(-bval*(bvec*Q*bvec)) Where Q is the quadratic form of the tensor and S0 is the diffusion signal measured with no diffusion weighting (baseline measured in B0 scans). M is represented as a Matlab 'sparse' matrix. This saves memory and speeds up calculations. However, it makes the computations required in building the matrix a little strange.
0001 function fe = feConnectomeBuildModel(fe) 0002 % Compute matrix to predict directional diffusion in each voxel from fibers 0003 % 0004 % fe = feConnectomeBuildModel(fe) 0005 % 0006 % INPUTS: fe - An fe structure, see feCreate.m 0007 % 0008 % See also: feFitModel.m, feComputePredictedSignal.m 0009 % 0010 % Example: 0011 % 0012 % Copyright (2013-2014), Franco Pestilli, Stanford University, pestillifranco@gmail.com. 0013 % 0014 % -- The LiFE Model -- 0015 % 0016 % LiFE stands for Linear Fascicle Evaluation. The quality of a white-matter 0017 % connectome is evaluated by using the connectome to model the diffusion 0018 % signal in the volume of the connectome. 0019 % 0020 % The LiFE model is set up in the following way. The diffusion signal in 0021 % the voxels comprising the volume of the connectome is modelled as a 0022 % linear equation of the following form: 0023 % 0024 % dSig = [M Miso] * w 0025 % 0026 % Where dSig is a nVoxels * nBvecs x 1 vector containing the predicted 0027 % diffusion signal in each voxel and each diffusion direction measured. M 0028 % is a block matrix that has the following shape: 0029 % 0030 % M = [fiber_component, isotropic_component] 0031 % 0032 % Where fiber_component is a nVoxels * nBvecs x nFibers matrix containing 0033 % the contribution of each fiber to the diffusion signal in each voxel and 0034 % direction. This matrix is sparse, because if a fiber does not pass 0035 % through a voxel, the contribution in that combination is set to 0. The 0036 % isotropic_component accounts for the part of the diffusion signal that is 0037 % not generated by the fiber, but rather by other elements of the tissue in 0038 % that voxel and freely diffusing water. 0039 % 0040 % The fiber_component is the *sum* over contributions from individual unique 0041 % fibers in the voxel. The isotropic_component is simply the mean diffusion 0042 % signal over all directions in the voxel (assumed to be spherical in 0043 % shape...). 0044 % 0045 % When a fiber passes through a voxel, the element of M at (voxel * bvec, fiber) 0046 % will contain the following the predicted diffusion signal as follows: 0047 % 0048 % S = S0 exp(-bval*(bvec*Q*bvec)) 0049 % 0050 % Where Q is the quadratic form of the tensor and S0 is the diffusion 0051 % signal measured with no diffusion weighting (baseline measured in B0 0052 % scans). 0053 % 0054 % M is represented as a Matlab 'sparse' matrix. This saves memory and 0055 % speeds up calculations. However, it makes the computations required in 0056 % building the matrix a little strange. 0057 0058 if notDefined('fe'), error('LiFE (fe = feCreate) struct needed'); end 0059 if ~isfield(fe,'life') 0060 error('LiFE - the field ''life'' is necessary in the fe structure.') 0061 end 0062 0063 disp('LiFE - Building the connectome model...'); 0064 tic 0065 % Indexes of actually used voxels, this can be less than the number of 0066 % voxels int he ROI in case some voxels have no fibers in them 0067 usedVoxels = feGet(fe,'usedVoxels'); 0068 nVoxels = length(usedVoxels); 0069 nBvecs = feGet(fe,'nBvecs'); 0070 vox_sparse_pSig = cell(nVoxels,1); 0071 0072 % Generate the signal predicted in each voxel by the connectome. 0073 % 0074 % This for-loop is written in such a way that matlab (ver later than 7.9) 0075 % will run it in parallel if parallel processing is allowed. 0076 fprintf('LiFE - Predicting diffusion signal in %i voxel...\n',nVoxels); 0077 0078 feOpenLocalCluster 0079 parfor vv = 1:nVoxels 0080 num_unique_fibers = feGet(fe,'unique f num',usedVoxels(vv)); 0081 0082 % This returns a matrix that is size nBvecs x num_unique_fibers 0083 voxelPSignal = feComputeVoxelSignal(fe,usedVoxels(vv)); 0084 0085 % Fibers in the connectome determine the directional signal in the voxel 0086 % signal, not the mean signal in the voxel. Here we first demean the 0087 % voxel signal we will predict. 0088 % 0089 % demeaned_pSig = voxelPSignal - repmat(mean(voxelPSignal, 1),nBvecs,1); 0090 % 0091 % The mean will be predicted by the Miso part of the matrix, not the 0092 % Mfiber part. 0093 % 0094 % Then we reorganize the demeaned signal into a column vector 0095 % 0096 % vox_sparse_pSig{vv} = reshape(demeaned_pSig', num_unique_fibers * nBvecs, 1) ; 0097 % 0098 % Somehow this column vector ends up occupying the right parts of the 0099 % Mfiber matrix when we are done. That miracle happens below. 0100 vox_sparse_pSig{vv} = reshape((voxelPSignal - repmat(mean(voxelPSignal, 1),nBvecs,1))', ... 0101 num_unique_fibers * nBvecs, 1) ; 0102 end 0103 fprintf('LiFE - Prediction computed in: %2.3fs.\n',toc); 0104 0105 tic 0106 fprintf('LiFE - Allocating the model prediction...') 0107 0108 % These variables will be used to generate the sparse matrix indices and 0109 % values (from the dw signal in each voxel-direction). We pre-allocate here 0110 % for speed. We assign the values in the loop below. More explanation 0111 % follows. 0112 0113 % Each row and column of the Mfiber matrix represents a unique contribution 0114 % to the signal in some voxel (group of rows) from a fiber (column). The 0115 % total number of elements of Mfiber is the number of unique fibers 0116 % (columns), times the number of voxels, times the number of directions 0117 % (rows). 0118 nUniqueFibersInEachVoxel = feGet(fe,'unique f num'); 0119 0120 % number of non-zero elements in the MFiber matrix 0121 M_siz = sum(nUniqueFibersInEachVoxel(1:nVoxels))*nBvecs; 0122 0123 % This is a vector that will contain the diffusion signals 0124 M_signal = zeros(M_siz, 1); 0125 0126 % These define the row/col entry in the Mfiber matrix where the diffusion 0127 % signal needs to be stored. 0128 M_rows = zeros(M_siz, 1); 0129 M_cols = zeros(M_siz, 1); 0130 0131 % Now, pre-allocate space for the measured diffusion signal. This will be 0132 % filled up in the following loop, too. 0133 fe = feSet(fe,'dSig measured',zeros(1,nVoxels*nBvecs)); 0134 0135 % The following lines create the *sparse-matrix* indexing into the M matrix. 0136 % Plus generate a vector of measured signal, full signal and demeaned. 0137 % 0138 % The following operations cannot be run in parallel because we address 0139 % vectors and matrices on the fly. 0140 end_idx = 0; 0141 for vv = 1:nVoxels 0142 num_unique_fibers = feGet(fe,'unique f num',usedVoxels(vv)); 0143 index_unique_fibers = cell2mat(feGet(fe,'unique f',usedVoxels(vv))); 0144 0145 % The first return is a binary vector of the locations of the rows of 0146 % Mfiber corresponding to the current voxel. We do the find to get the 0147 % indices of these rows. 0148 dense_rows = find(feGet(fe,'voxel rows',vv)); 0149 0150 % We determine the column/row combinations for the non-zero elements in 0151 % the part of M corresponding to the current voxel: 0152 sparse_rows = kron(dense_rows,ones(num_unique_fibers,1)); 0153 sparse_cols = repmat(index_unique_fibers, nBvecs, 1); 0154 0155 % We need to assign the sparse_* variables into the pre-allocated M_* 0156 % variables in the right place. 0157 start_idx = end_idx + 1; 0158 end_idx = end_idx + length(sparse_rows); 0159 0160 M_rows(start_idx:end_idx) = sparse_rows; 0161 M_cols(start_idx:end_idx) = sparse_cols; 0162 M_signal(start_idx:end_idx) = vox_sparse_pSig{vv}; 0163 0164 % Reorganize the diffusion data for each voxel into a long vector: 0165 % nDirs X nVoxels. 0166 fe.life.dSig(dense_rows) = feGet(fe,'diffusion signal in voxel',usedVoxels(vv)); 0167 end 0168 0169 % Install the matrix in the fe structure. 0170 fe = feSet(fe,'Mfiber',sparse(M_rows, M_cols, M_signal)); 0171 0172 fprintf('process completed in %2.3fs.\n',toc) 0173 disp('LiFE - DONE Building the connectome model.'); 0174 0175 return