Initial commit

.gitignore

@@ -0,0 +1,2 @@
+__pycache__
+*.pyc

README.md

@@ -0,0 +1,18 @@
+# espirit-python
+
+ESPIRiT implemented in Python.
+
+## Prerequisites
+
+ESPIRiT itself requires ```numpy```. ```matplotlib``` is needed to display images.
+
+## Usage
+
+Please see ```example.py``` for a usage example.
+
+## References
+
+- Uecker, M., Lai, P., Murphy, M. J., Virtue, P., Elad, M., Pauly, J. M., ... & Lustig, M. (2014). ESPIRiT—an
+  eigenvalue approach to autocalibrating parallel MRI: where SENSE meets GRAPPA. Magnetic resonance in medicine,
+  71(3), 990-1001.
+- Data from [mridata.org](mridata.org)

cfl.py

@@ -0,0 +1,41 @@
+# Copyright 2013-2015. The Regents of the University of California.
+# All rights reserved. Use of this source code is governed by 
+# a BSD-style license which can be found in the LICENSE file.
+#
+# Authors: 
+# 2013 Martin Uecker <[email protected]>
+# 2015 Jonathan Tamir <[email protected]>
+
+
+import numpy as np
+
+def readcfl(name):
+    # get dims from .hdr
+    h = open(name + ".hdr", "r")
+    h.readline() # skip
+    l = h.readline()
+    h.close()
+    dims = [int(i) for i in l.split( )]
+
+    # remove singleton dimensions from the end
+    n = np.prod(dims)
+    dims_prod = np.cumprod(dims)
+    dims = dims[:np.searchsorted(dims_prod, n)+1]
+
+    # load data and reshape into dims
+    d = open(name + ".cfl", "r")
+    a = np.fromfile(d, dtype=np.complex64, count=n);
+    d.close()
+    return a.reshape(dims, order='F') # column-major
+
+
+def writecfl(name, array):
+    h = open(name + ".hdr", "w")
+    h.write('# Dimensions\n')
+    for i in (array.shape):
+            h.write("%d " % i)
+    h.write('\n')
+    h.close()
+    d = open(name + ".cfl", "w")
+    array.T.astype(np.complex64).tofile(d) # tranpose for column-major order
+    d.close()

data/knee.hdr

@@ -0,0 +1,2 @@
+# Dimensions
+1 320 256 8 1 

espirit.py

