【信号处理】贝叶斯和分段IVIM模型拟合matlab代码

%% % %%% %%% %%% %%% %%% %%% %%% % % Use of IVIMmodelfit % %%% %%% %%% %%% %%% %%% %%% %% clc clear all close all im = "IVIMphantom_brain_SNR200.nii"; % can be 4 D matrix mask = "IVIMphantom_brain_mask.nii" ; % empty ([] or ’’ ) means no mask bvals = "bvals.txt" ; % can be vector matching 4 th dimension of im fittype = "seg" ; % ’seg’ or ’bayes’ , default ’seg’ % additional options specific for the fittype and if file input is used opts.outfile = "IVIMmaps" ; opts.its = 10 ; % should be larger opts.burns = 10 ; % should be larger opts.lim = [ 0 0 0 0 ; 3 e- 3 200 1 1 ]; opts.lim = opts.lim( : , 1 : 4 ); % 1 : 2 = ADC & S 0 , 1 : 3 = also f, 1 : 4 = also D* opts.bthr = 200 ; % threshold for segmented fitting = 200 s/mm2, for monoexp fit (ADC) = 0 s/mm2 % run model fit IVIMmodelfit(im,bvals,fittype,mask,opts); % display results maps = niftiread(opts.outfile+ ".nii" ); truepars = niftiread( "IVIMphantom_brain_truepars.nii" ); truepars = truepars( : , : , : ,[ 1 4 2 3 ]); displims = [ 0 0 0 0 ; 1.5 e- 3 200 0 . 2 0 . 05 ]; pars = { ’D’ , ’S_0’ , ’f’ , ’D*’ }; figure; for i = 1 :size (opts.lim, 2 ) subplot( 2 ,size(opts.lim, 2 ),i); imshow(rot9 0 (maps( : , : , 50 ,i), 3 ),displims( : ,i)); title(pars{i}) if i == 1 ylabel( ’Estimated parameter maps’ ); end subplot( 2 ,size(opts.lim, 2 ),i+size(opts.lim, 2 )); imshow(rot9 0 (truepars( : , : , 50 ,i), 3 ),displims( : ,i)); if i == 1 ylabel( ’True parameter maps’ ); end end colormap jet; %% % %%% %%% %%% %%% %%% %%% %%% %%% %%% %% % Generate data for examples % %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% % Input parameters n = 1 e3; % number of simulated data series D = 1 e- 3 *ones(n, 1 ); f = 0 . 3 *ones(n, 1 ); Dstar = 20 e- 3 *ones(n, 1 ); S 0 = 100 *ones(n, 1 ); b = [ 0 10 20 30 40 50 75 100 200 400 600 ]; % Generate data and add Rician noise Y_true = repmat(S 0 , 1 ,length(b)).*(( 1 -repmat(f, 1 ,length(b))).*exp(-D*b) + repmat(f, 1 ,length(b)).*exp(-Dstar*b)); SNR = 20 ; % difference between bayesian and segmented model fitting is more pronounced at low SNR (< 15 ) Y = abs(Y_true + repmat(S 0 , 1 ,length(b))/SNR .* ( randn(n,length(b)) + 1 i*randn(n,length(b)) )); % If you have data from an 4 D image (b-values in 4 th dimension) you can use % im2Y to get the correct format as: % Y = im2Y(im); % if you also have a 3 D mask (ROI) you can use it as a second argument: % Y = im2Y(im,mask); %% % %%% %%% %%% %%% %%% %%% %%% %%% %%% % % Example use of IVIM bayes % %%% %%% %%% %%% %%% %%% %%% %%% %%% %% lim = [0 0 0 0;1 3e-3 100e-3 200]; % order: f,D,D*,S 0 its = 1 e4; rician = false ; prior = { ’flat’ , ’lognorm’ , ’lognorm’ , ’flat’ }; burns = 1000 ; meanonly = false ; % Specifies all input arguments res_bayes1 = IVIM_bayes(Y,f,D,Dstar,S 0 ,b,lim,its,rician,prior,burns,meanonly); % Uses default arguments (including uniform priors) and scalar inputs for f, D, D* and S 0 res_bayes2 = IVIM_bayes(Y, 0 . 3 , 1 e- 3 , 20 e- 3 , 100 ,b,lim); %% % %%% %%% %%% %%% %%% %%% %%% %%% %% % Example use of IVIM_seg % %%% %%% %%% %%% %%% %%% %%% %%% %%% blim = 200 ; lim = [ 0 0 0 0 ; 3 e- 3 200 1 100 e- 3 ]; % NOTE! order here is D,S 0 ,f,D* disp_prog = true ; % Segmented fitting for biexponential IVIM, specifies all arguments res_seg1 = IVIM_seg(Y,b,lim,blim,disp_prog); % Estimation of on D and f through segmented fitting, using default arguments res_seg2 = IVIM_seg(Y,b,lim( : , 1 : 3 )); % Estimation of ADC (uses all b-values by default) res_seg3 = IVIM_seg(Y,b,lim( : , 1 : 2 )); %% % %%% %%% %%% %%% %%% % Plot results % %%% %%% %%% %%% %%% % figure; pars = { ’f’ , ’D’ , ’Dstar’ , ’S0’ }; for i = 1 :length (pars) subplot( 2 ,ceil(length(pars)/ 2 ),i); X = zeros(n, 7 ); X( : , 1 : 5 ) = ... [res_bayes1.(pars{i}).mode ’ res_bayes1.(pars{i}).mean’ ... res_bayes2.(pars{i}).mode ’ res_bayes2.(pars{i}).mean’ ... res_seg1.(pars{i}) ’]; if i ~= 3 X(:,end-1) = res_seg2.(pars{i})’ ; else X( : , end - 1 ) = nan; end if ~(i == 1 || i == 3 ) X( : , end ) = res_seg3.(pars{i}); else X( : , end ) = nan; end boxplot(X, ’labels’ ,{ ’logn mode’ , ’logn mean’ , ’flat mode’ ,... ’flat mean’ , ’seg biexp’ , ’seg D&f’ , ’seg ADC’ }); title(pars{i}); end

