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

Abstract

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.

About the Authors:

Huixuan Wu

Roles Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Software, Supervision, Writing – original draft, Writing – review & editing

Affiliation: Department of Aerospace Engineering, School of Engineering, University of Kansas, Lawrence, Kansas, United States of America

Pan Du

Roles Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review & editing

Affiliation: Department of Aerospace and Mechanical Engineering, University of Notre Dame, Notre Dame, Indiana, United States of America

Rohan Kokate

Roles Formal analysis

Affiliation: Department of Aerospace Engineering, School of Engineering, University of Kansas, Lawrence, Kansas, United States of America

ORCID logo https://orcid.org/0000-0003-1322-9747

Jian-Xun Wang

Roles Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Supervision, Writing – original draft, Writing – review & editing

* E-mail: [email protected]

Affiliation: Department of Aerospace and Mechanical Engineering, University of Notre Dame, Notre Dame, Indiana, United States of America

ORCID logo https://orcid.org/0000-0002-9030-1733

Abstract

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.

Figures

Table 2

Fig 1

Fig 2

Fig 3

Fig 4

Fig 5

Fig 6

Table 1

Fig 7

Fig 8

Table 2

Fig 1

Fig 2

Fig 3

Citation: Wu H, Du P, Kokate R, Wang J-X (2021) A semi-analytical solution and AI-based reconstruction algorithms for magnetic particle tracking. PLoS ONE 16(7): e0254051. https://doi.org/10.1371/journal.pone.0254051

Editor: Chi-Hua Chen, Fuzhou University, CHINA

Received: May 6, 2021; Accepted: June 18, 2021; Published: July 9, 2021

Copyright: © 2021 Wu et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: The data that support the fundings of this study are openly available in https://github.com/Jianxun-Wang/AIDenoising.

Funding: HW, National Science Foundation (Grant No. 1944187), https://www.nsf.gov/awardsearch/showAward?AWD_ID=1944187&HistoricalAwards=false JW, National Science Foundation (Grant No. 1934300), https://www.nsf.gov/awardsearch/showAward?AWD_ID=1934300&HistoricalAwards=false The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

List of symbols and abbreviations:

1. 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–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.

