function [ ] = UpdateGradient(TranMatrixfile, Schemefile, UpdatedFileName)
% UpdateGradient( TranMatrixfile, Schemefile, UpdatedFileName ) updates
% diffusion gradients from registration transformations.
%
% Input: 1) TranMatrixfile - Transformation Matrix File from mbalign -omat 
%           after registration, which contains all the transform matrics 
%           in ascii format.
%        2) Schemefile - Original scheme file for registered image.
% Output: UpdatedFileName - New scheme file with gradient directions updated 
%
% The result can be tested using 
% q_ori = [ 3-by-1 ];
% q_update = [ 3-by-1 ];
% degree3D=acos(sum(q_ori.*q_update)/(sqrt(sum(q_ori.^2))*sqrt(sum(q_update.^2))));
% degree3D/pi*180

% Author: Y.BAI
% Date: Sep.2007

% Read original scheme file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fid_1=fopen(Schemefile,'r','c');
tline = fgetl(fid_1); % 1st Line in scheme file : Diffusion time
Diffusion_Time=str2num(tline); 
tline = fgetl(fid_1); % 2nd Line in scheme file : Number of measurements
maxN=str2num(tline);
fclose(fid_1);

fid_1=fopen(Schemefile,'r','c');
no_line_scheme=(2+3*maxN);
qN=zeros(maxN*3,1);

for L=1:no_line_scheme
    tline = fgetl(fid_1);
    q_all(L)=str2num(tline);
end
fclose(fid_1);

% Read Transformation Matrix File
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
no_line_tran=4*maxN;
schemeMatrix=-100*ones(no_line_tran,4);

fid_0=fopen(TranMatrixfile,'r','c');
for L=1:no_line_tran
    tline = fgetl(fid_0); 
    schemeMatrix(L,:)=str2num(tline);
end
fclose(fid_0);

%  Extract Rotation Element
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F=-100*ones(4,4);
for N=1:maxN % N=14
    rowRange_start=1+4*(N-1);
    rowRange_end=4*N-1;
    F=schemeMatrix(rowRange_start:rowRange_end,1:3);
    RotMatrix(:,:,N)=(F*F')^(1/2)*F;
end

% Rotate q-Vector
%%%%%%%%%%%%%%%%%%
for N=1:maxN
    qx=q_all(3+3*(N-1));
    qy=q_all(4+3*(N-1));
    qz=q_all(5+3*(N-1));
    q_xyz=[qx;qy;qz];    
    R=RotMatrix(:,:,N);
    
%     outVec=R^(-1)*q_xyz;   % CAUTION: We are NOT 100 persents sure yet (16/07/2006)
    outVec=R*q_xyz;   % Comment: We think it should be R, NOT R^(-1). (28/08/2006)
    
    qN(1+3*(N-1))=outVec(1);
    qN(2+3*(N-1))=outVec(2);
    qN(3+3*(N-1))=outVec(3);
end

% Make new scheme file
%%%%%%%%%%%%%%%%%%%%%%%
fid_2=fopen(UpdatedFileName,'w','c');
for L=1:no_line_scheme
    if L==1
        fprintf(fid_2,'%s\n',num2str(Diffusion_Time));
    elseif L==2
        fprintf(fid_2,'%s\n',num2str(maxN));
    else
        fprintf(fid_2,'%s\n',num2str(qN(L-2)));
    end
end
fclose(fid_2);