Homec4science
Diffusion iCAPs public (master) Browse

iCAPs/packages/spm12/spm_dartel_integrate.m
31022b3a13e6master

Last Change
Raw File

spm_dartel_integrate.m
View Options

function [Phi,DPhi] = spm_dartel_integrate(U,t,K)
% Integrate a Dartel flow field
% FORMAT [Phi,DPhi] = spm_dartel_exp(U,t,K)
% U - name of flow field (nx x ny x nz x nt x 3)
% t - [t0 t1] Start and end time (values between 0 and 1)
% K - log2 of the Euler time steps to integrate the
% flow field.
%
% Phi - deformation field (nx x ny x nz x 3)
% DPhi - Jacobian determinant field (nx x ny x nz)
%
% The function integrates
% Phi(x,t) = \int_{t_0}^{t_1} U(Phi(x,t),t) dt
% where U is a piecewise constant flow field
%
% Note: this function is ready for LDDMM-style flow fields, even
% though the none of the official Dartel tools can generate them
% yet.
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
% John Ashburner
% $Id: spm_dartel_integrate.m 5506 2013-05-14 17:13:43Z john $
if isa(U,'char'), U = nifti(U); end;
if isa(U,'nifti'), U = U.dat; end;
% Figure out which bits of flow field to use, the numbers of
% Euler time steps, and the scales.
nt = size(U,4);
K = ceil(log2(2^K/nt));
big = 2^12;
t = min(max(t,0),1);
t0 = t(1);
t1 = t(2);
tt0 = t0*nt;
tt1 = t1*nt;
sc = zeros(1,nt);
if t1>t0,
for i=1:nt,
sc(i) = max(round((min(tt1,i)-max(tt0,i-1))*big)/big,0);
end
ind = find(sc~=0);
else
for i=1:nt,
sc(i) = max(round((min(tt0,i)-max(tt1,i-1))*big)/big,0);
end
ind = fliplr(find(sc~=0));
end
if isempty(ind), % t0==t1
ind = 1;
sc = 0;
ts = 0;
else
sc = sc(ind);
ts = zeros(1,numel(ind));
for i=1:numel(ind),
ts(i) = ceil(log2((2^K)*sc(i))-1/big);
sc(i) = sc(i)*2^(K-ts(i));
end
if t0>t1, sc = -sc; end
end
% Do the integrations
if nargout==1,
% Deformation field only
u = squeeze(single(U(:,:,:,ind(1),:)));
Phi = dartel3('Exp',u,[ts(1), sc(1)]);
for i=2:numel(ind),
u = squeeze(single(U(:,:,:,ind(i),:)));
Phi = dartel3('comp',Phi,dartel3('Exp',u,[ts(i), sc(i)]));
end
else
% Deformation and Jacobian determinant fields
u = squeeze(single(U(:,:,:,ind(1),:)));
[Phi,DPhi] = dartel3('Exp',u,[ts(1), sc(1), 1]);
for i=2:numel(ind),
u = squeeze(single(U(:,:,:,ind(i),:)));
[Phi1,DPhi1] = dartel3('Exp',u,[ts(i), sc(i), 1]);
[Phi,DPhi] = dartel3('comp',Phi,Phi1,DPhi,DPhi1);
clear Phi1 DPhi1
end
end
c4science · Help