Project: High Angular Resolution Diffusion Imaging Tools

Example to illustrate the Diffusion Orientation Transform (DOT) method using real spherical harmonic basis. The ODF is a contour of the propagator, evaluated at a given displacement length $$ l $$ :

$$ ODF(u) = {Propagator(u,r=l)}. $$

Contents

- Parameters for reconstruction: All this part is equal for each voxel

clear all; close all; clc;

% - Loading the measurement scheme

grad = load('/sampling_and_reconstruction_schemes/On_the_sphere/70_shell.txt');
% grad is a Nx3 matrix where grad(i,:) = [xi yi zi] and sqrt(xi^2 + yi^2 + zi^2) == 1 (unit vectors).
display(['The measurement scheme contains ' num2str(size(grad,1)) ' directions']);

% - Loading the reconstruction scheme (mesh where the ODF will be reconstructed)
%   ---------------------------------
% unit vectors on the sphere
disp(' ')
V = load('/sampling_and_reconstruction_schemes/On_the_sphere/724_shell.txt');
% facets
display(['The reconstruction scheme contains ' num2str(size(V,1)) ' directions']);
F = load('/sampling_and_reconstruction_schemes/On_the_sphere/724_sphere_facets.txt');
disp(' ');

%  --- Regularization parameter
Lambda = 0.006;
display(['The regularization parameter Lambda is ' num2str(Lambda) ',']);
display ('like in the Descoteaux paper "Regularized, Fast, and Robust Analytical Q-Ball Imaging".');

% --- real spherical harmonic reconstruction: parameters definition --- %
[Lmax Nmin] = obtain_Lmax(grad);
if Lmax >= 8
    Lmax = 8;
end
display(['The maximum order of the spherical harmonics decomposition is Lmax = ' num2str(Lmax) '.']);

[basisG, thetaG, phiG] = construct_SH_basis (Lmax, grad, 2, 'real');
[basisV, thetaV, phiV] = construct_SH_basis (Lmax, V, 2, 'real');
K = []; Laplac2 = []; Vector = [];
for L=0:2:Lmax
    for m=-L:L
        factor1 = (-1)^(L/2)*gamma((L+3)/2)/gamma(L + 3/2);
        K = [K; factor1];
        Laplac2 = [Laplac2; (L^2)*(L + 1)^2];
        Vector = [Vector L];
    end
end
Laplac = diag(Laplac2);

long = 0:2:Lmax; NL = length(long);
Ncoeff = (Lmax+1)*(Lmax+2)/2;
Matrix = zeros(Ncoeff,NL);

for i=1:length(long);
    ind = find(Vector==long(i));
    Matrix(ind,i) = 1;
end