function hsm = halfSampleMode(X) % calculates the half sample mode for each row in X if ~ismatrix(X) error(’X must be a matrix or vector’); end X = sort(X,2); n = size(X,2); hsm = HSM_rec(n,X); function hsm = HSM_rec(n,X) % special cases if size(X,2) == 1 hsm = X; return; elseif size(X,2) == 2 hsm = sum(X,2)/2; return; elseif size(X,2) == 3 hsm = zeros(size(X,1),1); low = (X(:,2) - X(:,1)) < (X(:,3) - X(:,2)); % use lower pair eq = (X(:,2) - X(:,1)) == (X(:,3) - X(:,2)); % use mid value if any(low) hsm(low) = sum(X(low,1:2),2)/2; end if any(eq) hsm(eq) = X(eq,2); end if any(~(low|eq)) hsm(~(low|eq)) = sum(X(~(low|eq),2:3),2)/2; % otherwise end return end % general case (n > 3) wmin = X(:,end) - X(:,1); N = ceil(n/2); j = ones(size(X,1),1); for i = 1:(n-N+1) w = X(:,i+N-1) - X(:,i); m = w < wmin; wmin(m) = w(m); j(m) = i; end rowsub = repmat((1:size(X,1))’,1,N); colsub = repmat(j,1,N) + repmat(0:N-1,size(X,1),1); Xsub = reshape(X(sub2ind(size(X),rowsub,colsub)),size(X,1),N); hsm = HSM_rec(N,Xsub);

function Y = im2Y(im,mask) % Transforms functional image data ( 4 D or 3 D array) into % data matrix with size VxT where V is the number of % voxels and T is the number of elements in the functional % dimensionon (e.g. time, TE, b-value). If a mask is supplied % only voxels within the mask are extracted. s = size(im); if numel(s) > 2 if nargin > 1 V = sum(mask( : )); else if numel(s) == 4 V = numel(im( : , : , : , 1 )); else V = numel(im( : , : , 1 )); end end Y = zeros(V,s( end )); for i = 1 :s ( end ) if nargin > 1 %mask_full = false ([size(mask) s( end )]);mask_full( : , : ,i) = mask; temp = im( : , : , : ,i); temp = temp(mask); else if numel(s) == 4 temp = im( : , : , : ,i); else temp = im( : , : ,i); end end Y( : ,i) = temp( : ); end else Y = im; end

