A semi-analytical solution and AI-based reconstruction algorithms for magnetic particle tracking

Magnetic particle tracking is a recently developed technology that can measure the translation and rotation of a particle in an opaque environment like a turbidity flow and fluidized-bed flow. The trajectory reconstruction usually relies on numerical optimization or filtering, which involve artificial parameters or thresholds. Existing analytical reconstruction algorithms have certain limitations and usually depend on the gradient of the magnetic field, which is not easy to measure accurately in many applications. This paper discusses a new semi-analytical solution and the related reconstruction algorithm. The new method can be used for an arbitrary sensor arrangement. To reduce the measurement uncertainty in practical applications, deep neural network (DNN)-based models are developed to denoise the reconstructed trajectory. Compared to traditional approaches such as wavelet-based filtering, the DNN-based denoisers are more accurate in the position reconstruction. However, they often over-smooth the velocity signal, and a hybrid method that combines the wavelet and DNN model provides a more accurate velocity reconstruction. All the DNN-based and wavelet methods perform well in the orientation reconstruction.


Introduction
Optical-based particle tracking technologies provide crucial knowledge and experimental guidance in the study of turbulence and complex flows [1,2]. However, the advanced optical approaches face severe problems and cannot be used in many opaque environments such as granular motion and dense particulate flows. Missing the experimental guidance, to some extent, retards the development of granular dynamics and multiphase flow theory. A group of non-optical tracking methods is therefore developed, e.g., radioactive particle tracking (RPT) and positron emission particle tracking (PEPT) [3][4][5][6][7][8][9][10], but they encounter new issues. For instance, they require expensive equipment and special expertise for radioactive material operation. The magnetic resonant imaging (MRI) has a limited temporal resolution [11], which is insufficient for high-speed motion measurements.

PLOS ONE
PLOS ONE | https://doi.org/10.1371/journal.pone.0254051 July 9, 2021 1 / 18 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 A new magnetic particle tracking (MPT) technology has recently been developed to address the issues facing the non-optical methods [12][13][14][15][16][17][18][19][20]. The working principle of MPT is to locate a magnetic source according to its field, which is modeled as a dipole. The MPT has the following advantages: 1) it is a safe technique involving no radiation; 2) it provides not only the translation but also the rotation information of a particle, which is critical to granular dynamics; 3) it has a sufficient temporal resolution to measure high-speed flow; 4) it is cost-efficient as the magnetometers and sensors are much less expensive than the equipment in other non-optical approaches.
The key to the MPT method lies in the reconstruction of the position and orientation of a magnetic dipole. In other words, given the measurements (B 1 , B 2 , . . .) at a few points (or field gradient rB), what is the position x and moment m of the magnetic source? The dipole field , where n is the unit vector in the x direction. In experimental applications, numerical optimization and filtering are widely used to calculate x and m, but these methods can be time-consuming and usually rely on artificial parameters or thresholds [12,18,20,21]. In contrast, if the dipole field equation can be inverted, we can find a set of analytical solutions x = x(B 1 , B 2 ,. . .) and m = m(B 1 , B 2 ,. . .), which provides an efficient way to directly calculate x and m without any parameters. In the early stage, a group of researchers developed the eigenvector method [22] and Nara method [23]. Later on, a scalar triangulation and ranging (STAR) method was proposed for real-time magnetic target localization [24], and this method was modified multiple times [25]. These analytical methods possess a clear physical meaning and have been used in practical problems. However, they involve the field gradient tensor rB, which casts special requirements on the magnetometer and the sensor arrangement. In addition, the STAR method has a larger asphericity error [26,27], and its modifications may involve a complex sensor setup [28]. In order to provide an accurate reconstruction method, this paper describes a new analytical solution that can be used for an arbitrarily arranged 3-axis magnetometer array. The reconstruction algorithm is accurate because it contains no assumptions other than the dipole model. However, for reconstruction in practical applications, denoising is an indispensable step since a real measurement contains uncertainty. The classic trajectory denoising method uses linear filtering [19,29]. Given the measured location y, the filtered position is where K is an integration kernel (e.g., a Gaussian kernel). Although this method is simple, the kernel used is fixed and cannot adapt to any temporal (or spatial) variation. To take care of the variation, more sophisticated denoising algorithms based on time-frequency decomposition have been developed. The general process of these approaches is twofold: (1) the original signal is transformed into the frequency domain, where the clean signal and noise can be represented by sparsely distributed spectral coefficients; (2) then a threshold scheme is applied to trim off noisy components, and the remaining are transferred back into the time domain for signal reconstruction. Wavelet transform (WT) is one classic example of time-frequency based methods and has been applied to a wide array of image denoising problems due to its advantage in removing random noise and improving the signalto-noise ratio (SNR), even when noise and signal frequencies overlap [30]. Unfortunately, an appropriate threshold scheme is often hard to determine in the presence of residual or drifting noise [31]. Recent attempts for denoising problems focus on AI-based algorithms, which are capable of suppressing drifting noise and capturing local features robustly given labeled training data. Vincent et al. [32] constructed a denoising auto-encoder (DAE) neural network, aiming to find robust representations of features from noisy input data. Subsequent works are dedicated in optimizing deep neural network (DNN) structures to achieve better performance in handling complex noise and interference [33][34][35][36]. Representative network structures include fully-connected multilayer perceptron (MLP), convolutional neural networks (CNN), and recurrent neural networks (RNN) such as long short-term memory (LSTM) net. The CNN-based methods are typically implemented in an encoding-decoding fashion, where latent features are first extracted by the encoder layers and details are then compensated by the decoder layers to recover a clean version of the original signal [37]. Another popular trend is to utilize RNN to preserve historical information and temporal coherence while denoising, which is effective when handling sequential data, e.g., time series [38]. In this work, we design a novel denoising algorithm by leveraging both unsupervised WT and supervised RNN models with gated recurrent units (GRU), aiming to reduce the noise of the particle trajectory and orientation time series. Using synthetic data, we evaluated the reconstruction performance by comparing it with pure WT, CNN and GRU denoising methods.
The rest of the paper is organized as follows. We first describe the proposed 2D analytical solution and generalize it to 3D in Section 2. Using synthetic data, we evaluate the reconstruction accuracy in Section 3 and then investigate the performance of the proposed denoising scheme in Section 4, where artificial noises are considered. Finally, Section 5 concludes the paper.

The 2D solution
In this section, we first discuss the analytical solution to a restricted problem. In the simple setup shown in Fig 1, there are two 3-axis magnetometers located at Point 0 and 1 with x 0 = (0, 0, 0) and x 1 = (L, 0, 0). They together provide six signals (B 0x , B 0y , B 0z ) and (B 1x , B 1y , B 1z ). The magnetic particle moves on the x-y plane (z = 0), but its magnetic moment is 3D, m = (m x , m y , m z ) 0 . According to the dipole field model, the magnetic field strengths at x 0 and x 1 are ½3n�ðm�nÞÀ m� jxÀ x 0; 1 j 3 , where x is the magnet's location, n is the normal vector in the x-x i direction (i = 0 or 1), m is the magnetic moment, and μ 0 is the magnetic permeability. The dipole field is exact for a uniform spherical magnetic bead, and it is a very good approximation for the far field of non-spherical beads [39]. Since x-x i lies on the x-y plane and n i = (n ix , n iy , 0) 0 , the magnetic field equation is where the index i is 0 or 1. To calculate m and x, the first step is to eliminate the nonlinear term |x − x i | 3 . Generally, the z components (B z and m z ) are not zero. Thus, we can normalize the x and y equations using the z component and define M = (M x , M y )' = (m x /m z , m y /m z )' and T i = (T ix , T iy )' = (B ix / B iz , B iy /B iz ) 0 . Hence, T i = 3n i (M � n i ) − M. The B z = 0 case will be discussed later. Define a unit vector t i in the x-y plane, t i = (−n iy , n ix , 0) 0 , so t i is normal to n i , i.e., t i � n i = 0. Now we decompose the vector T i to the t i and n i directions: as n is a unit vector, and Substitute the components of each vector, we obtain If this equation has a non-zero solution, the coefficient matrix must have a zero determinant. Hence, After rearranging, we get Note that the values of T components are known and fixed in a measurement, and M is unknown. Eq 3 describes two circles for i = 0 and 1. The radius is 3  These two circles have two joints, one of which is the solution M. Based on the geometric relationship, we can get two candidate solutions: S ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi Finally, we need to choose the correct solution from the two possible results. To proceed, selecting one solution M, we can obtain tan Here θ i is the angle between n i and the x-axis (Fig 1). Consequently, the magnet's position is Hence, the position x = (x, y, 0)', and the magnitude |x − x i | 3 can be determined. Thereafter, If the selected M is wrong, the two m z 's calculated using B 0z and B 1z do not match and this solution should be discarded. In addition, if the magnitude of m is known, we can also determine the correct solution by comparing the magnitude of reconstructed m with the known value. Now, we discuss the m z = 0 case. The above reconstruction algorithm can be written briefly as functions x = x(B 0 , B 1 ) and m = m(B 0 , B 1 ). They are finite for any m z = ± ε, where 0 < ε << 1. As continuous functions that describe a physical trajectory, the left and right limits must be equal, In other words, the B 0z = B 1z = 0 case is a removable singularity of the function. In real applications, the probability of measuring an exact zero m z is practically zero, so the above solution works for nearly all measurements.

From 2D to 3D
The above 2D analytical solution can be generalized to 3D with an arbitrary sensor arrangement. As shown in Fig 3, the magnet is located on the x'-y' plane, which has a β angle from the X-Y plane in the 3D space. The x' and X axes overlap. The angle β is regarded as a parameter to be determined later. If we rotate the frame of reference from XYZ to Xy'z', the 3D reconstruction reduces to the 2D problem. After the rotation, the magnetometer readings become

PLOS ONE
Semi-analytical and AI-based MPT reconstruction algorithm  g(B, β). Then, the position and moment in the original XYZ frame is Finally, we need to determine the parameter β by substituting x and m to the dipole equa- , where x and m are functions of B and the unknown parameter β, the dipole equation becomes On the right hand side, the field strength B is known. Recall that the 2D solution contains two possible results, so the ||B-h|| curves may include two branches, but only one of them reaches zero. The function h is highly non-linear, but it contains only one variable, as B are known. Therefore, it is actually efficient to solve Eq 7 using Newton's method in a practical problem, and the reconstructed position in a previous step (a known β) can be used as the initial estimation for the next time step.
For a 3D problem with an arbitrary 3-axis-sensor arrangement (each magnetometer should measure all three components of B), we can apply the above method to any pair of sensors and calculate the averaged reconstruction. In other words, if there are N magnetometers, they form

Trajectory generation and reconstruction
In this section, we will investigate the reconstruction accuracy using synthetic trajectories. The MPT method can be used to study, e.g., fluidization, where collisions among particles can be intense and hard to measure. Therefore, we design the trajectories to contain random turns

PLOS ONE
std(ε) = 0.05, the reconstructed positions deviate from the ground truth significantly. We define the position error as the mean deviation normalized by the measurement domain size, meanðx À x j jÞ=L. Fig 5d shows the position error's dependence on the measurement noise, which is nearly a straight line. For a typical sensor noise level of 3%, the reconstruction error is roughly 3-4% of the domain size, which is consistent with the filtering method [20]. The orientation reconstruction error, defined as meanðm À m j jÞ=jmj, is shown in Fig 5e.

Propagation of uncertainty
The reconstruction error is highly inhomogeneous in the spatial and orientation domain. This is caused by the propagation of error in the reconstruction. To provide a clear explanation, we differentiate Eq 1, dB i ¼ @B i @x j dx j þ @B i @m k dm k . Written in a matrix format, we have δB = PδV. δB = (δB 0x , δB 0y , δB 0z , δB 1x , δB 1y , δB 1z ) 0 can be regarded as the measurement noise and δV = (δx, δy, δz, δm x , δm y , δm z ) 0 is the reconstruction error, where prime means array transpose. The entry in the matrix P is P ij ¼ @B i @x j (for j � 3) or @B i @m j (for j � 4), which is a function of x and m reconstructed using the B. Define the invert of P as W = P -1 , we obtain δV = WδB. The position reconstruction error is therefore dx ¼ W 1 dB and the orientation error is where W 1 is the first three rows of W, and W 2 is the 4 th -6 th rows. The maximum ||δx|| can hence be estimated as, j dx j jj max ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where σ 1 is the square root of the maximum eigenvalue of W 0 1 W 1 . A similar argument shows that the max ||δm|| is j dm j jj max � s 2 j dB j jj ð8:2Þ where σ 2 is the square root of the maximum eigenvalue of W 0 2 W 2 . The coefficient σ 1 and σ 2 characterizes the propagation of error in the reconstruction.
To illustrate the spatial distribution of the error propagator, we select a specific case, in which the magnetic moment is |m| = 0.001Am 2 and it points at the (1, 1, 1) direction. The noise level of B measurements in a L = 0.1 m domain is in the order of micro Tesla (μT). Fig 6 shows σ 1 and σ 2 on a sample plane (z = 0 plane). Note that we present σ 1 and σ 2 with their metric units because the B vary significantly and normalization with one factor can be misleading. The results suggest that there are a few regions with high σ 1 and σ 2 values. Physically, in these regions, the magnetic field strength B are insensitive to at least one component of x or m. As a result, the matrix P becomes nearly singular and the eigenvalues σ 1 and σ 2 become very large. Hence, if the particle trajectory passes through these regions, the uncertainty is larger, which can be manifested as a burst of reconstruction error. In real applications, this error can be removed using more magnetometers in an array.

Denoising using wavelet transform and deep neural network
To examine the performance of WT and AI-based denoisers, we utilize them to process 2000 Lagrangian trajectories that are reconstructed from synthetic signals with 3% random noise. The reconstruction follows the two-magnetometer setting illustrated in Fig 3.

