Accurate Concentration Recovery for Quantitative Magnetic Particle Imaging Reconstruction via Nonconvex Regularization

Magnetic particle imaging (MPI) uses nonlinear response signals to noninvasively detect magnetic nanoparticles in space, and its quantitative properties hold promise for future precise quantitative treatments. In reconstruction, the system matrix based method necessitates suitable regularization terms, such as Tikhonov or non-negative fused lasso (NFL) regularization, to stabilize the solution. While NFL regularization offers clearer edge information than Tikhonov regularization, it carries a biased estimate of the $\mathbf {l}_{\mathbf {{1}}}$ penalty, leading to an underestimation of the reconstructed concentration and adversely affecting the quantitative properties. In this paper, a new nonconvex regularization method including min-max concave (MC) and total variation (TV) regularization is proposed. This method utilized MC penalty to provide nearly unbiased sparse constraints and adds the TV penalty to provide a uniform intensity distribution of images. By combining the alternating direction multiplication method (ADMM) and the two-step parameter selection method, a more accurate quantitative MPI reconstruction was realized. The performance of the proposed method was verified on the simulation data, the Open-MPI dataset, and measured data from a homemade MPI scanner. The results indicate that the proposed method achieves better image quality while maintaining the quantitative properties, thus overcoming the drawback of intensity underestimation by the NFL method while providing edge information. In particular, for the measured data, the proposed method reduced the relative error in the intensity of the reconstruction results from 28% to 8%.


I. INTRODUCTION
M AGNETIC particle imaging (MPI) is a novel imaging modality first introduced by Gleich and Weizenecker in 2005 [1].The spatial distribution of superparamagnetic iron oxide nanoparticles (SPIONs) can be reconstructed by utilizing the nonlinear magnetization response of SPIONs in alternating magnetic fields.As a noninvasive imaging method, MPI has highly sensitive and quantitative characteristics, which makes it promising for use in biological applications [2], [3], [4], [5], [6], [7].Examples include drug delivery [2], visualization of vascular interventions [3], in vivo tracking of minute numbers of cells [4], and quantitative monitoring of tumor immunotherapy [5].
In MPI reconstruction, the spatial distribution of SPIONs can be decoded using a system-matrix-based approach.However, because the system matrix-based inverse problem is ill posed, a direct solution is unstable [8], [9].Therefore, regularization, such as the Tikhonov penalty [10], nonnegative Garrote penalty [11], lasso penalty [12], nonnegative fused lasso (NFL) regularization [13], and elastic net regularization [14], is required to stabilize image reconstruction.Among these, the Tikhonov penalty is the most widely used regularization approach in MPI reconstruction [8], [10].Combined with the Kaczmarz method, this enables a fast solution using row iterations.However, Tikhonov regularization has a limited ability to suppress noise and over-smooth images.Considering the sparse and piecewise constant distribution of magnetic particles, Storath et al. introduced NFL regularization by combining the l 1 and total variation (TV) norms [13].This method can effectively suppress noise and maintain edge information in an image.Although the l 1 norm in NFL regularization is the most effective sparsity-induced term among convex regulars, it is a biased estimate that tends to underestimate the true solution [15], [16].This reduces the concentration of magnetic particles in the reconstructed images and weakens the advantage of MPI in terms of quantitative imaging.
To reduce the bias effect while maintaining sparsity, nonconvex regularizations, such as the smoothly clipped absolute deviation [17], log sum penalty [18], l q -norm penalty (0 < q < 1) [19], and min-max concave (MC) penalty [20], have been introduced.The optimization properties of nonconvex regularization are more similar to those of the l 1 -norm penalty than to those of the ideal l 0 -norm penalty.Therefore, nonconvex regularization can improve the reconstruction accuracy and suppress the bias effect.Several techniques can be adopted to solve nonconvex optimization problems, such as the Bayesian compressive sensing [21], focal underdetermined system solutions [22], and iterative reweighted least squares [23] methods.However, these methods are not suitable for solving large-scale problems because of their high computational costs.Selesnick proved that a nonconvex optimization problem has a simple closed-form solution for proximal operators when the penalty satisfies specific conditions (such as the MC penalty) [24], which can accelerate the solution of nonconvex optimization problems.Furthermore, when more than one penalty is included in the regular term, the alternating direction method of multipliers (ADMM) is suitable.It simplifies the solution by decomposing a complex problem into several simple subproblems.Recently, ADMM reconstruction has been applied to MPI reconstruction and has exhibited good performance [25].
In this paper, we propose a new regularization term for MPI reconstruction that enables accurate MPI concentration reconstruction, achieving better image quality while maintaining the quantitative advantages of MPI.The contributions of this work are as follows: (i) a regularization that better matches the quantitative reconstruction of MPI is proposed; (ii) a suitable solution method and a two-step parameter selection method are designed, which simplify the selection process of multiple parameters and realize the fast solution of the objective function; and (iii) the effectiveness and advantages of the proposed method are verified through experiments using Open-MPI data and measured data.
To evaluate the performance of the proposed method, we used different phantoms, including stenosis, overlapping ellipses with different intensities, and vascular trees.Nonnegative Tikhonov regularization utilizing the Kaczmarz method and NFL regularization employing the forward-backward (FB) method were applied for comparison.The results demonstrated that the proposed nonconvex regularization solved with ADMM could provide precise intensity estimation while maintaining edge information.

