Different basic undersampling methods.

.gitignore

@@ -0,0 +1,8 @@
+venv/
+*.cfl
+*.hdr
+fig/
+.ipynb_checkpoints/
+data/
+__pycache__/
+*.mat

README.md

@@ -0,0 +1,3 @@
+# ESMRMB Educational
+
+Testing different reconstruction methods.

cfl.py

@@ -0,0 +1,114 @@
+# Copyright 2013-2015. The Regents of the University of California.
+# Copyright 2021. Uecker Lab. University Center Göttingen.
+# 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 <uecker@eecs.berkeley.edu>
+# 2015 Jonathan Tamir <jtamir@eecs.berkeley.edu>
+
+from __future__ import print_function
+from __future__ import with_statement
+
+import numpy as np
+import mmap
+import os
+
+
+def readcfl(name):
+    # get dims from .hdr
+    with open(name + ".hdr", "rt") as h:
+        h.readline() # skip
+        l = h.readline()
+    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
+    with open(name + ".cfl", "rb") as d:
+        a = np.fromfile(d, dtype=np.complex64, count=n);
+    return a.reshape(dims, order='F') # column-major
+
+def readmulticfl(name):
+    # get dims from .hdr
+    with open(name + ".hdr", "rt") as h:
+        lines = h.read().splitlines()
+
+    index_dim = 1 + lines.index('# Dimensions')
+    total_size = int(lines[index_dim])
+    index_sizes = 1 + lines.index('# SizesDimensions')
+    sizes = [int(i) for i in lines[index_sizes].split()]
+    index_dims = 1 + lines.index('# MultiDimensions')
+
+    with open(name + ".cfl", "rb") as d:
+        a = np.fromfile(d, dtype=np.complex64, count=total_size)
+
+    offset = 0
+    result = []
+    for i in range(len(sizes)):
+        dims = ([int(i) for i in lines[index_dims + i].split()])
+        n = np.prod(dims)
+        result.append(a[offset:offset+n].reshape(dims, order='F'))
+        offset += n
+
+    if total_size != offset:
+        print("Error")
+
+    return result
+
+
+def writecfl(name, array):
+    with open(name + ".hdr", "wt") as h:
+        h.write('# Dimensions\n')
+        for i in (array.shape):
+                h.write("%d " % i)
+        h.write('\n')
+
+    size = np.prod(array.shape) * np.dtype(np.complex64).itemsize
+
+    with open(name + ".cfl", "a+b") as d:
+        os.ftruncate(d.fileno(), size)
+        mm = mmap.mmap(d.fileno(), size, flags=mmap.MAP_SHARED, prot=mmap.PROT_WRITE)
+        if array.dtype != np.complex64:
+            array = array.astype(np.complex64)
+        mm.write(np.ascontiguousarray(array.T))
+        mm.close()
+        #with mmap.mmap(d.fileno(), size, flags=mmap.MAP_SHARED, prot=mmap.PROT_WRITE) as mm:
+        #    mm.write(array.astype(np.complex64).tobytes(order='F'))
+
+def writemulticfl(name, arrays):
+    size = 0
+    dims = []
+
+    for array in arrays:
+        size += array.size
+        dims.append(array.shape)
+
+    with open(name + ".hdr", "wt") as h:
+        h.write('# Dimensions\n')
+        h.write("%d\n" % size)
+
+        h.write('# SizesDimensions\n')
+        for dim in dims:
+            h.write("%d " % len(dim))
+        h.write('\n')
+
+        h.write('# MultiDimensions\n')
+        for dim in dims:
+            for i in dim:
+                h.write("%d " % i)
+            h.write('\n')
+            
+    size = size * np.dtype(np.complex64).itemsize
+
+    with open(name + ".cfl", "a+b") as d:
+        os.ftruncate(d.fileno(), size)
+        mm = mmap.mmap(d.fileno(), size, flags=mmap.MAP_SHARED, prot=mmap.PROT_WRITE)
+        for array in arrays:
+            if array.dtype != np.complex64:
+                array = array.astype(np.complex64)
+            mm.write(np.ascontiguousarray(array.T))
+        mm.close()

plotBoMap.py

@@ -0,0 +1,48 @@
+import scipy.io as sio
+import matplotlib.pyplot as plt
+import numpy as np
+import matplotlib.gridspec as gridspec
+import matplotlib.cm as cm
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+
+# Bo Map
+mat_data_0=sio.loadmat('bo.mat')
+tmap = mat_data_0['boMap']
+print(np.min(tmap))
+print(np.max(tmap))
+vmin = -0.3
+vmax = 0.3
+
+def multiSlicePlot(data,nRow, nCol, nSlrep, vmin, vmax, sliceInit = 0):
+    """""
+    data = imageMatrix3D [sl,ph,rd]
+    nRow, nCol dimensions of the multislicePlot
+    nSlRep =  number of slices to represent
+    sliceInit = slice in wihich we start representing
+    """""
+    images = []
+    fig = plt.figure(figsize=(nRow, nCol), dpi=500)
+    gs1 = gridspec.GridSpec(nRow, nCol)
+    gs1.update(wspace=0.020, hspace=0.020)  # set the spacing between axes.
+    for i in range(nSlrep):
+        ii = i + sliceInit
+        # print(ii)
+        ax1 = plt.subplot(gs1[i])
+        # plt.axis('off')
+        ax1.set_xticklabels([])
+        ax1.set_yticklabels([])
+        ax1.set_aspect('equal')
+        dataAux = data[int(ii), :, :]
+        imgPlot = ax1.imshow(dataAux,vmin=vmin, vmax=vmax)
+        images.append(imgPlot)
+        ax1.axis('off')
+    return fig
+
+fig = multiSlicePlot(tmap*1e3,3,6,16,vmin,vmax,0)
+
+cbar_ax = fig.add_axes([0.92, 0.15, 0.01, 0.7]) 
+norm = plt.Normalize(vmin=vmin, vmax=vmax)
+cbar = plt.colorbar(cm.ScalarMappable(norm=norm, cmap=cm.viridis), cax=cbar_ax)
+cbar.set_label('Bo [mT]', fontsize=2)
+cbar.ax.tick_params(labelsize=2)
+plt.show()

plotComparativeTimes.py

@@ -0,0 +1,39 @@
+import scipy.io as sio
+import matplotlib.pyplot as plt
+import numpy as np
+import matplotlib.gridspec as gridspec
+import matplotlib.cm as cm
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+
+imagesPath=['T2_CPMG','T2_APCPMG', 'T2_APCP', 'T2_CP', 'T1']
+sliceRep = 12
+vmin = 0
+vmax = 200
+
+nMaps = len(imagesPath)
+fig = plt.figure(figsize=(1, 5), dpi=500)
+gs1 = gridspec.GridSpec(1, nMaps)
+gs1.update(wspace=0.020, hspace=0.020)  
+for i in range(nMaps):
+    if 'T1' in imagesPath[i]:
+        mat_data = sio.loadmat(imagesPath[i]+'.mat')
+        data = mat_data['t1map']
+        print(data.shape)
+    else:
+        mat_data = sio.loadmat(imagesPath[i]+'.mat')
+        data = mat_data['t2map']
+        print(data.shape)
+    ax1 = plt.subplot(gs1[i])
+    ax1.set_xticklabels([])
+    ax1.set_yticklabels([])
+    ax1.set_aspect('equal')
+    dataAux = data[sliceRep, :, :]
+    imgPlot = ax1.imshow(abs(dataAux),vmin=vmin, vmax=vmax)
+    ax1.axis('off')
+    ax1.set_title(imagesPath[i], fontsize=4)
+cbar_ax = fig.add_axes([0.92, 0.35, 0.01, 0.30]) 
+norm = plt.Normalize(vmin=vmin, vmax=vmax)
+cbar = plt.colorbar(cm.ScalarMappable(norm=norm, cmap=cm.viridis), cax=cbar_ax)
+cbar.set_label('Time [ms]', fontsize=2)
+cbar.ax.tick_params(labelsize=2)
+plt.show()

plotPDMap.py

