.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/tuto_06_b_dynamic.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_gallery_tuto_06_b_dynamic.py: 06.b. Dynamic acquisitions + reconstruction =========================================== .. _tuto_06b_dynamic: This tutorial shows how to simulate measurements of dynamic scenes using the :mod:`spyrit.core.meas` submodule. It also demonstrates how to reconstruct a clean image of the scene through motion compensation. It is based on the :class:`spyrit.core.meas.DynamicHadamSplit2d` class of the :mod:`spyrit.core.meas` submodule. .. image:: ../fig/tuto06b_network.png :width: 600 :align: center :alt: Overview of the dynamic pipeline We consider the inverse problem .. math:: y = \text{diag}(A x_{t=1,..., 2M}), where :math:`A \colon\, \mathbb{R}_+^{2M\times N}` is the acquisition matrix that contains positive DMD patterns, :math:`x_{t=1,..., 2M} \in \mathbb{R}^{N \times 2M}` is the temporal signal of interest, :math:`2M` is both the number of DMD patterns (positives and negatives) and the number of frames, :math:`N` is the dimension of the signal within the field of view, :math:`\text{diag}\colon\, \mathbb{R}^{2M \times 2M} \to \mathbb{R}^{2M}` extracts the diagonal of its input. Dynamic single-pixel imaging aims to reconstruct the temporal sequence :math:`x_{t=1,..., 2M}`. Our approach is based on motion-compensation. Assuming known motion during the acquisition, we build a new forward operator :math:`A_{\rm dyn}` that compensates the motion to a static reference frame :math:`x`, i.e. such that .. math:: \text{diag}(A x_{t=1,..., 2M}) = A_{\rm dyn} x. .. important:: The theory behind dynamic single-pixel imaging with motion compensation is detailed in [Maitre2024_1]_, [Maitre2024_2]_, [Maitre2026]_. .. warning:: This tutorial assumes the reader is already familiar with the static measurements operators. If not, we recommend starting with the tutorials :ref:`1a `, :ref:`1b `, and :ref:`1c `. Key concepts: - Simulation of a **dynamic single-pixel acquisition** accounting for motion during the measurement process - Creation of a **dynamic forward operator** :math:`A_{\text{dyn}}` using motion compensation with pattern warping or image warping - Numerical evidence that image warping is better suited for unbiased simulations - Regularized reconstruction with finite differences .. GENERATED FROM PYTHON SOURCE LINES 62-80 .. code-block:: Python import torch import torchvision import matplotlib.pyplot as plt import os import math import time from pathlib import Path from spyrit.misc.disp import torch2numpy from spyrit.misc.statistics import transform_norm import spyrit.misc.metrics as score from spyrit.misc.load_data import download_girder import spyrit.core.torch as spytorch from spyrit.core.prep import Unsplit from spyrit.core.warp import AffineDeformationField .. GENERATED FROM PYTHON SOURCE LINES 81-86 Set acquisition parameters: - img_size: Full image resolution (pixels) - n: Measurement pattern size (defines FOV) - und: Undersampling factor (1 = no undersampling) - M: Total number of measurements per frame .. GENERATED FROM PYTHON SOURCE LINES 86-119 .. code-block:: Python img_size = 88 # Full image side's size in pixels n = 64 # Measurement pattern side's size in pixels (Field of View) und = 1 # Undersampling factor (1 = full sampling) M = n**2 // und # Number of (positive, negative) measurements print(f"Image size: {img_size}x{img_size}") print(f"Measurement FOV: {n}x{n}") print(f"Measurements: {M}") # Dataset parameters i = 0 # Image index (modify to change the image) spyritPath = "../data/data_online/" imgs_path = os.path.join(spyritPath, "spyrit/") # Computation parameters device = torch.device("cuda" if torch.cuda.is_available() else "cpu") print(f"Using device: {device}") dtype = ( torch.float32 ) # Low precision for tutorial stability, feel free to use torch.float64 for better accuracy if memory allows simu_interp = "bilinear" # Interpolation for motion simulation reco_interp = "bilinear" # Interpolation for building forward operator time_dim = 1 # Time dimension index in tensors # Derived parameters meas_shape = (n, n) img_shape = (img_size, img_size) amp_max = (img_shape[0] - meas_shape[0]) // 2 # Border size for centering FOV .. rst-class:: sphx-glr-script-out .. code-block:: none Image size: 88x88 Measurement FOV: 64x64 Measurements: 4096 Using device: cpu .. GENERATED FROM PYTHON SOURCE LINES 120-121 Load an image from Tomoradio's warehouse. .. GENERATED FROM PYTHON SOURCE LINES 121-152 .. code-block:: Python # Download an RGB brain surface image. url_tomoradio = "https://tomoradio-warehouse.creatis.insa-lyon.fr/api/v1" data_root = Path("../data/data_online/2025_dynamic") # local path to data imgs_path = data_root / Path("images/") id_files = [ "69248e3204d23f6e964b16b7", # brain_surface_colorized.png ] try: download_girder(url_tomoradio, id_files, imgs_path) except Exception as e: print("Unable to download from the Tomoradio warehouse") print(e) # Create a transform for natural images to normalized image tensors transform = transform_norm(img_size=img_size) batch_size = 1 # Create dataset and loader dataset = torchvision.datasets.ImageFolder(root=data_root, transform=transform) dataloader = torch.utils.data.DataLoader(dataset, batch_size=batch_size) img, _ = dataloader.dataset[0] x = img.unsqueeze(0).to(dtype=dtype, device=device) print(f"Shape of input images: {x.shape}") x = (x - x.min()) / (x.max() - x.min()) n_wav = x.shape[1] .. rst-class:: sphx-glr-script-out .. code-block:: none File already exists at ../data/data_online/2025_dynamic/images/brain_surface_colorized.png Shape of input images: torch.Size([1, 3, 88, 88]) .. GENERATED FROM PYTHON SOURCE LINES 153-154 Plot the reference image .. GENERATED FROM PYTHON SOURCE LINES 154-165 .. code-block:: Python x_plot = x.moveaxis(1, -1).squeeze().cpu().numpy() plt.imshow(x_plot) if n_wav == 1: plt.colorbar(fraction=0.046, pad=0.04) plt.title("Ground truth") plt.axis("off") plt.show() .. image-sg:: /gallery/images/sphx_glr_tuto_06_b_dynamic_001.png :alt: Ground truth :srcset: /gallery/images/sphx_glr_tuto_06_b_dynamic_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 166-174 Define motion model and deformation fields ########################################## We simulate a pulsating motion using affine transformations (see :ref:`tutorial 6a `). Both forward and inverse deformation fields are needed for the tutorial: - Forward field: for motion simulation & reconstruction with image warping - Inverse field: for reconstruction with pattern warping .. GENERATED FROM PYTHON SOURCE LINES 174-225 .. code-block:: Python a = 0.2 # Scaling amplitude T = 1000 # Period of motion cycle print(f"Motion parameters:") print(f" Amplitude: {a:.2f}") print(f" Period: {T} ms") def s(t): return 1 + a * math.sin(t * 2 * math.pi / T) def f(t): return torch.tensor( [ [1 / s(t), 0, 0], [0, s(t), 0], [0, 0, 1], ], dtype=dtype, ) def f_inv(t): return torch.tensor( [ [s(t), 0, 0], [0, 1 / s(t), 0], [0, 0, 1], ], dtype=dtype, ) # Create time vector for 2*M frames (positive + negative measurements) n_frames = 2 * M time_vector = torch.linspace(0, 2 * T, n_frames) print(f"Time vector: {n_frames} frames from 0 to {2*T}ms") # Create instances of affine deformation fields def_field = AffineDeformationField( f, time_vector, img_shape, dtype=dtype, device=device ) def_field_inv = AffineDeformationField( f_inv, time_vector, img_shape, dtype=dtype, device=device ) print(f"Created deformation fields with shape: {def_field.field.shape}") .. rst-class:: sphx-glr-script-out .. code-block:: none Motion parameters: Amplitude: 0.20 Period: 1000 ms Time vector: 8192 frames from 0 to 2000ms Created deformation fields with shape: torch.Size([8192, 88, 88, 2]) .. GENERATED FROM PYTHON SOURCE LINES 226-227 We apply the deformation field to create a dynamic image sequence .. GENERATED FROM PYTHON SOURCE LINES 227-233 .. code-block:: Python x_motion = def_field(x, 0, n_frames, mode=simu_interp) x_motion = x_motion.moveaxis(time_dim, 1) print(f"Dynamic sequence shape: {x_motion.shape}") print(f"Generated {x_motion.shape[1]} frames for acquisition simulation") .. rst-class:: sphx-glr-script-out .. code-block:: none Dynamic sequence shape: torch.Size([1, 8192, 3, 88, 88]) Generated 8192 frames for acquisition simulation .. GENERATED FROM PYTHON SOURCE LINES 234-235 We display a few frames to visualize the motion. .. GENERATED FROM PYTHON SOURCE LINES 235-262 .. code-block:: Python # plot few frames n_frames_display = 1000 plt.figure(figsize=(12, 3)) n_rows, n_cols = 1, 4 for frame in range(n_frames): n_frame = n_frames_display * frame if n_frame >= n_frames or frame >= n_rows * n_cols: break plt.subplot(n_rows, n_cols, frame + 1) plt.imshow( x_motion[ 0, n_frame, :, amp_max : img_size - amp_max, amp_max : img_size - amp_max ] .moveaxis(0, -1) .view(*meas_shape, n_wav) .cpu() .numpy(), cmap="gray", ) # in X plt.title("frame %d" % (n_frame), fontsize=12) plt.axis("off") plt.tight_layout() plt.show() .. image-sg:: /gallery/images/sphx_glr_tuto_06_b_dynamic_002.png :alt: frame 0, frame 1000, frame 2000, frame 3000 :srcset: /gallery/images/sphx_glr_tuto_06_b_dynamic_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 263-277 Dynamic measurement simulation ############################## We simulate a dynamic single-pixel acquisition. .. note:: The dynamic measurement operator classes extend the static ones: - :class:`Linear` -> :class:`DynamicLinear` - :class:`LinearSplit` -> :class:`DynamicLinearSplit` - :class:`HadamSplit2d` -> :class:`DynamicHadamSplit2d` We used the :class:`DynamicHadamSplit2d` class (specialized for Hadamard patterns) for this tutorial. .. GENERATED FROM PYTHON SOURCE LINES 277-303 .. code-block:: Python from spyrit.core.meas import DynamicHadamSplit2d meas_op = DynamicHadamSplit2d( time_dim=time_dim, h=n, M=M, order=None, fast=True, reshape_output=False, img_shape=img_shape, noise_model=torch.nn.Identity(), white_acq=None, dtype=dtype, device=device, ) t1 = time.time() y1 = meas_op(x_motion) t2 = time.time() print(f"Computation time: {t2 - t1:.3f}s") print(f"Output shape: {y1.shape}") .. rst-class:: sphx-glr-script-out .. code-block:: none Computation time: 0.601s Output shape: torch.Size([1, 3, 8192]) .. GENERATED FROM PYTHON SOURCE LINES 304-306 We preprocess measurements for reconstruction using the differential strategy to combine the positive/negative measurements. .. GENERATED FROM PYTHON SOURCE LINES 306-315 .. code-block:: Python prep_op = Unsplit() y2 = prep_op(y1) print(f"\nMeasurement processing:") print(f" Raw measurements shape: {y1.shape}") print(f" Processed measurements shape: {y2.shape}") .. rst-class:: sphx-glr-script-out .. code-block:: none Measurement processing: Raw measurements shape: torch.Size([1, 3, 8192]) Processed measurements shape: torch.Size([1, 3, 4096]) .. GENERATED FROM PYTHON SOURCE LINES 316-321 Static reconstruction baseline ############################## Compute a static reconstruction for comparison. This ignores motion and treats all measurements as from a static scene. .. GENERATED FROM PYTHON SOURCE LINES 321-342 .. code-block:: Python from spyrit.core.meas import HadamSplit2d meas_op_stat = HadamSplit2d(h=n, M=M, order=None, dtype=dtype, device=device) print(f"\n=== Static Reconstruction (Baseline) ===") x_stat = meas_op_stat.fast_pinv(y2) print(f"Static reconstruction shape: {x_stat.shape}") # Quick quality check for static reconstruction x_ref_fov = x[0, :, amp_max : n + amp_max, amp_max : n + amp_max] static_psnr = score.psnr(torch2numpy(x_stat), torch2numpy(x_ref_fov)) print(f"Static reconstruction PSNR: {static_psnr:.2f} dB") plt.figure(figsize=(6, 6)) plt.imshow(torch2numpy(x_stat.moveaxis(1, -1).squeeze())) plt.title("Static Reconstruction \n (strong artifacts due to motion)") plt.axis("off") plt.show() .. image-sg:: /gallery/images/sphx_glr_tuto_06_b_dynamic_003.png :alt: Static Reconstruction (strong artifacts due to motion) :srcset: /gallery/images/sphx_glr_tuto_06_b_dynamic_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none === Static Reconstruction (Baseline) === Static reconstruction shape: torch.Size([1, 3, 64, 64]) Static reconstruction PSNR: 19.86 dB Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Got range [-0.21864055..1.0663859]. .. GENERATED FROM PYTHON SOURCE LINES 343-345 We set up regularization for the inverse problem with first-order finite differences (smoothing prior). We use Neumann boundary conditions. .. GENERATED FROM PYTHON SOURCE LINES 345-352 .. code-block:: Python Dx, Dy = spytorch.neumann_boundary(img_shape) D2 = Dx.T @ Dx + Dy.T @ Dy D2 = D2.to(dtype=dtype) print(f"Regularization matrix shape: {D2.shape}") .. rst-class:: sphx-glr-script-out .. code-block:: none Regularization matrix shape: torch.Size([7744, 7744]) .. GENERATED FROM PYTHON SOURCE LINES 353-363 Build dynamic system matrix ########################### Construct the dynamic forward operator :math:`H_{\rm dyn}` that accounts for motion. This matrix maps the reference image to measurements accounting for motion. Two approaches: 1. :attr:`warping='pattern'`: Warp the patterns. Need to pass the inverse deformation field. 2. :attr:`warping='image'`: Warp the image. Need to pass the forward deformation field. .. GENERATED FROM PYTHON SOURCE LINES 365-367 1. Pattern warping ------------------ .. GENERATED FROM PYTHON SOURCE LINES 367-379 .. code-block:: Python # Build the dynamic system matrix print("Building H_dyn using pattern warping (inverse deformation)...") meas_op.build_dynamic_forward( def_field_inv, warping="pattern", mode=reco_interp, verbose=False ) print(f"Dynamic system matrix shape: {meas_op.A_dyn.shape}") print( f"Dynamic system matrix (differential strategy AFTER motion compensation) shape: {meas_op.H_dyn.shape}" ) .. rst-class:: sphx-glr-script-out .. code-block:: none Building H_dyn using pattern warping (inverse deformation)... Be careful to use the inverse deformation field when warping patterns. Dynamic system matrix shape: torch.Size([8192, 7744]) Dynamic system matrix (differential strategy AFTER motion compensation) shape: torch.Size([4096, 7744]) .. GENERATED FROM PYTHON SOURCE LINES 380-384 .. note:: In order to use the differential strategy with dynamic scenes, it is important to build the dynamic system matrix after motion compensation [Maitre2026]_. .. GENERATED FROM PYTHON SOURCE LINES 387-393 Verify forward model accuracy ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Test the dynamic forward model by computing the residual of the forward model. Without measurement noise and using dtype=float64, the residual should be very small (:math:`\approx` 1e-12) [Maitre2024_2]_. When using the pattern warping approach, the forward model is non-zero, indicating a bias in the model. .. GENERATED FROM PYTHON SOURCE LINES 393-413 .. code-block:: Python A_dyn_x = meas_op.forward_A_dyn(x) residual_norm = torch.norm(y1 - A_dyn_x).item() relative_error = residual_norm / torch.norm(y1).item() print(f"\n=== Forward Model Verification ===") print(f" Predicted measurements shape: {A_dyn_x.shape}") print(f" Residual norm: {residual_norm:.2e}") print(f" Relative error: {relative_error:.2e}") # Visualize residual pattern (averaged over spectral channels) plt.figure(figsize=(4, 6)) residual_2d = abs(y1 - A_dyn_x).mean(dim=1).squeeze().cpu().numpy().reshape((2 * n, n)) plt.imshow(residual_2d, cmap="Spectral") plt.colorbar(fraction=0.046 * 6 / 3, pad=0.04) plt.title(f"Forward Model Residual |y - A_dyn·x| \n Max: {residual_2d.max():.2e}") plt.tight_layout() plt.show() .. image-sg:: /gallery/images/sphx_glr_tuto_06_b_dynamic_004.png :alt: Forward Model Residual |y - A_dyn·x| Max: 2.09e+01 :srcset: /gallery/images/sphx_glr_tuto_06_b_dynamic_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none === Forward Model Verification === Predicted measurements shape: torch.Size([1, 3, 8192]) Residual norm: 2.84e+02 Relative error: 1.78e-03 .. GENERATED FROM PYTHON SOURCE LINES 414-415 We display few dynamic patterns. .. GENERATED FROM PYTHON SOURCE LINES 415-438 .. code-block:: Python print(f"\nVisualizing dynamic matrix evolution...") H_dyn_diff_np = torch2numpy(meas_op.H_dyn) # plot few patterns n_frames_display = 500 plt.figure(figsize=(12, 3)) n_rows, n_cols = 1, 4 for frame in range(n_frames): n_frame = n_frames_display * frame if n_frame >= n_frames or frame >= n_rows * n_cols: break plt.subplot(n_rows, n_cols, frame + 1) plt.imshow( H_dyn_diff_np[n_frame].reshape(img_shape), cmap="gray", vmin=-1, vmax=1 ) # in X_{ext} plt.title("frame %d" % (2 * n_frame), fontsize=12) plt.axis("off") plt.tight_layout() plt.show() .. image-sg:: /gallery/images/sphx_glr_tuto_06_b_dynamic_005.png :alt: frame 0, frame 1000, frame 2000, frame 3000 :srcset: /gallery/images/sphx_glr_tuto_06_b_dynamic_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Visualizing dynamic matrix evolution... .. GENERATED FROM PYTHON SOURCE LINES 439-443 Analyze system conditioning ~~~~~~~~~~~~~~~~~~~~~~~~~~~ Move computations to CPU for optimized linear algebra. .. GENERATED FROM PYTHON SOURCE LINES 443-449 .. code-block:: Python print(f"\n=== Preparing for Reconstruction ===") print("Moving to CPU for optimized linear algebra...") H_dyn = meas_op.H_dyn.cpu() y2 = y2.cpu() .. rst-class:: sphx-glr-script-out .. code-block:: none === Preparing for Reconstruction === Moving to CPU for optimized linear algebra... .. GENERATED FROM PYTHON SOURCE LINES 450-452 Compute singular values to understand the inverse problem difficulty. High condition number indicates need for regularization. .. GENERATED FROM PYTHON SOURCE LINES 452-465 .. code-block:: Python print("Analyzing system matrix conditioning...") sing_vals = torch.linalg.svdvals(H_dyn) condition_number = (sing_vals[0] / sing_vals[-1]).item() sigma_max = sing_vals[0].item() sigma_min = sing_vals[-1].item() print(f"Singular value spectrum:") print(f" Maximum: {sigma_max:.2e}") print(f" Minimum: {sigma_min:.2e}") print(f" Condition number: {condition_number:.2e}") .. rst-class:: sphx-glr-script-out .. code-block:: none Analyzing system matrix conditioning... Singular value spectrum: Maximum: 1.11e+02 Minimum: 9.78e-02 Condition number: 1.13e+03 .. GENERATED FROM PYTHON SOURCE LINES 466-474 Solve the regularized least squares problem ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. math:: \arg \min_x \frac{1}{2} \| y - H_{\rm dyn}^{\rm wp} x \|^2 + \frac{\tilde{\eta}}{2} \| D x \|^2 where :math:`\tilde{\eta}` is the normalized regularization parameter scaled by the maximum singular value of :math:`H_{\rm dyn}` and :math:`D` is a first order finite difference operator. .. GENERATED FROM PYTHON SOURCE LINES 474-489 .. code-block:: Python eta = 1e-5 # Regularization parameter (adjust based on noise level) print(f"\n=== Dynamic Reconstruction ===") print(f"Regularization parameter: {eta:.1e}") start_time = time.time() x_dyn_wp = torch.linalg.solve( H_dyn.T @ H_dyn + eta * sigma_max**2 * D2, H_dyn.T @ y2.moveaxis(1, -1) ) solve_time = time.time() - start_time print(f"Reconstruction completed in {solve_time:.2f}s") print(f"Solution shape: {x_dyn_wp.shape}") .. rst-class:: sphx-glr-script-out .. code-block:: none === Dynamic Reconstruction === Regularization parameter: 1.0e-05 Reconstruction completed in 9.33s Solution shape: torch.Size([1, 7744, 3]) .. GENERATED FROM PYTHON SOURCE LINES 490-491 We plot the dynamic reconstruction .. GENERATED FROM PYTHON SOURCE LINES 491-504 .. code-block:: Python x_dyn_wp_plot = x_dyn_wp.view(*img_shape, n_wav) plt.figure(figsize=(6, 6)) plt.imshow(x_dyn_wp_plot.cpu().numpy(), cmap="gray") plt.title( "Dynamic Reconstruction \n (Motion Compensated with pattern warping)", fontsize=14 ) plt.axis("off") if n_wav == 1: plt.colorbar(fraction=0.046, pad=0.04) .. image-sg:: /gallery/images/sphx_glr_tuto_06_b_dynamic_006.png :alt: Dynamic Reconstruction (Motion Compensated with pattern warping) :srcset: /gallery/images/sphx_glr_tuto_06_b_dynamic_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Got range [-0.39748687..1.2340564]. .. GENERATED FROM PYTHON SOURCE LINES 505-507 2. Image warping ---------------- .. GENERATED FROM PYTHON SOURCE LINES 507-519 .. code-block:: Python # Build the dynamic system matrix print("Building H_dyn using image warping (forward deformation)...") meas_op.build_dynamic_forward( def_field, warping="image", mode=reco_interp, verbose=False ) print( f"Dynamic system matrix (differential strategy AFTER motion compensation) shape: {meas_op.H_dyn.shape}" ) .. rst-class:: sphx-glr-script-out .. code-block:: none Building H_dyn using image warping (forward deformation)... Dynamic system matrix (differential strategy AFTER motion compensation) shape: torch.Size([4096, 7744]) .. GENERATED FROM PYTHON SOURCE LINES 520-526 Verify forward model accuracy ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Test the dynamic forward model by computing the residual of the forward model. Without measurement noise and using dtype=float64, the residual should be very small (:math:`\approx` 1e-12) [Maitre2024_2]_. When using the image warping approach, the forward model is almost zero, indicating no bias in the model. .. GENERATED FROM PYTHON SOURCE LINES 526-546 .. code-block:: Python A_dyn_x = meas_op.forward_A_dyn(x) residual_norm = torch.norm(y1 - A_dyn_x).item() relative_error = residual_norm / torch.norm(y1).item() print(f"\n=== Forward Model Verification ===") print(f" Predicted measurements shape: {A_dyn_x.shape}") print(f" Residual norm: {residual_norm:.2e}") print(f" Relative error: {relative_error:.2e}") # Visualize residual pattern (averaged over spectral channels) plt.figure(figsize=(4, 6)) residual_2d = abs(y1 - A_dyn_x).mean(dim=1).squeeze().cpu().numpy().reshape((2 * n, n)) plt.imshow(residual_2d, cmap="Spectral") plt.colorbar(fraction=0.046 * 6 / 3, pad=0.04) plt.title(f"Forward Model Residual |y - A_dyn·x| \n Max: {residual_2d.max():.2e}") plt.tight_layout() plt.show() .. image-sg:: /gallery/images/sphx_glr_tuto_06_b_dynamic_007.png :alt: Forward Model Residual |y - A_dyn·x| Max: 1.59e-03 :srcset: /gallery/images/sphx_glr_tuto_06_b_dynamic_007.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none === Forward Model Verification === Predicted measurements shape: torch.Size([1, 3, 8192]) Residual norm: 4.04e-02 Relative error: 2.53e-07 .. GENERATED FROM PYTHON SOURCE LINES 547-548 We display few dynamic patterns .. GENERATED FROM PYTHON SOURCE LINES 548-571 .. code-block:: Python print(f"\nVisualizing dynamic matrix evolution...") H_dyn_diff_np = torch2numpy(meas_op.H_dyn) # plot few patterns n_frames_display = 500 plt.figure(figsize=(12, 3)) n_rows, n_cols = 1, 4 for frame in range(n_frames): n_frame = n_frames_display * frame if n_frame >= n_frames or frame >= n_rows * n_cols: break plt.subplot(n_rows, n_cols, frame + 1) plt.imshow( H_dyn_diff_np[n_frame].reshape(img_shape), cmap="gray", vmin=-1, vmax=1 ) # in X_{ext} plt.title("frame %d" % (2 * n_frame), fontsize=12) plt.axis("off") plt.tight_layout() plt.show() .. image-sg:: /gallery/images/sphx_glr_tuto_06_b_dynamic_008.png :alt: frame 0, frame 1000, frame 2000, frame 3000 :srcset: /gallery/images/sphx_glr_tuto_06_b_dynamic_008.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Visualizing dynamic matrix evolution... .. GENERATED FROM PYTHON SOURCE LINES 572-576 Analyze system conditioning ~~~~~~~~~~~~~~~~~~~~~~~~~~~ Move computations to CPU for optimized linear algebra. .. GENERATED FROM PYTHON SOURCE LINES 576-582 .. code-block:: Python print(f"\n=== Preparing for Reconstruction ===") print("Moving to CPU for optimized linear algebra...") H_dyn = meas_op.H_dyn.cpu() y2 = y2.cpu() .. rst-class:: sphx-glr-script-out .. code-block:: none === Preparing for Reconstruction === Moving to CPU for optimized linear algebra... .. GENERATED FROM PYTHON SOURCE LINES 583-585 Compute singular values to understand the inverse problem difficulty. High condition number indicates need for regularization. .. GENERATED FROM PYTHON SOURCE LINES 585-598 .. code-block:: Python print("Analyzing system matrix conditioning...") sing_vals = torch.linalg.svdvals(H_dyn) condition_number = (sing_vals[0] / sing_vals[-1]).item() sigma_max = sing_vals[0].item() sigma_min = sing_vals[-1].item() print(f"Singular value spectrum:") print(f" Maximum: {sigma_max:.2e}") print(f" Minimum: {sigma_min:.2e}") print(f" Condition number: {condition_number:.2e}") .. rst-class:: sphx-glr-script-out .. code-block:: none Analyzing system matrix conditioning... Singular value spectrum: Maximum: 1.10e+02 Minimum: 2.39e-01 Condition number: 4.61e+02 .. GENERATED FROM PYTHON SOURCE LINES 599-607 Solve the regularized least squares problem ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. math:: \arg \min_x \frac{1}{2} \| y - H_{\rm dyn}^{\rm wf} x \|^2 + \frac{\tilde{\eta}}{2} \| D x \|^2 where :math:`\tilde{\eta}` is the normalized regularization parameter scaled by the maximum singular value of :math:`H_{\rm dyn}` and :math:`D` is a first order finite difference operator. .. GENERATED FROM PYTHON SOURCE LINES 607-622 .. code-block:: Python eta = 1e-5 # Regularization parameter (adjust based on noise level) print(f"\n=== Dynamic Reconstruction ===") print(f"Regularization parameter: {eta:.1e}") start_time = time.time() x_dyn_wf = torch.linalg.solve( H_dyn.T @ H_dyn + eta * sigma_max**2 * D2, H_dyn.T @ y2.moveaxis(1, -1) ) solve_time = time.time() - start_time print(f"Reconstruction completed in {solve_time:.2f}s") print(f"Solution shape: {x_dyn_wf.shape}") .. rst-class:: sphx-glr-script-out .. code-block:: none === Dynamic Reconstruction === Regularization parameter: 1.0e-05 Reconstruction completed in 9.38s Solution shape: torch.Size([1, 7744, 3]) .. GENERATED FROM PYTHON SOURCE LINES 623-624 We plot the dynamic reconstruction .. GENERATED FROM PYTHON SOURCE LINES 624-637 .. code-block:: Python x_dyn_wf_plot = x_dyn_wf.view(*img_shape, n_wav) plt.figure(figsize=(6, 6)) plt.imshow(x_dyn_wf_plot.cpu().numpy(), cmap="gray") plt.title( "Dynamic Reconstruction \n (Motion Compensated with image warping)", fontsize=14 ) plt.axis("off") if n_wav == 1: plt.colorbar(fraction=0.046, pad=0.04) .. image-sg:: /gallery/images/sphx_glr_tuto_06_b_dynamic_009.png :alt: Dynamic Reconstruction (Motion Compensated with image warping) :srcset: /gallery/images/sphx_glr_tuto_06_b_dynamic_009.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Got range [-0.5651888..1.5256104]. .. GENERATED FROM PYTHON SOURCE LINES 638-645 Results overview ################ Compare the original image, static reconstruction, and dynamic reconstruction. This shows the improvement gained by accounting for motion. sphinx_gallery_thumbnail_number = 10 .. GENERATED FROM PYTHON SOURCE LINES 645-687 .. code-block:: Python print(f"\n=== Reconstruction Comparison ===") fig, ax = plt.subplots(2, 2, figsize=(10, 10)) # Original reference image ref_img = x.view(n_wav, *img_shape).moveaxis(0, -1).cpu().numpy() im0 = ax[0, 0].imshow(ref_img, cmap="gray") ax[0, 0].set_title("Reference Image", fontsize=18) ax[0, 0].axis("off") if n_wav == 1: fig.colorbar(im0, ax=ax[0, 0], fraction=0.046, pad=0.04) # Static reconstruction static_img = torch.zeros((*img_shape, n_wav), dtype=dtype, device=device) static_img[amp_max : n + amp_max, amp_max : n + amp_max] = x_stat.view( n_wav, *meas_shape ).moveaxis(0, -1) im1 = ax[0, 1].imshow(static_img.cpu().numpy(), cmap="gray") ax[0, 1].set_title("Static Reconstruction \n (Ignores Motion)", fontsize=18) ax[0, 1].axis("off") if n_wav == 1: fig.colorbar(im1, ax=ax[0, 1], fraction=0.046, pad=0.04) # Dynamic reconstruction with pattern warping im2 = ax[1, 0].imshow(x_dyn_wp_plot.cpu().numpy(), cmap="gray") ax[1, 0].set_title("Dynamic Reconstruction \n (pattern warping)", fontsize=18) ax[1, 0].axis("off") if n_wav == 1: fig.colorbar(im2, ax=ax[1, 0], fraction=0.046, pad=0.04) # Dynamic reconstruction with image warping im3 = ax[1, 1].imshow(x_dyn_wf_plot.cpu().numpy(), cmap="gray") ax[1, 1].set_title("Dynamic Reconstruction \n (image warping)", fontsize=18) ax[1, 1].axis("off") if n_wav == 1: fig.colorbar(im3, ax=ax[1, 1], fraction=0.046, pad=0.04) plt.tight_layout() plt.show() .. image-sg:: /gallery/images/sphx_glr_tuto_06_b_dynamic_010.png :alt: Reference Image, Static Reconstruction (Ignores Motion), Dynamic Reconstruction (pattern warping), Dynamic Reconstruction (image warping) :srcset: /gallery/images/sphx_glr_tuto_06_b_dynamic_010.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none === Reconstruction Comparison === Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Got range [-0.21864055..1.0663859]. Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Got range [-0.39748687..1.2340564]. Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Got range [-0.5651888..1.5256104]. .. GENERATED FROM PYTHON SOURCE LINES 688-690 Compute image quality metrics to quantify reconstruction improvement. Metrics are computed in the FOV region for fair comparison. .. GENERATED FROM PYTHON SOURCE LINES 690-734 .. code-block:: Python # Extract FOV region for metric calculation x_dyn_wp_in_X = x_dyn_wp_plot[amp_max : n + amp_max, amp_max : n + amp_max] x_dyn_wf_in_X = x_dyn_wf_plot[amp_max : n + amp_max, amp_max : n + amp_max] x_ref_in_X = x[0, :, amp_max : n + amp_max, amp_max : n + amp_max].moveaxis(0, -1) # Calculate metrics for both reconstructions psnr_static = score.psnr( torch2numpy(x_stat.view(n_wav, *meas_shape).moveaxis(0, -1)), torch2numpy(x_ref_in_X), ) ssim_static = score.ssim( torch2numpy(x_stat.view(n_wav, *meas_shape).moveaxis(0, -1)), torch2numpy(x_ref_in_X), ) psnr_dynamic_wp = score.psnr(torch2numpy(x_dyn_wp_in_X), torch2numpy(x_ref_in_X)) ssim_dynamic_wp = score.ssim(torch2numpy(x_dyn_wp_in_X), torch2numpy(x_ref_in_X)) psnr_dynamic_wf = score.psnr(torch2numpy(x_dyn_wf_in_X), torch2numpy(x_ref_in_X)) ssim_dynamic_wf = score.ssim(torch2numpy(x_dyn_wf_in_X), torch2numpy(x_ref_in_X)) print(f"\n=== Quantitative Results (in the SPC FOV X) ===") print(f"{'Method':<25} {'PSNR (dB)':<12} {'SSIM'}") print("-" * 48) print(f"{'Static':<25} {psnr_static:<12.2f} {ssim_static:.3f}") print( f"{'Dynamic (pattern warping)':<25} {psnr_dynamic_wp:<12.2f} {ssim_dynamic_wp:.3f}" ) print(f"{'Dynamic (image warping)':<25} {psnr_dynamic_wf:<12.2f} {ssim_dynamic_wf:.3f}") print("-" * 48 + "\n") improvement_wp = psnr_dynamic_wp - psnr_static improvement_wf = psnr_dynamic_wf - psnr_static print( f"Dynamic reconstruction (pattern warping) achieved a PSNR improvement of {improvement_wp:.2f} dB over static reconstruction." ) print( f"Dynamic reconstruction (image warping) achieved a PSNR improvement of {improvement_wf:.2f} dB over static reconstruction." ) .. rst-class:: sphx-glr-script-out .. code-block:: none === Quantitative Results (in the SPC FOV X) === Method PSNR (dB) SSIM ------------------------------------------------ Static 19.86 0.835 Dynamic (pattern warping) 25.16 0.943 Dynamic (image warping) 34.30 0.996 ------------------------------------------------ Dynamic reconstruction (pattern warping) achieved a PSNR improvement of 5.30 dB over static reconstruction. Dynamic reconstruction (image warping) achieved a PSNR improvement of 14.44 dB over static reconstruction. .. rst-class:: sphx-glr-timing **Total running time of the script:** (1 minutes 50.272 seconds) .. _sphx_glr_download_gallery_tuto_06_b_dynamic.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: tuto_06_b_dynamic.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: tuto_06_b_dynamic.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: tuto_06_b_dynamic.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_