.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "generated/autoexamples/example04_calibration.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_generated_autoexamples_example04_calibration.py: =========================================== Utils Tour: Coil Compression and NLINV =========================================== This example introduces two utility routines from :mod:`pygrog.utils`: 1. **Coil compression** (:func:`~pygrog.utils.coil_compress`) — reduces the number of receiver coils via PCA to speed up downstream processing without significant SNR loss. 2. **NLINV coil calibration** (:func:`~pygrog.utils.nlinv_calib`) — estimates coil sensitivity maps from undersampled k-space data using the nonlinear-inverse (NLINV) algorithm. All data use a **BrainWeb T1-weighted** phantom: a multi-coil dataset is constructed by multiplying it with synthetic coil sensitivity maps. .. GENERATED FROM PYTHON SOURCE LINES 18-33 .. code-block:: Python import matplotlib.colors as mcolors import matplotlib.pyplot as plt import numpy as np import torch from brainweb_dl import get_mri image = get_mri(0, "T1") .. rst-class:: sphx-glr-script-out .. code-block:: none Coil images shape: (16, 217, 181) K-space shape : (16, 217, 181) .. GENERATED FROM PYTHON SOURCE LINES 107-110 .. code-block:: Python fig, axes = plt.subplots(2, n_coils // 2, figsize=(11, 4)) .. image-sg:: /generated/autoexamples/images/sphx_glr_example04_calibration_001.png :alt: Synthetic multi-coil images (16 coils), Coil 1, Coil 2, Coil 3, Coil 4, Coil 5, Coil 6, Coil 7, Coil 8, Coil 9, Coil 10, Coil 11, Coil 12, Coil 13, Coil 14, Coil 15, Coil 16 :srcset: /generated/autoexamples/images/sphx_glr_example04_calibration_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 121-128 Coil Compression via PCA ========================= Use :func:`~pygrog.utils.coil_compress` to reduce the number of receiver coils via PCA on k-space data, reducing computational cost with minimal SNR loss. The function returns the compressed k-space and the compression matrix ``W`` of shape ``(n_virtual, n_coils)``. .. GENERATED FROM PYTHON SOURCE LINES 128-132 .. code-block:: Python # Prepare k-space data for compression kspace_flat = kspace_full.reshape(n_coils, -1) # (n_coils, n_samples) .. GENERATED FROM PYTHON SOURCE LINES 133-155 .. code-block:: Python n_virtual = 8 from pygrog.utils import coil_compress kspace_cc, W = coil_compress(kspace_flat, n_virtual) print(f"\nOriginal coils : {kspace_flat.shape[0]}") print(f"Virtual coils : {kspace_cc.shape[0]}") print(f"Compression matrix: {W.shape}") # Reconstruct coil images from compressed data kspace_cc_full = kspace_cc.reshape(n_virtual, *image_shape) coil_images_cc = np.fft.fftshift( np.fft.ifft2(np.fft.ifftshift(kspace_cc_full, axes=(-2, -1))), axes=(-2, -1) ) # RSS of compressed coils vs original rss_orig = np.sqrt((np.abs(coil_images) ** 2).sum(0)) rss_cc = np.sqrt((np.abs(coil_images_cc) ** 2).sum(0)) .. rst-class:: sphx-glr-script-out .. code-block:: none Original coils : 16 Virtual coils : 8 Compression matrix: (8, 16) .. GENERATED FROM PYTHON SOURCE LINES 156-160 .. code-block:: Python fig, axes = plt.subplots(1, 3, figsize=(11, 3.5)) .. image-sg:: /generated/autoexamples/images/sphx_glr_example04_calibration_002.png :alt: Coil compression, Reference, RSS (16 coils), RSS (8 virtual coils, PCA) :srcset: /generated/autoexamples/images/sphx_glr_example04_calibration_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 175-180 .. note:: Coil compression is lossless when ``n_coils >= original_n_coils`` and near-lossless when the retained energy fraction is high (e.g. > 0.99). Use a float threshold ``coil_compress(data, 0.99)`` for automatic rank selection. .. GENERATED FROM PYTHON SOURCE LINES 182-193 NLINV Sensitivity Estimation ============================= :func:`~pygrog.utils.nlinv_calib` estimates coil sensitivity maps jointly with the image using the nonlinear-inverse (NLINV) algorithm (Uecker et al.\ 2008). It requires a small fully-sampled ACS centre to bootstrap, but no explicit prior knowledge of the sensitivity maps. The function supports two modes: * **``cal_width=None``** - full k-space calibrationless reconstruction * **``cal_width=24``** - fast calibration mode with central k-space cropping .. GENERATED FROM PYTHON SOURCE LINES 195-223 .. code-block:: Python from pygrog.utils import nlinv_calib acs = 8 # fully-sampled ACS centre width (columns), as in MATLAB reference cal_width = 24 # NLINV internal calibration resolution for Step 2 # Cartesian undersampling mask: R=2 readout subsampling + 8-col ACS centre mask = np.zeros(image_shape, dtype=bool) mask[:, ::2] = True # R=2 subsampling on readout (column) direction cx = image_shape[1] // 2 cy = image_shape[0] // 2 mask[cy - acs // 2 : cy + acs // 2, cx - acs // 2 : cx + acs // 2] = True # ACS kspace_us = kspace_full * mask[np.newaxis] print(f"\nUndersampling factor : {mask.size / mask.sum():.2f}x") print(f"ACS region : all rows x {acs} cols") print(f"Undersampled k-space : {kspace_us.shape}") ncols_show = min(n_coils // 2, 4) # Zero-filled RSS image from undersampled data (baseline) coil_images_us = np.fft.fftshift( np.fft.ifft2(np.fft.ifftshift(kspace_us, axes=(-2, -1))), axes=(-2, -1) ) rss_us = np.sqrt((np.abs(coil_images_us) ** 2).sum(0)) .. rst-class:: sphx-glr-script-out .. code-block:: none Undersampling factor : 1.99x ACS region : all rows x 8 cols Undersampled k-space : (16, 217, 181) .. GENERATED FROM PYTHON SOURCE LINES 224-230 NLINV Step 1: Full k-space calibrationless reconstruction ========================================================== Call :func:`~pygrog.utils.nlinv_calib` with ``cal_width=None`` to solve jointly for image and sensitivity maps using the entire undersampled k-space. cal_width=None -> full k-space, no cropping, matching MATLAB reference .. GENERATED FROM PYTHON SOURCE LINES 230-249 .. code-block:: Python smaps_full, _, image_full = nlinv_calib( kspace_us, cal_width=None, ndim=2, mask=mask, ret_cal=True, ret_image=True, ) smaps_full_np = ( smaps_full.numpy() if isinstance(smaps_full, torch.Tensor) else smaps_full ) image_full_np = ( image_full.numpy() if isinstance(image_full, torch.Tensor) else image_full ) print(f"\n[Step 1] Smaps shape : {smaps_full_np.shape}") print(f"[Step 1] NLINV image shape: {image_full_np.shape}") .. rst-class:: sphx-glr-script-out .. code-block:: none [Step 1] Smaps shape : (16, 217, 181) [Step 1] NLINV image shape: (217, 181) .. GENERATED FROM PYTHON SOURCE LINES 250-253 .. code-block:: Python fig, axes = plt.subplots(3, ncols_show + 1, figsize=(3 * (ncols_show + 1), 8)) .. image-sg:: /generated/autoexamples/images/sphx_glr_example04_calibration_003.png :alt: NLINV - calibrationless image reconstruction (217x181), Reference, GT smap 1, GT smap 2, GT smap 3, GT smap 4, RSS (zero-filled), NLINV smap 1, NLINV smap 2, NLINV smap 3, NLINV smap 4, NLINV image, NLINV smap 1, NLINV smap 2, NLINV smap 3, NLINV smap 4, RGB key :srcset: /generated/autoexamples/images/sphx_glr_example04_calibration_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 287-293 NLINV Step 2: Calibration mode with central k-space cropping ============================================================ Call :func:`~pygrog.utils.nlinv_calib` with ``cal_width=24`` to crop to the central region and solve for low-resolution sensitivity maps and a synthesized calibration k-space patch (useful for GRAPPA/GROG training). .. GENERATED FROM PYTHON SOURCE LINES 293-313 .. code-block:: Python smaps_cal, grappa_train, image_cal = nlinv_calib( kspace_us, cal_width=cal_width, ndim=2, mask=mask, ret_cal=True, ret_image=True, ) smaps_cal_np = smaps_cal.numpy() if isinstance(smaps_cal, torch.Tensor) else smaps_cal grappa_train_np = ( grappa_train.numpy() if isinstance(grappa_train, torch.Tensor) else grappa_train ) image_cal_np = image_cal.numpy() if isinstance(image_cal, torch.Tensor) else image_cal print(f"\n[Step 2] Smaps shape : {smaps_cal_np.shape}") print(f"[Step 2] Low-res image shape: {image_cal_np.shape}") print(f"[Step 2] Cal k-space shape : {grappa_train_np.shape}") .. rst-class:: sphx-glr-script-out .. code-block:: none [Step 2] Smaps shape : (16, 217, 181) [Step 2] Low-res image shape: (24, 24) [Step 2] Cal k-space shape : (16, 24, 24) .. GENERATED FROM PYTHON SOURCE LINES 314-319 .. code-block:: Python # Low-res image and smaps at cal_width resolution fig, axes = plt.subplots(1, ncols_show + 1, figsize=(3 * (ncols_show + 1), 3)) .. image-sg:: /generated/autoexamples/images/sphx_glr_example04_calibration_004.png :alt: NLINV calibration mode - low-res result (24x24), NLINV image (24x24), NLINV smap 1, NLINV smap 2, NLINV smap 3, NLINV smap 4, RGB key :srcset: /generated/autoexamples/images/sphx_glr_example04_calibration_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 335-340 .. code-block:: Python # Zero-filled vs synthesized central k-space patch acr_zf = kspace_us[ .. image-sg:: /generated/autoexamples/images/sphx_glr_example04_calibration_005.png :alt: Calibration k-space (24x24) - zero-filled (top) vs NLINV-synthesized (bottom), Zero-filled - coil 1, Zero-filled - coil 2, Zero-filled - coil 3, Zero-filled - coil 4, Synthesized - coil 1, Synthesized - coil 2, Synthesized - coil 3, Synthesized - coil 4 :srcset: /generated/autoexamples/images/sphx_glr_example04_calibration_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 365-372 .. note:: In calibration mode (integer ``cal_width``) the sensitivity maps are estimated at low resolution and zero-padded to the full FOV - they are smoother but less accurate near the FOV boundary. Use ``cal_width=None`` when full image reconstruction quality matters; use an integer ``cal_width`` when only the synthesized k-space patch and fast coil estimates are needed (e.g.\ as input to GROG/GRAPPA kernel training). .. GENERATED FROM PYTHON SOURCE LINES 372-376 .. code-block:: Python plt.show() .. GENERATED FROM PYTHON SOURCE LINES 377-385 Multi-slice batched NLINV calibration ===================================== :func:`~pygrog.utils.nlinv_calib` accepts a leading batch axis and runs the calibration once per batch element. Sensitivity maps and image reconstructions are returned per-slice; the synthesized GRAPPA training k-space can optionally be averaged across the batch with ``train_reduce='mean'`` to produce a single shared kernel. .. GENERATED FROM PYTHON SOURCE LINES 385-402 .. code-block:: Python batch_slices = np.stack([kspace_us, kspace_us[:, ::-1, :]], axis=0) print(f"\n[Multi-slice] batched k-space shape : {batch_slices.shape}") # Per-slice smaps + shared (mean-reduced) GRAPPA training k-space. smaps_b, train_b_mean, image_b = nlinv_calib( batch_slices, cal_width=cal_width, ndim=2, ret_cal=True, ret_image=True, train_reduce="mean", ) print(f"[Multi-slice] smaps shape : {tuple(smaps_b.shape)}") print(f"[Multi-slice] image shape : {tuple(image_b.shape)}") print(f"[Multi-slice] mean-train shape : {tuple(train_b_mean.shape)}") .. rst-class:: sphx-glr-script-out .. code-block:: none [Multi-slice] batched k-space shape : (2, 16, 217, 181) [Multi-slice] smaps shape : (2, 16, 217, 181) [Multi-slice] image shape : (2, 24, 24) [Multi-slice] mean-train shape : (16, 24, 24) .. GENERATED FROM PYTHON SOURCE LINES 403-405 .. code-block:: Python fig, axes = plt.subplots(2, ncols_show, figsize=(3 * ncols_show, 5.5)) .. image-sg:: /generated/autoexamples/images/sphx_glr_example04_calibration_006.png :alt: Batched NLINV - per-slice smaps (shared GRAPPA training kernel), slice 0 - smap 1, slice 0 - smap 2, slice 0 - smap 3, slice 0 - smap 4, slice 1 - smap 1, slice 1 - smap 2, slice 1 - smap 3, slice 1 - smap 4, RGB key :srcset: /generated/autoexamples/images/sphx_glr_example04_calibration_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 9.896 seconds) .. _sphx_glr_download_generated_autoexamples_example04_calibration.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: example04_calibration.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: example04_calibration.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: example04_calibration.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_