A new magnetic particle tracking (MPT) technology has recently been developed to address the issues facing the non-optical methods [12–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 (B1, B2, …) at a few points (or field gradient ∇B), what is the position x and moment m of the magnetic source? The dipole field equation is , 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(B1, B2,…) and m = m(B1, B2,…), 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 ∇B, 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 x(t) = ∫ y(t − τ)K(τ)dτ, 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 signal-to-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–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.

2. Analytical solution

2.1 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 x0 = (0, 0, 0) and x1 = (L, 0, 0). They together provide six signals (B0x, B0y, B0z) and (B1x, B1y, B1z). The magnetic particle moves on the x-y plane (z = 0), but its magnetic moment is 3D, m = (mx, my, mz)′. According to the dipole field model, the magnetic field strengths at x0 and x1 are , where x is the magnet’s location, n is the normal vector in the x–xi 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–xi lies on the x-y plane and ni = (nix, niy, 0)′, the magnetic field equation is(1)where the index i is 0 or 1.

[Figure omitted. See PDF.]

Fig 1. The magnetometer arrangement.

https://doi.org/10.1371/journal.pone.0254051.g001

To calculate m and x, the first step is to eliminate the nonlinear term |x − xi|3. Generally, the z components (Bz and mz) are not zero. Thus, we can normalize the x and y equations using the z component and define M = (Mx, My)’ = (mx /mz, my /mz)’ and Ti = (Tix, Tiy)’ = (Bix/Biz, Biy/Biz)′. Hence, Ti = 3ni(M · ni) − M. The Bz = 0 case will be discussed later. Define a unit vector ti in the x-y plane, ti = (−niy, nix, 0)′, so ti is normal to ni, i.e., ti · ni = 0. Now we decompose the vector Ti to the ti and ni directions:(2.1)as n is a unit vector, and(2.2)

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(3)

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 and the center is at (–Tix/4, –Tiy/4), as illustrated in Fig 2.

[Figure omitted. See PDF.]

Fig 2. Two circles described by Eq 3.

The joints are candidate solutions to M.

https://doi.org/10.1371/journal.pone.0254051.g002

These two circles have two joints, one of which is the solution M. Based on the geometric relationship, we can get two candidate solutions:(4)where and . Here, and E = (T0xT1y − T1xT0y)/4.

Finally, we need to choose the correct solution from the two possible results. To proceed, selecting one solution M, we can obtain tan θ0 = (T0y − My)/(T0x − Mx) and tan θ1 = (T1y − My)/(T1x − Mx). Here θi is the angle between ni and the x-axis (Fig 1). Consequently, the magnet’s position is(5.1)(5.2)

Hence, the position x = (x, y, 0)’, and the magnitude |x − xi|3 can be determined. Thereafter, the z-component equations B0z = −μ0mz/4π|x − x0|3 and B1z = −μ0mz/4π|x − x1|3 provide two possible mz’s with corresponding (mx, my)’ = mz·M. If the selected M is wrong, the two mz’s calculated using B0z and B1z 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 mz = 0 case. The above reconstruction algorithm can be written briefly as functions x = x(B0, B1) and m = m(B0, B1). They are finite for any mz = ± ε, where 0 < ε << 1. As continuous functions that describe a physical trajectory, the left and right limits must be equal, . In other words, the B0z = B1z = 0 case is a removable singularity of the function. In real applications, the probability of measuring an exact zero mz is practically zero, so the above solution works for nearly all measurements.

2.2 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(5)

[Figure omitted. See PDF.]

Fig 3. The 3D setup and coordinate definition.

https://doi.org/10.1371/journal.pone.0254051.g003

Substituting B’ into the 2D solutions (Eqs 3 and 4), we can reconstruct the status of the magnet with the parameter β: x’ = x’(B’) = f(B, β) and m’ = m’(B’) = g(B, β). Then, the position and moment in the original XYZ frame is(6.1)(6.2)

Finally, we need to determine the parameter β by substituting x and m to the dipole equation. Define , where x and m are functions of B and the unknown parameter β, the dipole equation becomes(7)

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. Fig 4 shows a sample curve of ||B–h|| as a function of β.

[Figure omitted. See PDF.]

Fig 4. A sample curve illustrates the relationship between ||B–h|| and β.

https://doi.org/10.1371/journal.pone.0254051.g004

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 pairs. Each pair produces an x and m. The final reconstruction can be obtained by averaging x and m, which reduces errors in the results. Ideally, the uncertainty decreases as at a large N. This is better than the methods in literature, where the uncertainties scales as . However, the reconstruction speed of our method is lower.

3. Trajectory reconstruction

3.1 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 and sudden velocity changes. Fig 5a shows one sample trajectory section in a 3D space. The domain is a cube with the boundary size L = 0.1 m. One magnetometer is located as the origin and the other is at (L, 0, 0). The directions of the sensors are aligned with the frame of reference (XYZ). This setup has been plotted in Fig 3. Given the trajectory, we can simulate the magnetometer readings using the dipole model. Since real measurements always contain uncertainty, we model the field noise using a Gaussian distribution,where is the ground truth obtained from the synthetic trajectory, and ε is a random variable with zero-mean normal distribution. Here we investigate a series of noise levels, beginning with ε = 0 (no noise) to std(ε) = 0.01, 0.03, 0.05, 0.1, … 0.3, where std means the standard deviation.

[Figure omitted. See PDF.]

Fig 5.

a) A sample synthetic trajectory that contains sudden velocity changes and turns. The arrows indicate the magnetic moment direction. b) The original and reconstructed positions at various noise levels. The unit of x is meter. c) The ground truth and reconstructed orientations. Here the orientation is normalized using the magnitude |m|. d) The position errors at various noise levels. The error is normalized using the measurement domain size L. e) The orientation errors at different noise levels. The error is normalized using |m|.

https://doi.org/10.1371/journal.pone.0254051.g005

Fig 5b and 5c show a sample of reconstructed position x and orientation mx. At the zero noise level, our method in Section 2 can reconstruct the position and orientation of a magnet with no error. The zero-noise reconstructed trajectory overlaps the ground truth curve. As the measurement noise increases, the reconstruction error increases. For example, when 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, . 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 , is shown in Fig 5e.

3.2 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, . Written in a matrix format, we have δB = PδV. δB = (δB0x, δB0y, δB0z, δB1x, δB1y, δB1z)′ can be regarded as the measurement noise and δV = (δx, δy, δz, δmx, δmy, δmz)′ is the reconstruction error, where prime means array transpose. The entry in the matrix P is (for j ≤ 3) or (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 thereforeand the orientation error iswhere W1 is the first three rows of W, and W2 is the 4th-6th rows. The maximum ||δx|| can hence be estimated as,(8.1)where σ1 is the square root of the maximum eigenvalue of . A similar argument shows that the max ||δm|| is(8.2)where σ2 is the square root of the maximum eigenvalue of . 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.001Am2 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.

[Figure omitted. See PDF.]

Fig 6. The distributions of a) σ1 and b) σ2, which characterize the reconstruction error in position and orientation, respectively.

The unit of σ1 is [mm/μT] and that of σ2 is [Am2/μT]. Here L = 0.1 m, |m| = 0.001Am2 and pointing at the (1, 1, 1) direction.

https://doi.org/10.1371/journal.pone.0254051.g006

4. 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.

4.1 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.

[Figure omitted. See PDF.]

Table 1. DNN structure.

https://doi.org/10.1371/journal.pone.0254051.t001

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,(9)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 , 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),(10)where Xnoisy and Xclean 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,(11)where represents the forward finite difference of a time series X. Note that if X represents the velocity, is the acceleration. The performance of a denoising model can be characterized by the signal improvement ratio,(12)

Due to significant noise level in velocity, the signal is plotted in the symmetric pseudo log-scale to resize the continuous real value space of velocity axis,(13)where C determines the lowest resolution of the velocity axis.

4.2 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.

[Figure omitted. See PDF.]

Fig 7. Comparison of the performance of a) CNN, b) GRU, c) WT, d) WT+GRU (hybrid model) on Sample 1 randomly drawn from the testing dataset.

https://doi.org/10.1371/journal.pone.0254051.g007

[Figure omitted. See PDF.]

Fig 8. Comparison of the performance of a) CNN b) GRU c) WT d) WT+GRU (hybrid model) on Sample 2 randomly drawn from the testing dataset.

https://doi.org/10.1371/journal.pone.0254051.g008

[Figure omitted. See PDF.]

Table 2. Denoising models comparison.

https://doi.org/10.1371/journal.pone.0254051.t002

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, 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.

5. 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 non-linearity 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.

Acknowledgments

We thank Dr. Weishi Liu in the Mathematics Department at the University of Kansas for our discussion about the analytical solution. We are grateful for Mr. Mohit Nahar Prashanth, who proofreads the manuscript.

Citation: Wu H, Du P, Kokate R, Wang J-X (2021) A semi-analytical solution and AI-based reconstruction algorithms for magnetic particle tracking. PLoS ONE 16(7): e0254051. https://doi.org/10.1371/journal.pone.0254051

References

1. Ni R, Kramel S, Ouellette NT, Voth GA. Measurements of the coupling between the tumbling of rods and the velocity gradient tensor in turbulence. J Fluid Mech. 2015;766: 202–225.

2. Voth GA, Soldati A. Anisotropic particles in turbulence. Annual Review of Fluid Mechanics. 2017;49: 249–276.

3. Lin JS, Chen MM, Chao BT. A novel radioactive particle tracking facility for measurement of solids motion in gas fluidized beds. AIChE Journal. 1985;31: 465–473.

