123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960 |
- function def = spm_get_def(Bx,By,Bz,beta)
- % Calculating deformation field from coefficient vector.
- %
- % FORMAT [def] = get_def(Bx,By,Bz,beta);
- % Bx, By, Bz - Separable basis sets such that B=kron(Bz,kron(By,Bx));
- % beta - Coefficient vector for basis set.
- % or
- % FORMAT [def] = get_def(dim,order,beta);
- % dim - Dimensionality of image in x,y and z directions.
- % order - Order of DCT set in x,y and z directions.
- % beta - Coefficient vector for basis set.
- %_______________________________________________________________________
- %
- % Calculates the equivalent to kron(Bz,kron(By,Bx))*beta
- % in a faster and more efficient way. Note that this routine
- % can be used to calculate AtY efficiently as well since
- % A'*y = (diag(dfdy)*kron(Bz,kron(By,Bx)))'*y =
- % kron(Bz',kron(By',Bx'))*(diag(dfdy)*y) = get_def(Bx',By',Bz',dfdy.*y)
- %
- %_______________________________________________________________________
- % Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
- % Jesper Andersson
- % $Id: spm_get_def.m 1143 2008-02-07 19:33:33Z spm $
- if nargin == 4
- [nx,mx] = size(Bx);
- [ny,my] = size(By);
- [nz,mz] = size(Bz);
- if mx*my*mz ~= size(beta,1)
- warning('get_def: Size mismatch between beta and basis-set');
- return
- end
- elseif nargin == 3
- if numel(Bx) ~= 3 || numel(By) ~= 3
- warning('get_def: Wrong dimensionality on input');
- elseif prod(By) ~= size(Bz,1)
- warning('get_def: Size mismatch between beta and basis-set');
- return
- end
- nx = Bx(1); ny = Bx(2); nz = Bx(3);
- mx = By(1); my = By(2); mz = By(3);
- beta = Bz;
- Bx = spm_dctmtx(nx,mx);
- By = spm_dctmtx(ny,my);
- Bz = spm_dctmtx(nz,mz);
- end
- def = zeros(nx*ny*nz,1);
- for bf = 1:mz
- mbeta = reshape(beta((bf-1)*my*mx+1:bf*my*mx),mx,my);
- tmp = reshape(Bx*mbeta*By',nx*ny,1);
- for sl = 1:nz
- def((sl-1)*ny*nx+1:sl*ny*nx) = def((sl-1)*ny*nx+1:sl*ny*nx) + Bz(sl,bf)*tmp;
- end
- end
- return
|