A. Forward Problem
The forward problem of MPI involves establishing of a mathematical model of the process in which the nonlinear magnetization of SPIONs is excited by an alternating magnetic field and signal encoding within the signal chain.Assuming that the particles are always in thermal equilibrium, the nonlinear magnetization properties of SPIONs can be modeled by the Langevin function L(.), expressed as: where c(r ) is the concentration of SPIONs at position r , Ms is the magnetization of the particles at unit concentration, m is the magnetic moment, and e H = H (r, t) /||H (r, t)|| 2 is the normalized direction of the magnetic field strength H (r, t).µ 0 denotes the permeability of vacuum, k B denotes the Boltzmann constant and T is the temperature.Signal encoding is based on the MPI signal chain.The relationship between the particle distribution c and the induced voltage can be expressed as: where p (r ) is the sensitivity of receive coils.Using the Fourier transform and discretization, the model can be expressed as: where S ∈ C M×N denotes the system matrix, M is the number of frequencies, and N is the number of discrete points.c = [c 1 , c 2 , . . ., c N ] T ∈ R N ×1 denotes the particle concentration at N positions and u = [u 1 , u 2 , . . ., u N ] T ∈ C M×1 denotes the induced voltage signal in the frequency domain.

B. Reconstruction Methods
The reconstruction problem of the MPI is aimed at solving the linear problem (3).However, this approach can only acquire an approximate optimal solution because the system matrix of the MPI is ill-conditioned, which means that the singular values decay rapidly, and a small noise disturbance has a significant impact on the reconstruction results [26].Therefore, a regularization term is typically added when solving an inverse problem.The Tikhonov and NFL regularization are the two most well-known types of regularization.The former is the most widely used penalty in MPI reconstruction, whereas the latter is a state-of-the-art penalty.
1) Nonnegative Tikhonov Regularization: The objective function of the nonnegative Tikhonov regularization is: where • 2 2 denotes the Euclidean norm.Equation ( 4) is a typical least-squares optimization problem that can be solved using singular value decomposition [27] and conjugate gradient [28] methods.However, the computational cost is high when the system matrix is large.
Weizenecker introduced the Kaczmarz method to obtain a fast solution, which is calculated using row iterations and requires only a few iterations to reach convergence [29].
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
In addition, a weighting matrix was employed to weigh the rows of the system matrix by utilizing the row energy.Preprocessing can improve image quality and shorten the reconstruction time [10].Because of its short convergence time and good robustness, nonnegative Tikhonov regularization based on the Kaczmarz method is popular in practice [8], [30].
2) NFL Regularization: To improve the image quality, Storath introduced a near-isotropic TV norm to consider prior neighborhood structures while using the l 1 -norm to induce sparsity [13].The objective function can be expressed as: where • 1 denotes the l 1 -norm of c and T V near −iso (•) denotes the near-isotropic TV norm regularization.Compared to nonnegative Tikhonov regularization, this method suppresses more noise and preserves better edge information.Furthermore, the FB algorithm can avoid excessive computational costs, even if the objective function is complex.

III. PROPOSED MC REGULARIZATION METHOD
The two aforementioned popular MPI reconstruction methods have their own advantages and disadvantages in terms of reconstructed image quality.Tikhonov-based reconstruction can achieve fast and relatively high-quality reconstruction.However, the l 2 -norm penalty cannot provide sufficient noise suppression in the case of high signal noise, and it smooths the edges of magnetic particles, which is not conducive to obtain detailed information.The NFL method can achieve stronger noise suppression and edge recovery.However, its ability to predict particle concentration is not sufficiently accurate, which affects the quantitative properties of MPI imaging.This situation may be caused by the biased effect of the l 1 -norm.Accordingly, we propose the use of low-bias sparse regularization to improve the quality of MPI-reconstructed images and ensure quantitative accuracy while guaranteeing reconstruction quality.

A. Nonconvex Hybrid Penalty Combining MC and TV Regularization
The objective function based on the proposed method is as follows: where |•| T V denotes the isotropic TV norm, |•| MC denotes the MC penalty norm, and ε ≥ 0 depends on the noise variance.λ T V and λ MC are the regularization parameters.The MC function used in this study can be expressed as: where b is the scaling parameter and v is an intermediate parameter.
The MC regularization-based objective function is given as min λ is the regularization parameter.The objective function can be equally written in element form as 1 2 In this case, the objective function is convex when 1−λb 2 ≥ 0. Let θ = 1/b 2 .The firm threshold function is given as follows, with λ > 0 and θ > 0 [24]: The l 1 -norm problem is expressed as: The solution can be calculated using the soft-threshold function: where sign (•) is a symbolic function.
The derivation of the l 1 -norm and MC penalty is formulated as: A comparison of the l 0 -norm, l 1 -norm and MC penalty demonstrates the rationale for using the MC penalty (see Fig. 1).As shown in Fig. 1(a), a comparison of the penalties and threshold functions reveals that that the MC penalty is a more accurate estimator of the ideal l 0 -norm penalty than other penalties.In addition, a comparison of the derivatives shows that the MC penalty begins penalizing by applying the same rate as the l 1 -norm penalty.However, the l 1 -norm penalty tends to over-penalize which introduces estimation bias when it deviates from zero.In contrast, the MC penalty relaxes the penalization rate to zero and reduces the bias.Therefore, the reconstruction results obtained using the MC penalty can suppress the bias existing in the l 1 -norm penalty.Moreover, considering the continuity and uniformity of the SPION distribution in MPI images, TV regularization is also applied in our method to ensure a uniformly distributed intensity in the reconstructed results (see Fig. 1(b)).Based on the above prior knowledge, the proposed optimization method involving MC and TV regularization is reasonable.

B. Nonconvex Hybrid Penalty Based on the ADMM
Considering the advantages of the ADMM in solving large-scale sparse models, we apply the ADMM to solve the function (6) while utilizing the hybrid cost function form, LU decomposition, firm threshold, and Chambolle projections to speed up convergence.The detailed derivation process is as follows.
To solve the optimization problem of the two separable objective functions, we employ an ADMM with a hybrid cost function.The variable splitting scheme [31], [32] is employed to convert the function (6) into the following form: Here, z = [z 0 , z 1 , z 2 ] H , P = [S H I I] H , Q = −I and s = 0 are the introduced variables.[•] H is the conjugate transpose of the matrix, ι E(ε,I,u) z 0 accounts for the indicator function associated with the data fidelity constraint Sc − u 2 ≤ ε, and I is an identity matrix.
To improve the robustness of the algorithm, the augmented Lagrangian method is adopted.This method does not require strict convexity or finiteness of the objective function and can realize fast iterative convergence of the ADMM.The ADMM steps include the updating of variables c, z, and the Lagrange multiplier T , as shown in ( 15)-( 21).
Parameter β is the penalty coefficient of the Lagrange function [33].
By solving the above ADMM subproblems, ( 14) can be calculated.

C. Fast Solution of Subproblems
In (15), the direct calculation of the inversion term 2I + S H S −1 may be time-and memory-consuming owing Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
to the high matrix dimension.In particular, for 3D-MPI reconstruction, the number of pixels N is usually much greater than the frequency number M. Therefore, we adopted the LU decomposition scheme prior to inversion to avoid the direct inversion operation as: here 16) can be solved using the Moreau proximal map ι E(ε,I,u) (Sc k+1 − d 0 k ), and it can be easily calculated by the orthogonal projection of Gc k+1 − d 0 k on the closed ε-radius ball centered at u.The solution is as follows [33]: Similarly, the Moreau proximal form of ( 17) can be described as Here, we adopt an isotropic TV norm that can be quickly calculated using Chambolle projections as follows [34]: where p = p 1 , p 2 is the dual variable of z, ∇ denotes the gradient operator, div denotes the divergence operator, and τ represents the iteration step size, which is usually set to τ =0.248 [34].Equation ( 18) can be formulated as an MC-based optimization problem: To maintain the convexity of the subproblem, λ MC ≤ βθ is needed, and it can be solved using (10).

D. Parameter Selection
The complete procedure is shown in Algorithm 1.In the ADMM solution, five parameters exist: ε, β, θ, λ T V ,and λ MC .It is unrealistic to search for the optimal parameters directly; therefore, some simplifications are used in the selection of the parameters.Here, the parameter θ is set to 2 to maintain the convexity of the MC constrained subproblem and to enable the firm threshold function (10) to realize the fast solution.The data fidelity term ε is related to the noise level during scanning.The higher the noise, the larger ε should be set to provide greater noise suppression.The penalty parameter β Algorithm 1 ADMM with minimax concave & TV regularization for MPI 1: Output: Reconstructed data c ∈ C N ×1 .

3:
Repeat: 4: Calculate c k+1 by ( 21); 5: Until stop criterion reachedor k > K M AX ; 13: Return c affects the convergence speed, and a proper β can reduce the solution time, which is detrimental to the fast convergence of the algorithm.The parameters λ T V > 0 and λ MC > 0 are the coefficients of the two regularizations used.
The parameter ε can be selected roughly based on the standard deviation (STD) of the noise by ε = γ • ST D, γ > 0. For the parameter β, without the support of a priori information, a suitable value for the parameter is difficult to set before iteration.Therefore, based on the residual balancing strategy [35], we use an adaptive method to automatically update β in iterations and adaptively select the appropriate β value: here, τ incr > 1, τ decr > 1 and κ > 1 are parameters.
) are refer to the primal and dual residuals, respectively at iteration k.In general, the selection of multiple parameters is simplified into two steps in our work.First, the ε is selected according to the standard deviation of the noise, and the penalty parameter β is selected with the adaptive method as shown in (27).Second, λ T V and λ MC are screened in (0, 1), where λ T V +λ MC = 1 is satisfied.The algorithm is terminated when the maximum number of iterations or the stopping condition is reached.We set the stop condition as the relative change in the result between two iterations that was lower than the tolerance, described as c k+1 − c k / c k + 10 −4 < tol.All the aforementioned parameters were chosen to achieve the set tolerance or maximum number of iterations.

A. Simulation Experiment Setup
We implemented all the algorithms in MATLAB R2018b and simulated the data according to u = Sc + η, where c Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.denotes the ground truth of size N = 51 × 51, M = 9900 is the number of frequency points, and η is Gaussian noise.Three different phantoms were built in the simulation: a stenosis phantom, a vessel phantom with plaque in the main vessel, and an overlapping elliptical phantom with multiple intensities [13] (see Fig. 2(a)).The stenosis phantom has a simple structure, whereas the vascular and elliptical models are more complex.Additionally, the overlapping elliptical phantom has a more complex concentration distribution.Considering the reconstruction and time consumption, we preset the stop condition for the methods used.For the nonnegative Tikhonov regularization method, we performed 10 iterations.For the NFL and proposed methods, we stopped the iterations when the tolerance was less than 1×10 −3 or the number of iterations reached 40.The NFL method used the same parameter settings as those in the original study [13].
First, in order to verify the effectiveness of the proposed method under different structural complexity and concentration distribution complexity phantoms, reconstruction experiments were performed using the three designed phantoms and compared with the non-negative Tikhonov regularization and the NFL method.The data were corrupted by Gaussian noise with a 30dB signal-to-noise ratio (SNR) Subsequently, the profiles across the reconstructed results were extracted to compare the intensity distribution along the profiles.
Furthermore, for the quantitative reconstruction target in the paper, we focused on the overlapping elliptical phantom with multiple intensities.To compare the robustness of the algorithm under different noise levels, corrupted signals with SNR of 25 dB, 20 dB, and 15 dB were used.Error maps were used to visualize overall error, and the relative error between the reconstruction results and the ground truth in regions A1, A2, and A3 were calculated for quantitative comparison.

B. Open-MPI Dataset
To verify the feasibility of the proposed algorithm for the measured data, the system matrix and signal of phantoms in the Open-MPI dataset were used to verify our method [36].We selected a 2D MPI system matrix with 19 3 voxel scans using a 2D Lissajous trajectory.Because the measurements of the 19 planes are stitched together in the Open-MPI dataset, preprocessing is required when 2D sliced data are necessary.During data processing, frequency data with an SNR lower than three were discarded to obtain a system matrix of dimensions 824×19 3 .The image was then sliced into 19 2D system matrices of dimensions 824×19 3 , corresponding to the 19 planes scanned using a 2D Lissajous trajectory.Similarly, we sliced the measurement data to a size of 19,000 every 1000 times to obtain 19 signals, corresponding to 19 planes of data.Subsequently, a 2D system matrix of each plane with the corresponding phantom signals was obtained.
Experiments were conducted using shape and concentration phantoms from the dataset.After reconstruction using different methods, the absolute reconstructed intensities were compared using a shape phantom, and the accuracy of the concentration gradient was compared when multiple concentrations of particles were present in the field of view.

C. Measured Data
Furthermore, we validated the proposed algorithm by utilizing data measured using a prototype MPI scanner built in our laboratory [37].The scanner was encoded based on the field-free point, and the Cartesian trajectory was used to construct a 2D measurement sequence (selection field gradient: 1.5 T m −1 µ −1 0 in the x and y directions, drive field amplitude: 27 T m −1 µ −1 0 in the x and y directions, excitation frequency: 25 kHz, scanning frequency: 20 Hz).A delta sample with dimensions of 2 × 2 × 1 mm 3 was used to measure the system matrix S ∈ C M×N (frequency component M = 1.25 × 10 6 , position N = 121).The background signal was measured at each position to reduce the effect of dynamic noise and was utilized for frequency selection.A bar phantom containing two rectangular cavities was employed for the measurements and imaging tests.The edge-to-edge distance between the two cavities was 4 mm.A schematic diagram and the dimensional size are shown in Fig. 2(b).As mentioned above, Tikhonov regularization, NFL regularization and the proposed method were utilized to reconstruct the images, and the results were compared using qualitative and quantitative metrics.

D. Evaluation Indices
We used the normalized root mean square error (NRMSE), peak signal-to-noise ratio (PSNR), and structural similarity index (SSIM) to evaluate the reconstruction quality, and the reconstruction time were recorded.
The NRMSE is used to scale the RMSE to the range of (0, 1) and indicates the difference between the reconstructed image and the ground truth.The RMSE is defined as follows:    indicates better image quality.The PSNR is described as M AX I denotes the maximum intensity of the image.The SSIM, which considers brightness, contrast, and structure, was used to evaluate the similarity between two images and was calculated as follows: Here, µ I and σ 2 I are the mean and variance of an image, respectively.σ I R and σ I G are the covariance of I R and I G , respectively.

A. Reconstruction Comparison
The reconstruction results for the simulated phantoms are shown in Fig. 3, and the corresponding quantitative results are listed in Table I.Considerable background noise exists in the reconstruction results of the Tikhonov regularization.Moreover, the edges of the phantom were blurred, making it difficult to identify sharp edges.
For the NFL and proposed methods, the edges in the reconstructed images were sharper, and the intensity distribution inside the phantom was more uniform because of the TV norm.However, our proposed method can achieve better  reconstruction than the other methods, and it yields the most accurate recovery of the shape of the phantom among the three methods.Furthermore, as shown in Table I, the SSIM, PSNR, and NRSME of the reconstructed results obtained using the proposed method are superior to those of the results acquired by utilizing other methods.In addition, the comparison of the intensity distributions on the profile shows that the NFL regularization clearly underestimates the intensity compared to the ground truth, whereas the proposed method can could compensate for this drawback better, as shown in Fig. 3(b).

B. Error Comparison
To assess the ability to recover the intensity more accurately, we also reconstructed a multigray elliptical phantom with different intensities in each elliptical region, as shown in Fig. 2. The results are presented in Fig. 4 and Table II.In Fig. 4, the NFL method and proposed method both maintain the edges clearly.However, the error map in Fig. 4 reveals the inferiority of the NFL method.The intensity inside the ellipses is underestimated owing to the bias effect of the l 1 -norm.A clearer error map obtained using the proposed method indicates that this approach can provide a more accurate reconstruction.The relative errors of the reconstructed results and ground truths were calculated and are presented in II.The error values were reduced to within 5%.To evaluate the ability to recover the grayscale values more accurately, we focused on an elliptical phantom with different With SNR of 25 dB, 20 dB, 15 dB, we used Tikhonov, NFL, and the proposed regularization to reconstruct the signal and plot the error map, as depicted in Fig. 4. In addition, the average intensity in each elliptical region and the relative error were calculated and are presented in Table II.The error maps and the results in Table II reveal that although the NFL method can reconstruct the edge information more clearly, it underestimates the value inside each ellipse, which undermines the quantitative nature of MPI imaging, and its relative error is even greater than that of the Tikhonov method.The clearer error maps of the proposed method show that our method can provide more accurate reconstructions.The error can be suppressed to within 5%, thereby maintaining better quantitative properties.

C. Open-MPI Data Experiment
According to the results of our simulation experiments (Fig. 3 and Fig. 4), although the reconstruction results based on the Tikhonov penalty have the problem of over-smoothing, the intensity recovery is relatively accurate.Therefore, the intensity reconstructed using the Tikhonov penalty was used as a reference.
Fig. 5(a) shows the reconstruction results for the shape phantom in the Open-MPI dataset, and Fig. 5(b) displays the intensity distribution along the profile.For the Tikhonov regularization, the reconstruction results tend to converge to a maximum intensity of a few pixels and then become smooth (particularly in Profile 2 in Fig. 5(b)), which does not correspond to a uniform distribution of magnetic particles in the actual phantom cavity.For the NFL and proposed regularization method, as shown in Fig. 5(b), a stable region of intensity exists in the middle part along profile 2, which matches the actual situation.However, if the Tikhonov-based reconstruction results are taken as an approximation of the ground truth, Fig. 5(b) clearly shows that NFL regularization provides images with a significant intensity underestimation problem.6(b) shows that the reconstruction results of the NFL method significantly underestimate the reconstruction intensity at points P2 -P4.And, as shown in Table III, the signal intensities of P2 and P3 in the results of the NFL method are close, which does not match the concentration ratios in the Open-MPI dataset.In contrast, the results of the proposed method reduce the relative error (shown in Table III), and the concentrations of the P2 and P3 points show a gradient decreasing trend, which is consistent with the setup in the Open-MPI dataset.

D. Measurement Data Experiment
For the data measured by the homemade MPI scanner, frequency points were selected using the SNR of each frequency [38].Frequency points with an SNR > 10 dB were selected, and the number of filtered frequencies was 29.The images reconstructed using different regularizations are shown in Fig. 7(a).For Tikhonov regularization, the noise suppression is insufficient, leading to unclear reconstruction results.For the NFL and proposed regularization methods, two stripe distributions of magnetic particles can be clearly distinguished in the image.However, in the quantitative comparisons shown in Table IV, the concentration ratios of the two bars reconstructed by the NFL method deviate from the Ref value, whereas the ratios of the proposed method are closer to the Ref value and the Tikhonov reconstruction results, which illustrate that the proposed method enables quantitative reconstruction with better noise suppression and shape recovery.

