In [1]:
%pylab inline

Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.

In [2]:
from glob import glob
import os
import shutil
from nipype.interfaces.base import CommandLine
from nipype.interfaces.freesurfer import MRIConvert
from nipype.pipeline.engine import Node
DWIConvert = '/software/DTIPrep/dev/DWIConvert'
dicomfile = 'dicoms/125000-5-1.dcm'
dicomtemplate = 'dicoms/125000-5-%d.dcm'

Convert dicom to nifti using mri_convert from freesurfer

In [3]:
dcm2nii = Node(MRIConvert(), name='dcm2nii')
dcm2nii.config['execution'] = {'remove_unnecessary_outputs': False}
dcm2nii.base_dir = os.path.abspath('temp')
dcm2nii.inputs.in_file = os.path.abspath(dicomfile)
dcm2nii.inputs.out_file = 'mri_convert.nii.gz'
res = dcm2nii.run()
for filename in glob(os.path.join(dcm2nii.base_dir, dcm2nii.name, 'mri_convert*')):
    shutil.copy(filename, os.path.abspath('data'))
130103-19:36:44,399 workflow INFO:
	 Executing node dcm2nii in dir: /Users/satra/Dropbox/shares/share_kent/temp/dcm2nii
130103-19:36:44,400 workflow INFO:
	 Collecting precomputed outputs

Convert dicom to nrrd using DWIConvert

In [4]:
dcm2nrrd = Node(CommandLine(DWIConvert), name='dcm2nrrd')
dcm2nrrd.config['execution'] = {'remove_unnecessary_outputs': False}
dcm2nrrd.base_dir = os.path.abspath('temp')
dcm2nrrd.inputs.args = ' '.join(['--inputDicomDirectory', os.path.abspath('dicoms'), '--outputVolume', 'dwi_dicom2nrrd.nrrd'])
res = dcm2nrrd.run()
for filename in glob(os.path.join(dcm2nrrd.base_dir, dcm2nrrd.name, '*nrrd')):
    shutil.copy(filename, os.path.abspath('data'))
130103-19:36:44,668 workflow INFO:
	 Executing node dcm2nrrd in dir: /Users/satra/Dropbox/shares/share_kent/temp/dcm2nrrd
130103-19:36:44,669 workflow INFO:
	 Collecting precomputed outputs

Convert nifti to nrrd using DWIConvert

In [5]:
nii2nrrd = Node(CommandLine(DWIConvert), name='nii2nrrd')
nii2nrrd.config['execution'] = {'remove_unnecessary_outputs': False}
nii2nrrd.base_dir = os.path.abspath('temp')
nii2nrrd.inputs.args = ' '.join(['--inputVolume', os.path.abspath('data/mri_convert.nii.gz'),
                                 '--inputBVectors', os.path.abspath('data/mri_convert.mghdti.bvecs'), 
                                 '--inputBValues', os.path.abspath('data/mri_convert.mghdti.bvals'), 
                                 '--outputVolume', 'dwi_nifti2nrrd.nrrd', '--conversionMode', 'FSLToNrrd'])
res = nii2nrrd.run()
for filename in glob(os.path.join(nii2nrrd.base_dir, nii2nrrd.name, '*nrrd')):
    shutil.copy(filename, os.path.abspath('data'))
130103-19:36:45,105 workflow INFO:
	 Executing node nii2nrrd in dir: /Users/satra/Dropbox/shares/share_kent/temp/nii2nrrd
130103-19:36:45,106 workflow INFO:
	 Collecting precomputed outputs

Convert nrrd from dicom to nifti using DWIConvert

In [6]:
nrrd2nii1 = Node(CommandLine(DWIConvert), name='nrrd2nii1')
nrrd2nii1.config['execution'] = {'remove_unnecessary_outputs': False}
nrrd2nii1.base_dir = os.path.abspath('temp')
nrrd2nii1.inputs.args = ' '.join(['--inputVolume', os.path.abspath('data/dwi_dicom2nrrd.nrrd'),
                                 '--outputBVectors', 'dwi_dcmnrrd2nii.bvecs', 
                                 '--outputBValues', 'dwi_dcmnrrd2nii.bvals', 
                                 '--outputVolume', 'dwi_dcmnrrd2nii.nii.gz', '--conversionMode', 'NrrdToFSL'])
