Dual-Feature Frequency Component Compression Method for Accelerating Reconstruction in Magnetic Particle Imaging

The frequency component compression method (FCCM) has been widely used in magnetic particle imaging (MPI) technology to improve reconstruction efficiency. This method can reduce the reconstruction time by using the signal-to-noise ratio (SNR) feature to remove high noise frequency components. To further accelerate the reconstruction, a dual-feature frequency component compression method (DF-FCCM) was developed herein. A new energy spectral density (ESD) feature was introduced to describe the noise level of measurement signal. By using the SNR feature and the ESD feature, the DF-FCCM can select valuable frequency components that contain low noise both in the measurement signal and the system matrix. The reconstruction time can be reduced by using fewer frequency components that contain more valuable information. The efficiency and the robustness of the proposed method was validated through extensive simulation experiments. Further real experiments based on the OpenMPI data set verified that the DF-FCCM can be applied in real MPI reconstruction. Compared to the previous SNR-FCCM, the DF-FCCM can achieve similar or better reconstruction quality by using 25% reconstruction time. The proposed method can efficiently improve the reconstruction efficiency and potential to improve the online MPI imaging, which is essential for pre-clinical and clinical applications.


I. INTRODUCTION
M AGNETIC particle imaging (MPI) is a relatively novel molecular imaging technology [1]. It can image superparamagnetic iron oxide nanoparticles (SPIONs) non-invasively based on the movement of the field-free region and the non-linear magnetization response of SPIONs [2], [3]. Researchers have reported that MPI can in-vitro detect 156 μg/L iron concentrations [4], and provide positive contrast without background signals. It also supports sub-millimeter spatial resolution and millisecond temporal resolutions that allow for high quality real-time imaging [5], [6], [7]. MPI is a promising imaging technology and has been widely applied in various biomedical fields including cell tracking [8], cardiovascular imaging [9], blood pool imaging [10], and tumor detection [11].
In MPI technology, the desired particle concentration is recovered from the measurement signal. A variety of reconstruction methods have been proposed to reconstruct the particle distribution [12], [13], [14], [15]. System matrix reconstruction methods, such as the Kaczmarz method [16], provide high reconstruction quality and have been widely used. However, the large scale of the system matrix increases the calculation burden and reconstruction time [17]. This problem limited the efficiency of the system matrix reconstruction method.
To improve reconstruction efficiency, the signal-to-noise ratio based frequency components compression method (SNR-FCCM) was proposed to reduce the frequency component of the system matrix [18], [19]. This method can select hundreds of valuable frequency components, which means frequency components with high particle signal strength and low noise signal strength, for reconstruction. While the reconstruction time scales linearly with the number of frequency components, the SNR-FCCM can reduce the reconstruction time by using 5% to 10% frequency components of all frequencies for reconstruction. This method has been applied in many system matrix reconstruction researches and combined with different reconstruction methods [20], [21]. However, the reconstruction quality has to trade off against the reconstruction time. The image resolution can be rapidly decreased to 50% when using the SNR-FCCM to reduce the frequency components from 10% to 2.5% [22]. This problem limited the SNR-FCCM to compress the data to 10% or fewer.
To further improve the reconstruction efficiency, a dualfeature frequency component compression method (DF-FCCM) is proposed herein. A new energy spectral density (ESD) feature is introduced to supplement an evaluation of the measurement signal, which is essential to the reconstruction quality. In the DF-FCCM, the SNR feature and ESD feature are separately used to evaluate the quality of frequency components. The valuable frequency components, which contain more information and less noise, are selected and used for reconstruction. Thus, the DF-FCCM can accelerate the reconstruction by using less frequency components to achieve better reconstruction.
To evaluate the proposed compression method, extensive experiments are performed. The efficiency of the DF-FCCM was first evaluated. The performance of the DF-FCCM was compared with the DF-FCCM in different noise levels. Subsequently, the quantitative ability of MPI technology whether affected by the DF-FCCM or not was researched. Finally, to validate that the DF-FCCM can be applied in the real MPI data, we performed the real reconstruction experiments by using the OpenMPI data set [23]. These experiment results confirmed that the DF-FCCM is an efficiency compression method to accelerate the MPI reconstruction.
This paper is organized as follows: Section II presents the theory and the process of the proposed method. Section III discusses the experiments performed to validate the advantages of our method. Finally, the proposed method is discussed in Section IV, and concluding remarks have been provided in Section V.

