diff --git a/bin/fit_tensor b/bin/fit_tensor
index 8ffaa495f8..2ec744d897 100755
--- a/bin/fit_tensor
+++ b/bin/fit_tensor
@@ -11,7 +11,7 @@ from dipy.io.utils import nifti1_symmat
 from dipy.io.bvectxt import read_bvec_file, orientation_to_string
 from nibabel.trackvis import empty_header, write
 
-usage = """fit_tenor [options] dwi_images"""
+usage = """fit_tensor [options] dwi_images"""
 parser = OptionParser(usage)
 parser.add_option("-b","--bvec",help="text file with gradient directions")
 parser.add_option("-r","--root",help="root for files to be saved")
diff --git a/dipy/core/geometry.py b/dipy/core/geometry.py
index 2decc82ce6..fd745f15ef 100644
--- a/dipy/core/geometry.py
+++ b/dipy/core/geometry.py
@@ -24,7 +24,6 @@
 _TUPLE2AXES = dict((v, k) for k, v in _AXES2TUPLE.items())
 
 
-
 def sphere2cart(r, theta, phi):
     ''' Spherical to Cartesian coordinates
 
@@ -130,6 +129,18 @@ def cart2sphere(x, y, z):
     return r, theta, phi
 
 
+def sph2latlon(theta, phi):
+    """Convert spherical coordinates to latitude and longitude.
+
+    Returns
+    -------
+    lat, lon : ndarray
+        Latitude and longitude.
+
+    """
+    return np.rad2deg(theta - np.pi/2), np.rad2deg(phi - np.pi)
+
+
 def normalized_vector(vec, axis=-1):
     ''' Return vector divided by its Euclidean (L2) norm
 
diff --git a/dipy/core/meshes.py b/dipy/core/meshes.py
deleted file mode 100644
index bcb17c81bd..0000000000
--- a/dipy/core/meshes.py
+++ /dev/null
@@ -1,401 +0,0 @@
-''' Mesh analysis '''
-
-import numpy as np
-from scipy import sparse
-
-FLOAT64_EPS = np.finfo(np.float64).eps
-FLOAT_TYPES = np.sctypes['float']
-
-white = 0
-red = 1
-black = 2
-green = 3
-
-def sym_hemisphere(vertices,
-                   hemisphere='z',
-                   equator_thresh=None,
-                   dist_thresh=None):
-    """ Indices for hemisphere from an array of `vertices` on a sphere 
-
-    Selects the vertices from a sphere that lie in one hemisphere.
-    If there are pairs of symmetric points on the equator, we return only
-    the first occurring of each pair.
-
-    Parameters
-    ----------
-    vertices : (N,3) array-like
-       (x, y, z) Point coordinates of N vertices
-    hemisphere : str, optional
-       Which hemisphere to select.  Values of '-x', '-y', '-z' select,
-       respectively negative x, y, and z hemispheres; 'x', 'y', 'z'
-       select the positive x, y, and z hemispheres.  Default is 'z'
-    equator_thresh : None or float, optional
-       Threshold (+-0) to identify points as being on the equator of the
-       sphere.   If None, generate a default based on the data type
-    dist_thresh : None or float, optional
-       For a vertex ``v`` on the equator, if there is a vertex
-       ``v_dash`` in `vertices`, such that the Euclidean distance
-       between ``v * -1`` and ``v_dash`` is <= `dist_thresh`, then ``v``
-       is taken to be in the opposite hemisphere to ``v_dash``, and only
-       ``v``, not ``v_dash``, will appear in the output vertex indices
-       `inds`. None results in a threshold based on the input data type
-       of ``vertices``
-    
-    Returns
-    -------
-    inds : (P,) array
-       Indices into `vertices` giving points in hemisphere
-
-    Notes
-    -----
-    We expect the sphere to be symmetric, and so there may well be
-    points on the sphere equator that are both on the same diameter
-    line.  The routine returns the first of the two points in the
-    original order of `vertices`.
-    """
-    vertices = np.asarray(vertices)
-    assert vertices.shape[1] == 3
-    if len(hemisphere) == 2:
-        sign, hemisphere = hemisphere
-        if sign not in '+-':
-            raise ValueError('Hemisphere sign must be + or -')
-    else:
-        sign = '+'
-    try:
-        coord = 'xyz'.index(hemisphere)
-    except ValueError:
-        raise ValueError('Hemisphere must be (+-) x, y or z')
-    if equator_thresh is None or dist_thresh is None:
-        if not vertices.dtype.type in FLOAT_TYPES:
-            EPS = FLOAT64_EPS
-        else:
-            EPS = np.finfo(vertices.dtype.type).eps
-        if equator_thresh is None:
-            equator_thresh = EPS * 10
-        if dist_thresh is None:
-            dist_thresh = EPS * 20
-    # column with coordinates for selecting the hemisphere
-    sel_col = vertices[:,coord]
-    if sign == '+':
-        inds = sel_col > -equator_thresh
-    else:
-        inds = sel_col < equator_thresh
-    # find equator points
-    eq_inds, = np.where(
-        (sel_col < equator_thresh) & (sel_col > -equator_thresh))
-    # eliminate later points that are symmetric on equator
-    untested_inds = list(eq_inds)
-    out_inds = []
-    for ind in eq_inds:
-        untested_inds.remove(ind)
-        test_vert = vertices[ind,:] * -1
-        test_dists = np.sum(
-            (vertices[untested_inds,:] - test_vert)**2, axis=1)
-        sym_inds, = np.where(test_dists < dist_thresh)
-        for si in sym_inds:
-            out_ind = untested_inds[si]
-            untested_inds.remove(out_ind)
-            out_inds.append(out_ind)
-        if len(untested_inds) == 0:
-            break
-    inds[out_inds] = False
-    return np.nonzero(inds)[0]
-
-
-def vertinds_to_neighbors(vertex_inds, faces):
-    """ Return indices of neighbors of vertices given `faces`
-
-    Parameters
-    ----------
-    vertex_inds : sequence
-       length N.  Indices of vertices
-    faces : (F, 3) array-like
-       Faces given by indices of vertices for each of ``F`` faces
-
-    Returns
-    -------
-    adj : list
-       For each ``N`` vertex indicated by `vertex_inds`, the vertex
-       indices that are neighbors according to the graph given by
-       `faces`.  
-    """
-    full_adj = neighbors(faces)
-    adj = []
-    for i, n in enumerate(full_adj):
-        if i in vertex_inds:
-            adj.append(n)
-    return adj
-
-
-def neighbors(faces):
-    """ Return indices of neighbors for each vertex within `faces`
-
-    Parameters
-    ----------
-    faces : (F, 3) array-like
-       Faces given by indices of vertices for each of ``F`` faces
-
-    Returns
-    -------
-    adj : list
-       For each vertex found within `faces`, the vertex
-       indices that are neighbors according to the graph given by
-       `faces`.  We expand the list with empty lists in between
-       non-empty neighbors.  
-    """
-    faces = np.asarray(faces)
-    adj = {}
-    for face in faces:
-        a, b, c = face
-        if a in adj:
-            adj[a] += [b, c]
-        else:
-            adj[a] = [b, c]
-        if b in adj:
-            adj[b] += [a, c]
-        else:
-            adj[b] = [a, c]
-        if c in adj:
-            adj[c] += [a, b]
-        else:
-            adj[c] = [a, b]
-    N = max(adj.keys())+1
-    out = [[] for i in range(N)]
-    for i in range(N):
-        if i in adj:
-            out[i] = np.sort(np.unique(adj[i]))
-    return out
-
-
-def vertinds_faces(vertex_inds, faces):
-    """ Return faces containing any of `vertex_inds`
-
-    Parameters
-    ----------
-    vertex_inds : sequence
-       length N.  Indices of vertices
-    faces : (F, 3) array-like
-       Faces given by indices of vertices for each of ``F`` faces
-
-    Returns
-    ---------
-    less_faces : (P, 3) array
-       Only retaining rows in `faces` which contain any of `vertex_inds`
-    """
-    in_inds = []
-    vertex_inds = set(vertex_inds)
-    for ind, face in enumerate(faces):
-        if vertex_inds.intersection(face):
-            in_inds.append(ind)
-    return faces[in_inds]
-
-def vertinds_faceinds(vertex_inds, faces):
-    """ Return indices of faces containing any of `vertex_inds`
-
-    Parameters
-    ----------
-    vertex_inds : sequence
-       length N.  Indices of vertices
-    faces : (F, 3) array-like
-       Faces given by indices of vertices for each of ``F`` faces
-
-    Returns
-    ---------
-    in_inds : sequence
-       Indices of `faces` which contain any of `vertex_inds`
-    """
-    in_inds = []
-    vertex_inds = set(vertex_inds)
-    for ind, face in enumerate(faces):
-        if vertex_inds.intersection(face):
-            in_inds.append(ind)
-    return in_inds
-
-def edges(vertex_inds, faces):
-    r""" Return array of starts and ends of edges from list of faces
-        taking regard of direction.
-
-    Parameters
-    ----------
-    vertex_inds : sequence
-       length N.  Indices of vertices
-    faces : (F, 3) array-like
-       Faces given by indices of vertices for each of F faces
-
-    Returns
-    -------
-    edgearray : (E2, 2) array
-       where E2 = 2*E, twice the number of edges. If e= (a,b) is an
-       edge then [a,b] and [b,a] are included in edgearray.
-    """
-
-    edgedic = {}
-    for face in faces:
-        edgedic[(face[0],face[1])]=1
-        edgedic[(face[0],face[2])]=1
-        edgedic[(face[1],face[0])]=1
-        edgedic[(face[1],face[2])]=1
-        edgedic[(face[2],face[0])]=1
-        edgedic[(face[2],face[1])]=1
-    
-    start, end = zip(*edgedic)
-
-    edgearray = np.column_stack(zip(*edgedic))
-
-    return edgearray
- 
-def vertex_adjacencies(vertex_inds, faces):
-    """ Return matrix which shows the adjacent vertices
-        of each vertex
-        
-    Parameters
-    ----------
-    vertex_inds : sequence
-       length N.  Indices of vertices
-    
-    faces : (F, 3) array-like
-       Faces given by indices of vertices for each of F faces
-    
-    Returns
-    -------
-    """
-    edgearray = edges(vertex_inds, faces)
-    V = len(vertex_inds)
-    a = sparse.coo_matrix((np.ones(edgearray.shape[0]),
-                          (edgearray[:,0],edgearray[:,1])),
-                          shape=(V,V))
-    return a
-
-    
-def argmax_from_adj(vals, vertex_inds, adj_inds):
-    """ Indices of local maxima from `vals` given adjacent points
-
-    See ``reconstruction_performance`` for optimized versions of this
-    routine. 
-    
-    Parameters
-    ----------
-    vals : (N,) array-like
-       values at all vertices referred to in either of `vertex_inds` or
-       `adj_inds`'
-    vertex_inds : None or (V,) array-like
-       indices into `vals` giving vertices that may be local maxima.
-       If None, then equivalent to ``np.arange(N)``
-    adj_inds : sequence
-       For every vertex in ``vertex_inds``, the indices (into `vals`) of
-       the neighboring points
-
-    Returns
-    -------
-    inds : (M,) array
-       Indices into `vals` giving local maxima of vals, given topology
-       from `adj_inds`, and restrictions from `vertex_inds`.  Inds are
-       returned sorted by value at that index - i.e. smallest value (at
-       index) first.
-    """
-    vals = np.asarray(vals)
-    if vertex_inds is None:
-        vertex_inds = np.arange(vals.shape[0])
-    else:
-        vertex_inds = np.asarray(vertex_inds)
-    maxes = []
-    for i, adj in enumerate(adj_inds):
-        vert_ind = vertex_inds[i]
-        val = vals[vert_ind]
-        if np.all(val > vals[adj]):
-            maxes.append((val, vert_ind))
-    if len(maxes) == 0:
-        return np.array([])
-    maxes.sort(cmp=lambda x, y: cmp(x[0], y[0]))
-    vals, inds = zip(*maxes)
-    return np.array(inds)
-
-
-def peak_finding_compatible(vertices,
-                            hemisphere='z',
-                            equator_thresh=None,
-                            dist_thresh=None):
-    """ Check that a sphere mesh is compatible with ``peak_finding``
-
-    Parameters
-    ----------
-    vertices : (N,3) array-like
-       (x, y, z) Point coordinates of N vertices
-    hemisphere : str, optional
-       Which hemisphere to select.  Values of '-x', '-y', '-z' select,
-       respectively negative x, y, and z hemispheres; 'x', 'y', 'z'
-       select the positive x, y, and z hemispheres.  Default is 'z'
-    equator_thresh : None or float, optional
-       Threshold (+-0) to identify points as being on the equator of the
-       sphere.   If None, generate a default based on the data type
-    dist_thresh : None or float, optional
-       For a vertex ``v`` on the equator, if there is a vertex
-       ``v_dash`` in `vertices`, such that the Euclidean distance
-       between ``v * -1`` and ``v_dash`` is <= `dist_thresh`, then ``v``
-       is taken to be in the opposite hemisphere to ``v_dash``, and only
-       ``v``, not ``v_dash``, will appear in the output vertex indices
-       `inds`. None results in a threshold based on the input data type
-       of ``vertices``
-    
-    Returns
-    -------
-    compatible : bool
-       True if the sphere mesh is compatible with ``peak_finding``
-    """
-    inds = sym_hemisphere(vertices, hemisphere,
-                          equator_thresh, dist_thresh)
-    N = vertices.shape[0] // 2
-    return np.all(inds == np.arange(N))
-
-def euler_characteristic_check(vertices, faces, chi=2):
-    r'''
-    If $f$ = number of faces, $e$ = number_of_edges and $v$ = number of vertices,
-    the Euler formula says $f-e+v = 2$ for a mesh
-    on a sphere. Here, assuming we have a healthy triangulation every
-    face is a triangle, all 3 of whose edges should belong to exactly
-    two faces. So $2*e = 3*f$. To avoid integer division and consequential
-    integer rounding we test whether $2*f - 3*f + 2*v == 4$ or, more generally,
-    whether $2*v - f == 2*\chi$ where $\chi$ is the Euler characteristic of the mesh.
-
-    - Open chain (track) has $\chi=1$
-    - Closed chain (loop) has $\chi=0$
-    - Disk has $\chi=1$
-    - Sphere has $\chi=2$
-
-    Parameters
-    ----------
-    vertices : (N,3) array-like
-       (x, y, z) Point coordinates of N vertices
-    faces : (M,3) array-like of type int
-       (i1, i2, i3) Integer indices of the vertices of the (triangular) faces
-    chi : int, or None
-       The Euler characteristic of the mesh to be checked
-       
-    Returns
-    -------
-    check : bool
-       True if the mesh has Euler characteristic chi
-    '''
-    
-    v = vertices.shape[0]
-    f = faces.shape[0]
-    if 2*v-f==2*chi:
-        return True
-    else:
-        return False
-
-def adjacent_uncoloured(vertinds, vertex_colour, face_colour, faces):
-
-    adjacent_faces = np.array(vertinds_faceinds(vertinds, faces))
-
-    uncoloured_adjacent_faces = adjacent_faces[np.where(face_colour[adjacent_faces]==white)]
-
-    adjacent_vertices = np.array(list(set(faces[uncoloured_adjacent_faces].ravel())))
-    
-    l = list(adjacent_vertices)
-    
-    w = np.where(vertex_colour[l]==white)
-    
-    uncoloured_adjacent_vertices = adjacent_vertices[w]
-    
-    return(uncoloured_adjacent_vertices, uncoloured_adjacent_faces)
diff --git a/dipy/core/ndindex.py b/dipy/core/ndindex.py
new file mode 100644
index 0000000000..005be9c89c
--- /dev/null
+++ b/dipy/core/ndindex.py
@@ -0,0 +1,37 @@
+import numpy as np
+from numpy.lib.stride_tricks import as_strided
+
+def ndindex(shape):
+    """
+    An N-dimensional iterator object to index arrays.
+
+    Given the shape of an array, an `ndindex` instance iterates over
+    the N-dimensional index of the array. At each iteration a tuple
+    of indices is returned; the last dimension is iterated over first.
+
+    Parameters
+    ----------
+    shape : tuple of ints
+      The dimensions of the array.
+
+    Examples
+    --------
+    >>> from dipy.core.ndindex import ndindex
+    >>> shape = (3, 2, 1)
+    >>> for index in ndindex(shape):
+    ...     print index
+    (0, 0, 0)
+    (0, 1, 0)
+    (1, 0, 0)
+    (1, 1, 0)
+    (2, 0, 0)
+    (2, 1, 0)
+
+    """
+    if len(shape) == 0:
+        yield ()
+    else:
+        x = as_strided(np.zeros(1), shape=shape, strides=np.zeros_like(shape))
+        ndi = np.nditer(x, flags=['multi_index', 'zerosize_ok'], order='C')
+        for e in ndi:
+            yield ndi.multi_index
diff --git a/dipy/core/sphere.py b/dipy/core/sphere.py
index e7b3a0d402..56f52bc613 100644
--- a/dipy/core/sphere.py
+++ b/dipy/core/sphere.py
@@ -39,8 +39,11 @@ def faces_from_sphere_vertices(vertices):
 
     """
     from scipy.spatial import Delaunay
-    return Delaunay(vertices).convex_hull
-
+    faces = Delaunay(vertices).convex_hull
+    if len(vertices) < 2**16:
+        return np.asarray(faces, np.uint16)
+    else:
+        return faces
 
 def unique_edges(faces, return_mapping=False):
     """Extract all unique edges from given triangular faces.
