Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

bval/bvec files for Bruker enhanced DICOMs #265

Closed
isolovey opened this issue Feb 1, 2019 · 33 comments
Closed

bval/bvec files for Bruker enhanced DICOMs #265

isolovey opened this issue Feb 1, 2019 · 33 comments

Comments

@isolovey
Copy link
Contributor

isolovey commented Feb 1, 2019

I have another issue with Enhanced MR IOD DICOMs from Bruker. The bval/bvec files are not extracted from DTI series. The issue is likely that Bruker encodes the diffusion direction as a B-matrix. The bvecs need to be computed from that.

I put up two DTI series on Dropbox. There is also a Matlab script another lab member wrote to output bval/bvec files given a DICOM file as input.

@neurolabusc
Copy link
Collaborator

neurolabusc commented Feb 1, 2019

The Matlab script looks handy. The code suggests that Bruker Gradient matrix is reported with respect to the Freq (X), Phase (Y), Slice (Z) just like GE. In contrast, Siemens and Philips report gradient directions with respect to the scanner bore. Can you please confirm this, otherwise your script will provide incorrect results when slice angulations are applied. Also, since Bruker systems are often used in animals, I think we would want to make sure we have a generalized solution that works for axial, coronal and sagittal acquisitions. For humans, DWI scans are almost always axial as the S->I distance is shorter than full brain coronal A->P slices. This is not the case in rodents. The wiki links to a "dedicated document", e.g. Validate DTI Vectors that describes acquiring validation DWI datasets.

@isolovey
Copy link
Contributor Author

isolovey commented Feb 1, 2019

I'll ask around for a confirmation regarding the reference frame. I think your second link to the document is broken.

@neurolabusc
Copy link
Collaborator

Here is a link to a Dropbox zip archive that uses the Matlab code you provided as well as a C port by Connelly Barnes of the public domain Java Matrix library JAMA.

Note that the sign of bvec solutions differs for some of the vectors. (i=input matrix, e=eigenmatrix, v=3rd column of e). For example for the last vector in your 20 direction file:

C
m=[78.9072 269.576 -25.1249; 269.576 269.576 920.972; -25.1249 920.972 -85.8359];
e=[0.198185 -0.96076 0.194068; -0.639172 0.0234265 0.768707; 0.743089 0.276389 0.609448];
v = [0.194068 0.768707 0.609448];
Matlab
i=[78.9072 269.576 -25.1249; 269.576 269.576 920.972; -25.1249 920.972 -85.8359];
e=[0.198185 -0.96076 -0.194068; -0.639172 0.0234265 -0.768707; 0.743089 0.276389 -0.609448];
v = [-0.194068 -0.768707 -0.609448];

For computing preferential diffusion direction gradients with the precisely opposite polarity will sensitize the same orthogonal plane. However, I worry that Eddy uses this. I am not sure what the correct result is for the symmetric matrix provided by Bruker.

If you can resolve these tricky questions, the C code should be easy to plug in to dcm2niix (I added the function to the existing nifti1_io_core.cpp). Until then, I am afraid this is as far as I can go. Perhaps the engineers at Bruker can help.

@neurolabusc
Copy link
Collaborator

I'll ask around for a confirmation regarding the reference frame. I think your second link to the document is broken.

Fixed link

@neurolabusc
Copy link
Collaborator

According to the documentation, the rows and columns of the matrix are the X (right to left), Y (anterior to posterior) and Z (foot to head) patient-relative orthogonal axes.

@neurolabusc
Copy link
Collaborator

Just to clarify the situation with my C code and your Matlab code. You are suggesting that we compute [V,D] = eig(A) and the final column of the matrix V will be the gradient vector. The function [V,D] = eig(A) finds a solution where AV = VD. However, for this symmetric matrix we get a good solution if the final column is all positive or all negative. Therefore, the two solutions provide vectors that are parallel but point in precisely the opposite direction. How do we determine which is the correct direction that captures the spatial bias induced by the eddy currents?


C solution:


A=[78.9072 269.576 -25.1249; 269.576 269.576 920.972; -25.1249 920.972 -85.8359];
V=[0.198185 -0.96076 0.194068; -0.639172 0.0234265 0.768707; 0.743089 0.276389 0.609448];
D = [-884.715 0 0; 0 79.5619 0; 0 0 1067.8]


Matlab: Test C solution


format longg
A=[78.9072 269.576 -25.1249; 269.576 269.576 920.972; -25.1249 920.972 -85.8359];
V=[0.198185 -0.96076 0.194068; -0.639172 0.0234265 0.768707; 0.743089 0.276389 0.609448];
D = [-884.715 0 0; 0 79.5619 0; 0 0 1067.8]
AV
V
D

V =

              0.198185                  -0.96076                  0.194068
             -0.639172                 0.0234265                  0.768707
              0.743089                  0.276389                  0.609448

D =

              -884.715                         0                         0
                     0                   79.5619                         0
                     0                         0                    1067.8

ans =

       -175.3372444561            -76.4399052941            207.2260006664
         565.484650996          1.86391451200001             820.825576856
       -657.4226065856             21.9900509169            650.7691665276

ans =

        -175.337242275             -76.439891044               207.2258104
          565.48505598             1.86385685035               820.8253346
        -657.421984635             21.9900339791               650.7685744

Matlab: Generate and Test Matlab solution


A=[78.9072 269.576 -25.1249; 269.576 269.576 920.972; -25.1249 920.972 -85.8359];
[V,D] = eig(A) %where AV = VD
V
D
AV
V
D

V =

     0.198184949305482        -0.960760276207933        -0.194068074472953
    -0.639171914880749        0.0234264839528279        -0.768707007303358
     0.743089489292144         0.276389022051282        -0.609448208950585

D =

     -884.715203742877                         0                         0
                     0          79.5619197231669                         0
                     0                         0          1067.80058401971

ans =

     -175.337237803571         -76.4399319688634           -207.2260032618
      565.485110900447          1.86385603565098         -820.825791338569
      -657.42256891829          21.9900411848088         -650.769153447201 