res = nrrd2nii1.run()
for filename in glob(os.path.join(nrrd2nii1.base_dir, nrrd2nii1.name, 'dwi_dcmnrrd2nii*')):
    shutil.copy(filename, os.path.abspath('data'))
130103-19:36:45,582 workflow INFO:
	 Executing node nrrd2nii1 in dir: /Users/satra/Dropbox/shares/share_kent/temp/nrrd2nii1
130103-19:36:45,583 workflow INFO:
	 Collecting precomputed outputs

Convert nrrd from nifti back to nifti using DWIConvert

In [7]:
nrrd2nii2 = Node(CommandLine(DWIConvert), name='nrrd2nii2')
nrrd2nii2.config['execution'] = {'remove_unnecessary_outputs': False}
nrrd2nii2.base_dir = os.path.abspath('temp')
nrrd2nii2.inputs.args = ' '.join(['--inputVolume', os.path.abspath('data/dwi_nifti2nrrd.nrrd'),
                                  '--outputBVectors', 'dwi_niinrrd2nii.bvecs', 
                                  '--outputBValues', 'dwi_niinrrd2nii.bvals', 
                                  '--outputVolume', 'dwi_niinrrd2nii.nii.gz', '--conversionMode', 'NrrdToFSL'])
res = nrrd2nii2.run()
for filename in glob(os.path.join(nrrd2nii2.base_dir, nrrd2nii2.name, 'dwi_niinrrd2nii*')):
    shutil.copy(filename, os.path.abspath('data'))
130103-19:36:45,824 workflow INFO:
	 Executing node nrrd2nii2 in dir: /Users/satra/Dropbox/shares/share_kent/temp/nrrd2nii2
130103-19:36:45,825 workflow INFO:
	 Collecting precomputed outputs

In [8]:
import nibabel as nb
import numpy as np
np.set_printoptions(precision=8, suppress=True)

Geometry of converted nifti file

In [9]:
orignii = nb.load('data/mri_convert.nii.gz')
dcmnrrdnii = nb.load('data/dwi_dcmnrrd2nii.nii.gz')
niinrrdnii = nb.load('data/dwi_niinrrd2nii.nii.gz')
mricronnii = nb.load('data/20121221_145258DIFFUSIONHighRess005a001.nii.gz')

dicom to nifti: mri_convert

In [10]:
print orignii.get_affine()
print orignii.get_header().get_zooms()
[[  -1.99463236    0.11752418    0.08735005  119.13967133]
 [  -0.12566817   -1.98628283   -0.19720104  145.03933716]
 [   0.07516297   -0.20216042    1.98833561  -34.36771393]
 [   0.            0.            0.            1.        ]]
(2.0, 2.0, 1.9999992, 8.04)

dicom to nrrd to nifti: dwi_convert

In [11]:
print dcmnrrdnii.get_affine()
print dcmnrrdnii.get_header().get_zooms()
[[  -1.99808002   -0.02362024    0.08426505  126.69300079]
 [   0.008015     -1.96706998   -0.36133584  136.94200134]
 [   0.08714505   -0.36065215    1.96528006  -45.38460159]
 [   0.            0.            0.            1.        ]]
(1.9999956, 1.9999981, 1.9999975, 1.0)

dicom to nifti to nrrd to nifti: mri_convert + dwi_convert

In [12]:
print niinrrdnii.get_affine()
print niinrrdnii.get_header().get_zooms()
[[  -1.99807489   -0.00984556    0.04231431  124.40200043]
 [   0.00600271   -1.96706533   -0.18054311  142.48399353]
 [   0.04357169   -0.18049444    1.96528769  -48.8185997 ]
 [   0.            0.            0.            1.        ]]
(1.998559, 1.9753535, 1.9740167, 1.0)

dicom to nifti: mricron

In [13]:
print mricronnii.get_affine()
print mricronnii.get_header().get_zooms()
[[  -1.99463236   -0.11752418    0.08734999  134.06520081]
 [  -0.12566817    1.98628283   -0.19720103 -107.21861267]
 [   0.07516297    0.20216042    1.98833573  -60.04209137]
 [   0.            0.            0.            1.        ]]
