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
| 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 uen DCT(transformée en cosinus discrète)
basX = spm_dctmtx(VG(1).dim(1),k(1))*stabilise;
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;
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,%???
% 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)) +...
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; |