ans =

     -175.337237803571         -76.4399319688632           -207.2260032618
      565.485110900447          1.86385603565095         -820.825791338569
      -657.42256891829          21.9900411848087         -650.769153447201

@neurolabusc
Copy link
Collaborator

  1. Please extensively test latest commit on the development branch.
  2. Validation datasets would help greatly.
  3. The Bruker files are underspecified to generate FSL format bvec files. One can not determine the polarity of the bvec.
  4. The Matlab code illustrates the issue using Siemens XA10 data that provides the bmatrix, bvec and bval. The core issue an the eigenvector polarity is arbitrary (it is the same eigenvector, just a different eigenvalue). This explains why the estimated bvec often differs from the actual bvec by 180 degrees. Also note that code also estimates the bmatrix given the bvec and bval. Note the same bmatrix is generated regardless of whether the bvec is flipped.
  5. For most Diffusion computations (FA, tractography, etc) the polarity of the bvec does not matter. However, FSL's Eddy exploits the fact that while the diffusion signal is the same the spatial distortion is different.
  6. The latest commit includes a kludge: it enforces that the bvec X-value is positive. This makes Eddy see the acquisition as half shell (and therefore it should not attempt to correct differences that reflect noise).
  7. When converting these images dcm2niix will generate a warning Eddy users: B-matrix does not encode b-vector polarity (issue 265). The user should be aware that Eddy undistortion will not be as effective as if we could accurately determine the polarity of the b-vectors.
  8. I strongly suggest you work with the vendor to include the gradient vector (0018,9089) in future releases.
function bmat2bvec
%converts diffusion bmatrix to bvector, converts bvector+bvalue to bmatrix
% Rationale: 
%   Bruker only provides BMatrix+BVal, FSL expects BVec
%   https://github.com/rordenlab/dcm2niix/issues/265
% Issue:
%   The polarity of the BVec can not be extracted from the BMatrix
%    A vector and a flipped vector are the same eigen vector with different eigen value
%    "The choice of sign is arbitrary" 
%        http://nipy.org/dipy/theory/bmatrix.html
%    "Because the solution is a square root, the sign of the returned vector is arbitrary.  We set the vector to have a positive x component by convention."
%       https://raw.githubusercontent.com/matthew-brett/nibabel/master/nibabel/nicom/dwiparams.py
%    For FA/Diffusion the polarity of the gradient does not matter
%    However, it DOES matter for FSL's eddy
%      https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/eddy#eddy_--_a_tool_for_correcting_eddy_currents_and_movements_in_diffusion_data
%    Since the BMatrix is underspecified for estimating BVec polarity, 
%      we follow the convention of making the X vector positive
%      Eddy will see this as a half-shell acquisition
%    This script demonstrates the problem in two ways:
%       1. The estimated bvec is often 180-degrees different from actual bvec
%       2. The estimated bmat if identical for a flipped or non-flipped bvec
%
%The sample data comes from Siemens XA10 which saves both BVec (0018,9089) and BMatrix

% (0018,9089) FD 0.70429527759552002\0\-0.70990711450576782 #  24, 3 DiffusionGradientOrientation
% (0018,9087) FD 250                                      #   8, 1 DiffusionBValue
% (0018,9602) FD 125                                      #   8, 1 DiffusionBValueXX
% (0018,9603) FD 0                                        #   8, 1 DiffusionBValueXY
% (0018,9604) FD -126                                     #   8, 1 DiffusionBValueXZ
% (0018,9605) FD 0                                        #   8, 1 DiffusionBValueYY
% (0018,9606) FD 0                                        #   8, 1 DiffusionBValueYZ
% (0018,9607) FD 127                                      #   8, 1 DiffusionBValueZZ
bval = 250;
bvec = [0.70429527759552002, 0, -0.70990711450576782];
bmats = [125,0,-126,0,0,127];
testBMat(bmats, bvec, bval);

% (0018,9089) FD -0.70429527759552002\0\-0.70990711450576782 #  24, 3 DiffusionGradientOrientation
% (0018,9087) FD 250                                      #   8, 1 DiffusionBValue
% (0018,9602) FD 125                                      #   8, 1 DiffusionBValueXX
% (0018,9603) FD 0                                        #   8, 1 DiffusionBValueXY
% (0018,9604) FD 126                                      #   8, 1 DiffusionBValueXZ
% (0018,9605) FD 0                                        #   8, 1 DiffusionBValueYY
% (0018,9606) FD 0                                        #   8, 1 DiffusionBValueYZ
% (0018,9607) FD 127                                      #   8, 1 DiffusionBValueZZ
bval = 250;
bvec = [-0.70429527759552002, 0, -0.70990711450576782];
bmats = [125,0,126,0,0,127];
testBMat(bmats, bvec, bval);

% (0018,9089) FD 0\-0.70429527759552002\-0.70990711450576782 #  24, 3 DiffusionGradientOrientation
% (0018,9087) FD 250                                      #   8, 1 DiffusionBValue
% (0018,9602) FD 0                                        #   8, 1 DiffusionBValueXX
% (0018,9603) FD 0                                        #   8, 1 DiffusionBValueXY
% (0018,9604) FD 0                                        #   8, 1 DiffusionBValueXZ
% (0018,9605) FD 125                                      #   8, 1 DiffusionBValueYY
% (0018,9606) FD 126                                      #   8, 1 DiffusionBValueYZ
% (0018,9607) FD 127                                      #   8, 1 DiffusionBValueZZ
bval = 250;
bvec = [0, -0.70429527759552002, -0.70990711450576782];
bmats = [0,0,0,125,126,127];
testBMat(bmats, bvec, bval);

% (0018,9089) FD 0\-0.70995223522186279\0.70424985885620117 #  24, 3 DiffusionGradientOrientation
% (0018,9087) FD 250                                      #   8, 1 DiffusionBValue
% (0018,9602) FD 0                                        #   8, 1 DiffusionBValueXX
% (0018,9603) FD 0                                        #   8, 1 DiffusionBValueXY
% (0018,9604) FD 0                                        #   8, 1 DiffusionBValueXZ
% (0018,9605) FD 125                                      #   8, 1 DiffusionBValueYY
% (0018,9606) FD -124                                     #   8, 1 DiffusionBValueYZ
% (0018,9607) FD 123                                      #   8, 1 DiffusionBValueZZ
bval = 250;
bvec = [0, -0.70995223522186279, 0.70424985885620117];
bmats = [0,0,0,125,-124,123];
testBMat(bmats, bvec, bval);

