Updated code with working versions, cleaned up

create_mask.py

@@ -0,0 +1,37 @@
+import numpy as np
+from scipy.ndimage import measurements as m
+from scipy.ndimage import generate_binary_structure, binary_fill_holes
+from skimage.measure import regionprops
+
+from functions import imscale, smv_filter
+
+def create_mask(mag, level):
+    """ Creates mask using thresholding.
+    level: typically around 0.03 - 0.1
+    """
+    mag = imscale(np.abs(mag))
+    mask = mag > level
+    mask = polish_mask(mask)
+    return mask
+
+def polish_mask(mask):
+    """Chooses the largest connected structure and fills holes."""
+    labels, num_labels = m.label(mask, generate_binary_structure(3,3))
+    props = regionprops(labels, cache = False)
+    biggest = np.argmax([r.area for r in props]) + 1
+    mask = labels == biggest
+    mask = binary_fill_holes(mask)
+    return mask
+
+def create_mask_qsm(mag, voxel_size, level = 0.04):
+    """Creates mask for qsm using thresholding and spherical-mean-value
+    filtering. Use on echo closest to 30 ms.
+    """
+    mask = create_mask(mag, level)
+    mask = smv_filter(mask, 2, voxel_size)
+    mask = polish_mask(mask > 0.99)
+    mask = smv_filter(mask, 3, voxel_size)
+    mask = mask > 0.02
+    mask = smv_filter(mask, 1, voxel_size)
+    mask = mask > 0.97
+    return mask

functions.py

