From b105ecb7d17312759f42091804ccf16fb1a7c4f4 Mon Sep 17 00:00:00 2001 From: Brent McPherson Date: Wed, 11 Nov 2020 17:34:11 -0500 Subject: [PATCH 01/19] initial commit to new fork from old (?) brain-life repo. just comment note, pre changes --- life/fe/feBuildDictionaries.m | 6 ++++-- life/fe/feConnectomeEncoding.m | 3 +++ life/fe/feConnectomeInit.m | 6 +++--- life/fe/feGet.m | 2 +- 4 files changed, 11 insertions(+), 6 deletions(-) diff --git a/life/fe/feBuildDictionaries.m b/life/fe/feBuildDictionaries.m index 3f84152..db2f790 100755 --- a/life/fe/feBuildDictionaries.m +++ b/life/fe/feBuildDictionaries.m @@ -44,6 +44,7 @@ %DictTensors = zeros(9,Norient); % Initialize Tensors Dictionary matrix D = diag(fe.life.modelTensor); % diagonal matix with diffusivities +% THIS WILL GO IN THE LOOP BELOW %[ shells, ~, shellindex ] = unique(round(bvals)); % round to get around noise in shells - FRAGILE FIX @@ -54,11 +55,12 @@ for j=1:Norient [Rot,~, ~] = svd(orient(:,j)); % Compute the eigen vectors of the kernel orientation - Q = Rot*D*Rot'; + Q = Rot*D*Rot'; % THIS GETS PULLED BY SHELL %DictTensors(:,j) = Q(:); Dict(:,j) = exp(- bvals .* diag(bvecs*Q*bvecs')); % Compute the signal contribution of a fiber in the kernel orientation divided S0 + % THIS CAN BE MODIFIED - HARD-CODED TENSOR - COULD BE DKI? - % for every shell + % for every shell - THIS INCORPORATES EVERYTHING BUT ROT for k=1:size(ubv, 1) si = ubi == ubv(k); % find the indices for the shell DictSig(si,j) = Dict(si,j) - median(Dict(si,j)); % demedianed signal by shell (used to demean) diff --git a/life/fe/feConnectomeEncoding.m b/life/fe/feConnectomeEncoding.m index f4e3f38..e0d3463 100755 --- a/life/fe/feConnectomeEncoding.m +++ b/life/fe/feConnectomeEncoding.m @@ -144,6 +144,9 @@ roi_coords = feGet(fe,'roicoords'); roi_ind = sub2ind(imgsize,roi_coords(:,1)',roi_coords(:,2)',roi_coords(:,3)'); +% SIGNAL NORMALIZATION NEEDS TO HAPPEN BY SHELL IN THIS OPERATION +% IT DOES NOT CURRENTLY + % Restrict the Sparse Tensor and matrix A to the ROI voxels only Phi = Phi(:,roi_ind,:); % reduce tensor in 3rd dimension to roi voxels only A = A(roi_ind,:); diff --git a/life/fe/feConnectomeInit.m b/life/fe/feConnectomeInit.m index 6a36de7..3ee98d1 100755 --- a/life/fe/feConnectomeInit.m +++ b/life/fe/feConnectomeInit.m @@ -81,9 +81,9 @@ axialDiffusion = 1; radialDiffusion = 0; end -dParms(1) = axialDiffusion; -dParms(2) = radialDiffusion; -dParms(3) = radialDiffusion; +dParms(1) = axialDiffusion; % DEFINE THESE AS SEPARATE PER SHELL +dParms(2) = radialDiffusion; % FIGURE OUT HOW TO PASS A DIAGONAL ARRAY PER SHELL +dParms(3) = radialDiffusion; % MAKE DEFAULTS KIND TO ABSENCE, PASS nX3 ARRAY TO CONVERT TO DIAGONAL, USED IN feBuildDictionary Nphi = N; Ntheta = N; fe = feSet(fe,'model tensor',dParms); diff --git a/life/fe/feGet.m b/life/fe/feGet.m index 2194e17..8703a78 100755 --- a/life/fe/feGet.m +++ b/life/fe/feGet.m @@ -1310,7 +1310,7 @@ val = sqrt(mean((measured - predicted).^2,1)); val = val(feGet(fe,'voxelsindices',varargin)); - case {'voxrmses0norm'} + case {'voxrmses0norm'} % CAN THIS BE OPTIMIZED FOR A MULTISHELL? % A volume of RMSE normalized by the S0 value in each voxel. rmse = feGet(fe,'vox rmse'); s0 = feGet(fe,'b0signalimage')'; From 16babd75a1e60f7b39adf57b11dae353744ae984 Mon Sep 17 00:00:00 2001 From: Brent McPherson Date: Fri, 13 Nov 2020 02:23:23 -0500 Subject: [PATCH 02/19] added multishell support / partial implementation of kutosis forward model --- life/fe/feBuildDictionaries.m | 88 ++++++++++++++++++++++++++++++----- life/fe/feConnectomeInit.m | 23 ++++++--- 2 files changed, 92 insertions(+), 19 deletions(-) diff --git a/life/fe/feBuildDictionaries.m b/life/fe/feBuildDictionaries.m index db2f790..699ac65 100755 --- a/life/fe/feBuildDictionaries.m +++ b/life/fe/feBuildDictionaries.m @@ -14,7 +14,7 @@ % Cesar F. Caiafa ccaiafa@gmail.com tic -fprintf(['\n[%s] Computing demeaned and non-demeaned difussivities dictionaries in a (',num2str(Nphi),'x',num2str(Ntheta),')-grid', ' ...'],mfilename); +fprintf(['\n[%s] Computing demeaned and non-demeaned diffusivities dictionaries in a (',num2str(Nphi),'x',num2str(Ntheta),')-grid', ' ...'],mfilename); fprintf('took: %2.3fs.\n',toc) % Compute orientation vectors @@ -43,31 +43,95 @@ DictSig = zeros(nBvecs,Norient); % Initialize Signal Dictionary matrix %DictTensors = zeros(9,Norient); % Initialize Tensors Dictionary matrix -D = diag(fe.life.modelTensor); % diagonal matix with diffusivities -% THIS WILL GO IN THE LOOP BELOW +% catch Kurtosis estimates for debugging +akc = zeros(nBvecs,Norient); +KDict = zeros(nBvecs,Norient); %[ shells, ~, shellindex ] = unique(round(bvals)); % round to get around noise in shells - FRAGILE FIX +% pull the shell information ubv = feGet(fe, 'nshells'); ubi = feGet(fe, 'shellindex'); % Compute each dictionary column for a different kernel orientation for j=1:Norient - [Rot,~, ~] = svd(orient(:,j)); % Compute the eigen vectors of the kernel orientation - Q = Rot*D*Rot'; % THIS GETS PULLED BY SHELL - %DictTensors(:,j) = Q(:); - Dict(:,j) = exp(- bvals .* diag(bvecs*Q*bvecs')); % Compute the signal contribution of a fiber in the kernel orientation divided S0 - % THIS CAN BE MODIFIED - HARD-CODED TENSOR - COULD BE DKI? + % Compute the eigen vectors of the kernel orientation + [Rot,~, ~] = svd(orient(:,j)); - % for every shell - THIS INCORPORATES EVERYTHING BUT ROT - for k=1:size(ubv, 1) - si = ubi == ubv(k); % find the indices for the shell - DictSig(si,j) = Dict(si,j) - median(Dict(si,j)); % demedianed signal by shell (used to demean) + % for every shell + for k=1:size(ubv, 1) + + % find the indices for the shell + si = ubi == ubv(k); + + % create diagonal matix with diffusivities for current shell + % this assumes tensor fits for shell are entered in the order this will parse them in + D = diag(fe.life.modelTensor(k,:)); + + % estimate Q for tensor values in shell + Q = Rot*D*Rot'; + + % logical starts here: if kurtosis run these steps, else just fill in tensor + + % mean diffusivity of tensor? + md = mean(diag(D)); + + % pull tensor / kurtosis parameters for hard indexing of equations + dt = fe.life.modelTensor(k,:); + kt = fe.life.modelKurtosis; + + % apparent diffusion coefficient from dipy - + adc = bvecs(si,1) .* bvecs(si,1) * dt(1) + ... + 2 * bvecs(si,1) .* bvecs(si,2) * dt(2) + ... + bvecs(si,2) .* bvecs(si,2) .* dt(3) + ... + 2 * bvecs(si,1) .* bvecs(si,3) * 0 + ... + 2 * bvecs(si,2) .* bvecs(si,3) * 0 + ... + bvecs(si,3) .* bvecs(si,3) * 0; + + % apparent diffusion variance from dipy + % reordered to ensure match of kt params from mrtrix3 + adv = ... + bvecs(si,1) .* bvecs(si,1) .* bvecs(si,1) .* bvecs(si,1) .* kt(1) + ... % xxxx + bvecs(si,2) .* bvecs(si,2) .* bvecs(si,2) .* bvecs(si,2) .* kt(2) + ... % yyyy + bvecs(si,3) .* bvecs(si,3) .* bvecs(si,3) .* bvecs(si,3) .* kt(3) + ... % zzzz + 4 * bvecs(si,1) .* bvecs(si,1) .* bvecs(si,1) .* bvecs(si,2) .* kt(4) + ... % xxxy + 4 * bvecs(si,1) .* bvecs(si,1) .* bvecs(si,1) .* bvecs(si,3) .* kt(5) + ... % xxxz + 4 * bvecs(si,1) .* bvecs(si,2) .* bvecs(si,2) .* bvecs(si,2) .* kt(6) + ... % xyyy + 4 * bvecs(si,1) .* bvecs(si,3) .* bvecs(si,3) .* bvecs(si,3) .* kt(7) + ... % xzzz + 4 * bvecs(si,2) .* bvecs(si,2) .* bvecs(si,2) .* bvecs(si,3) .* kt(8) + ... % yyyz + 4 * bvecs(si,2) .* bvecs(si,3) .* bvecs(si,3) .* bvecs(si,3) .* kt(9) + ... % yzzz + 6 * bvecs(si,1) .* bvecs(si,1) .* bvecs(si,2) .* bvecs(si,2) .* kt(10) + ... % xxyy + 6 * bvecs(si,1) .* bvecs(si,1) .* bvecs(si,3) .* bvecs(si,3) .* kt(11) + ... % xxzz + 6 * bvecs(si,2) .* bvecs(si,2) .* bvecs(si,3) .* bvecs(si,3) .* kt(12) + ... % yyzz + 12 * bvecs(si,1) .* bvecs(si,1) .* bvecs(si,2) .* bvecs(si,3) .* kt(13) + ... % xxyz + 12 * bvecs(si,1) .* bvecs(si,2) .* bvecs(si,2) .* bvecs(si,3) .* kt(14) + ... % xyyz + 12 * bvecs(si,1) .* bvecs(si,2) .* bvecs(si,3) .* bvecs(si,3) .* kt(15); % xyzz + + % zero out noisy (bad) adc estimates + adc(adc < 0) = 0; + + % estimate apparent kutosis coefficient - store in temp var + takc = adv .* ((md / adc).^2)'; % always finds same param at 20/50, rest are zeros? + + % zero out noisy (bad) akc estimates - use tmp to avoid indexing of large matrix + takc(takc < -3/7) = -3/7; + + akc(si,j) = takc; + + % Compute the signal contribution of a fiber in the kernel orientation divided S0 + Dict(si,j) = exp(-bvals(si) .* diag(bvecs(si,:)*Q*bvecs(si,:)')); + KDict(si,j) = exp(-bvals(si) .* diag(bvecs(si,:)*Q*bvecs(si,:)') + (-bvals(si).^2 .* diag(bvecs(si,:)*Q*bvecs(si,:)').^2 .* akc(si,j))/6); + + % demedianed signal by shell (used to demean) + DictSig(si,j) = Dict(si,j) - median(Dict(si,j)); + end end +disp('Done.'); + fe = feSet(fe,'dictionary parameters',{Nphi,Ntheta,orient,Dict,DictSig}); end \ No newline at end of file diff --git a/life/fe/feConnectomeInit.m b/life/fe/feConnectomeInit.m index 3ee98d1..5941ec9 100755 --- a/life/fe/feConnectomeInit.m +++ b/life/fe/feConnectomeInit.m @@ -60,7 +60,7 @@ fe = feConnectomeSetDwi(fe,dwiFileRepeated,1); end -%% Anatomy Install the path tot he anatomical high-resolution file. +%% Anatomy Install the path to the anatomical high-resolution file. if ~notDefined('anatomyFile') fe = feSet(fe,'anatomy file',anatomyFile); end @@ -74,19 +74,28 @@ %tic if ~isempty(varargin) N = varargin{1}(1); - axialDiffusion = varargin{2}(1); - radialDiffusion = varargin{2}(2); + tensorParams = varargin{2}; + kurtosisPars = varargin{3}; + %axialDiffusion = varargin{2}(1); + %radialDiffusion = varargin{2}(2); + %radialDiffusion2 = varargin{2}(3); else % Default to stick and ball N = 360; + % very default tensor settings axialDiffusion = 1; radialDiffusion = 0; + tensorModel = [ axialDiffusion radialDiffusion radialDiffusion ]; + tensorParams = []; + kurtosisPars = []; + nshells = feGet(fe, 'nshells'); + for i=1:nshells + tensorParams = [ tensorParams; tensorModel ]; + end end -dParms(1) = axialDiffusion; % DEFINE THESE AS SEPARATE PER SHELL -dParms(2) = radialDiffusion; % FIGURE OUT HOW TO PASS A DIAGONAL ARRAY PER SHELL -dParms(3) = radialDiffusion; % MAKE DEFAULTS KIND TO ABSENCE, PASS nX3 ARRAY TO CONVERT TO DIAGONAL, USED IN feBuildDictionary Nphi = N; Ntheta = N; -fe = feSet(fe,'model tensor',dParms); +fe = feSet(fe,'model tensor',tensorParams); % TENSOR PARAMS NEEDS TO BE SANITIZED BETTER +fe.life.modelKurtosis = kurtosisPars; % add as feSet fxn? fe = feBuildDictionaries(fe,Nphi,Ntheta); %% Encode Connectome in multidimensional tensor From 9f498ef0d8ab0fdc90a10a394852be8523db5802 Mon Sep 17 00:00:00 2001 From: bcmcpher Date: Fri, 13 Nov 2020 02:30:17 -0500 Subject: [PATCH 03/19] added / modified files for new virtual lesion test --- life/compute/feComputeEvidence.m | 4 +- life/compute/feComputeEvidence_new.m | 277 ++++++++++++++++++++++++++ life/compute/feComputeEvidence_norm.m | 119 +++++++++-- 3 files changed, 386 insertions(+), 14 deletions(-) create mode 100644 life/compute/feComputeEvidence_new.m diff --git a/life/compute/feComputeEvidence.m b/life/compute/feComputeEvidence.m index d81890b..c317fb2 100755 --- a/life/compute/feComputeEvidence.m +++ b/life/compute/feComputeEvidence.m @@ -1,4 +1,4 @@ -function se = feComputeEvidence(rmse1,rmse2) +function se = feComputeEvidence(rmse1,rmse2,nbins) % Computes a series of distance metrics between two RMSE distribtions. % % Compute summary statistics on to characterize the lesion,: @@ -20,7 +20,7 @@ % Histograms se.xrange(1) = min([se.lesion.rmse.all se.nolesion.rmse.all]); se.xrange(2) = max([se.lesion.rmse.all se.nolesion.rmse.all]); -se.nbins = 60; +se.nbins = nbins; se.bins = linspace(se.xrange(1),se.xrange(2),se.nbins); [se.lesion.hist, se.lesion.xhist] = hist(se.lesion.rmse.all, se.bins); [se.nolesion.hist, se.nolesion.xhist] = hist(se.nolesion.rmse.all,se.bins); diff --git a/life/compute/feComputeEvidence_new.m b/life/compute/feComputeEvidence_new.m new file mode 100644 index 0000000..a5df0cf --- /dev/null +++ b/life/compute/feComputeEvidence_new.m @@ -0,0 +1,277 @@ +function se = feComputeEvidence_new(rmse1, rmse2, nbins, nmc, nboots) +% Computes a series of distance metrics between two RMSE distribtions. +% +% Compute summary statistics on to characterize the lesion,: +% We compute: +% - The strength of evidence, namely the effect size of the lesion +% - The Kullback-Leibler Divergence +% - Jeffrey's Divergence +% - The Earth Mover's distance +% +% Copyright (2020), Brent McPherson & Franco Pestilli, Indiana University +% + +% fill in default parameters + +if(~exist('nbins', 'var') || isempty(nbins)) + nbins = 128; +end + +if(~exist('nmc', 'var') || isempty(nmc)) + nmc = 5; +end + +if(~exist('nboots', 'var') || isempty(nboots)) + nboots = 250; % this almost certainly has to be higher, low for test +end + +% pull the size of a distribution +nobs1 = size(rmse1, 2); % nolesion +nobs2 = size(rmse2, 2); % lesion + +% build the combined null distribution +nrmse = [ rmse1, rmse2 ]; + +% the min/max observed values for SOE / d-prime histogram computation +minx = floor(mean(rmse1) - mean(rmse1)*.05); +maxx = ceil(mean(rmse2) + mean(rmse2)*.05); + +% prepare the distribution of errors and the histograms describing the errors. +se.nolesion.rmse.all = rmse1; +se.nolesion.rmse.mean = mean(se.nolesion.rmse.all); + +se.lesion.rmse.all = rmse2; +se.lesion.rmse.mean = mean(se.lesion.rmse.all); + +% compute histograms + +% define the bins based on the inputs / argument +se.nbins = nbins; +se.xrange = minmax([ se.lesion.rmse.all se.nolesion.rmse.all ]); +se.bins = linspace(se.xrange(1), se.xrange(2), se.nbins); + +% store params for reference +se.params.nmontecarlo = nmc; +se.params.nboots = nboots; + +% compute the histograms on the data +% the se.bins center have to be wrapped with inf to actually match hist +% https://www.mathworks.com/matlabcentral/answers/377930-histcounts-error-in-place-of-histc +[ se.lesion.hist, se.lesion.xhist ] = histcounts(se.lesion.rmse.all, [ se.bins inf ], 'Normalization', 'probability'); +[ se.nolesion.hist, se.nolesion.xhist ] = histcounts(se.nolesion.rmse.all, [ se.bins inf ], 'Normalization', 'probability'); + +% drop inf from end to keep bin centers the same size, b/c this is "fixed" +se.lesion.xhist = se.lesion.xhist(1:end-1); +se.nolesion.xhist = se.nolesion.xhist(1:end-1); + +%% compute the true differences between the rmse / histograms + +% dissimilarity between raw RMSE values +se.dis.name = 'Dissimilarity: https://nikokriegeskorte.org/category/representational-similarity-analysis/'; +se.dis.mean = pdist2(rmse1, rmse2, 'correlation'); + +% cosine dissimilarity between raw RMSE values +se.cos.name = 'Cosine Similarity: https://en.wikipedia.org/wiki/Cosine_similarity'; +se.cos.mean = pdist2(rmse1, rmse2, 'cosine'); + +% Kullback-Leibler Divergence +se.kl.name = sprintf('Kullback–Leibler divergence: http://en.wikipedia.org/wiki/Kullback-Leibler_divergence'); +tkl = se.nolesion.hist .* log2( (se.nolesion.hist) ./ (se.lesion.hist + eps) ); +se.kl.mean = nansum(tkl); +clear tkl + +% Jeffrey's divergence +se.j.name = sprintf('Jeffrey''s divergence: http://en.wikipedia.org/wiki/Divergence_(statistics)'); +tjd = se.nolesion.hist .* log2( (se.nolesion.hist) ./ ((se.lesion.hist + se.lesion.hist + eps)./2) ) + ... + se.lesion.hist .* log2( (se.lesion.hist) ./ ((se.lesion.hist + se.lesion.hist + eps)./2) ); +se.j.mean = nansum(tjd); clear tjd + +% Earth Mover's Distance: +% Note: This can be very slow and my require large amounts of memory for more than 1000 voxels +se.em.name = sprintf('Earth Mover''s distance: http://en.wikipedia.org/wiki/Earth_mover''s_distance'); + +try % Using Rubinov c-code fastest + + % estimate the distance between the bin centers + eDist = pdist2(se.nolesion.xhist', se.lesion.xhist', @(xi, xj) abs(xi-xj)); + + % compute the emd + tmp_em = emd_mex(se.nolesion.hist, se.lesion.hist, eDist); + +catch % else + + % compute EMD using slower MATLAB fxn + [ ~, tmp_em ] = emd(se.nolesion.xhist',se.lesion.xhist',se.nolesion.hist',se.lesion.hist',@gdf); + +end +se.em.mean = tmp_em; +clear tmp_emp + +%% add the bootstrapped null / se measures for all metrics + +% define null histograms for SOE estimate +mhist1 = nan(nbins, nmc); +mhist2 = nan(nbins, nmc); + +% preallocate test vectors for every measure (SE) +tcor = nan(nboots, nmc); +tcos = nan(nboots, nmc); +tkld = nan(nboots, nmc); +tjfd = nan(nboots, nmc); +temd = nan(nboots, nmc); +tdp1 = nan(nboots, nmc); +tdp2 = nan(nboots, nmc); + +% preallocate null vectors for every measure (pval) +ncor = nan(nboots, nmc); +ncos = nan(nboots, nmc); +nkld = nan(nboots, nmc); +njfd = nan(nboots, nmc); +nemd = nan(nboots, nmc); +ndpd = nan(nboots, nmc); + +for mc = 1:nmc + fprintf('.') + for boot = 1:nboots + + % pull random samples from each distribution for SE estimates + trmse1 = randsample(rmse1, nobs1, true); % nolesion + trmse2 = randsample(rmse2, nobs2, true); % lesion + + % pull random samples from combined distributions for p-value estimates + nrmse1 = randsample(nrmse, nobs1, true); % nolesion + nrmse2 = randsample(nrmse, nobs2, true); % lesion + + % build the boostrapped histograms for test and null + + % compute the test histograms w/ the estimated bins + [ t1hist, t1xhist ] = histcounts(trmse1, [ se.bins inf ], 'Normalization', 'probability'); + [ t2hist, t2xhist ] = histcounts(trmse2, [ se.bins inf ], 'Normalization', 'probability'); + + % fix the "fixed" bins + t1xhist = t1xhist(1:end-1); + t2xhist = t2xhist(1:end-1); + + % build the the combined null histogram + + % compute the null histograms w/ the estimated bins + [ n1hist, n1xhist ] = histcounts(nrmse1, [ se.bins inf ], 'Normalization', 'probability'); + [ n2hist, n2xhist ] = histcounts(nrmse2, [ se.bins inf ], 'Normalization', 'probability'); + + % fix the "fixed" bins + n1xhist = n1xhist(1:end-1); + n2xhist = n2xhist(1:end-1); + + % estimate each SE / null throw on raw RMSE values + + % SOE / d-prime test averages + tdp1(boot, mc) = mean(trmse1); + tdp2(boot, mc) = mean(trmse2); + + % SOE / d-prime null difference + ndp1 = mean(nrmse1); + ndp2 = mean(nrmse2); + ndpd(boot, mc) = ndp1 - ndp2; + + % estimate the test / null dissimilarity between rmse + tcor(boot, mc) = pdist2(trmse1, trmse2, 'correlation'); + ncor(boot, mc) = pdist2(nrmse1, nrmse2, 'correlation'); + + % estimate the test / null cosine distance between rmse + tcos(boot, mc) = pdist2(trmse1, trmse2, 'cosine'); + ncos(boot, mc) = pdist2(nrmse1, nrmse2, 'cosine'); + + % build the SE / null estimates on the resampled / null histograms + + % kl divergence + tkld(boot, mc) = nansum(t1hist .* log2( (t1hist) ./ (t2hist + eps) )); + nkld(boot, mc) = nansum(n1hist .* log2( (n1hist) ./ (n2hist + eps) )); + + % jeffery's divergence + tjfd(boot, mc) = nansum(t1hist .* log2( (t1hist) ./ ((t2hist + t2hist + eps)./2) ) + ... + t2hist .* log2( (t2hist) ./ ((t2hist + t2hist + eps)./2) )); + njfd(boot, mc) = nansum(n1hist .* log2( (n1hist) ./ ((n2hist + n2hist + eps)./2) ) + ... + n2hist .* log2( (n2hist) ./ ((n2hist + n2hist + eps)./2) )); + + % Earth Mover's Distance; cased to failover to matlab fxb if MEX fxn fails + try + + % estimate EMD on resampled / null + temd(boot, mc) = emd_mex(t1hist, t2hist, eDist); + nemd(boot, mc) = emd_mex(n1hist, n2hist, eDist); + + % pass to slower full matlab fxn if MEX fails + catch + + % slower MATLAB fxn + [ ~, temd(boot, mc) ] = emd(t1xhist', t2xhist', t1hist', t2hist', @gdf); + [ ~, nemd(boot, mc) ] = emd(n1xhist', n2xhist', n1hist', n2hist', @gdf); + + end + + end + + % pull the estimated bins for SOE histogram + mbins = linspace(minx, maxx, nbins); + + % compute and normalize histogram of unlesioned rmse + [ mhist1(:, mc), mxhist1 ] = histcounts(tdp1(:, mc), [ mbins inf ], 'Normalization', 'probability'); + [ mhist2(:, mc), mxhist2 ] = histcounts(tdp2(:, mc), [ mbins inf ], 'Normalization', 'probability'); + + % fix the "fixed" bins + mxhist1 = mxhist1(1:end-1); + mxhist2 = mxhist2(1:end-1); + +end +disp(' done.') + +%% store the permutation resluts in output struct + +% permuted dissimilarity measure +se.dis.se = mean(std(tcor, [], 1)); +se.dis.pval = min(sum(ncor > se.dis.mean) / nboots); + +% permuted cosine distance measure +se.cos.se = mean(std(tcos, [], 1)); +se.cos.pval = min(sum(ncos > se.cos.mean) / nboots); + +% permuted kl values +se.kl.se = mean(std(tkld, [], 1)); +se.kl.pval = min(sum(nkld > se.kl.mean) / nboots); + +% permuted jd values +se.j.se = mean(std(tjfd, [], 1)); +se.j.pval = min(sum(njfd > se.j.mean) / nboots); + +% permuted emd values +se.em.se = mean(std(temd, [], 1)); +se.em.pval = min(sum(nemd > se.em.mean) / nboots); + +% pull the tested SOE mean / std across MCs +smn = diff([ mean(tdp1, 1); mean(tdp2, 1) ]); +sse = sqrt(sum([ std(tdp1, [], 1); std(tdp2, [], 1) ].^2, 1)); + +% collapse SOE across MC runs +se.s.mean = mean(smn./sqrt(sse).^2); +se.s.std = std(smn./sqrt(sse).^2); +se.s.pval = min(sum(ndpd > smn) ./ nboots); + +% build mean histogram of SOE distributions +mn_mhist1 = mean(mhist1, 2); +mn_mhist2 = mean(mhist2, 2); + +% add 2 std around mean histogram estimates (confidence itervals) +ci_mhist1 = [mn_mhist1, mn_mhist1] + 2*[-std(mhist1, [], 2), std(mhist1, [], 2)]; +ci_mhist2 = [mn_mhist2, mn_mhist2] + 2*[-std(mhist2, [], 2), std(mhist2, [], 2)]; + +% store lesioned / unlesioned histograms +se.s.min_x = minx; +se.s.max_x = maxx; +se.s.unlesioned_mn = mn_mhist1; +se.s.unlesioned_ci = ci_mhist1; +se.s.unlesioned_xbins = mxhist1'; +se.s.lesioned_mn = mn_mhist2; +se.s.lesioned_ci = ci_mhist2; +se.s.lesioned_xbins = mxhist2'; + +end diff --git a/life/compute/feComputeEvidence_norm.m b/life/compute/feComputeEvidence_norm.m index 5398697..2f51eca 100644 --- a/life/compute/feComputeEvidence_norm.m +++ b/life/compute/feComputeEvidence_norm.m @@ -1,4 +1,4 @@ -function se = feComputeEvidence_norm(rmse1,rmse2) +function se = feComputeEvidence_norm(rmse1, rmse2, nbins) % Computes a series of distance metrics between two RMSE distribtions. % % Compute summary statistics on to characterize the lesion,: @@ -9,24 +9,33 @@ % - The Earth Mover's distance % % Copyright (2016), Franco Pestilli, Indiana University, frakkopesto@gmail.com. +% -% Prepare the distribution of errors and the histograms describing the -% erros. +% Prepare the distribution of errors and the histograms describing the errors. se.nolesion.rmse.all = rmse1; -se.lesion.rmse.all = rmse2; se.nolesion.rmse.mean = mean(se.nolesion.rmse.all); + +se.lesion.rmse.all = rmse2; se.lesion.rmse.mean = mean(se.lesion.rmse.all); % Histograms + +% define the bins based on the inputs / argument +se.nbins = nbins; se.xrange(1) = min([se.lesion.rmse.all se.nolesion.rmse.all]); se.xrange(2) = max([se.lesion.rmse.all se.nolesion.rmse.all]); -se.nbins = 60; se.bins = linspace(se.xrange(1),se.xrange(2),se.nbins); -[se.lesion.hist, se.lesion.xhist] = hist(se.lesion.rmse.all, se.bins); -[se.nolesion.hist, se.nolesion.xhist] = hist(se.nolesion.rmse.all,se.bins); -se.lesion.hist = se.lesion.hist ./ sum(se.lesion.hist); + +% compute the histograms on the data +[se.lesion.hist, se.lesion.xhist] = hist(se.lesion.rmse.all, se.bins); +[se.nolesion.hist, se.nolesion.xhist] = hist(se.nolesion.rmse.all, se.bins); + +% normalize the histograms +se.lesion.hist = se.lesion.hist ./ sum(se.lesion.hist); se.nolesion.hist = se.nolesion.hist ./ sum(se.nolesion.hist); +%% compute the true differences between the histograms + % Kullback-Leibler Divergence se.kl.name = sprintf('Kullback–Leibler divergence: http://en.wikipedia.org/wiki/Kullback-Leibler_divergence'); tmp = se.nolesion.hist .* log2( (se.nolesion.hist) ./ (se.lesion.hist + eps) ); @@ -62,7 +71,7 @@ clear tmp_emp % Strenght of evidence (effect size) -fprintf('[%s] Computing the Strength of Evidence... ',mfilename) +fprintf('[%s] Computing the Strength of Evidence... ', mfilename) se.s.name = sprintf('strength of evidence, d-prime: http://en.wikipedia.org/wiki/Effect_size'); se.s.nboots = 5000; se.s.nmontecarlo = 5; @@ -76,10 +85,33 @@ for inm = 1:se.s.nmontecarlo fprintf('.') for ibt = 1:se.s.nboots - nullDistributionW(ibt,inm) = mean(randsample(se.nolesion.rmse.all, sizeunlesioned,true)); - nullDistributionWO(ibt,inm) = mean(randsample(se.lesion.rmse.all,sizeunlesioned,true)); + nullDistributionW(ibt, inm) = mean(randsample(se.nolesion.rmse.all, sizeunlesioned, true)); + nullDistributionWO(ibt, inm) = mean(randsample(se.lesion.rmse.all, sizeunlesioned, true)); + end + + % build the combined null distribution + ndist = [ se.lesion.rmse.all, se.nolesion.rmse.all ]; + + % build the true difference to compare after null + empDiff = se.lesion.rmse.mean - se.nolesion.rmse.mean; + + % for every boot strap + for ibt = 1:se.s.nboots + + % pull a random sample from combined distributions + tnullDistributionW = mean(randsample(ndist, sizeunlesioned, true)); + tnullDistributionWO = mean(randsample(ndist, sizeunlesioned, true)); + + % pull the null difference for comparison + diffDist(ibt) = tnullDistributionW - tnullDistributionWO; + end + % get the count of nulls that are greater than observed + se.s.pval(inm) = sum(diffDist > empDiff) / se.s.nboots; + % find the code in CCA and copy it here to estimate this. + %pval(ii) = sum(grotRr(1, ii, :) > cca.cca.hocorrs(ii)) / Nperm; + % Distribution unlesioned [y(:,inm),xhis] = hist(nullDistributionW(:,inm),linspace(min_x,max_x,se.s.nbins)); y(:,inm) = y(:,inm)./sum(y(:,inm)); @@ -87,8 +119,9 @@ % Distribution lesioned [woy(:,inm),woxhis] = hist(nullDistributionWO(:,inm),linspace(min_x,max_x,se.s.nbins)); woy(:,inm) = woy(:,inm)./sum(woy(:,inm)); -end +end + se.s.mean = mean(diff([mean(nullDistributionW,1); ... mean(nullDistributionWO,1)])./sqrt(sum([std(nullDistributionW,[],1);std(nullDistributionWO,[],1)].^2,1))); se.s.std = std(diff([mean(nullDistributionW,1); ... @@ -109,4 +142,66 @@ se.s.min_x = min_x; se.s.max_x = max_x; +%% add the bootstrapped null / std measures for all metrics + +% pull the size of a distribution +nobs1 = size(rmse1, 2); +nobs2 = size(rmse2, 2); + +% build the combined null distribution +nrmse = [ rmse1, rmse2 ]; + +% preallocate null vector for every measure +ncor = nan(se.s.nboots, 1); +ncos = nan(se.s.nboots, 1); +nkld = nan(se.s.nboots, 1); + +for ibt = 1:se.s.nboots + + % pull random samples from each distribution + trmse1 = randsample(rmse1, nobs1, true); % nolesion + trmse2 = randsample(rmse2, nobs2, true); % lesion + + % pull random samples from combined distributions + nrmse1 = randsample(nrmse, nobs1, true); + nrmse2 = randsample(nrmse, nobs2, true); + + % estimate each null throw on raw RMSE values + ncor(ibt) = pdist2(nrmse1, nrmse2, 'correlation'); + ncos(ibt) = pdist2(nrmse1, nrmse2, 'cosine'); + + % build the null histogram + + % define the bins based on the inputs / argument + tse.xrange = minmax([ trmse1 trmse2 ]); + tse.bins = linspace(tse.xrange(1), tse.xrange(2), nbins); + + % compute the histograms on the data + [ tse.lesion.hist, tse.lesion.xhist ] = hist(trmse2, nbins); + [ tse.nolesion.hist, tse.nolesion.xhist ] = hist(trmse1, nbins); + + % normalize the histograms + tse.lesion.hist = tse.lesion.hist ./ sum(tse.lesion.hist); + tse.nolesion.hist = tse.nolesion.hist ./ sum(tse.nolesion.hist); + + % kl divergence estimated on null + tkld = tse.nolesion.hist .* log2( (tse.nolesion.hist) ./ (tse.lesion.hist + eps) ); + nkld(ibt) = nansum(tkld); + +end + +% dissimilarity between raw RMSE values +se.dis.name = 'Dissimilarity'; +se.dis.mean = pdist2(rmse1, rmse2, 'correlation'); +se.dis.pval = 1 - (sum(ncor > se.dis.mean) / se.s.nboots); + +% cosine dissimilarity between raw RMSE values +se.cos.name = 'Cosine Distance'; +se.cos.mean = pdist2(rmse1, rmse2, 'cosine'); +se.cos.pval = 1 - (sum(ncos > se.cos.mean) / se.s.nboots); + +se.nkl.name = "New KL Divergence"; +se.nkl.mean = se.kl.mean; +se.nkl.pval = 1 - (sum(nkld > se.kl.mean) / se.s.nboots); + end From 0d0f89da927f7b0d3e7ae8c74843e155112519c1 Mon Sep 17 00:00:00 2001 From: bcmcpher Date: Fri, 13 Nov 2020 21:07:15 -0500 Subject: [PATCH 04/19] modified to look for modelKurtosis in fe and estimate the apparent kurtosis coefficient (akc); modifications for multishell data - bookmark, not for release --- life/fe/feBuildDictionaries.m | 126 +++++++++++++++++++--------------- 1 file changed, 70 insertions(+), 56 deletions(-) diff --git a/life/fe/feBuildDictionaries.m b/life/fe/feBuildDictionaries.m index 699ac65..4770170 100755 --- a/life/fe/feBuildDictionaries.m +++ b/life/fe/feBuildDictionaries.m @@ -39,20 +39,75 @@ bvecs = feGet(fe,'bvecs'); % bvecs bvals = feGet(fe,'bvals'); % bvals +%bvals = bvals * 1000; + Dict = zeros(nBvecs,Norient); % Initialize Signal Dictionary matrix DictSig = zeros(nBvecs,Norient); % Initialize Signal Dictionary matrix %DictTensors = zeros(9,Norient); % Initialize Tensors Dictionary matrix % catch Kurtosis estimates for debugging -akc = zeros(nBvecs,Norient); -KDict = zeros(nBvecs,Norient); - -%[ shells, ~, shellindex ] = unique(round(bvals)); % round to get around noise in shells - FRAGILE FIX +%KDict = zeros(nBvecs,Norient); +akc = zeros(nBvecs,1); % pull the shell information ubv = feGet(fe, 'nshells'); ubi = feGet(fe, 'shellindex'); +% pull diffusion tensor +dt = feGet(fe, 'model tensor'); + +% pull kurtosis fit - needs feGet +kt = fe.life.modelKurtosis; + +% build apparent kurtosis coefficient - akc +for i=1:size(ubv,1) + + % pull tensor parameters for hard indexing of equations + sdt = dt(i,:); + + % find the indices for the shell + si = ubi == ubv(i); + + % mean diffusivity of tensor from forward model + md = mean(sdt); + + % apparent diffusion coefficient - matches dipy + adc = bvecs(si,1) .* bvecs(si,1) * sdt(1) + ... + 2 * bvecs(si,1) .* bvecs(si,2) * sdt(2) + ... + bvecs(si,2) .* bvecs(si,2) .* sdt(3) + ... + 2 * bvecs(si,1) .* bvecs(si,3) * 0 + ... % only store primary eigenvalues in + 2 * bvecs(si,2) .* bvecs(si,3) * 0 + ... % LiFE prediction. These are 0. + bvecs(si,3) .* bvecs(si,3) * 0; + + % apparent diffusion variance - matches dipy + adv = ... + bvecs(si,1) .* bvecs(si,1) .* bvecs(si,1) .* bvecs(si,1) .* kt(1) + ... % xxxx + bvecs(si,2) .* bvecs(si,2) .* bvecs(si,2) .* bvecs(si,2) .* kt(2) + ... % yyyy + bvecs(si,3) .* bvecs(si,3) .* bvecs(si,3) .* bvecs(si,3) .* kt(3) + ... % zzzz + 4 * bvecs(si,1) .* bvecs(si,1) .* bvecs(si,1) .* bvecs(si,2) .* kt(4) + ... % xxxy + 4 * bvecs(si,1) .* bvecs(si,1) .* bvecs(si,1) .* bvecs(si,3) .* kt(5) + ... % xxxz + 4 * bvecs(si,1) .* bvecs(si,2) .* bvecs(si,2) .* bvecs(si,2) .* kt(6) + ... % xyyy + 4 * bvecs(si,1) .* bvecs(si,3) .* bvecs(si,3) .* bvecs(si,3) .* kt(7) + ... % xzzz + 4 * bvecs(si,2) .* bvecs(si,2) .* bvecs(si,2) .* bvecs(si,3) .* kt(8) + ... % yyyz + 4 * bvecs(si,2) .* bvecs(si,3) .* bvecs(si,3) .* bvecs(si,3) .* kt(9) + ... % yzzz + 6 * bvecs(si,1) .* bvecs(si,1) .* bvecs(si,2) .* bvecs(si,2) .* kt(10) + ... % xxyy + 6 * bvecs(si,1) .* bvecs(si,1) .* bvecs(si,3) .* bvecs(si,3) .* kt(11) + ... % xxzz + 6 * bvecs(si,2) .* bvecs(si,2) .* bvecs(si,3) .* bvecs(si,3) .* kt(12) + ... % yyzz + 12 * bvecs(si,1) .* bvecs(si,1) .* bvecs(si,2) .* bvecs(si,3) .* kt(13) + ... % xxyz + 12 * bvecs(si,1) .* bvecs(si,2) .* bvecs(si,2) .* bvecs(si,3) .* kt(14) + ... % xyyz + 12 * bvecs(si,1) .* bvecs(si,2) .* bvecs(si,3) .* bvecs(si,3) .* kt(15); % xyzz + + % zero out noisy (bad) adc estimates + adc(adc < 0) = 0; + + % estimate apparent kutosis coefficient - matches dipy + akc(si) = adv .* ((md ./ adc).^2); + + % zero out noisy (bad) akc estimates + akc(akc < -3/7) = -3/7; + +end + % Compute each dictionary column for a different kernel orientation for j=1:Norient @@ -60,7 +115,7 @@ [Rot,~, ~] = svd(orient(:,j)); % for every shell - for k=1:size(ubv, 1) + for k=1:size(ubv,1) % find the indices for the shell si = ubi == ubv(k); @@ -72,65 +127,24 @@ % estimate Q for tensor values in shell Q = Rot*D*Rot'; - % logical starts here: if kurtosis run these steps, else just fill in tensor - - % mean diffusivity of tensor? - md = mean(diag(D)); - - % pull tensor / kurtosis parameters for hard indexing of equations - dt = fe.life.modelTensor(k,:); - kt = fe.life.modelKurtosis; - - % apparent diffusion coefficient from dipy - - adc = bvecs(si,1) .* bvecs(si,1) * dt(1) + ... - 2 * bvecs(si,1) .* bvecs(si,2) * dt(2) + ... - bvecs(si,2) .* bvecs(si,2) .* dt(3) + ... - 2 * bvecs(si,1) .* bvecs(si,3) * 0 + ... - 2 * bvecs(si,2) .* bvecs(si,3) * 0 + ... - bvecs(si,3) .* bvecs(si,3) * 0; - - % apparent diffusion variance from dipy - % reordered to ensure match of kt params from mrtrix3 - adv = ... - bvecs(si,1) .* bvecs(si,1) .* bvecs(si,1) .* bvecs(si,1) .* kt(1) + ... % xxxx - bvecs(si,2) .* bvecs(si,2) .* bvecs(si,2) .* bvecs(si,2) .* kt(2) + ... % yyyy - bvecs(si,3) .* bvecs(si,3) .* bvecs(si,3) .* bvecs(si,3) .* kt(3) + ... % zzzz - 4 * bvecs(si,1) .* bvecs(si,1) .* bvecs(si,1) .* bvecs(si,2) .* kt(4) + ... % xxxy - 4 * bvecs(si,1) .* bvecs(si,1) .* bvecs(si,1) .* bvecs(si,3) .* kt(5) + ... % xxxz - 4 * bvecs(si,1) .* bvecs(si,2) .* bvecs(si,2) .* bvecs(si,2) .* kt(6) + ... % xyyy - 4 * bvecs(si,1) .* bvecs(si,3) .* bvecs(si,3) .* bvecs(si,3) .* kt(7) + ... % xzzz - 4 * bvecs(si,2) .* bvecs(si,2) .* bvecs(si,2) .* bvecs(si,3) .* kt(8) + ... % yyyz - 4 * bvecs(si,2) .* bvecs(si,3) .* bvecs(si,3) .* bvecs(si,3) .* kt(9) + ... % yzzz - 6 * bvecs(si,1) .* bvecs(si,1) .* bvecs(si,2) .* bvecs(si,2) .* kt(10) + ... % xxyy - 6 * bvecs(si,1) .* bvecs(si,1) .* bvecs(si,3) .* bvecs(si,3) .* kt(11) + ... % xxzz - 6 * bvecs(si,2) .* bvecs(si,2) .* bvecs(si,3) .* bvecs(si,3) .* kt(12) + ... % yyzz - 12 * bvecs(si,1) .* bvecs(si,1) .* bvecs(si,2) .* bvecs(si,3) .* kt(13) + ... % xxyz - 12 * bvecs(si,1) .* bvecs(si,2) .* bvecs(si,2) .* bvecs(si,3) .* kt(14) + ... % xyyz - 12 * bvecs(si,1) .* bvecs(si,2) .* bvecs(si,3) .* bvecs(si,3) .* kt(15); % xyzz - - % zero out noisy (bad) adc estimates - adc(adc < 0) = 0; - - % estimate apparent kutosis coefficient - store in temp var - takc = adv .* ((md / adc).^2)'; % always finds same param at 20/50, rest are zeros? - - % zero out noisy (bad) akc estimates - use tmp to avoid indexing of large matrix - takc(takc < -3/7) = -3/7; - - akc(si,j) = takc; - % Compute the signal contribution of a fiber in the kernel orientation divided S0 Dict(si,j) = exp(-bvals(si) .* diag(bvecs(si,:)*Q*bvecs(si,:)')); - KDict(si,j) = exp(-bvals(si) .* diag(bvecs(si,:)*Q*bvecs(si,:)') + (-bvals(si).^2 .* diag(bvecs(si,:)*Q*bvecs(si,:)').^2 .* akc(si,j))/6); + %Dict(si,j) = exp(-bvals(si) .* diag(bvecs(si,:)*Q*bvecs(si,:)') + (-bvals(si).^2 .* diag(bvecs(si,:)*Q*bvecs(si,:)').^2 .* akc(si))/6); - % demedianed signal by shell (used to demean) - DictSig(si,j) = Dict(si,j) - median(Dict(si,j)); + % demeaned signal by shell + DictSig(si,j) = Dict(si,j) - mean(Dict(si,j)); end end -disp('Done.'); +% figure; hold on; +% for ii = 1:5000 +% plot(Dict(:, ii), KDict(:, ii), '.', 'color', 'k'); +% end +% title('Predicted Diffusion x Predicted Kurtosis'); +% xlabel('Diffusion'); +% ylabel('Kurtosis'); fe = feSet(fe,'dictionary parameters',{Nphi,Ntheta,orient,Dict,DictSig}); From 310cae53263c06f72ed6baf464d074fe546185d1 Mon Sep 17 00:00:00 2001 From: bcmcpher Date: Fri, 13 Nov 2020 21:09:06 -0500 Subject: [PATCH 05/19] changes for passing multishell / kurtosis parameters, fixed .gitignore for *.m~ files --- .gitignore | 2 +- life/fe/feConnectomeEncoding.m | 3 --- life/fe/feConnectomeSetDwi.m | 4 ++-- life/fe/feGet.m | 26 ++++++++++++++++++-------- 4 files changed, 21 insertions(+), 14 deletions(-) diff --git a/.gitignore b/.gitignore index 3084eda..6d725d9 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1 @@ -.m~ \ No newline at end of file +*.m~ diff --git a/life/fe/feConnectomeEncoding.m b/life/fe/feConnectomeEncoding.m index e0d3463..f4e3f38 100755 --- a/life/fe/feConnectomeEncoding.m +++ b/life/fe/feConnectomeEncoding.m @@ -144,9 +144,6 @@ roi_coords = feGet(fe,'roicoords'); roi_ind = sub2ind(imgsize,roi_coords(:,1)',roi_coords(:,2)',roi_coords(:,3)'); -% SIGNAL NORMALIZATION NEEDS TO HAPPEN BY SHELL IN THIS OPERATION -% IT DOES NOT CURRENTLY - % Restrict the Sparse Tensor and matrix A to the ROI voxels only Phi = Phi(:,roi_ind,:); % reduce tensor in 3rd dimension to roi voxels only A = A(roi_ind,:); diff --git a/life/fe/feConnectomeSetDwi.m b/life/fe/feConnectomeSetDwi.m index c876175..ff40399 100755 --- a/life/fe/feConnectomeSetDwi.m +++ b/life/fe/feConnectomeSetDwi.m @@ -48,13 +48,13 @@ dwiGet(dwi, 'diffusionimagenums')); % Handling multishell data. -% The original model was only able to handl single-shell. +% The original model was only able to handle single-shell. % This new version will handle also multishell data. % % To handle multishell data we will store the number of shells encoded % for the data (1,2,3, etc) and the indices to the bvecs for each shell % (this will NOT assume that the number of BVECS is equal for each shell -% (say 30 diffusion directiosn for each shell). +% (say 30 diffusion directions for each shell). % We extract the bvalues and find the unique of them and assign them to each shell. bvals = feGet(fe,'bvals'); diff --git a/life/fe/feGet.m b/life/fe/feGet.m index 8703a78..615e79a 100755 --- a/life/fe/feGet.m +++ b/life/fe/feGet.m @@ -566,7 +566,7 @@ voxelIndices = feGet(fe,'voxelsindices',varargin); val = fe.life.diffusion_signal_img(voxelIndices,:) - repmat(mean(fe.life.diffusion_signal_img(voxelIndices,:), 2),1,nBvecs); keyboard - % THis seems to be wrong + % THis seems to be wrong - DOES NOT DO MULTISHELL CORRECTLY case {'b0signalimage','b0vox'} % Get the diffusion signal at 0 diffusion weighting (B0) for this voxel @@ -929,11 +929,22 @@ % dSig = feGet(fe,'dsigdemeaned',coords); nVoxels = feGet(fe,'nVoxels'); nBvecs = feGet(fe,'nBvecs'); + shells = feGet(fe,'shellindex'); + nshell = feGet(fe,'nshells'); + + % THIS IS DEMEANING ACROSS SHELL + dSig = fe.life.diffusion_signal_img; + %dSig = reshape(fe.life.diffusion_signal_img',[1,nVoxels*nBvecs]); + + % loop over shells for demeaing + for shell = 1:size(nshell, 1) + s = nshell(shell); % index into listed shell for good logic + dSigMean = mean(dSig(:, shells == s), 2); + dSig(:, shells == s) = dSig(:, shells == s) - dSigMean; + end + %val = (dSig - reshape(repmat(mean(reshape(dSig, nBvecs, nVoxels),1),nBvecs,1), size(dSig)))'; + val = reshape(dSig, [nVoxels*nBvecs,1]); - dSig = reshape(fe.life.diffusion_signal_img',[1,nVoxels*nBvecs]); - val = (dSig - reshape(repmat( ... - mean(reshape(dSig, nBvecs, nVoxels),1),... - nBvecs,1), size(dSig)))'; % Return a subset of voxels if ~isempty(varargin) % voxelIndices = feGet(fe,'voxelsindices',varargin); @@ -1140,7 +1151,7 @@ % res = feGet(fe,'res sig full'); % res = feGet(fe, 'res sig full',coords); % res = feGet(fe, 'res sig full',voxelIndex); - val = (feGet(fe,'dsig full') - feGet(fe,'psig full')'); + val = (feGet(fe,'dsig full') - feGet(fe,'psig full')'); if ~isempty(varargin) % voxelIndices = feGet(fe,'voxelsindices',varargin); % voxelRowsToKeep = feGet(fe,'voxel rows',voxelIndices); @@ -1314,7 +1325,7 @@ % A volume of RMSE normalized by the S0 value in each voxel. rmse = feGet(fe,'vox rmse'); s0 = feGet(fe,'b0signalimage')'; - % Some voxels can have a S0=0. We replace the S0 in these voxles with + % Some voxels can have a S0=0. We replace the S0 in these voxels with % a NaN. idx = (s0 == 0); rmse(idx) = nan(size(find(idx))); @@ -1655,7 +1666,6 @@ otherwise help('feGet') fprintf('[feGet] Unknown parameter << %s >>...\n',param); - keyboard end end % END MAIN FUNCTION From b10fe31302546462d7f2bc6d5000e92264c361fe Mon Sep 17 00:00:00 2001 From: bcmcpher Date: Sat, 14 Nov 2020 11:44:56 -0500 Subject: [PATCH 06/19] adjusted refs to shell so demeaning actully happens --- life/fe/feConnectomeSetDwi.m | 2 +- life/fe/feGet.m | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/life/fe/feConnectomeSetDwi.m b/life/fe/feConnectomeSetDwi.m index ff40399..a4d50ce 100755 --- a/life/fe/feConnectomeSetDwi.m +++ b/life/fe/feConnectomeSetDwi.m @@ -59,7 +59,7 @@ % We extract the bvalues and find the unique of them and assign them to each shell. bvals = feGet(fe,'bvals'); [nshells, ~, shellindex] = unique(round(bvals)); -fe = feSet(fe, 'nshells', nshells); +fe = feSet(fe, 'nshells', length(nshells)); fe = feSet(fe, 'shellindex', shellindex); dim = dwi.nifti.dim; diff --git a/life/fe/feGet.m b/life/fe/feGet.m index 615e79a..e294a15 100755 --- a/life/fe/feGet.m +++ b/life/fe/feGet.m @@ -931,14 +931,14 @@ nBvecs = feGet(fe,'nBvecs'); shells = feGet(fe,'shellindex'); nshell = feGet(fe,'nshells'); - + ushell = unique(shells); % THIS IS DEMEANING ACROSS SHELL dSig = fe.life.diffusion_signal_img; %dSig = reshape(fe.life.diffusion_signal_img',[1,nVoxels*nBvecs]); % loop over shells for demeaing - for shell = 1:size(nshell, 1) - s = nshell(shell); % index into listed shell for good logic + for shell = 1:nshell + s = ushell(shell); % index into listed shell for good logic dSigMean = mean(dSig(:, shells == s), 2); dSig(:, shells == s) = dSig(:, shells == s) - dSigMean; end From 7555687260150036d14d9f5ed076f42de0a7f96f Mon Sep 17 00:00:00 2001 From: bcmcpher Date: Mon, 16 Nov 2020 17:44:37 -0500 Subject: [PATCH 07/19] add current version of feBuildDictionary that will fit kurtosis model to LiFE forward model --- life/fe/feBuildDictionaries.m | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/life/fe/feBuildDictionaries.m b/life/fe/feBuildDictionaries.m index 4770170..ced49d8 100755 --- a/life/fe/feBuildDictionaries.m +++ b/life/fe/feBuildDictionaries.m @@ -52,6 +52,7 @@ % pull the shell information ubv = feGet(fe, 'nshells'); ubi = feGet(fe, 'shellindex'); +ubl = unique(ubi); % pull diffusion tensor dt = feGet(fe, 'model tensor'); @@ -60,13 +61,13 @@ kt = fe.life.modelKurtosis; % build apparent kurtosis coefficient - akc -for i=1:size(ubv,1) +for i=1:ubv % pull tensor parameters for hard indexing of equations sdt = dt(i,:); % find the indices for the shell - si = ubi == ubv(i); + si = ubi == ubl(i); % mean diffusivity of tensor from forward model md = mean(sdt); @@ -115,10 +116,10 @@ [Rot,~, ~] = svd(orient(:,j)); % for every shell - for k=1:size(ubv,1) + for k=1:ubv % find the indices for the shell - si = ubi == ubv(k); + si = ubi == ubl(k); % create diagonal matix with diffusivities for current shell % this assumes tensor fits for shell are entered in the order this will parse them in @@ -128,8 +129,8 @@ Q = Rot*D*Rot'; % Compute the signal contribution of a fiber in the kernel orientation divided S0 - Dict(si,j) = exp(-bvals(si) .* diag(bvecs(si,:)*Q*bvecs(si,:)')); - %Dict(si,j) = exp(-bvals(si) .* diag(bvecs(si,:)*Q*bvecs(si,:)') + (-bvals(si).^2 .* diag(bvecs(si,:)*Q*bvecs(si,:)').^2 .* akc(si))/6); + %Dict(si,j) = exp(-bvals(si) .* diag(bvecs(si,:)*Q*bvecs(si,:)')); + Dict(si,j) = exp(-bvals(si) .* diag(bvecs(si,:)*Q*bvecs(si,:)') + (-bvals(si).^2 .* diag(bvecs(si,:)*Q*bvecs(si,:)').^2 .* akc(si))/6); % demeaned signal by shell DictSig(si,j) = Dict(si,j) - mean(Dict(si,j)); From 06fd3da618fd78e118ab20fa8bf30f1efe18e157 Mon Sep 17 00:00:00 2001 From: bcmcpher Date: Fri, 20 Nov 2020 16:45:22 -0500 Subject: [PATCH 08/19] build and write kurtosis dictionary in feBuildDictionaries, but use tensor model still - this is still locked on even when the kurtosis model isnt valid. Split one long print statement into 3 lines in feConnectomeEncoding. comments changed in feGet --- life/fe/feBuildDictionaries.m | 18 +++++++++++++----- life/fe/feConnectomeEncoding.m | 4 +++- life/fe/feGet.m | 3 ++- 3 files changed, 18 insertions(+), 7 deletions(-) diff --git a/life/fe/feBuildDictionaries.m b/life/fe/feBuildDictionaries.m index ced49d8..58dd940 100755 --- a/life/fe/feBuildDictionaries.m +++ b/life/fe/feBuildDictionaries.m @@ -39,6 +39,7 @@ bvecs = feGet(fe,'bvecs'); % bvecs bvals = feGet(fe,'bvals'); % bvals +%bvecs(:,1) = -bvecs(:,1); % better observed x predicted correlation %bvals = bvals * 1000; Dict = zeros(nBvecs,Norient); % Initialize Signal Dictionary matrix @@ -46,7 +47,8 @@ %DictTensors = zeros(9,Norient); % Initialize Tensors Dictionary matrix % catch Kurtosis estimates for debugging -%KDict = zeros(nBvecs,Norient); +KDict = zeros(nBvecs,Norient); +KDictSig = zeros(nBvecs,Norient); akc = zeros(nBvecs,1); % pull the shell information @@ -106,7 +108,7 @@ % zero out noisy (bad) akc estimates akc(akc < -3/7) = -3/7; - + end % Compute each dictionary column for a different kernel orientation @@ -123,17 +125,18 @@ % create diagonal matix with diffusivities for current shell % this assumes tensor fits for shell are entered in the order this will parse them in - D = diag(fe.life.modelTensor(k,:)); + D = diag(dt(k,:)); % estimate Q for tensor values in shell Q = Rot*D*Rot'; % Compute the signal contribution of a fiber in the kernel orientation divided S0 - %Dict(si,j) = exp(-bvals(si) .* diag(bvecs(si,:)*Q*bvecs(si,:)')); - Dict(si,j) = exp(-bvals(si) .* diag(bvecs(si,:)*Q*bvecs(si,:)') + (-bvals(si).^2 .* diag(bvecs(si,:)*Q*bvecs(si,:)').^2 .* akc(si))/6); + Dict(si,j) = exp(-bvals(si) .* diag(bvecs(si,:)*Q*bvecs(si,:)')); + KDict(si,j) = exp(-bvals(si) .* diag(bvecs(si,:)*Q*bvecs(si,:)') + (-bvals(si).^2 .* diag(bvecs(si,:)*Q*bvecs(si,:)').^2 .* akc(si))/6); % demeaned signal by shell DictSig(si,j) = Dict(si,j) - mean(Dict(si,j)); + KDictSig(si,j) = KDict(si,j) - mean(KDict(si,j)); end @@ -149,4 +152,9 @@ fe = feSet(fe,'dictionary parameters',{Nphi,Ntheta,orient,Dict,DictSig}); +% hard add kurtosis dictionaries / akc for debugging +fe.life.M.KDict = KDict; +fe.life.M.KDictSig = KDictSig; +fe.life.M.akc = akc; + end \ No newline at end of file diff --git a/life/fe/feConnectomeEncoding.m b/life/fe/feConnectomeEncoding.m index f4e3f38..e29bb66 100755 --- a/life/fe/feConnectomeEncoding.m +++ b/life/fe/feConnectomeEncoding.m @@ -62,7 +62,9 @@ nBatch = max([5,ceil(nTotalNodes/nNodesMax)]); % Number of batches nFib_Batch = ceil(nFibers/nBatch); % Number of fibers per batch -disp(['Max number of nodes per batch',num2str(nNodesMax),', Total number of nodes = ',num2str(nTotalNodes),'. Building the model in ',num2str(nBatch),' batches']) +disp(['Max number of nodes per batch: ',num2str(nNodesMax)]); +disp(['Total number of nodes = ',num2str(nTotalNodes)]); +disp(['Building the model in ',num2str(nBatch),' batches.']); fprintf('\n'); for n=1:nBatch diff --git a/life/fe/feGet.m b/life/fe/feGet.m index e294a15..ba8fd54 100755 --- a/life/fe/feGet.m +++ b/life/fe/feGet.m @@ -566,7 +566,8 @@ voxelIndices = feGet(fe,'voxelsindices',varargin); val = fe.life.diffusion_signal_img(voxelIndices,:) - repmat(mean(fe.life.diffusion_signal_img(voxelIndices,:), 2),1,nBvecs); keyboard - % THis seems to be wrong - DOES NOT DO MULTISHELL CORRECTLY + % THIS IS WRONG - DOES NOT DO MULTISHELL CORRECTLY + % BUT IS IT USED ANYWHERE? case {'b0signalimage','b0vox'} % Get the diffusion signal at 0 diffusion weighting (B0) for this voxel From 2ff8f772e28ae5baa7d9653146db62de3092a255 Mon Sep 17 00:00:00 2001 From: bcmcpher Date: Mon, 30 Nov 2020 18:22:59 -0500 Subject: [PATCH 09/19] fix feGet(fe, dsigdemeaned) by transposing input to match what was apparently the original orientation of the data for linear indexing --- life/fe/feGet.m | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/life/fe/feGet.m b/life/fe/feGet.m index ba8fd54..65726ca 100755 --- a/life/fe/feGet.m +++ b/life/fe/feGet.m @@ -933,19 +933,22 @@ shells = feGet(fe,'shellindex'); nshell = feGet(fe,'nshells'); ushell = unique(shells); - % THIS IS DEMEANING ACROSS SHELL - dSig = fe.life.diffusion_signal_img; - %dSig = reshape(fe.life.diffusion_signal_img',[1,nVoxels*nBvecs]); + + % HOPEFULLY THIS WORKS NOW - NEEDED TO BE TRANSPOSED + dSig = fe.life.diffusion_signal_img'; % loop over shells for demeaing for shell = 1:nshell s = ushell(shell); % index into listed shell for good logic - dSigMean = mean(dSig(:, shells == s), 2); - dSig(:, shells == s) = dSig(:, shells == s) - dSigMean; + dSigMean = mean(dSig(shells == s, :)); + dSig(shells == s, :) = dSig(shells == s, :) - repmat(dSigMean, [ sum(shells == s) 1 ]); end - %val = (dSig - reshape(repmat(mean(reshape(dSig, nBvecs, nVoxels),1),nBvecs,1), size(dSig)))'; val = reshape(dSig, [nVoxels*nBvecs,1]); + % % original demeaning operation + %dSig = reshape(fe.life.diffusion_signal_img',[1,nVoxels*nBvecs]); + %val = (dSig - reshape(repmat(mean(reshape(dSig, nBvecs, nVoxels),1),nBvecs,1), size(dSig)))'; + % Return a subset of voxels if ~isempty(varargin) % voxelIndices = feGet(fe,'voxelsindices',varargin); From 3a57b5e22c83d467ecd58403d2b4468c0b421e13 Mon Sep 17 00:00:00 2001 From: bcmcpher Date: Wed, 2 Dec 2020 15:05:49 -0500 Subject: [PATCH 10/19] fixed the problem of transposed imaging data during call to get demeaned data, updated comments to match / describe --- life/fe/feGet.m | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/life/fe/feGet.m b/life/fe/feGet.m index 65726ca..5b1fac3 100755 --- a/life/fe/feGet.m +++ b/life/fe/feGet.m @@ -934,7 +934,8 @@ nshell = feGet(fe,'nshells'); ushell = unique(shells); - % HOPEFULLY THIS WORKS NOW - NEEDED TO BE TRANSPOSED + % this works now - data needs to be transposed from how it's stored to + % match the unwound data used during prediction dSig = fe.life.diffusion_signal_img'; % loop over shells for demeaing From 5ca9042a1f7bd3eeb77249555c0629ceaaf1f898 Mon Sep 17 00:00:00 2001 From: bcmcpher Date: Wed, 2 Dec 2020 15:08:30 -0500 Subject: [PATCH 11/19] commented out all lines related to bvec flips and kurtosis fit. All the functions need to be updated to correctly parse inputs for these options - but that has to happen when the options work in the first place --- life/fe/feBuildDictionaries.m | 22 +++++++--------------- 1 file changed, 7 insertions(+), 15 deletions(-) diff --git a/life/fe/feBuildDictionaries.m b/life/fe/feBuildDictionaries.m index 58dd940..56d36d4 100755 --- a/life/fe/feBuildDictionaries.m +++ b/life/fe/feBuildDictionaries.m @@ -40,15 +40,15 @@ bvals = feGet(fe,'bvals'); % bvals %bvecs(:,1) = -bvecs(:,1); % better observed x predicted correlation -%bvals = bvals * 1000; +%bvals = bvals * 10; % does bad Dict = zeros(nBvecs,Norient); % Initialize Signal Dictionary matrix DictSig = zeros(nBvecs,Norient); % Initialize Signal Dictionary matrix %DictTensors = zeros(9,Norient); % Initialize Tensors Dictionary matrix % catch Kurtosis estimates for debugging -KDict = zeros(nBvecs,Norient); -KDictSig = zeros(nBvecs,Norient); +%KDict = zeros(nBvecs,Norient); +%KDictSig = zeros(nBvecs,Norient); akc = zeros(nBvecs,1); % pull the shell information @@ -132,29 +132,21 @@ % Compute the signal contribution of a fiber in the kernel orientation divided S0 Dict(si,j) = exp(-bvals(si) .* diag(bvecs(si,:)*Q*bvecs(si,:)')); - KDict(si,j) = exp(-bvals(si) .* diag(bvecs(si,:)*Q*bvecs(si,:)') + (-bvals(si).^2 .* diag(bvecs(si,:)*Q*bvecs(si,:)').^2 .* akc(si))/6); + %Dict(si,j) = exp(-bvals(si) .* diag(bvecs(si,:)*Q*bvecs(si,:)') + (-bvals(si).^2 .* diag(bvecs(si,:)*Q*bvecs(si,:)').^2 .* akc(si))/6); % demeaned signal by shell DictSig(si,j) = Dict(si,j) - mean(Dict(si,j)); - KDictSig(si,j) = KDict(si,j) - mean(KDict(si,j)); + %KDictSig(si,j) = KDict(si,j) - mean(KDict(si,j)); end end -% figure; hold on; -% for ii = 1:5000 -% plot(Dict(:, ii), KDict(:, ii), '.', 'color', 'k'); -% end -% title('Predicted Diffusion x Predicted Kurtosis'); -% xlabel('Diffusion'); -% ylabel('Kurtosis'); - fe = feSet(fe,'dictionary parameters',{Nphi,Ntheta,orient,Dict,DictSig}); % hard add kurtosis dictionaries / akc for debugging -fe.life.M.KDict = KDict; -fe.life.M.KDictSig = KDictSig; +%fe.life.M.KDict = KDict; +%fe.life.M.KDictSig = KDictSig; fe.life.M.akc = akc; end \ No newline at end of file From 5b6b18fd6e09248f2d5183bd788b2dce1299bfef Mon Sep 17 00:00:00 2001 From: bcmcpher Date: Fri, 11 Dec 2020 03:55:31 -0500 Subject: [PATCH 12/19] added feGet/feSet options for getting offset from input files, modified feConnectomeInit to pass warnings if a mismatch in file orientations is detected --- life/fe/feConnectomeInit.m | 37 ++++++++++++++++++++++++++++++++++++- life/fe/feGet.m | 5 +++++ life/fe/feSet.m | 5 +++++ 3 files changed, 46 insertions(+), 1 deletion(-) diff --git a/life/fe/feConnectomeInit.m b/life/fe/feConnectomeInit.m index 5941ec9..43b0e49 100755 --- a/life/fe/feConnectomeInit.m +++ b/life/fe/feConnectomeInit.m @@ -30,6 +30,9 @@ tempNi = niftiRead(dwiFile); fe = feSet(fe, 'img2acpc xform', tempNi.qto_xyz); fe = feSet(fe, 'acpc2img xform', inv(tempNi.qto_xyz)); + +% set the offset for files as well, then clear it +fe = feSet(fe, 'dwi offset', tempNi.qto_xyz(1:3,4)'); clear tempNi % Set up the fe name @@ -62,7 +65,39 @@ %% Anatomy Install the path to the anatomical high-resolution file. if ~notDefined('anatomyFile') - fe = feSet(fe,'anatomy file',anatomyFile); + + fe = feSet(fe,'anatomy file',anatomyFile); + + try + tempNi = niftiRead(anatomyFile); + catch + error('Anatomy file unable to be opened with niftiRead()'); + end + + fe = feSet(fe, 'anat offset', tempNi.qto_xyz(1:3,4)'); + clear tempNi + + % add check to file headers for consistent orientation between dwi / anat + + % load minimum data for check - ignore repeat dwi for now + doff = feGet(fe, 'dwi offset'); + aoff = feGet(fe, 'anat offset'); + + % check x dimension is consistently signed + if ~(all([ doff(1) aoff(1) ] < 0) || all([ doff(1) aoff(1) ] > 0)) + warning('POTENTIAL X-FLIP - qoffset_x is different between the dwi and anatomy files.'); + end + + % check y dimension is consistently signed + if ~(all([ doff(2) aoff(2) ] < 0) || all([ doff(2) aoff(2) ] > 0)) + warning('POTENTIAL Y-FLIP - qoffset_y is different between the dwi and anatomy files.'); + end + + % check z dimension is consistently signed + if ~(all([ doff(3) aoff(3) ] < 0) || all([ doff(3) aoff(3) ] > 0)) + warning('POTENTIAL Z-FLIP - qoffset_z is different between the dwi and anatomy files.'); + end + end %% Precompute Dictionary of orientations and canonical demeaned signals diff --git a/life/fe/feGet.m b/life/fe/feGet.m index 5b1fac3..6c0188f 100755 --- a/life/fe/feGet.m +++ b/life/fe/feGet.m @@ -1668,6 +1668,11 @@ s0 = fe.life.s0; val = ones(nTheta,1)*s0' + D*B; % with iso + case 'dwioffset' + val = fe.life.xform.dwi_offset; + case 'anatoffset' + val = fe.life.xform.anat_offset; + otherwise help('feGet') fprintf('[feGet] Unknown parameter << %s >>...\n',param); diff --git a/life/fe/feSet.m b/life/fe/feSet.m index 9d5b352..981b576 100755 --- a/life/fe/feSet.m +++ b/life/fe/feSet.m @@ -217,6 +217,11 @@ fe.life.M.tracts{Ntracts+1}.ind = find(val.index==0); fe.life.M.tracts{Ntracts+1}.name = 'not a tract'; + case 'dwioffset' + fe.life.xform.dwi_offset = val; + case 'anatoffset' + fe.life.xform.anat_offset = val; + otherwise error('Unknown parameter %s\n',param); end From 41dc560f17e3da13fec2b33858afb2b6a03ab68e Mon Sep 17 00:00:00 2001 From: bcmcpher Date: Sun, 13 Dec 2020 05:15:40 -0500 Subject: [PATCH 13/19] added flip correction as a fix if it detects one, streamlined storage of offset / header info --- life/fe/feBuildDictionaries.m | 121 +++++++++++++++++------------- life/fe/feConnectomeInit.m | 137 +++++++++++++++++++++------------- life/fe/feGet.m | 20 ++++- life/fe/feSet.m | 18 +++-- 4 files changed, 184 insertions(+), 112 deletions(-) diff --git a/life/fe/feBuildDictionaries.m b/life/fe/feBuildDictionaries.m index 56d36d4..2db5c9a 100755 --- a/life/fe/feBuildDictionaries.m +++ b/life/fe/feBuildDictionaries.m @@ -39,9 +39,6 @@ bvecs = feGet(fe,'bvecs'); % bvecs bvals = feGet(fe,'bvals'); % bvals -%bvecs(:,1) = -bvecs(:,1); % better observed x predicted correlation -%bvals = bvals * 10; % does bad - Dict = zeros(nBvecs,Norient); % Initialize Signal Dictionary matrix DictSig = zeros(nBvecs,Norient); % Initialize Signal Dictionary matrix %DictTensors = zeros(9,Norient); % Initialize Tensors Dictionary matrix @@ -49,7 +46,6 @@ % catch Kurtosis estimates for debugging %KDict = zeros(nBvecs,Norient); %KDictSig = zeros(nBvecs,Norient); -akc = zeros(nBvecs,1); % pull the shell information ubv = feGet(fe, 'nshells'); @@ -59,56 +55,70 @@ % pull diffusion tensor dt = feGet(fe, 'model tensor'); -% pull kurtosis fit - needs feGet -kt = fe.life.modelKurtosis; +% pull kurtosis fit +kt = feGet(fe, 'model Kurtosis'); -% build apparent kurtosis coefficient - akc -for i=1:ubv - - % pull tensor parameters for hard indexing of equations - sdt = dt(i,:); - - % find the indices for the shell - si = ubi == ubl(i); - - % mean diffusivity of tensor from forward model - md = mean(sdt); - - % apparent diffusion coefficient - matches dipy - adc = bvecs(si,1) .* bvecs(si,1) * sdt(1) + ... - 2 * bvecs(si,1) .* bvecs(si,2) * sdt(2) + ... - bvecs(si,2) .* bvecs(si,2) .* sdt(3) + ... - 2 * bvecs(si,1) .* bvecs(si,3) * 0 + ... % only store primary eigenvalues in - 2 * bvecs(si,2) .* bvecs(si,3) * 0 + ... % LiFE prediction. These are 0. - bvecs(si,3) .* bvecs(si,3) * 0; - - % apparent diffusion variance - matches dipy - adv = ... - bvecs(si,1) .* bvecs(si,1) .* bvecs(si,1) .* bvecs(si,1) .* kt(1) + ... % xxxx - bvecs(si,2) .* bvecs(si,2) .* bvecs(si,2) .* bvecs(si,2) .* kt(2) + ... % yyyy - bvecs(si,3) .* bvecs(si,3) .* bvecs(si,3) .* bvecs(si,3) .* kt(3) + ... % zzzz - 4 * bvecs(si,1) .* bvecs(si,1) .* bvecs(si,1) .* bvecs(si,2) .* kt(4) + ... % xxxy - 4 * bvecs(si,1) .* bvecs(si,1) .* bvecs(si,1) .* bvecs(si,3) .* kt(5) + ... % xxxz - 4 * bvecs(si,1) .* bvecs(si,2) .* bvecs(si,2) .* bvecs(si,2) .* kt(6) + ... % xyyy - 4 * bvecs(si,1) .* bvecs(si,3) .* bvecs(si,3) .* bvecs(si,3) .* kt(7) + ... % xzzz - 4 * bvecs(si,2) .* bvecs(si,2) .* bvecs(si,2) .* bvecs(si,3) .* kt(8) + ... % yyyz - 4 * bvecs(si,2) .* bvecs(si,3) .* bvecs(si,3) .* bvecs(si,3) .* kt(9) + ... % yzzz - 6 * bvecs(si,1) .* bvecs(si,1) .* bvecs(si,2) .* bvecs(si,2) .* kt(10) + ... % xxyy - 6 * bvecs(si,1) .* bvecs(si,1) .* bvecs(si,3) .* bvecs(si,3) .* kt(11) + ... % xxzz - 6 * bvecs(si,2) .* bvecs(si,2) .* bvecs(si,3) .* bvecs(si,3) .* kt(12) + ... % yyzz - 12 * bvecs(si,1) .* bvecs(si,1) .* bvecs(si,2) .* bvecs(si,3) .* kt(13) + ... % xxyz - 12 * bvecs(si,1) .* bvecs(si,2) .* bvecs(si,2) .* bvecs(si,3) .* kt(14) + ... % xyyz - 12 * bvecs(si,1) .* bvecs(si,2) .* bvecs(si,3) .* bvecs(si,3) .* kt(15); % xyzz - - % zero out noisy (bad) adc estimates - adc(adc < 0) = 0; +% if data is single shell +if ubv == 1 || all(isnan(kt)) - % estimate apparent kutosis coefficient - matches dipy - akc(si) = adv .* ((md ./ adc).^2); + % akc is invalid, set to nan + akc = nan(nBvecs,1); + +% if data is multishell +else - % zero out noisy (bad) akc estimates - akc(akc < -3/7) = -3/7; + % preallocate akc parameters to 0 + akc = zeros(nBvecs,1); + + % build apparent kurtosis coefficient - akc + for i=1:ubv + + % pull tensor parameters for hard indexing of equations + sdt = dt(i,:); + + % find the indices for the shell + si = ubi == ubl(i); + % mean diffusivity of tensor from forward model + md = mean(sdt); + + % apparent diffusion coefficient - matches dipy + adc = bvecs(si,1) .* bvecs(si,1) * sdt(1) + ... + 2 * bvecs(si,1) .* bvecs(si,2) * sdt(2) + ... + bvecs(si,2) .* bvecs(si,2) .* sdt(3) + ... + 2 * bvecs(si,1) .* bvecs(si,3) * 0 + ... % only store primary eigenvalues in + 2 * bvecs(si,2) .* bvecs(si,3) * 0 + ... % LiFE prediction. These are 0. + bvecs(si,3) .* bvecs(si,3) * 0; + + % apparent diffusion variance - matches dipy + adv = ... + bvecs(si,1) .* bvecs(si,1) .* bvecs(si,1) .* bvecs(si,1) .* kt(1) + ... % xxxx + bvecs(si,2) .* bvecs(si,2) .* bvecs(si,2) .* bvecs(si,2) .* kt(2) + ... % yyyy + bvecs(si,3) .* bvecs(si,3) .* bvecs(si,3) .* bvecs(si,3) .* kt(3) + ... % zzzz + 4 * bvecs(si,1) .* bvecs(si,1) .* bvecs(si,1) .* bvecs(si,2) .* kt(4) + ... % xxxy + 4 * bvecs(si,1) .* bvecs(si,1) .* bvecs(si,1) .* bvecs(si,3) .* kt(5) + ... % xxxz + 4 * bvecs(si,1) .* bvecs(si,2) .* bvecs(si,2) .* bvecs(si,2) .* kt(6) + ... % xyyy + 4 * bvecs(si,1) .* bvecs(si,3) .* bvecs(si,3) .* bvecs(si,3) .* kt(7) + ... % xzzz + 4 * bvecs(si,2) .* bvecs(si,2) .* bvecs(si,2) .* bvecs(si,3) .* kt(8) + ... % yyyz + 4 * bvecs(si,2) .* bvecs(si,3) .* bvecs(si,3) .* bvecs(si,3) .* kt(9) + ... % yzzz + 6 * bvecs(si,1) .* bvecs(si,1) .* bvecs(si,2) .* bvecs(si,2) .* kt(10) + ... % xxyy + 6 * bvecs(si,1) .* bvecs(si,1) .* bvecs(si,3) .* bvecs(si,3) .* kt(11) + ... % xxzz + 6 * bvecs(si,2) .* bvecs(si,2) .* bvecs(si,3) .* bvecs(si,3) .* kt(12) + ... % yyzz + 12 * bvecs(si,1) .* bvecs(si,1) .* bvecs(si,2) .* bvecs(si,3) .* kt(13) + ... % xxyz + 12 * bvecs(si,1) .* bvecs(si,2) .* bvecs(si,2) .* bvecs(si,3) .* kt(14) + ... % xyyz + 12 * bvecs(si,1) .* bvecs(si,2) .* bvecs(si,3) .* bvecs(si,3) .* kt(15); % xyzz + + % zero out noisy (bad) adc estimates + adc(adc < 0) = 0; + + % estimate apparent kutosis coefficient - matches dipy + akc(si) = adv .* ((md ./ adc).^2); + + % zero out noisy (bad) akc estimates + akc(akc < -3/7) = -3/7; + + end + end % Compute each dictionary column for a different kernel orientation @@ -130,9 +140,14 @@ % estimate Q for tensor values in shell Q = Rot*D*Rot'; + % STORE SEPARATE? THIS WILL GREATLY INCREASE THE SIZE + % Compute the signal contribution of a fiber in the kernel orientation divided S0 - Dict(si,j) = exp(-bvals(si) .* diag(bvecs(si,:)*Q*bvecs(si,:)')); - %Dict(si,j) = exp(-bvals(si) .* diag(bvecs(si,:)*Q*bvecs(si,:)') + (-bvals(si).^2 .* diag(bvecs(si,:)*Q*bvecs(si,:)').^2 .* akc(si))/6); + if all(kt == 0) || all(isnan(kt)) % if kurtosis is zeroed b/c it's invalid / requested + Dict(si,j) = exp(-bvals(si) .* diag(bvecs(si,:)*Q*bvecs(si,:)')); + else % if kurtosis parameters are passed + Dict(si,j) = exp(-bvals(si) .* diag(bvecs(si,:)*Q*bvecs(si,:)') + (-bvals(si).^2 .* diag(bvecs(si,:)*Q*bvecs(si,:)').^2 .* akc(si))/6); + end % demeaned signal by shell DictSig(si,j) = Dict(si,j) - mean(Dict(si,j)); @@ -143,10 +158,10 @@ end fe = feSet(fe,'dictionary parameters',{Nphi,Ntheta,orient,Dict,DictSig}); +fe = feSet(fe,'akc',akc); % hard add kurtosis dictionaries / akc for debugging %fe.life.M.KDict = KDict; %fe.life.M.KDictSig = KDictSig; -fe.life.M.akc = akc; end \ No newline at end of file diff --git a/life/fe/feConnectomeInit.m b/life/fe/feConnectomeInit.m index 43b0e49..9d2b738 100755 --- a/life/fe/feConnectomeInit.m +++ b/life/fe/feConnectomeInit.m @@ -1,4 +1,4 @@ -function fe = feConnectomeInit(dwiFile,fgFileName,feFileName,savedir,dwiFileRepeated,anatomyFile,varargin) +function fe = feConnectomeInit(dwiFile,fgFileName,dwiFileRepeated,anatomyFile,N,tensor,kurtosis,reset) % Initialize a new connectome (fe) structure. % % fe = feConnectomeInit(dwiFile,fgFileName,feFileName,savedir,dwiFileRepeated,anatomyFile,varargin); @@ -13,18 +13,35 @@ %feOpenLocalCluster +if(~exist('N', 'var') || isempty(N)) + N = 360; +end + +if(~exist('tensor', 'var') || isempty(tensor)) + axialDiffusion = 1; radialDiffusion = 0; + tensor = [ axialDiffusion radialDiffusion radialDiffusion ]; +end + +if(~exist('kurtosis', 'var') || isempty(kurtosis)) + kurtosis = zeros(15,1); +end + +if(~exist('reset', 'var') || isempty(reset)) + reset = true; +end + % Intialize the fe structure. fe = feCreate; -% Set the based dir for fe, this dire will be used -if notDefined('savedir') - if isstruct(fgFileName) - savedir = fullfile(fileparts(fgFileName.name),'life'); - else - savedir = fullfile(fileparts(fgFileName),'life'); - end -end -fe = feSet(fe,'savedir',savedir); +% Set the based dir for fe, this dir will be used +% if notDefined('savedir') +% if isstruct(fgFileName) +% savedir = fullfile(fileparts(fgFileName.name),'life'); +% else +% savedir = fullfile(fileparts(fgFileName),'life'); +% end +% end +% fe = feSet(fe,'savedir',savedir); % Set the xforms (transformations from diffusion data to acpc) tempNi = niftiRead(dwiFile); @@ -32,18 +49,18 @@ fe = feSet(fe, 'acpc2img xform', inv(tempNi.qto_xyz)); % set the offset for files as well, then clear it -fe = feSet(fe, 'dwi offset', tempNi.qto_xyz(1:3,4)'); +%fe = feSet(fe, 'dwi offset', tempNi.qto_xyz(1:3,4)'); clear tempNi -% Set up the fe name -if isstruct(fgFileName), n = fgFileName.name; -else [~,n] = fileparts(fgFileName); -end - -if notDefined('feFileName') - feFileName = sprintf('%s-%s', datestr(now,30),n); -end -fe = feSet(fe, 'name',feFileName); +% % Set up the fe name +% if isstruct(fgFileName), n = fgFileName.name; +% else [~,n] = fileparts(fgFileName); +% end +% +% if notDefined('feFileName') +% feFileName = sprintf('%s-%s', datestr(now,30),n); +% end +% fe = feSet(fe, 'name',feFileName); %% Load a connectome from disk if isstruct(fgFileName), fg = fgFileName; clear fgFileName @@ -70,13 +87,13 @@ try tempNi = niftiRead(anatomyFile); + fe = feSet(fe, 'anat2acpcxform', tempNi.qto_xyz); + fe = feSet(fe, 'acpc2anatxform', inv(tempNi.qto_xyz)); + clear tempNi catch error('Anatomy file unable to be opened with niftiRead()'); end - - fe = feSet(fe, 'anat offset', tempNi.qto_xyz(1:3,4)'); - clear tempNi - + % add check to file headers for consistent orientation between dwi / anat % load minimum data for check - ignore repeat dwi for now @@ -85,52 +102,70 @@ % check x dimension is consistently signed if ~(all([ doff(1) aoff(1) ] < 0) || all([ doff(1) aoff(1) ] > 0)) - warning('POTENTIAL X-FLIP - qoffset_x is different between the dwi and anatomy files.'); + if reset + warning('POTENTIAL X-FLIP - bvecs flipped on x-axis to correct.'); + bvecs = feGet(fe, 'bvecs'); + bvecs(:,1) = -bvecs(:,1); + fe = feSet(fe, 'bvecs', bvecs); + else + warning('POTENTIAL X-FLIP - qoffset_x is different between the dwi and anatomy files.'); + end end % check y dimension is consistently signed if ~(all([ doff(2) aoff(2) ] < 0) || all([ doff(2) aoff(2) ] > 0)) - warning('POTENTIAL Y-FLIP - qoffset_y is different between the dwi and anatomy files.'); + if reset + warning('POTENTIAL Y-FLIP - bvecs flipped on y-axis to correct.'); + bvecs = feGet(fe, 'bvecs'); + bvecs(:,2) = -bvecs(:,2); + fe = feSet(fe, 'bvecs', bvecs); + else + warning('POTENTIAL Y-FLIP - qoffset_y is different between the dwi and anatomy files.'); + end end % check z dimension is consistently signed if ~(all([ doff(3) aoff(3) ] < 0) || all([ doff(3) aoff(3) ] > 0)) - warning('POTENTIAL Z-FLIP - qoffset_z is different between the dwi and anatomy files.'); + if reset + warning('POTENTIAL Z-FLIP - bvecs flipped on z-axis to correct.'); + bvecs = feGet(fe, 'bvecs'); + bvecs(:,3) = -bvecs(:,3); + fe = feSet(fe, 'bvecs', bvecs); + else + warning('POTENTIAL Z-FLIP - qoffset_z is different between the dwi and anatomy files.'); + end end end +%% check / fix the tensor model for multishell input + +nshells = feGet(fe, 'nshells'); +if (size(tensor, 1) == 1) && nshells > 1 + tensorParams = repmat(tensor, [ nshells 1 ]); +else + tensorParams = tensor; +end + +% if data is single shell, zero out kurtosis model regardless of user input +if nshells == 1 + kurtosisParams = nan(15,1); +else + kurtosisParams = kurtosis; +end + %% Precompute Dictionary of orientations and canonical demeaned signals % Define discretization steps for building Dictionary of canonical difussivities % These numbers represents the steps in which the range [0,pi] is divided % for spherical coordinates phi: azimuth and theta: pi/2-elevetion -% Set model for canonical tensor -%tic -if ~isempty(varargin) - N = varargin{1}(1); - tensorParams = varargin{2}; - kurtosisPars = varargin{3}; - %axialDiffusion = varargin{2}(1); - %radialDiffusion = varargin{2}(2); - %radialDiffusion2 = varargin{2}(3); -else % Default to stick and ball - N = 360; - % very default tensor settings - axialDiffusion = 1; - radialDiffusion = 0; - tensorModel = [ axialDiffusion radialDiffusion radialDiffusion ]; - tensorParams = []; - kurtosisPars = []; - nshells = feGet(fe, 'nshells'); - for i=1:nshells - tensorParams = [ tensorParams; tensorModel ]; - end -end Nphi = N; Ntheta = N; -fe = feSet(fe,'model tensor',tensorParams); % TENSOR PARAMS NEEDS TO BE SANITIZED BETTER -fe.life.modelKurtosis = kurtosisPars; % add as feSet fxn? +fe = feSet(fe,'model tensor',tensorParams); +fe = feSet(fe,'model kurtosis',kurtosisParams); +% if single shell data is passed (impossible to fit kurtosis), +% feBuildDictionaries will zero out any submitted kurtosis parameters + fe = feBuildDictionaries(fe,Nphi,Ntheta); %% Encode Connectome in multidimensional tensor diff --git a/life/fe/feGet.m b/life/fe/feGet.m index 6c0188f..3fe14e8 100755 --- a/life/fe/feGet.m +++ b/life/fe/feGet.m @@ -1667,11 +1667,25 @@ s0 = fe.life.s0; val = ones(nTheta,1)*s0' + D*B; % with iso - + + case {'modelkurtosis'} + val = fe.life.modelKurtosis; + case {'akc'} + val = fe.life.M.akc; + case {'dwi2acpc'} + val = fe.life.xform.img2acpc; + case {'dwi2img'} + val = fe.life.xform.acpc2img; + case {'anat2acpc'} + val = fe.life.xform.anat2acpc; + case {'anat2img'} + val = fe.life.xform.anat2img; case 'dwioffset' - val = fe.life.xform.dwi_offset; + tmp = feGet(fe, 'dwi2acpc'); + val = tmp(1:3,4)'; case 'anatoffset' - val = fe.life.xform.anat_offset; + tmp = feGet(fe, 'anat2acpc'); + val = tmp(1:3,4)'; otherwise help('feGet') diff --git a/life/fe/feSet.m b/life/fe/feSet.m index 981b576..e15a8c6 100755 --- a/life/fe/feSet.m +++ b/life/fe/feSet.m @@ -217,11 +217,19 @@ fe.life.M.tracts{Ntracts+1}.ind = find(val.index==0); fe.life.M.tracts{Ntracts+1}.name = 'not a tract'; - case 'dwioffset' - fe.life.xform.dwi_offset = val; - case 'anatoffset' - fe.life.xform.anat_offset = val; - + case {'xformanat2acpc','anat2acpc','anat2acpcxform'} + fe.life.xform.anat2acpc = val; + case {'xformacpc2anat','acpc2anat','acpc2anatxform'} + fe.life.xform.acpc2anat = val; +% case 'dwioffset' +% fe.life.xform.dwi_offset = val; +% case 'anatoffset' +% fe.life.xform.anat_offset = val; + case {'modelkurtosis'} + fe.life.modelKurtosis = val; + case {'akc'} + fe.life.M.akc = val; + otherwise error('Unknown parameter %s\n',param); end From be97e7b751d31daaf87d60c8223a4d4691b80991 Mon Sep 17 00:00:00 2001 From: bcmcpher Date: Thu, 11 Mar 2021 17:07:05 -0500 Subject: [PATCH 14/19] added checks to loading multishell info to fill in old single shell structures --- life/fe/feGet.m | 18 ++++++++++++++++-- life/fe/feSet.m | 4 ++-- 2 files changed, 18 insertions(+), 4 deletions(-) diff --git a/life/fe/feGet.m b/life/fe/feGet.m index 3fe14e8..f14b5fd 100755 --- a/life/fe/feGet.m +++ b/life/fe/feGet.m @@ -506,13 +506,27 @@ % This is the number of unique shells for multishell data % % val = feGet(fe,'nshells'); - val = fe.life.shells_n; + + % if nshells isn't present, it's an old single shell structure + if ~isfield(fe.life, 'nshells') + val = 1; % fill in structure and set the value to 1 + feSet(fe, 'nshells', val); % add it to the structure + else % otherwise pull the value set during encoding + val = fe.life.nshells; + end case {'shellindex','indextoeachshellbvecs'} % This is the index to the bvecs associated to each unique shell % % val = feGet(fe,'nshshellindexells'); - val = fe.life.shells_index; + + % if shellsindex isn't present, it's an old single shell structure + if ~isfield(fe.life, 'shells_index') + val = ones(feGet(fe, 'nbvals'), 1); % fill in all ones for the index + feSet(fe, 'shellindex', val); % add it to the structure + else % otherwise pull the value set during encoding + val = fe.life.shells_index; + end case {'bvecs'} % Diffusion directions. diff --git a/life/fe/feSet.m b/life/fe/feSet.m index e15a8c6..ad6e2f5 100755 --- a/life/fe/feSet.m +++ b/life/fe/feSet.m @@ -90,7 +90,7 @@ fe.life.bvals = val; case {'nshells','numberofshells'} % This is the number of unique shells for multishell data - fe.life.shells_n = val; + fe.life.nshells = val; case {'shellindex','indextoeachshellbvecs'} % This is the index to the bvecs associated to each unique shell fe.life.shells_index = val; @@ -170,7 +170,7 @@ fe.rep.bvals = val; case {'nshellsrepeat','numberofshellsrepeat'} % This is the number of unique shells for multishell data - fe.rep.shells_n = val; + fe.rep.nshells = val; case {'shellindexrepeat','indextoeachshellbvecsrepeat'} % This is the index to the bvecs associated to each unique shell fe.rep.shells_index = val; From 0d1d222d551d2d6b11cafcf77a05b66f03a12647 Mon Sep 17 00:00:00 2001 From: bcmcpher Date: Fri, 12 Mar 2021 11:27:13 -0500 Subject: [PATCH 15/19] add an isdeployed check around mex fxns so they can be compiled --- life/compute/feFitModel.m | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/life/compute/feFitModel.m b/life/compute/feFitModel.m index 54f7edb..a9b4539 100755 --- a/life/compute/feFitModel.m +++ b/life/compute/feFitModel.m @@ -56,12 +56,14 @@ % Check for mexfiles and generate them if necessary -% Mtransp_times_b function -checkMexCompiled('-largeArrayDims', '-output', 'Mtransp_times_b', '-DNDEBUG','Mtransp_times_b.c', 'Mtransp_times_b_sub.c') -% M_times_w function -checkMexCompiled('-largeArrayDims', '-output', 'M_times_w', '-DNDEBUG', 'M_times_w.c', 'M_times_w_sub.c') -% compute_diag function -checkMexCompiled('-largeArrayDims', '-output', 'compute_diag', '-DNDEBUG', 'compute_diag.c', 'compute_diag_sub.c') +if ~isdeployed + % Mtransp_times_b function + checkMexCompiled('-largeArrayDims', '-output', 'Mtransp_times_b', '-DNDEBUG','Mtransp_times_b.c', 'Mtransp_times_b_sub.c') + % M_times_w function + checkMexCompiled('-largeArrayDims', '-output', 'M_times_w', '-DNDEBUG', 'M_times_w.c', 'M_times_w_sub.c') + % compute_diag function + checkMexCompiled('-largeArrayDims', '-output', 'compute_diag', '-DNDEBUG', 'compute_diag.c', 'compute_diag_sub.c') +end if nargin <6 % no initial w0 is provided [nFibers] = size(M.Phi,3); From d9c5d94fbd01778612d4a0e6a6f9421dc653406f Mon Sep 17 00:00:00 2001 From: bcmcpher Date: Mon, 29 Mar 2021 05:25:11 -0400 Subject: [PATCH 16/19] renamed fxn for sanity in debugging --- .../{feComputeEvidence_new.m => feComputeEvidence_fixed.m} | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) rename life/compute/{feComputeEvidence_new.m => feComputeEvidence_fixed.m} (99%) diff --git a/life/compute/feComputeEvidence_new.m b/life/compute/feComputeEvidence_fixed.m similarity index 99% rename from life/compute/feComputeEvidence_new.m rename to life/compute/feComputeEvidence_fixed.m index a5df0cf..ebf53b4 100644 --- a/life/compute/feComputeEvidence_new.m +++ b/life/compute/feComputeEvidence_fixed.m @@ -1,4 +1,4 @@ -function se = feComputeEvidence_new(rmse1, rmse2, nbins, nmc, nboots) +function se = feComputeEvidence_fixed(rmse1, rmse2, nbins, nmc, nboots) % Computes a series of distance metrics between two RMSE distribtions. % % Compute summary statistics on to characterize the lesion,: From b2275dac36d14cf6c22cf392cffaf946b02088b7 Mon Sep 17 00:00:00 2001 From: bcmcpher Date: Fri, 2 Apr 2021 15:17:36 -0400 Subject: [PATCH 17/19] added function to convert 3d tensor to 2d matrix for optimization --- life/utility/feSparseMatrixCoords.m | 83 +++++++++++++++++++++++++++++ 1 file changed, 83 insertions(+) create mode 100644 life/utility/feSparseMatrixCoords.m diff --git a/life/utility/feSparseMatrixCoords.m b/life/utility/feSparseMatrixCoords.m new file mode 100644 index 0000000..ff1391c --- /dev/null +++ b/life/utility/feSparseMatrixCoords.m @@ -0,0 +1,83 @@ +function [ out ] = feSparseMatrixCoords(M) +%[ mat ] = feSparseMatrixCoords(Phi) +% Convert the ENCODE model object (Phi) and the dictionary signal (D) to +% a set of coordinates to be used in the 2d sparse matrix version of +% LiFE that incorporates streamline groupings during optimization. +% +% INPUTS: +% M - The model object from an initialized ENCODE object +% example: M = feGet(fe, 'model'; +% +% OUTPUT: +% out - The conversion of the 3d sparse tensor and dictionary to the +% coordinates / values required for the 2d sparse matrix LiFE model +% +% Brent McPherson put this function together, but +% Andy Womack and Dan McDonald did most of the driving. +% +% Brent McPherson, Indiana University (c) 2021 +% + +% pull the number of streamlines +nfib = size(M.Phi, 3); + +% Each slice (streamline) of Phi needs to be multiplied by the dictionary +% signal to create the 2d matrix for of each streamlines prediction. The +% streamlines can be stacked into a single 2d matrix by adding a column +% indicating which streamline the entries (subscripts, values) they +% correspond to. + +disp([ 'Converting ' num2str(nfib) ' streamlines in Phi to sparse matrix coordinates...' ]); + +% initialize the output +out = []; +%out = nan(size(M.Phi.vals,1), 4); % preallocate output to final size + +% for every streamline +for fib = 1:nfib + + % pull the slice of the tensor + slc = M.Phi(:,:,fib); + + % use the sptensor toolbox to multiply the slice by the dictionary for + % making each streamlines prediction + mm = ttm(slc, M.DictSig, 1); + + % create the index to identify the streamline in the new matrix + idx = ones(size(mm.vals)) * ii; + + % for this streamline, store the values from the sparse product: + % - x (dictionary) / y (voxel) subscripts into the sparse model matrix + % - the index of the streamline being for this value + % - the value stored at that entry + iter = [ mm.subs, idx, mm.vals ]; + + % append this streamline to the output + out = [ out; iter ]; + %out(zzz, :) = iter; % need zzz to be indices of out for the streamline + +end + +% +% concerns / things to note +% +% mm.subs is type double from ttm - why? +% - this is dumb, these should be int - ought to fix +% - indices are integers, why store them as doubles? +% - BE CAREFUL WITH THE PRECISION WHEN SAVING THIS OUTPUT +% -- because they're saved as doubles, the integer indices will round if +% the precision isn't high enough, resulting in errors. +% +% In theory, the sptensor toolbox should be able to transform the tensor +% to a matrix with a function distributed in the tool. However, that +% doesn't appear to be the case (their documentation / notation is bad). +% +% Further, the tensor-matrix multiplication doesn't work for these types +% (handling sparse objects sparsely). So many attempts to work on the 3d +% tensor directly resulted in comically large memory overflows. Which is +% weird, because they have a @sptensor/ttm function that ought to handle +% sparse tensors correctly. Maybe it worked in an older version of matlab? +% + +end + From aa44cb8bcec88c8dc3ac8637ccc8ad380927f723 Mon Sep 17 00:00:00 2001 From: bcmcpher Date: Fri, 2 Apr 2021 16:40:06 -0400 Subject: [PATCH 18/19] updated documentation --- life/utility/feSparseMatrixCoords.m | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/life/utility/feSparseMatrixCoords.m b/life/utility/feSparseMatrixCoords.m index ff1391c..6eec134 100644 --- a/life/utility/feSparseMatrixCoords.m +++ b/life/utility/feSparseMatrixCoords.m @@ -6,7 +6,7 @@ % % INPUTS: % M - The model object from an initialized ENCODE object -% example: M = feGet(fe, 'model'; +% example: M = feGet(fe, 'model'); % % OUTPUT: % out - The conversion of the 3d sparse tensor and dictionary to the @@ -31,12 +31,12 @@ % initialize the output out = []; -%out = nan(size(M.Phi.vals,1), 4); % preallocate output to final size +%out = nan(xxx, 4); % preallocate output for speed, but I don't know what that is... % for every streamline for fib = 1:nfib - % pull the slice of the tensor + % pull the slice of the model tensor for the streamline slc = M.Phi(:,:,fib); % use the sptensor toolbox to multiply the slice by the dictionary for @@ -44,15 +44,15 @@ mm = ttm(slc, M.DictSig, 1); % create the index to identify the streamline in the new matrix - idx = ones(size(mm.vals)) * ii; + idx = ones(size(mm.vals)) * fib; % for this streamline, store the values from the sparse product: % - x (dictionary) / y (voxel) subscripts into the sparse model matrix - % - the index of the streamline being for this value - % - the value stored at that entry + % - the index of the streamline for these values + % - the value stored at each entry iter = [ mm.subs, idx, mm.vals ]; - - % append this streamline to the output + + % append this streamline to the output matrix out = [ out; iter ]; %out(zzz, :) = iter; % need zzz to be indices of out for the streamline @@ -73,10 +73,10 @@ % doesn't appear to be the case (their documentation / notation is bad). % % Further, the tensor-matrix multiplication doesn't work for these types -% (handling sparse objects sparsely). So many attempts to work on the 3d -% tensor directly resulted in comically large memory overflows. Which is -% weird, because they have a @sptensor/ttm function that ought to handle -% sparse tensors correctly. Maybe it worked in an older version of matlab? +% (handling sparse objects sparsely). So attempts to work on the 3d tensor +% directly resulted in comically large memory overflows. Which is weird, +% because they have a @sptensor/ttm function that ought to handle sparse +% tensors correctly. Maybe it worked in an older version of matlab? % end From 0a8af27f735b3ee54c4edd557f63f8310e5ad3e4 Mon Sep 17 00:00:00 2001 From: Brent McPherson Date: Fri, 2 Apr 2021 17:20:43 -0400 Subject: [PATCH 19/19] preallocate and print useful updates in feSparseMatrixCoords --- life/utility/feSparseMatrixCoords.m | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/life/utility/feSparseMatrixCoords.m b/life/utility/feSparseMatrixCoords.m index 6eec134..adf1ecf 100644 --- a/life/utility/feSparseMatrixCoords.m +++ b/life/utility/feSparseMatrixCoords.m @@ -27,15 +27,20 @@ % indicating which streamline the entries (subscripts, values) they % correspond to. -disp([ 'Converting ' num2str(nfib) ' streamlines in Phi to sparse matrix coordinates...' ]); +fprintf('Converting %d streamlines in Phi to sparse matrix coordinates...\n', nfib); % initialize the output -out = []; -%out = nan(xxx, 4); % preallocate output for speed, but I don't know what that is... +%out = []; +out = cell(nfib, 1); % for every streamline for fib = 1:nfib + % print out every 1000th streamline as an update + if mod(fib, 1000) == 0 + fprintf('Running streamline %d\n', fib); + end + % pull the slice of the model tensor for the streamline slc = M.Phi(:,:,fib); @@ -53,11 +58,16 @@ iter = [ mm.subs, idx, mm.vals ]; % append this streamline to the output matrix - out = [ out; iter ]; - %out(zzz, :) = iter; % need zzz to be indices of out for the streamline + %out = [ out; iter ]; + out{fib} = iter; end +fprintf('Stacking the streamline data into sparse matrix coordinates...\n'); + +% merge the output +out = cat(1, out{:}); + % % concerns / things to note %