If there's a script for updating the bvec file using the output of eddy, (the .eddy_parameters file), I would appreciate if someone would direct me to it.
If there isn't would someone mind checking out the sign conventions and the matrix ordering in the matlab code pasted below?
I checked the sign conventions by checking the output of eddy as applied to a dataset with artificially injected rotations.
I'm pretty sure about the sign conventions for the individual rotation matrices.
However, I couldn't get a clear conclusion about the composite rotation.
If it would help motivate you, I would be willing to write a bash version of this.
Thanks,
Ken
function update_bvec_eddy(filename_eddy_parameters, filename_bvecs, filename_updated_bvecs)
%% ------------------------------------------------------------------------
%% WE ASSUME THAT THE IMAGE ORIENTATION IS RPI.
%% FOR EXAMPLE, THE OUTPUT OF
%%
%% fslhd datafile.nii | grep qform | grep orient
%%
%% IS:
%%
%% qform_xorient Right-to-Left
%% qform_yorient Posterior-to-Anterior
%% qform_zorient Inferior-to-Superior
%%
%% WE ASSUME THAT THE BVECS ARE ORIENTED IN AGREEMENT WITH THE DATA
%% ------------------------------------------------------------------------
%% load text files
% Gradient vectors. size(bvecs) = [3, number_of_volumes]
bvecs = load(filename_bvecs);
% Output of eddy. size(eddy_parameters) = [number_of_volumes, 16]
eddy_parameters = load(filename_eddy_parameters);
%% Extract rotation angles
% Column 4-6 are rotation angles, as mentioned in:
% http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/EDDY/Faq
% "Will eddy rotate my bevcs [sic] for me?"
theta_x = eddy_parameters(:,4);
theta_y = eddy_parameters(:,5);
theta_z = eddy_parameters(:,6);
%% Construct rotation matrices
%% HERE IS ONE PLACE WHERE THINGS MIGHT GO WRONG BASED ON SIGN CONVENTION
% rotation around x-axis
c = cos(theta_x);
s = sin(theta_x);
Rx = [[1 0 0]; ...
[0 c -s]; ...
[0 s c]];
% rotation around y-axis
c = cos(theta_y);
s = sin(theta_y);
Ry = [[c 0 s]; ...
[0 1 0]; ...
[-s 0 c]];
% rotation around z-axis
c = cos(theta_z);
s = sin(theta_z);
Rz = [[c -s 0]; ...
[s c 0]; ...
[0 0 1]];
%% Construct overall, composite correction
%% HERE IS THE OTHER PLACE WHERE THINGS MIGHT GO WRONG
%% ONE ALTERNATIVE IS
%% R = Rx * Ry * Rz
R = inv(Rx * Ry * Rz);
%% Apply correction
bvecs_corrected = R * bvecs;
%% Save output
fid = fopen(filename_updated_bvecs, 'w');
for i = 1:3
fprintf(fid, '%.06f ', bvecs_corrected(i,:));
fprintf(fid, '\n');
end
fclose(fid);
|