1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92
| function Tr = snbasis(VG,VF,VWG,VWF,Affine,fwhm,cutoff,nits,reg)
% 3D Basis Function Normalization
% FORMAT Tr = snbasis(VG,VF,VWG,VWF,Affine,fwhm,cutoff,nits,reg)
% VG - Template volumes (see spm_vol).
% VF - Volume to normalize.
% VWG - weighting Volume - for template.
% VWF - weighting Volume - for object.
% Affine - A 4x4 transformation (in voxel space).
% fwhm - smoothness of images.
% cutoff - frequency cutoff of basis functions.
% nits - number of iterations.
% reg - regularisation.
% Tr - Discrete cosine transform of the warps in X, Y & Z.
%
% snbasis performs a spatial normalization based upon a 3D
% discrete cosine transform.
%
%______________________________________________________________________
fwhm = [fwhm 30];
% Number of basis functions for x, y & z
%-----------------------------------------------------------------------
tmp = sqrt(sum(VG(1).mat(1:3,1:3).^2));
k = max(round((VG(1).dim(1:3).*tmp)/cutoff),[1 1 1]);
% Scaling is to improve stability.
%-----------------------------------------------------------------------
stabilise = 8; % on fait une DCT(transformée en cosinus discrète)
basX = spm_dctmtx(VG(1).dim(1),k(1))*stabilise; %matrice carée de la taille
basY = spm_dctmtx(VG(1).dim(2),k(2))*stabilise;
basZ = spm_dctmtx(VG(1).dim(3),k(3))*stabilise;
% j'obtiens les dérivées de la DCT
dbasX = spm_dctmtx(VG(1).dim(1),k(1),'diff')*stabilise;
dbasY = spm_dctmtx(VG(1).dim(2),k(2),'diff')*stabilise;
dbasZ = spm_dctmtx(VG(1).dim(3),k(3),'diff')*stabilise;
vx1 = sqrt(sum(VG(1).mat(1:3,1:3).^2));
vx2 = vx1;
% %définition des trois directions orthogonales
kx = (pi*((1:k(1))'-1)/VG(1).dim(1)/vx1(1)).^2; ox=ones(k(1),1);%??
ky = (pi*((1:k(2))'-1)/VG(1).dim(2)/vx1(2)).^2; oy=ones(k(2),1);%??
kz = (pi*((1:k(3))'-1)/VG(1).dim(3)/vx1(3)).^2; oz=ones(k(3),1);%??
if 1,%du quick end dirty, la condition va toujours être vérifié!!!
% BENDING ENERGY REGULARIZATION
% Estimate a suitable sparse diagonal inverse covariance matrix for
% the parameters (IC0).
%-----------------------------------------------------------------------
IC0 = (1*kron(kz.^2,kron(ky.^0,kx.^0)) +... %kron:coeiff de kroneker
1*kron(kz.^0,kron(ky.^2,kx.^0)) +...
1*kron(kz.^0,kron(ky.^0,kx.^2)) +...
2*kron(kz.^1,kron(ky.^1,kx.^0)) +...
2*kron(kz.^1,kron(ky.^0,kx.^1)) +...
2*kron(kz.^0,kron(ky.^1,kx.^1)) );
IC0 = reg*IC0*stabilise^6;
IC0 = [IC0*vx2(1)^4 ; IC0*vx2(2)^4 ; IC0*vx2(3)^4 ; zeros(prod(size(VG))*4,1)];
IC0 = sparse(1:length(IC0),1:length(IC0),IC0,length(IC0),length(IC0));
else
% MEMBRANE ENERGY (LAPLACIAN) REGULARIZATION
%-----------------------------------------------------------------------
IC0 = kron(kron(oz,oy),kx) + kron(kron(oz,ky),ox) + kron(kron(kz,oy),ox);
IC0 = reg*IC0*stabilise^6;
IC0 = [IC0*vx2(1)^2 ; IC0*vx2(2)^2 ; IC0*vx2(3)^2 ; zeros(prod(size(VG))*4,1)];
IC0 = sparse(1:length(IC0),1:length(IC0),IC0,length(IC0),length(IC0));
end;
% Generate starting estimates.
%-----------------------------------------------------------------------
s1 = 3*prod(k);
s2 = s1 + numel(VG)*4;
T = zeros(s2,1);
T(s1+(1:4:numel(VG)*4)) = 1;
pVar = Inf;
for iter=1:nits,
fprintf(' iteration %2d: ', iter);
[Alpha,Beta,Var,fw] = spm_brainwarp(VG,VF,Affine,basX,basY,basZ,dbasX,dbasY,dbasZ,T,fwhm,VWG, VWF);
if Var>pVar, scal = pVar/Var ; Var = pVar; else scal = 1; end;
pVar = Var;
T = (Alpha + IC0*scal)\(Alpha*T + Beta);
fwhm(2) = min([fw fwhm(2)]);
fprintf(' FWHM = %6.4g Var = %g\n', fw,Var);
end;
% Values of the 3D-DCT - for some bizarre reason, this needs to be done
% as two seperate statements in Matlab 6.5...
%-----------------------------------------------------------------------
Tr = reshape(T(1:s1),[k 3]);
drawnow;
Tr = Tr*stabilise.^3;
return; |
Partager