@@ -0,0 +1,117 @@
+import numpy as np
+
+fft  = lambda x, ax : np.fft.fftshift(np.fft.fftn(np.fft.ifftshift(x, axes=ax), axes=ax, norm='ortho'), axes=ax) 
+ifft = lambda X, ax : np.fft.fftshift(np.fft.ifftn(np.fft.ifftshift(X, axes=ax), axes=ax, norm='ortho'), axes=ax) 
+
+def espirit(X, k, r, t, c):
+    """
+    Derives the ESPIRiT operator.
+
+    Arguments:
+      X: Multi channel k-space data. Expected dimensions are (sx, sy, sz, nc), where (sx, sy, sz) are volumetric 
+         dimensions and (nc) is the channel dimension.
+      k: Parameter that determines the k-space kernel size. If X has dimensions (1, 256, 256, 8), then the kernel 
+         will have dimensions (1, k, k, 8)
+      r: Parameter that determines the calibration region size. If X has dimensions (1, 256, 256, 8), then the 
+         calibration region will have dimensions (1, r, r, 8)
+      t: Parameter that determines the rank of the auto-calibration matrix (A). Singular values below t times the
+         largest singular value are set to zero.
+      c: Crop threshold that determines eigenvalues "=1".
+    Returns:
+      maps: This is the ESPIRiT operator. It will have dimensions (sx, sy, sz, nc, nc) with (sx, sy, sz, :, idx)
+            being the idx'th set of ESPIRiT maps.
+    """
+
+    sx = np.shape(X)[0]
+    sy = np.shape(X)[1]
+    sz = np.shape(X)[2]
+    nc = np.shape(X)[3]
+
+    sxt = (sx//2-r//2, sx//2+r//2) if (sx > 1) else (0, 1)
+    syt = (sy//2-r//2, sy//2+r//2) if (sy > 1) else (0, 1)
+    szt = (sz//2-r//2, sz//2+r//2) if (sz > 1) else (0, 1)
+
+    # Extract calibration region.    
+    C = X[sxt[0]:sxt[1], syt[0]:syt[1], szt[0]:szt[1], :].astype(np.complex64)
+
+    # Construct Hankel matrix.
+    p = (sx > 1) + (sy > 1) + (sz > 1)
+    A = np.zeros([(r-k+1)**p, k**p * nc]).astype(np.complex64)
+
+    idx = 0
+    for xdx in range(max(1, C.shape[0] - k + 1)):
+      for ydx in range(max(1, C.shape[1] - k + 1)):
+        for zdx in range(max(1, C.shape[2] - k + 1)):
+          # numpy handles when the indices are too big
+          block = C[xdx:xdx+k, ydx:ydx+k, zdx:zdx+k, :].astype(np.complex64) 
+          A[idx, :] = block.flatten()
+          idx = idx + 1
+
+    # Take the Singular Value Decomposition.
+    U, S, VH = np.linalg.svd(A, full_matrices=True)
+    V = VH.conj().T
+
+    # Select kernels.
+    n = np.sum(S >= t * S[0])
+    V = V[:, 0:n]
+
+    kxt = (sx//2-k//2, sx//2+k//2) if (sx > 1) else (0, 1)
+    kyt = (sy//2-k//2, sy//2+k//2) if (sy > 1) else (0, 1)
+    kzt = (sz//2-k//2, sz//2+k//2) if (sz > 1) else (0, 1)
+
+    # Reshape into k-space kernel, flips it and takes the conjugate
+    kernels = np.zeros(np.append(np.shape(X), n)).astype(np.complex64)
+    kerdims = [(sx > 1) * k + (sx == 1) * 1, (sy > 1) * k + (sy == 1) * 1, (sz > 1) * k + (sz == 1) * 1, nc]
+    for idx in range(n):
+        kernels[kxt[0]:kxt[1],kyt[0]:kyt[1],kzt[0]:kzt[1], :, idx] = np.reshape(V[:, idx], kerdims)
+
+    # Take the iucfft
+    axes = (0, 1, 2)
+    kerimgs = np.zeros(np.append(np.shape(X), n)).astype(np.complex64)
+    for idx in range(n):
+        for jdx in range(nc):
+            ker = kernels[::-1, ::-1, ::-1, jdx, idx].conj()
+            kerimgs[:,:,:,jdx,idx] = fft(ker, axes) * np.sqrt(sx * sy * sz)/np.sqrt(k**p)
+
+    # Take the point-wise eigenvalue decomposition and keep eigenvalues greater than c
+    maps = np.zeros(np.append(np.shape(X), nc)).astype(np.complex64)
+    for idx in range(0, sx):
+        for jdx in range(0, sy):
+            for kdx in range(0, sz):
+
+                Gq = kerimgs[idx,jdx,kdx,:,:]
+
+                u, s, vh = np.linalg.svd(Gq, full_matrices=True)
+                for ldx in range(0, nc):
+                    if (s[ldx]**2 > c):
+                        maps[idx, jdx, kdx, :, ldx] = u[:, ldx]
+
+    return maps
+
+def espirit_proj(x, esp):
+    """
+    Construct the projection of multi-channel image x onto the range of the ESPIRiT operator. Returns the inner
+    product, complete projection and the null projection.
+
+    Arguments:
+      x: Multi channel image data. Expected dimensions are (sx, sy, sz, nc), where (sx, sy, sz) are volumetric 
+         dimensions and (nc) is the channel dimension.
+      esp: ESPIRiT operator as returned by function: espirit
+
+    Returns:
+      ip: This is the inner product result, or the image information in the ESPIRiT subspace.
+      proj: This is the resulting projection. If the ESPIRiT operator is E, then proj = E E^H x, where H is 
+            the hermitian.
+      null: This is the null projection, which is equal to x - proj.
+    """
+    ip = np.zeros(x.shape).astype(np.complex64)
+    proj = np.zeros(x.shape).astype(np.complex64)
+    for qdx in range(0, esp.shape[4]):
+        for pdx in range(0, esp.shape[3]):
+            ip[:, :, :, qdx] = ip[:, :, :, qdx] + x[:, :, :, pdx] * esp[:, :, :, pdx, qdx].conj()
+
+    for qdx in range(0, esp.shape[4]):
+        for pdx in range(0, esp.shape[3]):
+          proj[:, :, :, pdx] = proj[:, :, :, pdx] + ip[:, :, :, qdx] * esp[:, :, :, pdx, qdx]
+
+    return (ip, proj, x - proj)

example.py

@@ -0,0 +1,58 @@
+import cfl
+from espirit import espirit, espirit_proj, ifft
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+# Load data
+X = cfl.readcfl('data/knee')
+x = ifft(X, (0, 1, 2))
+
+# Derive ESPIRiT operator
+esp = espirit(X, 6, 24, 0.01, 0.9925)
+# Do projections
+ip, proj, null = espirit_proj(x, esp)
+
+# Figure code
+
+esp  = np.squeeze(esp)
+x    = np.squeeze(x)
+ip   = np.squeeze(ip)
+proj = np.squeeze(proj)
+null = np.squeeze(null)
+
+print("Close figures to continue execution...")
+
+# Display ESPIRiT operator
+for idx in range(8):
+    for jdx in range(8):
+        plt.subplot(8, 8, (idx * 8 + jdx) + 1)
+        plt.imshow(np.abs(esp[:,:,idx,jdx]), cmap='gray')
+        plt.axis('off')
+plt.show()
+
+dspx = np.power(np.abs(np.concatenate((x[:, :, 0], x[:, :, 1], x[:, :, 2], x[:, :, 3], x[:, :, 4], x[:, :, 5], x[:, :, 6], x[:, :, 7]), 0)), 1/3)
+dspip = np.power(np.abs(np.concatenate((ip[:, :, 0], ip[:, :, 1], ip[:, :, 2], ip[:, :, 3], ip[:, :, 4], ip[:, :, 5], ip[:, :, 6], ip[:, :, 7]), 0)), 1/3)
+dspproj = np.power(np.abs(np.concatenate((proj[:, :, 0], proj[:, :, 1], proj[:, :, 2], proj[:, :, 3], proj[:, :, 4], proj[:, :, 5], proj[:, :, 6], proj[:, :, 7]), 0)), 1/3)
+dspnull = np.power(np.abs(np.concatenate((null[:, :, 0], null[:, :, 1], null[:, :, 2], null[:, :, 3], null[:, :, 4], null[:, :, 5], null[:, :, 6], null[:, :, 7]), 0)), 1/3)
+
+print("NOTE: Contrast has been changed")
+
+# Display ESPIRiT projection results 
+plt.subplot(1, 4, 1)
+plt.imshow(dspx, cmap='gray')
+plt.title('Data')
+plt.axis('off')
+plt.subplot(1, 4, 2)
+plt.imshow(dspip, cmap='gray')
+plt.title('Inner product')
+plt.axis('off')
+plt.subplot(1, 4, 3)
+plt.imshow(dspproj, cmap='gray')
+plt.title('Projection')
+plt.axis('off')
+plt.subplot(1, 4, 4)
+plt.imshow(dspnull, cmap='gray')
+plt.title('Null Projection')
+plt.axis('off')
+plt.show()

index.html

@@ -0,0 +1,12 @@
+<center><h1 id="espirit-python">espirit-python</h1></center>
+<h2 id="description">Description</h2>
+<p>ESPIRiT implemented in Python.</p>
+<h2 id="prerequisites">Prerequisites</h2>
+<p>ESPIRiT itself requires <code>numpy</code>. <code>matplotlib</code> is needed to display images.</p>
+<h2 id="usage">Usage</h2>
+<p>Please see <code>example.py</code> for a usage example.</p>
+<h2 id="references">References</h2>
+<ul>
+<li>Uecker, M., Lai, P., Murphy, M. J., Virtue, P., Elad, M., Pauly, J. M., ... &amp; Lustig, M. (2014). ESPIRiT—an eigenvalue approach to autocalibrating parallel MRI: where SENSE meets GRAPPA. Magnetic resonance in medicine, 71(3), 990-1001.</li>
+<li>Data from <a href="http://mridata.org">mridata.org</a></li>
+</ul>