Dear FSL people,
I am trying to generate a DWI from the parameters estimated in Bedpostx using the ball and sticks model equation:
Sj = S0[(1-sum(fn))*exp(-bj*d) + sum(fn*exp(-bj*d*(xjT * vn).^2))]
I know this is not something usual, but in theory I should be able to create a new 4D volume (for any xj and bj) calculating the Sj value in each of the voxels, since after Bedpostx I have an estimation from the unknowns in the above equation for each of the voxels.
I have written a matlab code which reads the nifti files and calculates the value of Sj in each voxel in this form (using the same bvecs and bvals):
% File paths
dataFileName = fullfile('~/100408/T1w/Diffusion/data.nii.gz');
b0FileName = fullfile(subj_dir,'mean_S0samples.nii.gz');
dFileName = fullfile(subj_dir,'mean_dsamples.nii.gz');
f1FileName = fullfile(subj_dir,'mean_f1samples.nii.gz');
f2FileName = fullfile(subj_dir,'mean_f2samples.nii.gz');
f3FileName = fullfile(subj_dir,'mean_f3samples.nii.gz');
dyad1FileName = fullfile(subj_dir,'dyads1.nii.gz');
dyad2FileName = fullfile(subj_dir,'dyads2.nii.gz');
%dyad2FileName = fullfile(subj_dir,'dyads2_thr0.05.nii.gz');
dyad3FileName = fullfile(subj_dir,'dyads3.nii.gz');
%dyad3FileName = fullfile(subj_dir,'dyads3_thr0.05.nii.gz');
% Read all Nifti files
data = niftiRead(dataFileName);
b0 = niftiRead(b0FileName);
d = niftiRead(dFileName);
f1 = niftiRead(f1FileName);
f2 = niftiRead(f2FileName);
f3 = niftiRead(f3FileName);
dyad1 = niftiRead(dyad1FileName);
dyad2 = niftiRead(dyad2FileName);
dyad3 = niftiRead(dyad3FileName);
% Load bval and bvec
bvecs = load(fullfile(subj_dir,'bvecs'));
bvals = load(fullfile(subj_dir,'bvals'));
% Normalize bvec
r = bvecs(:,w)/ sqrt(bvecs(:,w)' * bvecs(:,w));
b = bvals(w);
S0 = b0.data(x,y,z);
% Ball compartment
ball = (1 - (f1.data(x,y,z) + f2.data(x,y,z) + f3.data(x,y,z))) * exp(-b*d.data(x,y,z));
%Stick compartment
stick = f1.data(x,y,z) * exp(-b*d.data(x,y,z)*(r'*[dyad1.data(x,y,z,1);dyad1.data(x,y,z,2);dyad1.data(x,y,z,3)]).^2) + ...
f2.data(x,y,z) * exp(-b*d.data(x,y,z)*(r'*[dyad2.data(x,y,z,1);dyad2.data(x,y,z,2);dyad2.data(x,y,z,3)]).^2) + ...
f3.data(x,y,z) * exp(-b*d.data(x,y,z)*(r'*[dyad3.data(x,y,z,1);dyad3.data(x,y,z,2);dyad3.data(x,y,z,3)]).^2);
Sj = S0*(ball + stick);
However if I compare the value obtained from this formula with the measured one (the one from data.nii.gz) they are very different (MSE = 6.0865e+11).
Should I not be able to recover approximately the same value in each voxel?
I can not seem to find what I am doing wrong, so any ideas are much appreciated!
Best,
Maria
########################################################################
To unsubscribe from the FSL list, click the following link:
https://www.jiscmail.ac.uk/cgi-bin/WA-JISC.exe?SUBED1=FSL&A=1
This message was issued to members of www.jiscmail.ac.uk/FSL, a mailing list hosted by www.jiscmail.ac.uk, terms & conditions are available at https://www.jiscmail.ac.uk/policyandsecurity/
|