more doc, cleanup

Author: poolio

python/thunder/imgprocessing/regmethods/utils.py

@@ -1,6 +1,5 @@
 """ Shared utilities for registration methods """
 
-import numpy as np
 from numpy import ndarray
 
 from thunder.rdds.images import Images
@@ -110,6 +109,7 @@ def zeroBorder(vol, border=1, cval=0.0):
         border: scalar or tuple containing borders for each dimension
         cval: constant value to replace boundary voxels with, default is 0.0
     """
+    import numpy as np # sorry
     vol = vol.copy() # create new copy so we don't overwrite vol
     dims = np.array(vol.shape)
     if np.size(border) == 1:
@@ -130,62 +130,6 @@ def zeroBorder(vol, border=1, cval=0.0):
     return vol
 
 
-def imageGradients(im, sigma=None):
-    """Compute gradients of volume in each dimension using a Sobel filter.
-
-    Parameters
-    ----------
-    im : ndarray
-        single volume
-    sigma : float or tuple, optional, default = None
-        smoothing amount to apply to volume before computing gradients
-    """
-    from scipy.ndimage.filters import gaussian_filter, sobel
-    if sigma is not None:
-       im = gaussian_filter(im, sigma)
-    grads = [sobel(im, axis=dim, mode='constant') / 8.0 for dim in xrange(im.ndim)]
-    return grads
-
-
-def imageJacobian(vol, tfm, grid=None, sigma=None, normalize=True, border=1, order=1):
-    """Compute Jacobian of volume w.r.t. transformation parameters
-
-    Args:
-        I: volume
-        tfm: Transform object
-        sigma: smoothing bandwidth for gradients (None for no smoothing)
-        normalize: Whether to normalize images before aligning.
-        border: Number or tuple of border sizes to zero after transforming.
-    Returns:
-        It: vectorized version of transformed volume
-        J: nvoxels x nparameters Jacobian matrix
-    """
-    # Compute gradients
-    # Transform image and gradients
-    # Note: does not work if we compute gradients after transforming (borders?)
-    #It = zeroBorder(tfm.apply(I), border)
-    #TODO(ben): check that tfm then grad works
-#        It = transform_vol(I, tfm, grid, order=order)
-#        grads = [transform_vol(grad, tfm, grid, order=order) for grad in grads]
-    if grid is None:
-        grid = GridTransformer(vol.shape)
-    grads = imageGradients(vol, sigma)
-    It = tfm.apply(vol, grid, order=order)
-    It = zeroBorder(It, border)
-    grads = [tfm.apply(grad, grid, order=order) for grad in grads]
-    grads = [zeroBorder(grad, border) for grad in grads]
-
-    if normalize:
-        normI = np.linalg.norm(It.ravel())
-        if normI == 0.0:
-          raise ValueError('Transform yields volume of zeroes.')
-        grads = [(1. / normI) * I - (I * It).sum() / (normI**3) * It for I in grads]
-        It /= normI
-
-    jacobianVols = tfm.jacobian(grads, grid.homo_points)
-    return It, jacobianVols
-
-
 def solveLinearized(vec, jacobian, reference, robust=False):
     """Solve linearized registration problem for change in transformation parameters and weighting on reference basis.
 
@@ -214,22 +158,100 @@ def solveLinearized(vec, jacobian, reference, robust=False):
         model = QuantReg(vec, A).fit(q=quantile)
         params = model.params
     else:
-        model = np.linalg.lstsq(A, vec)
+        from numpy.linalg import lstsq
+        model = lstsq(A, vec)
         params = model[0]
-    deltaTransform, coeff = np.split(params, [jacobian.shape[1]])
+    from numpy import split
+    deltaTransform, coeff = split(params, [jacobian.shape[1]])
     return deltaTransform, coeff
 
 
-def volumesToMatrix(vols):
-    if not isinstance(vols, list):
-        return v2v(vols)
-    else:
-        return np.column_stack([v2v(v) for v in vols]).squeeze()
-
 def v2v(v, dims=None):
+    """Convert vector to volume or volume to vector.
+
+    Parameters
+    ----------
+    v : array
+        volume or vectorized version of volume
+    dims : tuple, optional
+        shape of volume, must be set to reshape vectors into volumes
+
+    Returns
+    -------
+    volume if input is a vector, and vector if input is a volume
+    """
     if v.ndim == 1:
+        from numpy import prod
         assert dims
         assert v.size == np.prod(dims)
         return v.reshape(dims)
     else:
-        return v.ravel()
+        return v.ravel()
+
+
+def volumesToMatrix(vols):
+    """Convert list of volumes to a matrix.
+
+    Parameters
+    ----------
+    vols : list of arrays
+
+    Returns
+    -------
+    array with size nvoxels by number of volumes
+    """
+    from numpy import column_stack
+    if not isinstance(vols, list):
+        return v2v(vols)
+    else:
+        return column_stack([v2v(v) for v in vols]).squeeze()
+
+
+def imageGradients(im, sigma=None):
+    """Compute gradients of volume in each dimension using a Sobel filter.
+
+    Parameters
+    ----------
+    im : ndarray
+        single volume
+    sigma : float or tuple, optional, default = None
+        smoothing amount to apply to volume before computing gradients
+    """
+    from scipy.ndimage.filters import gaussian_filter, sobel
+    if sigma is not None:
+        im = gaussian_filter(im, sigma)
+    grads = [sobel(im, axis=dim, mode='constant') / 8.0 for dim in xrange(im.ndim)]
+    return grads
+
+
+def imageJacobian(vol, tfm, grid=None, sigma=None, normalize=True, border=1, order=1):
+    """Compute Jacobian of volume w.r.t. transformation parameters
+
+    Args:
+        vol: volume
+        tfm: Transform object
+        sigma: smoothing bandwidth for gradients (None for no smoothing)
+        normalize: Whether to normalize images before aligning.
+        border: Number or tuple of border sizes to zero after transforming.
+        order: interpolation order used by map_coordinates
+    Returns:
+        tvol : array
+            transformed volume
+        jacobianVols : list of arrays
+            list of volume Jacobians, one for each parameter of the transformation
+    """
+    if grid is None:
+        from thunder.imgprocessing.transformation import GridTransformer
+        grid = GridTransformer(vol.shape)
+    grads = imageGradients(vol, sigma)
+    tvol = zeroBorder(tfm.apply(vol, grid, order=order))
+    grads = [zeroBorder(tfm.apply(grad, grid, order=order)) for grad in grads]
+    if normalize:
+        norm = np.linalg.norm(tvol.ravel())
+        if norm == 0.0:
+            raise ValueError('Transform yields volume of zeroes.')
+        # Update gradients to reflect normalization
+        grads = [grad / norm - (grad * tvol).sum() / (norm**3) * tvol for grad in grads]
+        tvol /= norm
+    jacobianVols = tfm.jacobian(grads, grid.homo_points)
+    return tvol, jacobianVols