Accelerated susceptibility-driven positive contrast MRI reconstruction based on primal-dual optimization with minimal parameter tuning

Introduction

MR compatible metallic devices, such as brachytherapy seeds, stents and so on, are widely utilized in clinical settings (1-4). However, these interventional metallic implants often appear as signal voids (dark holes) in conventional magnetic resonance (MR) images due to susceptibility-induced artifacts (5-7). To overcome this limitation, a susceptibility-driven positive contrast magnetic resonance imaging (PC-MRI) technique has been developed to estimate arbitrary magnetic susceptibility distributions of metallic devices, thereby facilitating high-contrast visualization of these implants.

The reconstruction problem associated with this technique is inherently a non-smooth convex optimization task characterized by a known saddle-point structure, which is typically addressed using the nonlinear conjugate gradient (CG) algorithm (8-10). This algorithm employs a sequential iterative optimization approach to approximate the locally optimal solution with adjustable relaxation parameters. However, the traditional CG method encounters significant limitations in the context of PC-MRI. It requires smoothing the non-differentiable ℓ₁-norm using a conditioning parameter (ϵ), which alters the original problem. This results in intrinsic drawbacks, including sensitivity to initial conditions and parameter selection (particularly ϵ), slow convergence rates when high accuracy is desired, and difficulties in achieving optimal imaging reconstruction due to the ill-posed nature of the inversion (11). Consequently, there is an urgent need for an optimization framework capable of effectively and robustly addressing the specific challenges of PC-MRI reconstruction.

As we know, the first-order primal-dual (PD) algorithm, also known as the Chambolle-Pock algorithm, is a well-established framework for a class of structured convex optimization problems and offers several notable advantages, including the ability to handle non-differentiable objective functions, ease of implementation, and optimal convergence rates across multiple sub-classes (12,13). This algorithm has been employed in various applications, such as, reconstructing general MR images (14), computing generalized total variation in Bayesian MRI (15), and generating images in quantitative susceptibility mapping (QSM) (16,17). Given the reconstruction challenges inherent in PC-MRI, the first-order PD algorithm holds potential for application in this context.

So, leveraging the advantages of the first-order PD method could effectively address those issues of ill-posed inversion and parameters selection in the current PC-MRI methodologies (8,9,18-20). We propose a susceptibility-driven PC-MRI technique that leverages the PD framework to solve the exact ℓ₁-minimization problem without smoothing approximations. The proposed first-order PD algorithm exhibits simplified parameter tuning for solving ℓ₁-minimization problems in optimal reconstruction. Simulations, phantom data and patient experiments validated the PD method’s performance in susceptibility-driven positive contrast imaging using central-difference gradient operators, and compare this with the conventional CG method. The results have shown that the PD method can improve the reconstruction image quality with less reconstruction time than before.

However, PD method involves extensive matrix-vector operations during the image reconstruction. Sequential processing based on a central processing unit (CPU) fails to exploit the inherent parallelism in these computations, resulting in low efficiency. Furthermore, CPU-based methods require setting up a cycle to reconstruct one image per cycle, which further prolongs the computation time (21). These limitations lead to reconstruction times of several minutes for the PD method, even though it significantly outperforms conventional CG methods in speed (13). Therefore, we further accelerate PD reconstruction using graphics processing unit (GPU) parallel computing. GPUs leverage the parallelism in matrix-vector operations and enable multi-slice data reconstruction simultaneously (22,23). Under the same experimental conditions, this approach achieves comparable image quality while drastically reducing reconstruction times.

Theory

Susceptibility-driven positive contrast imaging

The susceptibility-driven PC method utilizes a modified two-dimensional (2D) fast spin echo (FSE) sequence to expedite data acquisition (8). The technique employs an equivalent short echo time (TE) by introducing a gradient shift on the readout gradient (T_shift, typically 0.2–0.7 ms) during MR data acquisition. Thus, a phase difference ∆Φ can be obtained through the phase difference between Φ(T_shift) and Φ(0), where Φ(T_shift) and Φ(0) are respectively obtained by the modified FSE sequence with T_shift and without T_shift (8). The local field variation ΔB(r) caused by MR-compatible metallic devices can be determined through, ∆B(r)=∆Φ/γB₀T_shift. Once ∆B(r) is determined, the susceptibility map χ(r), which serves as a positive-contrast MR image, is reconstructed by solving an ℓ1-norm optimization problem as shown in the following Eq. [1]. This approach leverages the sparsity of χ(r) in the spatial image gradient domain to enable efficient and accurate susceptibility mapping.

χ=argminχ||W(Dχ−ΔB)||22+λ||MGχ||1[1]

In this formulation, D denotes the dipole kernel convolution operator, while W represents the weighting matrix that normalizes the magnitude image. The spatial gradient operator G=[Gx,Gy,Gz]T obtains susceptibility gradients across all three dimensions. To preserve fine structural details in susceptibility distributions, M is a masking matrix that mitigates excessive smoothing of small-scale structures. The regularization parameter λ adjusts the balance between data consistency and the desired smoothness of the reconstructed susceptibility map. To solve the optimization problem of Eq. [1], the conventional nonlinear CG algorithm and the first-order PD algorithm can be utilized to derive the final susceptibility mapping χ, and thereby the PC images of the metallic devices are subsequently derived from χ.

Nonliear CG algorithm

In the conventional nonlinear CG algorithm, the ℓ1-norm ‖x‖1=∑i=1n|xi| is smoothly approximated by ∑i=1nxi2+ϵ for computational purpose, ϵ is small and positive. Then, the Euler-Lagrange equation corresponding to the energy functional can be derived as follows:

λ(WD)∗(W(Dχ−ΔB))+(MG)∗MGχ(MGχ)2+ϵ=0[2]

where ∗ denotes the adjoint operation. As a consequence, χ is solved using the conventional CG with a slightly different objective that involves the condition parameter ϵ.

First-order PD algorithm

The first-order PD algorithm is a class of iterative optimization methods designed to solve the non-smooth convex optimization problems, particularly those that can be formulated as saddle-point problems. It has been proven that convergence exists in mathematics (13). The algorithm operates by iteratively updating the primal and its dual variables to converge to the optimal solution. The PD model is obtained as:

minx∈Xmaxy∈Y{⟨Kx,y⟩+g(x)−f∗(y)}[3]

Here, our reconstruction problem of PC-MRI is not explicitly covered in (13). Inspired by (24), the overall objective function for our reconstruction problem can be expressed as f(Kx). Since both terms in Eq. [1] involve linear transformations, the objective function for our reconstruction problem can be formulated as follows:

f(y,z)=‖y−n‖22+λ‖z‖1andg(χ)=0[4]

While

y=WDχ,n=WDΔB,z=MGχK=(WDMG)[5]

Using the Legendre transform of (15), we obtain the convex conjugate of f:

f1*(y)=12λ‖y‖22+⟨y,n⟩,f2*(z)=δBOX(1)(z)[6]

Therefore, the minimization of Eq. [1] can be written to its dual formulation:

maxy,z{−12λ‖y‖22−⟨y,n⟩−δBOX(1)(z)−δ0S(−D*W*y−G*M*z)}[7]

and

