Parallel Magnetic Resonance imaging (MRI)¶
Standard MRI is described by the Fourier transform \(\mathcal{F}\) as forward operator (here in two dimensions). To accelerate data acquisition, parallel MRI uses simultaneous measurements by \(N\) receiver coils. This allows undersampling of the Fourier domain leading to speed-ups. Parallel MRI is described by the forward operator
\[\begin{split}F\left(\begin{array}{c}\rho \\ c_1\\ \vdots \\ c_N\end{array}\right)
= \left(\begin{array}{c}M\cdot\mathcal{F}(c_1 \cdot \rho)\\ \vdots \\
M\cdot\mathcal{F}(c_N \cdot \rho)\end{array}\right).\end{split}\]
Here \(\rho\) describes the hydrogen density and is the main quantity of interest. To take into account effects such as motion artifacts, \(\rho\) has to be modeled as a complex-valued function. \(c_1,\dots, c_N\) describe complex-valued coil profile, which may be assumed to be smooth. As they depend on the sample \(\rho\), they must be reconstructed together with \(\rho\). \(M\) is a 0-1-mask describing the undersampling pattern.
[1]:
import logging
import matplotlib.pyplot as plt
import matplotlib as mplib
from matplotlib.colors import hsv_to_rgb
import numpy as np
from scipy.io import loadmat
import regpy.stoprules as rules
import regpy.util as uti
from examples.mri.mri import cartesian_sampling, normalize, parallel_mri, sobolev_smoother, estimate_sampling_pattern
from regpy.operators import PtwMultiplication
from regpy.solvers import RegularizationSetting
from regpy.solvers.nonlinear.irgnm import IrgnmCG
from regpy.vecsps import UniformGridFcts
from regpy.hilbert import L2
import os
logging.basicConfig(
level=logging.INFO,
format='%(asctime)s %(levelname)s %(name)-40s :: %(message)s'
)
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[1], line 11
9 import regpy.stoprules as rules
10 import regpy.util as uti
---> 11 from examples.mri.mri import cartesian_sampling, normalize, parallel_mri, sobolev_smoother, estimate_sampling_pattern
12 from regpy.operators import PtwMultiplication
13 from regpy.solvers import RegularizationSetting
ModuleNotFoundError: No module named 'examples'
Complex to rgb conversion¶
Converts array of complex numbers into array of RGB color values for plotting. The hue corresponds to the argument. The brighntess corresponds to the absolute value.
Parameters >z : numpy.ndarray >array of complex numbers
Returns > numpy.ndarray > Array that contains three values for each value in z containing the RGB representation of this value.
[2]:
def complex_to_rgb(z):
HSV = np.dstack( (np.mod(np.angle(z)/(2.*np.pi),1), 1.0*np.ones(z.shape), np.abs(z)/np.max((np.abs(z[:]))), ))
return hsv_to_rgb(HSV)
Load data from file and estimate sampling pattern¶
[3]:
data = loadmat('data/ksp3x2.mat')['Y']
data = np.transpose(data,(2,0,1))*(100/np.linalg.norm(data))
# normalize and transpose data
nrcoils,n1,n2 = data.shape
grid = UniformGridFcts((-1, 1, n1), (-1, 1, n2), dtype=complex)
mask = estimate_sampling_pattern(data)
plt.imshow(mask.T); plt.title('Undersampling pattern of data')
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
File /usr/local/lib/python3.12/site-packages/scipy/io/matlab/_mio.py:39, in _open_file(file_like, appendmat, mode)
38 try:
---> 39 return open(file_like, mode), True
40 except OSError as e:
41 # Probably "not found"
FileNotFoundError: [Errno 2] No such file or directory: 'data/ksp3x2.mat'
During handling of the above exception, another exception occurred:
FileNotFoundError Traceback (most recent call last)
Cell In[3], line 1
----> 1 data = loadmat('data/ksp3x2.mat')['Y']
2 data = np.transpose(data,(2,0,1))*(100/np.linalg.norm(data))
3 # normalize and transpose data
File /usr/local/lib/python3.12/site-packages/scipy/io/matlab/_mio.py:233, in loadmat(file_name, mdict, appendmat, spmatrix, **kwargs)
88 """
89 Load MATLAB file.
90
(...) 230 3.14159265+3.14159265j])
231 """
232 variable_names = kwargs.pop('variable_names', None)
--> 233 with _open_file_context(file_name, appendmat) as f:
234 MR, _ = mat_reader_factory(f, **kwargs)
235 matfile_dict = MR.get_variables(variable_names)
File /usr/local/lib/python3.12/contextlib.py:137, in _GeneratorContextManager.__enter__(self)
135 del self.args, self.kwds, self.func
136 try:
--> 137 return next(self.gen)
138 except StopIteration:
139 raise RuntimeError("generator didn't yield") from None
File /usr/local/lib/python3.12/site-packages/scipy/io/matlab/_mio.py:17, in _open_file_context(file_like, appendmat, mode)
15 @contextmanager
16 def _open_file_context(file_like, appendmat, mode='rb'):
---> 17 f, opened = _open_file(file_like, appendmat, mode)
18 try:
19 yield f
File /usr/local/lib/python3.12/site-packages/scipy/io/matlab/_mio.py:45, in _open_file(file_like, appendmat, mode)
43 if appendmat and not file_like.endswith('.mat'):
44 file_like += '.mat'
---> 45 return open(file_like, mode), True
46 else:
47 raise OSError(
48 'Reader needs file name or open file-like object'
49 ) from e
FileNotFoundError: [Errno 2] No such file or directory: 'data/ksp3x2.mat'
Set up forward operator¶
[4]:
sobolev_index = 32
full_mri_op = parallel_mri(grid=grid, ncoils=nrcoils,centered=True)
sampling = PtwMultiplication(full_mri_op.codomain,(1.+0j)* mask)
#cartesian_sampling(full_mri_op.codomain, mask=mask)
smoother = sobolev_smoother(full_mri_op.domain, sobolev_index, factor=220.)
parallel_mri_op = sampling * full_mri_op * smoother
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[4], line 3
1 sobolev_index = 32
----> 3 full_mri_op = parallel_mri(grid=grid, ncoils=nrcoils,centered=True)
4 sampling = PtwMultiplication(full_mri_op.codomain,(1.+0j)* mask)
5 #cartesian_sampling(full_mri_op.codomain, mask=mask)
NameError: name 'parallel_mri' is not defined
Set up initial guess¶
We use constant density and zero coil profiles as initial guess.
[5]:
init = parallel_mri_op.domain.zeros()
init_density, _ = parallel_mri_op.domain.split(init)
init_density[...] = 1
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[5], line 1
----> 1 init = parallel_mri_op.domain.zeros()
2 init_density, _ = parallel_mri_op.domain.split(init)
3 init_density[...] = 1
NameError: name 'parallel_mri_op' is not defined
Set up regularization method¶
[6]:
setting = RegularizationSetting(op=parallel_mri_op, penalty=L2, data_fid=L2)
solver = IrgnmCG(
setting=setting,
data=data,
regpar=1,
regpar_step=1/3.,
init=init
)
stoprule = rules.CountIterations(max_iterations=5)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[6], line 1
----> 1 setting = RegularizationSetting(op=parallel_mri_op, penalty=L2, data_fid=L2)
3 solver = IrgnmCG(
4 setting=setting,
5 data=data,
(...) 8 init=init
9 )
11 stoprule = rules.CountIterations(max_iterations=5)
NameError: name 'RegularizationSetting' is not defined
Run solver by hand and plot iterates¶
Get an iterator from the solver
[7]:
it = iter(solver)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[7], line 1
----> 1 it = iter(solver)
NameError: name 'solver' is not defined
Iterate this cell by hand about 5 times
[8]:
"""
reco, reco_data = next(it)
rho, coils = smoother.codomain.split(smoother(reco))
#rho, coils = normalize(rho,coils)
fig = plt.figure(figsize = (15,9))
gs = fig.add_gridspec(3,7)
axs = [fig.add_subplot(gs[0:3, 0:3])]
im0 = axs[0].imshow(np.abs(rho),cmap=mplib.colormaps['Greys_r'],origin='lower')
axs[0].xaxis.set_ticklabels([])
axs[0].yaxis.set_ticklabels([])
for j in range(3):
for k in range(3,7):
axs.append(fig.add_subplot(gs[j,k]))
axs[-1].xaxis.set_ticklabels([])
axs[-1].yaxis.set_ticklabels([])
for j in range(nrcoils):
axs[1+j].imshow(complex_to_rgb(coils[j,:,:]),origin='lower')
plt.show()
"""
[8]:
"\nreco, reco_data = next(it)\n\nrho, coils = smoother.codomain.split(smoother(reco))\n#rho, coils = normalize(rho,coils)\n\nfig = plt.figure(figsize = (15,9))\n\ngs = fig.add_gridspec(3,7)\naxs = [fig.add_subplot(gs[0:3, 0:3])]\nim0 = axs[0].imshow(np.abs(rho),cmap=mplib.colormaps['Greys_r'],origin='lower')\naxs[0].xaxis.set_ticklabels([])\naxs[0].yaxis.set_ticklabels([])\nfor j in range(3):\n for k in range(3,7):\n axs.append(fig.add_subplot(gs[j,k]))\n axs[-1].xaxis.set_ticklabels([])\n axs[-1].yaxis.set_ticklabels([])\nfor j in range(nrcoils):\n axs[1+j].imshow(complex_to_rgb(coils[j,:,:]),origin='lower')\nplt.show()\n"
… or run solver by stopping rule¶
[9]:
for reco, reco_data in solver.while_(stoprule):
rho, coils = smoother.codomain.split(smoother(reco))
#rho, coils = normalize(rho,coils)
fig = plt.figure(figsize = (15,9))
gs = fig.add_gridspec(3,7)
axs = [fig.add_subplot(gs[0:3, 0:3])]
axs[0].imshow(np.abs(rho),cmap=mplib.colormaps['Greys_r'],origin='lower')
axs[0].xaxis.set_ticklabels([])
axs[0].yaxis.set_ticklabels([])
for j in range(3):
for k in range(3,7):
axs.append(fig.add_subplot(gs[j,k]))
axs[-1].xaxis.set_ticklabels([])
axs[-1].yaxis.set_ticklabels([])
for j in range(nrcoils):
axs[1+j].imshow(complex_to_rgb(coils[j,:,:]),origin='lower')
plt.show()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[9], line 1
----> 1 for reco, reco_data in solver.while_(stoprule):
2 rho, coils = smoother.codomain.split(smoother(reco))
3 #rho, coils = normalize(rho,coils)
NameError: name 'solver' is not defined