% (0018,9089) FD 0.7070954442024231\-0.70709550380706787\-0.0056565827690064907 #  24, 3 DiffusionGradientOrientation
% (0018,9087) FD 250                                      #   8, 1 DiffusionBValue
% (0018,9602) FD 125                                      #   8, 1 DiffusionBValueXX
% (0018,9603) FD -125                                     #   8, 1 DiffusionBValueXY
% (0018,9604) FD -1                                       #   8, 1 DiffusionBValueXZ
% (0018,9605) FD 125                                      #   8, 1 DiffusionBValueYY
% (0018,9606) FD 1                                        #   8, 1 DiffusionBValueYZ
% (0018,9607) FD 0                                        #   8, 1 DiffusionBValueZZ
bval = 250;
bvec = [0.7070954442024231, -0.70709550380706787, -0.0056565827690064907];
bmats = [125,-125,-1,125,1,0];
testBMat(bmats, bvec, bval);

% (0018,9089) FD -0.7070954442024231\-0.70709550380706787\-0.0056565827690064907 #  24, 3 DiffusionGradientOrientation
% (0018,9087) FD 250                                      #   8, 1 DiffusionBValue
% (0018,9602) FD 125                                      #   8, 1 DiffusionBValueXX
% (0018,9603) FD 125                                      #   8, 1 DiffusionBValueXY
% (0018,9604) FD 1                                        #   8, 1 DiffusionBValueXZ
% (0018,9605) FD 125                                      #   8, 1 DiffusionBValueYY
% (0018,9606) FD 1                                        #   8, 1 DiffusionBValueYZ
% (0018,9607) FD 0                                        #   8, 1 DiffusionBValueZZ
bval = 250;
bvec = [-0.7070954442024231, -0.70709550380706787, -0.0056565827690064907];
bmats = [125,125,1,125,1,0];
testBMat(bmats, bvec, bval);

% (0018,9089) FD 0.70587193965911865\0\-0.7083393931388855 #  24, 3 DiffusionGradientOrientation
% (0018,9087) FD 2000                                     #   8, 1 DiffusionBValue
% (0018,9602) FD 1000                                     #   8, 1 DiffusionBValueXX
% (0018,9603) FD 0                                        #   8, 1 DiffusionBValueXY
% (0018,9604) FD -1003                                    #   8, 1 DiffusionBValueXZ
% (0018,9605) FD 0                                        #   8, 1 DiffusionBValueYY
% (0018,9606) FD 0                                        #   8, 1 DiffusionBValueYZ
% (0018,9607) FD 1007                                     #   8, 1 DiffusionBValueZZ
bval = 2000;
bvec = [0.70587193965911865, 0, -0.7083393931388855];
bmats = [1000,0,-1003,0,0,1007];
testBMat(bmats, bvec, bval);

% (0018,9089) FD -0.70587193965911865\0\-0.7083393931388855 #  24, 3 DiffusionGradientOrientation
% (0018,9087) FD 2000                                     #   8, 1 DiffusionBValue
% (0018,9602) FD 1000                                     #   8, 1 DiffusionBValueXX
% (0018,9603) FD 0                                        #   8, 1 DiffusionBValueXY
% (0018,9604) FD 1003                                     #   8, 1 DiffusionBValueXZ
% (0018,9605) FD 0                                        #   8, 1 DiffusionBValueYY
% (0018,9606) FD 0                                        #   8, 1 DiffusionBValueYZ
% (0018,9607) FD 1007                                     #   8, 1 DiffusionBValueZZ
bval = 2000;
bvec = [-0.70587193965911865, 0, -0.7083393931388855];
bmats = [1000,0,1003,0,0,1007];
testBMat(bmats, bvec, bval);

% (0018,9089) FD 0\-0.70587193965911865\-0.7083393931388855 #  24, 3 DiffusionGradientOrientation
% (0018,9087) FD 2000                                     #   8, 1 DiffusionBValue
% (0018,9602) FD 0                                        #   8, 1 DiffusionBValueXX
% (0018,9603) FD 0                                        #   8, 1 DiffusionBValueXY
% (0018,9604) FD 0                                        #   8, 1 DiffusionBValueXZ
% (0018,9605) FD 1000                                     #   8, 1 DiffusionBValueYY
% (0018,9606) FD 1003                                     #   8, 1 DiffusionBValueYZ
% (0018,9607) FD 1007                                     #   8, 1 DiffusionBValueZZ
bval = 2000;
bvec = [0, -0.70587193965911865, -0.7083393931388855];
bmats = [0,0,0,1000,1003,1007];
testBMat(bmats, bvec, bval);

% (0018,9089) FD 0\-0.70834684371948242\0.70586460828781128 #  24, 3 DiffusionGradientOrientation
% (0018,9087) FD 2000                                     #   8, 1 DiffusionBValue
% (0018,9602) FD 0                                        #   8, 1 DiffusionBValueXX
% (0018,9603) FD 0                                        #   8, 1 DiffusionBValueXY
% (0018,9604) FD 0                                        #   8, 1 DiffusionBValueXZ
% (0018,9605) FD 1000                                     #   8, 1 DiffusionBValueYY
% (0018,9606) FD -997                                     #   8, 1 DiffusionBValueYZ
% (0018,9607) FD 993                                      #   8, 1 DiffusionBValueZZ
bval = 2000;
bvec = [0, -0.70834684371948242, 0.70586460828781128];
bmats = [0,0,0,1000,-997,993];
testBMat(bmats, bvec, bval);

