In [1]:

    

import numpy as np
import nibabel as nib
import dsiadapt as dsi
from dipy.viz import fvtk
import dipy.core.gradients as grad
from dipy.data import get_sphere

%pylab inline
np.set_printoptions(threshold=numpy.nan)


    

    




Populating the interactive namespace from numpy and matplotlib

In [2]:

    

# ODF sphere directions
sphere = get_sphere('symmetric724')
sphere = sphere.subdivide(2)


    

In [3]:

    

fnArr = np.array(['DSI11_exvivo', 'DSI11_exvivo', 'DSI11_invivo_b10k', 'DSI11_invivo_b7k'])
pdfwgtArr = np.array([8, 2, 2, 2]);

qgridsz = 257; # Q-space grid size
qgridcenter = qgridsz//2;

for ii in np.arange(fnArr.shape[0]):
    fn = fnArr[ii];
    pdfwgt = pdfwgtArr[ii];
    data = nib.load('data/' + fn + '_roi.nii.gz').get_data();
    gtab = grad.gradient_table('data/' + fn + '_bvals.txt', 'data/' + fn + '_bvecs_dipy.txt');
    mdd = np.loadtxt('data/' + fn + '_stats.txt')[0];
    fov = np.loadtxt('data/' + fn + '_stats.txt')[1];

    rend = mdd / fov * qgridsz;
    dsimodel = dsi.DiffusionSpectrumModel(gtab, qgrid_size=qgridsz, filter_width=np.inf, 
               filter_type='none', pdfwgt=pdfwgt, r_start=0, r_end=rend, r_step=0.2);
    dsifit = dsimodel.fit(data);
    odf = dsifit.odf(sphere);

    # Save ODF
    r = fvtk.ren();
    fvtk.add(r, fvtk.sphere_funcs(odf, sphere));
    fvtk.camera(r, pos=(0, 1, 0), viewup=(0, 0, 1))
    odfname = 'figure5/' + fn + '_roi_odf_wgt_' + str(pdfwgt) + '.png';
    fvtk.record(r, n_frames=1, out_path=odfname, size=(1600,1000))
    #fvtk.show(r)


    

    




Camera Position (0.00,0.00,1.00)
Camera Focal Point (0.00,0.00,0.00)
Camera View Up (0.00,1.00,0.00)
-------------------------------------
Camera New Position (0.00,1.00,0.00)
Camera New Focal Point (0.00,0.00,0.00)
Camera New View Up (0.00,0.00,1.00)






    




---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-3-9fe006687360> in <module>()
     17                filter_type='none', pdfwgt=pdfwgt, r_start=0, r_end=rend, r_step=0.2);
     18     dsifit = dsimodel.fit(data);
---> 19     odf = dsifit.odf(sphere);
     20 
     21     # Save ODF

/usr/local/anaconda/lib/python2.7/site-packages/dipy/reconst/multi_voxel.pyc in __call__(self, *args, **kwargs)
     97             item = self[ijk]
     98             if item is not None:
---> 99                 result[ijk] = item(*args, **kwargs)
    100         return _squash(result)

/data/qytian/QTResearchResult/DSIQspace/DSIQspace/dsiadapt.pyc in odf(self, sphere)
    317             self.model.cache_set('interp_coords', sphere, interp_coords)
    318 
--> 319         Pr = self.pdf()
    320 
    321         # calculate the orientation distribution function

/data/qytian/QTResearchResult/DSIQspace/DSIQspace/dsiadapt.pyc in pdf(self, normalized, clipped)
    168         # apply fourier transform
    169         Pr = fftshift(np.real(fftn(ifftshift(Sq),
--> 170                                    3 * (self.qgrid_sz, ))))
    171         # clipping negative values to 0 (ringing artefact)
    172         if clipped:

/usr/local/anaconda/lib/python2.7/site-packages/scipy/fftpack/basic.pyc in fftn(x, shape, axes, overwrite_x)
    591 
    592     """
--> 593     return _raw_fftn_dispatch(x, shape, axes, overwrite_x, 1)
    594 
    595 

/usr/local/anaconda/lib/python2.7/site-packages/scipy/fftpack/basic.pyc in _raw_fftn_dispatch(x, shape, axes, overwrite_x, direction)
    606 
    607     overwrite_x = overwrite_x or _datacopied(tmp, x)
--> 608     return _raw_fftnd(tmp,shape,axes,direction,overwrite_x,work_function)
    609 
    610 

/usr/local/anaconda/lib/python2.7/site-packages/scipy/fftpack/basic.pyc in _raw_fftnd(x, s, axes, direction, overwrite_x, work_function)
    510             x, copy_made = _fix_shape(x, s[i], i)
    511             overwrite_x = overwrite_x or copy_made
--> 512         return work_function(x,s,direction,overwrite_x=overwrite_x)
    513 
    514     # We ordered axes, because the code below to push axes at the end of the

KeyboardInterrupt: