Solving the MPI reconstruction problem with automatically tuned regularization parameters

In the field of medical imaging, magnetic particle imaging (MPI) poses a promising non-ionizing tomographic technique with high spatial and temporal resolution. In MPI, iterative solvers are used to reconstruct the particle distribution out of the measured voltage signal based on a system matrix. The amount of regularization needed to reconstruct an image of good quality differs from measurement to measurement, depending on the MPI system and the measurement settings. Finding the right choice for the three major parameters controlling the regularization is commonly done by hand and requires time and experience. In this work, we study the reduction to a single regularization parameter and propose a method that enables automatic reconstruction. The method is qualitatively and quantitatively validated on several MPI data sets showing promising results.


Introduction
Magnetic Particle Imaging (MPI) is a tomographic imaging technique with huge potential in non-invasive diagnostic imaging (Gleich and Weizenecker 2005).Using various magnetic fields, it directly measures the distribution of superparamagnetic tracer in a three-dimensional region of interest.The method is sensitive with respect to tracer concentration and enables imaging at a very high temporal and spatial resolution without ionizing radiation (Gleich andWeizenecker 2005, Knopp et al 2017).After measuring a voltage signal in several receive coils, the particle distribution can be determined by solving an inverse problem based on a measured system matrix (Knopp et al 2017).This system is commonly solved using a regularized least squares approach and an iterative solver (Knopp et al 2010).Regularization is an important ingredient of this process, since it strongly influences the quality of the solution (as it can be seen in figure 1).Using an iterative solver, the amount of regularization mainly depends on three parameters: the weighting of the regularization term itself, the frequency selection of the voltage signal and the number of total iterations (Boberg et al 2021).For the reconstruction of a noisy measurement, the regularization needs to be strong.This is achieved by using only frequencies with good signal-to-noise ratio (SNR), a highly weighted regularization term and only few iterations.This comes at the cost of a lower resolution.For clear data, only a weak or no regularization at all is needed.Thus, the vast majority of frequency components should be used for a reconstruction with more iterations to ensure a sharp image.Generally speaking, the choice of these reconstruction parameters depends on the quality of the measurement itself, i.e. the scanner topology and the measurement settings, as well as on shape and iron concentration of the phantom.When considering time-dependent MPI data, e.g. the circulation of a tracer bolus through the cardiovascular system, the concentration may change over time.Thus, different parameters are needed for different stages of the bolus.Moreover, in the clinical context, there is a maximum iron dose for safe usage of superparamagnetic tracers in humans (Southern and Pankhurst 2018).Thus, a high concentrated tracer circulating constantly through the body for a longer time period is not suitable.New concepts like negative contrast perfusion try to tackle this problem (Mohn et al 2023).This underlines the importance of getting good reconstruction results on preferably low iron concentrations and thus the relevance of finding a good amount of regularization at each stage of a bolus.
Algorithms on the choice of regularization strength for general imaging inverse problems and especially magnetic resonance imaging are broadly studied in academic research.There are different strategies to emphasize: the classical methods based on the L-curve and the condition of the system (Hansen and O'Leary 1993, Lin et al 2004, Ying et al 2004), SURE-based methods (Ramani et al 2012), Bayesian approaches (Saquib et al 1998, Chaabene et al 2020) and learning-based methods (Toma and Weller 2020).However, all of these methods tackle only the regularization term and not the choice of the frequency selection.Up to now, only the L-curve approach has been adapted to MPI (Knopp et al 2008).Considerations on an automatically chosen SNR threshold have been done for single sided MPI (Volkens et al 2022).Furthermore, there has been proposed a more efficient frequency selection method based on the energy spectral density recently (Zhang et al 2023).
However, it remains a major challenge for each single MPI reconstruction to identify a good choice of parameters, leading to an ideal image.This is due to the complexity of the iterative solving process of the illposed inverse MPI reconstruction problem with different parameters not only controlling the regularization term but through the frequency selection also the system itself.The current state-of-the-art modus operandi is still trial-and-error as described in Weizenecker et al (2007), requiring much experience in MPI reconstruction.In order to establish MPI in medical imaging and, in particular, taking into account medical feasibility, an easy reconstruction leading to good results would be most beneficial.The methods presented in Knopp and Hofmann (2016) enable an online reconstruction for MPI, however all settings for the regularization must be done by hand.
Motivated from this practical point of view, we present a way to reduce the MPI reconstruction parameter set to a single parameter by adapting the SNR threshold and regularization weight in each iteration.Based on this, we propose a hands-free reconstruction method using a stopping criterion based on the quality of the data.A validation of the method is shown on various MPI data: dilution series with two different phantoms and an in vivo bolus experiment on the cardiovascular system of a mouse.
This paper is the full and extended version of our conference abstract (Scheffler et al., hands-free reconstruction for MPI, IWMPI 2023).The method has been further developed and the focus has been broadened to include sound theoretical descriptions of the proposed methods.The evaluation of additional data, including in vivo data, and the careful discussion of the method are also substantial additions to this manuscript.

Basics of MPI image reconstruction
In MPI several magnetic fields are used to derive the three-dimensional density distribution of superparamagnetic iron-oxid particles (SPIOs) inside a region of interest.The particle magnetization gets excited using an excitation field, which is spatially homogeneous but sinusoidal in time.A spatial resolution is achieved by adding a time-constant spatial gradient field bringing the SPIOs in saturation everywhere but a small low-field region.The nonlinear magnetization response of the SPIOs in the low-field region induces a voltage is the system function with support W Ì  . 3The system function is typically measured on a spatial grid with equidistantly placed grid-points r n for n ä I N = {1, K, N} by using a delta-probe.In combination with the discrete sampling of the voltage signal leading to a finite number of frequency components k ä I K = {1, K, K}, the discrete MPI equation can be given in matrix-vector form as with h Î  K is the measured voltage vector with additive background noise.Solving (1) for the particle concentration vector poses an ill-conditioned inverse problem, which is commonly adressed by formulating the Tikhonov regularized least squares problem with the parameter l Î +  controlling the influence of the regularization term.Additional real and non- negativity constraints are beneficial (Weizenecker et al 2009, Knopp et al 2017).Using an iterative solver (e.g. the Kaczmarz method), one ends up with the reconstructed particle distribution.
In the course of the reconstruction in MPI, the regularization is controlled by three main parameters as described in Boberg et al (2021):  , controlling the strength of the regularization term in (2).Bigger values for λ lead to a stronger regularization.
(ii) The frequency selection, determined by a minimum frequency k min and an SNR threshold Q Î +  , such that only frequency components in where SNR k denotes the SNR of frequency component k, are considered for reconstruction.By excluding noisy frequencies, a high threshold Θ leads to a small frequency selection and a weak noise amplification.Thus, the results get less noisy at the cost of loosing resolution.Due to the distribution of high-SNR frequencies around higher harmonics of the excitation frequency (Rahmer et al 2012, Szwargulski andKnopp 2017), the connection between Θ and the number of rows in the linear system is not linear.
(iii) The number of iterations i Î  used by the iterative solver.Stopping the reconstruction early is equivalent to a further regularization and less noise amplification (Fleming 1990).
The minimum frequency is mainly for excluding frequencies close to the excitation frequency, which are hardly measured correctly due to hardware limitations (Them et al 2016, Thieben et al 2023).In the following, we will drop further considerations on the minimum frequency and consider the frequency selection only to be determined by the SNR threshold Θ.The effect of Θ to the system can be described using projections from K to K Θ frequency components.Thus, we end up with the final reconstruction problem as The regularization of this reconstruction process is described by the parameter set , , .

Projection to a single parameter
For the standard Kaczmarz reconstruction, the parameter set  is defined a priori and has to match the noise level of the used MPI scanner as well as the measurement settings.Important scanner-dependent factors are properties of the receive coils as well as the overall quality of the receive path, which includes the receive coils and all electronic components involved in voltage processing and digitization, e.g. a low-noise amplifier (Graeser et al 2020, von Gladiss et al 2020).Moreover, the amount and concentration of the tracer directly determines the SNR of the measurement and thus the amount of needed regularization.To achieve a reconstruction result of good quality without any user input, an appropriate parameter choice has to be derived automatically from the measurement data.To this end, the first goal is to condense  into a single parameter, i.e. mapping the regularization strength λ and the SNR threshold Θ to ι.Thus, in each iteration there are different but consistent values for λ and Θ. Starting with high values in the first iterations, these reconstruction steps only work on few high-SNR frequencies with a strong regularization.For measurements with low signal intensity, the reconstruction can already stop at this point to achieve an artifact-free image with low noise.For measurements with high signal intensity, more iterations with a higher amount of used frequencies and less regularization are needed.Thus, low values for Θ and λ are chosen for later iterations.The transition in between should be a steady and smooth decrease.These properties are satisfied by functions of the form and j ä {Θ, λ}.Here, α j controls the starting value in the first iteration, and β j as well as γ j control the behavior and slope of decrease for the transition of the parameters in later iterations.For h Θ also a minimum SNR threshold Q Î +  min is needed, to exclude pure noise frequencies also for late iterations.It is important to note, that the values of α λ , α Θ and Q min are dependent on the total noise level and therefore on the receive path of the MPI scanner.Thus, it is sufficient to determine these parameters once for a certain receive path.Afterwards they can be used for all measurements using this receive path.
When altering the scanner hardware (e.g.changing the receive coil), an empty measurement u f empty over several frames empty empty can be used to calculate the mean noise h via Then, the parameters α λ , α Θ and Q min , which depend on the noise level of the scanner, can be estimated by a scaling of h.Since short empty measurements before each signal measurement are commonly done for background subtraction, auto-tuning the parameters from this empty measurement (as described in section 3) does not pose additional work.

Stopping criterion
After the described parameter projection, the number of iterations ι of the iterative solver remains to control the amount of regularization.For an automated reconstruction a suitable stopping criterion is needed.When the measurement is noisy, this stopping criterion should take effect already after few iterations.For less noisy measurements, less regularization is needed and the stopping criterion should not take effect too early, allowing more iterations.Typical stopping criteria for iterative solvers are build upon the L-curve between the 2-norms of the residuum r = Sc − u meas and c.The usage of the regularized Kaczmarz method poses an advantage, since the auxiliary vector Î  . In order to take into account the change in system size and regularization in each iteration, it is important to normalize r λ with respect to the amount of considered frequencies in each iteration.Then, a stopping criterion based on a jump in the curvature Y Î  of the L-curve (Knopp et al 2008) can be used.It is meaningful that the curvature should exceed a certain level dependent on the size of the solution norm, to hinder the reconstruction from stopping too early.In conclusion, we end up with the following stopping criterion at iteration step Î  i : , where c 1 is the particle concentration vector after the first iteration.The values for δ and ε describe the lower thresholds for curvature and change of curvature at which the stopping criterion takes effect.

Bolus reconstructions
Because of the high frame rate of typical MPI measurements, it can be assumed, that the tracer concentration does not change abruptly between two frames of time-dependent MPI data, e.g.bolus experiments.Thus, the number of iterations of the hands-free reconstruction should also not change too rapidly between two frames.
To exploit this additional form of prior knowledge and enforce this behavior, we propose the following two-step flattening approach.
(i) The number of iterations is not only dependent on the stopping criterion but also on the number of iterations of the previous frames.This value defines for each frame Î  f a number of iterations to be expected i Î  f exp .The actual number of iterations must not exceed a certain margin around the expected number i f exp .When the stopping criterion is triggered before the margin around the expected number of iterations, more iterations are performed until this margin is reached.When the stopping criterion is not triggered as the highest number of iterations in the margin is reached, the reconstruction stops at this point.The margin should be big enough to allow a gradually adjustment of the number of iterations over the course of the frames, but low enough to hinder the number of iterations to jump back and forth.
(ii) After a first reconstruction using the stopping criterion as well as the additional constraints of step (i), a second reconstruction is performed to handle remaining peaks.In this second round, the number of iterations over the course of the frames gets smoothed using a rolling average of the number of iterations from the first step (i).

Grid search reconstruction
For validation of the proposed method, a quantitative comparison to a reference solution is needed next to subjective statements on the reconstructed images.To this end, a three-dimensional parameter grid ( ) 1 2 3 1 2 3 1 2 3  is introduced.Note that as mentioned before, the connection between the parameter values and the amount of regularization is nonlinear.To span over a wide enough space of regularization strength, the parameter grid needs to have higher order exponential spacing in each component.A costly brute-force approach with 8 3 reconstructions is used to find the optimal parameter sets compared to a ground-truth phantom in two different measures.The normalized root mean squared error (NRMSE) was chosen because of its establishment and interpretability.Because the NRMSE does have issues when it comes to noise, we also consider the structural similarity index measure (SSIM) as a perception based measure (Wang et al 2004).For both measures, the reconstruction results get normalized to the interval [0, 1].We ranked all reconstructions from the parameter grid in both measures, respectively.The reference reconstruction is then chosen as the reconstruction with the parameter set with the lowest additive rank over both measures.

Study design
The proposed hands-free reconstruction was implemented by adapting the regularized Kaczmarz method from the MPIReco.jlframework in Julia (Knopp et al 2019).The general structure of the reconstruction is given as pseudo code in algorithm 1.Here, dcdr i is the approximate derivative of the solution norm at the iteration Î  i 0 as an intermediate step to calculate the approximate curvature Ψ i .We suggest a zero-initialization for the initial values.The function getParameter() calculates λ i and Θ i at iteration Î  i by using the mapping functions h λ and h Θ .
The full code is published under a Creative Commons Attribution 4.0 International license and can be found at https://github.com/IBIResearch/HandsfreeReco.Algorithm 1. Hands-free reconstruction.
, , dcdr , , dcdr , Phys.Med.Biol.69 (2024) 045024 K Scheffler et al To verify the method on a broad variety of MPI data with different iron concentrations, it was tested on the dilution series described in Boberg et al (2021).Moreover, the bolus study describing the blood-flow of a mouse heart from Graeser et al (2017) was used to test the hands-free reconstruction on in vivo data.All studies have been performed on the preclinical MPI system 25/20FF (Bruker Corporation, Ettlingen, Germany), however the bolus experiment used a dedicated receive coil with better SNR.Based on experience on MPI reconstruction with the Kaczmarz method, the parameters describing the slope of the mapping functions were chosen as β λ = 0.2, γ λ = 5 and β Θ = 0.28, γ Θ = 2.The remaining parameters dependent on the noise level of the receive path, were calculated from the mean noise of the data sets h (derived from empty measurements, as described in 2.2), such that they work for all data sets on all concentrations: .

Dilution series
The dilution series is taken from (Boberg et al 2021) and was recorded using a three dimensional Lissajous trajectory and a gradient strength of (−0.75, − 0.75, 1.5) T m −1 .The drive-field amplitudes were set to 12 mT in each direction, resulting in a FOV of 32 × 32 × 16 mm.The system matrices were measured using a 2 × 2 × 1mm delta-sample filled with perimag (micromod Partikeltechnologie GmbH, Rostock, Germany) diluted to 300 mmol l −1 on a 21 × 21 × 24 grid.The single dot phantom series consists of 12 measurements where the iron concentration is halved for each new measurement starting with k = 400 1 dot mmol l −1 in a small glass capillary with an inner diameter of 2.4 mm.The dot phantom was only filled with a small amount (20 μl) of tracer.Furthermore, we consider the dilution series on the rat kidney phantom.Two kidneys from a realistic rat phantom (Exner et al 2019) are each filled with 781 μl tracer.The iron concentration is halved 5 times, starting with k = 3.05 1 kidney mmol l −1 .For the reference reconstruction as well as the error analysis a ground truth was voxelized out of the CAD files of the dot phantom and the kidney phantom.The 8 3 sized parameter grid for the reference reconstruction was derived proportional to higher order exponential functions giving 1.4 10 , 8.8 10 , 4.5 10 , 2 10 , 7.6 10 , 0.28, 1. 1, 7.3 , 1.32, 1.64, 2.24, 3.4, 5.88, 11.71, 27.33, 76 , 1, 2, 3, 5, 7, 11, 16, 25 .
The grid searches were performed in parallel on 8 kernels and took about 6 to 8 h for each concentration on both phantoms resulting in a total net reconstruction time of about 5.25 d.

In vivo bolus experiment
In addition to phantom measurements we consider the in vivo measurements from (Graeser et al 2017), which were measured with a dedicated high-sensitivity receive chain that shows a very different frequency response and noise characteristic than the standard receive coils used in the previous phantom experiments.The measurements were performed using a finer system matrix dedicated to the fine structures in the cardiovascular system of a mouse.The cylindrical delta probe had a diameter as well as a height of 0.7 mm, filled with undiluted LS-008 ((Khandhar et al 2017), 91.7 mmol l −1 ).The system matrix was measured on a 46 × 36 × 19 grid.The measurement of the mouse itself used three separate boli of LS-008, sequentially injected into the tail vain.In this work, we only consider the second bolus with a volume of 10 μl concentrated at 9.17 mmol l −1 .The concentration strongly changes in the course of the bolus throughout the cardiovascular system of the mouse.Thus, different regularization strength in different stages of the bolus are needed for good image reconstruction.
The measurement was performed on a 32 × 32 × 16 mm 3 sized FOV covering the mouse heart and took 21.54 ms per frame, resulting in a temporal resolution of 46 Hz.For evaluation we consider a passage of the measurement over 115 frames in total.
No averaging on the data was performed to preserve the temporal resolution.This leads to a stronger influence of noise and thus may lead to an erratic number of iterations in the course of the frames.Since we know, that the bolus does only need a gradually changing strength of regularization, we want to enforce a smooth transition in the reconstruction results between the frames.To this end, the number of iterations of the previous frames is also considered next to the plain stopping criterion.Let i Î  f be the number of iterations of the hands-free reconstruction at frame f ä {4, K, 115}.To get the number of expected iterations for ι f , we weight the latest numbers of iterations via Then, the number of iterations ι f for each frame is restricted to fulfill This rather large margin enables the algorithm to easily adjust the number of iterations over the course of the frames, but still allows large jumps in the number of iterations between two frames.Therefore, the course of the number of iterations over the frames is further smoothed in a second reconstruction step, by using a rolling average.In other words, we set the maximum number of iterations for each frame in the second step as the average of the two previous and two subsequent frames from the first round.The result of this two-step flattening is exemplarily shown in figure 3.For anatomical reference, the reconstructed MPI data is overlayed on magnetic resonance imaging (MRI) data that was acquired with an isotropic gradient echo sequence using the preclinical 7 T device BioSpec 70/30 (Bruker Corporation, Ettlingen, Germany) as described in Graeser et al (2017).

Dot phantom
The reconstruction results of the dilution series with the dot phantom are shown in figure 4. On the left side (framed by the blue box), the reconstructions using the parameter projection are given, with different values for ι  as the only remaining parameter.It can be seen over the course of the dilution series, that the optimal number of iterations strongly differs with the concentration.For high concentrations, the best results are achieved after many iterations (ι = 22).Over the course of decreasing concentrations, the result after too many iterations (with the corresponding low regularization strength) starts to get noisy.Thus the best results are achieved after a decreasing number of iterations.Eventually, for the lowest concentrations, the best results are obtained after only one to three iterations.
The stopping criterion is triggered in a similar manner over the course of decreasing concentrations.Starting with 24 iterations for the highest concentration, the hands-free method stops after fewer iterations for each further dilution step.For the lowest concentrations, the method stops after 2 iterations.A comparison to the reference reconstruction resulting from the grid search is shown on the right side of figure 4. By using the grid search regularization parameters we are able to produce even sharper images, however, this comes with slightly enhanced noise.All grid search parameters are given in table 1 next to the resulting NRMSD and SSIM values compared to the voxelized phantom.
Considering the NRMSD results of the hands-free method, the lowest error values are achieved for the highest concentrations and are very close to the NRSMD values of the grid search reconstruction.Over the course of decreasing concentrations, the NRMSD values get worse and are significantly higher than the NRMSD values of the reference images resulting from the grid search.For the lowest concentrations the NRMSD values corresponding to hands-free method and grid search are closer to each other again.The SSIM value of both methods are similar to each other for all concentrations.On the intermediate concentrations, the SSIM values of the hands-free method are-in contrast to the NRMSD results-even higher than the SSIM values from the reconstruction resultung from the grid search.

Kidney phantom
The reconstruction results of the kidney phantom from the hands-free method as well as the results of the grid search are shown in figure 5. Similar to the dot-phantom, the stopping criterion is triggered at earlier iterations for lower concentrations.
In comparison to the reference image resulting from the regularization parameter set from the grid search (as given in table 2), the results from the hands-free method are qualitatively similar for each iteration.This is confirmed when considering the NRMSD as well as the SSIM of both methods in table 2. The deviation to the voxelized phantom is close between both methods in both measures over all concentrations.

In vivo bolus
For the evaluation of the hands-free method on the in vivo bolus, we consider figure 6. Shown is the hands-free reconstruction result column-wise compared to reconstructions using a strong and a weak regularization constantly over all frames.The bolus is displayed row-wise for 8 frames with equidistant time difference of 0.32 s in a single xz-slice.The total time distance between the first and the last shown frame is 2.25 s.By laying the reconstructed MPI signal on top of the corresponding MRI data, the flow of the bolus through the mouse can be followed.We expect to see the bolus entering the vena cava inferior in the first frames and evolve to the heart by entering the right atrium.Following the results of Graeser et al (2017), the bolus will then fill the right ventricle and lung.However, since we show only the slice, where the vena cava and the right atrium are visible, we expect the MPI signal to vanish.
We see that the strong regularized solution lacks of accuracy in the first frames (frames 5, 20, and 35), when the bolus enters the vena cava and flows to the atrium at high concentration.When the bolus gets weaker at later frames, its position is still identifiable but the reconstructed signal is very weak.
In contrast, the solution using weak regularization has high accuracy but also shows large artifacts and noise.For later frames, the bolus is even not distinguishable from the noise anymore.
The hands-free method is able to adapt the regularization strength frame-wise (with the number of iterations after which the stopping criterion is triggered superposed in white).It generates images with high accuracy and low noise for all shown frames.After entering the vena cava inferior in frames 5 and 20, the tracer enters the right atrium of the mouse heart, as it can be seen in the frames 35, 50 and 65.Furthermore, the vanishing of the tracer in the mouse heart in the later frames is clearly visible.For more details, a comparison to the images shown in Graeser et al (2017) is meaningful.However, it has to be noted, that Graeser et al use the bolus experiment with a ten-fold higher concentration.The hands-free method is able to give good results on the bolus progression at a lower concentration.
The total course of number of iterations over the frames is shown in the lower part of figure 6.It can be seen, that the number of iterations climbs (and therefore the amount of regularization shrinks), when the bolus enters the vena cava inferior at the first frames.Then, the number of iterations gradually decreases and the amount of regularization increases, as the bolus is weaker at later frames.

Discussion and outlook
Considering the presented reconstruction results of the dilution series of the dot phantom, we see that the presented regularization parameter condensation works as intended.The number of iterations is inversely proportional to the amount of regularization and poses a trade-off between blurred images and sharp but noisy images.Thus, for high concentrations, a higher number of iterations and for low concentrations a low number of iterations gives the best results.This itself is a very promising step towards simplifying the problem of choosing the right amount of regularization for system-matrix-based MPI reconstruction.
The choice of the parameters of the mapping functions α λ , α Θ and Q min is dependent on the scanner hardware.The step from choosing different regularization parameter setups for each measurement to choosing different mapping function parameters only one time for an MPI system is huge.Furthermore, a first idea on auto-leveling these parameters out of the noise data was given and tested on different receive paths.
In combination with the proposed stopping criterion, the hands-free method is able to find a reasonable amount of regularization for all 12 different concentrations of the dot phantom as well as the 6 concentrations of the kidney phantom and thus over a vast range of MPI data.It selects a good compromise between resolution and noise.With enough time and computation power, it is possible to find better solutions for each measurement by a traditional reconstruction with hand-chosen regularization parameters.However, this does not weaken the huge potential of the hands-free reconstruction, since it aims to find a good reconstruction result without any user input in a short time.This goal is achieved and the hands-free reconstruction enables automatic reconstructions without experience and prior knowledge on the data.Further validation of the method requires future application to different data sets and different MPI systems.A runtime analysis of the Shown is a solution using strong respectively weak regularization as well as the hands-free method on several frames with equidistant time difference.The number of iterations after which the stopping criterion was triggered is superposed in white.The results are normalized to the maximum of the hands-free solution c max HF , with the common colorbar given below.Furthermore, the number of iterations of the hands-free method over the course of the frames is given on the right side.
hands-free reconstruction was not part of this work.However, since the underlying iterative solver remains, the reconstruction times are comparable to the regular Tikhonov regularized Kaczmarz method, for which online reconstruction has been shown to be possible (Knopp and Hofmann 2016).
The hands-free methods also shows its benefits on the time-dependent in vivo data with varying concentration.Whereas the reconstructions using either strong or weak (static) regularization do not give good results over all frames as they are either good for the frames with low concentration or respectively with high concentration.The hands-free method is able to adapt the regularization strength over the course of the frames and does give superior results.In comparison to Graeser et al (2017), the hands-free method is able to give good results on a ten-fold weaker concentration.This is especially important with regard of the usability of MPI for human imaging, since too high tracer concentrations are a medical issue (Southern and Pankhurst 2018).
To hinder jumps in the course of reconstructed images, a dual approach with two rounds of reconstructions as described in 3.2 was used.This procedure is efficient to produce a smooth transition between the single frames.However, for an instant online reconstruction a method using only one round of reconstruction may be beneficial.In this case, the second round can be used for post-processing the data.
Combining the results of the hands-free method on the different data sets, the hands-free method poses a very promising step for MPI application in the medical context, where fast and user friendly solutions are important, especially considering time-dependent bolus measurements.The hands-free method enables good reconstruction results on a low concentrated bolus, which is important when considering the maximum amount of human compatible iron particles.Furthermore, the method can also be used to get a fast first view on new data.
All data were derived from the same MPI scanner, but different receive coils with differing frequency responses and noise characteristics.Thus, it was shown, that the method is capable to adapt to different MPI receive chains.However, it remains to show that the method can be directly transferred to other MPI systems with different exciting frequencies and gradient strengths.In doing so, the auto-leveling of the mappingparameters dependent on the receive path can get further investigated and improved.
Furthermore, it is meaningful to adapt different methods on hyperparameter tuning for inverse problems, like bayesian approaches or learning-based methods to the MPI setting.Adaptions on such methods would be necessary, to include the choice of the frequency selection.But especially combinations with the method presented in this paper are conceivable, for instance an advanced estimation of the parameters of the mapping functions.

Conclusion
In this work we investigated a possible simplification of the problem of finding the right amount of regularization for system matrix based MPI reconstruction.By condensing the regularization parameter triplet to a single parameter (the number of iterations of the iterative solver) and deriving a fitting stopping criterion, a method was presented, that is able to produce good reconstruction results in a fast manner.This can be achieved without additional user input, prior knowledge on the data nor additional measurements.The method was validated on several MPI measurements using different phantoms on a broad band of iron concentrations.Furthermore it was shown, that the method is very promising for in vivo bolus experiments with timedependent concentration.This hands-free reconstruction represents a crucial step towards medical usage of MPI systems.

Figure 1 .
Figure 1.The amount of regularization strongly influences the quality of the reconstructed images in MPI.
The resulting mapping functions are shown in figure 2. The parameters for the stopping criterion were set to d

Figure 3 .
Figure3.Enforcing a smooth transition between the frames of the in vivo data using the hands-free method is done by a two-step flattening.The number of iterations from the plain stopping criterion is shown in blue, the restriction by an equalization following (3) is shown in orange and the final flattened course is shown in green.

Figure 2 .
Figure2.Mapping functions for the parameter condensation on the different data sets (calculated from the mean noise which is derived from empty measurements).The left image shows the course of the regularization pre-factor λ over the iterations, the right image the course of the SNR threshold Θ. Note, that the y-axes are logarithmic.

Figure 4 .
Figure 4. Reconstruction results for the dot phantom at each concentration (k = 400 1 dot mmol l −1 ) in the xy-plane with a maximum projection on the central slices in z-direction.On the left side (blue box), the reconstructions resulting from the parameter projection with various number of iterations are shown.On the right side (yellow box), the hands-free reconstruction (the number of total iterations dedicated by the stopping criterion is superposed) and the reconstructions with an optimal parameter set resulting from the grid search are shown.All results are normalized to the range [0, 1], the colorbar is given next to the phantom.

Figure 5 .
Figure 5. Reconstruction results for the kidney phantom.Shown are the hands-free method (top, the number of total iterations dedicated by the stopping criterion is superposed) and the optimal parameter set resulting from the grid search (bottom) for each concentration (k = 3.05 1 kid mmol l −1 ) in the xy-plane on slice 10 in z-direction.All results are normalized to the range [0, 1], the colorbar is given next to the phantom.

Figure 6 .
Figure6.MPI results of the in vivo bolus visualized on top of corresponding MRI data.Shown is a solution using strong respectively weak regularization as well as the hands-free method on several frames with equidistant time difference.The number of iterations after which the stopping criterion was triggered is superposed in white.The results are normalized to the maximum of the hands-free solution c max HF , with the common colorbar given below.Furthermore, the number of iterations of the hands-free method over the course of the frames is given on the right side.

Table 1 .
Error results for the dot phantom.

Table 2 .
Error results for the kidney phantom.