AI-and WT-based denoising method
Two types of DNN structures (i.e., CNN and GRU) are studied for developing the AI-based denoiser. The first proposed denoising method is based on a CNN autoencoder, which contains two 1-D convolutional layers encoding the noisy signal into the latent space and two 1-D deconvolution layers decoding the latent signal to the denoised one after a linear layer. Each convolution/deconvolution layer has a kernel size of 3×1 and stride size of one, followed by a rectified linear unit (ReLU) activation function to account for nonlinearity (more details can be found in Table 1). The convolutional and deconvolutional layers have symmetric parameters such that the output has the same dimension as the input. The noisy trajectory data are fed into the CNN autoencoder with a size of 500×3. The denoised signals can be obtained by the forward evaluation of the network after sufficient training. The second AI-based model uses GRU to preserve the long-term memory of the sequence data. After conducting our preliminary studies on the synthetic data, the GRU is chosen over the LSTM because GRU achieves a  close performance to LSTM with less trainable parameters, hence higher efficiency. The DNN within the denoiser starts with two GRU layers stacked together, which convert three channels into nine channels in the hidden layer, and ends with a linear layer that reduces the dimension back to three channels. The input data structure and implementation of the GRU-based model are identical to that in the aforementioned CNN-based model. More details about the GRU parameters can be found in Table 1. Both DNN-based denoisers are implemented in PyTorch, which is an open-source python platform for machine/deep learning. As mentioned, a synthetic dataset of 2000 trajectories is built for training and testing the DNN-based denoiser. The whole dataset contains 2000 samples, each of which has reconstructed results with artificial noise introduced by the synthetic sensor and the corresponding clean signal as ground truth. The datasets are divided into 1800 samples for training and 200 samples for testing. The training dataset is split into mini batches with a size of 200, which is randomly shuffled, for the neural network to update the model parameters in each epoch. Both CNN and GRU models are trained on an NVIDIA GTX2080Ti GPU using the Adam optimizer with a constant learning rate of 0.001. Each training needs at least 1200 epochs to converge, which takes roughly 0.5 hours on the single GPU. For example, the CNN can achieve a root mean square error (RMSE) of 0.0021 after 1200 epochs of training which takes about 0.58 hours.
In addition, a WT denoising method is also developed for comparison purpose, which is based on the discrete wavelet transformation (DWT) described as follow, where X(k) are the DWT coefficients, x(j) and g(j) represent the input signal and wavelet filters respectively, k is the shift coefficient, and C is the scaling factor (normally chosen as 2). The WT model uses VisuShrink [40], which applies a global threshold defined as g ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 2logQ p , where γ and Q denote the noise variance and number of signal elements (or image pixel), respectively. VisuShrink threshold is renowned for its capability of removing Gaussian noise with high probability and therefore widely applied in image denoising problems. Our WT model adopts Coiflet wavelet basis function with wavelet level set to five and uses soft threshold mode. The WT transformation is implemented by using scikit-image [41] in python and applied on the synthetic data regarding each time series as a 1D image. Note that, different from the AI-based denoisers, which require the clean signals as labeled data to learn from, the WT-based denoiser directly rejects noise in the given signal without any labels. However, as mentioned above, the WT-based approaches have difficulties dealing with residual noises, and it is hard to specify appropriate filter banks and suitable hyperparameters that work for all samples.
We define a metric to describe the relative noise strength, referred to as noise level (N), where X noisy and X clean represent the noisy and clean signals, respectively. Considering that the noisy signal has large fluctuations of high frequency, a second metric is introduced to evaluate the extent of fluctuation in the following form, where _ X represents the forward finite difference of a time series X. Note that if X represents the velocity, _ X is the acceleration. The performance of a denoising model can be characterized by the signal improvement ratio, Due to significant noise level in velocity, the signal is plotted in the symmetric pseudo logscale to resize the continuous real value space of velocity axis, where C determines the lowest resolution of the velocity axis.

Denoising results
We utilize the above methods to denoise the position, velocity, and orientation of the magnetic particle. The denoiser performance is illustrated using two typical samples in the testing dataset (Figs 7 and 8). The ground truth is also plotted for comparison. The signals in both time and frequency domains are studied to evaluate the degrees of noise reduction and oversmoothing. For the position, all the proposed models provide results in good agreement with the ground truth. AI-based denoising methods display better performance than WT models when handling noise with large magnitude. However, large deviations from the ground truth still exist at some time intervals. The noise level and performance metrics are listed in Table 2. The CNN-and GRU-based algorithms outperform the WT model in terms of noise level, fluctuation level, as well as signal improvement ratio. Moreover, the GRU-based model suppresses both the noise and fluctuation better than the CNN-based model, which can be attributed to the inherent advantage of capturing long-term memory of sequential data. The Fast Fourier Transform (FFT) plots show that both the CNN and GRU-based models can significantly reduce the high-frequency noise of the trajectory signals; however, notable noisy components across the entire frequency domain still remain. In comparison, the WT method suppresses much less high-frequency noises. To further enhance the performance, a hybrid method is developed: first the WT filtering is used to preprocess the signal and reduce the noise, and then the GRU is applied on the filtered results to obtain the final velocity. By combining WT and GRU methods, the hybrid model yields the best denoised trajectory signal with a frequency distribution nearly identical to the ground truth. Hence the hybrid model stands out for denoising the position signal. The magnitude of velocity obtained from the noisy signal is 100 times larger than that of the clean signal. This poses a great challenge to velocity denoising, especially for the AI-based models because the magnitude of the true velocity signal is too small to be distinguished from numerical error. Surprisingly both AI-based models have lower noise and fluctuation levels as shown in Table 2. Figs 7 and 8 rectify our observation by showing that the AI-based models fail to capture the shape of the clean velocity signal in both samples. This indicates that the DNN model alone may not be able to reduce the noise without trimming off the velocity signal, which is a limitation of the neural network approach when directly applied on highly noisy data without any preprocessing (e.g., wavelet filtering). On the other hand, the WT method correctly captures the shape of the clean signal but fails to reduce large fluctuations in certain regions (Figs 7c and 8c). According to Figs 7d and 8d and Table 2, the hybrid model successfully attenuates the velocity fluctuations while retaining the shape of the signal. Note that the trajectories studied here are reconstructed using the 2-sensor setup. The burst of error issue as explained in Section 3.2 exists in these trajectories. Using more magnetic sensors can help completely remove the large fluctuation in certain regions and improve the accuracy. In addition, there still exist discrepancies between the shape of the denoised velocity and the clean signal in certain locations (e.g., the left end of velocity plot in Fig 7d), partially due to the fact that the velocity of synthetic trajectories contains Heaviside step functions. The discrepancies can be improved by modifying the preprocessing algorithm and increasing the depth of the DNN, which requires more training data and longer time. Moreover, it can be observed that the GRU-based model significantly oversmoothes the velocity signals, while the performance of the WT method is case-dependent. For sample 1, the WT model well recovers the frequency distribution at low frequency (<60 Hz), but it significantly deviates from the ground truth at higher frequencies (>60 Hz), indicating poor denoising performance. For sample 2, the WT-denoised result shows a good agreement with the ground truth. The hybrid model,