% Creating the kernel for the reconstruction
A = recon_matrix(basisG,Laplac,Lambda);
% where A = inv(Ymatrix'*Ymatrix + Lambda*Laplacian)*Ymatrix';
% and Ymatrix = basisG.
% The matrix, A, is the same for all the voxels

The measurement scheme contains 70 directions
 
The reconstruction scheme contains 724 directions
 
The regularization parameter Lambda is 0.006,
like in the Descoteaux paper "Regularized, Fast, and Robust Analytical Q-Ball Imaging".
The maximum order of the spherical harmonics decomposition is Lmax = 8.

- Simulated diffusion MRI signal using a multi-tensor model

% --- Experimental parameters --- %
% b-value
b = 3*10^3;% s/mm^2
% diffusivities
lambda1 = 1.5*10^-3;% mm^2/s
lambda2 = 0.3*10^-3;% mm^2/s
% S0 is the image at b-value = 0
S0 = 100;
% Signal to noise ratio
SNR = 20;
display(['The signal to noise ratio in this experiment is SNR = ' num2str(SNR) '.']);
% Add a measurement with b-value = 0 into the gradient matrix (standard format in which the signal is measured in the scanner)
grad = [0 0 0; grad];
% Note that more measurements for b-value = 0 can be added...

% --- Intravoxel properties --- %
% Using a mixture of two fibers (more fibers can be used)
% volume fraction for each fiber (the sum of this elemenst is 1 by definition)
f1 = 0.5; f2 = 0.5;
% angular orientation corresponding to each fiber (main eigenvector for each diffusion tensor)
anglesF = [0  0;
           70 0];
% [ang1 ang2]: -> ang1 is the rotation in degrees from x to y axis (rotation using the z-axis);
%                 ang2 is the rotation in degrees from the x-y plane to the z axis (rotation using the y axis)
% anglesF(i,:) defines the orientation of the fiber with volume fraction fi.

%  --- Computation of the diffusion signal using the above experimental parameters and intravoxel properties %
% *** Signal with realistic rician noise
add_rician_noise = 1;
S = create_signal_multi_tensor(anglesF, [f1  f2],...
    [lambda1, lambda2, lambda2], b, grad, S0, SNR, add_rician_noise);

% *** Ideal signal (no noise)
add_rician_noise = 0;
Snonoise = create_signal_multi_tensor(anglesF, [f1  f2],...
    [lambda1, lambda2, lambda2], b, grad, S0, SNR, add_rician_noise);

The signal to noise ratio in this experiment is SNR = 20.

- Isolate the S0 component from the signal

% --- find zeros
grad_sum = sum(abs(grad),2);
ind_0 = find( grad_sum == 0);
S0_est = max(S(ind_0));
S(ind_0) = []; Snonoise(ind_0) = [];
grad(ind_0,:) = [];

- DOT-parameters definition

long = 0:2:Lmax;
Ro = 15*10^-3; % mm
display (['The diffusion propagator will be evaluated at a displacement length of ' num2str(Ro) ' mm']);
t = 20*10^-3; % s
display (['The experimental diffusion time is assumed to be ' num2str(t) ' seconds']);

The diffusion propagator will be evaluated at a displacement length of 0.015 mm
The experimental diffusion time is assumed to be 0.02 seconds

- DOT reconstruction from the signal without noise

% ---
ADC = abs(-log(Snonoise/S0)/b)';
Dto = ADC.*t;
% ---
A0 = recon_matrix(basisG,Laplac,0);
% Note that in this case we created an additional matrix A0,
% this matrix is different to A because the regularization parameter is set
% to zero (Lambda = 0).
ODFDot_nonoise = ReconstructionDOT(A0, Lmax, basisV, K, Ro, Dto, Matrix);
ODFDot_nonoise(ODFDot_nonoise<0) = 0; % discarding negative values, which also appears in the noise-free signal!
ODFDot_nonoise = ODFDot_nonoise/sum(ODFDot_nonoise);

- Regularized DOT reconstruction from the noised signal

display('This is the real/effective time required for computing the ODF at each voxel:')
tic
ADC = abs(-log(S/S0_est)/b)';
% --- ADC denoising
ADCd = (basisG*(A*ADC(:))); % denoised ADC using the spherical harmonic series
% ---
Dt = ADCd.*t;
ODFDot = ReconstructionDOT(A0, Lmax, basisV, K, Ro, Dt, Matrix);
ODFDot(ODFDot<0) = 0; % discarding negative values
ODFDot = ODFDot/sum(ODFDot);
toc

This is the real/effective time required for computing the ODF at each voxel:
Elapsed time is 0.015737 seconds.

- Plots

% --- opengl('software');
feature('UseMesaSoftwareOpenGL',1);

set(gcf, 'color', 'white');
screen_size = get(0, 'ScreenSize');
f1 = figure(1); hold on;
set(f1, 'Position', [0 0 screen_size(3) screen_size(4) ] );
cameramenu;
% -----
Origen = [0 0 0.25]; % Voxel coordinates [x y z] where the ODF will be displayed
angle = 0:.01:2*pi;
R = ones(1,length(angle));
R1 = R.*f2/f1;
subplot(1,2,1)
polarm(angle,R,'k'); hold on; polarm(angle, R1,'k');
% the ODF is scaled only for the display
plot_ODF(ODFDot_nonoise./max(ODFDot_nonoise), V, F, Origen);
title({['\fontsize{14} Propagator evaluated at Ro = ' num2str(mean(Ro)) 'mm.'];'(Signal without noise)'},'FontWeight','bold');
axis equal; axis off;

subplot(1,2,2)
polarm(angle,R,'k'); hold on; polarm(angle, R1,'k');
% the ODF is scaled only for the display
plot_ODF(ODFDot./max(ODFDot),V, F, Origen);
title({['\fontsize{14} Propagator evaluated at Ro = ' num2str(mean(Ro)) 'mm.'];...
       ['SNR = ' num2str(SNR) '; Regularization parameter {\lambda}  = ' num2str(Lambda)]},'FontWeight','bold');
axis equal; axis off;


Published with MATLAB®