g*(χ)=δ0S(χ)={0χ∈S+∞otherwiseδBOX(a)(η)={0‖η‖∞≤a∞‖η‖∞>a[8]

While ‖·‖∞ norm corresponds to the maximum absolute value of the variable, under the premise that at least one feasible solution exists satisfying the Eq. [1] and Eq. [7], the proximal operator is leveraged to compute the optimal susceptibility mapping by iteratively determining descent directions for the convex objective function.

proxσ[H](z)=argminz′{H(z′)+‖z−z′‖222σ}[9]

The detailed pseudocode for both the conventional CG algorithm and the proposed PD algorithm was provided in Appendix 1 for clarity and reproducibility.

GPU-accelerated PD algorithm for positive contrast reconstruction (GPU-based PD)

The PD algorithm for positive contrast reconstruction traditionally relies on CPU-based computation. However, CPU-driven iterative methods fail to fully exploit the parallelism inherent in vector/matrix operations, resulting in the reconstruction times ranging from tens of seconds to several minutes. To accelerate reconstruction, we implemented the PD algorithm on GPUs, leveraging their massive parallelism to significantly reduce processing time without compromising image quality. The PD algorithm is inherently well-suited for GPU acceleration due to its iterative matrix-vector operations, which can be efficiently parallelized across multiple cores.

Specifically, the iterative components of the PD algorithm were implemented on a GPU during the reconstruction of the positive contrast image. Because the data across different slices exhibit no inter-dependency in the PC-MR image reconstruction, the GPU approach employs a dynamic parallelism-based multi-data concurrent computing strategy (25,26) based on this characteristic, and executes computations across multiple slice-layer data simultaneously. This approach fundamentally eliminates the single-layer loop in CPU-based reconstruction methods across all operational scenarios. Thus, the complex matrix-vector calculation can be efficiently processed through the dynamic parallelism-based multi-data concurrent computing module. For simple vector-matrix calculations, GPU thread invocation is straightforward. The Fourier Transform is computed using the GPU-accelerated computation model, which is provided by the CUDA library. Compared with the CPU-based implementation, the GPU-accelerated approach significantly improves reconstruction efficiency and reduces time consumption by minimizing both data slice loops and the vector-matrix operation loops. The detailed processing workflow is illustrated in Figure 1.

Figure 1 GPU-accelerated PD framework for positive contrast MR image reconstruction. CUDA, Compute Unified Device Architecture; FFT, fast Fourier transform; GPU, graphics processing unit; PD, primal-dual.

Methods

Data set descriptions

In the simulated experiments, synthetic MRI data with a known susceptibility map were generated through computer simulations and subsequently processed using the susceptibility-driven PC-MRI. To create the simulated images, a high-resolution susceptibility map containing a single seed was created as a reference (referred to as the true image), based on the configuration of brachytherapy seeds. For the single-seed simulations, the seed was oriented perpendicular to the B₀ field.

To further evaluate the imaging performance and convergence behavior of the two iterative algorithms in different metal devices, a phantom study was conducted using two geometrically distinct interventional devices with differing dimensions and morphologies. All of the experimental datasets were acquired on a 3 Tesla SIEMENS Tim Trio MRI system (SIEMENS AG, Erlangen, Germany) equipped with an eight-channel phased-array receiver coil. In both phantom experiments, data were collected using the modified FSE sequence, with and without a T_shift applied to the readout gradient (8). A turbo factor (L) of 7 was chosen to balance signal-to-noise ratio (SNR) and acquisition speed, as L was constrained by the T2 relaxation time of the imaged object.

The first phantom experiment was designed to evaluate the two algorithms for imaging brachytherapy seeds. The phantom was constructed following the method described in (8), with five dummy brachytherapy seeds (Seeds Biological Pharmacy Ltd., Tianjin, China) placed in a porcine tissue. Similarly, a plastic stick was inserted into the tissue to simulate a cavity, and a bamboo toothpick was included to mimic a capillary within human tissue. Additionally, a small animal bone was incorporated into the phantom to replicate a human bone structure. In this experiment, the seeds were oriented perpendicular to the B₀ magnetic field. Coronal image slices were acquired orthogonal to the longitudinal axis of the seed (8). The sequence parameters for the first tissue phantom were also the same as (8) as follows: field of view (FOV) =120×120×15 mm³, repetition time (TR) =2,000 ms, TE =18 ms, matrix size =192×192×10, slice thickness =1.5 mm (no slice gap), in-plane resolution =0.625×0.625 mm², bandwidth =134 Hz/pixel, and T_shift values of 0.6 ms (with echo shift) and 0 ms (without echo shift) (8). The sequence parameters of the second tissue phantom were identical to the first tissue phantom, except that the spectrally selective attenuated inversion recovery (SPAIR) was used for the fat suppression.

In the second phantom experiment, conducted using a tracheal stent [manufactured by Micro-Tech (Nanjing) Co., Ltd., length × diameter =60 mm × 20 mm]. The stent was perpendicularly placed into a gelatin phantom doped with 1.0 g/L copper sulfate solution, and it was made of 0.24 mm diameter Ni-Ti alloy wire (18). In this experiment, two datasets were acquired with or without a T_shift of 0.6 ms using the modified 2D FSE sequence of PC-MRI (8) with the same scan parameters of reference (18): TR =2,000 ms, TE =18 ms, matrix size =192×192, FOV =128×128 mm², resolution is 0.67×0.67 mm², bandwidth =134 Hz/Pixel, slice thickness is 1.5 mm with no gap, and slice/partition number was 20 (18).

During the in vivo radiotherapy experiments, imaging was performed on two prostate tumor patients with permanently implanted brachytherapy seeds. The scans were conducted using a 24-channel flexible receiver coil (SIEMENS AG, Erlangen, Germany), and fat suppression was achieved using a spectral pre-saturation with inversion recovery (SPIR) preparation pulse. The scan parameters were as follows: TE/TR =17/3,200 ms, ETL =7, T_shift =0.6 ms, bandwidth =134 Hz/pixel, matrix size =320×320, FOV =295×295 mm², spatial resolution =0.92×0.92 mm², slice thickness =2.5 mm with no gap, and partition number was 44. All clinical data were obtained from the Shenzhen Institute of Advanced Technology, Chinese Academy of Sciences. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by the Ethics Committee of Shenzhen Institute of Advanced Technology, Chinese Academy of Sciences, and written informed consent was obtained from all patients. Throughout the study, data were acquired using the positive contrast imaging sequence, both with and without T_shift (8).

Parameter settings

In the conventional CG, the numerical conditioning parameter was set to ∈=10⁻⁸, the maximum number of inner iterations was 100, the maximum number of outer iterations was 20, and the tolerance for the inner CG loop was set to 0.01 (the detailed justification for the selection of these parameters was provided in Appendix 2). The PD method implemented in our study employed empirical parameter selection for the all experiments, we set σ=τ=1/L. While the theoretical convergence criteria of this approach assumes L²στ<1, our practical implementation demonstrates consistent convergence behavior with the parameter choice σ=τ=1/L, as pointed out in (24). The relative error norm ‖uk+1−uk‖/‖uk‖ with u^k was the solution of the k-th iteration served as the stopping criterion for both the conventional CG method and the proposed PD algorithm. All of the reconstruction algorithms were executed on a workstation equipped with an Intel Xeon E5-2660 v3 CPU, 128 GB of memory, and an NVIDIA Tesla V100 SXM2 GPU. The software environment consisted of MATLAB R2014a and CUDA Toolkit 10.0 running on Ubuntu 18.04 LTS operating system. Additionally, to enable standardized visual comparison, all reconstructed images were displayed individually using identical grayscale mapping parameters during qualitative assessment.

Quantitative evaluation methods

To objectively evaluate the performance of the reconstruction algorithms, quantitative metrics including the SNR and full width at half maximum (FWHM) were calculated from the reconstructed positive contrast images.

SNR Measurement: SNR was used to assess the image quality and noise level in the reconstructed susceptibility maps. It was calculated using the following formula:

SNR=SIseed/NoiseSD[10]

Where, SI_seed represents the mean signal intensity within a region of interest (ROI) placed over a metallic seed, and Noise_SD is the standard deviation of the signal intensity within a large ROI drawn in a homogeneous background area devoid of any visible artifacts or structures. The average SNR (SNR_Avg) was reported for multiple seeds (e.g., the five seeds in the phantom experiments, three seeds in the in vivo experiments) to ensure statistical reliability.

FWHM Measurement: FWHM was used to quantify the spatial sharpness and accuracy of the seed’s localization. The measurement procedure was as follows:

• Region selection: a linear profile was drawn perpendicularly across the center of a metallic seed in the positive contrast image;
• Intensity profile extraction: the signal intensity values along this line profile were plotted;
• Width calculation: the FWHM was computed as the distance between the two points on the intensity profile where the signal intensity drops to 50% of the peak value. A smaller FWHM indicates a sharper and more precisely localized representation of the seed. The average FWHM (FWHM_Avg) was calculated from multiple seeds to provide a representative measure of sharpness for each method and dataset.

These quantitative analyses were performed on both phantom and in vivo experimental results to provide robust evidence supporting the visual observations.

Convergence behavior

The conventional CG and proposed PD algorithm were run on the simulated susceptibility map using a single seed. The root mean square error (RMSE) was computed to illustrate the convergence behavior of each algorithm. The RMSE was defined as

RMSE=‖uk−u*‖22N[11]

Where N denotes the number of pixels in the reconstructed susceptibility map,u^k denotes the reconstructed susceptibility map obtained at the k-th iteration, and u* denotes the reference image.

Results

Results were organized based on computer simulations, phantom experiment, in vivo experiments with convergence analysis, and quantitative analysis of the running time.

Computer simulations

Figure 2 shows the positive contrast images obtained by the conventional CG and PD algorithms. A comparison with the PD algorithm was conducted using the conventional CG algorithm under different numerical conditioning settings: ϵ=10⁻⁶ and ϵ=10⁻⁸. The results indicated that both the conventional CG (ϵ=10⁻⁸) and PD algorithms could get comparable results with the true positive contrast image. In principle, the PD method was not influenced by variations in the parameters ϵ.

Figure 2 The positive contrast images were acquired by the conventional CG iterative algorithm with different numerical conditioning ϵ=10⁻⁶ and ϵ=10⁻⁸, and by the PD algorithm. CG, conjugate gradient; PD, primal-dual.

In addition, Figure 3 shows the quantitative comparison of the reconstruction accuracy between CG and PD optimization methods in simulated data. Figure 3A displays the simulated susceptibility distribution employed for method validation, which features a uniform square object. The red dotted line indicates the profile along which quantitative comparisons were conducted between the true susceptibility and the reconstructed values.

Figure 3 Quantitative comparison of the reconstruction accuracy between CG and PD optimization methods using simulated data. (A) Schematic illustration of the simulated susceptibility distribution: a uniform square object (white bounding box). The red dotted line indicates the profile line where quantitative measurements were obtained. (B) Profile plots along the red dotted line comparing the true susceptibility distribution (red solid line) with reconstructions from CG under two numerical conditioning (ϵ=10⁻⁸, light blue line with inverted triangles; ϵ=10⁻⁶, yellow line) and PD (purple solid line). The x-axis represents spatial coordinates along the profile (in pixels units), and the y-axis represents susceptibility values (×10⁻⁵ ppm). CG, conjugate gradient; PD, primal dual; ppm, parts per million.

Figure 3B shows a one-dimensional cross-section intersecting the seed, further demonstrating that the PD iteration algorithm achieves higher positive contrast than conventional CG iteration. Although the conventional CG with conditioning ∈=10⁻⁸ could obtain the comparable results, the susceptibility values obtained with the PD method were closer to the true value than those from the conventional CG method.

Phantom experiment

Figure 4 shows experimental results from two pork tissue phantoms, with seeds oriented perpendicular to the B₀ field. The scan parameters of tissue 2 were identical to those of tissue 1, except that the spectrally selective attenuated inversion recovery was used for the fat suppression. Figure 4A shows a photograph of the tissue phantom, while Figure 4B,4C displays the magnitude images of the center slice.

Figure 4 Imaging results from the two-tissue phantom experiment. (A) Photograph of the phantom, with a red arrow indicating metallic seeds, a green arrow indicating a bone fragment, and a blue arrow indicating a plastic object. The seeds were placed perpendicular to the B0 field. (B,C) Magnitude images of Tissue 1 and Tissue 2, respectively. (D,E) Positive contrast images of Tissue 1 and Tissue 2 were reconstructed using the conventional CG algorithm. (F,G) Positive contrast images of the same tissues reconstructed using the PD algorithm. Zoomed-in views of the areas within the white boxes are provided, allowing for a detailed comparison of the improved sharpness achieved by the PD algorithm. CG, conjugate gradient; PD, primal-dual.

Notably, small regions surrounding the five seeds exhibited some signal enhancement due to distortion at TE. Nevertheless, dark spots appeared around each seed, covering up to 25 pixels (approximately 10 mm²), resulting from susceptibility artifacts. Consequently, accurately pinpointing the exact location of the seeds proved challenging. Figure 4D,4E shows the positive contrast image of the two tissue phantoms generated by the conventional CG iteration algorithm, while Figure 4F,4G presents those produced by the PD iteration algorithm. The results indicate that PD achieves comparable or superior results to conventional CG while requiring significantly shorter reconstruction time.

Figure 5 demonstrates the positive contrast visualization of the stent-graft using both methods across different image slices. Figure 5A provides a photograph of the stent-graft phantom. Figure 5B displays the magnitude image of the normal MR images.

Figure 5 Images from the stent phantom: (A) photograph of the phantom; (B) magnitude image; (C,D) positive contrast images generated by the conventional CG iteration algorithm; (E,F) positive contrast images generated by the PD iteration algorithm. Note that the stent was placed in the phantom such that they were parallel to the B0 field. CG, conjugate gradient; PD, primal-dual.

Figure 5C,5D shows the PC MR images of different stent slices reconstructed by the conventional CG method. Figure 5E,5F presents the results obtained by the PD iteration algorithm. The PD method achieved higher contrast visualization of the stent compared to the conventional CG method.

In vivo experiment

The results of two in vivo radiotherapy experiments are shown in Figure 6. Compared to the homogeneous gelatin and porcine muscle tissue phantoms, human tissue exhibited greater heterogeneity, resulting in significantly more inhomogeneous and noisier images in PC-MR. This increased noise posed greater challenges for visualizing and precisely localizing the seeds. However, most seeds (indicated by the white arrow) remained detectable in the positive contrast images generated by both CG and PD methods. While the PD method yields brachytherapy metal seeds that are both sharper and greater in number, demonstrating superior performance than the CG method.

Figure 6 The positive contrast images of brachytherapy metal seeds for two radiotherapy patients with prostate tumors. (A,D) Magnitude images, where the red box indicated the prostate region containing a tumor, and the white arrow pointed to a radioactive particle (e.g., a brachytherapy seed); (B,E) axial PC images were generated by the CG method; and (C,F) axial PC images were generated by the PD method. Zoomed-in views of the areas within the white boxes were provided, allowing for a detailed comparison of the enhanced sharpness achieved by the PD algorithm. CG, conjugate gradient; PC, positive contrast; PD, primal-dual.

Quantitative analysis

To substantiate the visual observations of improved image quality, a quantitative analysis was conducted by calculating the FWHM and SNR near the metallic seeds for both the CG and PD methods, as detailed in the Methods section. The results are summarized in Figure 7. In the phantom experiments, the PD method consistently yielded lower FWHM values compared to the CG method (e.g., Tissue 2: 1.21 for PD vs. 1.26 for CG), indicating sharper delineation of the seeds. Correspondingly, the PD method achieved higher SNR values (e.g., Tissue 2: 190.21 for PD vs. 189.56 for CG), demonstrating improved noise suppression. The superiority of the PD method was more pronounced in the in vivo experiments. For Patient 2, the PD method reduced the FWHM from 1.15 (CG) to 1.02 (PD) and increased the SNR from 110.34 (CG) to 131.65 (PD), representing a significant enhancement in both image sharpness and quality. These quantitative metrics provided strong evidence supporting the visual findings that the PD method achieved superior positive contrast imaging.

Figure 7 Quantitative comparison of image sharpness and quality. (A) A representative positive contrast image showing a metallic seed and the line profile (blue line) used for FWHM measurement. The inset shows a zoomed-in view of the seed. (B) The susceptibility variance along the blue line in (A), illustrating the FWHM calculation. The x-axis represents sliding window width (pixels), and the y-axis represents variance of susceptibility values (×10⁻⁵). (C) The average FWHM (FWHM_Avg) and average SNR (SNR_Avg) for the CG and PD methods across phantom and in vivo experiments. The SNR_Avg/FWHM_Avg respectively represent the average SNR/FWHM of the five seeds (the white square box in Figure 4) or the average SNR/FWHM of the three seeds (the white square box in Figure 6). CG, conjugate gradient; FWHM, full width at half maximum; PD, primal-dual; SNR, signal-to-noise ratio.

Convergence analysis

The convergence behavior of conventional CG and PD methods on simulated data one is shown in Figure 8. The Y-axis represents the RMSE between the reconstructed images and the reference, while the X-axis represents either the iteration number or the computation time. As seen in Figure 8A, the PD method (red solid line) converged much faster than conventional CG with a numerical conditioning of 10⁻⁶ (blue line), and as rapidly as conventional CG with a numerical conditioning of 10⁻⁸ (green line) in the first couple of iterations (iterations <14). However, CG with 10⁻⁸ became less accurate in later iterations (iterations >14), although it nearly converged to the same RMSE as PD in the final iteration (iterations >1,400). After 2,000 iterations, the relative errors were 1.005% for the proposed PD method and 1.006% for conventional CG. As shown in Figure 8, the PD algorithm reached a relative error of 1.005% at 400 iterations faster than the conventional CG algorithm. The PD algorithm demonstrates superior convergence speed in terms of both iteration count and actual computation time, highlighting its strong potential for time-sensitive clinical inline reconstruction.

Figure 8 Convergence behavior of the conventional CG and the PD algorithm using the simulated data (shown in Figure 2). (A) The convergence plot versus the number of iterations. Although the conventional CG with ϵ=10⁻⁸ performed as rapidly as the proposed PD in the initial few iterations, the conventional CG did not achieve the accuracy that the PD does (as indicated by the gap between the red and green lines and the zoom-in plot of the first 300 iterations). (B) Convergence plot versus computation time, displayed for the first 450 seconds (CPU time in seconds). The PD algorithm exhibits superior convergence speed in terms of both iteration count and actual computation time, highlighting its strong potential for time-sensitive clinical inline reconstruction. CG, conjugate gradient; CPU, central processing unit; PD, primal-dual; RMSE, root mean square error.

Running time analysis

Quantitative analysis of the running time for the conventional CG method, the PD iteration technique, and the GPU-based PD method using simulated data, phantom experiments and in vivo experiments is shown in Table 1. The total number of iterations for both the conventional CG and PD was 2,000 in the simulated data. In addition, the number of inner iterations per outer loop for the conventional CG method was 100 and the 20 outer iterations were used with ϵ=10⁻⁸. The running times of each algorithm were approximately 2,693.9 seconds for the conventional CG and 962.4 seconds for the PD method with the simulated data, respectively. Therefore, the reconstruction computational time of the PD method was 2–3 times faster than that of the conventional CG iteration method when the same total number of iterations was used.

All of the phantom experiments, the total number of iterations for the conventional CG was set to 2,000, with 20 outer iterations. The relative error between two adjacent iterations (10⁻² for the conventional CG and 10⁻³ for the PD algorithm) was used to balance the image quality and reconstruction speed. Additionally, in all phantom experiments, the conventional CG iterations stopped at 2,000 iterations without ever reaching the target relative error of 10⁻². In contrast, the PD algorithm successfully reached the given relative error of 10⁻³ and stopped at 455 iterations (tissue 1 experiment), 355 iterations (tissue 2 experiment), and 225 iterations (stent phantom experiment), respectively. Comparing to the conventional CG iteration method, the PD method demonstrated 2–4 times faster convergence while still maintaining accurate visualization and localization of interventional devices with various shapes and sizes.

For the patient experiment, the iteration counts for PD were 190, and the outer iteration counts for CG were 6. The relative error between two adjacent iterations was simultaneously set to 10⁻² for both methods to balance the image quality and reconstructing speed. The running times for each algorithm were approximately 265.64 seconds for the conventional CG method, and 128.14 seconds for the PD on the patient 1 data and 289.53 seconds for the conventional CG method, and 100.11 seconds for the PD on the patient 2 data, respectively. After GPU acceleration, the PD reconstruction time was reduced to 64.25 seconds on simulated data, 5.22 seconds on the tissue 1 phantom data, 5.91 seconds on the tissue 2 phantom data, 5.65 seconds on the stent phantom data, 32.04 seconds on the patient 1 data, and 20.03 seconds on the patient 2 data. The experimental results demonstrated that GPU-based PD achieved significant acceleration in PC-MR reconstruction compared to PD reconstruction without GPU acceleration. The proposed architecture achieved 4–15 times speedup while maintaining virtually identical image quality.

Discussion

This work demonstrates that the accuracy and convergence of the conventional CG method highly depend on the conditioning parameter ϵ. The accuracy of conventional CG can be improved by reducing the tolerance parameter ϵ from 10⁻⁶ to 10⁻⁸, however, this adjustment also results in a substantial slowdown in the algorithm’s convergence rate. This occurs because the linear system solved by the nested CG loop becomes ill-conditioned. More precisely, the convergence rate of the conventional CG iteration is governed by the distribution (spectrum) of the eigenvalues of the linear system (27), and the value of the conditioning parameter ϵ influences the location of the spectrum. Therefore, from a computational perspective, the algorithm becomes impractical when ϵ is set to a very small value (e.g., 10⁻¹⁰) to achieve higher accuracy. Moreover, this limitation makes the conventional CG less desirable compared to the proposed PD approach, as the numerical conditioning of the Euler-Lagrange equation deteriorates, potentially leading to an inaccurate suboptimal solution. However, when addressing non-smooth convex optimization problems, the continuation method provides a solution (28). Specifically, at each outer iteration, the tolerance parameter ϵ is set according to the rule ϵ=γ.ϵ, where γ represents a decreasing factor.

In the first couple of iterations, the conventional CG algorithm appears faster than the PD method, as shown in Figure 8. However, when high accuracy is required, the conventional CG method becomes less preferable compared to the PD method. This is because the nested CG loop requires more iterations to achieve a tighter tolerance, in addition to needing more outer iterations. In contrast, the PD is the fastest and demonstrates better accuracy even with early termination. Unlike the CG method, the PD algorithm directly addresses the exact non-smooth optimization problem and attains a sufficiently optimal solution within 500–2,000 iterations, enabling both rapid and precise reconstruction. Since the PD algorithm solely involves matrix-vector multiplications and coordinate-wise projections, it can reduce reconstruction time by up to 2 to 4 times compared to the conventional CG iterations, and even by as much as 15 times when using GPUs. Additionally, the PD method offers more flexible parameter selection (29,30), which further enhances its convergence speed.

While λ is a critical parameter that balances data consistency and spatial smoothness, and thus requires careful consideration in the work. Initial experience indicates that a proper parameter λ can better suppress noise and streaking artifacts commonly found in susceptibility estimation. Other parameters in both the conventional CG and the PD algorithms are fixed, as they can provide comparable quality. Therefore, there are no free parameters except λ. Similar to other minimization problems, the performance of the algorithm depends on the choice of the regularization parameter λ, which can be set empirically. Figure 9 presents the reconstruction results on the computer simulation experiment using different λ. The RMSE of the 5×5 pixels region of the seeds is used to evaluate the reconstruction. The results in Figure 9 show that the RMSE values of the PD method are lower than those of the conventional CG method across all tested λ values. Therefore, better positive contrast visualization can be achieved using the PD iteration algorithm. Moreover, the PD method allows a wider range of suitable parameters (λ=2−16). Based on these findings, we set λ=10 for the computer simulation.

Figure 9 The RMSE of the two iteration methods with different reconstruction parameter λ. CG, conjugate gradient; PD, primal-dual; RMSE, root mean square error.

In addition, we incorporate a first-order central difference scheme to promote sparsity in our approach. However, the central difference scheme’s insensitivity to checkerboard-like noise may potentially introduce artifacts into the final susceptibility map. We indeed observed such artifacts in our phantom experiments, though the checker-board artifacts did not s adversely affect the final positive contrast image quality, as shown in Figure 4. Nevertheless, the quality of the PC MRI can be further improved by employing alternative gradient schemes, such as the forward difference method or the second-order total generalized variation [TGV (15,31)], as well as by incorporating other sparse transformation techniques, such as the wavelet transform (32,33).

In the in vivo radiotherapy experiments, the 128-second reconstruction time for PD algorithm represented a 52% reduction compared to conventional CG’s 265.64 seconds on the same dataset. It represented a significant advancement towards feasible inline reconstruction in clinical settings. With GPU implementation, reconstruction times were further reduced to 32.04 seconds for patient 1 data. So, the GPU-based accelerated implementation (20–32 seconds) established genuine feasibility for inline reconstruction in clinical practice.

While AI-based methods would have the potential for faster reconstruction times, which would be meaningful to compare with the model-based PD method. However, AI-based methods require extensive training datasets and may suffer from performance degradation, especially in scenarios with limited training data or when the imaging conditions significantly differ from the training set. Our future research could explore hybrid approaches that combine the strengths of both model-based optimization and AI-based techniques to achieve even faster and more accurate MRI reconstruction.

Conclusions

A novel susceptibility-driven positive contrast imaging technique is proposed based on the PD iteration technique. The results from simulations, phantom experiments, and in vivo experiments demonstrate that the PD method has shown advantages in convergence and accuracy over the conventional CG technique. It not only reduces the reconstruction time by up to 2–4 times but also achieves better image quality and provides easily adjustable parameters for PC-MRI of the interventional devices. The reconstruction speed can be increased to 4–15 times after GPU acceleration.

Acknowledgments

None.

Footnote

Data Sharing Statement: Available at https://qims.amegroups.com/article/view/10.21037/qims-2025-2018/dss

Funding: This study was partially supported by grants from the National Key R&D Program of China (No. 2021YFF0501402), Guangzhou Medical University Research Capability Enhancement Project (No. 02-410-2405125), the Science and Technology Plan Program of Guangzhou (No. 2025A04J4360), and the Natural Science Foundation of Guangdong Province (No. 2026A050002626), as well as by the Guang Dong Engineering Technology Research Center of Respiratory Rehabilitation Medicine.

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://qims.amegroups.com/article/view/10.21037/qims-2025-2018/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by the Ethics Committee of Shenzhen Institute of Advanced Technology, Chinese Academy of Sciences, and written informed consent was obtained from all patients.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. D'Amico A, Cormack R, Kumar S, Tempany CM. Real-time magnetic resonance imaging-guided brachytherapy in the treatment of selected patients with clinically localized prostate cancer. J Endourol 2000;14:367-70. [Crossref] [PubMed]
2. Lammer J, Bosiers M, Zeller T, Schillinger M, Boone E, Zaugg MJ, Verta P, Peng L, Gao X, Schwartz LB. First clinical trial of nitinol self-expanding everolimus-eluting stent implantation for peripheral arterial occlusive disease. J Vasc Surg 2011;54:394-401. [Crossref] [PubMed]
3. Raval AN, Telep JD, Guttman MA, Ozturk C, Jones M, Thompson RB, Wright VJ, Schenke WH, DeSilva R, Aviles RJ, Raman VK, Slack MC, Lederman RJ. Real-time magnetic resonance imaging-guided stenting of aortic coarctation with commercially available catheter devices in Swine. Circulation 2005;112:699-706. [Crossref] [PubMed]
4. Dolidze DD, Kovantsev SD, Bagatelia ZA, Bumbu AV, Barinov YV, Chechenin GM, Pichugina NV, Gogolashvili DG. Ultrasound-guided core-needle biopsy for diagnosis of thyroid cancer. Khirurgiia (Mosk) 2025;87-95. [Crossref] [PubMed]
5. Wachowicz K, Thomas SD, Fallone BG. Characterization of the susceptibility artifact around a prostate brachytherapy seed in MRI. Med Phys 2006;33:4459-67. [Crossref] [PubMed]
6. Frahm C, Gehl HB, Melchert UH, Weiss HD. Visualization of magnetic resonance-compatible needles at 1.5 and 0.2 Tesla. Cardiovasc Intervent Radiol 1996;19:335-40. [Crossref] [PubMed]
7. Glowinski A, Adam G, Bücker A, van Vaals J, Günther RW. A perspective on needle artifacts in MRI: an electromagnetic model for experimentally separating susceptibility effects. IEEE Trans Med Imaging 2000;19:1248-52. [Crossref] [PubMed]
8. Shi C, Xie G, Zhang Y, Zhang X, Chen M, Su S, Dong Y, Liu X, Ji J. Accelerated susceptibility-based positive contrast imaging of MR compatible metallic devices based on modified fast spin echo sequences. Phys Med Biol 2017;62:2505-20. [Crossref] [PubMed]
9. Nosrati R, Paudel M, Ravi A, Pejovic-Milic A, Morton G, Stanisz GJ. Potential applications of the quantitative susceptibility mapping (QSM) in MR-guided radiation therapy. Phys Med Biol 2019;64:145013. [Crossref] [PubMed]
10. Dong Y, Chang Z, Ji J. Imaging and localizing interventional devices by susceptibility mapping using MRI. Annu Int Conf IEEE Eng Med Biol Soc 2014;2014:1541-4. [Crossref] [PubMed]
11. Daubechies I, Defrise M, De Mol C. An iterative thresholding algorithm for linear inverse problems with a sparsity constraint. Commun Pure Appl Math 2004;57:1413-57.
12. Zhao Y, Allen-Zhao Z, Wang L, He X, Mao Q. Inertial primal-dual projection neurodynamic approaches for constrained convex optimization problems and application to sparse recovery. Neural Netw 2025;186:107274. [Crossref] [PubMed]
13. Chambolle A, Pock T. A first-order primal-dual algorithm for convex problems with applications to imaging. J Math Imaging Vis 2011;40:120-45.
14. Cheng J, Jia S, Wang H, Liang D. MR image reconstruction using the Chambolle-Pock algorithm. In Proceedings of the Joint Annual Meeting of ISMRM 2018;3552.
15. Knoll F, Bredies K, Pock T, Stollberger R. Second order total generalized variation (TGV) for MRI. Magn Reson Med 2011;65:480-91. [Crossref] [PubMed]
16. Langkammer C, Bredies K, Poser BA, Barth M, Reishofer G, Fan AP, Bilgic B, Fazekas F, Mainero C, Ropele S. Fast quantitative susceptibility mapping using 3D EPI and total generalized variation. Neuroimage 2015;111:622-30. [Crossref] [PubMed]
17. Chatnuntawech I, McDaniel P, Cauley SF, Gagoski BA, Langkammer C, Martin A, Grant PE, Wald LL, Setsompop K, Adalsteinsson E, Bilgic B. Single-step quantitative susceptibility mapping with variational penalties. NMR Biomed 2017;
18. Shi C, Xie G, Liang D, Wang H, Huang Y, Ren Y, Xue Y, Chen H, Su S, Liu X. Positive visualization of MR-compatible nitinol stent using a susceptibility-based imaging technique. Quant Imaging Med Surg 2019;9:477-90. [Crossref] [PubMed]
19. Vafay Eslahi S, Ji J. Accelerated positive contrast MRI of interventional devices using parallel compressed sensing imaging. Magn Reson Imaging 2019;60:130-6. [Crossref] [PubMed]
20. Eslahi SV, Dhulipala PV, Shi Caiyun, Xie Guoxi, Ji JX. Parallel compressive sensing in a hybrid space: Application in interventional MRI. Annu Int Conf IEEE Eng Med Biol Soc 2017;2017:3260-3. [Crossref] [PubMed]
21. Valkonen T. A primal-dual hybrid gradient method for nonlinear operators with applications to MRI. Inverse Probl 2014;30:055012.
22. Kirimtat A, Krejcar O. GPU-Based Parallel Processing Techniques for Enhanced Brain Magnetic Resonance Imaging Analysis: A Review of Recent Advances. Sensors (Basel) 2024.
23. Gai ND. Highly Efficient and Accurate Deep Learning-Based Classification of MRI Contrast on a CPU and GPU. J Digit Imaging 2022;35:482-95. [Crossref] [PubMed]
24. Sidky EY, Jorgensen JH, Pan X. Convex optimization problem prototyping for image reconstruction in computed tomography with the ChambollePock algorithm. Phys Med Biol 2012;57:3065-91. [Crossref] [PubMed]
25. Wang H, Peng H, Chang Y, Liang D. A survey of GPU-based acceleration techniques in MRI reconstructions. Quant Imaging Med Surg 2018;8:196-208. [Crossref] [PubMed]
26. Quelhas KN, Henn MA, Farias R, Tew WL, Woods SI. GPU-accelerated parallel image reconstruction strategies for magnetic particle imaging. Phys Med Biol 2024;
27. Yahaya J, Kumam P, Salisu S, Sitthithakerngkiet K. Spectral-like conjugate gradient methods with sufficient descent property for vector optimization. PLoS One 2024;19:e0302441. [Crossref] [PubMed]
28. Pan S, Dai Q, Chen H. Global optimality analysis and solution of the ℓ(0) total variation signal denoising model. Math Biosci Eng 2023;20:6932-46. [Crossref] [PubMed]
29. Loris I, Verhoeven C. On a generalization of the iterative soft-thresholding algorithm for the case of non-separable penalty. Inverse Probl 2011;27:125007.
30. Chen P, Huang J, Zhang X. A primal-dual fixed point algorithm for minimization of the sum of three convex separable functions. Fixed Point Theory Appl 2016;54.
31. Liu Z, Li Y, Wang W, Liu L, Chen R. Mesh Total Generalized Variation for Denoising. IEEE Trans Vis Comput Graph 2022;28:4418-33. [Crossref] [PubMed]
32. Kumar A, Tomar H, Mehla VK, Komaragiri R, Kumar M. Stationary wavelet transform based ECG signal denoising method. ISA Trans 2021;114:251-62. [Crossref] [PubMed]
33. Ding Y, Xu X, Peng L, Zhang L, Li W, Cao W, Gao X. Wavelet transform-based frequency self-adaptive model for functional brain network. Cereb Cortex 2023;33:11181-94. [Crossref] [PubMed]