4. Parker DJ, Allen DA, Benton DM, Fowles P, McNeil PA, Tan M, et al. Developments in particle tracking using the Birmingham positron camera. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment. North-Holland; 1997;392: 421–426.

5. Parker DJ, Forster RN, Fowles P, Takhar PS. Positron emission particle tracking using the new Birmingham positron camera. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment. North-Holland; 2002;477: 540–545.

6. Parker DJ, Hawkesworth MR, Broadbent CJ, Fowles P, Fryer TD, McNeil PA. Industrial positron-based imaging: Principles and applications. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment. North-Holland; 1994;348: 583–592.

7. Parker DJ, Broadbent CJ, Fowles P, Hawkesworth MR, McNeil P. Positron emission particle tracking—a technique for studying flow within engineering equipment. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment. North-Holland; 1993;326: 592–607.

8. Wildman RD, Parker DJ. Coexistence of Two Granular Temperatures in Binary Vibrofluidized Beds. Physical Review Letters. 2002;88: 064301. pmid:11863810

9. Heindel TJ. A review of X-Ray flow visualization with applications to multiphase flows. Journal of Fluids Engineering. 2011;133: 074001.

10. Nadeem H, Heindel TJ. Review of noninvasive methods to characterize granular mixing. Powder Technology. Elsevier B.V.; 2018;332: 331–350.

11. Stannarius R. Magnetic Resonance Imaging of Granular Materials. Review of Scientific Instruments. 2017;88: 051806. pmid:28571451

12. Buist KA, van Erdewijk TW, Deen NG, Kuipers JAM. Determination and comparison of rotational velocity in a pseudo 2-D fluidized bed using magnetic particle tracking and discrete particle modeling. AlChE. 2015;61: 3198–3207.

13. Tao X, Wu H. A Comparison of the sequential quadratic programing algorithm and extended kalman filter method in the magnetic particle tracking reconstruction. AIAA Scitech 2019 Forum. 2019. p. 272. https://arc.aiaa.org/doi/abs/10.2514/6.2019-0272

14. Köhler A, Pallarès D, Johnsson F. Magnetic tracking of a fuel particle in a fluid-dynamically down-scaled fluidised bed. Fuel Processing Technology. Elsevier B.V.; 2017;162: 147–156.

15. Köhler A, Rasch A, Pallarès D, Johnsson F. Experimental characterization of axial fuel mixing in fluidized beds by magnetic particle tracking. Powder Technology. Elsevier B.V.; 2017;316: 492–499.

16. Zhang L, Weigler F, Idakiev V, Jiang Z, Mörl L, Mellmann J, et al. Experimental study of the particle motion in flighted rotating drums by means of Magnetic Particle Tracking. Powder Technology. Elsevier B.V.; 2018;339: 817–826.

17. Buist KA, van der Gaag AC, Deen NG, Kuipers JAM. Improved magnetic particle tracking technique in dense gas fluidized beds. AlChE. 2014;60: 3133–3142.

18. Buist KA, Jayaprakash P, Kuipers JAM, Deen NG, Padding JT. Magnetic particle tracking for nonspherical particles in a cylindrical fluidized bed. AIChE Journal. 2017;63: 5335–5342. pmid:29213144

19. Tao X, Wu H. The translational and rotational motions of a cylindrical particle in a granular shear flow inside a split bottom Couette cell. Physics of Fluids. AIP Publishing, LLC; 2020;32: 073310.

20. Tao X, Tu X, Wu H. A new development in magnetic particle tracking technology and its application in a sheared dense granular flow. Review of Scientific Instruments. 2019;90: 065116. pmid:31255032

21. Mema I, Buist KA, Kuipers JAM, Padding JT. Fluidization of spherical versus elongated particles: Experimental investigation using magnetic particle tracking. AIChE Journal. 2019;66: e16895.