A. System Matrix Reconstruction Method in MPI
The reconstruction method is necessary for MPI to recover particle concentrations from the measured voltage. The relationship between the particle concentration c(r) at position r ∈ Ω and the measurement signal u(t) provided in previous research [19], can be described as follows: In the linear system, s(r, t) is the integral kernel describing the mapping between the particle concentration and voltage signal. As the signal u(t) and kernel s(r, t) are both periodic continuous sequences, they can be expanded into Fourier series and discretized using a rectilinear rule with N integration points r n , n = 1, 2, . . . , N. The discrete frequency equation can be obtained:û The discrete (2) can be compactly written in a matrix sequence form as follows: where S ∈ C M ×N is the system matrix, u ∈ C M ×1 is the frequency measurement sequence and c ∈ R N ×1 is the particle concentration sequence. N represents the number of calibration points, which determines the number of pixels in the reconstructed image, and M represents the number of frequency components for the dimensions of the system matrix.
To solve the linear equation system (3), the Kaczmarz method is used in the work of Knopp et al. [24]. The iteration equation is the following: where η denotes a temporary variable, i = 1, 2, 3, . . . , I denotes the iteration number, z is the residual item, c is the particle distribution, s j is the jth row in the system matrix S, u j is the jth element in the measurement signal u and λ is the regularization parameter.

B. SNR Based Frequency Component Compression Method
In previous studies [18], [19], the SNR feature was proposed to quantify the noise in the frequency components of the system matrix. The SNR-FCCM was applied to reduce the system matrix after the calibration step.
In the SNR-FCCM, the SNR feature of each frequency component is determined by the calibration scans with and without a sample probe [19]. The SNR of a specific frequency f i can be calculated using the following equation: where i = 1, 2, . . . , M, S δ and S 0 denote the frequency component of the calibration signal with and without the sample probe, respectively, and N is the number of calibration points. After acquiring the feature sequence SNR ∈ R M ×1 through calibration measurement, the threshold value σ SNR is applied. The low noise frequency components are selected using the following equation: The dimension of the frequency component is reduced from M to M SNR and the selected frequency component is f SNR ∈ R M SNR ×1 . The system matrix S and the measurements signal u can be reduced to a small system matrix S SNR = S(:, f SNR ) and small measurement signal u SNR = u(f SNR ) respectively.
Further, a reduced linear equation system can be described as: The reconstruction time can be decreased for fast MPI imaging by solving this reduced linear equation system.

C. ESD Feature of the Measurement Signal
Because the SNR feature sequence is directly calculated using the system matrix S with Eq. (5), the noise contained in the measurement signal is not considered by this feature. Thus, the SNR-FCCM can not select the frequency components that contain more valuable information both in the system matrix and measurement signal. To address this problem and further improve the efficiency of the FCCM, the ESD feature is introduced herein to describe the noise contained in the measurement signal.
Based on previous research regarding discrete data sequences and Parseval's theorem [25], the distribution of sequence energy can be represented by the ESD feature as follows: is the discrete-time Fourier transform (DTFT) of the measurement signal. The measured frequency signal is practically acquired by coupling various types of noise signals into the particle signal as follows: whereû p andû n represents the particle and other noise signals, respectively. Thus, the ESD feature of the measurement signal can be described as follows: In our theory, the ESD feature will focus on the random noise. Thus, when the noise in (10) is assumed to the random noise, the energy of the noiseû n is uniform at each frequency component [26]. The ESD feature of noise can be described as follows: where L is the noise level of the measurement signal. From (10) and (11), the noise level of the measurement signal can be described with the following equation: Based on this equation, the ESD s can be used to evaluate the random noise level in each frequency component, while the exact noise value is not necessary. Thus, the frequency components can be directly selected based on the ESD feature using the following equation: This feature is beneficial to select valuable frequency components which contain low noise and valuable information in the corresponding measurement signal.

