Merge pull request dipy#1212 from aaossa/master

Follow PEP8 in reconst (part 2)

dipy/reconst/dsi.py

@@ -180,8 +180,8 @@ def rtop_signal(self, filtering=True):
 
         Parameters
         ----------
-        filtering : boolean
-            default true, perform the hanning filtering
+        filtering : boolean, optional
+            Whether to perform Hanning filtering. Default: True
 
         Returns
         -------
@@ -200,13 +200,15 @@ def rtop_signal(self, filtering=True):
 
     def rtop_pdf(self, normalized=True):
         r""" Calculates the return to origin probability from the propagator, which is
-        the propagator evaluated at zero (see Descoteaux et Al. [1]_, Tuch [2]_, Wu et al. [3]_)
+        the propagator evaluated at zero (see Descoteaux et Al. [1]_,
+        Tuch [2]_, Wu et al. [3]_)
         rtop = P(0)
 
         Parameters
         ----------
-        normalized : boolean
-            default true, normalize the propagator by its sum in order to obtain a pdf
+        normalized : boolean, optional
+            Whether to normalize the propagator by its sum in order to obtain a
+            pdf. Default: True.
 
         Returns
         -------
@@ -243,12 +245,14 @@ def msd_discrete(self, normalized=True):
                     MSD:{DSI}=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} P(\hat{\mathbf{r}}) \cdot \hat{\mathbf{r}}^{2} \ dr_x \ dr_y \ dr_z
                 \end{equation}
 
-        where $\hat{\mathbf{r}}$ is a point in the 3D Propagator space (see Wu et. al [1]_).
+        where $\hat{\mathbf{r}}$ is a point in the 3D Propagator space
+        (see Wu et. al [1]_).
 
         Parameters
         ----------
-        normalized : boolean
-            default true, normalize the propagator by its sum in order to obtain a pdf
+        normalized : boolean, optional
+            Whether to normalize the propagator by its sum in order to obtain a
+            pdf. Default: True
 
         Returns
         -------
@@ -573,7 +577,7 @@ def pdf(self):
         self.model.cache_set('deconv_psf', self.model.gtab, DSID_PSF)
         # apply fourier transform
         Pr = fftshift(np.abs(np.real(fftn(ifftshift(Sq),
-                      3 * (self.qgrid_sz, )))))
+                                          3 * (self.qgrid_sz, )))))
         # threshold propagator
         Pr = threshold_propagator(Pr)
         # apply LR deconvolution
@@ -645,8 +649,9 @@ def LR_deconv(prop, psf, numit=5, acc_factor=1):
         reBlurred = np.real(np.fft.ifftn(otf * np.fft.fftn(prop_deconv)))
         reBlurred[reBlurred < eps] = eps
         # Update the estimate
-        prop_deconv = prop_deconv * (np.real(np.fft.ifftn(otf *
-                                                          np.fft.fftn((prop / reBlurred) + eps)))) ** acc_factor
+        prop_deconv = prop_deconv * (
+            np.real(np.fft.ifftn(
+                otf * np.fft.fftn((prop / reBlurred) + eps)))) ** acc_factor
         # Enforce positivity
         prop_deconv = np.clip(prop_deconv, 0, np.inf)
     return prop_deconv / prop_deconv.sum()

dipy/reconst/dti.py

@@ -1373,12 +1373,12 @@ def wls_fit_tensor(design_matrix, data, return_S0_hat=False):
     log_s = np.log(data)
     w = np.exp(np.einsum('...ij,...j', ols_fit, log_s))
     fit_result = np.einsum('...ij,...j',
-                  pinv(design_matrix * w[..., None]),
-                  w * log_s)
+                           pinv(design_matrix * w[..., None]),
+                           w * log_s)
     if return_S0_hat:
         return (eig_from_lo_tri(fit_result,
                                 min_diffusivity=tol / -design_matrix.min()),
-                                np.exp(-fit_result[:, -1]))
+                np.exp(-fit_result[:, -1]))
     else:
         return eig_from_lo_tri(fit_result,
                                min_diffusivity=tol / -design_matrix.min())
@@ -1437,7 +1437,7 @@ def ols_fit_tensor(design_matrix, data, return_S0_hat=False):
     if return_S0_hat:
         return (eig_from_lo_tri(fit_result,
                                 min_diffusivity=tol / -design_matrix.min()),
-                                np.exp(-fit_result[:, -1]))
+                np.exp(-fit_result[:, -1]))
     else:
         return eig_from_lo_tri(fit_result,
                                min_diffusivity=tol / -design_matrix.min())

dipy/reconst/fwdti.py

@@ -638,8 +638,8 @@ def nls_iter(design_matrix, sig, S0, Diso=3e-3, mdreg=2.7e-3,
         # Process water volume fraction f
         f = this_tensor[7]
         if f_transform:
-            f = 0.5 * (1 + np.sin(f - np.pi/2))        
-        
+            f = 0.5 * (1 + np.sin(f - np.pi/2))
+
         params = np.concatenate((evals, evecs[0], evecs[1], evecs[2],
                                  np.array([f])), axis=0)
     return params

dipy/reconst/gqi.py

@@ -37,13 +37,13 @@ def __init__(self,
 
         .. [2] Garyfallidis E, "Towards an accurate brain tractography", PhD
         thesis, University of Cambridge, 2012.
-        
+
         Notes
         -----
         As of version 0.9, range of the sampling length in GQI2 has changed
-	to match the same scale used in the 'standard' method [1]_. This 
-        means that the value of `sampling_length` should be approximately 
-        1 - 1.3 (see [1]_, pg. 1628). 
+        to match the same scale used in the 'standard' method [1]_. This
+        means that the value of `sampling_length` should be approximately
+        1 - 1.3 (see [1]_, pg. 1628).
 
         Examples
         --------
@@ -53,7 +53,7 @@ def __init__(self,
 
         >>> from dipy.data import dsi_voxels
         >>> data, gtab = dsi_voxels()
-        >>> from dipy.core.subdivide_octahedron import create_unit_sphere 
+        >>> from dipy.core.subdivide_octahedron import create_unit_sphere
         >>> sphere = create_unit_sphere(5)
         >>> from dipy.reconst.gqi import GeneralizedQSamplingModel
         >>> gq = GeneralizedQSamplingModel(gtab, 'gqi2', 1.1)
@@ -75,7 +75,7 @@ def __init__(self,
         scaling = np.sqrt(self.gtab.bvals * 0.01506)
         tmp = np.tile(scaling, (3, 1))
         gradsT = self.gtab.bvecs.T
-        b_vector = gradsT * tmp # element-wise product
+        b_vector = gradsT * tmp  # element-wise product
         self.b_vector = b_vector.T
 
     @multi_voxel_fit
@@ -109,18 +109,22 @@ def odf(self, sphere):
         self.gqi_vector = self.model.cache_get('gqi_vector', key=sphere)
         if self.gqi_vector is None:
             if self.model.method == 'gqi2':
-                H=squared_radial_component
-                #print self.gqi_vector.shape
-                self.gqi_vector = np.real(H(np.dot(self.model.b_vector,                                         sphere.vertices.T) * self.model.Lambda))
+                H = squared_radial_component
+                # print self.gqi_vector.shape
+                self.gqi_vector = np.real(H(np.dot(
+                    self.model.b_vector, sphere.vertices.T) *
+                    self.model.Lambda))
             if self.model.method == 'standard':
-                self.gqi_vector = np.real(np.sinc(np.dot(self.model.b_vector,                                   sphere.vertices.T) * self.model.Lambda / np.pi))
+                self.gqi_vector = np.real(np.sinc(np.dot(
+                    self.model.b_vector, sphere.vertices.T) *
+                    self.model.Lambda / np.pi))
             self.model.cache_set('gqi_vector', sphere, self.gqi_vector)
 
         return np.dot(self.data, self.gqi_vector)
 
 
 def normalize_qa(qa, max_qa=None):
-    """ Normalize quantitative anisotropy. 
+    """ Normalize quantitative anisotropy.
 
     Used mostly with GQI rather than GQI2.
 
@@ -139,8 +143,8 @@ def normalize_qa(qa, max_qa=None):
     Notes
     -----
     Normalized quantitative anisotropy has the very useful property
-    to be very small near gray matter and background areas. Therefore, 
-    it can be used to mask out white matter areas. 
+    to be very small near gray matter and background areas. Therefore,
+    it can be used to mask out white matter areas.
 
     """
     if max_qa is None:
@@ -165,31 +169,39 @@ def npa(self, odf, width=5):
 
     Nimmo-Smith et. al  ISMRM 2011
     """
-    #odf = self.odf(s)
+    # odf = self.odf(s)
     t0, t1, t2 = triple_odf_maxima(self.odf_vertices, odf, width)
     psi0 = t0[1] ** 2
     psi1 = t1[1] ** 2
     psi2 = t2[1] ** 2
-    npa = np.sqrt((psi0 - psi1) ** 2 + (psi1 - psi2) ** 2 + (psi2 - psi0) ** 2) / np.sqrt(2 * (psi0 ** 2 + psi1 ** 2 + psi2 ** 2))
-    #print 'tom >>>> ',t0,t1,t2,npa
+    npa = (np.sqrt(
+        (psi0 - psi1) ** 2 +
+        (psi1 - psi2) ** 2 +
+        (psi2 - psi0) ** 2) /
+        np.sqrt(2 * (psi0 ** 2 + psi1 ** 2 + psi2 ** 2)))
+    # print 'tom >>>> ',t0,t1,t2,npa
 
-    return t0,t1,t2,npa
+    return t0, t1, t2, npa
 
 
 def equatorial_zone_vertices(vertices, pole, width=5):
     """
     finds the 'vertices' in the equatorial zone conjugate
     to 'pole' with width half 'width' degrees
     """
-    return [i for i,v in enumerate(vertices) if np.abs(np.dot(v,pole)) < np.abs(np.sin(np.pi*width/180))]
+    return [i
+            for i, v in enumerate(vertices)
+            if np.abs(np.dot(v, pole)) < np.abs(np.sin(np.pi * width / 180))]
 
 
 def polar_zone_vertices(vertices, pole, width=5):
     """
     finds the 'vertices' in the equatorial band around
     the 'pole' of radius 'width' degrees
     """
-    return [i for i,v in enumerate(vertices) if np.abs(np.dot(v,pole)) > np.abs(np.cos(np.pi*width/180))]
+    return [i
+            for i, v in enumerate(vertices)
+            if np.abs(np.dot(v, pole)) > np.abs(np.cos(np.pi * width / 180))]
 
 
 def upper_hemi_map(v):
@@ -201,9 +213,10 @@ def upper_hemi_map(v):
 
 def equatorial_maximum(vertices, odf, pole, width):
     eqvert = equatorial_zone_vertices(vertices, pole, width)
-    #need to test for whether eqvert is empty or not
+    # need to test for whether eqvert is empty or not
     if len(eqvert) == 0:
-        print('empty equatorial band at %s  pole with width %f' % (np.array_str(pole), width))
+        print('empty equatorial band at %s  pole with width %f' %
+              (np.array_str(pole), width))
         return None, None
     eqvals = [odf[i] for i in eqvert]
     eqargmax = np.argmax(eqvals)
@@ -213,18 +226,21 @@ def equatorial_maximum(vertices, odf, pole, width):
     return eqvertmax, eqvalmax
 
 
-def patch_vertices(vertices,pole, width):
+def patch_vertices(vertices, pole, width):
     """
     find 'vertices' within the cone of 'width' degrees around 'pole'
     """
-    return [i for i,v in enumerate(vertices) if np.abs(np.dot(v,pole)) > np.abs(np.cos(np.pi*width/180))]
+    return [i
+            for i, v in enumerate(vertices)
+            if np.abs(np.dot(v, pole)) > np.abs(np.cos(np.pi * width / 180))]
 
 
 def patch_maximum(vertices, odf, pole, width):
     eqvert = patch_vertices(vertices, pole, width)
-    #need to test for whether eqvert is empty or not
+    # need to test for whether eqvert is empty or not
     if len(eqvert) == 0:
-        print('empty cone around pole %s with with width %f' % (np.array_str(pole), width))
+        print('empty cone around pole %s with with width %f' %
+              (np.array_str(pole), width))
         return np.Null, np.Null
     eqvals = [odf[i] for i in eqvert]
     eqargmax = np.argmax(eqvals)
@@ -239,26 +255,27 @@ def odf_sum(odf):
 
 def patch_sum(vertices, odf, pole, width):
     eqvert = patch_vertices(vertices, pole, width)
-    #need to test for whether eqvert is empty or not
+    # need to test for whether eqvert is empty or not
     if len(eqvert) == 0:
-        print('empty cone around pole %s with with width %f' % (np.array_str(pole), width))
+        print('empty cone around pole %s with with width %f' %
+              (np.array_str(pole), width))
         return np.Null
     return np.sum([odf[i] for i in eqvert])
 
 
 def triple_odf_maxima(vertices, odf, width):
 
-    indmax1 = np.argmax([odf[i] for i,v in enumerate(vertices)])
+    indmax1 = np.argmax([odf[i] for i, v in enumerate(vertices)])
     odfmax1 = odf[indmax1]
     pole = vertices[indmax1]
     eqvert = equatorial_zone_vertices(vertices, pole, width)
-    indmax2, odfmax2 = equatorial_maximum(vertices,\
-                                              odf, pole, width)
-    indmax3 = eqvert[np.argmin([np.abs(np.dot(vertices[indmax2],vertices[p])) for p in eqvert])]
+    indmax2, odfmax2 = equatorial_maximum(vertices, odf, pole, width)
+    indmax3 = eqvert[np.argmin([np.abs(np.dot(vertices[indmax2], vertices[p]))
+                                for p in eqvert])]
     odfmax3 = odf[indmax3]
     """
     cross12 = np.cross(vertices[indmax1],vertices[indmax2])
     cross12 = cross12/np.sqrt(np.sum(cross12**2))
     indmax3, odfmax3 = patch_maximum(vertices, odf, cross12, 2*width)
     """
-    return [(indmax1, odfmax1),(indmax2, odfmax2),(indmax3, odfmax3)]
+    return [(indmax1, odfmax1), (indmax2, odfmax2), (indmax3, odfmax3)]

dipy/reconst/interpolate.py

@@ -1,4 +1,5 @@
-"""Interpolators wrap arrays to allow the array to be indexed in continuous coordinates
+"""Interpolators wrap arrays to allow the array to be indexed in
+continuous coordinates
 
 This module uses the trackvis coordinate system, for more information about
 this coordinate system please see dipy.tracking.utils
@@ -10,15 +11,18 @@
 from numpy import array
 from dipy.reconst.recspeed import trilinear_interp
 
+
 class OutsideImage(Exception):
     pass
 
+
 class Interpolator(object):
     """Class to be subclassed by different interpolator types"""
     def __init__(self, data, voxel_size):
         self.data = data
         self.voxel_size = array(voxel_size, dtype=float, copy=True)
 
+
 class NearestNeighborInterpolator(Interpolator):
     """Interpolates data using nearest neighbor interpolation"""
 
@@ -31,6 +35,7 @@ def __getitem__(self, index):
         except IndexError:
             raise OutsideImage
 
+
 class TriLinearInterpolator(Interpolator):
     """Interpolates data using trilinear interpolation
 

dipy/reconst/mapmri.py

@@ -923,7 +923,6 @@ def predict(self, qvals_or_gtab, S0=100.):
         E = S0 * np.dot(M, self._mapmri_coef)
         return E
 
-
     def pdf(self, r_points):
         """ Diffusion propagator on a given set of real points.
         if the array r_points is non writeable, then intermediate

dipy/reconst/odf.py

@@ -5,6 +5,7 @@
 # Classes OdfModel and OdfFit are using API ReconstModel and ReconstFit from
 # .base
 
+
 class OdfModel(ReconstModel):
 
     """An abstract class to be sub-classed by specific odf models
@@ -67,4 +68,3 @@ def minmax_normalize(samples, out=None):
     out -= sample_mins
     out /= (sample_maxes - sample_mins)
     return out
-