% (0018,9089) FD 0.70710510015487671\-0.70710533857345581\-0.0021213062573224306 #  24, 3 DiffusionGradientOrientation
% (0018,9087) FD 2000                                     #   8, 1 DiffusionBValue
% (0018,9602) FD 1000                                     #   8, 1 DiffusionBValueXX
% (0018,9603) FD -1000                                    #   8, 1 DiffusionBValueXY
% (0018,9604) FD -3                                       #   8, 1 DiffusionBValueXZ
% (0018,9605) FD 1000                                     #   8, 1 DiffusionBValueYY
% (0018,9606) FD 3                                        #   8, 1 DiffusionBValueYZ
% (0018,9607) FD 0                                        #   8, 1 DiffusionBValueZZ
bval = 2000;
bvec = [0.70710510015487671, -0.70710533857345581, -0.0021213062573224306];
bmats = [1000,-1000,-3,1000,3,0];
testBMat(bmats, bvec, bval);

% (0018,9089) FD -0.70710510015487671\-0.70710533857345581\-0.0021213062573224306 #  24, 3 DiffusionGradientOrientation
% (0018,9087) FD 2000                                     #   8, 1 DiffusionBValue
% (0018,9602) FD 1000                                     #   8, 1 DiffusionBValueXX
% (0018,9603) FD 1000                                     #   8, 1 DiffusionBValueXY
% (0018,9604) FD 3                                        #   8, 1 DiffusionBValueXZ
% (0018,9605) FD 1000                                     #   8, 1 DiffusionBValueYY
% (0018,9606) FD 3                                        #   8, 1 DiffusionBValueYZ
% (0018,9607) FD 0                                        #   8, 1 DiffusionBValueZZ
bval = 2000;
bvec = [-0.70710510015487671, -0.70710533857345581, -0.0021213062573224306];
bmats = [1000,1000,3,1000,3,0];
testBMat(bmats, bvec, bval);
%end bmat2bvec()

function testBMat(bmats, bvec, bval)
fprintf('STEP1: estimate bvec from known bmat;\n')
bm = [bmats(1), bmats(2), bmats(3);
    bmats(2), bmats(4), bmats(5);
    bmats(3), bmats(5), bmats(6)];
[V,D,W] = eig(bm);
printMat('bmat', bm);
fprintf('bvec = [%g %g %g];\n', bvec(1), bvec(2), bvec(3) );
fprintf('bval = %g;\n', bval );
fprintf('[V,D,W] = eig(bmat);\n')
printMat('D', D);
printMat('W', W);
printMat('V', V);
bvecEst= V(:,3);
fprintf('estimated_bvec = [%g %g %g];\n', bvecEst(1), bvecEst(2), bvecEst(3) );
deg = atan2d(norm(cross(bvec,bvecEst)),dot(bvec,bvecEst));
if deg > 1
    fprintf(2,'estimated_bvec_error_deg %g\n', deg);
end
fprintf('STEP2: estimate bmat from known bval/bvec;\n')
bmEst = bvec2bmatSiemens(bvec, bval);
printMat('estimated_bmat', bmEst);
%fprintf('max_bmat_error = %g\n', max(abs(bm(:)-bmEst(:))) );
%flip polarity of bvec
bvec = -bvec;
bmEst = bvec2bmatSiemens(bvec, bval);
printMat('estimated_bmat_flipped_bvec', bmEst);
fprintf('\n');
%end testBMat()

function bmEst = bvec2bmatSiemens(bvec, bval)
bmats = bval*[bvec(1).^2 bvec(1).*bvec(2) bvec(1).*bvec(3) bvec(2).^2 bvec(2).*bvec(3) bvec(3).^2];
bmEst = [bmats(1), bmats(2), bmats(3);
    bmats(2), bmats(4), bmats(5);
    bmats(3), bmats(5), bmats(6)];
%end bvec2bmatSiemens()

function bmEst = bvec2bmat(bvec, bval)
% https://cmiclab.cs.ucl.ac.uk/CMIC/NiftyFit-Release/blob/master/fit-lib/_dwiDti.cpp
% https://github.com/stijnimaging/ExploreDTI_scripts
% https://raw.githubusercontent.com/nipy/dipy/master/dipy/reconst/qtdmri.py (cite Basser et al. 1994)
bmats = bval*[bvec(1).^2 2*bvec(1).*bvec(2) 2*bvec(1).*bvec(3) bvec(2).^2 2*bvec(2).*bvec(3) bvec(3).^2];
bmEst = [bmats(1), bmats(2), bmats(3);
    bmats(2), bmats(4), bmats(5);
    bmats(3), bmats(5), bmats(6)];
%end bvec2bmat()

function printMat(str, m)
fprintf('%s=[%g %g %g; %g %g %g; %g %g %g];\n', str, ...
    m(1,1), m(1,2), m(1,3), ...
    m(2,1), m(2,2), m(2,3), ...
    m(3,1), m(3,2), m(3,3));
%end printMat()

@mharms
Copy link
Collaborator

mharms commented Feb 3, 2019

Regarding (7) above ( "Eddy undistortion will not be as effective"), I would argue that the problem is more severe than that, and that attempting to use eddy in the situation where the polarity of the bvecs is unknown would be outright inaccurate. I wonder if the 'bvecs' file in this situation shouldn't be named differently, so as to directly indicate the unknown polarity problem.

@tashrifbillah
Copy link
Contributor

Okay, so I shall try to come up with a mathematical formulation:

If we compare the 1st and 2nd expression for Ai , you should reach B=ggT which is also stated there. Then, the problem reduces to finding eigen decomposition of B in terms of g and gT. The math is as follows:

B= U A UT
where each column of U is a right eigenvector of B corresponding to each eigenvalue in the diagonal matrix A.

Then,
B= U A1/2 A1/2 UT
= (U A1/2) (U A1/2)T
= g gT

(since A is a diagonal eigenvalue matrix, its square root is straightforward)

If I follow my above math, then the first step would be computing eigenvalue of B, one of which is negative (surprising). Consequently, when you compute square root, you have a complex number.