22. Wynn W, Frahm C, Carroll P, Clark R, Wellhoner J, Wynn M. Advanced superconducting gradiometer/magnetometer arrays and a novel signal processing technique. IEEE Transactions on Magnetics. 1975;11: 701–707.

23. Nara T, Suzuki S, Ando S. A closed-form formula for magnetic dipole localization by measurement of its magnetic field and spatial gradients. IEEE Transactions on Magnetics. 2006;42: 3291–3293.

24. Wiegert RF. Magnetic STAR technology for real-time localization and classification of unexploded ordnance and buried mines. Detection and Sensing of Mines, Explosive Objects, and Obscured Targets XIV, International Society for Optics and Photonics. 2009. p. 73031U.

25. Yin G, Li P, Wei Z, Liu G, Yang Z, Zhao L. Magnetic dipole localization and magnetic moment estimation method based on normalized source strength. Journal of Magnetism and Magnetic Materials. Elsevier B.V.; 2020;502: 166450.

26. Sui Y, Li G, Wang S, Lin J. Asphericity errors correction of magnetic gradient tensor invariants method for magnetic dipole localization. IEEE Transactions on Magnetics. 2012;48: 4701–4706.

27. Wang C, Zhang X, Qu X, Pan X, Fang G, Chen L. A modified magnetic gradient contraction based method for ferromagnetic target localization. Sensors. 2016;16: 2168. pmid:27999322

28. Jin HH, Zhuang ZH, Wang HB. None-asphericity-error method for magnetic dipole target detection. IEEE Geoscience and Remote Sensing Letters. 2018;15: 1294–1298.

29. Romano GP, Ouellette NT, Xu H, Bodenschatz E, Steinberg V, Meneveau C, et al. Chapter 10 Measurements of turbulent flows. In: Tropea C, Yarin AL, Foss JF, editors. Handbook of Experimental Fluid Mechanics. Springer; 2007. pp. 754–856.

30. Tajane K, Pitale R, Umale J. Comparative analysis of mother wavelet functions with the ecg signals. Journal of Engineering Research and Applications. 2014;4: 38–41.

31. Xiong P, Wang H, Liu M, Zhou S, Hou Z, Liu X. ECG signal enhancement based on improved denoising auto-encoder. Engineering Applications of Artificial Intelligence. 2016;52: 194–202.

32. Vincent P, Larochelle H, Bengio Y, Manzagol PA. Extracting and composing robust features with denoising autoencoders. Proceedings of the 25th international conference on Machine learning. 2008. pp. 1096–1103.

33. Antczak K. Deep recurrent neural networks for ECG signal denoising. arXiv. 2018; 1807.11551.

34. Zhu W, Mousavi SM, Beroza GC. Seismic signal denoising and decomposition using deep neural networks. IEEE Transactions on Geoscience and Remote Sensing. 2019;57: 9476–9488.

35. Arsene CT, Hankins R, Yin H. Deep learning models for denoising ECG signals. 2019 27th European Signal Processing Conference (EUSIPCO). 2019. pp. 1–5.

36. Tian C, Xu Y, Zuo W. Image denoising using deep CNN with batch renormalization. Neural Networks. Neural Networks. 2020;121: 461–473. pmid:31629201

37. Fotiadou E, Konopczyński T, Hesser J, Vullings R. End-to-end trained encoded-decoder convolutional neural network for fetal electrocardiogram signal denoising. Physiological measurement. 2020;41: 015005. pmid:31918422

38. Haque KN, Yousuf MA, Rana R. Image denoising and restoration with CNN-LSTM encoder decoder with direct attention. arXiv. 2018; 1801.05141.

39. Jackson JD. Classical Electrodynamics. 3rd ed. Wiley; 1998.

40. Donoho DL, Johnstone IM. Ideal spatial adaptation by wavelet shrinkage. Biometrica. 1994;81: 425–455.

41. Van der Walt S, Schönberger JL, Nunez-Iglesias J, Boulogne F, Warner JD, Yager N, et al. scikit-image: image processing in Python. PeerJ. 2014;2: e453. pmid:25024921