@@ -0,0 +1,46 @@
+import scipy.io as sio
+import matplotlib.pyplot as plt
+import numpy as np
+import matplotlib.gridspec as gridspec
+import matplotlib.cm as cm
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+
+# PD
+mat_data_0=sio.loadmat('PD.mat')
+img = mat_data_0['image3D']
+vmin = np.min(np.abs(img))
+vmax = np.max(np.abs(img))
+
+def multiSlicePlot(data,nRow, nCol, nSlrep, vmin, vmax, sliceInit = 0):
+    """""
+    data = imageMatrix3D [sl,ph,rd]
+    nRow, nCol dimensions of the multislicePlot
+    nSlRep =  number of slices to represent
+    sliceInit = slice in wihich we start representing
+    """""
+    images = []
+    fig = plt.figure(figsize=(nRow, nCol), dpi=500)
+    gs1 = gridspec.GridSpec(nRow, nCol)
+    gs1.update(wspace=0.020, hspace=0.020)  # set the spacing between axes.
+    for i in range(nSlrep):
+        ii = i + sliceInit
+        # print(ii)
+        ax1 = plt.subplot(gs1[i])
+        # plt.axis('off')
+        ax1.set_xticklabels([])
+        ax1.set_yticklabels([])
+        ax1.set_aspect('equal')
+        dataAux = data[int(ii), :, :]
+        imgPlot = ax1.imshow(dataAux,vmin=vmin, vmax=vmax)
+        images.append(imgPlot)
+        ax1.axis('off')
+    return fig
+
+fig = multiSlicePlot(np.abs(img),3,6,16,vmin,vmax,0)
+
+cbar_ax = fig.add_axes([0.92, 0.15, 0.01, 0.7]) 
+norm = plt.Normalize(vmin=vmin, vmax=vmax)
+cbar = plt.colorbar(cm.ScalarMappable(norm=norm, cmap=cm.viridis), cax=cbar_ax)
+cbar.set_label('PD [a.u.]', fontsize=2)
+cbar.ax.tick_params(labelsize=2)
+plt.show()

plotT1Map.py

@@ -0,0 +1,12 @@
+import scipy.io as sio
+import matplotlib.pyplot as plt
+import numpy as np
+from bart import bart
+
+# T1 dict
+T1_dict = sio.loadmat('T1.mat')
+
+rawName = 'T2_CPMG.mat'
+mat_data_0=sio.loadmat(rawName)
+
+print(mat_data_0.keys())

plotTimesMap.py

@@ -0,0 +1,52 @@
+import scipy.io as sio
+import matplotlib.pyplot as plt
+import numpy as np
+import matplotlib.gridspec as gridspec
+import matplotlib.cm as cm
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+
+# Map
+rawName = 'T2_CPMG.mat'
+mat_data_0=sio.loadmat(rawName)
+if 'T1' in rawName:
+    tmap = mat_data_0['t1map']
+else:
+    tmap = mat_data_0['t2map']
+print(np.min(tmap))
+print(np.max(tmap))
+vmin = 0
+vmax = 200
+
+def multiSlicePlot(data,nRow, nCol, nSlrep, vmin, vmax, sliceInit = 0):
+    """""
+    data = imageMatrix3D [sl,ph,rd]
+    nRow, nCol dimensions of the multislicePlot
+    nSlRep =  number of slices to represent
+    sliceInit = slice in wihich we start representing
+    """""
+    images = []
+    fig = plt.figure(figsize=(nRow, nCol), dpi=500)
+    gs1 = gridspec.GridSpec(nRow, nCol)
+    gs1.update(wspace=0.020, hspace=0.020)  # set the spacing between axes.
+    for i in range(nSlrep):
+        ii = i + sliceInit
+        # print(ii)
+        ax1 = plt.subplot(gs1[i])
+        # plt.axis('off')
+        ax1.set_xticklabels([])
+        ax1.set_yticklabels([])
+        ax1.set_aspect('equal')
+        dataAux = data[int(ii), :, :]
+        imgPlot = ax1.imshow(dataAux,vmin=vmin, vmax=vmax)
+        images.append(imgPlot)
+        ax1.axis('off')
+    return fig
+
+fig = multiSlicePlot(tmap,3,6,16,vmin,vmax,0)
+
+cbar_ax = fig.add_axes([0.92, 0.15, 0.01, 0.7]) 
+norm = plt.Normalize(vmin=vmin, vmax=vmax)
+cbar = plt.colorbar(cm.ScalarMappable(norm=norm, cmap=cm.viridis), cax=cbar_ax)
+cbar.set_label(rawName+ '[ms]', fontsize=2)
+cbar.ax.tick_params(labelsize=2)
+plt.show()

rawDatasExplanation.txt