B=[[78.9072,269.576,-25.1249],[269.576,269.576,920.972],[-25.1249,920.972,-85.8359]]
import numpy as np
lam, v= np.linalg.eig(i)
lam
array([1067.80058402,   79.56191972, -884.71520374])
# eigenvalues ordered greatest to smallest (opposite to MATLAB ordering)
lam= np.array([complex(x) for x in lam])
sq= np.sqrt(lam)
sq
array([32.67721812 +0 j,  8.91974886 +0 j, 0 +29.74416252j])

With complex eigenvalue, the resultant g has unique polarity, however complex.

g= v @ sq
g
array([14.91134518 +5.89484534j, 
24.91024819-19.01163331j,
17.44975139+22.10257453j])

I am trying to establish the way we can reach unequivocal gradient direction from B matrix. In this particular case, the non positive-semi-definiteness of B matrix, is giving us complex number, which should not be the case.

So, @isolovey-robarts , would it be possible to double check correctness of B matrix you provided?

@tashrifbillah
Copy link
Contributor

As an aside, I would love to know the theory about how the third eigenvector of a B matrix gives gradient direction.

@neurolabusc
Copy link
Collaborator

A simple example demonstrates that the B_matrix is underspecified to store the polarity of the B_vector. Consider the B_vector [1, 0,0] and its polar opposite [-1,0,0]: both create the B_matrix [1,0,0; 0, 0,0; 0,0,0]. Therefore, given only this B_matrix it is impossible to determine the polarity.

I agree with @tashrifbillah that using one vector of the SVD is unintuitive. The concept is described here. While this may not be the only solution, pragmatically it appears to work.

Note that many compute the B-Matrix by pre-multiplying the off diagonal entries by 2 (e.g. here, here and here). While not specified in the DICOM standard it seems that this is not done by Siemens or Bruker. In this way the DICOM B_matrix is similar to the NRRD format.

Jesper Andersson agrees with @mharms that we should not use the term "bvec" for files that can not determine polarity. Therefore, the latest commit will create "mvec" files, which have the same format as the "bvec" files but alert the user to the fact that the values come from the matrix and therefore the polarity is unknown. He also provides the following advice:

Normally when data is acquired on the half-shell we recommend using --slm=linear. That essentially does a linear regression of the estimated EC parameters on the x-, y- and z- components of the bvecs, and then mean-corrects those components. So instead of ending up in a “non-zero average distortions sitting somewhere inside the half-shell, is is extrapolated to the zero point in the Q-coordinate system. It is actually a lot simpler to do it than it is to explain it. What it sounds like with the Bruker data is that you don’t really know if you are on the half- or whole-shell, so those linear fits would fail. So one would therefore strongly advice to not use --slm=linear for the Bruker data.

This issue has been resolved to the full extent possible with the provided data. Validation would require images from a phantom or brain that exhibits tract-like diffusion signal. It would really help to have different series that validate axial, coronal and sagittal slices. Further, it woul be nice to see the slice angulation is handled correctly as the FSL bvec format reports the vectors with respect to the image, while the DICOM format is oriented relative to the participant. Without new data, I think we should close this issue.

@mharms
Copy link
Collaborator

mharms commented Feb 6, 2019

@tashrifbillah

As an aside, I would love to know the theory about how the third eigenvector of a B matrix gives gradient direction.

When I run one of the @neurolabusc examples above, the bvec (i.e., DiffusionGradientOrientation value, (0018,9089)) is the first (not third) eigenvector returned by Matlab's svd, which is what I would expect.

@mharms
Copy link
Collaborator

mharms commented Feb 6, 2019

Also, in this era in which accurate polarity of the bvec's matter, this is probably a good point for a historical reminder that prior to Siemens VB17 software line, the polarity of the bvec reported in the DICOM (DiffusionGradientOrientation value, (0018,9089)) was occasionally wrong. This is therefore an issue in VB15 and earlier data. And to make matters worse, the direction is itself sometimes simply wrong for VB13 and earlier data.

@neurolabusc
Copy link
Collaborator

Others have pointed out the bizarre nature of the B-Matrices provided in the Bruker DICOM header. One generally expects the B-Matrix to be a positive semidefinite having non-negative eigenvalues only. Further, the sum of the diagonals doesn’t equal the scalar b-value. The Matlab code below shows the B-matrices provided by @isolovey-robarts sample Bruker DICOM data with attempts to estimate b-vector and b-value from these. I am closing this issue until the issue is clarified. Until this is resolved, stable releases of dcm2niix will attempt to resolve the B-vector from the B-matrix. Any user who wants to explore this can work with the developmental branch code.

function brukerMat
%converts diffusion bmatrix to bvector, converts bvector+bvalue to bmatrix
% Rationale: 
%   Bruker only provides BMatrix+BVal, FSL expects BVec
%   https://github.com/rordenlab/dcm2niix/issues/265
%   Values as reported by Bruker BioSpec 94/30

bmats=[12.826 12.826 -13.8393 12.826 12.826 -13.8393];
bvec=[0.721174 0.692192 -0.0279017];
bval=40.6283;
testBMat(bmats, bvec, bval);

bmats=[212.982 -66.2239 -81.6229 -66.2239 20.7722 25.7217];
bvec=[0.916828 -0.208834 -0.340316];
bval=265.726;
testBMat(bmats, bvec, bval);

bmats=[68.8207 145.803 93.3369 145.803 308.898 197.743];
bvec=[0.335607 0.653464 0.678493];
bval=504.387;
testBMat(bmats, bvec, bval);

bmats=[18.564 -75.6026 108.966 -75.6026 316.969 -457.449];
bvec=[0.187981 -0.871057 -0.453787];
bval=995.805;
testBMat(bmats, bvec, bval);

bmats=[12.826 12.826 -13.8393 12.826 12.826 -13.8393];
bvec=[0.721174 0.692192 -0.0279017];
bval=40.6283;
testBMat(bmats, bvec, bval);

bmats=[17.5466 61.6668 -19.5919 61.6668 226.936 -69.6628];
bvec=[0.171555 0.798707 0.576745];
bval=266.465;
testBMat(bmats, bvec, bval);

bmats=[184.866 -306.533 -240.962 -306.533 508.275 399.548];
bvec=[0.476145 -0.453951 -0.753136];
bval=1007.3;
testBMat(bmats, bvec, bval);