@@ -196,8 +199,6 @@ def z(self):
     @auto_attr
     def faces(self):
         faces = faces_from_sphere_vertices(self.vertices)
-        if len(self.theta) < 2**16:
-            faces = np.asarray(faces, dtype='uint16')
         return faces
 
     @auto_attr
@@ -299,7 +300,8 @@ def __init__(self, x=None, y=None, z=None,
         """Create a HemiSphere from points"""
 
         sphere = Sphere(x=x, y=y, z=z, theta=theta, phi=phi, xyz=xyz)
-        uniq_vertices, mapping = remove_similar_vertices(sphere.vertices, tol)
+        uniq_vertices, mapping = remove_similar_vertices(sphere.vertices, tol,
+                                                         return_mapping=True)
         uniq_vertices *= 1 - 2*(uniq_vertices[:, -1:] < 0)
         if faces is not None:
             faces = np.asarray(faces)
@@ -486,6 +488,47 @@ def interp_rbf(data, sphere_origin, sphere_target,
     return rbfi(sphere_target.x, sphere_target.y, sphere_target.z)
 
 
+def euler_characteristic_check(sphere, chi=2):
+    r"""Checks the euler characteristic of a sphere
+
+    If $f$ = number of faces, $e$ = number_of_edges and $v$ = number of
+    vertices, the Euler formula says $f-e+v = 2$ for a mesh on a sphere. More
+    generally, whether $f -e + v == \chi$ where $\chi$ is the Euler
+    characteristic of the mesh.
+
+    - Open chain (track) has $\chi=1$
+    - Closed chain (loop) has $\chi=0$
+    - Disk has $\chi=1$
+    - Sphere has $\chi=2$
+    - HemiSphere has $\chi=1$
+
+    Parameters
+    ----------
+    sphere : Sphere
+        A Sphere instance with vertices, edges and faces attributes.
+    chi : int, optional
+       The Euler characteristic of the mesh to be checked
+
+    Returns
+    -------
+    check : bool
+       True if the mesh has Euler characteristic $\chi$
+
+    Examples
+    --------
+    >>> euler_characteristic_check(unit_octahedron)
+    True
+    >>> hemisphere = HemiSphere.from_sphere(unit_icosahedron)
+    >>> euler_characteristic_check(hemisphere, chi=1)
+    True
+
+    """
+    v = sphere.vertices.shape[0]
+    e = sphere.edges.shape[0]
+    f = sphere.faces.shape[0]
+    return (f - e + v) == chi
+
+
 octahedron_vertices = np.array(
     [[ 1.0 , 0.0,  0.0],
      [-1.0,  0.0,  0.0],
diff --git a/dipy/core/tests/test_ndindex.py b/dipy/core/tests/test_ndindex.py
new file mode 100644
index 0000000000..43090f024a
--- /dev/null
+++ b/dipy/core/tests/test_ndindex.py
@@ -0,0 +1,14 @@
+from dipy.core.ndindex import ndindex
+
+import numpy as np
+from numpy.testing import assert_array_equal
+
+def test_ndindex():
+    x = list(ndindex((1, 2, 3)))
+    expected = [ix for ix, e in np.ndenumerate(np.zeros((1, 2, 3)))]
+    assert_array_equal(x, expected)
+
+def test_ndindex_0d():
+    x = list(ndindex(np.array(1).shape))
+    expected = [()]
+    assert_array_equal(x, expected)
diff --git a/dipy/data/__init__.py b/dipy/data/__init__.py
index 758757e8e8..c0d3fbc7f3 100644
--- a/dipy/data/__init__.py
+++ b/dipy/data/__init__.py
@@ -6,6 +6,7 @@
 import cPickle
 import gzip
 from dipy.core.gradients import gradient_table
+from dipy.core.sphere import Sphere
 import numpy as np
 
 from ..utils.arrfuncs import as_native_array
@@ -20,6 +21,7 @@
 class DataError(Exception):
     pass
 
+
 def get_sim_voxels(name='fib1'):
     """ provide some simulated voxel data
 
@@ -106,16 +108,14 @@ def get_sphere(name='symmetric362'):
 
     Returns
     -------
-    vertices : ndarray
-        vertices for sphere
-    faces : ndarray
-        faces
+    sphere : a dipy.core.sphere.Sphere class instance
 
     Examples
     --------
     >>> import numpy as np
     >>> from dipy.data import get_sphere
-    >>> verts, faces = get_sphere('symmetric362')
+    >>> sphere = get_sphere('symmetric362')
+    >>> verts, faces = sphere.vertices, sphere.faces
     >>> verts.shape
     (362, 3)
     >>> faces.shape
@@ -131,8 +131,7 @@ def get_sphere(name='symmetric362'):
     res = np.load(fname)
     # Set to native byte order to avoid errors in compiled routines for
     # big-endian platforms, when using these spheres.
-    return (as_native_array(res['vertices']),
-            as_native_array(res['faces']))
+    return Sphere(xyz=as_native_array(res['vertices']))
 
 
 def get_data(name='small_64D'):
@@ -194,8 +193,9 @@ def get_data(name='small_64D'):
     if name=='grad514':
         return pjoin(THIS_DIR,'grad_514.txt')
 
+
 def dsi_voxels():
-    fimg,fbvals,fbvecs = get_data('small_101D')
+    fimg, fbvals, fbvecs = get_data('small_101D')
     bvals = np.loadtxt(fbvals)
     bvecs = np.loadtxt(fbvecs).T
     img = load(fimg)
diff --git a/dipy/reconst/dti.py b/dipy/reconst/dti.py
index 4b1882381c..e532df5485 100644
--- a/dipy/reconst/dti.py
+++ b/dipy/reconst/dti.py
@@ -1,96 +1,91 @@
 #!/usr/bin/python
-""" Classes and functions for fitting tensors """
-# 5/17/2010
-
+import warnings
 import numpy as np
-
 from dipy.reconst.maskedview import MaskedView, _makearray, _filled
 from dipy.reconst.modelarray import ModelArray
 from dipy.data import get_sphere
+from ..core.geometry import vector_norm
+from dipy.core.onetime import auto_attr
 
-class Tensor(ModelArray):
-    """ Fits a diffusion tensor given diffusion-weighted signals and gradient info
 
-    Tensor object that when initialized calculates single self diffusion
-    tensor [1]_ in each voxel using selected fitting algorithm
-    (DEFAULT: weighted least squares [2]_)
-    Requires a given gradient table, b value for each diffusion-weighted
-    gradient vector, and image data given all as arrays.
+class TensorModel(object):
+    """ Diffusion Tensor
+    """
+    def __init__(self, gtab, fit_method="WLS", *args, **kwargs):
+        """ A Diffusion Tensor Model [1]_, [2]_.
 
-    Parameters
-    ----------
-    data : array ([X, Y, Z, ...], g)
-        Diffusion-weighted signals. The dimension corresponding to the
-        diffusion weighting must be the last dimenssion
-    bval : array (g,)
-        Diffusion weighting factor b for each vector in gtab.
-    gtab : array (g, 3)
-        Diffusion gradient table found in DICOM header as a array.
-    mask : array, optional
-        The tensor will only be fit where mask is True. Mask must must
-        broadcast to the shape of data and must have fewer dimensions than data
-    thresh : float, default = None
-        The tensor will not be fit where data[bval == 0] < thresh. If multiple
-        b0 volumes are given, the minimum b0 signal is used.
-    fit_method : funciton or string, default = 'WLS'
-        The method to be used to fit the given data to a tensor. Any function
-        that takes the B matrix and the data and returns eigen values and eigen
-        vectors can be passed as the fit method. Any of the common fit methods
-        can be passed as a string.
-    *args, **kargs :
-        Any other arguments or keywards will be passed to fit_method.
-
-    common fit methods:
-        'WLS' : weighted least squares
-            dti.wls_fit_tensor
-        'LS' : ordinary least squares
-            dti.ols_fit_tensor
-
-    Attributes
-    ----------
-    D : array (..., 3, 3)
-        Self diffusion tensor calculated from cached eigenvalues and 
-        eigenvectors.
-    mask : array
-        True in voxels where a tensor was fit, false if the voxel was skipped
-    B : array (g, 7)
-        Design matrix or B matrix constructed from given gradient table and
-        b-value vector.
-    evals : array (..., 3) 
-        Cached eigenvalues of self diffusion tensor for given index. 
-        (eval1, eval2, eval3)
-    evecs : array (..., 3, 3)
-        Cached associated eigenvectors of self diffusion tensor for given 
-        index. Note: evals[..., j] is associated with evecs[..., :, j]
-
-
-    Methods
-    -------
-    fa : array
-        Calculates fractional anisotropy [2]_.
-    md : array
-        Calculates the mean diffusivity [2]_.
-        Note: [units ADC] ~ [units b value]*10**-1
+        Parameters
+        ----------
+        gtab : GradientTable
+        fit_method : str or callable 
+            str can be one of the following: 
+            'WLS' for weighted least squares
+                dti.wls_fit_tensor
+            'LS' for ordinary least squares
+                dti.ols_fit_tensor
+        
+            callable has to have the signature:
+              fit_method(design_matrix, data, *args, **kwargs)
 
-    See Also
-    --------
-    dipy.io.bvectxt.read_bvec_file, dipy.core.qball.ODF
+        args, kwargs : arguments and key-word arguments passed to the
+           fit_method. See dti.wls_fit_tensor, dti.ols_fit_tensor for details 
 
-    References
-    ----------
-    .. [1] Basser, P.J., Mattiello, J., LeBihan, D., 1994. Estimation of
-       the effective self-diffusion tensor from the NMR spin echo. J Magn
-       Reson B 103, 247-254.
-    .. [2] Basser, P., Pierpaoli, C., 1996. Microstructural and physiological
-       features of tissues elucidated by quantitative diffusion-tensor MRI.
-       Journal of Magnetic Resonance 111, 209-219.
-
-    Examples
-    ----------
-    For a complete example have a look at the main dipy/examples folder
-    """
+        References
+        ----------
+        .. [1] Basser, P.J., Mattiello, J., LeBihan, D., 1994. Estimation of
+           the effective self-diffusion tensor from the NMR spin echo. J Magn
+           Reson B 103, 247-254.
+        .. [2] Basser, P., Pierpaoli, C., 1996. Microstructural and
+           physiological features of tissues elucidated by quantitative
+           diffusion-tensor MRI.  Journal of Magnetic Resonance 111, 209-219.
+
+        """
+        if not callable(fit_method):
+            try:
+                self.fit_method = common_fit_methods[fit_method]
+            except KeyError:
+                raise ValueError('"'+str(fit_method)+'" is not a known fit '
+                                 'method, the fit method should either be a '
+                                 'function or one of the common fit methods')
+        self.bvec = gtab.bvecs
+        self.bval = gtab.bvals
+        self.design_matrix = design_matrix(self.bvec.T, self.bval)
+        self.args = args
+        self.kwargs = kwargs
+
+    def fit(self, data):
+        """
+        Fit method of the DTI model class
+
+        Parameters
+        ----------
+        data : array
+            The measured signal from one voxel.
+
+        """
+        dti_params = self.fit_method(self.design_matrix, data,
+                                     *self.args, **self.kwargs)
+        return TensorFit(self, dti_params)
+
+
+class TensorFit(object):
+    def __init__(self, model, model_params):
+        """
+        Initialize a TensorFit class instance.
+        """
+        self.model_params = model_params
+
+    @property
+    def shape(self):
+        return self.model_params.shape[:-1]
+
+    @property
+    def directions(self):
+        """
+        For tracking - return the primary direction in each voxel
+        """
+        return self.evecs[0,0]
 
-    ### Eigenvalues Property ###
     @property
     def evals(self):
         """
@@ -98,7 +93,6 @@ def evals(self):
         """
         return _filled(self.model_params[..., :3])
 
-    ### Eigenvectors Property ###
     @property
     def evecs(self):
         """
@@ -108,85 +102,22 @@ def evecs(self):
         evecs = _filled(self.model_params[..., 3:])
         return evecs.reshape(self.shape + (3, 3))
 
-    def __init__(self, data, b_values, grad_table, mask=True, thresh=None,
-                 fit_method='WLS', verbose=False, *args, **kargs):
-        """
-        Fits a tensors to diffusion weighted data.
-
-        """
-
-        if not callable(fit_method):
-            try:
-                fit_method = common_fit_methods[fit_method]
-            except KeyError:
-                raise ValueError('"'+str(fit_method)+'" is not a known fit '+
-                                 'method, the fit method should either be a '+
-                                 'function or one of the common fit methods')
-
-        #64 bit design matrix makes for faster pinv
-        B = design_matrix(grad_table.T, b_values)
-        self.B = B
-
-        mask = np.atleast_1d(mask)
-        if thresh is not None:
-            #Define total mask from thresh and mask
-            #mask = mask & (np.min(data[..., b_values == 0], -1) >
-            #thresh)
-            #the assumption that the lowest b_value is always 0 is
-            #incorrect the lowest b_value could also be higher than 0
-            #this is common with grid q-spaces
-            min_b0_sig = np.min(data[..., b_values == b_values.min()], -1)
-            mask = mask & (min_b0_sig > thresh)
-
-        #if mask is all False
-        if not mask.any():
-            raise ValueError('between mask and thresh, there is no data to '+
-            'fit')
-
-        #and the mask is not all True
-        if not mask.all():
-            #leave only data[mask is True]
-            data = data[mask]
-            data = MaskedView(mask, data, fill_value=0)
-
-        #Perform WLS fit on masked data
-        dti_params = fit_method(B, data, *args, **kargs)
-        self.model_params = dti_params
-
-    ### Self Diffusion Tensor Property ###
-    def _getD(self):
+    @property
+    def quadratic_form(self):
         """Calculates the 3x3 diffusion tensor for each voxel"""
-        params, wrap = _makearray(self.model_params)
-        evals = params[..., :3]
-        evecs = params[..., 3:]
-        evals_flat = evals.reshape((-1, 3))
-        evecs_flat = evecs.reshape((-1, 3, 3))
-        D_flat = np.empty(evecs_flat.shape)
-        for ii in xrange(len(D_flat)):
-            Q = evecs_flat[ii]
-            L = evals_flat[ii]
-            D_flat[ii] = np.dot(Q*L, Q.T)
-        D = _filled(wrap(D_flat))
-        D.shape = self.shape + (3, 3)
-        return D
-
-    D = property(_getD, doc = "Self diffusion tensor")
+        evecs = self.evecs
+        evals = self.evals
+        # use einsum to do `evecs * evals * evecs.T` where * is matrix multiply
+        return np.einsum('...ij,...j,...kj->...ik', evecs, evals, evecs)
 
     def lower_triangular(self, b0=None):
-        D = self._getD()
-        return lower_triangular(D, b0)
+        return lower_triangular(self.quadratic_form, b0)
 
-    def fa(self, fill_value=0, nonans=True):
+    @auto_attr
+    def fa(self):
         r"""
         Fractional anisotropy (FA) calculated from cached eigenvalues.
 
-        Parameters
-        ----------
-        fill_value : float
-            value of fa where self.mask == True.
-        nonans : Bool
-            When True, fa is 0 when all eigenvalues are 0, otherwise fa is nan
-
         Returns
         ---------
         fa : array (V, 1)
@@ -208,15 +139,17 @@ def fa(self, fill_value=0, nonans=True):
         ev2 = evals[..., 1]
         ev3 = evals[..., 2]
 
-        if nonans:
-            all_zero = (ev1 == 0) & (ev2 == 0) & (ev3 == 0)
-        else:
-            all_zero = 0.
+        # Make sure not to get nans:
+        all_zero = (ev1 == 0) & (ev2 == 0) & (ev3 == 0)
+
         fa = np.sqrt(0.5 * ((ev1 - ev2)**2 + (ev2 - ev3)**2 + (ev3 - ev1)**2)
                       / (ev1*ev1 + ev2*ev2 + ev3*ev3 + all_zero))
+
         fa = wrap(np.asarray(fa))
-        return _filled(fa, fill_value)
+        # Fill with zeros outside of the mask
+        return _filled(fa, 0)
 
+    @auto_attr
     def md(self):
         r"""
         Mean diffusitivity (MD) calculated from cached eigenvalues.
@@ -232,24 +165,24 @@ def md(self):
 
         .. math::
 
-            ADC = \frac{\lambda_1+\lambda_2+\lambda_3}{3}
+            MD = \frac{\lambda_1+\lambda_2+\lambda_3}{3}
         """
-        #adc/md = (ev1+ev2+ev3)/3
         return self.evals.mean(-1)
 
-        
-    def ind(self):
-        ''' Quantizes eigenvectors with maximum eigenvalues  on an
-        evenly distributed sphere so that the can be used for tractography.
+    def odf(self, sphere):
+        lower = 4 * np.pi * np.sqrt(np.prod(self.evals, -1))
+        projection = np.dot(sphere.vertices, self.evecs)
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore")
+            projection /=  np.sqrt(self.evals)
+            odf = (vector_norm(projection) ** -3) / lower
+        # Zero evals are non-physical, we replace nans with zeros
+        any_zero = (self.evals == 0).any(-1)
+        odf = np.where(any_zero, 0, odf)
+        # Move odf to be on the last dimension
+        odf = np.rollaxis(odf, 0, odf.ndim)
+        return odf
 
-        Returns
-        ---------
-        IN : array, shape(x,y,z) integer indices for the points of the
-        evenly distributed sphere representing tensor  eigenvectors of
-        maximum eigenvalue
-    
-        '''
-        return quantize_evecs(self.evecs,odf_vertices=None)
 
 def wls_fit_tensor(design_matrix, data, min_signal=1):
     r"""
@@ -332,6 +265,7 @@ def wls_fit_tensor(design_matrix, data, min_signal=1):
     dti_params = wrap(dti_params)
     return dti_params
 
+
 def _wls_iter(ols_fit, design_matrix, sig, min_signal=1):
     '''
     Function used by wls_fit_tensor for later optimization.
@@ -343,6 +277,7 @@ def _wls_iter(ols_fit, design_matrix, sig, min_signal=1):
     tensor = from_lower_triangular(D)
     return decompose_tensor(tensor)
 
+
 def _ols_iter(inv_design, sig, min_signal=1):
     '''
     Function used by ols_fit_tensor for later optimization.
@@ -425,6 +360,7 @@ def ols_fit_tensor(design_matrix, data, min_signal=1):
     dti_params = wrap(dti_params)
     return dti_params
 
+
 def _ols_fit_matrix(design_matrix):
     """
     Helper function to calculate the ordinary least squares (OLS)
@@ -444,9 +380,12 @@ def _ols_fit_matrix(design_matrix):
     U,S,V = np.linalg.svd(design_matrix, False)
     return np.dot(U, U.T)
 
+
 _lt_indices = np.array([[0, 1, 3],
                         [1, 2, 4],
                         [3, 4, 5]])
+
+
 def from_lower_triangular(D):
     """
     Returns a tensor given the six unique tensor elements
@@ -467,8 +406,11 @@ def from_lower_triangular(D):
     """
     return D[..., _lt_indices]
 
+
 _lt_rows = np.array([0, 1, 1, 2, 2, 2])
 _lt_cols = np.array([0, 0, 1, 0, 1, 2])
+
+
 def lower_triangular(tensor, b0=None):
     """
     Returns the six lower triangular values of the tensor and a dummy variable
@@ -497,7 +439,8 @@ def lower_triangular(tensor, b0=None):
         D[..., :6] = tensor[..., _lt_rows, _lt_cols]
         return D
 
-def tensor_eig_from_lo_tri(B, data):
+
+def tensor_eig_from_lo_tri(data):
     """Calculates parameters for creating a Tensor instance
 
     Calculates tensor parameters from the six unique tensor elements. This
@@ -506,8 +449,6 @@ def tensor_eig_from_lo_tri(B, data):
 
     Parameters:
     -----------
-    B :
-        not currently used
     data : array_like (..., 6)
         diffusion tensors elements stored in lower triangular order
 
@@ -531,6 +472,7 @@ def tensor_eig_from_lo_tri(B, data):
     dti_params = wrap(dti_params)
     return dti_params
 
+
 def decompose_tensor(tensor):
     """
     Returns eigenvalues and eigenvectors given a diffusion tensor
@@ -574,6 +516,7 @@ def decompose_tensor(tensor):
 
     return eigenvals, eigenvecs
 
+
 def design_matrix(gtab, bval, dtype=None):
     """
     Constructs design matrix for DTI weighted least squares or least squares 
@@ -608,6 +551,7 @@ def design_matrix(gtab, bval, dtype=None):
     B[:, 6] = np.ones(bval.size)
     return -B
 
+
 def quantize_evecs(evecs, odf_vertices=None):
     ''' Find the closest orientation of an evenly distributed sphere
 
@@ -624,70 +568,224 @@ def quantize_evecs(evecs, odf_vertices=None):
     '''
     max_evecs=evecs[...,:,0]
     if odf_vertices==None:
-        odf_vertices, _ = get_sphere('symmetric362')
+        odf_vertices = get_sphere('symmetric362').vertices
     tup=max_evecs.shape[:-1]
     mec=max_evecs.reshape(np.prod(np.array(tup)),3)
     IN=np.array([np.argmin(np.dot(odf_vertices,m)) for m in mec])
     IN=IN.reshape(tup)
     return IN
 
-class TensorStepper(object):
-    """Used for tracking diffusion tensors, has a next_step method"""
-    def _get_angel_limit(self):
-        return np.arccos(self.dot_limit)*180/pi
-    def _set_angel_limit(self, angle):
-        if angle >= 0 and angle <= 90:
-            self.dot_limit = cos(angle*pi/180)
-        else:
-            raise ValueError("angle should be between 0 and 180")
-    angel_limit = property(_get_angel_limit, _set_angel_limit)
-
-    def _get_fa_limit(self):
-        return self._fa_limit
-    def _set_fa_limit(self, arg):
-        self._fa_limit = arg
-        mask = self.fa_vol > arg
-        self._interp_inst = _interpolator(self.evec1_vol, self.voxe_size, mask)
-    fa_limit = property(_get_fa_limit, _set_fa_limit)
-
-    def __init__(self, fa_vol, evec1_vol, voxel_size, interpolator,
-                 fa_limit=None, angle_limit=None):
-        self.voxel_size = voxel_size
-        self.angle_limit = angle_limit
-        if fa_vol.shape != evec1_vol.shape[:-1]:
-            msg = "the fa and eigen vector volumes are not the same shape"
-            raise ValueError(msg)
-        if evec1_vol.shape[-1] != 3:
-            msg = "eigen vector volume should have vecetors of length 3 " + \
-                  "along the last dimmension"
-            raise ValueError(msg)
-        self.evec1_vol = evec1_vol
-        self.fa_vol = fa_vol
-        self._interpolator = interpolator
-        #self._interp_inst is created when fa_limit is set
-        self.fa_limit = fa_limit
-
-    def next_step(location, prev_step):
-        """Returns the nearest neighbor tensor for location"""
-        step = self._interp_inst[location]
-        angle_dot = dot(step, prev_step)
-        if np.abs(angle_dot) < self.dot_limit:
-            raise StopIteration
-        if angle_dot > 0:
-            return step
-        else:
-            return -step
-
-def stepper_from_tensor(tensor, *args, **kargs):
-    """stepper_from_tensor(tensor, fa_vol, evec1_vol, voxel_size, interpolator)
-    """
-    fa_vol = tensor.fa()
-    evec1_vol = tensor.evec[..., 0]
-    stepper = TensorStepper(fa_vol, evec1_vol, *args, **kargs)
-    return stepper
-
 common_fit_methods = {'WLS': wls_fit_tensor,
                       'LS': ols_fit_tensor,
-                      'from lower triangular': tensor_eig_from_lo_tri,
+                      'OLS': ols_fit_tensor,
                      }
 
+
+
+# For backwards compatibility:
+class Tensor(TensorFit, ModelArray):
+    """
+    For backwards compatibility, we continue to support this form of the Tensor
+    fitting.
+
+    """
+    def __init__(self, data, b_values, b_vectors, mask=True, thresh=None,
+                 fit_method='WLS', verbose=False, *args, **kargs):
+        """ Fits tensors to diffusion weighted data.
+
+        Fits a diffusion tensor given diffusion-weighted signals and gradient
+        info. Tensor object that when initialized calculates single self
+        diffusion tensor [1]_ in each voxel using selected fitting algorithm
+        (DEFAULT: weighted least squares [2]_) Requires a given b-vector table,
+        b value for each diffusion-weighted gradient vector, and image data
+        given all as arrays.
+
+        Parameters
+        ----------
+        data : array ([X, Y, Z, ...], g)
+            Diffusion-weighted signals. The dimension corresponding to the
+            diffusion weighting must be the last dimension
+
+        bval : array (g,)
+            Diffusion weighting factor b for each vector in gtab.
+
+        bvec : array (g, 3)
+            Diffusion gradient table found in DICOM header as a array.
+
+        mask : array, optional
+            The tensor will only be fit where mask is True. Mask must must
+            broadcast to the shape of data and must have fewer dimensions than
+            data
+
+        thresh : float, default = None
+            The tensor will not be fit where data[bval == 0] < thresh. If
+            multiple b0 volumes are given, the minimum b0 signal is used.
+
+        fit_method : funciton or string, default = 'WLS'
+            The method to be used to fit the given data to a tensor. Any
+            function that takes the B matrix and the data and returns eigen
+            values and eigen vectors can be passed as the fit method. Any of
+            the common fit methods can be passed as a string.
+
+        *args, **kargs :
+            Any other arguments or keywards will be passed to fit_method.
+
+        common fit methods:
+
+            'WLS' : weighted least squares
+
+                dti.wls_fit_tensor
+
+            'LS' : ordinary least squares
+
+                dti.ols_fit_tensor
+
+        Attributes
+        ----------
+        D : array (..., 3, 3)
+            Self diffusion tensor calculated from cached eigenvalues and
+            eigenvectors.
+        mask : array
+            True in voxels where a tensor was fit, false if the voxel was skipped
+        B : array (g, 7)
+            Design matrix or B matrix constructed from given gradient table and
+            b-value vector.
+        evals : array (..., 3)
+            Cached eigenvalues of self diffusion tensor for given index.
+            (eval1, eval2, eval3)
+        evecs : array (..., 3, 3)
+            Cached associated eigenvectors of self diffusion tensor for given
+            index. Note: evals[..., j] is associated with evecs[..., :, j]
+
+        Methods
+        -------
+        fa : array
+            Calculates fractional anisotropy [2]_.
+        md : array
+            Calculates the mean diffusivity [2]_.
+            Note: [units ADC] ~ [units b value]*10**-1
+
+        Examples
+        ----------
+        For a complete example have a look at the main dipy/examples folder
+
+        """
+        warnings.warn("This implementation of DTI will be deprecated in a future release, consider using TensorModel", DeprecationWarning)
+        if not callable(fit_method):
+            try:
+                fit_method = common_fit_methods[fit_method]
+            except KeyError:
+                raise ValueError('"'+str(fit_method)+'" is not a known fit '+
+                                 'method, the fit method should either be a '+
+                                 'function or one of the common fit methods')
+
+        #64 bit design matrix makes for faster pinv
+        B = design_matrix(b_vectors.T, b_values)
+        self.B = B
+
+        mask = np.atleast_1d(mask)
+        if thresh is not None:
+            #Define total mask from thresh and mask
+            #mask = mask & (np.min(data[..., b_values == 0], -1) >
+            #thresh)
+            #the assumption that the lowest b_value is always 0 is
+            #incorrect the lowest b_value could also be higher than 0
+            #this is common with grid q-spaces
+            min_b0_sig = np.min(data[..., b_values == b_values.min()], -1)
+            mask = mask & (min_b0_sig > thresh)
+
+        #if mask is all False
+        if not mask.any():
+            raise ValueError('between mask and thresh, there is no data to '+
+            'fit')
+
+        #and the mask is not all True
+        if not mask.all():
+            #leave only data[mask is True]
+            data = data[mask]
+            data = MaskedView(mask, data, fill_value=0)
+
+        #Perform WLS fit on masked data
+        dti_params = fit_method(B, data, *args, **kargs)
+        self.model_params = dti_params
+
+    # For backwards compatibility:
+    D = TensorFit.quadratic_form
+
+    def ind(self):
+        """
+        Quantizes eigenvectors with maximum eigenvalues  on an
+        evenly distributed sphere so that the can be used for tractography.
+
+        Returns
+        ---------
+        IN : array, shape(x,y,z) integer indices for the points of the
+        evenly distributed sphere representing tensor  eigenvectors of
+        maximum eigenvalue
+
+        """
+        return quantize_evecs(self.evecs, odf_vertices=None)
+
+    def fa(self, fill_value=0, nonans=True):
+        r"""
+        Fractional anisotropy (FA) calculated from cached eigenvalues.
+
+        Parameters
+        ----------
+        fill_value : float
+            value of fa where self.mask == True.
+
+        nonans : Bool
+            When True, fa is 0 when all eigenvalues are 0, otherwise fa is nan
+
+        Returns
+        ---------
+        fa : array (V, 1)
+            Calculated FA. Note: range is 0 <= FA <= 1.
+
+        Notes
+        --------
+        FA is calculated with the following equation:
+
+        .. math::
+
+            FA = \sqrt{\frac{1}{2}\frac{(\lambda_1-\lambda_2)^2+(\lambda_1-
+                        \lambda_3)^2+(\lambda_2-lambda_3)^2}{\lambda_1^2+
+                        \lambda_2^2+\lambda_3^2} }
+
+        """
+        evals, wrap = _makearray(self.model_params[..., :3])
+        ev1 = evals[..., 0]
+        ev2 = evals[..., 1]
+        ev3 = evals[..., 2]
+
+        if nonans:
+            all_zero = (ev1 == 0) & (ev2 == 0) & (ev3 == 0)
+        else:
+            all_zero = 0.
+
+        fa = np.sqrt(0.5 * ((ev1 - ev2)**2 + (ev2 - ev3)**2 + (ev3 - ev1)**2)
+                      / (ev1*ev1 + ev2*ev2 + ev3*ev3 + all_zero))
+
+        fa = wrap(np.asarray(fa))
+        return _filled(fa, fill_value)
+
+
+    def md(self):
+        r"""
+        Mean diffusitivity (MD) calculated from cached eigenvalues.
+
+        Returns
+        ---------
+        md : array (V, 1)
+            Calculated MD.
+
+        Notes
+        --------
+        MD is calculated with the following equation:
+
+        .. math::
+
+            MD = \frac{\lambda_1+\lambda_2+\lambda_3}{3}
+        """
+        return self.evals.mean(-1)
diff --git a/dipy/reconst/multi_voxel.py b/dipy/reconst/multi_voxel.py
new file mode 100644
index 0000000000..cebd3ea7f2
--- /dev/null
+++ b/dipy/reconst/multi_voxel.py
@@ -0,0 +1,98 @@
+"""Tools to easily make multi voxel models"""
+import numpy as np
+from numpy.lib.stride_tricks import as_strided
+from numpy import ndindex
+
+
+def multi_voxel_model(SingleVoxelModel):
+    """Class decorator to turn a single voxel model into a multi voxel model
+
+    See Also
+    --------
+    dipy/docs/examples/multiVoxelModel.py
+
+    """
+    class MultiVoxelModel(SingleVoxelModel):
+        """A subclass of SingleVoxelModel that fits many voxels"""
+
+        def fit(self, data, mask=None):
+            """Fit model for every voxel in data"""
+            # If only one voxel just return a normal fit
+            if data.ndim == 1:
+                return SingleVoxelModel.fit(self, data)
+
+            # Make a mask if mask is None
+            if mask is None:
+                shape = data.shape[:-1]
+                strides = (0,) * len(shape)
+                mask = as_strided(np.array(True), shape=shape, strides=strides)
+            # Check the shape of the mask if mask is not None
+            elif mask.shape != data.shape[:-1]:
+                raise ValueError("mask and data shape do not match")
+
+            # Fit data where mask is True
+            fit_array = np.empty(data.shape[:-1], dtype=object)
+            for ijk in ndindex(*data.shape[:-1]):
+                if mask[ijk]:
+                    fit_array[ijk] = SingleVoxelModel.fit(self, data[ijk])
+            return MultiVoxelFit(self, fit_array, mask)
+
+    return MultiVoxelModel
+
+
+class MultiVoxelFit(object):
+    """Holds an array of fits and allows access to their attributes and
+    methods"""
+    def __init__(self, model, fit_array, mask):
+        self.model = model
+        self.fit_array = fit_array
+        self.mask = mask
+
+    @property
+    def shape(self):
+        return self.fit_array.shape
+
+    def __getattribute__(self, attr):
+        try:
+            return object.__getattribute__(self, attr)
+        except AttributeError:
+            result = CallableArray(self.fit_array.shape, dtype=object)
+            for ijk in ndindex(*result.shape):
+                if self.mask[ijk]:
+                    result[ijk] = getattr(self.fit_array[ijk], attr)
+            return _squash(result, self.mask)
+
+    # I leave this in hesitantly, I'm not sure this will be be easy to support
+    # for all models
+    def __getitem__(self, index):
+        item = self.fit_array[index]
+        if isinstance(item, np.ndarray):
+            return MultiVoxelFit(self.model, item, self.mask[index])
+        else:
+            return item
+
+
+class CallableArray(np.ndarray):
+    """An array which can be called like a function"""
+    def __call__(self, *args, **kwargs):
+        result = np.empty(self.shape, dtype=object)
+        for ijk in ndindex(*self.shape):
+            item = self[ijk]
+            if item is not None:
+                result[ijk] = item(*args, **kwargs)
+        return _squash(result)
+
+
+def _squash(arr, mask=None):
+    """Makes a prettier array"""
+    if mask is None:
+        mask = arr != np.array(None)
+    not_none = arr[mask]
+    not_none = not_none.tolist()
+    tmp = np.array(not_none)
+    if tmp.dtype == object:
+        return arr
+    shape = arr.shape + tmp.shape[1:]
+    result = np.zeros(shape, tmp.dtype)
+    result[mask] = tmp
+    return result
diff --git a/dipy/reconst/odf.py b/dipy/reconst/odf.py
index b7cf0a1cc2..1699ebd204 100644
--- a/dipy/reconst/odf.py
+++ b/dipy/reconst/odf.py
@@ -90,13 +90,14 @@ def peak_directions(odf, sphere, relative_peak_threshold,
     ----------
     odf : 1d ndarray
         The odf function evaluated on the vertices of `sphere`
-    sphere :
-        The sphere on which odf was evaluated
+    sphere : Sphere
+        The Sphere providing discrete directions for evaluation.
     relative_peak_threshold : float
-        A relative threshold for excluding small peaks
-    min_separation_angle : float
-        An angle threshold in degrees. Peaks too close to a larger peak are
-        excluded.
+        Only return peaks greater than ``relative_peak_threshold * m`` where m
+        is the largest peak.
+    min_separation_angle : float in [0, 90] The minimum distance between
+        directions. If two peaks are too close only the larger of the two is
+        returned.
 
     Returns
     -------
@@ -118,16 +119,40 @@ def peak_directions(odf, sphere, relative_peak_threshold,
         indices = indices[:first_too_small]
 
     directions = sphere.vertices[indices]
-    directions, mappiing = remove_similar_vertices(directions,
-                                                   min_separation_angle)
+    directions = remove_similar_vertices(directions, min_separation_angle)
     return directions
 
 class PeaksAndMetrics(object):
     pass
 
-def peaks_from_model(model, data, mask=None, return_odf=False, gfa_thr=0.02, 
-                     normalize_peaks=False):
-    """Fits the model to data and computes peaks and metrics"""
+def peaks_from_model(model, data, sphere, relative_peak_threshold,
+                     min_separation_angle, mask=None, return_odf=False,
+                     gfa_thr=0.02, normalize_peaks=False):
+    """Fits the model to data and computes peaks and metrics
+
+    Parameters
+    ----------
+    model : a model instance
+        `model` will be used to fit the data.
+    sphere : Sphere
+        The Sphere providing discrete directions for evaluation.
+    relative_peak_threshold : float
+        Only return peaks greater than ``relative_peak_threshold * m`` where m
+        is the largest peak.
+    min_separation_angle : float in [0, 90] The minimum distance between
+        directions. If two peaks are too close only the larger of the two is
+        returned.
+    mask : array, optional
+        If `mask` is provided, voxels that are False in `mask` are skipped and
+        no peaks are returned.
+    return_odf : bool
+        If True, the odfs are returned.
+    gfa_thr : float
+        Voxels with gfa less than `gfa_thr` are skipped, no peaks are returned.
+    normalize_peaks : bool
+        If true, all peak values are calculated relative to `max(odf)`.
+
+    """
 
     data_flat = data.reshape((-1, data.shape[-1]))
     size = len(data_flat)
@@ -146,13 +171,13 @@ def peaks_from_model(model, data, mask=None, return_odf=False, gfa_thr=0.02,
     peak_indices.fill(-1)
 
     if return_odf:
-        odf_array = np.zeros((size, len(model.sphere.vertices)))
+        odf_array = np.zeros((size, len(sphere.vertices)))
 
     global_max = -np.inf
     for i, sig in enumerate(data_flat):
         if not mask[i]:
             continue
-        odf = model.fit(sig).odf()
+        odf = model.fit(sig).odf(sphere)
         if return_odf:
             odf_array[i] = odf
 
@@ -160,14 +185,22 @@ def peaks_from_model(model, data, mask=None, return_odf=False, gfa_thr=0.02,
         if gfa_array[i] < gfa_thr:
             global_max = max(global_max, odf.max())
             continue
-        pk, ind = local_maxima(odf, model.sphere.edges)
-        """
-        # Will update this later when filter_peaks is nailed down
-        pk, ind = _filter_peaks(pk, ind,
-                                model._distance_matrix,
-                                model.relative_peak_threshold,
-                                model._cos_distance_threshold)
-        """
+        pk, ind = local_maxima(odf, sphere.edges)
+
+        # Remove small peaks.
+        gt_threshold = pk >= (relative_peak_threshold * pk[0])
+        pk = pk[gt_threshold]
+        ind = ind[gt_threshold]
+
+        # Keep peaks which are unique, which means remove peaks that are too
+        # close to a larger peak.
+        _, where_uniq = remove_similar_vertices(sphere.vertices[ind],
+                                                min_separation_angle,
+                                                return_index=True)
+        pk = pk[where_uniq]
+        ind = ind[where_uniq]
+
+        # Calculate peak metrics
         global_max = max(global_max, pk[0])
         n = min(npeaks, len(pk))
         qa_array[i, :n] = pk[:n] - odf.min()
diff --git a/dipy/reconst/recspeed.pyx b/dipy/reconst/recspeed.pyx
index d0c9dd0ced..bfe1bfd949 100644
--- a/dipy/reconst/recspeed.pyx
+++ b/dipy/reconst/recspeed.pyx
@@ -75,45 +75,12 @@ cdef double wght(int i, double r) nogil:
         return 1.-r
 
 
-@cython.wraparound(False)
-def _filter_peaks(cnp.ndarray[cnp.float_t, ndim=1, mode='c'] odf_value,
-                  cnp.ndarray[cnp.int_t, ndim=1, mode='c'] odf_ind,
-                  cnp.ndarray[cnp.float_t, ndim=2, mode='c'] sep_matrix,
-                  float relative_threshold, float isolation):
-    """Filters peaks based on odf_value and angular distance
-
-    Assumes that odf_value is sorted in descending order. Looks up odf_ind in
-    sep_matrix to determine the angular separation between two points. Returns
-    a subset of the peaks that pass the relative_threshold and isolation
-    criterion.
-    """
-    cdef:
-        int i, j, pass_all
-        int count = 1
-        float threshold = relative_threshold * odf_value[0]
-        cnp.ndarray[cnp.int_t, ndim=1, mode='c'] find = odf_ind.copy()
-        cnp.ndarray[cnp.float_t, ndim=1, mode='c'] fvalue = odf_value.copy()
-
-    for i in range(odf_value.shape[0]):
-        if odf_value[i] < threshold:
-            break
-        pass_all = 1
-        for j in range(count):
-            if sep_matrix[odf_ind[i], find[j]] >= isolation:
-                pass_all = 0
-                break
-        if pass_all:
-            find[count] = odf_ind[i]
-            fvalue[count] = odf_value[i]
-            count += 1
-
-    return fvalue[:count].copy(), find[:count].copy()
-
-
 @cython.boundscheck(False)
 @cython.wraparound(False)
 def remove_similar_vertices(cnp.ndarray[cnp.float_t, ndim=2, mode='strided'] vertices,
-                            double theta):
+                            double theta,
+                            bint return_mapping=False,
+                            bint return_index=False):
     """remove_similar_vertices(vertices, theta)
 
     Returns vertices that are separated by at least theta degrees from all
@@ -135,12 +102,14 @@ def remove_similar_vertices(cnp.ndarray[cnp.float_t, ndim=2, mode='strided'] ver
     mapping : (N,) ndarray
         Indices into unique_vertices. For each vertex in `vertices` the index
         of a vertex in `unique_vertices` that is less than theta degrees away.
+
     """
     if vertices.shape[1] != 3:
         raise ValueError()
     cdef:
         cnp.ndarray[cnp.float_t, ndim=2, mode='c'] unique_vertices
         cnp.ndarray[cnp.uint16_t, ndim=1, mode='c'] mapping
+        cnp.ndarray[cnp.uint16_t, ndim=1, mode='c'] index
         char pass_all
         size_t i, j
         size_t count = 0
@@ -150,7 +119,14 @@ def remove_similar_vertices(cnp.ndarray[cnp.float_t, ndim=2, mode='strided'] ver
     if n > 2**16:
         raise ValueError("too many vertices")
     unique_vertices = np.empty((n, 3), dtype=np.float)
-    mapping = np.empty(n, dtype=np.uint16)
+    if return_mapping:
+        mapping = np.empty(n, dtype=np.uint16)
+    else:
+        mapping = None
+    if return_index:
+        index = np.empty(n, dtype=np.uint16)
+    else:
+        index = None
 
     for i in range(n):
         pass_all = 1
@@ -163,16 +139,28 @@ def remove_similar_vertices(cnp.ndarray[cnp.float_t, ndim=2, mode='strided'] ver
                        c * unique_vertices[j, 2])
             if sim > cos_similarity:
                 pass_all = 0
-                mapping[i] = j
+                if return_mapping:
+                    mapping[i] = j
                 break
         if pass_all:
             unique_vertices[count, 0] = a
             unique_vertices[count, 1] = b
             unique_vertices[count, 2] = c
-            mapping[i] = count
+            if return_mapping:
+                mapping[i] = count
+            if return_index:
+                index[count] = i
             count += 1
 
-    return unique_vertices[:count].copy(), mapping
+    if return_mapping and return_index:
+        return unique_vertices[:count].copy(), mapping, index[:count].copy()
+    elif return_mapping:
+        return unique_vertices[:count].copy(), mapping
+    elif return_index:
+        return unique_vertices[:count].copy(), index[:count].copy()
+    else:
+        return unique_vertices[:count].copy()
+
 
 #@cython.boundscheck(False)
 @cython.wraparound(False)
@@ -242,103 +230,6 @@ def local_maxima(cnp.ndarray[cnp.float64_t, ndim=1, mode='c'] codf,
     order = peakvalues.argsort()[::-1]
     return peakvalues[order], peakidx[order]
 
-@cython.boundscheck(False)
-@cython.wraparound(False)
-def peak_finding(odf, odf_faces):
-    ''' Hemisphere local maxima from sphere values and faces
-
-    Return local maximum values and indices. Local maxima (peaks) are
-    given in descending order.
-
-    The sphere mesh, as defined by the vertex coordinates ``vertices``
-    and the face indices ``odf_faces``, has to conform to the check in
-    ``dipy.core.meshes.peak_finding_compatible``.  If it does not, then
-    the results from peak finding routine will be unpredictable.
-
-    Parameters
-    ------------
-    odf : (N,) array of dtype np.float64
-       function values on the sphere, where N is the number of vertices
-       on the sphere
-    odf_faces : (M,3) array of dtype np.uint16
-       faces of the triangulation on the sphere, where M is the number
-       of faces on the sphere
-
-    Returns
-    ---------
-    peaks : (L,) array, dtype np.float64
-       peak values, shape (L,) where L can vary and is the number of
-       local moximae (peaks).  Values are sorted, largest first
-    inds : (L,) array, dtype np.uint16
-       indices of the peak values on the `odf` array corresponding to
-       the maxima in `peaks`
-    
-    Notes
-    -----
-    In summary this function does the following:
-
-    Where the smallest odf values in the vertices of a face put
-    zeros on them. By doing that for the vertices of all faces at the
-    end you only have the peak points with nonzero values.
-
-    For precalculated odf_faces look under
-    dipy/data/evenly*.npz to use them try numpy.load()['faces']
-    
-    Examples
-    ----------
-    This is called from GeneralizedQSampling or QBall and other models with orientation
-    distribution functions.
-
-    See also
-    -----------
-    dipy.core.meshes
-    '''
-    cdef:
-        cnp.ndarray[cnp.uint16_t, ndim=2] cfaces = np.ascontiguousarray(odf_faces)
-        cnp.ndarray[cnp.float64_t, ndim=1] codf = np.ascontiguousarray(odf)
-        cnp.ndarray[cnp.float64_t, ndim=1] cpeak = odf.copy()
-        int i=0
-        int test=0
-        int lenfaces = len(cfaces)
-        double odf0,odf1,odf2
-        int find0,find1,find2
-    
-    for i in range(lenfaces):
-
-        find0 = cfaces[i,0]
-        find1 = cfaces[i,1]
-        find2 = cfaces[i,2]        
-        
-        odf0=codf[find0]
-        odf1=codf[find1]
-        odf2=codf[find2]       
-
-        if odf0 >= odf1 and odf0 >= odf2:
-            cpeak[find1] = 0
-            cpeak[find2] = 0
-            continue
-
-        if odf1 >= odf0 and odf1 >= odf2:
-            cpeak[find0] = 0
-            cpeak[find2] = 0
-            continue
-
-        if odf2 >= odf0 and odf2 >= odf1:
-            cpeak[find0] = 0
-            cpeak[find1] = 0
-            continue
-
-    peak=np.array(cpeak)
-    peak=peak[0:len(peak)/2]
-
-    #find local maxima and give fiber orientation (inds) and magnitude
-    #peaks in a descending order
-
-    inds=np.where(peak>0)[0]
-    pinds=np.argsort(peak[inds])
-    peaks=peak[inds[pinds]][::-1]
-
-    return peaks, inds[pinds][::-1]
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
diff --git a/dipy/reconst/tests/test_dsi.py b/dipy/reconst/tests/test_dsi.py
index d4559e395e..b84de8cd6a 100644
--- a/dipy/reconst/tests/test_dsi.py
+++ b/dipy/reconst/tests/test_dsi.py
@@ -5,8 +5,8 @@
 from dipy.reconst.odf import gfa
 from dipy.sims.voxel import SticksAndBall
 from dipy.core.sphere import Sphere
-from dipy.utils.spheremakers import sphere_vf_from
 from dipy.core.gradients import gradient_table
+from dipy.data import get_sphere
 from numpy.testing import assert_equal
 from dipy.core.subdivide_octahedron import create_unit_sphere
 from dipy.core.sphere_stats import angular_similarity
@@ -14,8 +14,7 @@
 
 def test_dsi():
     #load symmetric 724 sphere
-    vertices, faces = sphere_vf_from('symmetric724')
-    sphere = Sphere(xyz=vertices)
+    sphere = get_sphere('symmetric724')
     #load icosahedron sphere
     sphere2 = create_unit_sphere(5)
     btable = np.loadtxt(get_data('dsi515btable'))    
diff --git a/dipy/reconst/tests/test_dti.py b/dipy/reconst/tests/test_dti.py
index 07233f72b4..f5f39e3a71 100644
--- a/dipy/reconst/tests/test_dti.py
+++ b/dipy/reconst/tests/test_dti.py
@@ -3,19 +3,78 @@
 """
 
 import numpy as np
-from nose.tools import assert_true, assert_false, \
-     assert_equal, assert_almost_equal, assert_raises
+from nose.tools import (assert_true, assert_equal,
+                        assert_almost_equal, assert_raises)
 from numpy.testing import assert_array_equal, assert_array_almost_equal
-from dipy.testing import parametric
-import os
-
 import dipy.reconst.dti as dti
 from dipy.reconst.dti import lower_triangular, from_lower_triangular
 from dipy.reconst.maskedview import MaskedView
-import nibabel as nib
 from dipy.io.bvectxt import read_bvec_file
-from dipy.data import get_data
+from dipy.data import get_data, dsi_voxels
+from dipy.core.subdivide_octahedron import create_unit_sphere
+from dipy.reconst.odf import gfa
+import dipy.core.gradients as grad
+
+def test_TensorModel():
+    data, gtab = dsi_voxels()
+    dm = dti.TensorModel(gtab, 'LS')
+    dtifit = dm.fit(data[0, 0, 0])
+    assert_equal(dtifit.fa < 0.5, True)
+    dm = dti.TensorModel(gtab, 'WLS')
+    dtifit = dm.fit(data[0, 0, 0])
+    assert_equal(dtifit.fa < 0.5, True)
+    sphere = create_unit_sphere(4)
+    assert_equal(len(dtifit.odf(sphere)), len(sphere.vertices))
+    assert_almost_equal(dtifit.fa, gfa(dtifit.odf(sphere)), 1)
+
+    # Check that the multivoxel case works: 
+    dtifit = dm.fit(data)
+    assert_equal(dtifit.fa.shape, data.shape[:3])
+                 
+    # Make some synthetic data
+    b0 = 1000.
+    bvecs, bvals = read_bvec_file(get_data('55dir_grad.bvec'))
+    gtab = grad.gradient_table_from_bvals_bvecs(bvals, bvecs.T)
+    # The first b value is 0., so we take the second one:
+    B = bvals[1]
+    #Scale the eigenvalues and tensor by the B value so the units match
+    D = np.array([1., 1., 1., 0., 0., 1., -np.log(b0) * B]) / B
+    evals = np.array([2., 1., 0.]) / B
+    md = evals.mean()
+    tensor = from_lower_triangular(D)
+    evecs = np.linalg.eigh(tensor)[1]
+    #Design Matrix
+    X = dti.design_matrix(bvecs, bvals)
+    #Signals
+    Y = np.exp(np.dot(X,D))
+    assert_almost_equal(Y[0], b0)
+    Y.shape = (-1,) + Y.shape
 
+    # Test fitting with different methods: #XXX Add NNLS methods! 
+    for fit_method in ['OLS', 'WLS']:
+        tensor_model = dti.TensorModel(gtab,
+                                       fit_method=fit_method)
+        
+        tensor_fit = tensor_model.fit(Y)
+        assert_equal(tensor_fit.shape, Y.shape[:-1])
+        assert_array_almost_equal(tensor_fit.evals[0], evals)
+
+        assert_array_almost_equal(tensor_fit.quadratic_form[0], tensor,
+                                  err_msg =\
+        "Calculation of tensor from Y does not compare to analytical solution")
+
+        assert_almost_equal(tensor_fit.md[0], md)
+
+        assert_array_almost_equal(tensor_fit.directions.shape[0], 3)
+        
+    # Test error-handling:
+    assert_raises(ValueError,
+                  dti.TensorModel,
+                  gtab,
+                  fit_method='crazy_method')
+
+    
+    
 def test_tensor_scalar_attributes():
     """
     Tests that the tensor class scalar attributes (FA, ADC, etc...) are
@@ -54,6 +113,7 @@ def test_tensor_scalar_attributes():
     #assert_array_equal(n_list % 2, 0)
     #assert_raises(ValueError, qball.sph_harm_ind_list, 1)
 
+
 def test_fa_of_zero():
     dummy_gtab = np.zeros((10,3))
     dummy_bval = np.zeros((10,))
@@ -62,6 +122,7 @@ def test_fa_of_zero():
     assert_equal(ten.fa(), 0)
     assert_true(np.isnan(ten.fa(nonans=False)))
 
+
 def test_WLS_and_LS_fit():
     """
     Tests the WLS and LS fitting functions to see if they returns the correct
@@ -114,6 +175,7 @@ def test_WLS_and_LS_fit():
     assert_almost_equal(tensor_est.md(), md)
     assert_array_almost_equal(tensor_est.lower_triangular(b0), D)
 
+
 def test_masked_array_with_Tensor():
     data = np.ones((2,4,56))
     mask = np.array([[True, False, False, True],
@@ -144,6 +206,7 @@ def test_masked_array_with_Tensor():
     assert_equal(tensor.evecs.shape, (3,3))
     assert_equal(type(tensor.model_params), np.ndarray)
 
+
 def test_passing_maskedview():
     data = np.ones((2,4,56))
     mask = np.array([[True, False, False, True],
@@ -177,6 +240,7 @@ def test_passing_maskedview():
     assert_equal(tensor.evecs.shape, (3,3))
     assert_equal(type(tensor.model_params), np.ndarray)
 
+
 def test_init():
     data = np.ones((2,4,56))
     mask = np.ones((2,4),'bool')
@@ -193,6 +257,7 @@ def test_init():
     assert_raises(ValueError, dti.Tensor, data, bval, gtab.T,
                         fit_method=0)
 
+
 def test_lower_triangular():
     tensor = np.arange(9).reshape((3,3))
     D = lower_triangular(tensor)
@@ -212,6 +277,7 @@ def test_lower_triangular():
     result[:] = [0, 3, 4, 6, 7, 8, 0]
     assert_array_equal(D, result)
 
+
 def test_from_lower_triangular():
     result = np.array([[0, 1, 3],
                        [1, 2, 4],
diff --git a/dipy/reconst/tests/test_gqi.py b/dipy/reconst/tests/test_gqi.py
index 130bd70777..6aa6d3b8ff 100644
--- a/dipy/reconst/tests/test_gqi.py
+++ b/dipy/reconst/tests/test_gqi.py
@@ -4,7 +4,7 @@
 from dipy.core.gradients import gradient_table
 from dipy.sims.voxel import SticksAndBall
 from dipy.reconst.gqi import GeneralizedQSamplingModel
-from dipy.utils.spheremakers import sphere_vf_from
+from dipy.data import get_sphere
 from numpy.testing import (assert_equal, 
                            assert_almost_equal, 
                            run_module_suite)
@@ -16,8 +16,7 @@
 
 def test_gqi():
     #load symmetric 724 sphere
-    vertices, faces = sphere_vf_from('symmetric724')
-    sphere = Sphere(xyz=vertices)
+    sphere = get_sphere('symmetric724')
     #load icosahedron sphere
     sphere2 = create_unit_sphere(5)
     btable = np.loadtxt(get_data('dsi515btable'))    
diff --git a/dipy/reconst/tests/test_multi_voxel.py b/dipy/reconst/tests/test_multi_voxel.py
new file mode 100644
index 0000000000..d98a36d8c3
--- /dev/null
+++ b/dipy/reconst/tests/test_multi_voxel.py
@@ -0,0 +1,100 @@
+import numpy as np
+import numpy.testing as npt
+from dipy.reconst.multi_voxel import _squash, multi_voxel_model, CallableArray
+from dipy.reconst.shm import QballOdfModel
+from dipy.core.sphere import unit_icosahedron
+
+
+def test_squash():
+    A = np.ones((3, 3), dtype=float)
+    B = np.asarray(A, object)
+    npt.assert_array_equal(A, _squash(B))
+
+    B[2, 2] = None
+    A[2, 2] = 0
+    npt.assert_array_equal(A, _squash(B))
+
+    for ijk in np.ndindex(*B.shape):
+        B[ijk] = np.ones((2,))
+    A = np.ones((3, 3, 2))
+    npt.assert_array_equal(A, _squash(B))
+
+    B[2, 2] = None
+    A[2, 2] = 0
+    npt.assert_array_equal(A, _squash(B))
+
+
+def test_CallableArray():
+    callarray = CallableArray((2, 3), dtype=object)
+
+    # Test without Nones
+    callarray[:] = range
+    expected = np.empty([2, 3, 4])
+    expected[:] = range(4)
+    npt.assert_array_equal(callarray(4), expected)
+
+    # Test with Nones
+    callarray[0, 0] = None
+    expected[0, 0] = 0
+    npt.assert_array_equal(callarray(4), expected)
+
+
+def test_multi_voxel_model():
+
+    class SillyModel(object):
+
+        def fit(self, data, mask=None):
+            return SillyFit(model)
+
+    class SillyFit(object):
+
+        def __init__(self, model):
+            self.model = model
+
+        model_attr = 2.
+
+        def odf(self, sphere):
+            return np.ones(len(sphere.phi))
+
+        @property
+        def directions(self):
+            n = np.random.randint(0, 10)
+            return np.zeros((n, 3))
+
+    # Wrap the SillyModel
+    MultiVoxelSillyModel = multi_voxel_model(SillyModel)
+
+    # Test the single voxel case
+    model = MultiVoxelSillyModel()
+    single_voxel = np.zeros(64)
+    fit = model.fit(single_voxel)
+    npt.assert_equal(type(fit), SillyFit)
+
+    # Test without a mask
+    many_voxels = np.zeros((2, 3, 4, 64))
+    fit = model.fit(many_voxels)
+    expected = np.empty((2, 3, 4))
+    expected[:] = 2.
+    npt.assert_array_equal(fit.model_attr, expected)
+    expected = np.ones((2, 3, 4, 12))
+    npt.assert_array_equal(fit.odf(unit_icosahedron), expected)
+    npt.assert_equal(fit.directions.shape, (2, 3, 4))
+
+    # Test with a mask
+    mask = np.eye(3).astype('bool')
+    data = np.zeros((3, 3, 64))
+    fit = model.fit(data, mask)
+    npt.assert_array_equal(fit.model_attr, np.eye(3)*2)
+    odf = fit.odf(unit_icosahedron)
+    npt.assert_equal(odf.shape, (3, 3, 12))
+    npt.assert_array_equal(odf[~mask], 0)
+    npt.assert_array_equal(odf[mask], 1)
+
+    # Test fit.shape
+    npt.assert_equal(fit.shape, (3, 3))
+
+    # Test indexing into a fit
+    npt.assert_equal(type(fit[0, 0]), SillyFit)
+    npt.assert_equal(fit[:2, :2].shape, (2, 2))
+
+
diff --git a/dipy/reconst/tests/test_odf.py b/dipy/reconst/tests/test_odf.py
index 4b88eb7926..3721f65f0f 100644
--- a/dipy/reconst/tests/test_odf.py
+++ b/dipy/reconst/tests/test_odf.py
@@ -85,7 +85,7 @@ def test_peaksFromModel():
     # Test basic case
     model = SimpleOdfModel()
     odf_argmax = _odf.argmax()
-    pam = peaks_from_model(model, data, normalize_peaks=True)
+    pam = peaks_from_model(model, data, _sphere, .5, 45, normalize_peaks=True)
 
     assert_array_equal(pam.gfa, gfa(_odf))
     assert_array_equal(pam.peak_values[:, 0], 1.)
@@ -97,21 +97,22 @@ def test_peaksFromModel():
     assert_array_equal(pam.peak_indices[:, 1:], -1)
 
     # Test that odf array matches and is right shape
-    pam = peaks_from_model(model, data, return_odf=True)
+    pam = peaks_from_model(model, data, _sphere, .5, 45, return_odf=True)
     expected_shape = (len(data), len(_odf))
     assert_equal(pam.odf.shape, expected_shape)
     assert_true((_odf == pam.odf).all())
     assert_array_equal(pam.peak_values[:, 0], _odf.max())
-   
+
     # Test mask
     mask = (np.arange(10) % 2) == 1
 
-    pam = peaks_from_model(model, data, mask=mask, normalize_peaks=True)
+    pam = peaks_from_model(model, data, _sphere, .5, 45, mask=mask,
+                           normalize_peaks=True)
     assert_array_equal(pam.gfa[~mask], 0)
     assert_array_equal(pam.qa[~mask], 0)
     assert_array_equal(pam.peak_values[~mask], 0)
     assert_array_equal(pam.peak_indices[~mask], -1)
-    
+
     assert_array_equal(pam.gfa[mask], gfa(_odf))
     assert_array_equal(pam.peak_values[mask, 0], 1.)
     assert_array_equal(pam.peak_values[mask, 1:], 0.)
diff --git a/dipy/reconst/tests/test_peak_finding.py b/dipy/reconst/tests/test_peak_finding.py
index 1c19e05f39..d01427187d 100644
--- a/dipy/reconst/tests/test_peak_finding.py
+++ b/dipy/reconst/tests/test_peak_finding.py
@@ -1,14 +1,13 @@
 import numpy as np
-from nose.tools import assert_true, assert_false, assert_equal, assert_almost_equal, assert_raises
-from numpy.testing import assert_array_equal
-from dipy.reconst.recspeed import (peak_finding, local_maxima, _filter_peaks,
-                                   remove_similar_vertices)
+import numpy.testing as npt
+from dipy.reconst.recspeed import local_maxima, remove_similar_vertices
 from dipy.data import get_sphere, get_data
 from dipy.core.sphere import unique_edges, HemiSphere
 from dipy.sims.voxel import all_tensor_evecs, multi_tensor_odf
 
 def test_local_maxima():
-    vertices, faces=get_sphere('symmetric724')
+    sphere = get_sphere('symmetric724')
+    vertices, faces = sphere.vertices, sphere.faces
     edges = unique_edges(faces)
     odf = abs(vertices.sum(-1))
     odf[1] = 10.
@@ -16,8 +15,8 @@ def test_local_maxima():
     odf[505] = 505
 
     peak_values, peak_index = local_maxima(odf, edges)
-    assert_array_equal(peak_values, [505, 143, 10])
-    assert_array_equal(peak_index, [505, 143, 1])
+    npt.assert_array_equal(peak_values, [505, 143, 10])
+    npt.assert_array_equal(peak_index, [505, 143, 1])
 
     hemisphere = HemiSphere(xyz=vertices, faces=faces)
     vertices_half, edges_half = hemisphere.vertices, hemisphere.edges
@@ -26,45 +25,12 @@ def test_local_maxima():
     odf[143] = 143.
 
     peak_value, peak_index = local_maxima(odf, edges_half)
-    assert_array_equal(peak_value, [143, 10])
-    assert_array_equal(peak_index, [143, 1])
+    npt.assert_array_equal(peak_value, [143, 10])
+    npt.assert_array_equal(peak_index, [143, 1])
 
     odf[20] = np.nan
-    assert_raises(ValueError, local_maxima, odf, edges_half)
+    npt.assert_raises(ValueError, local_maxima, odf, edges_half)
 
-def test_peak_finding():
-
-    vertices, faces=get_sphere('symmetric724')
-    odf=np.zeros(len(vertices))
-    odf = np.abs(vertices.sum(-1))
-
-    odf[1] = 10.
-    odf[505] = 505.
-    odf[143] = 143.
-
-    peaks, inds=peak_finding(odf.astype('f8'), faces.astype('uint16'))
-    print peaks, inds
-    edges = unique_edges(faces)
-    peaks, inds = local_maxima(odf, edges)
-    print peaks, inds
-
-    hemisphere = HemiSphere(xyz=vertices, faces=faces)
-    vertices_half, edges_half = hemisphere.vertices, hemisphere.edges
-    n = len(vertices_half)
-    peaks, inds = local_maxima(odf[:n], edges_half)
-    print peaks, inds
-    mevals=np.array(([0.0015,0.0003,0.0003],
-                    [0.0015,0.0003,0.0003]))
-    e0=np.array([1,0,0.])
-    e1=np.array([0.,1,0])
-    mevecs=[all_tensor_evecs(e0),all_tensor_evecs(e1)]
-    odf = multi_tensor_odf(vertices, [0.5,0.5], mevals, mevecs)
-    peaks, inds=peak_finding(odf, faces)
-    print peaks, inds
-    peaks2, inds2 = local_maxima(odf[:n], edges_half)
-    print peaks2, inds2
-    assert_equal(len(peaks), 2)
-    assert_equal(len(peaks2), 2)
 
 def test_remove_similar_peaks():
     vertices = np.array([[1., 0., 0.],
@@ -77,76 +43,41 @@ def test_remove_similar_peaks():
     norms = np.sqrt((vertices*vertices).sum(-1))
     vertices = vertices/norms[:, None]
 
-    uv, mapping = remove_similar_vertices(vertices, .01)
-    assert_array_equal(uv, vertices[:6])
-    assert_array_equal(mapping, range(6) + [0])
-    uv, mapping = remove_similar_vertices(vertices, 30)
-    assert_array_equal(uv, vertices[:4])
-    assert_array_equal(mapping, range(4) + [1, 0, 0])
-    uv, mapping = remove_similar_vertices(vertices, 60)
-    assert_array_equal(uv, vertices[:3])
-    assert_array_equal(mapping, range(3) + [0, 1, 0, 0])
-
-
-def test_filter_peaks():
-    # The setup
-    peak_values = np.array([1, .9, .8, .7, .6, .2, .1])
-    peak_points = np.array([[1., 0., 0.],
-                            [0., 0., 1.],
-                            [0., .9, .1],
-                            [0., 0., 1.],
-                            [0., 1., 0.],
-                            [0., 0., 1.],
-                            [.9, .1, 0.],
-                            [0., 0., 1.],
-                            [0., 0., 1.],
-                            [0., 0., 1.],
-                            [1., 1., 0.],
-                            [0., 0., 1.],
-                            [0., 1., 1.]])
-    norms = np.sqrt((peak_points*peak_points).sum(-1))
-    peak_points = peak_points/norms[:, None]
-
-    # Filter above peaks down to thre peaks
-    copy_peak_values = peak_values.copy()
-    ind = np.arange(0, len(peak_points), 2, dtype='int')
-    copy_ind = ind.copy()
-    print ind
-    sep_mat = abs(np.dot(peak_points, peak_points.T))
-    fvalues, find = _filter_peaks(peak_values, ind, sep_mat, .5, .9)
-    assert_array_equal(find, [0,2,8])
-    assert_array_equal(fvalues, [1., .9, .6])
-    # Check that the arguments have not been modified by _filter_peaks
-    assert_array_equal(ind, copy_ind)
-    assert_array_equal(peak_values, copy_peak_values)
-    # Test on a larger set of peaks
-    v, faces=get_sphere('symmetric724')
-    sep_mat = np.dot(v, v.T)
-    values = np.arange(len(v), 0., -1.)
-    ind = np.arange(len(values), dtype='int')
-    # Filter nothing if thresholds are 0. and 1.
-    fvalues, find = _filter_peaks(values, ind, sep_mat, 0., 1.)
-    assert_array_equal(find, ind)
-    assert_array_equal(fvalues, values)
-    # Return only the largest peak.
-    fvalues, find = _filter_peaks(values, ind, sep_mat, 0., -1.1)
-    assert_array_equal(find, [0])
-    assert_array_equal(fvalues, [len(values)])
-    fvalues, find = _filter_peaks(values, ind, sep_mat, 1., 0.)
-    assert_array_equal(find, [0])
-    assert_array_equal(fvalues, [len(values)])
-
-    fvalues, find = _filter_peaks(values, ind, sep_mat, .5, 1.)
-    assert_array_equal(fvalues, values[values >= .5*values[0]])
-    assert_array_equal(find, ind[values >= .5*values[0]])
-    values = values[1:].copy()
-    ind = ind[1:].copy()
-    sep_mat = sep_mat[1:, 1:].copy()
-    fvalues, find = _filter_peaks(values, ind, sep_mat, .5, 1.)
-    assert_array_equal(fvalues, values[values >= .5*values[0]])
-    assert_array_equal(find, ind[values >= .5*values[0]])
+    # Return unique vertices
+    uv = remove_similar_vertices(vertices, .01)
+    npt.assert_array_equal(uv, vertices[:6])
+
+    # Return vertices with mapping and indices
+    uv, mapping, index = remove_similar_vertices(vertices, .01,
+                                                 return_mapping=True,
+                                                 return_index=True)
+    npt.assert_array_equal(uv, vertices[:6])
+    npt.assert_array_equal(mapping, range(6) + [0])
+    npt.assert_array_equal(index, range(6))
+
+    # Test mapping with different angles
+    uv, mapping = remove_similar_vertices(vertices, .01, return_mapping=True)
+    npt.assert_array_equal(uv, vertices[:6])
+    npt.assert_array_equal(mapping, range(6) + [0])
+    uv, mapping = remove_similar_vertices(vertices, 30, return_mapping=True)
+    npt.assert_array_equal(uv, vertices[:4])
+    npt.assert_array_equal(mapping, range(4) + [1, 0, 0])
+    uv, mapping = remove_similar_vertices(vertices, 60, return_mapping=True)
+    npt.assert_array_equal(uv, vertices[:3])
+    npt.assert_array_equal(mapping, range(3) + [0, 1, 0, 0])
+
+    # Test index with different angles
+    uv, index = remove_similar_vertices(vertices, .01, return_index=True)
+    npt.assert_array_equal(uv, vertices[:6])
+    npt.assert_array_equal(index, range(6))
+    uv, index = remove_similar_vertices(vertices, 30, return_index=True)
+    npt.assert_array_equal(uv, vertices[:4])
+    npt.assert_array_equal(index, range(4))
+    uv, index = remove_similar_vertices(vertices, 60, return_index=True)
+    npt.assert_array_equal(uv, vertices[:3])
+    npt.assert_array_equal(index, range(3))
 
 
 if __name__ == '__main__':
-    test_peak_finding()
-
+    import nose
+    nose.runmodule()
diff --git a/dipy/reconst/tests/test_sphere_max.py b/dipy/reconst/tests/test_sphere_max.py
deleted file mode 100644
index df73e87b8a..0000000000
--- a/dipy/reconst/tests/test_sphere_max.py
+++ /dev/null
@@ -1,196 +0,0 @@
-""" Testing sphere maxima finding and associated routines
-"""
-
-from os.path import join as pjoin, dirname
-import numpy as np
-from dipy.data import get_sphere
-
-from dipy.core.meshes import (
-    sym_hemisphere,
-    neighbors,
-    vertinds_to_neighbors,
-    vertinds_faces,
-    argmax_from_adj,
-    peak_finding_compatible,
-    edges,
-    vertex_adjacencies)
-
-import dipy.reconst.recspeed as dcr
-
-from nose.tools import assert_true, assert_false, \
-     assert_equal, assert_raises
-
-from numpy.testing import assert_array_equal, assert_array_almost_equal
-
-# 8 faces (two square pyramids)
-VERTICES = np.array([
-        [0, 0, 1],
-        [1, 0, 0],
-        [0, 1, 0],
-        [-1, 0, 0],
-        [0, -1, 0],
-        [0, 0, -1]], dtype=np.float)
-FACES = np.array([
-        [0, 1, 2],
-        [0, 2, 3],
-        [0, 3, 4],
-        [0, 4, 1],
-        [5, 1, 2],
-        [5, 2, 3],
-        [5, 3, 4],
-        [5, 4, 1]])
-N_VERTICES = VERTICES.shape[0]
-VERTEX_INDS = np.array([0,1,2,3,4,5])
-
-DATA_PATH = pjoin(dirname(__file__), '..', 'matrices')
-# vertex, face tuple
-SPHERE_DATA = get_sphere('symmetric362')
-
-def test_sym_hemisphere():
-    assert_raises(ValueError, sym_hemisphere,
-                        VERTICES, 'k')
-    assert_raises(ValueError, sym_hemisphere,
-                        VERTICES, '%z')
-    for hem in ('x', 'y', 'z', '-x', '-y', '-z'):
-        vert_inds = sym_hemisphere(VERTICES, hem)
-        assert_equal(vert_inds.shape, (3,))
-        verts = VERTICES[vert_inds]
-        # there are no symmetrical points remanining in vertices
-        for vert in verts:
-            assert_false(np.any(np.all(
-                    vert * -1 == verts)))
-    # Test the sphere mesh data
-    vertices, _ = SPHERE_DATA
-    n_vertices = vertices.shape[0]
-    vert_inds = sym_hemisphere(vertices)
-    assert_array_equal(vert_inds,
-                       np.arange(n_vertices / 2))
-
-
-def test_vertinds_neighbors():
-    adj = neighbors(FACES)
-    assert_array_equal(adj,
-                             [[1, 2, 3, 4],
-                              [0, 2, 4, 5],
-                              [0, 1, 3, 5],
-                              [0, 2, 4, 5],
-                              [0, 1, 3, 5],
-                              [1, 2, 3, 4]])
-    adj = vertinds_to_neighbors(np.arange(6),
-                                FACES)
-    assert_array_equal(adj,
-                             [[1, 2, 3, 4],
-                              [0, 2, 4, 5],
-                              [0, 1, 3, 5],
-                              [0, 2, 4, 5],
-                              [0, 1, 3, 5],
-                              [1, 2, 3, 4]])
-    # subset of inds gives subset of faces
-    adj = vertinds_to_neighbors(np.arange(3),
-                                      FACES)
-    assert_array_equal(adj,
-                             [[1, 2, 3, 4],
-                              [0, 2, 4, 5],
-                              [0, 1, 3, 5]])
-    # can be any subset
-    adj = vertinds_to_neighbors(np.arange(3,6),
-                                      FACES)
-    assert_array_equal(adj,
-                             [[0, 2, 4, 5],
-                              [0, 1, 3, 5],
-                              [1, 2, 3, 4]])
-    # just test right size for the real mesh
-    vertices, faces = SPHERE_DATA
-    n_vertices = vertices.shape[0]
-    adj = vertinds_to_neighbors(np.arange(n_vertices),
-                                      faces)
-    assert_equal(len(adj), n_vertices)
-    assert_equal(len(adj[1]), 6)
-
-
-def test_vertinds_faces():
-    # routines to strip out faces
-    f2 = vertinds_faces(range(6), FACES)
-    assert_array_equal(f2, FACES)
-    f2 = vertinds_faces([0, 5], FACES)
-    assert_array_equal(f2, FACES)
-    f2 = vertinds_faces([0], FACES)
-    assert_array_equal(f2, FACES[:4])
-
-
-def test_neighbor_max():
-    # test ability to find maxima on sphere using neighbors
-    vert_inds = sym_hemisphere(VERTICES)
-    adj_inds = vertinds_to_neighbors(vert_inds, FACES)
-    # test slow and fast routine
-    for func in (argmax_from_adj, dcr.argmax_from_adj):
-        # all equal, no maxima
-        vert_vals = np.zeros((N_VERTICES,))
-        inds = func(vert_vals,
-                    vert_inds,
-                    adj_inds)
-        assert_equal(inds.size, 0)
-        # just ome max
-        for max_pos in range(3):
-            vert_vals = np.zeros((N_VERTICES,))
-            vert_vals[max_pos] = 1
-            inds = func(vert_vals,
-                        vert_inds,
-                        adj_inds)
-            assert_array_equal(inds, [max_pos])
-        # maxima outside hemisphere don't appear
-        for max_pos in range(3,6):
-            vert_vals = np.zeros((N_VERTICES,))
-            vert_vals[max_pos] = 1
-            inds = func(vert_vals,
-                        vert_inds,
-                        adj_inds)
-            assert_equal(inds.size, 0)
-        # use whole mesh, with two maxima
-        w_vert_inds = np.arange(6)
-        w_adj_inds = vertinds_to_neighbors(w_vert_inds, FACES)
-        vert_vals = np.array([1.0, 0, 0, 0, 0, 2])
-        inds = func(vert_vals, w_vert_inds, w_adj_inds)
-        assert_array_equal(inds, [0, 5])
-        # check too few vals raises sensible error.  For the Cython
-        # version of the routine, the test below causes odd errors and
-        # segfaults with numpy SVN vintage June 2010 (sometime after
-        # 1.4.0 release) - see
-        # http://groups.google.com/group/cython-users/browse_thread/thread/624c696293b7fe44?pli=1
-        # assert_raises(IndexError, func, vert_vals[:3],
-        # w_vert_inds, w_adj_inds)
-
-
-def test_performance():
-    # test this implementation against Frank Yeh implementation
-    vertices, faces = SPHERE_DATA
-    n_vertices = vertices.shape[0]
-    vert_inds = sym_hemisphere(vertices)
-    adj = vertinds_to_neighbors(vert_inds, faces)
-    np.random.seed(42)
-    vert_vals = np.random.uniform(size=(n_vertices,))
-    maxinds = argmax_from_adj(vert_vals, vert_inds, adj)
-    maxes, pfmaxinds = dcr.peak_finding(vert_vals, faces)
-    assert_array_equal(maxinds, pfmaxinds[::-1])
-
-
-def test_sym_check():
-    assert_true(peak_finding_compatible(VERTICES))
-    vertices, faces = SPHERE_DATA
-    assert_true(peak_finding_compatible(vertices))
-    assert_false(peak_finding_compatible(vertices[::-1]))
-
-
-def test_adjacencies():
-    faces = FACES
-    vertex_inds = VERTEX_INDS
-    edgearray = edges(vertex_inds, faces)
-    assert_array_equal(edgearray.shape,(24,2))
-    assert_array_equal(edgearray,
-                             [[3, 0], [5, 4], [2, 1], [5, 1],
-                              [2, 5], [0, 3], [4, 0], [1, 2],
-                              [1, 5], [0, 4], [5, 3], [4, 1],
-                              [3, 2], [4, 5], [1, 4], [2, 3],
-                              [1, 0], [3, 5], [0, 1], [5, 2],
-                              [2, 0] ,[4, 3], [3, 4], [0, 2]])
-    assert_array_equal(vertex_adjacencies(vertex_inds, faces).shape,(6,6))
diff --git a/dipy/sims/tests/test_voxel.py b/dipy/sims/tests/test_voxel.py
index cf32fb6339..020af964f8 100644
--- a/dipy/sims/tests/test_voxel.py
+++ b/dipy/sims/tests/test_voxel.py
@@ -46,7 +46,8 @@ def test_single_tensor():
 
 
 def test_multi_tensor():
-    vertices, faces = get_sphere('symmetric724')
+    sphere = get_sphere('symmetric724')
+    vertices, faces = sphere.vertices, sphere.faces
     mevals=np.array(([0.0015, 0.0003, 0.0003],
                      [0.0015, 0.0003, 0.0003]))
     e0 = np.array([1, 0, 0.])
diff --git a/dipy/sims/voxel.py b/dipy/sims/voxel.py
index cd8f94b8c1..66bd306a4f 100644
--- a/dipy/sims/voxel.py
+++ b/dipy/sims/voxel.py
@@ -2,7 +2,6 @@
 
 import numpy as np
 from dipy.core.geometry import sphere2cart
-from dipy.reconst.dti import design_matrix, lower_triangular
 from dipy.core.geometry import vec2vec_rotmat
 
 
@@ -249,7 +248,8 @@ def multi_tensor_odf(odf_verts, mf, mevals=None, mevecs=None):
     >>> import numpy as np
     >>> from dipy.sims.voxel import multi_tensor_odf, all_tensor_evecs
     >>> from dipy.data import get_sphere
-    >>> vertices, faces = get_sphere('symmetric724')
+    >>> sphere = get_sphere('symmetric724')
+    >>> vertices, faces = sphere.vertices, sphere.faces
     >>> mevals=np.array(([0.0015, 0.0003, 0.0003],[0.0015, 0.0003, 0.0003]))
     >>> e0 = np.array([1, 0, 0.])
     >>> e1 = np.array([0., 1, 0])
diff --git a/dipy/tracking/eudx.py b/dipy/tracking/eudx.py
index 428fc6accd..4d802d76de 100644
--- a/dipy/tracking/eudx.py
+++ b/dipy/tracking/eudx.py
@@ -119,8 +119,9 @@ def __init__(self, a, ind,
         x,y,z,g=self.a.shape
         self.Np=g
         if odf_vertices==None:
-            vertices, faces = get_sphere('symmetric362')
-            self.odf_vertices = np.ascontiguousarray(vertices, dtype='f8')
+            sphere = get_sphere('symmetric362')
+            vertices, faces = sphere.vertices, sphere.faces
+            self.odf_vertices = vertices
         else:
             self.odf_vertices = np.ascontiguousarray(odf_vertices, dtype='f8')
         try:
diff --git a/dipy/utils/spheremakers.py b/dipy/utils/spheremakers.py
index 5704ce7d86..0e8b98fa38 100644
--- a/dipy/utils/spheremakers.py
+++ b/dipy/utils/spheremakers.py
@@ -1,6 +1,8 @@
 """ Factory function(s) for spheres """
 
 from dipy.data import get_sphere
+from dipy.core.sphere import Sphere
+
 
 def sphere_vf_from(input):
     """ Return sphere vertices and faces from a variety of inputs
@@ -20,7 +22,7 @@ def sphere_vf_from(input):
         Indices into `vertices`
     """
     if hasattr(input, 'keys'):
-        return input['vertices'], input['faces']
+        return Sphere(xyz=input['vertices'])
     if isinstance(input, basestring):
         return get_sphere(input)
     return input
diff --git a/dipy/utils/tests/test_spheremakers.py b/dipy/utils/tests/test_spheremakers.py
index ecb25ef6d9..a26737fc97 100644
--- a/dipy/utils/tests/test_spheremakers.py
+++ b/dipy/utils/tests/test_spheremakers.py
@@ -13,7 +13,8 @@
 def test_spheremakers():
     # Test inputs to spheremakers
     # Example data given string
-    v, f = sphere_vf_from('symmetric362')
+    sphere = sphere_vf_from('symmetric362')
+    v, f = sphere.vertices, sphere.faces
     assert_equal(f.shape[1], 3)
     assert_equal(v.shape[1], 3)
     # Given tuple
@@ -21,7 +22,7 @@ def test_spheremakers():
     assert_array_equal(vdash, v)
     assert_array_equal(fdash, f)
     # Given dict
-    vdash, fdash = sphere_vf_from({'vertices': v,
-                                   'faces': f})
-    assert_array_equal(vdash, v)
-    assert_array_equal(fdash, f)
+    sphere = sphere_vf_from({'vertices': v})
+    vdash, fdash = sphere.vertices, sphere.faces
+    assert_array_almost_equal(vdash, v)
+    assert_array_almost_equal(fdash, f)
diff --git a/dipy/viz/__init__.py b/dipy/viz/__init__.py
index 542cec1701..fc9beb65ed 100644
--- a/dipy/viz/__init__.py
+++ b/dipy/viz/__init__.py
@@ -30,4 +30,4 @@
     from ._show_odfs import show_odfs
 
 if has_mpl:
-    from projections import *
+    import projections
diff --git a/dipy/viz/projections.py b/dipy/viz/projections.py
index 1539441df3..d3fb80aafe 100644
--- a/dipy/viz/projections.py
+++ b/dipy/viz/projections.py
@@ -12,13 +12,13 @@
 matplotlib, has_mpl, setup_module = optional_package("matplotlib")
 plt, _, _ = optional_package("matplotlib.pyplot")
 tri, _, _ = optional_package("matplotlib.tri")
+bm, _, _ = optional_package("mpl_toolkits.basemap")
 
 import dipy.core.geometry as geo
 
 
-def sph_project(vertices, val, ax=None, vmin=None, vmax=None,
-                cmap=None, cbar=True, triang=False):
-
+def sph_project(vertices, val, ax=None, vmin=None, vmax=None, cmap=None,
+                cbar=True, tri=False, boundary=False, **basemap_args):
     """Draw a signal on a 2D projection of the sphere.
 
     Parameters
@@ -42,26 +42,49 @@ def sph_project(vertices, val, ax=None, vmin=None, vmax=None,
 
     triang: Whether to display the plot triangulated as a pseudo-color plot.
 
+    boundary: Whether to draw the boundary around the projection in a black line
+
     Returns
     -------
-    fig : figure
-        Matplotlib figure
+    ax : axis
+        Matplotlib figure axis
 
     Examples
     --------
     >>> from dipy.data import get_sphere
-    >>> verts,faces=get_sphere('symmetric724')
-    >>> ax = sph_project(verts,np.random.rand(len(verts)))
+    >>> verts = get_sphere('symmetric724').vertices
+    >>> ax = sph_project(verts.T, np.random.rand(len(verts.T)))
 
     """
+    if ax is None:
+        fig, ax = plt.subplots(1)
+
     if cmap is None:
         cmap = matplotlib.cm.hot
-    if ax is None:
-        _, ax = plt.subplots(1)
-    fig = ax.get_figure()
 
-    x = vertices[:, 0]
-    y = vertices[:, 1]
+    basemap_args.setdefault('projection', 'ortho')
+    basemap_args.setdefault('lat_0', 0)
+    basemap_args.setdefault('lon_0', 0)
+    basemap_args.setdefault('resolution', 'c')
+
+    from mpl_toolkits.basemap import Basemap
+
+    m = Basemap(**basemap_args)
+    if boundary:
+        m.drawmapboundary()
+
+    # Rotate the coordinate system so that you are looking from the north pole:
+    verts_rot = np.array(np.dot(np.matrix([[0,0,-1],[0,1,0],[1,0,0]]), vertices))
+
+    # To get the orthographic projection, when the first coordinate is positive:
+    neg_idx = np.where(verts_rot[0]>0)
+
+    # rotate the entire bvector around to point in the other direction:
+    verts_rot[:, neg_idx] *= -1
+
+    _, theta, phi = geo.cart2sphere(verts_rot[0], verts_rot[1], verts_rot[2])
+    lat, lon = geo.sph2latlon(theta, phi)
+    x, y = m(lon, lat)
 
     my_min = np.nanmin(val)
     if vmin is not None:
@@ -71,34 +94,30 @@ def sph_project(vertices, val, ax=None, vmin=None, vmax=None,
     if vmax is not None:
         my_max = vmax
 
-    r = (val - my_min)/float(my_max-my_min)
-
-    # Enforce the maximum and minumum boundaries, if there are values
-    # outside those boundaries:
-    r[r<0]=0
-    r[r>1]=1
+    if tri:
+        m.pcolor(x, y, val, vmin=my_min, vmax=my_max, tri=True, cmap=cmap)
 
-    if triang:
-        triang = tri.Triangulation(x, y)
-        ax.tripcolor(triang, r, cmap=cmap)
     else:
         cmap_data = cmap._segmentdata
         red_interp, blue_interp, green_interp = (
-                interp.interp1d(np.array(cmap_data[gun])[:,0],
-                                np.array(cmap_data[gun])[:,1]) for gun in
-                                ['red', 'blue','green'])
+        interp.interp1d(np.array(cmap_data[gun])[:,0],
+                        np.array(cmap_data[gun])[:,1]) for gun in
+                                                  ['red', 'blue','green'])
 
+        r = (val - my_min)/float(my_max-my_min)
+
+        # Enforce the maximum and minumum boundaries, if there are values
+        # outside those boundaries:
+        r[r<0]=0
+        r[r>1]=1
 
         for this_x, this_y, this_r in zip(x,y,r):
             red = red_interp(this_r)
             blue = blue_interp(this_r)
             green = green_interp(this_r)
-            ax.plot(this_x, this_y, 'o',
-                    c=[red.item(), green.item(), blue.item()])
-
+            m.plot(this_x, this_y, 'o',
+                   c=[red.item(), green.item(), blue.item()])
 
-    ax.set_aspect('equal')
-    ax.set_axis_off()
     if cbar:
         mappable = matplotlib.cm.ScalarMappable(cmap=cmap)
         mappable.set_array([my_min, my_max])
@@ -109,7 +128,5 @@ def sph_project(vertices, val, ax=None, vmin=None, vmax=None,
         cax = fig.add_axes([l+w+0.075, b, 0.05, h], frameon=False)
         fig.colorbar(mappable, cax=cax) # draw colorbar
 
-    ax.set_xlim([-1.1, 1.1])
-    ax.set_ylim([-1.1, 1.1])
-
     return ax
+
diff --git a/doc/examples/multiVoxelModel.py b/doc/examples/multiVoxelModel.py
new file mode 100644
index 0000000000..de5426a8b6
--- /dev/null
+++ b/doc/examples/multiVoxelModel.py
@@ -0,0 +1,85 @@
+"""The `multi_voxel_model` is a class decorator to help easily write multi voxel
+models. A developer simply needs to write a description of the single voxel
+case and wrap it using `multi_voxel_model` like bellow."""
+
+import numpy as np
+from dipy.core.sphere import unit_icosahedron
+from dipy.reconst.multi_vox import multi_voxel_model
+
+"""First the developer should write out the single voxel as bellow and either
+wrap or decorate the Model Class"""
+
+class SingleVoxelModel(object):
+    def fit(self, data, mask=None):
+        n = np.max(data)
+        return SingleVoxelFit(self, n)
+
+class SingleVoxelFit(object):
+    model_attr = 1.0
+    def __init__(self, model, n):
+        self.model = model
+        self.n = n
+
+    def odf(self, sphere):
+        return np.ones(len(sphere.phi))
+
+    @property
+    def directions(self):
+        return np.zeros((self.n, 3))
+
+MultiVoxelModel = multi_voxel_model(SingleVoxelModel)
+
+@multi_voxel_model
+class DecoratedModel(SingleVoxelModel):
+    """Now to show how all this works
+
+    To show how the single voxel case works
+    ---------------------------------------
+    >>> model = SingleVoxelModel()
+    >>> fit = model.fit(4)
+    >>> fit.model_attr
+    1.0
+    >>> fit.directions.shape
+    (4, 3)
+    >>> fit.odf(unit_icosahedron).shape
+    (12,)
+
+
+    Now we use the MultiVoxelModel
+    ------------------------------
+    >>> model = MultiVoxelModel()
+    >>> data = np.arange(1, 6).reshape((2, 3, 1))
+    >>> fit = model.fit(data)
+    >>> fit.model_attr.shape
+    (2, 3)
+    >>> np.all(fit.model_attr == 1.)
+    True
+    >>> fit.directions.shape
+    (2, 3)
+    >>> fit.directions[0, 0].shape
+    (1, 3)
+    >>> fit.odf(unit_icosahedron).shape
+    (2, 3, 12)
+
+
+    Of course using using `multi_voxel_model` as a decorator or as a wrapper
+    function is exactly the same
+    -------------------------------------------------------------------------
+    >>> model = DecoratedModel()
+    >>> data = np.arange(1, 6).reshape((2, 3, 1))
+    >>> fit = model.fit(data)
+    >>> fit.model_attr.shape
+    (2, 3)
+    >>> np.all(fit.model_attr == 1.)
+    True
+    >>> fit.directions.shape
+    (2, 3)
+    >>> fit.directions[0, 0].shape
+    (1, 3)
+    >>> fit.odf(unit_icosahedron).shape
+    (2, 3, 12)
+
+
+    """
+    pass
+
diff --git a/doc/faq.rst b/doc/faq.rst
index fdb8fba10a..e4fa1d7f32 100644
--- a/doc/faq.rst
+++ b/doc/faq.rst
@@ -17,27 +17,28 @@ Theoretical
   pulsed gradient spin-echo (PGSE) sequence, at the time of readout
   $b=\gamma^{2}G^{2}\delta^{2}\left(\Delta-\frac{\delta}{3}\right)$
   where $\gamma$ is the gyromagnetic radio, $\delta$ denotes the pulse
-  width, $G$ is the gradient amplitude and $\Delta$ the centre to
-  centre spacing. $\gamma$ is a constant, but we can change the other
+  width, $G$ is the gradient amplitude and $\Delta$ the centre-to-centre 
+  spacing. $\gamma$ is a constant, but we can change the other
   three parameters and in that way control the b-value.
 
 2. **What is q-space?**
 
   Q-space is the space of one or more 3D spin displacement wave vectors
-  $\mathbf{q}$ as shown in equation \ref{eq:fourier}. The vector $\mathbf{q}$
+  $\mathbf{q}$ as shown in equation $\ref{eq:fourier}$. The vector $\mathbf{q}$
   parametrises the space of diffusion gradients. It is related to the
-  applied magnetic gradient $\mathbf{g}$ by the formula $\mathbf{q}=(2\pi)^{-1}\gamma\delta\mathbf{g}$.
+  applied magnetic gradient $\mathbf{g}$ by the formula 
+  $\mathbf{q}=(2\pi)^{-1}\gamma\delta\mathbf{g}$.
   Every single vector $\mathbf{q}$ has the same orientation as the
   direction of diffusion gradient $\mathbf{g}$ and length proportional
   to the strength $g$ of the gradient field. Every single point in
   q-space corresponds to a possible 3D volume of the MR signal for a specific
   gradient direction and strength. Therefore if, for example, we have
-  programmed the scanner to apply 60 gradient directions then our data
-  should have 60 diffusion volumes with each volume obtained for a specific
+  programmed the scanner to apply 60 gradient directions, then our data
+  should have 60 diffusion volumes, with each volume obtained for a specific
   gradient. A Diffusion Weighted Image (DWI) is the volume acquired
   from only one direction gradient.
 
-3. **What DWI stands for?**
+3. **What does DWI stand for?**
 
   Diffusion Weighted Imaging (DWI) is MRI imaging designed to be sensitive
   to diffusion. A diffusion weighted image is a volume of voxel data gathered 
@@ -46,81 +47,95 @@ Theoretical
   should be low if there is greater mobility of water molecules along
   the specified gradient direction and it should be high if there is
   less movement in that direction. Yes, it is counterintuitive but correct!
-  However greater mobility gives greater opportunity for the proton spins to be dephased
-  producing a smaller RF signal.
+  However, greater mobility gives greater opportunity for the proton spins to
+  be dephased, producing a smaller RF signal.
 
 4. **Why dMRI and not DTI?**
 
-  Diffusion MRI (dMRI or dwMRI) are prefered terms if you want to speak about diffusion weighted MRI in general. 
-  DTI (diffusion tensor imaging) is just one of the many ways you can reconstruct the voxel from your measured signal. 
-  There are plenty of others for example DSI, GQI, QBI etc.     
+  Diffusion MRI (dMRI or dwMRI) are the preferred terms if you want to speak 
+  about diffusion weighted MRI in general. DTI (diffusion tensor imaging) is 
+  just one of the many ways you can reconstruct the voxel from your measured 
+  signal. There are plenty of others, for example DSI, GQI, QBI, etc.
 
 5. **What is the recommended practice for registration of diffusion datasets?**
 
-  Registration can be tricky. But this is what usually works for us for normal healthy adult subjects. 
-  We register the FA (fractional anisotropy) images to the FMRIB_FA_1mm template which is in MNI space 
-  using ``flirt`` and ``fnirth`` from FSL. Then we can apply the warping displacements in any other scalar volumes 
-  that we have to register that scalar volume into the MNI space. We need the corresponding inverse displacements 
-  to map a tractography into MNI space. 
+  Registration can be tricky. But this is what usually works for us for normal 
+  healthy adult subjects. We register the FA (fractional anisotropy) images 
+  to the FMRIB_FA_1mm template which is in MNI space using ``flirt`` and 
+  ``fnirt`` from FSL. Then we can apply the warping displacements in any 
+  other scalar volumes that we have to register that scalar volume into the 
+  MNI space. We need the corresponding inverse displacements 
+  to map a tractography into MNI space.
 
 6. **What is the difference between Image coordinates and World coordinates?**
 
-  Image coordinates have positive integer values and represent the centres $(i, j, k)$ of the voxels. There is an affine transform 
-  (stored in the nifti file) that takes the image coordinates and transforms them to millimeter (mm) in real world space. 
-  World coordinates have floating point precision and your dataset have 3 real dimensions e.g. $(x, y, z)$.
+  Image coordinates have positive integer values and represent the centres 
+  $(i, j, k)$ of the voxels. There is an affine transform (stored in the 
+  nifti file) that takes the image coordinates and transforms them to 
+  millimeter (mm) in real world space. World coordinates have floating point 
+  precision and your dataset has 3 real dimensions e.g. $(x, y, z)$.
   
 7. **Why 'tracks' and not 'tracts'?**
 
-  Tractography is only an approximation or simulation - if you prefer - of the real tracts (brain neural fiber pathways 
-  or brain nerves). Therefore we prefer to call these simulated tracts as tracks (trajectories or curves represented as sequences of
-  points joined by line segments) so that others will be clear 
-  that they are not the real tracts (fibers) but only an estimate or suggestion. 
-  We hope that in the future tractography could reach a point that what you see on
-  your screen is a very faithful representation of what is actually in the white matter of the brain. 
-  However the field is not yet at this level of detail.    
+  Tractography is only an approximation or simulation - if you prefer - of 
+  the real tracts (brain neural fiber pathways or brain nerves). Therefore 
+  we prefer to call these simulated tracts as tracks (trajectories or curves 
+  represented as sequences of points joined by line segments) so that 
+  others will be clear that they are not the real tracts (fibers), but 
+  only an estimate or suggestion. We hope that in the future tractography 
+  could reach a point that what you see on your screen is a very faithful 
+  representation of what is actually in the white matter of the brain. 
+  However the field is not yet at this level of detail.
 
 8. **Why use 'deterministic' and not 'probabilistic' tractography?**
 
-  We wanted to create at the outset a tractographic method which will help us and you to get closer to datasets 
-  in a very efficient way. Therefore, we created first the ``EuDX`` (Euler Delta Crossings) algorithm which is a tracking method 
-  which can work both with model or model-free input and resolve also
-  crossing fibers with a high order of crossings. Also it is very fast to calculate (~2 minutes for 1 million tracks ). 
-  We hope that at a later stage we will be able to incorporate and test more methods e.g. probabilistic, global and graph-theoretic.
+  We wanted to create at the outset a tractographic method which will help us 
+  and you to get closer to datasets in a very efficient way. Therefore, we 
+  created first the ``EuDX`` (Euler Delta Crossings) algorithm which is a 
+  tracking method which can work both with model or model-free input and 
+  resolve also crossing fibers with a high order of crossings. Also it is 
+  very fast to calculate (~2 minutes for 1 million tracks). We hope that at 
+  a later stage we will be able to incorporate and test more methods e.g. 
+  probabilistic, global and graph-theoretic.
   
-9. **We made the mistake in our lab of generating datasets with nonisotropic voxel sizes wusehat do we do?**
+9. **We made the mistake in our lab of generating datasets with nonisotropic voxel sizes. What do we do?**
   
-  You need to resample your raw data to an isotropic size. Have a look at the module ``dipy.align.noniso2iso``. 
-  (We think it is a mistake to acquire nonisotropic data because the directional resolution of the data will depend on
-  the orientation of the gradient with respect to the voxels, being lower when aligned with a longer voxel dimension.)
+  You need to resample your raw data to an isotropic size. Have a look at 
+  the module ``dipy.align.noniso2iso``. (We think it is a mistake to 
+  acquire nonisotropic data because the directional resolution of the data 
+  will depend on the orientation of the gradient with respect to the 
+  voxels, being lower when aligned with a longer voxel dimension.)
   
-10. **Why nonisotropic voxel sizes are a bad idea in diffusion?**
+10. **Why are nonisotropic voxel sizes a bad idea in diffusion?**
   
-  If for example you have $2 \times 2 \times 4 \textrm{mm}^3$ voxels, the last dimension will
-  be averaged over the double distance and less detail will be captured compared
-  to the other two dimensions. Furthermore, with very nonisotropic voxels 
-  the uncertainty on orientation estimates will depend on the position of 
-  the subject in the scanner.
+  If, for example, you have $2 \times 2 \times 4\ \textrm{mm}^3$ voxels, the 
+  last dimension will be averaged over the double distance and less detail 
+  will be captured compared to the other two dimensions. Furthermore, with 
+  very nonisotropic voxels the uncertainty on orientation estimates will 
+  depend on the position of the subject in the scanner.
 
 ---------
 Practical
 ---------
 
-1. **Why python and not matlab or some other language?**
+1. **Why Python and not MATLAB or some other language?**
 
-  python is free, batteries included, very well designed,  painless to read and easy to use. 
+  Python is free, batteries included, very well-designed, painless to read 
+  and easy to use. 
   There is nothing else like it. Give it a go. 
-  Once with python always with python. 
+  Once with Python, always with Python. 
   
-2. **Isn't python slow?**
+2. **Isn't Python slow?**
 
-  True, some times python can be slow if you are using for example multiple nested for loops. 
-  In that case we use cython which takes execution up to C speed.
+  True, sometimes Python can be slow, if you are using multiple nested 
+  ``for`` loops, for example. 
+  In that case, we use Cython, which takes execution up to C speed.
 
-3. **What numerical libraries do you use in python?**
+3. **What numerical libraries do you use in Python?**
 
-  The best ever designed numerical library - numpy.
+  The best ever designed numerical library - NumPy.
 
-2. **Which python console do your recommend?**
+2. **Which Python console do you recommend?**
 
   ``ipython``
 
@@ -132,29 +147,34 @@ Practical
 
 4. **What about interactive visualization?**
 
-  There is already interaction in the ``fvtk`` module but we have started a new project 
-  only for visualization which we plan to integrate in ``dipy`` in the near future for more information 
-  have a look at http://fos.me
+  There is already interaction in the ``fvtk`` module, but we have started a 
+  new project only for visualization which we plan to integrate in ``dipy`` 
+  in the near future.  For more information, have a look at http://fos.me
 
 5. **Which file formats do you support?**
 
-  Nifti (.nii), Dicom (Siemens(read-only)), Trackvis (.trk), Dipy (.dpy), Numpy (.npy, ,npz), text 
-  and any other formats supported by nibabel and pydicom.
+  Nifti (.nii), Dicom (Siemens(read-only)), Trackvis (.trk), Dipy (.dpy), 
+  Numpy (.npy, ,npz), text and any other formats supported by nibabel and 
+  pydicom.
 
-  You can also read/save in Matlab version v4 (Level 1.0), v6 and v7 to 7.2 using scipy.io.loadmat. For higher versions >= 7.3
-  you can use pytables or any other python to hdf5 library e.g. h5py .
+  You can also read/save in Matlab version v4 (Level 1.0), v6 and v7 to 7.2, 
+  using `scipy.io.loadmat`. For higher versions >= 7.3, you can use pytables 
+  or any other python-to-hdf5 library e.g. h5py.
 
-  For object serialization you can used dipy.io.pickles function load_pickle, save_pickle.
+  For object serialization you can use `dipy.io.pickles` functions 
+  `load_pickle`, `save_pickle`.
 
 6. **What is dpy**?
 
-  ``dpy`` is an ``hdf5`` file format which we use in dipy to store tractography and other information. 
-  This allows us to store huge tractographies and load different parts of the datasets 
+  ``dpy`` is an ``hdf5`` file format which we use in dipy to store 
+  tractography and other information. This allows us to store huge 
+  tractographies and load different parts of the datasets 
   directly from the disk as if it were in memory.
 
 7. **Which python editor should I use?**
 
-  Any text editor would do the job but we prefer the following Aptana, Emacs, Vim and Eclipse (with PyDev).
+  Any text editor would do the job but we prefer the following: Aptana, 
+  Emacs, Vim and Eclipse (with PyDev).
 
 8. **I have problems reading my dicom files using nibabel, what should I do?**
 
diff --git a/doc/index.rst b/doc/index.rst
index da28cccb4d..628aa47da9 100644
--- a/doc/index.rst
+++ b/doc/index.rst
@@ -4,7 +4,7 @@
 Dipy
 ####
 
-Dipy_ is an *international*, **free** and **open soure** software project for
+Dipy_ is an *international*, **free** and **open source** software project for
 **diffusion** *magnetic resonance imaging* **analysis**.
 
 Depends on a few standard libraries: python_ (the core language), numpy_ (for