VI. DISCUSSION
The ill condition of the system matrix of the MPI system renders the direct inverse of the linear equation unstable, and several regularizations have been proposed to solve this problem.The commonly used regularization approach, the Tikhonov penalty, provides a fast solution but leads to an oversmoothing problem.NFL regularization provides clear edge information; however, the l 1 -norm and l 1 -based near-isotropic TV norm used in NFL are biased estimates.Our experimental findings demonstrate that the reconstruction results obtained using NFL regularization significantly underestimate the intensity.
We propose a novel, nonconvex regularization method that can compensate for the concentration underestimation caused by NFL regularization and can precisely recover the intensity of a reconstructed image.In the proposed hybrid regularization, the nonconvex MC penalty was used to suppress the intensity bias produced by the l 1 penalty, and the TV penalty was added to ensure that the concentration distribution in the results was uniform.To avoid complex nonconvex objective problems, we chose MC regularization, which can be solved quickly using the firm threshold when the convexity condition is satisfied.Faced with the hybrid regularization of nonconvex MC and the isotropic TV penalty, the ADMM algorithm was employed to solve the problem; this algorithm transformed the problem into several simple subproblems for iteration.LU decomposition and the Chambolle projection method were adopted for acceleration in the subproblem solving.
We verified the validity of our proposed regularization method using simulated data, Open-MPI data, and homemade MPI data.A stenosis phantom, a vessel phantom with plaque in the main vessel, and overlapping ellipses of different intensities were used in the simulation.Reconstruction experiments were conducted using three phantoms with an SNR of 30 dB.The effectiveness of the method was verified through qualitative and quantitative comparisons.The accuracy and stability of the intensity recovery method were verified using overlapping elliptical phantom signals at different noise levels.The phantoms in the Open-MPI dataset and the data measured by a homemade MPI scanner were utilized to verify the generality of the proposed regularization on real-world data.
Although the proposed method enables more accurate image reconstruction, it may also introduce larger errors at the edges (see the error maps in Fig. 4).This may be due to the greater sparsity of the MC regularization we used, which is similar to concentrating a large range of errors into some pixel points at the edges.
The two-step method was used in our approach to simplify the selection of multiple parameters for the objective function with hybrid constraints.The noise related parameter ε was selected first, while β was selected using the adaptive updating strategy [35].In the update strategy shown in (27), the smaller the update threshold κ, the smaller the value of β.In our experiments, the data acquired under the same conditions used the same κ: κ was set to 10 for the simulation data, and κ was set to 80 for the Open-MPI data.This non-specific setting may lead to the selection of non-optimal parameters and adversely affect the reconstruction results, such as the larger error at the edge pixels in Fig. 4 and the reconstruction results of the proposed method in Fig. 6, where some background noise exists near P1.While the quality of the reconstruction results can be improved by fine tuning, this approach is incompatible with the intent of simplifying the selection of parameters.A more refined multiparameter selection method may optimize this contradiction.
The point P4 in Fig. 6, which has the lowest concentration, cannot be reconstructed by either the proposed method or the NFL method.The proposed method suppresses the weak signals from low-concentration points while providing sufficient noise suppression, reducing the dynamic range of the reconstruction method.Therefore, in practice, the appropriate regularization term and reconstruction method need to be selected according to the specific needs.Meanwhile, highquality reconstruction optimization methods that incorporate existing high dynamic range reconstruction methods remain to be investigated [39].
In Table VI, a comparison of the reconstruction times shows that although the proposed method achieved highquality reconstruction, the reconstruction speed is still slower than that of the NFL method, which is a limitation of the proposed method.We consider that the NFL method achieves faster reconstruction because it converts the 2D problem into a 1D problem when solving the objective function [13].However, we consider the runtime of the proposed method to be acceptable in non-real-time reconstruction scenarios.In addition, the current work is applicable to 2D reconstruction.In future research, the proposed method can be extended to 3D conditions, which will mainly involve improving the Chambolle projection algorithm presented in (25).Hence, the proposed methodology still needs to be optimized.Recently, deep learning-based methods have made progress in MPI imaging [40], [41], [42], [43], [44], [45], [46], [47].Through the use of substantial volumes of training data, learning-based methods can obtain accurate data distributions and priori information to achieve high-quality signal denoising [44], [47], system matrix optimization [42], [45], and image reconstruction [41], [46].The design of regularizations in physics-driven reconstruction can optimize the loss function via deep learning [48].Therefore, in future work, the proposed nonconvex regularization can be applied to optimize the reconstruction results of deep learning methods.In addition, direct comparisons between deep learning methods and physics-driven methods have yet to be performed, including the quality of the reconstruction results, the computation time, and the generalizability and robustness of the algorithms.

VII. CONCLUSION
In this study, we propose a regularization method that contains a nonconvex MC function and uses the ADMM algorithm to solve the optimization problem.During the iterative solution process, a subproblem containing nonconvex functions can maintain convexity and can be solved quickly by employing a threshold function.Reconstruction experiments on the simulated data, Open-MPI dataset, and the data measured by the homemade MPI scanner show that the proposed method outperforms the traditional Tikhonov and NFL regularizations and can alleviate the intensity underestimation problem existing in NFL regularization.The results demonstrate that the developed approach can achieve superior image quality while maintaining the quantitative advantages of MPI.