bmats=[261.364 -207.139 136.949 -207.139 164.164 -108.537];
bvec=[0.93904 -0.296907 0.173352];
bval=497.369;
testBMat(bmats, bvec, bval);

bmats=[12.826 12.826 -13.8393 12.826 12.826 -13.8393];
bvec=[0.721174 0.692192 -0.0279017];
bval=40.6283;
testBMat(bmats, bvec, bval);

bmats=[28.031 20.5705 78.8185 20.5705 16.3182 58.0581];
bvec=[0.626191 0.230151 0.744926];
bval=266.056;
testBMat(bmats, bvec, bval);

bmats=[175.229 65.3131 -232.109 65.3131 24.3443 -86.5161];
bvec=[0.858803 0.175133 -0.481441];
bval=507.107;
testBMat(bmats, bvec, bval);

bmats=[806.061 368.024 165.275 368.024 168.029 75.4597];
bvec=[0.839318 0.495448 0.223778];
bval=1008.34;
testBMat(bmats, bvec, bval);

bmats=[12.826 12.826 -13.8393 12.826 12.826 -13.8393];
bvec=[0.721174 0.692192 -0.0279017];
bval=40.6283;
testBMat(bmats, bvec, bval);

bmats=[17.289 -51.2599 -35.5735 -51.2599 160.011 111.616];
bvec=[0.266501 -0.518918 -0.81222];
bval=255.24;
testBMat(bmats, bvec, bval);

bmats=[34.0121 -19.5627 126.818 -19.5627 14.0057 -73.2678];
bvec=[0.834704 -0.062723 0.547115];
bval=520.955;
testBMat(bmats, bvec, bval);

bmats=[725.612 -411.351 178.406 -411.351 233.195 -101.138];
bvec=[0.952174 -0.286507 0.106202];
bval=1002.75;
testBMat(bmats, bvec, bval);

bmats=[12.826 12.826 -13.8393 12.826 12.826 -13.8393];
bvec=[0.721174 0.692192 -0.0279017];
bval=40.6283;
testBMat(bmats, bvec, bval);

bmats=[41.2232 41.1216 -81.8416 41.1216 41.0204 -81.64];
bvec=[0.84358 0.450714 -0.29194];
bval=244.808;
testBMat(bmats, bvec, bval);

bmats=[473.904 109.475 -71.9371 109.475 25.2893 -16.6179];
bvec=[0.959963 0.253802 -0.118555];
bval=514.657;
testBMat(bmats, bvec, bval);

bmats=[78.9072 269.576 -25.1249 269.576 920.972 -85.8359];
bvec=[0.194068 0.768707 0.609448];
bval=1014.78;
testBMat(bmats, bvec, bval);
%end bruker()

function testBMat(bmats, bvec, bval)
fprintf('STEP1: estimate bvec from known bmat;\n')
bm = [bmats(1), bmats(2), bmats(3);
    bmats(2), bmats(4), bmats(5);
    bmats(3), bmats(5), bmats(6)];
[V,D,W] = eig(bm);
fprintf('bvec = [%g %g %g];\n', bvec(1), bvec(2), bvec(3) );
fprintf('bval = %g;\n', bval );
fprintf('[V,D,W] = eig(bmat);\n')
printMat('D', D);
printMat('W', W);
printMat('V', V);
bvecEst= V(:,3);
fprintf('estimated_bvec = [%g %g %g];\n', bvecEst(1), bvecEst(2), bvecEst(3) );
deg = atan2d(norm(cross(bvec,bvecEst)),dot(bvec,bvecEst));
if deg > 1
    fprintf(2,'estimated_bvec_error_deg %g\n', deg);
end
fprintf('STEP2: estimate bmat from known bval/bvec;\n')
bmEst = bvec2bmatSiemens(bvec, bval);
printMat('bmat', bm);
printMat('estimated_bmat', bmEst);
%fprintf('max_bmat_error = %g\n', max(abs(bm(:)-bmEst(:))) );
%flip polarity of bvec
bvec = -bvec;
bmEst = bvec2bmatSiemens(bvec, bval);
printMat('estimated_bmat_flipped_bvec', bmEst);
fprintf('estimated_bval=%g;\n', sum(diag(bm)));
%isPositiveDefinite(bm);
fprintf('\n');
%end testBMat()

function isPositiveDefinite(mat)
p = all(eig(mat) > eps);
if (~p)
    error('Not positive definite!')
end

function bmEst = bvec2bmatSiemens(bvec, bval)
%Like slicer, the B-matrix off-diagonal entries are NOT pre-multiplied by 2
% https://www.na-mic.org/wiki/NAMIC_Wiki:DTI:Nrrd_format
bmats = bval*[bvec(1).^2 bvec(1).*bvec(2) bvec(1).*bvec(3) bvec(2).^2 bvec(2).*bvec(3) bvec(3).^2];
bmEst = [bmats(1), bmats(2), bmats(3);
    bmats(2), bmats(4), bmats(5);
    bmats(3), bmats(5), bmats(6)];
%end bvec2bmatSiemens()

function printMat(str, m)
fprintf('%s=[%g %g %g; %g %g %g; %g %g %g];\n', str, ...
    m(1,1), m(1,2), m(1,3), ...
    m(2,1), m(2,2), m(2,3), ...
    m(3,1), m(3,2), m(3,3));
%end printMat()

xiangruili added a commit to xiangruili/dicm2nii that referenced this issue Feb 8, 2019
@xiangruili
Copy link

I tried the bvec related stuff. The private tag (0177,1101) stores a complete version of bmatrix (9 numbers each direction). Those 9 numbers seem to give correct bval. Based on this, my guess is that the 6-number bmatrix in PerFrameSeq may not be right.

Based on my guess, I updated dicm2nii.m. If there is a real dataset, please test the bvec to verify/deny my guess.

@neurolabusc
Copy link
Collaborator

