The first thing we need to do is to setup the G matrix given the diffusion encoding directions (and the b-value, if you want quantitative measures such as ADC).

[cc lang="matlab"]
%% First we need to find out how many of the gradient
%% directions are all zero…
zeroi = find( sum(gradients’)’ == 0 );
nonzeroi = setdiff( 1:nrows(gradients), zeroi );

gradients = gradients( nonzeroi, : );

%% Setup G
for ii=1:nrows(gradients)

G(ii,:) = b_value * [ gradients(ii,:).*gradients(ii,:) ...
2*gradients(ii,1).*gradients(ii,2) ...
2*gradients(ii,1).*gradients(ii,3) ...
2*gradients(ii,2).*gradients(ii,2) ];
end

Ginv = pinv(G);
[/cc]

Now, given the G matrix, we want to calculate the diffusion tensor values given the gradient diffusion encoding, G, and the acquired signal y. tv is the tensor values (in vector form) and TV is the tensor matrix. Then we use singular value decomposition (svd) to compute the eigenvalues and eigenvectors of the tensor (per voxel).

[cc lang="matlab"]
for ii=1:nrows(d)
for jj=1:ncols(d)
for sl=1:size(d,3)
 
data = squeeze( d(ii,jj,sl,:) );
data = data(nonzeroi) ./ mean(data(zeroi));
data = – log( data );
 
% tensor values
tv = Ginv * data;
 
TV=zeros(3,3);
TV(1,1) = tv(1); TV(2,2) = tv(2); TV(3,3) = tv(3);
TV(1,2) = tv(4); TV(2,1) = tv(4);
TV(1,3) = tv(5); TV(3,1) = tv(5);
TV(2,3) = tv(6); TV(3,2) = tv(6);
 
[U,S,V] = svd(TV);
 
dt(ii,jj,sl).evalues = diag(S);
dt(ii,jj,sl).evector = U(:,1);
end
end
end
[/cc]
Now, given the eigenvalues values, evalues in Matlab, and corresponding eigenvectors, evector in Matlab, we can compute lots of different possible measures.