D. Dual-Feature Frequency Component Compression Method
Based on the ESD feature, a dual-feature frequency component compression method (DF-FCCM) is introduced. In the DF-FCCM, the SNR feature is used in the first selection step, and the ESD feature is used in the second selection step. The flow chart of DF-FCCM is performed in Fig. 1.(a) and the detail process of DF-FCCM is illustrated in Fig. 1.(b).
The first compression is applied after acquiring the system matrix in the calibration step. With the system matrix S, the feature sequence SNR ∈ R M ×1 can be calculated by (5). Then the feature sequence SNR is rearranged to form a new descending feature sequence SNR d which is satisfied by the following relationship: where I 1 is the index value of the descending sequence, M is the length of frequency component sequence and m 1 = 1, 2, . . . , M is the index number. Subsequently, the selected frequency sequence f r1 can be acquired by using the following equation: The parameter K 1 is the used frequency percentage in the first step. With the selected frequency f r1 , the system matrix S can be reduced to a new smaller system matrix S r1 = {S(:, f) | f ∈ f r1 } to save the storage space and improve the reconstruction efficiency.
The second compression is applied after acquiring the measurement signal. First, the signal is reduced using the selected frequency sequence f r1 to obtain the reduced signal u r1 = {u(f ) | f ∈ f r1 } corresponding to the system matrix S r1 . Then the ESD feature sequence of the signal u r1 is calculated using the following equation: The feature sequence ESD is rearranged to form a new descending sequence ESD d which is satisfied by the following relationship: where I 2 is the index value of the ESD feature sequence and m 2 = 1, 2, . . . , M × K 1 is the index number. Based on the feature sequence, the selected frequency f r2 can be acquired using the following equation: where parameter K 2 is the used frequency percentage in the second step. Using the selected frequency component f r2 , the system matrix is reduced to S r2 = S r1 (:, f r2 ) and the signal is reduced to u r2 = u r1 (f r2 ).
After the frequency components is compressed, the reduction linear equation system can be formulated as follows: where S r2 and u r2 is the compressed data. A proper reconstruction method, such as the Kaczmarz method, can be used to solve this reduced linear equation system with higher speed.

III. EXPERIMENTS & RESULTS
To evaluate the performance of the proposed method, extensive experiments were performed herein. First, the efficiency of the proposed method was evaluated by comparing it with the SNR-FCCM at different noise levels. Then, we validated that the DF-FCCM did not negatively influence the quantitative ability of MPI technology. Finally, the performance of the DF-FCCM in real data was evaluated. The real signals from the OpenMPI dataset were used to compare the performance of the DF-FCCM and the SNR-FCCM.
All reconstructions were performed on a computer with an Intel Core T M i7-6700 K CPU (4.00 GHz) and 64 GB RAM. All programming was coded by Matlab under Windows 10 environment.
A. Experiments Set 1) Calibration Data: Calibration data from the OpenMPI data set were used with the three-dimensional Lissajous scanning sequence and Bruker Preclinical MPI Scanner [23]. The field-of-view had a size of 37 mm × 37 mm × 18.5 mm, and the size of the gridding space was 37 × 37 × 37. A sample probe was fabricated using the Perimag tracer with a concentration of 100 mmol/l. The size of the sample probe was 2 mm × 2 mm × 1 mm.
2) SNR Feature: In the following experiments, both the simulation and real experiments, the SNR feature from the OpenMPI data set was used. While the SNR feature was acquired during system matrix calibration, the size of the SNR feature was same to the frequency component number of the system matrix.
3) Simulation Signal: To make the simulation closer to the real physics, we used the calibration data from the OpenMPI data set and following the work of Von et al. [6] to synthesize the signal. For the particle distribution, three different numerical phantoms were constructed: the stripe phantom, the stenosis phantom, and the letter phantom. The phantoms were shown in Fig. 2. The original size of phantoms was 37 × 37 voxels. For using the calibration data from the OpenMPI data set with the size of 37 × 37 × 37, the phantoms were enlarged to the same size by filling the void voxels with 0. The Gaussian noise was added to the simulation measurement signal before reconstruction.

4) Real Signal:
The real MPI signals from the OpenMPI dataset were used in real data experiments [23]. The measurement signals were converted into frequency signals using the fast Fourier transform and averaged over 1000 times. The frequency component was also processed using a band-pass filter to filter out the direct feed through. Only frequencies larger than 70 kHz were selected.

5) Reconstruction Parameters:
To fairly compare the efficiency of compression methods, we used regularization parameter λ = 0.001 and iteration number I = 10 with the Kaczmarz reconstruction method in all experiments [24].
6) DF-FCCM Parameters: There are two essential parameters: K 1 and K 2 in the DF-FCCM method. K 1 represents the used frequency percentage in the first step of DF-FCCM. It determines the frequency number for the following step. To balance the time cost on data compression and the loss of valuable frequency components, we applied K 1 = 10% as a fixed parameter in the following experiments. In this scenario, the exact SNR value is 0.1279. K 2 represents the used frequency percentage in the second step, which also indicates the final frequency number for reconstruction. In the following experiments, we used the 'Threshold' to describe the percentage of the frequency component for reconstruction, which was equal to the K 2 parameter in DF-FCCM. As an example, when the whole frequency number M is 76263, 76263 × 1% ≈ 762 frequency components were used for the threshold of 1%.

B. Evaluation Indicators
To evaluate the reconstruction results, the peak signal-to-noise ratio (PSNR) indicator, structural similarity (SSIM) indicator, full-width at half-maximum (FWHM) indicator and reconstruction time were used herein.
The PSNR was commonly used to evaluate the quality of images, and it can be calculated by the equation: where x a and x r represent the actual and reconstruction results, respectively. The SSIM indicator was used to evaluate the similarity between actual result and reconstruction result, and the equation is as follows: where c 1 and c 2 are the constant values. The FWHM indicator was used to evaluate the spatial resolution of results [27]. It was defined as the width at which a function decays to half of its maximum. For the reconstruction results, the profile along the center line of the image was used to calculate the FWHM. In general, a lower FWHM indicates a better spatial resolution.
Additionally, the reconstruction time was calculated to quantify the efficiency of the compression methods. We recorded the reconstruction time from reading the raw signal and system matrix to finish the final results. Both the time spending on processing data and reconstruction were taken into account.
To make the evaluation more robust, all reconstruction in the experiments were repeated 1000 times. In each reconstruction, the noises were different, and the evaluation was individual. To acquire reliable indicators, the evaluation indicators of all results were averaged as the final value.

C. Validating the Efficiency of the DF-FCCM
To validate the efficiency of the proposed method, the SNR-FCCM and DF-FCCM were applied to the letter phantom. The threshold was set from 0.5% to 4.5% with a 0.25% interval. Three Gaussian noise levels: 30 dB, 40 dB, and 50 dB were tested, respectively [28]. The statistical lines of all the indicators were performed in Fig. 3.
According to the statistical lines of indicators, it can be concluded that the DF-FCCM can achieve better reconstruction quality than the SNR-FCCM when using the same number of frequency components. In different noise levels, the DF-FCCM can robustly acquire better reconstruction quality both in the FWHM, PSNR and SSIM indicators. In 40 dB noise levels, the FWHM indicators increased from 5.56 to 8.33 when using the SNR-FCCM to compress the frequency from 4% to 1%. These results represented that the spatial resolution decreased nearly 50% by the SNR-FCCM. On contrary, both indicators were changed minimally when using the DF-FCCM to compress the frequency from 4% to 1%. These results proved that the DF-FCCM can efficiently compress the frequency while maintaining high reconstruction quality.
Further, the visualization of three specific thresholds (1%, 2%, and 4%) were performed in Fig. 4. Comparing the results of SNR-FCCM and DF-FCCM, the DF-FCCM performed better visualization. When use DF-FCCM compressed the frequency from 4% to 1%, the quality of results decreased minimally. As for the SNR-FCCM, the shape of results were distorted and the spatial resolution obviously decreased when compressing the frequency from 4% to 1%. The indicators of results with 1%, 2%, and 4% threshold were performed in Table I to quantify the efficiency of the DF-FCCM. In 30 dB noise level, the SNR-FCCM only acquired a low quality result with 4% frequency components. The SSIM indicator of this result was 0.37. In contrast, the DF-FCCM spent 25% time of SNR-FCCM to accomplish a higher quality result than the SNR-FCCM. Similar results can be found in other noise levels.
To sum up, it can be concluded that the DF-FCCM can efficiency accelerate the reconstruction in different noise levels. The DF-FCCM can use 1% frequency components to acquire better results than 4% frequency component results of SNR-FCCM. This improvement can efficiently accelerate the MPI reconstruction while maintaining the quality of results.

D. Validating the Quantitative Ability of the DF-FCCM
After validating the efficiency of the DF-FCCM, we further explored the quantitative ability of the DF-FCCM herein.
For MPI technology, the concentration of the SPIO tracer was generally consistent and accurate across the scans. Thus, the concentration is linearity with the reconstruction intensity. As a new signal processing method, we tested whether the linear quantitative ability of MPI results would be affected by the DF-FCCM or not.
Based on the stenosis phantom, we simulated five different concentrations by setting the value of c to 1 × 10 7 , 2 × 10 7 , 3 ×   According to the visualization of results, it can be concluded that both the DF-FCCM and SNR-FCCM had no affection on the linear relationship between reconstruction intensity and the concentration. The intensity of reconstruction results from both methods were changed across the concentration of SPIO. To further quantify the linear relationship, we calculated the mean values of regions, and performed the statistical line in Fig. 6. It can be concluded that the linearity between intensity and concentration were not affected by both methods.
Subsequently, we validated that the proposed method did not affect the dynamic range of reconstruction results. The dynamic range scenarios was an important context in MPI technology [29], which means the ability to reconstruct a low intensity target while a high intensity target exist. To quantify the dynamic range of the proposed method, the stripe phantom was used to simulation. The concentration of left target was set to c = 5 × 10 7 [Particles/m 3 ], and the concentration of right target was changed. Five different concentration ratios of two targets : 1:0.25, 1:0.5, 1:1, 1:2, and 1:4 were set. We applied 40 dB Gaussian noise to each signal and the data was compressed to 1526 frequency components by the DF-FCCM and the SNR-FCCM, respectively. The reconstruction results of five concentration ratios were performed in Fig. 7.  According to the visualization of results, the DF-FCCM can reconstruction two targets when the concentration ratio was 1:0.25 and 1:4. Thus, the dynamic range of the DF-FCCM can be seen as 1:0.25 and 1:4. This range was the same to the dynamic range of SNR-FCCM, which mean the DF-FCCM did not affect the dynamic range scenarios of results.
Thus, the DF-FCCM did not affect the quantitative ability of MPI technology, which was an important characteristic. Both the linear relationship and the dynamic range of reconstruction results can be maintained after compression by using the DF-FCCM.

E. Validating the DF-FCCM Can Be Applied in Real Experiments
After validating the DF-FCCM with simulation experiments, the performance of the DF-FCCM in real MPI data was further explored. With the calibration data and the measurement signal from OpenMPI data set, we reconstructed the resolution and shape phantom with the SNR-FCCM and DF-FCCM, respectively. The images of resolution phantom and shape phantom were performed in Figs. 8 and 9, respectively.
Comparing the results of DF-FCCM and SNR-FCCM in 4% threshold, the visualization of results was approximate. However, when the frequency components compressed to 1%, the quality of SNR-FCCM decreased. The SNR-FCCM leads to  the boundaries of the phantoms were blurred and the shape of results were distorted. In contrast, the quality of results decreased minimally when compressing the frequency by using the DF-FCCM. The tubes in the resolution phantom can be easily distinguished and the boundaries of the shape phantom were changed minimally.
To further quantify the quality of results, we applied the reconstruction result with 7630 frequency number as the real result and evaluated other results. The statistical lines were shown in Fig. 10 and the evaluation indicators of specific number frequency components were presented in Table II. According to the tables and the lines, the quality indicators of DF-FCCM results were approximate or better than the SNR-FCCM results. The PSNR value was apparently improved by the DF-FCCM. These results proved that the proposed method can reduce nearly 75% reconstruction time while keep high quality reconstruction.
Thus, the DF-FCCM was shown to well perform when applying to real MPI data. In spite of the real data containing more complex noise, such as the background noise and the distortion  [26], the DF-FCCM can acquire high quality results with 25% reconstruction time. This proved that the DF-FCCM was an advanced signal processing method in MPI technology and can be widely used.

