Homec4science
Diffusion Slepians (master) Browse

Slepi/Code/Utilities/spm12/spm_uw_apply.m
8a186cd82567master

Last Change
Raw File

spm_uw_apply.m
View Options

function varargout = spm_uw_apply(ds,flags)
% Reslices images volume by volume
% FORMAT spm_uw_apply(ds,[flags])
% or
% FORMAT P = spm_uw_apply(ds,[flags])
%
%
% ds - a structure created by spm_uw_estimate.m containing the fields:
% ds can also be an array of structures, each struct corresponding
% to one sesssion (it hardly makes sense to try and pool fields across
% sessions since there will have been a reshimming). In that case each
% session is unwarped separately, unwarped into the distortion space of
% the average (default) position of that series, and with the first
% scan on the series defining the pahse encode direction. After that each
% scan is transformed into the space of the first scan of the first series.
% Naturally, there is still only one actual resampling (interpolation).
% It will be assumed that the same unwarping parameters have been used
% for all sessions (anything else would be truly daft).
%
% .P - Images used when estimating deformation field and/or
% its derivative w.r.t. modelled factors. Note that this
% struct-array may contain .mat fields that differ from
% those you would observe with spm_vol(P(1).fname). This
% is because spm_uw_estimate has an option to re-estimate
% the movement parameters. The re-estimated parameters are
% not written to disc (in the form of .mat files), but rather
% stored in the P array in the ds struct.
%
% .order - Number of basis functions to use for each dimension.
% If the third dimension is left out, the order for
% that dimension is calculated to yield a roughly
% equal spatial cut-off in all directions.
% Default: [8 8 *]
% .sfP - Static field supplied by the user. It should be a
% filename or handle to a voxel-displacement map in
% the same space as the first EPI image of the time-
% series. If using the FieldMap toolbox, realignment
% should (if necessary) have been performed as part of
% the process of creating the VDM. Note also that the
% VDM mut be in undistorted space, i.e. if it is
% calculated from an EPI based field-map sequence
% it should have been inverted before passing it to
% spm_uw_estimate. Again, the FieldMap toolbox will
% do this for you.
% .regorder - Regularisation of derivative fields is based on the
% regorder'th (spatial) derivative of the field.
% Default: 1
% .lambda - Fudge factor used to decide relative weights of
% data and regularisation.
% Default: 1e5
% .jm - Jacobian Modulation. If set, intensity (Jacobian)
% deformations are included in the model. If zero,
% intensity deformations are not considered.
% .fot - List of indexes for first order terms to model
% derivatives for. Order of parameters as defined
% by spm_imatrix.
% Default: [4 5]
% .sot - List of second order terms to model second
% derivatives of. Should be an nx2 matrix where
% e.g. [4 4; 4 5; 5 5] means that second partial
% derivatives of rotation around x- and y-axis
% should be modelled.
% Default: []
% .fwhm - FWHM (mm) of smoothing filter applied to images prior
% to estimation of deformation fields.
% Default: 6
% .rem - Re-Estimation of Movement parameters. Set to unity means
% that movement-parameters should be re-estimated at each
% iteration.
% Default: 0
% .noi - Maximum number of Iterations.
% Default: 5
% .exp_round - Point in position space to do Taylor expansion around.
% 'First', 'Last' or 'Average'.
% .p0 - Average position vector (three translations in mm
% and three rotations in degrees) of scans in P.
% .q - Deviations from mean position vector of modelled
% effects. Corresponds to deviations (and deviations
% squared) of a Taylor expansion of deformation fields.
% .beta - Coeffeicents of DCT basis functions for partial
% derivatives of deformation fields w.r.t. modelled
% effects. Scaled such that resulting deformation
% fields have units mm^-1 or deg^-1 (and squares
% thereof).
% .SS - Sum of squared errors for each iteration.
%
% flags - a structure containing various options. The fields are:
%
% mask - mask output images (1 for yes, 0 for no)
% To avoid artifactual movement-related variance the realigned
% set of images can be internally masked, within the set (i.e.
% if any image has a zero value at a voxel than all images have
% zero values at that voxel). Zero values occur when regions
% 'outside' the image are moved 'inside' the image during
% realignment.
%
% mean - write mean image
% The average of all the realigned scans is written to
% mean*.<ext>.
%
% interp - the interpolation method (see e.g. spm_bsplins.m).
%
% which - Values of 0 or 1 are allowed.
% 0 - don't create any resliced images.
% Useful if you only want a mean resliced image.
% 1 - reslice all the images.
%
% udc - Values 1 or 2 are allowed
% 1 - Do only unwarping (not correcting
% for changing sampling density).
% 2 - Do both unwarping and Jacobian correction.
%
%
% prefix - Filename prefix for resliced image files. Defaults to 'u'.
%
% The spatially realigned images are written to the orginal
% subdirectory with the same filename but prefixed with an 'u'.
% They are all aligned with the first.
%__________________________________________________________________________
% Copyright (C) 2003-2015 Wellcome Trust Centre for Neuroimaging
% Jesper Andersson
% $Id: spm_uw_apply.m 6301 2015-01-12 17:23:08Z guillaume $
SVNid = '$Rev: 6301 $';
%-Say hello
%--------------------------------------------------------------------------
SPMid = spm('FnBanner',mfilename,SVNid);
%-Parameters
%--------------------------------------------------------------------------
tiny = 5e-2;
def_flags = spm_get_defaults('realign.write');
def_flags.udc = 1;
def_flags.prefix = 'u';
defnames = fieldnames(def_flags);
if nargin < 1 || isempty(ds)
ds = load(spm_select(1,'.*uw\.mat$','Select Unwarp result file'),'ds');
ds = ds.ds;
end
%
% Default to using Jacobian modulation for the reslicing if it was
% used during the estimation phase.
%
if ds(1).jm ~= 0
def_flags.udc = 2;
end
%
% Replace defaults with user supplied values for all fields
% defined by user.
%
if nargin < 2 || isempty(flags)
flags = def_flags;
end
for i=1:length(defnames)
if ~isfield(flags,defnames{i})
flags.(defnames{i}) = def_flags.(defnames{i});
end
end
if numel(flags.which) == 2
flags.mean = flags.which(2);
flags.which = flags.which(1);
end
ntot = 0;
for i=1:length(ds)
ntot = ntot + length(ds(i).P);
end
hold = [repmat(flags.interp,1,3) flags.wrap];
%
% Create empty sfield for all structs.
%
[ds.sfield] = deal([]);
%
% Make space for output P-structs if required
%
if nargout > 0
oP = cell(length(ds),1);
end
%-Create mask, if required
%--------------------------------------------------------------------------
if flags.mask || flags.mean,
fprintf('%-40s: %30s','Computing mask...',''); %-#
spm_progress_bar('Init',ntot,'Computing available voxels',...
'volumes completed');
[x,y,z] = ndgrid(1:ds(1).P(1).dim(1),1:ds(1).P(1).dim(2),1:ds(1).P(1).dim(3));
xyz = [x(:) y(:) z(:) ones(prod(ds(1).P(1).dim(1:3)),1)]; clear x y z;
if flags.mean
Count = zeros(prod(ds(1).P(1).dim(1:3)),1);
Integral = zeros(prod(ds(1).P(1).dim(1:3)),1);
end
% if flags.mask
msk = zeros(prod(ds(1).P(1).dim(1:3)),1);
% end
tv = 1;
for s=1:length(ds)
def_array = zeros(prod(ds(s).P(1).dim(1:3)),size(ds(s).beta,2));
Bx = spm_dctmtx(ds(s).P(1).dim(1),ds(s).order(1));
By = spm_dctmtx(ds(s).P(1).dim(2),ds(s).order(2));
Bz = spm_dctmtx(ds(s).P(1).dim(3),ds(s).order(3));
if isfield(ds(s),'sfP') && ~isempty(ds(s).sfP)
T = ds(s).sfP.mat\ds(1).P(1).mat;
txyz = xyz * T';
c = spm_bsplinc(ds(s).sfP,ds(s).hold);
ds(s).sfield = spm_bsplins(c,txyz(:,1),txyz(:,2),txyz(:,3),ds(s).hold);
ds(s).sfield = ds(s).sfield(:);
clear c txyz;
end
for i=1:size(ds(s).beta,2)
def_array(:,i) = spm_get_def(Bx,By,Bz,ds(s).beta(:,i));
end
sess_msk = zeros(prod(ds(1).P(1).dim(1:3)),1);
for i = 1:numel(ds(s).P)
T = inv(ds(s).P(i).mat) * ds(1).P(1).mat;
txyz = xyz * T';
txyz(:,2) = txyz(:,2) + spm_get_image_def(ds(s).P(i),ds(s),def_array);
tmp = false(size(txyz,1),1);
if ~flags.wrap(1), tmp = tmp | txyz(:,1) < (1-tiny) | txyz(:,1) > (ds(s).P(i).dim(1)+tiny); end
if ~flags.wrap(2), tmp = tmp | txyz(:,2) < (1-tiny) | txyz(:,2) > (ds(s).P(i).dim(2)+tiny); end
if ~flags.wrap(3), tmp = tmp | txyz(:,3) < (1-tiny) | txyz(:,3) > (ds(s).P(i).dim(3)+tiny); end
sess_msk = sess_msk + real(tmp);
spm_progress_bar('Set',tv);
tv = tv+1;
end
msk = msk + sess_msk;
if flags.mean, Count = Count + repmat(length(ds(s).P),prod(ds(s).P(1).dim(1:3)),1) - sess_msk; end % Changed 23/3-05
%
% Include static field in estmation of mask.
%
if isfield(ds(s),'sfP') && ~isempty(ds(s).sfP)
T = inv(ds(s).sfP.mat) * ds(1).P(1).mat;
txyz = xyz * T';
tmp = false(size(txyz,1),1);
if ~flags.wrap(1), tmp = tmp | txyz(:,1) < (1-tiny) | txyz(:,1) > (ds(s).sfP.dim(1)+tiny); end
if ~flags.wrap(2), tmp = tmp | txyz(:,2) < (1-tiny) | txyz(:,2) > (ds(s).sfP.dim(2)+tiny); end
if ~flags.wrap(3), tmp = tmp | txyz(:,3) < (1-tiny) | txyz(:,3) > (ds(s).sfP.dim(3)+tiny); end
msk = msk + real(tmp);
end
if isfield(ds(s),'sfield') && ~isempty(ds(s).sfield)
ds(s).sfield = [];
end
end
if flags.mask, msk = find(msk ~= 0); end
fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...done'); %-#
end
%-Reslicing images
%--------------------------------------------------------------------------
fprintf('%-40s: %30s','Reslicing images...',''); %-#
spm_progress_bar('Init',ntot,'Reslicing','volumes completed');
jP = ds(1).P(1);
jP = rmfield(jP,{'fname','descrip','n','private'});
jP.dim = jP.dim(1:3);
jP.dt = [spm_type('float64'), spm_platform('bigend')];
jP.pinfo = [1 0]';
tv = 1;
for s=1:length(ds)
def_array = zeros(prod(ds(s).P(1).dim(1:3)),size(ds(s).beta,2));
Bx = spm_dctmtx(ds(s).P(1).dim(1),ds(s).order(1));
By = spm_dctmtx(ds(s).P(1).dim(2),ds(s).order(2));
Bz = spm_dctmtx(ds(s).P(1).dim(3),ds(s).order(3));
if isfield(ds(s),'sfP') && ~isempty(ds(s).sfP)
T = ds(s).sfP.mat\ds(1).P(1).mat;
txyz = xyz * T';
c = spm_bsplinc(ds(s).sfP,ds(s).hold);
ds(s).sfield = spm_bsplins(c,txyz(:,1),txyz(:,2),txyz(:,3),ds(s).hold);
ds(s).sfield = ds(s).sfield(:);
clear c txyz;
end
for i=1:size(ds(s).beta,2)
def_array(:,i) = spm_get_def(Bx,By,Bz,ds(s).beta(:,i));
end
if flags.udc > 1
ddef_array = zeros(prod(ds(s).P(1).dim(1:3)),size(ds(s).beta,2));
dBy = spm_dctmtx(ds(s).P(1).dim(2),ds(s).order(2),'diff');
for i=1:size(ds(s).beta,2)
ddef_array(:,i) = spm_get_def(Bx,dBy,Bz,ds(s).beta(:,i));
end
end
for i = 1:length(ds(s).P)
%
% Read undeformed image.
%
T = inv(ds(s).P(i).mat) * ds(1).P(1).mat;
txyz = xyz * T';
if flags.udc > 1
[def,jac] = spm_get_image_def(ds(s).P(i),ds(s),def_array,ddef_array);
else
def = spm_get_image_def(ds(s).P(i),ds(s),def_array);
end
txyz(:,2) = txyz(:,2) + def;
if flags.udc > 1
jP.dat = reshape(jac,ds(s).P(i).dim(1:3));
jtxyz = xyz * T';
c = spm_bsplinc(jP.dat,hold);
jac = spm_bsplins(c,jtxyz(:,1),jtxyz(:,2),jtxyz(:,3),hold);
end
c = spm_bsplinc(ds(s).P(i),hold);
ima = spm_bsplins(c,txyz(:,1),txyz(:,2),txyz(:,3),hold);
if flags.udc > 1
ima = ima .* jac;
end
%
% Write it if so required.
%
if flags.which
PO = ds(s).P(i);
PO.fname = spm_file(PO.fname, 'prefix',flags.prefix);
PO.mat = ds(1).P(1).mat;
PO.descrip = 'spm - undeformed';
ivol = ima;
if flags.mask
ivol(msk) = NaN;
end
ivol = reshape(ivol,PO.dim(1:3));
PO = spm_write_vol(PO,ivol);
if nargout > 0
oP{s}(i) = PO;
end
end
%
% Build up mean image if so required.
%
if flags.mean
Integral = Integral + nan2zero(ima);
end
spm_progress_bar('Set',tv);
tv = tv+1;
end
if isfield(ds(s),'sfield') && ~isempty(ds(s).sfield)
ds(s).sfield = [];
end
end
spm_progress_bar('Clear');
fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...done'); %-#
if flags.mean
fprintf('%-40s: %30s','Writing mean image...',''); %-#
sw = warning('off','MATLAB:divideByZero');
Integral = Integral./Count;
warning(sw);
PO = ds(1).P(1);
PO.fname = spm_file(ds(1).P(1).fname,'prefix',['mean' flags.prefix]);
PO.pinfo = [max(max(max(Integral)))/32767 0 0]';
PO.descrip = 'spm - mean undeformed image';
PO.dt = [spm_type('int16') spm_platform('bigend')];
PO.n = [1 1];
ivol = reshape(Integral,PO.dim);
spm_write_vol(PO,ivol);
fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...done'); %-#
end
if nargout > 0
varargout{1} = oP;
end
fprintf('%-40s: %30s\n','Completed',spm('time')) %-#
%==========================================================================
function v = nan2zero(v)
v(~isfinite(v)) = 0;
c4science · Help