@@ -0,0 +1,81 @@
+import numpy as np
+from scipy.interpolate import RegularGridInterpolator as rgi
+from scipy.ndimage.interpolation import rotate
+
+def ifftnc(arr):
+    scaling = np.sqrt(arr.size)
+    return np.fft.fftshift(np.fft.ifftn(np.fft.fftshift(arr))) * scaling
+
+def fftnc(arr):
+    scaling = np.sqrt(arr.size)
+    return np.fft.fftshift(np.fft.fftn(np.fft.fftshift(arr))) / scaling
+
+def interp3(xx, yy, zz, V, points, shape):
+    interp = rgi((xx, yy, zz), V, bounds_error = False, fill_value = 0)
+    values = interp(points)
+    values = values.reshape((int(shape[1]),int(shape[0]),int(shape[2])))
+    values = rotate(values, -90)
+    values = np.fliplr(values)
+    return values
+
+def imscale(image):
+    max_v = np.max(image)
+    min_v = np.min(image)
+    return (image - min_v) / (max_v - min_v)
+
+def ball_3d(radius, voxel_size):
+    xx, yy, zz = np.meshgrid(np.arange(-radius, radius + 1),
+                             np.arange(-radius, radius + 1),
+                             np.arange(-radius, radius + 1))
+    xx = xx*voxel_size[1]
+    yy = yy*voxel_size[0]
+    zz = zz*voxel_size[2]
+    sq_dist = np.square(xx) + np.square(yy) + np.square(zz)
+    index = np.logical_and(sq_dist <= (radius**2 + radius / 3), sq_dist >= 0.1)
+    return index
+
+def total_field(shape, radius, voxel_size):
+    kernel_3d = ball_3d(radius, voxel_size)
+    kernel_3d = kernel_3d / np.sum(kernel_3d)
+    f0 = np.zeros(shape)
+    f0[f0.shape[0]//2 - radius : f0.shape[0]//2 + radius + 1,
+       f0.shape[1]//2 - radius : f0.shape[1]//2 + radius + 1,
+       f0.shape[2]//2 - radius : f0.shape[2]//2 + radius + 1] = kernel_3d
+    f0 = fftnc(f0 * np.sqrt(f0.size))
+    return f0
+
+def smv_filter(phase, radius, voxel_size):
+    f0 = total_field(phase.shape, radius, voxel_size)
+    out_phase = ifftnc(f0 * fftnc(phase))
+    return np.real(out_phase)
+
+def pad(arr, pad_size):
+    pad_size = tuple(pad_size)
+    while len(pad_size) < len(arr.shape):
+        pad_size = pad_size + (0,)
+    pad_size = tuple([(size, size) for size in pad_size])
+    return np.pad(arr, pad_size, 'constant')
+
+def left_pad(arr, axis):
+    pad_size = tuple([(int(i == axis), 0) for i, a in enumerate(arr.shape)])
+    return np.pad(arr, pad_size, 'constant')
+
+def unpad(arr, pad_size):
+    pad_size = tuple(pad_size)
+    while len(pad_size) < len(arr.shape):
+        pad_size = pad_size + (0,)
+    pad_size = tuple([(size, size) for size in pad_size])
+    slices = []
+    for p in pad_size:
+        if p[1] == 0:
+            slices.append(slice(p[0], None))
+        else:
+            slices.append(slice(p[0], -p[1]))
+    return arr[slices]
+
+def bbox_slice(bbox):
+    front, back = bbox[:len(bbox)//2], bbox[len(bbox)//2:]
+    slices = []
+    for i in range(len(front)):
+        slices.append(slice(front[i], back[i]))
+    return slices

laplacian_unwrap.py

@@ -0,0 +1,45 @@
+import numpy as np
+
+from functions import ifftnc, fftnc, pad, left_pad, unpad
+
+def laplacian_unwrap(phase, voxel_size, pad_size = (12,12,12)):
+    """
+    Implements Laplacian-based unwrapping 
+    (http://onlinelibrary.wiley.com/doi/10.1002/nbm.3056/full)
+    
+    Steven Cao, Hongjiang Wei, Wei Li, Chunlei Liu
+    University of California, Berkeley
+    """
+    # add 4th singleton dimension if phase is 3d
+    if len(phase.shape) == 3:
+        phase = phase[:,:,:,None]
+    assert len(phase.shape) == 4, 'Phase must be 3d or 4d'
+    phase = pad(phase, pad_size)
+
+    if phase.shape[2] % 2 == 1:
+        phase = left_pad(phase, 2)
+
+    yy, xx, zz = np.meshgrid(np.arange(1, phase.shape[1]+1),
+                             np.arange(1, phase.shape[0]+1),
+                             np.arange(1, phase.shape[2]+1))
+    field_of_view = np.multiply(voxel_size, phase.shape[0:3])
+    xx, yy, zz = ((xx - phase.shape[0]/2 - 1) / field_of_view[0],
+                  (yy - phase.shape[1]/2 - 1) / field_of_view[1],
+                  (zz - phase.shape[2]/2 - 1) / field_of_view[2])
+    k2 = np.square(xx) + np.square(yy) + np.square(zz)
+    xx, yy, zz, field_of_view = None, None, None, None
+
+    laplacian = np.zeros(phase.shape, dtype = np.complex)
+    phi = np.zeros(phase.shape, dtype = np.complex)
+    for i in range(0, phase.shape[3]):
+        laplacian[:,:,:,i] = (np.cos(phase[:,:,:,i]) * 
+                 ifftnc(k2 * fftnc(np.sin(phase[:,:,:,i]))) - 
+                 np.sin(phase[:,:,:,i]) * 
+                 ifftnc(k2 * fftnc(np.cos(phase[:,:,:,i]))))
+        phi[:,:,:,i] = fftnc(laplacian[:,:,:,i]) / k2
+        phi[phase.shape[0]//2,phase.shape[1]//2,phase.shape[2]//2,i] = 0
+        phi[:,:,:,i] = ifftnc(phi[:,:,:,i])
+    laplacian = unpad(laplacian, pad_size)
+    phi = unpad(phi, pad_size)
+    phi = np.real(phi)
+    return phi, laplacian

qsm_star.py

@@ -0,0 +1,144 @@
+import numpy as np
+from functions import ifftnc, fftnc, pad, unpad
+
+def qsm_star(phase, mask, voxel_size, TE, pad_size = (8,8,20), 
+             B0 = 3, B0_dir = (0,0,1), tau = 1e-6, d2_thresh = .065):
+    d2 = calc_d2_matrix(phase.shape, voxel_size, B0_dir)
+    # sample = np.abs(d2) > d2_thresh
+    d2 = pad(d2, pad_size)
+    phase, mask = pad(phase, pad_size), pad(mask, pad_size)
+    A = lambda x: ifftnc(d2 * fftnc(x * mask))
+    AT = A
+    susc = sparsa(phase, np.zeros(phase.shape), A, AT, tau)
+    susc *= mask
+    susc = unpad(susc, pad_size)
+    susc *= scaling_factor(B0, TE)
+    return np.real(susc)
+
+def scaling_factor(B0, TE_ms):
+    gamma = 42.575 * 2 * np.pi
+    TE = TE_ms * 1e-3
+    return 1 / (gamma * B0 * TE)
+
+def calc_d2_matrix(shape, voxel_size, B0_dir):
+    shape = np.array(shape)
+    voxel_size = np.array(voxel_size)
+    B0_dir = np.array(B0_dir)
+    field_of_view = shape * voxel_size
+    ry, rx, rz = np.meshgrid(np.arange(-shape[1]//2, shape[1]//2),
+                             np.arange(-shape[0]//2, shape[0]//2),
+                             np.arange(-shape[2]//2, shape[2]//2))
+    rx, ry, rz = rx/field_of_view[0], ry/field_of_view[1], rz/field_of_view[2]
+    sq_dist = rx**2 + ry**2 + rz**2
+    sq_dist[sq_dist==0] = 1e-6
+    d2 = ((B0_dir[0]*rx + B0_dir[1]*ry + B0_dir[2]*rz)**2)/sq_dist
+    d2 = (1/3 - d2)
+    return d2
+
+def soft(x, tau):
+    if np.sum(np.abs(tau)) == 0:
+        return x
+    else:
+        y = np.abs(x) - tau
+        y[y < 0] = 0
+        y = y / (y + tau) * x
+        return y
+
+def sparsa(y, x, A, AT, tau, verbose = False, min_iter = 5, alpha = 1):
+    """
+    SpaRSA solver modified from the Matlab version by Mario Figueiredo, 
+    Robert Nowak, and Stephen Wright (http://www.lx.it.pt/~mtf/SpaRSA/)
+    
+    Inputs:
+    y: 1D vector or 2D array (image) of observations
+    x: initial guess for the solution x
+    A: matrix to be applied to x (in function form)
+    AT: transpose of A (in function form)
+    tau: regularization parameter (scalar)
+    
+    Compared to the Matlab code, it uses the following options:
+    psi = soft
+    phi = l1 norm
+    StopCriterion = 0 - algorithm stops when the relative 
+                        change in the number of non-zero 
+                        components of the estimate falls 
+                        below 'tolerance'
+    debias = 0 (no)
+    BB_variant = 1 - standard BB choice  s'r/s's
+    monotone = 0 (don't enforce monotonic decrease of f)
+    safeguard = 0 (don't enforce a "sufficient decrease" over the largest
+                   objective value of the past M iterations.)
+    continuation = 0 (don't do the continuation procedure)
+    
+    Steven Cao
+    University of California, Berkeley
+    
+    ==========================================================================    
+    SpaRSA version 2.0, December 31, 2007
+     
+    This function solves the convex problem 
+    
+    arg min_x = 0.5*|| y - A x ||_2^2 + tau phi(x)
+    
+    d'*A'*(A x - y) + 0.5*alpha*d'*d 
+    
+    where alpha is obtained from a BB formula. In a monotone variant, alpha is
+    increased until we see a decreasein the original objective function over
+    this step. 
+    
+    -----------------------------------------------------------------------
+    Copyright (2007): Mario Figueiredo, Robert Nowak, Stephen Wright
+    ----------------------------------------------------------------------
+    """
+    tolerance = 0.01 #min amount of changes needed to keep going
+    max_iter = 100
+    min_alpha = 1e-30
+    max_alpha = 1e30
+
+    psi_function = soft
+    phi_function = lambda x: np.sum(np.abs(x)) #phi is l1
+
+    nonzero_x = x != 0
+
+    resid = A(x) - y
+    gradient = AT(resid)
+    alpha = 1
+    f = (0.5 * np.vdot(resid.flatten(), resid.flatten())
+        + tau * phi_function(x))
+
+    iterations = 0
+    while True:
+        gradient = AT(resid)
+        prev_x = x
+        prev_resid = resid
+
+        x = psi_function(prev_x - gradient * (1/alpha),
+                         tau/alpha)
+        dx = x - prev_x
+        Adx = A(dx)
+        resid = prev_resid + Adx
+        f = (0.5 * np.vdot(resid.flatten(), resid.flatten())
+            + tau * phi_function(x))
+
+        dd = np.vdot(dx.flatten(), dx.flatten())
+        dGd = np.vdot(Adx.flatten(), Adx.flatten())
+        alpha = dGd / (np.finfo(np.double).tiny + dd)
+        alpha = min(max_alpha, max(min_alpha, alpha))
+        iterations += 1
+
+        # compute the stopping criterion based on the change
+        # in the number of non-zero components of the estimate
+        prev_nonzero_x = nonzero_x
+        nonzero_x = x != 0
+        changes = np.sum(nonzero_x != prev_nonzero_x) / np.sum(nonzero_x)
+        if verbose:
+            print('-----------------------',
+                  '\nObjective f:', f,
+                  '\nAlpha:', alpha,
+                  '\nPercent change in x:', changes,
+                  '\nChange tolerance:', tolerance, 
+                  '\nIterations:', iterations)
+        if (changes <= tolerance and
+           iterations >= min_iter and
+           iterations < max_iter):
+           return x