13. Calibrationless CS-pMR image reconstruction from undersampled Cartesian data¶
In this tutorial we will reconstruct a 2D MR image from multicoil Cartesian under-sampled kspace measurements.
We use the toy datasets available in pysap, more specifically a 2D brain slice and under-sampled Cartesian acquisition over 32 channels. We compare zero-order image reconstruction with calibrationless multi-coil Compressed sensing reconstructions (analysis vs synthesis formulation) using the FISTA algorithm for the synthesis formulation and the Condat-Vu algorithm for the analysis formulation. Structured sparsity will be promoted in the wavelet domain, using either Symmlet-8 (analysis and synthesis) or undecimated bi-orthogonal wavelets (analysis only) considering group-LASSO or OSCAR-based regularization. The multicoil data $(y_\ell)_\ell$ is collected across multiple, say $L$, channels.
We remind that the synthesis formulation of the Calibrationless CS-PMRI problem reads (minimization in the sparsifying domain):
$$ \widehat{Z} = \text{arg}\,\min_{Z\in C^{n_\Psi\times L}} \frac{1}{2} \sum_{\ell=1}^L\|y_\ell - \Omega F \Psi^*z_\ell \|_2^2 + \lambda {\cal R}(Z) $$
where $Z= [z_1, \ldots, z_L]$ and $X = [x_1,\ldots, x_L]\in C^{n\times L}$ such that $x_l = \Psi^* z_l$. The image solution is given by $\widehat{x} = \Psi^*\widehat{z}$. For an orthonormal wavelet transform, we have $n_\Psi=n$ while for a frame we may have $n_\Psi > n$. The regularization term promotes structured sparsity. For instance when one chooses group-LASSO regularization ${\cal R}(Z) = \sum_{i=1}^{n_\Psi} \|z_i\|_2$, where the L2 norm involves the $L$ channels per wavelet coefficient $z_i$.
The analysis formulation consists in minimizing the following cost function (min. in the image domain):
$$ \widehat{X} = \text{arg}\,\min_{X\in C^{n\times L}} \frac{1}{2} \sum_{\ell=1}^L \|y_\ell - \Omega F x_\ell\|_2^2 + \lambda {\cal R}( \Psi X)\, . $$
- Author: Chaithya G R & Philippe Ciuciu
- Date: 01/07/2021
- Target: ATSI MSc students, Paris-Saclay University
In [37]:
# Package import
from mri.operators import FFT, WaveletN, OWL
'''
FFT ：傅里叶变换算子，用于将 MRI 图像转换到 k-space（频域）。
WaveletN ：正交小波变换，用于稀疏表示 MRI 图像。
OWL：重叠加权 L1 规范（Overlapping Weighted L1），用于稀疏约束。
'''
from mri.reconstructors import CalibrationlessReconstructor
'''
CalibrationlessReconstructor ：用于 MRI 自校准（Calibrationless）重建，不依赖于灵敏度图（SMaps）
'''
from pysap.data import get_sample_data
'''
从 PySAP（Python Sparse Analysis Package）中获取 示例 MRI 数据。
'''
# Third party import
from modopt.opt.proximity import GroupLASSO
'''
Group LASSO 正则化算子，在 稀疏优化（Sparse Optimization）中用于 组稀疏性约束。
'''
from modopt.math.metrics import ssim
'''
计算 结构相似性指数（SSIM），用于衡量重建图像与原始图像之间的相似度。
'''
import numpy as np
import matplotlib.pyplot as plt
In [38]:
# Loading input data
cartesian_ref_image = get_sample_data('2d-pmri').data
image = np.linalg.norm(cartesian_ref_image, axis=0)
# Obtain MRI cartesian mask
mask = get_sample_data("cartesian-mri-mask").data
In [39]:
print(get_sample_data('2d-pmri').data.shape)
print(get_sample_data("cartesian-mri-mask").data.shape)
print(np.where((get_sample_data("cartesian-mri-mask").data)[:, 0] == 1 ))
(32, 512, 512)
(512, 512)
(array([ 8, 21, 35, 38, 66, 74, 95, 100, 130, 134, 142, 146, 166,
173, 177, 180, 182, 183, 185, 195, 197, 198, 203, 206, 207, 210,
215, 219, 220, 221, 222, 223, 225, 226, 232, 235, 237, 238, 239,
241, 244, 245, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256,
257, 258, 259, 260, 261, 262, 264, 265, 266, 267, 270, 271, 272,
277, 278, 283, 285, 290, 292, 293, 294, 295, 297, 301, 302, 304,
306, 320, 326, 339, 341, 345, 358, 387, 388, 393, 396, 400, 411,
436, 447, 462, 465, 488, 489, 506]),)
In [40]:
# View Input
plt.subplot(1, 2, 1)
plt.imshow(np.abs(image), cmap='gray')
plt.title("MRI Data")
plt.subplot(1, 2, 2)
plt.imshow(mask, cmap='gray')
plt.title("K-space Sampling Mask")
plt.show()
$\leadsto$ The code above is to generate the mask, the key part is
mask = get_sample_data("cartesian-mri-mask").data
where we use the package to import the cartesian-mri-mask already existing.
Generate the kspace¶
From the 2D brain slice and the acquisition mask, we retrospectively undersample the k-space using a cartesian acquisition mask We then reconstruct the zero order solution as a baseline
In [41]:
# 2. Define the Fourier operator (FFT)
# Get the locations of the kspace samples and the associated observations
fourier_op = FFT(mask=mask, shape=image.shape,
n_coils=cartesian_ref_image.shape[0])
# 3. Generate the undersampled k-space data
'''
.op(cartesian_ref_image) applies the forward Fourier transform (FFT) only at the sampled locations corresponding to the mask
kspace_obs is a subsampled k-space representation of the image with size (32, 512, 512) .
'''
kspace_obs = fourier_op.op(cartesian_ref_image)
In [42]:
# Zero order solution
# Perform the inverse Fourier transform
zero_soln = np.linalg.norm(fourier_op.adj_op(kspace_obs), axis=0)
base_ssim = ssim(zero_soln, image)
plt.imshow(np.abs(zero_soln), cmap='gray')
plt.title('Zero Order Solution : SSIM = ' + str(np.around(base_ssim, 3)))
plt.show()
$\leadsto$ Here we apply the Fourier operator and its ajoint to perform the zero order solution :
fourier_op = FFT(samples=kspace_loc, shape=image.shape)fourier_op = FFT(mask=mask, shape=image.shape, n_coils=cartesian_ref_image.shape[0])where we define the Fourier transform (FFT) operator only for the sampled locations corresponding to the
mask,then we apply this operator to the
image:kspace_obs = fourier_op.op(cartesian_ref_image)to obtain the undersampled k-space obs data.To compute the zero order solution, we just apply the ajoint of the operator
zero_soln = np.linalg.norm(fourier_op.adj_op(kspace_obs), axis=0), which serves as a basic reconstruction:
Synthesis formulation: FISTA vs POGM optimization¶
We now want to refine the zero order solution using a FISTA optimization. The cost function is set to Proximity Cost + Gradient Cost
In [6]:
# Setup the operators
linear_op = WaveletN(
wavelet_name='sym8',
nb_scale=4,
n_coils=cartesian_ref_image.shape[0],
)
coeffs = linear_op.op(cartesian_ref_image)
regularizer_op = GroupLASSO(weights=6e-8)
Setup reconstructor:¶
In [7]:
# Setup Reconstructor
reconstructor = CalibrationlessReconstructor(
fourier_op=fourier_op,
linear_op=linear_op,
regularizer_op=regularizer_op,
gradient_formulation='synthesis',
verbose=1,
)
WARNING: Making input data immutable.
Lipschitz constant is 1.0999999344348907
In [8]:
# Run the FISTA reconstruction and view results
image_rec, costs, metrics = reconstructor.reconstruct(
kspace_data=kspace_obs,
optimization_alg='fista',
num_iterations=100,
)
image_rec = np.linalg.norm(image_rec, axis=0)
recon_ssim = ssim(image_rec, image)
plt.imshow(np.abs(image_rec), cmap='gray')
plt.title('Iterative FISTA Reconstruction : SSIM = ' + str(np.around(recon_ssim, 3)))
plt.show()
WARNING: Making input data immutable.
- mu: 6e-08 - lipschitz constant: 1.0999999344348907 - data: (512, 512) - wavelet: <mri.operators.linear.wavelet.WaveletN object at 0x7f4bd4ba5790> - 4 - max iterations: 100 - image variable shape: (512, 512) - alpha variable shape: (32, 291721) ---------------------------------------- Starting optimization...
100%|██████████| 100/100 [01:55<00:00, 1.16s/it]
- final iteration number: 100 - final log10 cost value: 6.0 - converged: False Done. Execution time: 115.74814624199644 seconds ----------------------------------------
POGM optimization¶
In [9]:
# Run the POGM reconstruction and view results
image_rec2, costs, metrics = reconstructor.reconstruct(
kspace_data=kspace_obs,
optimization_alg='pogm',
num_iterations=100,
)
image_rec2 = np.linalg.norm(image_rec2, axis=0)
recon2_ssim = ssim(image_rec2, image)
plt.imshow(np.abs(image_rec2), cmap='gray')
plt.title('Iterative POGM Reconstruction : SSIM = ' + str(np.around(recon2_ssim, 3)))
plt.show()
- mu: 6e-08 - lipschitz constant: 1.0999999344348907 - data: (512, 512) - wavelet: <mri.operators.linear.wavelet.WaveletN object at 0x7f4bd4ba5790> - 4 - max iterations: 100 - image variable shape: (32, 512, 512) ---------------------------------------- Starting optimization...
100%|██████████| 100/100 [02:12<00:00, 1.33s/it]
- final iteration number: 100 - final log10 cost value: 6.0 - converged: False Done. Execution time: 132.82787203100452 seconds ----------------------------------------
In [10]:
# Setup the operators
linear_op = WaveletN(
wavelet_name='sym8',
nb_scale=4,
n_coils=cartesian_ref_image.shape[0],
)
coeffs = linear_op.op(cartesian_ref_image)
regularizer_op = OWL(
alpha=1.05e-8,
beta=0,
mode='band_based',
n_coils=cartesian_ref_image.shape[0],
bands_shape=linear_op.coeffs_shape,
)
# Setup Reconstructor
reconstructor = CalibrationlessReconstructor(
fourier_op=fourier_op,
linear_op=linear_op,
regularizer_op=regularizer_op,
gradient_formulation='synthesis',
verbose=1,
)
WARNING: Making input data immutable.
Lipschitz constant is 1.1
In [11]:
# Run the FISTA reconstruction and view results
image_rec, costs, metrics = reconstructor.reconstruct(
kspace_data=kspace_obs,
optimization_alg='fista',
num_iterations=100,
)
image_rec = np.linalg.norm(image_rec, axis=0)
recon_ssim = ssim(image_rec, image)
plt.imshow(np.abs(image_rec), cmap='gray')
plt.title('Iterative Reconstruction : SSIM = ' + str(np.around(recon_ssim, 2)))
plt.show()
- mu: [<modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7f4bd4bb05d0>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7f4bd4bcd450>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7f4bd420b8d0>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7f4bd4209790>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7f4bd4208ad0>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7f4bd4bd5c50>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7f4bd423fcd0>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7f4bd423ee50>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7f4bd423fc50>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7f4bd4bd45d0>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7f4bd423f490>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7f4bd423f8d0>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7f4bd423d610>] - lipschitz constant: 1.1 - data: (512, 512) - wavelet: <mri.operators.linear.wavelet.WaveletN object at 0x7f4bd4b99c10> - 4 - max iterations: 100 - image variable shape: (512, 512) - alpha variable shape: (32, 291721) ---------------------------------------- Starting optimization...
100%|██████████| 100/100 [02:53<00:00, 1.74s/it]
- final iteration number: 100 - final log10 cost value: 6.0 - converged: False Done. Execution time: 173.98121633499977 seconds ----------------------------------------
In [12]:
linear_op = WaveletN(
wavelet_name='sym8',
nb_scale=4,
n_coils=cartesian_ref_image.shape[0],
)
#padding_mode="periodization"
regularizer_op = GroupLASSO(6e-8)
In [13]:
reconstructor = CalibrationlessReconstructor(
fourier_op=fourier_op,
linear_op=linear_op,
regularizer_op=regularizer_op,
gradient_formulation='analysis',
verbose=1,
)
WARNING: Making input data immutable.
Lipschitz constant is 1.1
In [15]:
import pysap
In [16]:
x_final, costs, metrics = reconstructor.reconstruct(
kspace_data=kspace_obs,
optimization_alg='condatvu',
num_iterations=100,
)
image_rec = pysap.Image(data=np.sqrt(np.sum(np.abs(x_final)**2, axis=0)))
recon_ssim = ssim(image_rec, image)
plt.imshow(np.abs(image_rec), cmap='gray')
plt.title('Condat-Vu Reconstruction\nSSIM = ' + str(recon_ssim))
plt.show()
WARNING: <class 'mri.operators.linear.wavelet.WaveletN'> does not inherit an operator parent.
- mu: 6e-08 - lipschitz constant: 1.1 - tau: 0.9523809433107514 - sigma: 0.5 - rho: 1.0 - std: None - 1/tau - sigma||L||^2 >= beta/2: True - data: (512, 512) - wavelet: <mri.operators.linear.wavelet.WaveletN object at 0x7f4bd425cf10> - 4 - max iterations: 100 - number of reweights: 0 - primal variable shape: (32, 512, 512) - dual variable shape: (32, 291721) ---------------------------------------- Starting optimization...
100%|██████████| 100/100 [02:08<00:00, 1.29s/it]
- final iteration number: 100 - final cost value: 1000000.0 - converged: False Done. Execution time: 129.45200876399758 seconds ----------------------------------------
In [22]:
from mri.operators import NonCartesianFFT, WaveletUD2, WaveletN
In [23]:
# Undecimated Wavelets
'''
linear_op = WaveletUD2(
wavelet_id=24,
nb_scale=4,
n_coils=cartesian_ref_image.shape[0],
)
'''
Out[23]:
'\nlinear_op = WaveletUD2(\n wavelet_id=24,\n nb_scale=4,\n n_coils=cartesian_ref_image.shape[0],\n)\n'
In [18]:
#regularizer_op = GroupLASSO(6e-8)
coeffs = linear_op.op(cartesian_ref_image)
regularizer_op = OWL(
alpha=1.05e-8,
beta=0,
mode='band_based',
n_coils=cartesian_ref_image.shape[0],
bands_shape=linear_op.coeffs_shape,
)
In [19]:
reconstructor = CalibrationlessReconstructor(
fourier_op=fourier_op,
linear_op=linear_op,
regularizer_op=regularizer_op,
gradient_formulation='analysis',
verbose=1,
)
WARNING: Making input data immutable.
Lipschitz constant is 1.1000001311302186
In [20]:
x_final, costs, metrics = reconstructor.reconstruct(
kspace_data=kspace_obs,
optimization_alg='condatvu',
num_iterations=100,
)
image_rec = pysap.Image(data=np.sqrt(np.sum(np.abs(x_final)**2, axis=0)))
recon_ssim = ssim(image_rec, image)
plt.imshow(np.abs(image_rec), cmap='gray')
plt.title('Condat-Vu Reconstruction\nSSIM = ' + str(recon_ssim))
plt.show()
WARNING: <class 'mri.operators.linear.wavelet.WaveletN'> does not inherit an operator parent.
- mu: [<modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7f4bd4bbf210>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7f4bd4bd5390>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7f4bcb43e2d0>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7f4bcb459390>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7f4bf20626d0>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7f4bd411c090>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7f4dac12e110>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7f4dac12f2d0>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7f4bd42e8790>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7f4bd411ea10>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7f4bd42ebdd0>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7f4bd42e9910>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7f4bd411b690>] - lipschitz constant: 1.1000001311302186 - tau: 0.9523808838412695 - sigma: 0.5 - rho: 1.0 - std: None - 1/tau - sigma||L||^2 >= beta/2: True - data: (512, 512) - wavelet: <mri.operators.linear.wavelet.WaveletN object at 0x7f4bd425cf10> - 4 - max iterations: 100 - number of reweights: 0 - primal variable shape: (32, 512, 512) - dual variable shape: (32, 291721) ---------------------------------------- Starting optimization...
100%|██████████| 100/100 [03:26<00:00, 2.06s/it]
- final iteration number: 100 - final cost value: 1000000.0 - converged: False Done. Execution time: 207.19260900700465 seconds ----------------------------------------
Q1
Observe and comment the image quality using calibrationless reconstruction with either the Group-LASSO or the Ordered Weighted L1-norm (OWL) regularization.
Which one provides the best image quality? Why?
$\leadsto$ First, we try using linear_op = WaveletUD2, gradient_formulation='analysis' and optimization_alg='condatvu', in num_iterations=100, the reconstruction takes approximately 1h30 , having the lowest image quality and the computation time is the longest.
In [28]:
import pysap
In [ ]:
linear_op = WaveletUD2(
wavelet_id=24,
nb_scale=4,
n_coils=cartesian_ref_image.shape[0],
)
coeffs = linear_op.op(cartesian_ref_image)
regularizer_op = OWL(
alpha=1.05e-8,
beta=0,
mode='band_based',
n_coils=cartesian_ref_image.shape[0],
bands_shape=linear_op.coeffs_shape,
)
reconstructor = CalibrationlessReconstructor(
fourier_op=fourier_op,
linear_op=linear_op,
regularizer_op=regularizer_op,
gradient_formulation='analysis',
verbose=1,
)
x_final, costs, metrics = reconstructor.reconstruct(
kspace_data=kspace_obs,
optimization_alg='condatvu',
num_iterations=100,
)
image_rec = pysap.Image(data=np.sqrt(np.sum(np.abs(x_final)**2, axis=0)))
recon_ssim = ssim(image_rec, image)
plt.imshow(np.abs(image_rec), cmap='gray')
plt.title('Condat-Vu Reconstruction\nSSIM = ' + str(recon_ssim))
plt.show()
$\leadsto$ We decide to use linear_op = WaveletN, gradient_formulation='synthesis' and optimization_alg='fista'
In [30]:
# Setup the operators
linear_op = WaveletN(
wavelet_name='sym8',
nb_scale=4,
n_coils=cartesian_ref_image.shape[0],
)
coeffs = linear_op.op(cartesian_ref_image)
regularizer_op = GroupLASSO(weights=6e-8)
# Setup Reconstructor
reconstructor = CalibrationlessReconstructor(
fourier_op=fourier_op,
linear_op=linear_op,
regularizer_op=regularizer_op,
gradient_formulation='synthesis',
verbose=1,
)
# Run the FISTA reconstruction and view results
image_rec, costs, metrics = reconstructor.reconstruct(
kspace_data=kspace_obs,
optimization_alg='fista',
num_iterations=100,
)
image_rec = np.linalg.norm(image_rec, axis=0)
recon_ssim = ssim(image_rec, image)
plt.imshow(np.abs(image_rec), cmap='gray')
plt.title('Iterative FISTA Reconstruction : SSIM = ' + str(np.around(recon_ssim, 3)))
plt.show()
WARNING: Making input data immutable.
Lipschitz constant is 1.0999999344348907 - mu: 6e-08 - lipschitz constant: 1.0999999344348907 - data: (512, 512) - wavelet: <mri.operators.linear.wavelet.WaveletN object at 0x7ff5c1dfea10> - 4 - max iterations: 100 - image variable shape: (512, 512) - alpha variable shape: (32, 291721) ---------------------------------------- Starting optimization...
100%|██████████| 100/100 [02:11<00:00, 1.32s/it]
- final iteration number: 100 - final log10 cost value: 6.0 - converged: False Done. Execution time: 131.70454503200017 seconds ----------------------------------------
In [29]:
# Setup the operators
linear_op = WaveletN(
wavelet_name='sym8',
nb_scale=4,
n_coils=cartesian_ref_image.shape[0],
)
coeffs = linear_op.op(cartesian_ref_image)
regularizer_op = OWL(
alpha=1.05e-8,
beta=0,
mode='band_based',
n_coils=cartesian_ref_image.shape[0],
bands_shape=linear_op.coeffs_shape,
)
# Setup Reconstructor
reconstructor = CalibrationlessReconstructor(
fourier_op=fourier_op,
linear_op=linear_op,
regularizer_op=regularizer_op,
gradient_formulation='synthesis',
verbose=1,
)
# Run the FISTA reconstruction and view results
image_rec, costs, metrics = reconstructor.reconstruct(
kspace_data=kspace_obs,
optimization_alg='fista',
num_iterations=100,
)
image_rec = np.linalg.norm(image_rec, axis=0)
recon_ssim = ssim(image_rec, image)
plt.imshow(np.abs(image_rec), cmap='gray')
plt.title('Iterative Reconstruction : SSIM = ' + str(np.around(recon_ssim, 2)))
plt.show()
WARNING: Making input data immutable.
Lipschitz constant is 1.1 - mu: [<modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7ff5c3b871d0>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7ff5bfc7e7d0>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7ff5bfcf2750>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7ff5bfcf2a50>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7ff5bfcf2c10>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7ff5c1d98d10>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7ff5c1d99250>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7ff5c1e1aa90>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7ff5c1e18cd0>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7ff5c1d9bdd0>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7ff5bfcfb790>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7ff5bfcf9650>, <modopt.opt.proximity.OrderedWeightedL1Norm object at 0x7ff5bfcfbfd0>] - lipschitz constant: 1.1 - data: (512, 512) - wavelet: <mri.operators.linear.wavelet.WaveletN object at 0x7ff5bfcb0610> - 4 - max iterations: 100 - image variable shape: (512, 512) - alpha variable shape: (32, 291721) ---------------------------------------- Starting optimization...
100%|██████████| 100/100 [03:03<00:00, 1.84s/it]
- final iteration number: 100 - final log10 cost value: 6.0 - converged: False Done. Execution time: 183.59860118600773 seconds ----------------------------------------
In [53]:
# 导入 MRI 重建运算符
# Package import
from mri.operators import FFT, WaveletN
'''
WaveletN : 正交小波变换（如 Symmlet-8 或 Daubechies），用于 合成公式（Synthesis Formulation） 压缩感知重建。
WaveletUD2 : 未降采样双正交小波变换（Undecimated Bi-Orthogonal Wavelet Transform），用于 分析公式（Analysis Formulation） 压缩感知重建。
'''
# 导入 MRI 处理的实用工具
from mri.operators.utils import convert_mask_to_locations
'''
convert_mask_to_locations ：2D binary mask -> a list of sampled k-space coordinates (x, y)
'''
# 导入 MRI 图像重建器
from mri.reconstructors.utils.extract_sensitivity_maps \
import get_Smaps, extract_k_space_center_and_locations
from mri.reconstructors import SelfCalibrationReconstructor
'''
get_Smaps ：从 k-space 提取 灵敏度图（Sensitivity Maps），用于 多通道 MRI 重建（pMRI）。
extract_k_space_center_and_locations ：获取 K-space 采样中心，用于自校准（Self-Calibration）。
SelfCalibrationReconstructor ：用于 基于自校准（Self-Calibrated）策略 进行 MRI 图像重建。
'''
# 导入 MRI 样本数据
from pysap.data import get_sample_data
# Third party import
from modopt.math.metrics import ssim
from modopt.opt.linear import Identity
from modopt.opt.proximity import SparseThreshold
'''
ssim : 结构相似性指数（SSIM），用于 评估重建图像的质量，与原始 MRI 进行比较。
Identity : 单位变换（Identity Transformation），即 不对数据做任何变换，在优化过程中常用于保持数据结构。
SparseThreshold : 稀疏阈值运算（L1 正则化 / 软阈值收缩），用于 压缩感知重建中的稀疏性约束。
'''
import numpy as np
import matplotlib.pyplot as plt
In [54]:
# 2. Define the Fourier operator (FFT)
# Get the locations of the kspace samples and the associated observations
fourier_op = FFT(mask=mask, shape=image.shape,
n_coils=cartesian_ref_image.shape[0])
# 3. Generate the undersampled k-space data
'''
.op(cartesian_ref_image) applies the forward Fourier transform (FFT) only at the sampled locations corresponding to the mask
kspace_obs is a subsampled k-space representation of the image with size (32, 512, 512) .
'''
kspace_obs = fourier_op.op(cartesian_ref_image)
In [62]:
import time
In [72]:
def reconstruction (regularizer_op_, num_iterations_) :
if regularizer_op_ in ["GroupLASSO", "OWL"] :
time_0 = time.time()
# Setup the operators
linear_op = WaveletN(
wavelet_name='sym8',
nb_scale=4,
n_coils=cartesian_ref_image.shape[0],
)
coeffs = linear_op.op(cartesian_ref_image)
if regularizer_op_ == "GroupLASSO" :
regularizer_op = GroupLASSO(weights=6e-8)
elif regularizer_op_ == "OWL" :
regularizer_op = OWL(
alpha=1.05e-8,
beta=0,
mode='band_based',
n_coils=cartesian_ref_image.shape[0],
bands_shape=linear_op.coeffs_shape,
)
# Setup Reconstructor
reconstructor = CalibrationlessReconstructor(
fourier_op=fourier_op,
linear_op=linear_op,
regularizer_op=regularizer_op,
gradient_formulation='synthesis',
verbose=1,
)
# Run the FISTA reconstruction and view results
image_rec, costs, metrics = reconstructor.reconstruct(
kspace_data=kspace_obs,
optimization_alg='fista',
num_iterations=num_iterations_,
)
image_rec = np.linalg.norm(image_rec, axis=0)
recon_ssim = ssim(image_rec, image)
time_1 = time.time()
elif regularizer_op_ == "SparseThreshold" :
time_0 = time.time()
linear_op = WaveletN(
wavelet_name='sym8',
nb_scale=4,
)
regularizer_op = SparseThreshold(Identity(), 1.5e-8, thresh_type="soft")
# Setup Reconstructor
reconstructor = SelfCalibrationReconstructor(
fourier_op=fourier_op,
linear_op=linear_op,
regularizer_op=regularizer_op,
gradient_formulation='synthesis',
kspace_portion=0.01,
verbose=1,
)
image_rec, costs, metrics = reconstructor.reconstruct(
kspace_data=kspace_obs,
optimization_alg='fista',
num_iterations=num_iterations_,
)
recon_ssim = ssim(image_rec, image)
time_1 = time.time()
time_computation = time_1 - time_0
return image_rec, recon_ssim, time_computation
In [ ]:
vector_regularizer_op=["GroupLASSO", "OWL", "SparseThreshold"]
vector_num_iterations = [20, 100, 400, 500]
fig, axs = plt.subplots(len(vector_num_iterations), 3, figsize=(12, 12))
matrix_time_computation = np.zeros((len(vector_num_iterations), 3))
for i in range(len(vector_num_iterations)) :
for j in range(3):
num_iterations_ = vector_num_iterations[i]
regularizer_op_ = vector_regularizer_op[j]
image_rec, recon_ssim, time_computation = reconstruction(regularizer_op_, num_iterations_)
# 显示重建的图像
axs[i, j].imshow(np.abs(image_rec), cmap='gray')
axs[i, j].set_xlabel(f"{regularizer_op_} : SSIM = {recon_ssim:.3f}")
axs[i, j].xaxis.set_label_position('top') # 移动 xlabel 到顶部
axs[i, j].set_ylabel(f"num_iter = {num_iterations_}")
matrix_time_computation[i][j] = time_computation
axs[0, 0].set_title("'WaveletN' + 'fista' : calibrationless")
axs[0, 1].set_title("'WaveletN' + 'fista' : calibrationless")
axs[0, 2].set_title("'WaveletN' + 'fista' : self-calibrated ")
# Add a main title (suptitle) to the figure
fig.suptitle("Iterative Reconstruction Results")
# Display the figure with adjusted layout
plt.tight_layout(rect=[0, 0, 1, 0.96]) # Adjust layout to fit suptitle
plt.show()
In [ ]:
vector_regularizer_op=["GroupLASSO", "OWL", "SparseThreshold"]
vector_num_iterations = np.concatenate([
np.linspace(1, 80, 20),
np.linspace(81, 100, 2)
])
matrix_recon_ssim = np.zeros((len(vector_num_iterations), 3))
matrix_time_computation = np.zeros((len(vector_num_iterations), 3))
for i in range(len(vector_num_iterations)) :
for j in range(3):
num_iterations_ = vector_num_iterations[i]
regularizer_op_ = vector_regularizer_op[j]
image_rec, recon_ssim, time_computation = reconstruction(regularizer_op_, int(num_iterations_))
matrix_recon_ssim[i][j] = recon_ssim
matrix_time_computation[i][j] = time_computation
# Add a main title (suptitle) to the figure
fig.suptitle("Iterative Reconstruction Results")
# Display the figure with adjusted layout
plt.tight_layout(rect=[0, 0, 1, 0.96]) # Adjust layout to fit suptitle
plt.show()
plt.figure(figsize=(10, 6))
# WaveletN Synthesis FISTA plot "GroupLASSO"
plt.plot(vector_num_iterations, matrix_recon_ssim[:, 0], linestyle='-', color='red',
markersize=8, markerfacecolor='red', label='WaveletN Synthesis FISTA : regularizer_op = GroupLASSO | calibrationless')
# WaveletN Synthesis FISTA plot "OWL"
plt.plot(vector_num_iterations, matrix_recon_ssim[:, 1], linestyle='-', color='blue',
markersize=8, markerfacecolor='blue', label='WaveletN Synthesis FISTA : regularizer_op = OWL | calibrationless')
# WaveletN Synthesis FISTA plot "OWL"
plt.plot(vector_num_iterations, matrix_recon_ssim[:, 2], linestyle='-', color='green',
markersize=8, markerfacecolor='green', label='WaveletN Synthesis FISTA : regularizer_op = SparseThreshold | self-calibrated')
# Labels and title
plt.xlabel('num_iterations')
plt.ylabel('recon_ssim')
plt.title('SSIM Score vs Number of Iterations')
# Grid and legend
# plt.grid(True)
plt.legend()
# Show the plot
plt.show()
$\leadsto$ The calibrationless reconstruction using the Ordered Weighted L1-norm (OWL) regularization leads to a better reconstructed image quality compared to using the Group-LASSO regularization :
OWL applies the non-uniform shrinkage to the wavelet coefficients however Group-LASSO regularization enforces block sparsity. Thus the OWL is more accurate than the Group-LASSO regularization, it effectively handles the high frenquency components and suppressing of noise.
Q2
Compare calibrationless reconstruction to the self-calibrated version you run previously in the Cartesian setting (in 3.), both in terms of image quality and computing time.
$\leadsto$ Self-calibrated reconstruction results in both improved reconstructed image quality and reduced computing time compared to the calibrationless reconstruction.
In [78]:
# Add a main title (suptitle) to the figure
fig.suptitle("Computing Time")
# Display the figure with adjusted layout
plt.tight_layout(rect=[0, 0, 1, 0.96]) # Adjust layout to fit suptitle
plt.show()
plt.figure(figsize=(10, 6))
# WaveletN Synthesis FISTA plot "GroupLASSO"
plt.plot(vector_num_iterations, matrix_time_computation[:, 0], linestyle='-', color='red',
markersize=8, markerfacecolor='red', label='WaveletN Synthesis FISTA : regularizer_op = GroupLASSO | calibrationless')
# WaveletN Synthesis FISTA plot "OWL"
plt.plot(vector_num_iterations, matrix_time_computation[:, 1], linestyle='-', color='blue',
markersize=8, markerfacecolor='blue', label='WaveletN Synthesis FISTA : regularizer_op = OWL | calibrationless')
# WaveletN Synthesis FISTA plot "OWL"
plt.plot(vector_num_iterations, matrix_time_computation[:, 2], linestyle='-', color='green',
markersize=8, markerfacecolor='green', label='WaveletN Synthesis FISTA : regularizer_op = SparseThreshold | self-calibrated')
# Labels and title
plt.xlabel('num_iterations')
plt.ylabel('matrix_time_computation')
plt.title('Computing Time vs Number of Iterations')
# Grid and legend
# plt.grid(True)
plt.legend()
# Show the plot
plt.show()
<Figure size 640x480 with 0 Axes>
$\leadsto$