@@ -0,0 +1,36 @@
+
+# T2 Maps
+t2maps_dict = {
+    'kSpaces3D': sampled, # kSpace [kRd, kPh, kSl, kSpace_echo_1, kSpace_echo_2, ..., kSpace_echo_nETL]
+    'images3D': imgMSE,   # images [nSl, nPh, nETL, nRD]
+    'echoSpacing': echoSpacing, # Echo spacing
+    'FOV': fov,           # FOV [sl,ph,rd]
+    't2map': t2Map        # t2Map (ms) [nSl, nPh, nRD]
+}
+
+# T1 Map
+t1map_dict = {
+    'kSpaces3D': sampled, # kSpace [kRd, kPh, kSl, kSpace_echo_1, kSpace_echo_2, ..., kSpace_echo_nETL]
+    'images3D': imgAll,   # images [nSl, nPh, nImg, nRD]
+    'inversionTimes': tI_vector, # Inversion times corresponding each image
+    'FOV': fov,           # FOV [sl,ph,rd]
+    't1map': t1Map        # t2Map (ms) [nSl, nPh, nRD]
+}
+
+# Bo Map
+bomap_dict = {
+    'kSpace3D_0': sampled_0, # kSpace delay 0 [kRd, kPh, kSl, kSpace]
+    'image3D_0': img0,   # images delay 0 [nSl, nPh, nRD]
+    'kSpace3D_delay': sampled_delay, # kSpace delay [kRd, kPh, kSl, kSpace]
+    'image3D_delay': img80,   # images delay [nSl, nPh, nRD]
+    'delay': delay,         # delay (s)
+    'FOV': fov,             # FOV [sl,ph,rd]
+    'boMap': boCeros        # boMap (T) [nSl, nPh, nRD]
+}
+
+# PD 
+pd_dict = {
+    'kSpace3D': sampled_0, # kSpace [kRd, kPh, kSl, kSpace]
+    'image3D': img0,   # images [nSl, nPh, nRD]
+    'FOV': fov,             # FOV [sl,ph,rd]
+}

reco.py

@@ -0,0 +1,78 @@
+import scipy.io as sio
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.optimize import curve_fit
+import sys,os
+from cfl import writecfl
+
+os.environ['TOOLBOX_PATH'] = '/home/jpfitzer/bart-0.9.00'
+os.environ['BART_TOOLBOX_PATH'] = '/home/jpfitzer/bart-0.9.00'
+sys.path.append('/home/jpfitzer/bart-0.9.00/bart/python')
+
+rawName = 'T2_CPMG.mat'
+mat_data_0=sio.loadmat(rawName)
+
+# Format: sampled, # kSpace [kRd, kPh, kSl, kSpace_echo_1, kSpace_echo_2, ..., kSpace_echo_nETL] (102400, 23)
+# So the first three are the coordinates of the kspace, and the rest are the echoes
+print(mat_data_0['kSpaces3D'].shape)
+
+kSpaces3D = mat_data_0['kSpaces3D']
+# self.mapVals['sampled'] = np.concatenate((kRD, kPH, kSL, dataAll_sampled), axis=1)
+
+# nReadout, nPhase, nSlice
+nPoints = (80, 80, 16)
+
+echo_train_length = mat_data_0['kSpaces3D'].shape[1] - 3 # Because the first 3 are kRD, kPH, kSL -> should give 20
+print(f"Echo train length: {echo_train_length}")
+
+echo_spacing = mat_data_0['echoSpacing'][0][0]
+print(f"Echo spacing: {echo_spacing}")
+
+k_readout = kSpaces3D[:, 0]
+k_phase = kSpaces3D[:, 1]
+k_slice = kSpaces3D[:, 2]
+
+# The rest of the data is the echoes
+echos = kSpaces3D[:, 3:]
+
+# Reshape the kspace data for bart
+kSpace = echos.reshape(nPoints[2], nPoints[1], nPoints[0], echo_train_length)
+
+
+
+print(kSpace.shape)
+
+# cfl = writecfl('kSpace', kSpace)
+
+# Create the image with bart fft -i 7 kSpace fft -> three dimensional
+# Put the echos on the fifth dimension:
+# bart transpose 3 5
+# Put the slices on the correct dimension
+# bart transpose 0 2
+
+# traj = writecfl('traj', np.stack((k_readout, k_phase, k_slice), axis=1))
+
+# Echo times with echo spacing
+TE = np.linspace(echo_spacing, echo_spacing * echo_train_length, echo_train_length, endpoint=True)
+print("TE: ", TE)
+
+# Create the echotimes file:
+# bart vec ... echo_times
+# bart scale 0.001 echo_times echo_times_scaled
+# Move the echo_times to the correct dimension:
+# bart transpose 0 5 echo_times_scaled echo_times_final
+
+# Fit the model:
+# bart mobafit -T echo_times_final fft_transposed fit
+
+# Now some values will be very large so we can apply a threshold to obtain a mask
+# bart threshold -M 1000 reco/fit reco/mask
+
+# Multiply the fit with the mask
+# bart fmac reco/fit reco/fit reco/fit_mask
+
+# Select slice
+# bart slice 6 1 reco/fit_mask reco/R2_map
+
+# Invert the data to get T2
+# bart invert reco/R2_map reco/T2_map