PLOS ONE
though slightly oversmoothes the velocity signals, outperforms the method by GRU or WT model alone.
In terms of the orientation signal, all proposed models suppress the noise well, among which the CNN-based model has the best denoising performance. Interestingly, the denoising performance of the GRU-based model has been surpassed by the WT model. This is because the overall noise magnitude is small and large noise occasionally spikes up in random locations, which might confuse the GRU from a time sequence point of view. It is worth noting that the hybrid model has the best denoising performance for all variables (positions, velocity, and orientation) compared to AI-based models or the WT model, especially for velocity where a large noise magnitude is present.

Summary
This paper presents an analytical solution to the magnet particle positioning problem. Specifically, we can calculate the 3D position and orientation of a magnetic source based on the field vectors B at a few points, which are measured using 3-axis magnetometers. The solution allows us to develop a reconstruction algorithm for 3D particle tracking, which provides a novel approach to detect the flow in an opaque environment such as fluidized bed or turbidity flow. With zero noise, the reconstruction is exact. As the measurement noise increases, the mean reconstruction error grows linearly. For the two-magnetometer setup, the average position error is 3-4%, and orientation error is 5%, given that the measurement noise is 3%. The investigation of uncertainty propagation shows that at the same level of noise, the reconstruction error is highly non-homogeneous in the space and orientation domain, due to the high nonlinearity of the magnetic dipole equation. The burst of error in certain regions can result in severe problems in denoising. Using more sensors can reduce the error.
In order to improve the reconstruction accuracy in practical scenarios where measurement noise is inevitable, we employ the CNN, GRU and WT-based methods to denoise the position, velocity, and orientation signals. All these methods lead to a significant improvement in the reconstruction accuracy. In the regions where the error is extremely large, the AI-based models provide better trajectory reconstruction. However, pure AI-based models fail to capture the shape of the velocity signals because the noise is 100 times larger than the true velocity magnitude. In contrast, the wavelet method can roughly capture the trend of the true velocity signal with a few exceptions. Therefore, we develop a hybrid approach that preprocesses the sequence by filtering out high-frequency noise with WT and then denoises the signal using GRU. The performance of the hybrid method outperforms the other denoisers in velocity denoising. Finally, all the methods perform well in the orientation denoising, and the CNN and hybrid models are slightly better. From the FFT analysis, we observe that the AI-based methods (GRU and CNN) tend to oversmooth the velocity signals, while the hybrid model can better capture the frequency distribution for most cases. In general, the hybrid method combining WT and GRU shows the best performance.