Merge pull request #21 from Garyfallidis/dsi-with-b0s

BF: make dsi to work with data with multiple b0s etc.

dipy/reconst/dsi.py

@@ -114,8 +114,9 @@ def __init__(self,
         self.qgrid_size = qgrid_size
         # necessary shifting for centering
         self.origin = self.qgrid_size // 2
+
         # hanning filter width
-        self.filter = hanning_filter(gtab, filter_width)
+        self.filter = hanning_filter(gtab, filter_width, self.origin)
         # odf sampling radius
         self.qradius = np.arange(r_start, r_end, r_step)
         self.qradiusn = len(self.qradius)
@@ -159,6 +160,7 @@ def pdf(self, normalized=True):
         # create the signal volume
         Sq = np.zeros((self.qgrid_sz, self.qgrid_sz, self.qgrid_sz))
         # fill q-space
+
         for i in range(len(values)):
             qx, qy, qz = self.model.qgrid[i]
             Sq[qx, qy, qz] += values[i]
@@ -222,7 +224,7 @@ def rtop_pdf(self, normalized=True):
          PhD Thesis, 2002.
 
         .. [3] Wu Y. et. al, "Computation of Diffusion Function Measures
-        in q -Space Using Magnetic Resonance Hybrid Diffusion Imaging", 
+        in q -Space Using Magnetic Resonance Hybrid Diffusion Imaging",
         IEEE TRANSACTIONS ON MEDICAL IMAGING, vol. 27, No. 6, p. 858-865, 2008
 
         """
@@ -309,31 +311,53 @@ def create_qspace(gtab, origin):
     ----------
     gtab : GradientTable
     origin : (3,) ndarray
-        center of the qspace
+        center of qspace
 
     Returns
     -------
     qgrid : ndarray
         qspace coordinates
     """
     # create the q-table from bvecs and bvals
-    qtable = create_qtable(gtab)
+    qtable = create_qtable(gtab, origin)
+
     # center and index in qspace volume
     qgrid = qtable + origin
     return qgrid.astype('i8')
 
 
-def create_qtable(gtab):
+def create_qtable(gtab, origin):
     """ create a normalized version of gradients
+
+    Parameters
+    ----------
+    gtab : GradientTable
+    origin : (3,) ndarray
+        center of qspace
+
+    Returns
+    -------
+    qtable : ndarray
     """
+
     bv = gtab.bvals
-    bmin = np.sort(bv)[1]
+    bsorted = np.sort(bv[np.bitwise_not(gtab.b0s_mask)])
+    for i in range(len(bsorted)):
+        bmin = bsorted[i]
+        try:
+            if bv.max() / bmin > origin * 2:
+                continue
+            else:
+                break
+        except ZeroDivisionError:
+            continue
+
     bv = np.sqrt(bv / bmin)
     qtable = np.vstack((bv, bv, bv)).T * gtab.bvecs
     return np.floor(qtable + .5)
 
 
-def hanning_filter(gtab, filter_width):
+def hanning_filter(gtab, filter_width, origin):
     """ create a hanning window
 
     The signal is premultiplied by a Hanning window before
@@ -344,14 +368,16 @@ def hanning_filter(gtab, filter_width):
     ----------
     gtab : GradientTable
     filter_width : int
+    origin : (3,) ndarray
+        center of qspace
 
     Returns
     -------
     filter : (N,) ndarray
         where N is the number of non-b0 gradient directions
 
     """
-    qtable = create_qtable(gtab)
+    qtable = create_qtable(gtab, origin)
     # calculate r - hanning filter free parameter
     r = np.sqrt(qtable[:, 0] ** 2 + qtable[:, 1] ** 2 + qtable[:, 2] ** 2)
     # setting hanning filter width and hanning

dipy/reconst/tests/test_dsi.py

@@ -20,6 +20,7 @@
 def test_dsi():
     #load symmetric 724 sphere
     sphere = get_sphere('symmetric724')
+
     #load icosahedron sphere
     sphere2 = create_unit_sphere(5)
     btable = np.loadtxt(get_data('dsi515btable'))
@@ -30,18 +31,19 @@ def test_dsi():
 
     ds = DiffusionSpectrumModel(gtab)
 
-    #symmetric724
+    # symmetric724
     dsfit = ds.fit(data)
     odf = dsfit.odf(sphere)
-    directions, _, _ = peak_directions(odf, sphere, .35, 25)
+
+    directions, _, _ = peak_directions(odf, sphere)
     assert_equal(len(directions), 2)
     assert_almost_equal(angular_similarity(directions, golden_directions),
                         2, 1)
 
-    #5 subdivisions
+    # 5 subdivisions
     dsfit = ds.fit(data)
     odf2 = dsfit.odf(sphere2)
-    directions, _, _ = peak_directions(odf2, sphere2, .35, 25)
+    directions, _, _ = peak_directions(odf2, sphere2)
     assert_equal(len(directions), 2)
     assert_almost_equal(angular_similarity(directions, golden_directions),
                         2, 1)
@@ -51,7 +53,7 @@ def test_dsi():
     for sbd in sb_dummies:
         data, golden_directions = sb_dummies[sbd]
         odf = ds.fit(data).odf(sphere2)
-        directions, _, _ = peak_directions(odf, sphere2, .35, 25)
+        directions, _, _ = peak_directions(odf, sphere2)
         if len(directions) <= 3:
             assert_equal(len(directions), len(golden_directions))
         if len(directions) > 3:
@@ -78,13 +80,18 @@ def test_multib0_dsi():
     new_bvecs = np.concatenate([gtab.bvecs, np.zeros((1, 3))])
     new_bvals = np.concatenate([gtab.bvals, [0]])
     new_gtab = gradient_table(new_bvals, new_bvecs)
-    DS = DiffusionSpectrumModel(new_gtab)
-    sphere = get_sphere('symmetric724')
-    DSfit = DS.fit(new_data)
-    PDF = DSfit.pdf()
-    assert_equal(new_data.shape[:-1] + (17, 17, 17), PDF.shape)
-    assert_equal(np.alltrue(np.isreal(PDF)), True)
+    ds = DiffusionSpectrumModel(new_gtab)
+    sphere = get_sphere('repulsion724')
+    dsfit = ds.fit(new_data)
+    pdf = dsfit.pdf()
+    odf = dsfit.odf(sphere)
+    assert_equal(new_data.shape[:-1] + (17, 17, 17), pdf.shape)
+    assert_equal(np.alltrue(np.isreal(pdf)), True)
+    from dipy.viz import fvtk
 
+    ren = fvtk.ren()
+    fvtk.add(ren, fvtk.sphere_funcs(odf, sphere))
+    fvtk.show(ren)
 
 def sticks_and_ball_dummies(gtab):
     sb_dummies={}