function out = IVIM_bayes(Y,f,D,Dstar,S 0 ,b,lim,n,rician,prior,burns,meanonly) % out = IVIM_bayes(Y,f,D,Dstar,S 0 ,b,N,lim,rician,prior,burns,meanonly) % % Input: % Y: a VxB matrix containing all image data where V is the number of % voxels in the analysis ROI and B is the number of b-values % f: a Vx1 vector containg the perfusion fraction parameters for all % voxels from a least squares fit % D: a Vx1 vector containg the diffusion parameters for all voxels % from a least squares fit % Dstar: a Vx1 vector containg the pseudo diffusion parameters for % all voxels from a least squares fit % b: a 1 xB vector containing all b-values % lim: a 2 x4 or 2 x5 matrix with lower ( 1 st row) and upper ( 2 nd row) % limits of all parameters in the order f,D,D*,S 0 ,noisevar % n: number of iterations after "burn in" ( default: 10000 ) % rician: scalar boolean. If true Rician noise distribution is used, % else Gaussian noise distribution is used ( default: false ) % prior: cell of size equal to number of estimated parameters % the cells should contain string equal to ’flat’ , ’reci’ or % ’lognorm’ to define the prior disitribution for the % corresponding parameter ( default: { ’flat’ , ’flat’ , ’flat’ , ’flat’ , ’reci’ } % ’flat’ = uniform, ’reci’ = reciprocal, ’lognorm’ = lognormal % lognormal is only available for D and D* % burns: number of burn- in steps ( default: 1000 ) % meanonly: scalar boolean. if true , only the posterior mean and std % are estimated (substantially more memory efficient and % therefore faster) ( default: false ) % % Output: % out: a struct with the fields D, f, Dstar and S 0 (Vx1) containing the % voxelwise mean, median, mode and standard deviation of each parameter % % By Oscar Jalnefjord 2018 -08- 27 % % If you use this function in research, please cite: % Gustafsson et al. 2017 Impact of prior distributions and central tendency % measures on Bayesian intravoxel incoherent motion model fitting, MRM %%% %%% %%% %%% %%% %%% % Error handling % %%% %%% %%% %%% %%% %%% if ~isrow(b) b = b ’; if ~isrow(b) error(’ b must be a row vector ’); end end B = length(b); if size(Y,2) ~= B Y = Y’ ; if size(Y, 2 ) ~= B error( ’Y must have same number of columns as b’ ); end end V = size(Y, 1 ); if isvector(f) && length(f) > 1 if ~isequal([V 1 ],size(f)) f = f ’; if ~isequal([V 1],size(f)) error(’ f must be a column vector with the same number of rows as Y or a scalar ’); end end elseif isscalar(f) f = f*ones(V,1); else error(’ f must be a column vector or a scalar ’); end if isvector(D) && length(D) > 1 if ~isequal([V 1],size(D)) D = D’ ; if ~isequal([V 1 ],size(D)) error( ’D must be a column vector with the same number of rows as Y or a scalar’ ); end end elseif isscalar(D) D = D*ones(V, 1 ); else error( ’D must be a column vector or a scalar’ ); end if isvector(Dstar) && length(Dstar) > 1 if ~isequal([V 1 ],size(Dstar)) Dstar = Dstar ’; if ~isequal([V 1],size(Dstar)) error(’ Dstar must be a column vector with the same number of rows as Y or a scalar ’); end end elseif isscalar(Dstar) Dstar = Dstar*ones(V,1); else error(’ Dstar must be a column vector or a scalar ’); end if isvector(S0) && length(S0) > 1 if ~isequal([V 1],size(S0)) S0 = S0’ ; if ~isequal([V 1 ],size(S 0 )) error( ’S0 must be a column vector with the same number of rows as Y or a scalar’ ); end end elseif isscalar(S 0 ) S 0 = S 0 *ones(V, 1 ); else error( ’S0 must be a column vector or a scalar’ ); end % N = number of iterations if nargin < 8 n = 10000 ; else if ~(isscalar(n) && n > 0 && (round(n) == n) && isreal(n)) error( ’N must be a positive scalar integer’ ); end end if nargin < 9 rician = false ; % use rician noise distribution end if nargin < 10 prior = { ’flat’ , ’flat’ , ’flat’ , ’flat’ , ’reci’ }; % use flat prior distributions end % lim if rician limsz = [ 2 , 5 ]; else limsz = [ 2 , 4 ]; end if ~isequal(size(lim),limsz) error( ’lim must be 2x%d’ ,limsz( 2 )); end if rician && ~isequal(size(lim),[ 2 , 5 ]) error( ’noise variance must be estimated for Rician noise distribution’ ); end % burn- in steps if nargin < 11 burns = 1000 ; end % mean only if nargin < 12 meanonly = false ; end %%% %%% %%% %%% %%% %%% %%% %%% % % Parameter preparation % %%% %%% %%% %%% %%% %%% %%% %%% % M = length(f); nbs = length(b); if rician % startvalue for variance is2 = nbs./sum((Y-repmat(S 0 , 1 ,nbs).*(( 1 -repmat(f, 1 ,nbs)).*exp(-D*b) + ... repmat(f, 1 ,nbs).*exp(-(D+Dstar)*b))).^ 2 , 2 ); % 1 /s2 end if meanonly voxelsPerRun = V; else voxelsPerRun = min(V, 1500 ); % to avoid filling the memory for N > 40 , 000 end if meanonly thetasum = zeros(M, 4 +rician); theta2sum = zeros(M, 4 +rician); end % burn- in parameters burnUpdateInterval = 100 ; burnUpdateFraction = 1 / 2 ; %%% %%% %%% %%% %%% %%% %%% %%% % Parameter estimation % %%% %%% %%% %%% %%% %%% %%% %%% for i = 1 :ceil (M/voxelsPerRun) % partition voxels usedVoxels = (i- 1 )*voxelsPerRun + 1 :min (i*voxelsPerRun,M); fpart = f(usedVoxels); Dpart = D(usedVoxels); Dstarpart = Dstar(usedVoxels); S0part = S 0 (usedVoxels); if rician is2part = is2(usedVoxels); end Ypart = Y(usedVoxels, : ); Mpart = min(i*voxelsPerRun,M) - (i- 1 )*voxelsPerRun; % initialize parameter vector if meanonly theta = zeros(Mpart, 4 +rician, 2 ); else theta = zeros(Mpart, 4 +rician,n+burns); end theta( : , 1 : 4 , 1 ) = [fpart, Dpart, Dstarpart,S0part]; if rician theta( : , 5 , 1 ) = is2part; end % step length parameter w = zeros(Mpart, 4 +rician); w( : , 1 : 4 ) = [fpart Dpart Dstarpart S0part]/ 10 ; if rician w( : , 5 ) = 0 . 01 *ones(Mpart, 1 ); end N = zeros(Mpart, 4 +rician); % number of accepted samples % iterate for j = 2 , 3 ,...,n for j = 2 :n + burns % initialize theta(j) if meanonly theta( : , : , 2 ) = theta( : , : , 1 ); thetanew = theta( : , : , 2 ); thetaold = theta( : , : , 1 ); else theta( : , : ,j) = theta( : , : ,j- 1 ); thetanew = theta( : , : ,j); thetaold = theta( : , : ,j- 1 ); end % sample each parameter for k = 1 : 4 +rician % sample s and r and update s = thetaold( : ,k) + randn(Mpart, 1 ).*w( : ,k); r = rand(Mpart, 1 ); thetas = thetanew; thetas( : ,k) = s; alpha = acc_MH(thetas,thetanew,Ypart,b,lim,rician,prior{k}); sample_ok = r < alpha; thetanew(sample_ok,k) = thetas(sample_ok,k); thetanew(~sample_ok,k) = thetaold(~sample_ok,k); % reject samples N( : ,k) = N( : ,k) + sample_ok; end % prepare for next iteration if meanonly theta( : , : , 1 ) = thetanew; else theta( : , : ,j) = thetanew; end % save parameter value after burn- in phase if meanonly && j > burns thetasum = thetasum + thetanew; theta2sum = theta2sum + thetanew.^ 2 ; end % adapt step length if j <= burns*burnUpdateFraction && mod(j,burnUpdateInterval) == 0 w = w*(burnUpdateInterval+ 1 )./( 2 *((burnUpdateInterval+ 1 )-N)); N = zeros(Mpart, 4 +rician); end % Display iteration every 500 th iteration if ~mod(j, 500 ) && j > burns disp([ ’Iterations: ’ num2str(j-burns)]); elseif ~mod(j, 100 ) && j < burns disp([ ’Burn in-steps: ’ num2str(j)]); elseif j == burns disp([ ’Burn in complete: ’ num2str(j)]); end end % Saves distribution measures if meanonly %mean out.f.mean(usedVoxels) = thetasum( : , 1 )/n; out.D.mean(usedVoxels) = thetasum( : , 2 )/n; out.Dstar.mean(usedVoxels) = thetasum( : , 3 )/n; out.S 0 .mean(usedVoxels) = thetasum( : , 4 )/n; % standard deviation out.f.std(usedVoxels) = sqrt(theta2sum( : , 1 )/n-(thetasum( : , 1 )/n).^ 2 ); out.D.std(usedVoxels) = sqrt(theta2sum( : , 2 )/n-(thetasum( : , 2 )/n).^ 2 ); out.Dstar.std(usedVoxels) = sqrt(theta2sum( : , 3 )/n-(thetasum( : , 3 )/n).^ 2 ); out.S 0 .std(usedVoxels) = sqrt(theta2sum( : , 4 )/n-(thetasum( : , 4 )/n).^ 2 ); else %mean out.f.mean(usedVoxels) = mean(squeeze(theta( : , 1 ,burns + 1 :n+burns )), 2 ); out.D.mean(usedVoxels) = mean(squeeze(theta( : , 2 ,burns + 1 :n+burns )), 2 ); out.Dstar.mean(usedVoxels) = mean(squeeze(theta( : , 3 ,burns + 1 :n+burns )), 2 ); out.S 0 .mean(usedVoxels) = mean(squeeze(theta( : , 4 ,burns + 1 :n+burns )), 2 ); %median out.f.median(usedVoxels) = median(squeeze(theta( : , 1 ,burns + 1 :n+burns )), 2 ); out.D.median(usedVoxels) = median(squeeze(theta( : , 2 ,burns + 1 :n+burns )), 2 ); out.Dstar.median(usedVoxels) = median(squeeze(theta( : , 3 ,burns + 1 :n+burns )), 2 ); out.S 0 .median(usedVoxels) = median(squeeze(theta( : , 4 ,burns + 1 :n+burns )), 2 ); % mode out.f.mode(usedVoxels) = halfSampleMode(squeeze(theta( : , 1 ,burns + 1 :n+burns ))); out.D.mode(usedVoxels) = halfSampleMode(squeeze(theta( : , 2 ,burns + 1 :n+burns ))); out.Dstar.mode(usedVoxels) = halfSampleMode(squeeze(theta( : , 3 ,burns + 1 :n+burns ))); out.S 0 .mode(usedVoxels) = halfSampleMode(squeeze(theta( : , 4 ,burns + 1 :n+burns ))); % standard deviation out.f.std(usedVoxels) = std(squeeze(theta( : , 1 ,burns + 1 :n+burns )), 1 , 2 ); out.D.std(usedVoxels) = std(squeeze(theta( : , 2 ,burns + 1 :n+burns )), 1 , 2 ); out.Dstar.std(usedVoxels) = std(squeeze(theta( : , 3 ,burns + 1 :n+burns )), 1 , 2 ); out.S 0 .std(usedVoxels) = std(squeeze(theta( : , 4 ,burns + 1 :n+burns )), 1 , 2 ); end end function alpha = acc_MH(thetas,thetaj,Y,b,lim,rician,prior) % theta = [f, D, Dstar,S 0 , 1 /s2]; M = size(thetas, 1 ); N = length(b); q = zeros(M, 1 ); % p(theta |lim) pts = min((thetas >= repmat(lim(1,:),M,1)) & (thetas <= repmat(lim(2,:),M,1)),[],2); % D < D* pts = pts & (thetas(:,2) < thetas(:,3)); % signal model Ss = repmat(thetas(pts,4),1,N).*((1-repmat(thetas(pts,1),1,N)).*exp(-thetas(pts,2)*b) + repmat(thetas(pts,1),1,N).*exp(-thetas(pts,3)*b)); Sj = repmat(thetaj(pts,4),1,N).*((1-repmat(thetaj(pts,1),1,N)).*exp(-thetaj(pts,2)*b) + repmat(thetaj(pts,1),1,N).*exp(-thetaj(pts,3)*b)); ptsptj = double(pts); if strcmp(prior,’reci’) diffpar = find(thetas(1,:) ~= thetaj(1,:)); ptsptj(pts) = thetaj(pts,diffpar)./thetas(pts,diffpar); % rejects samples outside the limits % ~pts already == 0 elseif strcmp(prior,’lognorm’) diffpar = find(thetas(1,:) ~= thetaj(1,:)); % ptsptj = pts; if diffpar == 2 mu = -6; s = 1; elseif diffpar == 3 mu = -3.5; s = 1; else error(’lognorm prior not available’); % only for D and D* end ptsptj(pts) = lognormprior(thetas(pts,diffpar),mu,s)./lognormprior(thetaj(pts,diffpar),mu,s); % ~pts already == 0 elseif ~strcmp(prior,’flat’) error(’unknown prior’); end if rician % rician noise distribution q(pts) = exp(-N*log(thetaj(pts,5))+sum(Y(pts,:).^2,2).*0.5.*thetaj(pts,5)+sum(Sj.^2,2).*0.5.*thetaj(pts,5)-sum(logIo(Y(pts,:).*Sj.*repmat(thetaj(pts,5),1,N)),2)... +N*log(thetas(pts,5))-sum(Y(pts,:).^2,2).*0.5.*thetas(pts,5)-sum(Ss.^2,2).*0.5.*thetas(pts,5)+sum(logIo(Y(pts,:).*Ss.*repmat(thetas(pts,5),1,N)),2)) .* ptsptj(pts); else % gaussian noise distribution q(pts) = (sum((Y(pts,:)-Sj).^2,2)./sum((Y(pts,:)-Ss).^2,2)).^(N/2) .* ptsptj(pts); end alpha = min(1,q); function p = lognormprior(x,mu,s) p = 1./(s*sqrt(2*pi)*x).*exp(-(log(x)-mu).^2/(2*s^2)); function y = logIo(x) % Returns the natural log of the 0th order modified Bessel function of first kind for an argument x % Follows the exponential implementation of the Bessel function in Numerical Recipes, Ch. 6 % % Translated to MATLAB from C++ from the FSL source code and vectorized for % faster computations in MATLAB b = abs(x); y = zeros(size(x)); a1 = (x(b < 3.75)/3.75).^2; a2 = 3.75./b(b >= 3.75); y(b < 3.75) = log(1.0 + a1.*(3.5156229 + a1.*(3.0899424 + a1.*(1.2067492 + a1.*(0.2659732 + a1.*(0.0360768 + a1.*0.0045813)))))); y(b >= 3.75) = b(b >= 3.75) + log((0.39894228+a2.*(0.01328592+a2.*(0.00225319+a2.*(-0.00157565+a2.*(0.00916281+a2.*(-0.02057706+a2.*(0.02635537+a2.*(-0.01647633+a2.*0.00392377))))))))./sqrt(b(b>=3.75)));

function [pars,mask, gof] = IVIM_seg(Y,b,lim,blim,disp_prog) % function [pars,mask, gof] = IVIM_seg(Y,b,lim,blim,disp_prog) % % Function for monoexponential or stepwise biexponential fitting to % diffusion weighted MR data. The size of the input variable lim determines % the type of fit (monoexponential if size(lim) = [ 2 2 ], stepwise % biexponential if size(lim) = [ 2 4 ], estimation of D and f if % size(lim) = [ 2 3 ]) % % Input % - b is a column vector with B elements. It contains the b values % - Y is a BxV-matrix where V is the number of voxels. It contains the data % - lim is a 2 x2-, 2 x3 or 2 x4-matrix where the first row gives the lower % limit and the second row the upper limits of the model parameters. The % order of the parameters is [D,S 0 ], [D,S 0 ,f] or [D,S 0 ,f,D*] depending on % the fit % - blim is a scalar that determines the b-values used in the first of the % fits (b == blim is included) % - disp_prog is a scalar boolean. If it is set to "true" the progress of % the model fit is printed to the command window % % Output % - pars is a struct with fields for all model parameters % - mask is a struct with fields D and Dstar specifying the voxels that % converged to a solution within the specified limits % - gof is a struct with fields SSE (sum of squared errors) and R2 % (coefficient of determination) % % By Oscar Jalnefjord 2018 -08- 27 % % If you use this function in research, please cite: % Jalnefjord et al. 2018 Comparison of methods for estimation of the % intravoxel incoherent motion (IVIM) diffusion coefficient (D) and % perfusion fraction (f), MAGMA %%% %%% %%% %%% %%% %%% %%% %%% %%% % % Error handling %%% %%% %%% %%% %%% %%% %%% %%% %%% % % b if ~iscolumn(b) if isrow(b) b = b ’; else error(’ b must be a vector ’); end end % Y if size(Y,1) ~= length(b) if size(Y,2) == length(b) Y = Y’ ; else error( ’Dimensions of Y and b must agree’ ); end end % lim if isequal(size(lim),[ 2 , 2 ]) % Limits [D,S 0 ] step2 = false ; % fit monoexponential model estf = false ; elseif isequal(size(lim),[ 2 , 3 ]) % Limits [D,S 0 ,f] step2 = false ; estf = true ; % calculate f if ~any(b== 0 ) error( ’b=0 is need to calculate f’ ); end elseif isequal(size(lim),[ 2 , 4 ]) % Limits [D,S 0 ,f,D*] step2 = true ; % fit biexponential model estf = true ; else error( ’lim must be 2x2, 2x3 or 2x4’ ); end % blim if nargin < 4 if step2 == false && estf == false blim = 0 ; else blim = 200 ; end end % Display progress if nargin < 5 disp_prog = true ; else if ~isscalar(disp_prog) || ~islogical(disp_prog) error( ’disp_prog must be a scalar boolean’ ); end end if disp_prog fprintf( ’ Starting optization process ’ ); fprintf( ’--------------------------- ’ ); end % Common parameters optlim = 1 e- 6 ; n = size(Y, 2 ); %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% % % Fitting D and A (S 0 ) % %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% if disp_prog fprintf( ’Fitting D and A ’ ); end if estf % prepare variables for fitting at high b-values Yred = Y(b>=blim, : ); bred = b(b>=blim); if length(bred) < 2 error( ’blim is too large. At least two b-value must be >=blim’ ) end else % monoexponential fit Yred = Y; bred = b; end % Mask to remove background mask_back = sum(Y, 1 ) > 0 ; if disp_prog fprintf( ’Discarding %d voxels as background ’ ,sum(~mask_back)); end % Allocation D = zeros( 1 ,n); maskD = mask_back; % Optimization [D(mask_back),maskD(mask_back)] = optimizeD(Yred( : ,mask_back),bred,optlim,lim( : , 1 ),disp_prog); % Calculates A based on D A = sum(Yred.*exp(-bred*D))./sum(exp(- 2 *bred*D)); % Assures that A is within limits D(A < lim( 1 , 2 )) = lim( 1 , 1 ); D(A > lim( 2 , 2 )) = lim( 2 , 1 ); maskD(A < lim( 1 , 2 ) | A > lim(2,2)) = false ; A(A < lim(1,2)) = lim(1,2); A(A > lim(2,2)) = lim(2,2); % Calculate f if estf && ~step2 if disp_prog fprintf(’ Calculating f ’); end f = 1-A./mean(Y(b==0,:),1); f(f < lim(1,3)) = lim(1,3); f(f > lim(2,3)) = lim(2,3); S0 = mean(Y(b==0,:),1); end if step2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Fitting of D* and f ( and S0) % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if disp_prog fprintf(’ Fitting D* and f ’); end % calculation of f given D* and the previously estimated parameters (D and A) fcalc = @(Dstar,mask) sum(exp(-b*Dstar).*(Y(:,mask)-repmat(A(mask),length(b),1).*exp(-b*D(mask))))./sum(exp(-b*Dstar).*(Y(:,mask)-repmat(A(mask),length(b),1).*(exp(-b*D(mask))-exp(-b*Dstar)))); % calculate possible range of f to remove voxels where f cannot % result in a fit within limits f1 = zeros(1,n); f1(mask_back) = fcalc(lim(1,4)*ones(1,sum(mask_back)),mask_back); f2 = zeros(1,n); f2(mask_back) = fcalc(lim(2,4)*ones(1,sum(mask_back)),mask_back); maskf = maskD & (((f1 > lim(1,3)) & (f1 < lim(2,3))) | ((f2 > lim( 1 , 3 )) & (f2 < lim( 2 , 3 )))); if disp_prog fprintf( ’Discarding %d voxels due to f out of bounds ’ ,sum(maskD) - sum(maskf)); end % Prepares output Dstar = zeros( 1 ,n); Dstar(f1 < lim( 1 , 3 )) = lim( 1 , 4 ); Dstar(f2 > lim( 2 , 3 )) = lim( 2 , 4 ); maskDstar = maskf; % Optimization [Dstar(maskf),maskDstar(maskf)] = optimizeD(Y( : ,maskf)-repmat(A(maskf),length(b), 1 ).*exp(-b*D(maskf)),b,optlim,lim( : , 4 ),disp_prog); % Calculates f f = lim( 1 , 1 )*ones( 1 ,n); f(f1 < lim( 1 , 3 )) = lim( 1 , 3 ); f(f2 > lim( 2 , 3 )) = lim( 2 , 3 ); f(maskDstar) = fcalc(Dstar(maskDstar),maskDstar); % Checks for f out of bounds maskf = maskf & (f > lim( 1 , 3 )) & (f < lim( 2 , 3 )); Dstar(f < lim( 1 , 3 )) = lim( 1 , 4 ); Dstar(f > lim( 2 , 3 )) = lim( 2 , 4 ); f(f < lim( 1 , 3 )) = lim( 1 , 3 ); f(f > lim( 2 , 3 )) = lim( 2 , 3 ); maskDstar = maskDstar & maskf; % Calculate S 0 S 0 = Y( 1 , : ); S 0 (f < . 9 ) = A(f < . 9 )./( 1 - f(f < . 9 )); % risk for Inf if f = 1 is used end if nargout == 3 nbs = length(b); if step2 gof.SSE = sum((Y - repmat(S 0 ,nbs, 1 ).*(( 1 -repmat(f,nbs, 1 )).*exp(-b*D)+repmat(f,nbs, 1 ).*exp(-b*(D+Dstar)))).^ 2 , 1 ); else gof.SSE = sum((Y - repmat(A,nbs, 1 ).*exp(-b*D)).^ 2 , 1 ); end SStot = sum((Y-repmat(mean(Y, 1 ),nbs, 1 )).^ 2 , 1 ); gof.R2 = 1 -gof.SSE./SStot; % R2 can be negative end % Saves output variables pars.D = D; if step2 pars.A = A; end mask.D = maskD; if estf pars.f = f; pars.S 0 = S 0 ; if step2 pars.Dstar = Dstar; mask.Dstar = maskDstar; end else pars.S 0 = A; end if disp_prog fprintf( ’ Done! ’ ); end function [D,mask_background] = optimizeD(Y,b,optlim,Dlim,disp_prog) % Prepares variables n = size(Y, 2 ); D = zeros( 1 ,n); yb = Y .* repmat(b, 1 ,n); % precalculated for speed %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% % Control if a minimum is within the interval %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% % checks that all diff < 0 for Dlow Dlow = Dlim( 1 )*ones( 1 ,n); difflow = Ddiff(Y,yb,b,Dlow,n); low_check = difflow < 0 ; % difflow must be < 0 if the optimum is within the interval % checks that all diff > 0 for Dhigh Dhigh = Dlim( 2 )*ones( 1 ,n); diffhigh = Ddiff(Y,yb,b,Dhigh,n); high_check = diffhigh > 0 ; % diffhigh must be > 0 if the optimum if within the interval % sets parameter value with optimum out of bounds D(~low_check) = Dlim( 1 ); % difflow > 0 means that the mimimum has been passed D(~high_check) = Dlim( 2 ); % diffhigh < 0 means that the minium is beyond the interval % Only the voxels with a possible minimum should be optimized mask = low_check & high_check; if disp_prog fprintf( ’Discarding %d voxels due to parameters out of bounds ’ ,sum(~mask)); end % Saves a mask to know which voxels that has been optimized mask_background = mask; % Allocates all variables sz = size(Dlow); D_lin = zeros(sz); diff_lin = zeros(sz); D_mid = zeros(sz); diff_mid = zeros(sz); q1_lin = zeros(sz); q2_lin = zeros(sz); q1_mid = zeros(sz); q2_mid = zeros(sz); %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% % Iterative method for finding the point where diff = 0 % %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% k = 0 ; while any(mask) % Continues is there are voxels left to optimize % Assumes diff is linear within the search interval [Dlow Dhigh] D_lin(mask) = Dlow(mask) - difflow(mask) .* (Dhigh(mask) - Dlow(mask)) ./ (diffhigh(mask) - difflow(mask)); % Calculates diff in the point of intersection given by the previous expression [diff_lin(mask),q1_lin(mask),q2_lin(mask)] = Ddiff(Y( : ,mask),yb( : ,mask),b,D_lin(mask),sum(mask)); % As a protential speed up, the mean of Dlow and Dhigh is also calculated D_mid(mask) = (Dlow(mask) + Dhigh(mask))/ 2 ; [diff_mid(mask),q1_mid(mask),q2_mid(mask)] = Ddiff(Y( : ,mask),yb( : ,mask),b,D_mid(mask),sum(mask)); % If diff < 0 , then the point of intersection or mean is used as the % new Dlow. Only voxels with diff < 0 are updated at this step. Linear % interpolation or the mean is used depending of which method that % gives the smallest diff updatelow_lin = (diff_lin < 0 ) & ((diff_mid > 0 ) | ((D_lin > D_mid) & (diff_mid < 0))); updatelow_mid = (diff_mid < 0) & ((diff_lin > 0) | ((D_mid > D_lin) & (diff_lin < 0 ))); Dlow(updatelow_lin) = D_lin(updatelow_lin); Dlow(updatelow_mid) = D_mid(updatelow_mid); % If diff > 0 , then the point of intersection or mean is used as the % new Dhigh. Only voxels with diff > 0 are updated at this step. % Om diff 锟絩 > 0 ska sk锟絩ningspunkten anv锟絥das som ny Dhigh. Linear % interpolation or the mean is used depending of which method that % gives the smallest diff updatehigh_lin = (diff_lin > 0 ) & ((diff_mid < 0 ) | ((D_lin < D_mid) & (diff_mid > 0))); updatehigh_mid = (diff_mid > 0) & ((diff_lin < 0) | ((D_mid < D_lin) & (diff_lin > 0 ))); Dhigh(updatehigh_lin) = D_lin(updatehigh_lin); Dhigh(updatehigh_mid) = D_mid(updatehigh_mid); % Updates the mask to exclude voxels that fulfills the optimization % limit from the mask opt_lin = abs( 1 -q1_lin./q2_lin) < optlim; opt_mid = abs( 1 -q1_mid./q2_mid) < optlim; D(opt_lin) = D_lin(opt_lin); D(opt_mid) = D_mid(opt_mid); % Not optimal if both D_lin and D_mean fulfills the optimization limit, but has a small impact on the result as long as optlim is small % Updates the mask mask = mask & (~(opt_lin | opt_mid)); % Calculates diff for the new bounds if any(mask) difflow(mask) = Ddiff(Y(:,mask),yb(:,mask),b,Dlow(mask),sum(mask)); diffhigh(mask) = Ddiff(Y(:,mask),yb(:,mask),b,Dhigh(mask),sum(mask)); end k = k + 1; if disp_prog fprintf(’Iteration %2d: %6d voxels left ’,k,sum(mask)) end end function [diff,q1,q2] = Ddiff(Y,yb,b,D,n) q1 = sum(exp(-2*b*D),1) .* sum(yb.*exp(-b*D),1); q2 = sum(Y.*exp(-b*D),1).*sum(repmat(b,1,n).*exp(-2*b*D),1); diff = q1 - q2;

function maps = IVIMmodelfit(im,bvals,fittype,roi,opts) % function maps = IVIMmodelfit(im,bvals,fittype,roi,opts) % maps = IVIMmodelfit(im,bvals,fittype,roi,opts) % IVIMmodelfit(im,bvals) % % Wrapper for monoexponential or stepwise biexponential least-squares % fitting, and biexponential bayesian fitting to diffusion weighted MR % data. % % Input % - im is a 4 D matrix with b-values representing the 4 th dimension, or % the path to a nifti file containing the corresponding data. It contains % the diffusion weighted image data % - bvals is vector or a path to an FSL-style (space separated) text file % containing the b-values used to obtain the data in im % - fittype is a string the determines which type of fitting that is used. % "seg" ( default ) gives least-squares fitting while "bayes" gives % bayesian fitting % - roi is a 3 D matrix matching im or a path to a nifti file with the % corresponding data. It is use to mask out data of interest from im % (e.g. a brain mask or a tumor roi). This is recommended for % computational speed. roi = [] or roi = ’’ corresponds to no mask ( i.e. % fitting in all voxels ) % - opts is a struct used to set additional options ( suitable default % values are provided for all these ). % * Common: lim ( parameter limits in a 2 x4 matrix with order [D,S0,f,D*], % see IVIM_seg for advanced use ), outfile ( name of nifti file % containing the parameter maps ) % * seg: bthr, dispprog ( see IVIM_seg ) % * bayes: its, burns, rician, meanonly, prior ( see IVIM_bayes ) % % - blim is a scalar that determines the b-values used in the first of the % fits ( b == blim is included ) % - disp_prog is a scalar boolean. If it is set to " true " the progress of % the model fit is printed to the command window % % Output % - maps is a 4D matrix containing the estimated parameter maps. If im % points out a nifti file, the results are also written to a nifti file % % Example % % im = "IVIMdata.nii.gz" ; % mask = "IVIMdata_mask.nii.gz" ; % bvals = "bvals.txt" ; % fittype = "seg" ; % opts.bthr = 200 ; % opts.dispprog = false ; % opts.outfile = "IVIMmaps" ; % IVIMmodelfit(im,bvals,fittype,mask,opts); % % % By Oscar Jalnefjord 2020 -08 -04 % % If you use this function in research, please cite ref 1 if you use % least-squares fitting and ref 2 if you use bayesian fitting: % [ 1 ] Jalnefjord et al. 2018 Comparison of methods for estimation of the % intravoxel incoherent motion ( IVIM ) diffusion coefficient ( D ) and % perfusion fraction ( f ), MAGMA % [2] Gustafsson et al. 2017 Impact of prior distributions and central % tendency measures on Bayesian intravoxel incoherent motion model % fitting, MRM %%%%%%%%%%%%%%%%%%%%% % First input check % %%%%%%%%%%%%%%%%%%%%% fittypes = { ’seg’ , ’bayes’ }; if nargin < 3 fittype = fittypes{ 1 }; else if isstring ( fittype ) fittype = char (fittype); end if ischar ( fittype ) if ~ any ( strcmp(fittype,fittypes )) error ( ’Unknown fit type: %s’ ,fittype ) ; end else error ( ’fittype must be a character string’ ) ; end end %%%%%%%%%%%%%% % Check opts % %%%%%%%%%%%%%% options = { ’bthr’ , ’lim’ , ’dispprog’ , ’its’ , ’burns’ ,... ’rician’ , ’meanonly’ , ’prior’ , ’outfile’ }; if isstruct ( opts ) opts_fields = fieldnames(opts); for i = 1 :length(opts_fields) if ~any(strcmp(opts_fields{i},options)) error( ’Unknown field "%s" in struct opts’ ,opts_fields{i}); end end else error ( ’opts must be a struct’ ) ; end %%%%%%%%%%%%%%%%%%% % Read from files % %%%%%%%%%%%%%%%%%%% w2f = false ; % write results to file if isstring ( im ) im = char (im); end if ischar ( im ) w2f = true ; % read b- values if isstring ( bvals ) bvals = char (bvals); end if ischar ( bvals ) try fid = fopen(bvals); b = fscanf(fid, ’%f ’ ); fclose(fid); catch error ( ’Unable to read b-values from file: %s’ ,bvals ) ; end else error ( ’bvals must be a character string’ ) ; end % read image file imfile = im; try info = niftiinfo(im); im = niftiread(info); catch error ( ’Unable to read image from file: %s’ ,imfile ) ; end % read roi file if nargin > 3 roifile = roi; try roi = logical(niftiread(roifile)); catch error ( ’Unable to read roi from file: %s’ ,roifile ) ; end end else % should be array if ~isnumeric(im) error( ’im must be a character string or a numeric array’ ); end if isnumeric ( bvals ) b = bvals; else error( ’b must be a numeric vector if im is’ ); end if nargin < 3 if ~isnumeric(roi) error( ’roi must be a numeric array if im is’ ); end end end %%%%%%%%%%%%%%%%%%%%%% % Check input arrays % %%%%%%%%%%%%%%%%%%%%%% % im should be 4 D array d = size(im); if length ( d ) ~ = 4 error( ’The image must have 4 dimensions’ ); end % b should be vector with length equal to 4 th dimension of im if isvector ( b ) if length ( b ) ~ = d( 4 ) error( ’Number of b-values in b must equal the size of the 4th dimension of the image’ ); end else error ( ’The b-values must be contained in a vector’ ) ; end % roi should be if nargin < 4 roi = true (d( 1 : 3 )); else if ~isequal(size(roi),d( 1 : 3 )) error( ’The ROI must have the same size as the first three dimensions of the image’ ); end end %%%%%%%%%%%%%%%%% % Run model fit % %%%%%%%%%%%%%%%%% % Turn image into data matrix Y = im2Y(im,roi); % Define which parameters to estimate pars = { ’D’ , ’S0’ , ’f’ , ’Dstar’ }; lim = [ 0 0 0 0 ; 3e-3 2 *max(Y(:)) 1 1 ]; blim = 200 ; dispprog = true ; n = 10000 ; burns = 2000 ; rician = false ; meanonly = false ; prior = { ’flat’ , ’lognorm’ , ’lognorm’ , ’flat’ , ’reci’ }; fstart = 0.1 *ones(size(Y(:, 1 ))); Dstart = 1e-3 *ones(size(Y(:, 1 ))); Dstarstart = 20e-3 *ones(size(Y(:, 1 ))); % Run model fit if any ( strcmp( ’lim’ ,opts_fields )) lim = opts.lim; end switch fittype case fittypes{ 1 } % segmented if any ( strcmp( ’bthr’ ,opts_fields )) blim = opts.bthr; end if any ( strcmp( ’dispprog’ ,opts_fields )) dispprog = opts.dispprog; end if ( size(lim, 1 ) ~ = 2 ) || (size(lim, 2 ) < 2 ) || (size(lim, 2 ) > 4 ) error( ’lim must be of size 2xP where P is in range 2-4’ ); end res = IVIM_seg(Y,b,lim,blim,dispprog); case fittypes{ 2 } % Bayesian if any ( strcmp( ’its’ ,opts_fields )) n = opts.its; end if any ( strcmp( ’burns’ ,opts_fields )) burns = opts.burns; end if any ( strcmp( ’rician’ ,opts_fields )) rician = opts.rician; end if any ( strcmp( ’meanonly’ ,opts_fields )) meanonly = opts.meanonly; end if any ( strcmp( ’prior’ ,opts_fields )) prior = opts.prior; end if ~isequal(size(lim),[ 2 , 4 ]) error( ’lim must be of size [2,4] for Bayesian fitting’ ); end res = IVIM_bayes(Y,fstart,Dstart,Dstarstart,mean(Y(:,b==min(b)), 2 ),b,lim(:,[ 3 1 4 2 ]),n,rician,prior,burns,meanonly); end %%%%%%%%%%%%%%%%%% % Prepare output % %%%%%%%%%%%%%%%%%% pars = pars( 1 :size(lim, 2 )); maps = nan([size(roi) length(pars)]); for i = 1 :length(pars) temp = maps(:,:,:,i); switch fittype case fittypes{ 1 } % Segmented temp ( roi ) = res.(pars{i}); case fittypes{ 2 } % Bayesian if meanonly temp ( roi ) = res.(pars{i}).mean; else temp(roi) = res.(pars{i}).mode; end end maps ( :,:,:,i ) = temp; end if w2f % check if output if strcmp ( imfile(end -2 :end ),’.gz’) comp = true ; else comp = false ; end info.Description = ’IVIM parameter maps’ ; info.ImageSize( 4 ) = length(pars); outfile = ’IVIMmaps’ ; if any ( strcmp( ’outfile’ ,opts_fields )) outfile = opts.outfile; end try niftiwrite ( maps,outfile,info, ’Compressed’ ,comp ) ; catch error ( ’Failed to write parameter maps to file: %s’ ,outfile ) ; end end

function y = logIo(x) % Returns the natural log of the 0th order modified Bessel function of first kind for an argument x % Follows the exponential implementation of the Bessel function in Numerical Recipes, Ch. 6 % % Translated to MATLAB from C++ from the FSL source code and vectorized for % faster computations in MATLAB b = abs(x); % if b < 3.75 % a = x/3.75; % a = a.^2; % % Bessel function evaluation % y = 1.0+a*(3.5156229+a*(3.0899424+a*(1.2067492+a*(0.2659732+a*(0.0360768+a*0.0045813))))); % y = log (y); % else % a = 3.75/b; % %Bessel function evaluation % %y=(exp(b)/sqrt(b))*(0.39894228+a*(0.01328592+a*(0.00225319+a*(-0.00157565+a*(0.00916281+a*(-0.02057706+a*(0.02635537+a*(-0.01647633+a*0.00392377)))))))); % %Logarithm of Bessel function % y=b+ log ((0.39894228+a*(0.01328592+a*(0.00225319+a*(-0.00157565+a*(0.00916281+a*(-0.02057706+a*(0.02635537+a*(-0.01647633+a*0.00392377))))))))/sqrt(b)); % end y = zeros(size(x)); a1 = (x(b < 3.75)/3.75).^2; a2 = 3.75./b(b >= 3.75); y(b < 3.75) = log(1.0 + a1.*(3.5156229 + a1.*(3.0899424 + a1.*(1.2067492 + a1.*(0.2659732 + a1.*(0.0360768 + a1.*0.0045813)))))); y(b >= 3.75) = b(b >= 3.75) + log((0.39894228+a2.*(0.01328592+a2.*(0.00225319+a2.*(-0.00157565+a2.*(0.00916281+a2.*(-0.02057706+a2.*(0.02635537+a2.*(-0.01647633+a2.*0.00392377))))))))./sqrt(b(b>=3.75)));