Bonjour,

j'essaye de comprendre ce code, mais j'y arrive pas du tout;
j'ai compris quelques lignes du code(j'ai mis des commentaires en français), mais pour le reste, j'ai beau cherché je ne trouve pas ce qu'on y fait exactement.
Merci beaucoup pour votre aide

Code : Sélectionner tout - Visualiser dans une fenêtre à part
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;