@xiangruili you are correct. Looking at the VisuAcqDiffusionBMatrix stored in (0177,1101) reveals the error in the public tags. Specifically, the public tags are supposed to store [xx,xy,xz,yy,yz,zz], but Bruker has saved [xx,xy,xz,yx,yy,yz]. Since we have the b-value (0018,9087) and bvalue=xx+yy+zz, we can detect the error and recover zz. Your solution of reading the 0177,1101 works well, but is laborious to implement in C. Therefore, I have decided to rely on the public tags to detect this error and recover the value.

Be aware that my prior concern regarding the vectors [1,0,0] and [-1,0,0] creating the same BMatrix still holds: you can not determine BVector polarity from the BMatrix. Therefore, I would suggest you have dicm2nii use the extension ".mvec" instead of ".bvec" when the vectors are inferred from the matrix (as dcm2niix does).

@neurolabusc
Copy link
Collaborator

@xiangruili - I note that (0177,1101) also includes a field with DwDir that appears to be vectors with polarity and unit length. Better yet, they seem to report polarity. They are plausible on their own, but their direction seems in conflict with the values reported in the b-matrix. At first, I thought this might reflect a different frame of reference, but the slices are acquired aligned to the scanner bore (ImageOrientationPatient (0020,0037) DS [1\0\0\0\1\0]). Any thoughts about these?

##$PVM_DwDir=( 20, 3 )
1 0 0 -0.923015691400132 0.254004318110112 -0.289004913125285  -0.371098172012216 -0.786206934706051 0.494130349755311 0.106054844535604 
-0.563291296920237 -0.819423751647734 1 0 0 0.189954606273003 
0.952772314621957 0.23694337729843 0.429960230517931 -0.712934056649499 
0.553948762109148 0.722421713983863 -0.57254165515527 -0.387689205777878 1 0 
0 -0.334855039145283 -0.159930764964911 0.928598004077517 0.592151048844286 
0.220055549982116 0.775198484437812 0.897809235803934 0.409912902761262 
-0.160965798401374 1 0 0 -0.184009476732084 0.800041203182973 
-0.571029408771847 -0.260883045433759 0.0589730962176642 0.963567439532755 
-0.851828799615187 0.482902946260722 0.202959209297985 1 0 0 
-0.406076755760944 -0.405076566707345 -0.819154834897076 0.973810111526766 
0.224956099752581 0.0329942399872185 0.280904227984603 0.959672807349532 
-0.0109962509175467

@xiangruili
Copy link

xiangruili commented Feb 9, 2019

Using PVM_DwDir would avoid the sign issue from bmatrix, I think. There is another parameter PVM_DwGradVec. The two are quite similar with some difference. The latter has dim of (13, 3) for DTI_EPI_6dir_2bval, while the former has (6,3), omitting the copy for another bval. The latter is not normalized though. The vectors for bval of zero (or small bval) also show difference, but that would not matter. I think the latter is more convenient.

There are multiple parameters related to bmatrix:

PVM_DwBMat=( 20, 3, 3 )
PVM_DwBMatImag=( 20, 3, 3 )
PVM_DwBMatMag=( 20, 3, 3 )
PVM_DwBMatPat=( 20, 3, 3 )

We may not bother their meaning, but simply use PVM_DwGradVec or corrected bmatrix in public tag if the sign issue is figured out.

Users will need to test the bvec to make sure the sign is correct.

@neurolabusc neurolabusc reopened this Feb 9, 2019
@neurolabusc
Copy link
Collaborator

@isolovey-robarts can you confirm that your system has the ParaVison 6.0.1 patch ( PCISW-2290-2337-DCMExport.linux.tar.gz) installed. I you contact Bruker they can send you a copy and install instructions.

@xiangruili
Copy link

If Bruker agrees, the best solution is to store the real bvec into (0018,9089), following dicom standard. That will avoid all the confusion.

Another piece of information is the slice timing. I tried in the private tags, but with no success.

@neurolabusc
Copy link
Collaborator

I am going to close this issue. This is based on comments sent to me from an experienced Bruker user who noted they found bruker-dicom too inconsistent to start the conversion from the bruker-dicom. It may still be feasible, but I found that there were incongruencies between ParaVision version 5.1 and version 6 in the dicom as well. As I do not have this equipment, and I have only seen a limited set of Bruker DWI data (all of which has clear errors in the public DICOM tags), there is no way to reliably resolve this.

@neurolabusc
Copy link
Collaborator

neurolabusc commented Feb 15, 2023 via email

@nie-xingju
Copy link

thanks so much , the only thing I could find for VisuAcqDiffusionBMatrix is as below,

VisuAcqDiffusionBMatrix - Diffusion B-Matrix(es) in subject coordinate system.

The parameter VisuAcqDiffusionBMatrix is set to the PVM_DwBMat diffusion b-matrix and made dependent to the diffusion loop.

VisuAcqDiffusionBMatrix - Diffusion b-matrix(es) in subject coordinate system. First dimension is the number of b-matrices and second dimension is 9 for a 3x3 matrix.

@neurolabusc
Copy link
Collaborator

neurolabusc commented Feb 15, 2023

We can not use the b-matrix, as polar opposite b-vectors generate the identical b-matrix (e.g. consider the b-vectors [1,0,0] and [-1,0,0] which both generate the b-matrix [1,0,0; 0,1,0; 0,0,1]). We need to distinguish these for tools like Eddy to work effectively. In theory, the vector PVM_DwGradVec does encode this, and it looks like dicm2nii leverages this. However, it remains unclear to me if these vectors are in the frame of reference of the world or the image. In other words, if a scan is acquired with an angulation relative to the the scanner bore, what is the frame of reference for the vector. Perhaps @xiangruili has resolved this.

@xiangruili
Copy link

The private tag (0177,1100) contains PVM_DwGradRead, PVM_DwGradPhase, and PVM_DwGradSlice. Based on the names, I think those are in the frame image. For simplicity, I am using PVM_DwGradVec, which contains all 3 dimensions. This is based on the only testing dataset without angulation relative to the scanner bore. If someone can provide another dataset with angulation, we will be able to verify the solution.

@neurolabusc
Copy link
Collaborator