Fig. 1 .
Fig. 1.Comparison of the l 0 -norm, l 1 -norm, and MC penalty functions.(a) Comparison of the l 0 -norm, l 1 -norm, and MC penalty; (b) comparison of the 1D reconstruction results for l 1 , the MC and the corresponding TV hybrid penalty.The simulated data were used in experiments and the phantom have concentration varied in the middle region.The ordinate indicates the particle concentration in linear and arbitrary scales.The reconstruction results show that the MC penalty is effective at suppressing the bias introduced by the l 1 penalty.Furthermore, by adding TV regularization, the concentration distribution of the results becomes more uniform.For both hybrid reconstructions in (b), the parameter of the TV regularization is 0.02.

Fig. 2 .
Fig. 2. Phantoms used in our study.(a) Phantoms used in the simulation study, including the stenosis phantom (left), the vessel phantom with a plaque in the main vessel (middle), and the ellipses with different intensities (right, and the intensity in A1:1; A2:0.8;A3:0.6).(b) Phantom used in the real-world experiment.Bar phantom (left) and its schematic (right) used for the measurements (unit: mm).
) I R and I G denote the reconstructed and ground-truth images, respectively.N R M S E ∈ [0, 1].A smaller value Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

Fig. 3 .
Fig. 3. Reconstruction results for phantoms.(a) Reconstructed images of the stenosis phantom (top), vessel phantom (middle), and elliptical phantom (bottom).(b) Intensity distribution along the profile.The profile of the stenosis phantom is column25, that of the vessel phantom in column10, and that of the elliptical phantom in column 21. (SeeTable I for a quantitative comparison, Table V for regularization parameters and Table VI for the computation time.)

Fig. 4 .
Fig. 4. Visualization reconstructions and error maps of the elliptical phantom with different intensities.(SeeTable II for a quantitative comparison, Table V for the regularization parameters and Table VI for the computation time.)

Fig. 5 .Fig. 6 .
Fig. 5. Reconstruction of phantom shape using Open-MPI data.(a) Visualization results of the phantom shape and (b) intensity distribution along the profile given in (a).The data in (b) were interpolated.(SeeTable V for the regularization parameters and Table VI for the computation time.)

Furthermore
, although NFL regularization provides clear edges, the reconstructed image has only a central part, with the pixels filled by particles.Pixels that were not completely filled with particles were discarded during reconstruction.Our proposed regularization can suppress the intensity underestimation problem and reduce the over-sparsity problem of NFL regularization.Moreover, the central part of the results reconstructed using our regularization method exhibited a uniform intensity distribution (Profile 2 in Fig.5(b)).Fig.6(a)shows the reconstruction results of the concentration phantom in the Open-MPI dataset, and Fig.6(b) displays the intensity distribution along the profile that passes through four different concentrations (P1, P2, P3, and P4).To compare the concentration gradients reconstructed using different penalties more clearly, the intensity along the profiles was

Fig. 7 .
Fig. 7. Reconstruction was performed the data measured by a homemade MPI scanner.The intensity of the reconstructed images was normalized to [0,1].(a) Reconstructed results, where the dashed red boxes in the images indicate the reconstructed rods (b) Intensity distribution along the profile given in (a), where the arrow shows the direction.(SeeTable IV for the quantitative comparison, Table V for the parameters used and Table VI for the computation time).
Table I for a quantitative comparison, Table V for regularization parameters and Table VI for the computation time.)

TABLE I QUANTITATIVE
RESULTS OF THE PHANTOMS USED Table II for a quantitative comparison, Table V for the regularization parameters and Table VI for the computation time.)

TABLE II MEAN
VALUE AND RELATIVE ERROR OF THE DIFFERENT AREAS IN THE ELLIPTICAL PHANTOM Table III for the quantitative results, Table V for the regularization parameters and Table VI for the computation time).
Table IV for the quantitative comparison, Table V for the parameters used and Table VI for the computation time).

TABLE IV QUANTITATIVE
RESULTS OF THE BAR PHANTOM TABLE REGULARIZATION PARAMETERS FOR THE RECONSTRUCTED RESULTS