IV. DISCUSSION
In this study, we developed the compression method in MPI by considering the measurement while selecting frequency components. A second selection step with the ESD feature is introduced in the DF-FCCM to select more valuable frequency components, which contain more information and less noise. With this new compression method, high-quality reconstruction can be achieved with less time. Simulation experiments and real experiments validated that the reconstruction can be accelerated nearly 75% with the DF-FCCM.
This improvement is beneficial to the online MPI imaging. Because online imaging can provide feedback for position adjustment [30], achieving online imaging is important to apply technology in biomedical research, such as visual catheter navigation [31], [32]. Previously, the 3D online MPI technology has been accomplished by using the SNR-FCCM [22]. However, because the reconstruction quality has to trade off against the reconstruction time, the update rate of imaging was limited to 1 Hz. The proposed method has potential to be used to improve the online MPI imaging to a 4 Hz or higher update rate, which is essential for the pre-clinical and clinical applications.
In MPI technology, there are three types of noise: the random noise, the background noise and the distortion noise [26]. The background noise and the distortion noise can be removed by using the SNR feature during the first step of DF-FCCM. Thus, in the second selection step, the ESD feature is selected to against the random noise. Based on this two-step approach, we can assume the noise is uniform in (11). Therefore, Gaussian noise is applied in the simulation experiments to evaluate the efficiency of the DF-FCCM.
The compression methods, both the SNR-FCCM and the DF-FCCM, are trade-off between the imaging quality and reconstruction time. When the reconstruction is accelerated, the loss of information and the decrease of quality is unavoidable. In the proposed DF-FCCM, a part of the high SNR but low ESD frequency component was replaced by the high ESD frequency component, which contained more valuable information. Thus, the information loss of DF-FCCM is lower, and it can acquire better results with less time than the SNR-FCCM.
In our experiments, the improvement of DF-FCCM in real experiments is not as clear as in the simulation experiments. The DF-FCCM lead to the center tube separated by other tubes in resolution phantom, which could be an artifact. This phenomenon can be explained by the difference of noise levels. Based on the theory of DF-FCCM, its efficiency is relatively to the random noise level. In the noise simulation experiments, the improvement of DF-FCCM is more apparent when the noise level was 30 dB or 40 dB. However, the real signals from the OpenMPI data set are high quality data. It contains few random noises and might reach 50 dB noise level. In the future, more real experiments with different random noise levels will be used to further validate the efficiency of the DF-FCCM.
The DF-FCCM is a preliminary explore of the frequency component selection method. It can be combined with other compression methods to further accelerate reconstruction. The compressed sensing based methods were proposed to reduce the calibration points dimension because of the highly sparse system matrix in several known transform domains [33], [34], [35], [36]. Moreover, the coded calibration scene method [37], low-rank tensors method [38], and neural network methods [39], [40] are also developed to reduce the calibration point dimension. These methods can be combined with the proposed method to further reduce the data from multi-dimensions.
We introduced the ESD feature in the DF-FCCM to evaluate the measurement signal in this manuscript. More features can be applied to further improve the efficiency of FCCM. The energy of the system matrix is a potential feature that can be incorporated into the FCCM. In previous studies, the energy feature of system matrix was applied as a weight parameter in iteration while the SNR-FCCM was applied to remove frequency components with a low SNR feature [24]. Because the energy of the system matrix is connected to the resolution of result imaging, the energy feature can be added in FCCM to select frequency components which can achieve high-resolution results. This possibility will be explored in future work.
Finally, there are some limitations in this study. Recent research has proved that distortion noise is contained in the MPI signal [26]. Though this noise can be removed by the first selection step by using the SNR feature, it can still affect the performance of the ESD feature when appearing on high SNR frequency components. More research about the influence of the distortion noise will be explored in the future. Furthermore, we used K 1 = 10% as a fixed parameter in our research. An adaptive method for selecting this parameter can be developed to promote the efficiency of the DF-FCCM. Additionally, the proposed method is only validated on the OpenMPI data set. Future work will expand the proposed method to other data sets and MPI instruments to validate the robustness and efficiency of the proposed method.

V. CONCLUSION
The DF-FCCM was developed to accelerate reconstruction in MPI imaging. By adding a second selection step with the ESD feature, the frequency components which contribute more to the reconstruction can be selected by the DF-FCCM. This method can efficiently reduce the reconstruction time while achieving high reconstruction quality. The simulation and real experiments validated that compared to the SNR-FCCM, the DF-FCCM can accelerate more than 75% reconstruction efficiency. The proposed method has shown potential to further accelerate the three-dimensional online imaging, which is beneficial to the pre-clinical and clinical applications of MPI.