@nie-xingju can you provide two series of DTI DICOM images that use the same sequence, but where one is aligned orthogonal to the scanner bore and the other has an angulation applied? Ideally, a phantom or animal with anisotropic diffusion would help. This would test @xiangruili's hypothesis that vectors are applied in image space.

@nie-xingju
Copy link

sure, I could collect some data in the near future. My system is Bruker 94/20 with Paravision 360 v2.0 or v3.3

thanks,

@nie-xingju can you provide two series of DTI DICOM images that use the same sequence, but where one is aligned orthogonal to the scanner bore and the other has an angulation applied? Ideally, a phantom or animal with anisotropic diffusion would help. This would test @xiangruili's hypothesis that vectors are applied in image space.

@nie-xingju
Copy link

nie-xingju commented Mar 10, 2023

Hello Dr. Rorden,
sorry for the delay. I acquired some 30 direction DTI data with fixed mouse brain with Bruker default DTIEPI protocol.
Bruker 94/20 with Paravision 360 v2.0 , raw, dicom and nifti data are saved in each subfolder.

Width: 19.2000 mm (96)
Height: 19.2000 mm (96)
Depth: 157.5 mm (315)
Voxel size: 0.2x0.2x0.5 mm^3

all slice orientation: axial
#33: no angulation
#37: slice 5 degree angulation
#40: slice -5 degree angulation
#42: phase -175 degree angulation
#43: phase 175 degree angulation
#45: read 5 degree angulation
#46: read -5 degree angulation

best and please let me know if you need to test more

1.zip
33.zip
37.zip
40.zip
42.zip
43.zip
45.zip
46.zip

@neurolabusc
Copy link
Collaborator

  1. The rationale for the angulation is to decouple the scanner's world space frame of reference from the image space frame of reference. This is described here . Specifically, see figures 2 an image acquired with angulation and figure 3 for a correct and interpretation. The 5 and 175 degree angulation are far to close to the bore alignment to visualize errors. I think you want to acquire around 30 degrees orientation.

angulation

  1. Patient Position (0018,5100) reports Feet First Prone ('FFP'). Is this accurate? One of my concerns is the quadrupeds are often scanned without correctly specifying the orientation.

mouse

  1. I have been recently assigned a lot of new administrative responsibilities, and have limited time to support features that are outside my expertise. I do think a comprehensive solution will require not only validation images, but also painstaking work to validate the data. I think the best solution will be for engineers from the manufacturer to contribute, either with a stand alone solution or with pull requests to dicm2nii (Matlab) or dcm2niix (C).

@xiangruili
Copy link

xiangruili commented Mar 11, 2023

This is my observation from the new dataset, mainly from '45: read 5 degree angulation'.

bval and bvec are stored in the enhanced dicom format with public tags, which is nice, but see below for possible issue.

It seems the new dicom does not contain the private tag anymore. I found the file method contains those information seen in Private_0177_0011 in older dicom. Based on this file, DwGradVec is simply DwGradRead, DwGradPhase and DwGradSlice putting together, even with 5-deg angulation. If the Read/Phase/Slice names are meaningful, these suggest image reference as I guessed. Image reference is fine, although it is not as definitive as real world reference.

The problem is that those stored in the enhanced dicom are the same as DwGradVec, except normalized. Those dicom public tags should be in Patient reference, according to dicom standard.

There is a possibility that all of these bvec are in Patient reference, although the names suggest otherwise. This can be verified or denied with a larger angulation, as @neurolabusc has pointed out.

If the currently stored bvec are in Patient reference, we will need to fix the code, so the early data with private tag gives correct output. Otherwise, Bruker needs to fix the bvec in enhanced dicom, so they are in Patient reference.

EDIT on 3/13/2023:
Checked the scan "46: read -5 degree angulation", which has opposite rotation to scan 45. The method files give exactly the same DwGrad numbers for the 5 and -5 angulation. The same for both enhanced dicom files. I had thought I could figure out something, but realized the reference can be either way: either all are in Patient reference (then the names are misleading), or all are in image reference (then the dicom files have wrong values).

@neurolabusc
Copy link
Collaborator

@xiangruili well spotted. I also notice that the images with a low Diffusion b-value (0018 9087) omit the Diffusion Gradient Orientation (0018 9089). These are clearly intended as the B=0 images. However, DICOM convention is to record the b-vector for these images regardless of the low amplitude.

@neurolabusc
Copy link
Collaborator

These images are also very anisotropic (0.2mm in plane, 0.8mm slice thickness). For fitting diffusion tensors it helps to have isotropic data. Perhaps a smaller matrix in plane (64x64 instead of 96x96), decreased field of view and thinner slices would help.

@nie-xingju
Copy link

nie-xingju commented May 25, 2023

hello all, I apologize for the delay and collected more data with positive and minus 35 degree rotation in slice, phase and read direction.
Bruker 94/20 with Paravision 360 v2.0, raw, enhanced multi-frame DICOM data are saved in each subfolder.
30 direction DTI data with fixed mouse brain by Bruker 2D DTIEPI protocol.

head first, supine position if animal alive (data uploaded on March 10th is foot first prone)

Width: 19.2000 mm (64)
Height: 19.2000 mm (64)
Depth: 5.7 mm (19 slices)
Voxel size: 0.3x0.3x0.3 mm^3

all slice orientation: axial

3: no angulation
4: slice -35 degree angulation
5: slice 35 degree angulation
7: phase -145 degree angulation
8: phase 145 degree angulation
9: read -35 degree angulation
10: read 35 degree angulation

it seems to me the major changes are PVM_DwBMatMag and PVM_DwBMatPat after changing the angulation.

best and please let me know if you need to test more,
DTI_EPI_30dir_sat (E3) -no angulation.zip
DTI_EPI_30dir_read_35_degree (E10).zip
DTI_EPI_30dir_read_minus_35_degree (E9).zip
DTI_EPI_30dir_phase_145_degree (E8).zip
DTI_EPI_30dir_phase_minus_145_degree (E7).zip
DTI_EPI_30dir_slice_35_degree (E5).zip
DTI_EPI_30dir_slice_minus_35_degree (E4).zip

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

6 participants