(2.0, 2.0, 1.9999992, 8.04)

Gradient vectors from the various converters

In [14]:
dcm2nrrd2nii_bvecs = np.genfromtxt('data/dwi_dcmnrrd2nii.bvecs')
freesurfbvecs = np.genfromtxt('data/mri_convert.mghdti.bvecs')
mricronbvecs = np.genfromtxt('data/20121221_145258DIFFUSIONHighRess005a001.bvec')
In [15]:
def plot_angles(vector1, vector2):
    plot([vector1[idx, :].dot(vector2[idx, :].T) for idx in range(vector1.shape[0])])
    plot([vector1[idx, :].dot(vector1[idx, :].T) for idx in range(vector1.shape[0])], ':')
    plot([vector2[idx, :].dot(vector2[idx, :].T) for idx in range(vector1.shape[0])], '.-')
In [16]:
subplot(311);plot_angles(dcm2nrrd2nii_bvecs, freesurfbvecs)
subplot(312);plot_angles(dcm2nrrd2nii_bvecs, mricronbvecs.T)
subplot(313);plot_angles(freesurfbvecs, mricronbvecs.T)

from dicom -> nrrd -> nifti via DWIConvert

In [17]:
dcm2nrrd2nii_bvecs[10:14, ...]
Out[17]:
array([[-0.999965  , -0.00722049, -0.00433271],
       [-0.588614  , -0.382054  , -0.712439  ],
       [ 0.896099  , -0.0168485 , -0.443534  ],
       [ 0.72873   ,  0.418467  , -0.542067  ]])

from dicom to nifti via mri_convert

In [18]:
freesurfbvecs[10:14, ...]
Out[18]:
array([[ 1.      ,  0.      ,  0.      ],
       [ 0.586792, -0.376583, -0.716841],
       [ 0.894195,  0.022917,  0.447091],
       [ 0.729281, -0.41097 ,  0.547041]])

from dicom to nifti via mricron

In [19]:
mricronbvecs[..., 10:14].T
Out[19]:
array([[ 0.99718219,  0.06522338, -0.03706211],
       [ 0.63686717, -0.27377155,  0.72072834],
       [ 0.87571853,  0.03634979, -0.48145169],
       [ 0.73273391, -0.42719257, -0.52972406]])
In [20]:
from nibabel.nicom import dicomwrappers as dw
-c:1: UserWarning: The DICOM readers are highly experimental, unstable, and only work for Siemens time-series at the moment
Please use with caution.  We would be grateful for your help in improving them

In [21]:
from nibabel.nicom import dicomreaders as dr
In [22]:
data, aff, bvals, bvecs = dr.read_mosaic_dwi_dir('dicoms', force=True)
/software/nipy-repo/nibabel/nibabel/nicom/dicomreaders.py:101: RuntimeWarning: invalid value encountered in divide
  g = q / b

from dicom to numpy via nibabel

In [23]:
bvecs[10:14, ...]
Out[23]:
array([[ 0.98130271,  0.14653902, -0.12478504],
       [        nan,         nan,         nan],
       [ 0.06472833,  0.09414854, -0.99345171],
       [ 0.69413   , -0.06984238,  0.71645348]])

from dicom via nibabel

In [24]:
for i in range(11, 15):
    print dw.csar.get_g_vector(dw.wrapper_from_file(dicomtemplate % i, force=True).csa_header)
[ 0.99995732 -0.00577336  0.00722241]
[ 0.58759272  0.38297507  0.71278679]
[ 0.8965317  -0.02854692 -0.44205895]
[ 0.7288003   0.41807222 -0.54227829]

In [25]:
for i in range(11, 15):
    nw = dw.wrapper_from_file(dicomtemplate % i, force=True)
    print nw.csa_header['tags']['DiffusionGradientDirection']['items']
[0.99995732, -0.00577336, 0.00722241]
[0.58759272, 0.38297507, 0.71278679]
[0.8965317, -0.02854692, -0.44205895]
[0.7288003, 0.41807222, -0.54227829]

In [25]: