Particle importance analysis: A fast localization method under indeterminate magnetic sources

Abstract

For magnetically actuated robots driven by external magnetic sources (EMSs), rapid and stable localization of the internal permanent magnet (IPM) is essential for closed-loop control. Existing sensor-array-based localization methods require a fixed model with an offline-calibrated dipole moment to separate the magnetic field generated by the IPM from the measured total field. However, model accuracy varies with the distance to the sensor array, leading to inaccurate separation and unstable localization. Additionally, nonlinear optimization is typically used for localization with sensor arrays, which is computationally inefficient. We propose a sampling-based method that jointly estimates the dipole moments and the IPM position, thereby eliminating the need for magnetic field separation. Given the EMS positions and the measured total field, our particle importance analysis efficiently identifies the sampled particle that most closely approximates the dipole position and concurrently estimates the optimal dipole moments. Adaptive particle sampling enables high-accuracy localization using sparsely sampled particles, thereby improving computational efficiency. We validated the approach experimentally with EMSs actuated by manipulators and electromagnetic systems. With a 72-sensor array and the IPM located 20 cm above the sensor plane, the method achieved an average computation time of 1.2 ms and localization errors of 4.70 mm and 4.17°. Compared with an optimization-based method, our method reduced the maximum localization error by 59% and improved computational efficiency by a factor of 30. The method reliably localized in permanent magnet, electromagnetic, and hybrid magnetic fields. Finally, we applied the method to the simultaneous localization and actuation of a magnetic endoscope, validating its feasibility in a colon phantom.

1. Introduction

Magnetically actuated robotic systems have promising applications in medical examination and minimally invasive surgeries, with research demonstrating their feasibility in magnetic endoscopes (Taddese et al., 2018) and endovascular procedures (Von Arx et al., 2023). Such systems typically use external magnetic sources (EMSs) to actuate robots that contain an internal permanent magnet (IPM). The IPM is localized either by embedded sensors integrated into the robot or by external sensor arrays surrounding the workspace. Compared with embedded sensors, external sensor arrays enable more compact robot designs (Xu et al., 2022), improve scalability (Sun et al., 2025), and reduce on-board power consumption. However, external arrays demand computationally intensive processing, offer a limited effective sensing range, and rely on accurate calibration of the EMSs, which hinders their widespread adoption.

External sensor arrays typically measure the combined magnetic field generated by both the IPM and the EMSs, necessitating the separation of the IPM field for localization. Existing approaches often employ a dipole model with offline-calibrated magnetic moments (Son et al., 2018) to estimate the EMS field for this separation, as illustrated in Figure 1(a). However, since the accuracy of the model varies with distance (Petruska and Abbott, 2012), using a fixed dipole moment leads to inaccurate separation. Due to the substantial size difference between the EMS and the IPM, minor separation inaccuracies lead to significant localization errors. When the IPM is far from the sensor array, the separation error combined with disturbances leads to unstable localization performance. Moreover, with the expansion of sensor arrays and the increase in the number of sensors, commonly used nonlinear-optimization-based localization methods become computationally inefficient.

Figure 1.

The illustration of internal permanent magnet (IPM) localization under the effect of external magnetic sources (EMSs). (a) EMSs can be generated by permanent magnets (PM) or electromagnets (EM). Most existing methods require offline segmentation of the IPM magnetic field for localization, but inaccuracies in magnetic dipole moment calibration lead to segmentation errors. (b) Our method avoids the prior knowledge of the EMS moments, using particle sampling and importance analysis to simultaneously calibrate the dipole moments of magnets, and estimate the position and heading of the IPM.

We propose a fast sampling-based localization method that does not require prior knowledge of dipole moments, thereby reducing errors caused by inaccurate field separation. As shown in Figure 1(b), using only the measured magnetic field and the EMS positions, our importance analysis effectively identifies the optimal particle for IPM localization while simultaneously estimating the dipole moments of magnets. The main contributions of this article are as follows:

Importance Analysis: To eliminate reliance on prior knowledge of dipole moments during localization, we treat the dipole moments of both the IPM and EMSs as variables. Sampled particles with unknown moments represent the possible magnet positions in the workspace, and the total magnetic field is modeled as the sum of contributions from all particles. An error reduction ratio approach is then used to identify the most significant particle for IPM localization. Subsequently, for fast dipole moment calibration, we neglect insignificant particles in the model and apply orthogonal least squares to estimate the optimal dipole moments.

Adaptive Particle Sampling: The accuracy of sampling-based methods depends on sampling density and range, but dense sampling over large areas reduces computational efficiency. We propose an adaptive sampling method that begins with sparsely sampled particles across a wide area and dynamically adjusts the sampling range to maintain accuracy while improving efficiency. The method adapts by leveraging overall particle importance from the importance analysis to determine both the convergence direction and rate of the sampling range.

We evaluated the localization performance under the influence of EMSs at multiple distances between the IPM and the sensors. Our method achieves an average accuracy of 4.70 mm and 4.17° when the IPM is 20 cm from the sensor board, with an average computation time of only 1.2 ms. Additionally, we tested the method under permanent magnets, electromagnets, and hybrid magnetic fields. Our method can stably localize the IPM under these conditions while simultaneously estimating the dipole moments of the EMSs. Finally, phantom experiments were conducted to evaluate its performance in simultaneous localization and actuation for a soft endoscope.

The remainder of this article is organized as follows: Section 2 introduces the related work on magnetic localization. Section 3 presents our method, including details on importance analysis and adaptive sampling. In Section 4, we validate the method through simulations and experiments and evaluate its performance. Section 5 concludes the article and discusses future work.

2. Related work

2.1. Embedded and external sensor arrays

Embedded sensors and external sensor arrays are two commonly used sensing configurations in magnetic actuation systems. In the embedded sensor scheme, Hall sensors are placed near the internal permanent magnet (IPM) to measure its relative position with respect to external magnetic sources (EMSs) (Ciuti et al., 2011). To prevent the IPM field from dominating the sensor measurements, the sensor layout must be carefully designed. Some studies place compensating magnets near the sensors to offset the IPM’s magnetic field (Yim and Sitti, 2013), while more commonly, multiple single-axis sensors are deployed at specific locations where the internal magnetic field can be considered negligible (Popek et al., 2016). However, due to the concentrated distribution of sensors, embedded configurations often encounter regions with magnetic field singularities, which significantly degrade localization accuracy. To address this issue, Taddese et al. (2018) superimpose a time-varying magnetic field generated by electromagnets onto the static magnetic field to achieve nonsingular localization. In addition, to obtain global orientation information, Huang et al. (2024) employ Helmholtz coils to generate a uniform directional magnetic field that provides a determined heading angle. In scenarios involving multiple magnetic sources, the above methods do not generalize directly. da Veiga et al. (2023) and Pittiglio et al. (2022) employ a Kalman filter to localize the IPM when two permanent magnets actuate it. In navigation systems with multiple electromagnets, Von Arx et al. (2023) add an oscillating magnetic field onto the actuation field, enabling six-DoF localization using only Hall sensors.

Embedded sensors have been extensively studied and widely applied, yet they still require additional magnetic fields to resolve singularities. Moreover, embedded sensors fail to localize the IPM when the EMS field cannot be measured. External sensor arrays naturally provide global IPM localization regardless of EMS measurability. By collecting measurements from multiple points, they also avoid singularity issues. Sensor arrays further offer high scalability, enabling coverage of the entire human body (Xu et al., 2020) or integration into flexible substrates for in-body capsule tracking (Su et al., 2024b, 2024c). Despite these advantages, external sensor arrays have notable limitations. They use Hall sensors (Son et al., 2015) or magneto-inductive sensors (Watson and Morimoto, 2020) to measure the total magnetic field and rely on the IPM field for localization. The IPM field decays rapidly with distance (Su et al., 2017), making localization sensitive to environmental disturbances. Im et al. (2023) use differential measurements from two-layer arrays to mitigate the disturbances, and Shao and Guo (2019) optimize sensor placement based on signal-to-noise ratios. However, such specialized hardware designs reduce generality. Additionally, computational complexity grows as more sensors are involved (Xu et al., 2022). Optimized placement can reduce bulkiness (Song et al., 2017), but the number of sensors remains far higher than in embedded schemes. Zhang et al. (2023a) mount sensor arrays at the ends of external driving coils, allowing a small number of sensors to localize the IPM in large workspaces. However, this approach is unsuitable for permanent magnets due to sensor saturation. We therefore focus on generalizable localization methods for sensor arrays that achieve fast and stable IPM localization across various types of EMSs.

2.2. Magnetic localization algorithms

Gradient-based methods, including the Levenberg–Marquardt algorithm (LMA), are widely used for magnet localization with sensor arrays, where the magnet position is estimated by minimizing the discrepancy between measured and modeled fields (Wang et al., 2006). These methods are straightforward to construct and implement, but their stability depends on the accuracy of the initial guess (Song et al., 2016a, 2016b). Due to the high nonlinearity of the magnetic field model, a poor initial guess leads to local minima and significant localization errors. Inaccurate initialization also increases the number of iterations, causing variable computation times. Compared with gradient-based methods, heuristic optimization approaches such as particle swarm optimization (PSO) and whale optimization have stronger global search capabilities (Hu et al., 2009). However, when applied to sensor arrays, these methods are computationally expensive due to the large number of sensors and the high-dimensional search space. A common compromise is a hybrid strategy in which heuristics perform a global search for initialization, followed by refinement using LMA (Hu et al., 2010; Lv et al., 2021). While this balances global optimality and computation speed, it remains inefficient for large-scale sensor arrays and is sensitive to model errors and disturbances. Sampling-based methods, such as particle filtering, offer strong global search capabilities and robustness to disturbances (Taddese et al., 2018). However, they also struggle with sensor arrays and are better suited to embedded sensor schemes with fewer sensors. Selecting an optimal subset of sensors can reduce computation and improve efficiency (Xu et al., 2022), but this is ineffective when multiple external magnetic sources are present. Overall, existing methods face computational limitations for large sensor arrays, and errors in the magnetic field model lead to separation errors, resulting in unstable localization.

Magnetic field models based on the point-dipole assumption are widely used due to their computational efficiency, but their accuracy varies with sensor distance. Multipole models can significantly improve modeling accuracy, but they increase parameter complexity and require training on specially designed datasets (Charreyron et al., 2021). Analytical models derived from the Biot–Savart Law provide high precision, but the multiple integrals involved make real-time computation difficult (Zhang et al., 2023b). Using the Jacobian matrix of the magnetic field avoids the inaccuracy introduced by the dipole approximation, but the resulting large-scale matrix reduces computational efficiency when applied to sensor arrays (Di Natali et al., 2016). Learning-based methods have also been proposed to address these challenges. For example, backpropagation neural networks can directly map the relationship between the magnetic field and the magnet position. Fu et al. (2024) and Su et al. (2023) use two neural networks for magnet localization and theoretical model correction, while Su et al. (2024a) propose a self-supervised framework that separates the IPM-generated field under unpredictable interference. However, these approaches typically require offline training, making direct transfer to systems with different configurations difficult.

Table 1 summarizes the localization performance of existing magnetically actuated robots. Most existing methods rely on fixed models obtained through offline calibration. While this is sufficient for embedded sensors that only measure the EMS field, sensor arrays require accurate separation of the IPM field from the total magnetic field. Distance-dependent modeling errors make these methods highly sensitive to sensor height. In addition, larger EMS fields introduce greater magnetic field errors, limiting both EMS size and driving capability. Moreover, limited computational efficiency restricts the number of sensors that can be utilized. We aim to propose a fast localization method that simultaneously estimates the IPM position and calibrates the EMS model parameters. This avoids the instability caused by field separation and allows stable localization at different heights and with different types and sizes of EMSs.

Table 1.

Summary of localization methods in magnetically actuated robots.

Method (sensor config.)	Model	Knowledge of dipole moments	Comp. time (sensor num.)	Accuracy (dist. to sensor)	EMS (size/num.)
Taddese et al. (2018) (embedded)	Dipole model	Predefined/fixed	10 ms (6)	≤5 mm/6°	PM (1048 cm3)
Von Arx et al. (2023) (embedded)	Multipole model	Predefined/fixed	100 ms (1)	≤1 mm/2°	Three EMs
Song et al. (2016b) (external)	Dipole model	Predefined/fixed	142 ms (8)	1.7 mm/2.1° ( $\le 10$ cm)	PM (1.6 cm3)
Xu et al. (2020) (external)	Dipole model	Predefined/fixed	1 s (16)	5.8 mm/7.9° ( $\le 10$ cm)	PM (8.1 cm3)
Xu et al. (2022) (external)	Dipole model	Predefined/fixed	20 ms (8 in 80) a	4.4 mm/6.6° ( $\le 10$ cm) 15.2 mm/78.5° (14 cm)	PM (65 cm3)
Our method (external)	Dipole model	Real-time calibrated	1.2 ms (72)	≤2 mm/3° ( $\le 10$ cm) 4.70 mm/4.17° (20 cm)	PM (672 cm3)/six EMs/hybrid

^aUsing 8 out of the 80 magnetic sensors as a subarray for localization.

3. Method

This section introduces the problem and presents the proposed method, whose algorithmic framework is illustrated in Figure 2. ¹ A summary of the principal symbols used in this section is provided in Table 2.

Figure 2.

The algorithm consists of adaptive sampling and importance analysis. It constructs a state matrix P using the positions of the sampled particles x _i (Section 3.2). The proposed importance analysis evaluates the particle importance values g _c based on the measured magnetic field, and then localizes the IPM position and estimates the dipole moments of both the IPM and EMS (Sections 3.2 and 3.4). The adaptive sampling module adjusts the sampling range according to the standard deviation σ and the maximum particle importance value $g_{\text{c}}^{*}$ (Section 3.3).

Table 2.

List of symbols.

Variables	Descriptions
xc, xa	The dipole positions of the internal permanent magnet and external magnetic sources
mc, ma	The magnetic dipole moment of the internal permanent magnet and external magnetic sources
B, B j	The magnetic field and the measurement of the jth Hall sensor
T, R	Homogeneous transformation matrix and rotation matrix
x i	Position of a particle
r i,j	Vector that points from the ith particle to the jth sensor
P	State matrix for importance analysis
W, A	Orthogonal matrix and triangular matrix decomposed from P
Ω	Auxiliary vector
θ, θ*	Parameter vector related to dipole moments and the vector of calibrated magnetic dipole moment
gc, $g_{\mathrm{c}}^{*}$	Importance vector of all particles and its maximum element
Σ	Standard deviation of the importance values of all particles
ηc, ηs	Thresholds of importance and standard deviation
Ρ	Range matrix of adaptive sampling
Λ	Matrix for adjusting the sampling direction
Q1, Q2	Diagonal weight matrices

3.1. Problem formulation

The considered scenario consists of multiple external magnetic sources (EMSs), an internal permanent magnet (IPM), and an external sensor board, as illustrated in Figure 3. The EMSs generate magnetic fields using either permanent magnets (PMs) or electromagnets (EMs). In PM-based systems, field modulation is achieved through manipulator motion, whereas in EM-based systems it is realized by adjusting the coil current. The magnetization direction of each magnet is aligned with the x-axis of its body frame O _a . The IPM, defined by the body frame O _c , moves in response to the magnetic fields generated by the EMSs. A Hall-sensor-based board is placed beneath the IPM workspace to measure the total magnetic field B produced by both the EMSs and the IPM in real time.

Figure 3.

The coordinate frames of the EMSs, IPM, and sensor board are defined, and the positions and headings of both the IPM and EMS are expressed in the sensor frame O_s. The x-axes of the magnets are defined as their heading directions, which align with their respective magnetic dipole moments.

This article assumes that the magnetic fields generated by the IPM and EMSs follow a point-dipole model, which has a form as

$\mathit{B}\left(\mathit{r},\mathit{m}\right)=\frac{{\mu }_0}{4\pi {‖\mathit{r}‖}^3}\left(3\widehat{\mathit{r}}{\widehat{\mathit{r}}}^{\text{T}}-\mathit{I}\right)\mathit{m}.$

(1)

r is a vector point from a magnet to a sensor, m is the dipole moment of the magnet, $\mathit{I}\in {\mathbb{R}}^{3\times 3}$ is an identity matrix, μ₀ is the permeability of free space, and the superscript hat denotes a unit vector. We also assume that the IPM position can be approximated by the position of the dipole, which is reasonable given the small size of the IPM. At last, the magnetic field measured by the sensor array is within the range of the Hall sensors, which can be achieved by adjusting the distance between the sensors and magnets.

The relationship between the body frames of the IPM, EMS, and the sensor frame can be represented using a homogeneous transformation matrix,

${}^s\mathit{T}_{a,c}=\left[\begin{array}{cc}\hfill {}^s\mathit{R}_{a,c}\hfill & \hfill {}^s\mathit{x}_{a,c}\hfill \\ \hfill \mathbf{0}\hfill & \hfill 1\hfill \end{array}\right]\in {\mathbb{R}}^{4\times 4},$

(2)

where ${}^s\mathit{R}_a,{}^s\mathit{R}_c\in {\mathbb{R}}^{3\times 3}$ are the rotation matrices of EMS and IPM, and ${}^s\mathit{x}_a,{}^s\mathit{x}_c\in {\mathbb{R}}^3$ are their dipole position. The dipole moment of EMS and IPM are ${\mathit{m}}_a,{\mathit{m}}_c\in {\mathbb{R}}^3$ . We define a vector as ${\mathit{e}}_1={\left[\mathrm{1,0,0}\right]}^{\top }$ , the relationship between the rotation matrix and dipole moment is described by (3),

${\widehat{\mathit{m}}}_{a,c}\approx {}^s\mathit{R}^{a,c}{\mathit{e}}_1,$

(3)

The positions of the EMSs, denoted by x _a , are obtained from the robot forward kinematics, while their dipole moments m _a are unknown. The objective is to estimate the IPM position x _c and dipole moment ${\widehat{\mathit{m}}}_c$ in real time within the sensor frame O _s , while simultaneously calibrating the dipole moments of the EMSs.

3.2. Importance analysis

We propose a localization method that rapidly evaluates the contributions of randomly sampled particles to the measured magnetic field and estimates the IPM position based on their relative importance. We first describe the method for the single-magnet case. Random sampling is performed over the feasible space to generate n candidate dipole positions, and importance analysis is then used to evaluate their contributions to the measured magnetic field.

This process is shown in Figure 4(a), the black, red, and blue points are sampled within circular regions with radii of 5 mm, 10 mm, and 15 mm, respectively, from the actual dipole. Their corresponding importance values are shown in Figure 4(b). The importance of particles decreases as their distance from the actual dipole increases, and conversely, it increases as they move closer. Therefore, we treat the particle with the highest importance as an estimated position of the magnet.

Figure 4.

Illustration of importance analysis. (a) Random particles and the estimated dipole position obtained through importance analysis. (b) The importance of particles in different ranges of distance from the actual dipole position. Particles farther from the dipole exhibit lower importance.

The relationship between the particles and the magnetic field is formulated as a matrix equation. The vector from the ith candidate position to the jth sensor is denoted by r _ij , and the corresponding dipole moment is assumed to be m _i . Thereafter, the magnetic field generated by the n particles at a Hall sensor j is

${\mathit{B}}_j=\sum _{i=1}^n\frac{{\mu }_0}{4\pi {‖{\mathit{r}}_{ij}‖}^3}\left(3{\widehat{\mathit{r}}}_{ij}{\widehat{\mathit{r}}}_{ij}^{\top }-\mathit{I}\right){\mathit{m}}_i.$

(4)

The magnetic field ${\mathit{B}}_j\in {\mathbb{R}}^3$ is a three-dimensional vector. To simplify the computation, only the z-component of the magnetic field is used for magnet position estimation. Accordingly, (4) is simplified to the form as

$B_{z,j}=\sum _{i=1}^n\frac{{\mu }_0{\mathit{d}}_{ij}^{\top }}{4\pi {‖{\mathit{r}}_{ij}‖}^3}{\mathit{m}}_i,$

(5)

where ${\mathit{d}}_{ij}=\left(3{\widehat{r}}_{1,ij}{\widehat{r}}_{3,ij},3{\widehat{r}}_{2,ij}{\widehat{r}}_{3,ij},3{\widehat{r}}_{3,ij}^2-1\right)\in {\mathbb{R}}^3$ is a coefficient vector related to the magnetic field in the z-axis.

For a multi-sensor array, the magnetic field model is expressed by

${\mathit{B}}_z=\mathit{P}\mathit{\theta },$

(6)

where $\mathit{\theta }=\left({\mathit{m}}_1,\dots ,{\mathit{m}}_n\right)\in {\mathbb{R}}^{3n}$ is a parameter vector that includes the dipole moment of particles. ${\mathit{B}}_z=\left(B_{z,1},\dots ,B_{z,m}\right)\in {\mathbb{R}}^m$ is the measurement of m sensors in the array. The matrix P is related to the position of potential particles, which has a form as

$\mathit{P}=\frac{{\mu }_0}{4\pi }\left[\begin{array}{ccc}\hfill {\mathit{p}}_{11}\hfill & \hfill \cdots \hfill & \hfill {\mathit{p}}_{n1}\hfill \\ \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill {\mathit{p}}_{1m}\hfill & \hfill \cdots \hfill & \hfill {\mathit{p}}_{nm}\hfill \end{array}\right]\in {\mathbb{R}}^{m\times 3n}.$

(7)

According to the model (5), the elements in P are calculated by ${\mathit{p}}_{ij}={\mathit{d}}_{ij}^{\top }/{‖{\mathit{r}}_{ij}‖}^3$ . The matrix P depends only on the positions of the sampled particles and is therefore known once these positions are specified. In principle, the parameter vector θ and the dipole moments m _i can be estimated using a least-squares method. However, the large number of particles leads to mutual interference among nearby candidates, and the required matrix inversion is computationally expensive. The proposed importance analysis efficiently identifies the most influential particles and improves computational efficiency by discarding those with negligible contributions.

From the construction of P, each column represents the influence of a sampled particle on the sensor measurements. Particles closer to the true dipole position exert greater influence and thus correspond to more important columns. Consequently, the dipole position can be estimated by identifying the columns with the highest importance.

We define P as being decomposed into the product of an orthogonal matrix $\mathit{W}=\left[{\mathit{w}}_1,\dots ,{\mathit{w}}_{3n}\right]\in {\mathbb{R}}^{m\times 3n}$ and an upper triangular matrix $\mathit{A}\in {\mathbb{R}}^{3n\times 3n}$ ,

$\mathit{P}=\mathit{W}\mathit{A}.$

(8)

The importance of a column ${\mathit{w}}_i\in {\mathbb{R}}^m$ with respect to the measurement B _z is evaluated using (9),

$g({\mathit{w}}_i)=\frac{{\left({\mathit{w}}_i^{\top }{\mathit{B}}_z\right)}^2}{\left({\mathit{w}}_i^{\top }{\mathit{w}}_i\right)\left({\mathit{B}}_z^{\top }{\mathit{B}}_z\right)}.$

(9)

This metric is also known as the error reduction ratio (ERR) of an orthogonal vector (Billings and Wei, 2005). Columns with higher ERR values contribute more to the measurement, and the corresponding sampled particles provide estimates of the dipole position.

This process is derived as follows. An auxiliary vector $\mathit{\omega }=\mathit{A}\mathit{\theta }\in {\mathbb{R}}^m$ is defined, and the magnetic model is reformulated as the product of an orthogonal matrix W and the vector ω ,

${\mathit{B}}_z=\mathit{W}\mathit{\omega }.$

(10)

Since W is obtained from the orthogonal decomposition of P, its columns correspond to those of P and are associated with the sampled particles.

Let ω _i denote the elements of the vector ω , which are computed by the following equation:

${\omega }_i=\frac{{\mathit{w}}_i^{\top }{\mathit{B}}_z}{{\mathit{w}}_i^{\top }{\mathit{w}}_i}.$

(11)

By combining (10) and (11), a new equation is obtained and expressed as

${\mathit{B}}_z=\sum _{i=1}^m{\omega }_i{\mathit{w}}_i+\mathit{\xi },$

(12)

where $\mathit{\xi }\in {\mathbb{R}}^m$ is the error vector. By multiplying both sides of the equation by B _z , (13) is obtained,

${\mathit{\xi }}^{\top }\mathit{\xi }-2{\mathit{B}}_z^{\top }\mathit{\xi }=\sum _{i=1}^m{\omega }_i^2{\mathit{\omega }}_i^{\top }{\mathit{w}}_i-{\mathit{B}}_z^{\top }{\mathit{B}}_z.$

(13)

The left-hand side of the equation is related to the error, while the right-hand side shows that the error decreases with each additional column w _i . The corresponding error reduction ratio for column w _i is computed using (14).

$g({\mathit{w}}_i)=\frac{{\omega }_i^2{\mathit{\omega }}_i^{\top }{\mathit{w}}_i}{{\mathit{B}}_z^{\top }{\mathit{B}}_z}.$

(14)

Equation (14) can be further simplified to (9) by substituting (11). It can be observed that ${\sum }_{i=1}^{3n}g({\mathit{w}}_i)\approx 1$ . In our simulations, the cumulative g ( w _i ) of significant vectors exceeds 0.95, indicating that vectors with negligible contributions to error reduction can be discarded.

To improve computational efficiency, importance analysis is integrated into the orthogonal decomposition, enabling dipole position estimation without completing the full decomposition of P in (8).

Let us define the ith column in P as P _i . We use the Gram–Schimit method to decomposite the matrix P, and the orthogonal vector in the matrix W has a form as

${\mathit{w}}_i={\mathit{P}}_i-\sum _{j=1}^{k-1}\frac{{\mathit{P}}_i^{\top }{\mathit{w}}_j}{{\mathit{w}}_j^{\top }{\mathit{w}}_j}{\mathit{w}}_j.$

(15)

Each time a new vector P _i in P is introduced, this formula orthogonalizes the vector with respect to the existing orthogonal vectors w _j . After that, the elements in triangular matrix A are calculated by

$A_{j,i}=\left\{\begin{array}{c}\begin{array}{cc}\hfill 1,\hfill & \hfill j=i\hfill \end{array}\hfill \\ \begin{array}{cc}\hfill \frac{{\mathit{P}}_i^{\top }{\mathit{w}}_j}{{\mathit{w}}_j^{\top }{\mathit{w}}_j},\hfill & \hfill j<i\hfill \end{array}\hfill \\ \begin{array}{cc}\hfill 0,\hfill & \hfill j>i\hfill \end{array}\hfill \end{array}\right..$

(16)

The above method evaluates the importance of the vectors in P. To assess the importance of a given particle, we compute the combined significance of its associated vectors. According to (7), each sampled particle corresponds to three adjacent column vectors in P. Thus, particle k is associated with P_3(k−1)+1, P_3(k−1)+2, and P_3(k−1)+3. Since the orthogonal vectors in (15) are decomposed sequentially, their ERR values can be computed during the decomposition process using (9). For each k ∈ [1, n], we therefore perform orthogonal decomposition on the three corresponding vectors in sequence and sum their ERR values to quantify the importance of particle k. This procedure is applied to all sampled particles, and the particle with the highest importance is selected as the estimated dipole position.

The importance analysis algorithm is presented in Algorithm 1. The inputs include the parameter matrix P, the measurement B _z , the number of dipoles n _d , and the number of particles n. The outputs are the index set of important dipoles Ω_s and their estimated dipole moments ${\mathit{\theta }}^{*}\in {\mathbb{R}}^{3n_d}$ . In our case, a single magnet is considered, that is, n _d = 1. From Line 9 to Line 11, the matrix P is factorized using the Gram–Schmidt method defined in (15). The importance of each sampled particle is evaluated by summing the ERR values of its three associated orthogonal vectors according to (9). In Line 13, the particle with the highest importance is selected as the dipole position estimate. The corresponding orthogonal vectors are retained, while the remaining ones are discarded, as shown in Line 14. Subsequently, the upper triangular matrix A and the auxiliary vector ω are computed using (11) and (16). Given the relationship ω = Aθ, the dipole moment is obtained via backward substitution from ${\theta }_{3n_{\text{d}}}^{*}$ to ${\theta }_1^{*}$ , as shown in Line 19. This avoids matrix inversion and reduces computational complexity. The algorithm estimates the dipole position during decomposition rather than requiring a full decomposition of P, thereby improving computational efficiency. Moreover, by monitoring the norm of the orthogonal vectors ${\tilde{\mathit{w}}}_i$ during decomposition, the method remains effective even when P is rank-deficient. When the norm falls below a predefined threshold, the corresponding particle is neglected, thereby avoiding singularity issues in subsequent orthogonalization using (15).

3.3. Adaptive particle sampling

In importance analysis, the dipole position is selected from sampled particles. Particle density affects localization accuracy, while the sampling range influences global optimality. However, dense sampling significantly increases computational cost. Therefore, instead of uniformly sampling the entire workspace at high density, we dynamically adjust the sampling range to balance accuracy and computational efficiency. The randomly sampled particles x _i are expressed by

${\mathit{x}}_i={\mathit{x}}_0+\mathit{\rho }\left[\begin{array}{c}\hfill u_i\sin \left({\varphi }_i\right)\cos \left({\phi }_i\right)\hfill \\ \hfill u_i\sin \left({\varphi }_i\right)\sin \left({\phi }_i\right)\hfill \\ \hfill u_i\cos \left({\varphi }_i\right)\hfill \end{array}\right],i=1,\dots ,n,$

(17)

where x₀ is the center of the sampling range, ${\varphi }_i\sim U\left(0,\pi \right)$ , ${\phi }_i\sim U\left(\mathrm{0,2}\pi \right)$ , $u_i\sim U\left(0,\pi \right)$ are three randomly sampled variables that describe the particles in the polar coordinate system. $\mathit{\rho }\in {\mathbb{R}}^{3\times 3}$ is a diagonal matrix that adjusts the sampling range. By tuning ρ , the method supports sparse sampling over large regions and dense sampling within smaller regions.

The most straightforward strategy is to begin with a small number of particles sampled over a large range to obtain a coarse estimate. The sampling range is then reduced, and the same number of particles are resampled around the previous estimate, effectively increasing sampling density without increasing computational cost. Importance analysis is repeated to refine the estimate. This process continues until the desired accuracy is achieved. As shown in Figure 5(a), the sampling range progressively contracts while the particle center converges toward the dipole position.

Figure 5.

Illustration of two particle sampling methods. (a) Performing multiple sampling iterations with the shrinking range in one loop. (b) Real-time iteration of the sampling process with varying ranges. (c) Position errors of two sampling methods. (d) The importance of all particles under different sampling ranges. The value of the highest importance and the variance of particle importance can be used as indicators to adaptively adjust the sampling range.

Another strategy adjusts the sampling range in real time. As illustrated in Figure 5(b), the range is continuously updated rather than waiting for convergence within each sampling period, resulting in a gradual narrowing over time. This real-time iteration concept is also employed in optimization-based control. For example, in nonlinear model predictive control using sequential quadratic programming (SQP), fully converged solutions are computationally expensive. Therefore, real-time iteration performs a single QP update per control cycle and uses the result to initialize the next cycle (Minniti et al., 2019), significantly improving computational efficiency.

The real-time iteration sampling consists of two processes: prediction and update. In the prediction process of the kth step, we estimate the current sampling center ${\mathit{x}}_0\left(k\right)$ based on the difference between the estimated position ${\mathit{x}}^{*}\left(k-1\right)$ and the center position of the sampled particles from the previous period. The prediction process is expressed by

${\mathit{x}}_0\left(k\right)={\mathit{x}}^{*}\left(k-1\right)+\mathrm{\Delta }{\mathit{x}}_0\left(k-1\right),$

(18)

where $\mathrm{\Delta }{\mathit{x}}_0\left(k\right)$ is the rate of change of the sampling center, and is calculated by $\mathrm{\Delta }{\mathit{x}}_0\left(k\right)={\mathit{x}}^{*}\left(k\right)-{\mathit{x}}^{*}\left(k-1\right)$ . To prevent sudden changes in the sampling center when the target moves significantly, we limit the maximum value of $\mathrm{\Delta }{\mathit{x}}_0\left(k\right)$ . The prediction provides a coarse estimate, allowing the sampling process to cover the actual position within a smaller range.

The update process adjusts the sampling range by parameter ρ in (17). The results of the importance analysis are used to determine the parameter and assess whether the current error can be further reduced. Recall the error reduction ratio (ERR) defined in (14). For a multi-magnet system with n _d dipoles, let ${\mathit{w}}_{i,1}^{*}$ , ${\mathit{w}}_{i,2}^{*}$ , and ${\mathit{w}}_{i,3}^{*}$ denote the three orthogonal vectors associated with the ith dipole. The sum of their error reduction ratios approximately equals one,

$\sum _{i=1}^{n_d}\sum _{j=1}^{3}g({\mathit{w}}_{i,j}^{*})\approx 1.$

(19)

Let the n _d th dipole denote the IPM, and its importance in each iteration is calculated by

$g_{\mathrm{c}}=\frac{{\sum }_{j=1}^{3}g({\mathit{w}}_{n_d,j}^{*})}{1-{\sum }_{i=1}^{n_d-1}{\sum }_{j=1}^{3}g({\mathit{w}}_{i,j}^{*})}\times 100\%.$

(20)

As the iterations proceed and the sampling range narrows, particle importance is computed using the importance analysis and (20). In each iteration, 50 particles are used for estimation, and their importance values are shown in Figure 5(d). The maximum importance in each iteration is marked with a star, and the corresponding particle is selected to estimate the IPM position. As the sampling range contracts, the maximum importance accounts for an increasing proportion, while the variance among particle importance values gradually decreases.

We define the importance of the estimated dipole acquired by (20) in step k as $g_c^{*}\left(k\right)$ , and the standard deviation of the importance of all the n particles as $\sigma \left(k\right)$ . The update process adjust the sampling range by an expression (21),

$\begin{array}{c}\mathit{\rho }\left(k+1\right)=\hfill \\ \left\{\begin{array}{cc}\mathit{\rho }\left(k\right)+s\mathit{\lambda }\left(k\right)\mathit{\rho }\left(k\right),\hfill & g_{\mathrm{c}}^{*}\left(k\right)⩽{\eta }_{\mathrm{c}}\hfill \\ \left(1-s\right)\mathit{\rho }\left(k\right)+s\mathit{\lambda }\left(k\right)\mathit{\rho }\left(k\right),\hfill & g_{\mathrm{c}}^{*}\left(k\right)>{\eta }_{\mathrm{c}},\sigma \left(k\right)>{\eta }_{\mathrm{s}}\hfill \\ \mathit{\rho }\left(k\right),\hfill & g_{\mathrm{c}}^{*}\left(k\right)>{\eta }_{\mathrm{c}},\sigma \left(k\right)⩽{\eta }_{\mathrm{s}}\hfill \end{array}\right.\hfill \end{array}$

(21)

In this expression, s ∈ [0, 1] is a constant that increases or decreases the range, η _c and η _s are two thresholds. If the importance of the estimate $g_c^{*}$ is below η _c , (21) increases the sampling range to better cover the actual dipole position. When $g_c^{*}$ is above η _c and σ is below η _s , (21) keeps the current sampling range. If both $g_c^{*}$ and σ are higher than the threshold, (21) reduces the sampling range to focus on a more precise area. $\mathit{\lambda }\in {\mathbb{R}}^{3\times 3}$ is a diagonal matrix used to differentiate the sampling range across the three-dimensional directions. We incorporate the deviation of the estimated position x* from the sampling center x₀ to adjust the magnitudes of the adjustments in each direction. The expression of λ expressed by

$\mathit{\lambda }\left(k\right)=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\mathit{x}}^{*}\left(k\right)-{\mathit{x}}_0\left(k\right)\right){\mathit{\rho }}^{-1}\left(k\right).$

(22)

When the sampling range is expanded, if the estimated value deviates significantly from the sampling center in a certain direction, that direction will undergo more aggressive adjustment. Conversely, when the deviation is small, the direction will tend to maintain the original range. When the sampling range is reduced, the trend is reversed.

We compared these two strategies. The first samples over an extensive range in each period and theoretically converge faster. Still, it typically requires three to five iterations to obtain an estimate, resulting in lower computational efficiency. In contrast, real-time iteration converges over multiple periods but is computationally more efficient, and its accuracy is comparable when the magnet motion is continuous and stable. This analysis is confirmed by simulation using the same magnet trajectory. The tracking errors are shown in Figure 5(c). Both methods achieve similar average errors, with real-time iteration achieving an accuracy of 0.2 mm after three iterations.

3.4. Localization under external magnetic sources

To estimate the pose of the IPM under the influence of the EMSs, we need to incorporate the EMS positions into the model formulation.

We define m as the number of sensors, n_c as the number of particles for IPM localization, and n_a as the number of EMSs. Based on the sampled particles and the known EMS position, we can compute two parameter matrices ${\mathit{P}}_{\text{c}}\in {\mathbb{R}}^{m\times 3n_{\text{c}}}$ and ${\mathit{P}}_{\text{a}}\in {\mathbb{R}}^{m\times 3n_{\text{a}}}$ by (7). Their corresponding parameter vectors are ${\mathit{\theta }}_{\text{c}}\in {\mathbb{R}}^{3n_{\text{c}}}$ and ${\mathit{\theta }}_{\text{a}}\in {\mathbb{R}}^{3n_{\text{a}}}$ , respectively. Thereafter, the model (6) is written as

${\mathit{B}}_{\text{z}}=\left[{\mathit{P}}_{\text{c}}{\mathit{P}}_{\text{a}}\right]\left[\begin{array}{c}\hfill {\mathit{\theta }}_{\text{c}}\hfill \\ \hfill {\mathit{\theta }}_{\text{a}}\hfill \end{array}\right]+{\mathit{e}}_0,$

(23)

where e ₀ is the error of the model. In the importance analysis, we extract the subset of the matrix P_c, P_a by ERR, and then we use the least square method to get their moment ${\mathit{\theta }}^{*}\in {\mathbb{R}}^{3\left(1+n_{\text{a}}\right)}$ , as described in Algorithm 1. Let us define the subset as ${\mathit{P}}_{\text{c},s^{*}}\in {\mathbb{R}}^{m\times 3}$ , and s* is the index of the important subset. The importance analysis process is equivalent to solving an optimization problem expressed by

$\begin{array}{c}\underset{{\mathit{\theta }}^{*},s^{*}}{\min }\frac{1}{2}{\mathit{e}}_0^{\top }{\mathit{e}}_0,\hfill \\ \begin{array}{cc}s.t.\hfill & {\mathit{B}}_z=\left[{\mathit{P}}_{\text{c},s^{*}}{\mathit{P}}_{\text{a}}\right]{\mathit{\theta }}^{*}+{\mathit{e}}_0,\hfill \\ \hfill & s^{*}\subseteq \left\{1,\dots ,n_{\text{c}}\right\},\hfill \end{array}\hfill \end{array}$

(24)

Note that directly solving the position using model (23) makes it difficult to ensure continuity. When the magnet is far from the sensors, the measurements become weak and highly sensitive to noise. To mitigate this issue, the previous magnet position is incorporated into the importance analysis. The particle positions are denoted by ${\mathit{x}}_{\text{c},i}\in {\mathbb{R}}^3$ , with corresponding moments ${\mathit{\theta }}_{\text{c},i}\in {\mathbb{R}}^3$ . Their initial guesses, x_c0 and θ _c0, are set to the solution from the previous step. To enhance continuity, the current estimate is constrained to remain close to the previous one. Two constraints are introduced to limit variation between consecutive steps,

$\sum _{i=1}^{n_{\text{c}}}{\mathit{\theta }}_{\text{c},i}\cdot \left({\mathit{x}}_{\text{c},i}-{\mathit{x}}_{\text{c}0}\right)={\mathit{e}}_1,$

(25a)

$\sum _{i=1}^{n_{\text{c}}}{\mathit{\theta }}_{\text{c},i}-{\mathit{\theta }}_{\text{c}0}={\mathit{e}}_2.$

(25b)

In these equations, e ₁ and e ₂ are error vectors that also need to be minimized with different weights. (25b) limits the difference between the sum of the particle moments θ _c,i and the solved moment from the previous step. (25a) makes those particles closer to the previous solution more likely to be assigned larger moments, while points farther away are assigned a moment of zero to minimize e ₁. By combining these constraints with (23), we will obtain a new formula as given in (26a),

$\left[\begin{array}{c}\hfill {\mathit{B}}_{\text{z}}\hfill \\ \hfill \mathbf{0}\hfill \\ \hfill {\mathit{Q}}_2{\mathit{\theta }}_{\text{c}0}\hfill \end{array}\right]=\left[\begin{array}{cc}\hfill {\mathit{P}}_{\text{c}}\hfill & \hfill {\mathit{P}}_{\text{a}}\hfill \\ \hfill {\mathit{Q}}_1\mathit{C}\hfill & \hfill \mathbf{0}\hfill \\ \hfill {\mathit{Q}}_2\mathit{D}\hfill & \hfill \mathbf{0}\hfill \end{array}\right]\left[\begin{array}{c}\hfill {\mathit{\theta }}_{\text{c}}\hfill \\ \hfill {\mathit{\theta }}_{\text{a}}\hfill \end{array}\right]+\mathit{E},$

(26a)

where C and D are two matrices that constrain the changes in solved position and moment, as represented by

$\mathit{C}=\left[\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\mathit{x}}_{\text{c},i}-{\mathit{x}}_{\text{c}0}\right),\dots ,\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\mathit{x}}_{\text{c},n_{\text{c}}}-{\mathit{x}}_{\text{c}0}\right)\right]\in {\mathbb{R}}^{3\times 3n_{\text{c}}},$

(26b)

$\mathit{D}=\left[\mathit{I},\dots ,\mathit{I}\right]\in {\mathbb{R}}^{3\times 3n_{\text{c}}}.$

(26c)

In these equations, ${\mathit{Q}}_1,{\mathit{Q}}_2\in {\mathbb{R}}^{3\times 3}$ are diagonal weight matrices, and E is an error vector that is related to e ₀, e ₁, and e ₂. Based on the equality constraint (26a), we can construct a similar optimization problem as (24), and use the same importance analysis method to solve the estimated positions.

4. Experiments

This section evaluates the feasibility and performance of PIA through five sets of simulations and experiments. Simulations verify the robustness of PIA to model parameter errors and sensor noise. Experiments demonstrate that, compared with conventional optimization-based methods, PIA achieves faster computation, a longer sensing range, and accurate estimation of the IPM position and heading under different EMS configurations. Finally, PIA is applied to simultaneous localization and actuation in magnetic endoscopy using a colon phantom.

The experimental platform consists of manipulators, electromagnetic systems, a sensor board, and an NDI camera. The EMSs include permanent magnets, electromagnets, and hybrid magnetic sources. Permanent magnets are actuated by manipulators, while electromagnets are controlled by regulating coil currents. A computer running Ubuntu 24.04 on an Intel i7-11850 processor with 32 GB RAM handles manipulator control and magnet localization. The sensor board contains 72 Hall sensors, which are processed by an STM32 microcontroller and transmitted to the computer via serial communication. Both the manipulators and the NDI camera communicate with the computer through EtherCAT. ROS2 is used to publish and subscribe to sensor data, camera measurements, and manipulator joint states.

The hardware and algorithm parameters used in the simulations and experiments are summarized in Table 3. Three permanent magnets of different sizes were tested, with nominal dipole moment magnitudes of 4.16, 16.6, and 780 Am². These values are used only for analysis, since PIA does not require prior knowledge of the magnetic moment. Among the algorithm parameters, Q ₁ and Q₂ weight the deviation of the current estimates from the previous position and heading. Smaller values reduce magnetic field errors. The parameters η _c and η _s control adaptive sampling and correspond to the percentage error and the standard deviation of particle importance, respectively.

Table 3.

Experimental parameters.

​Parameters
Hardware Parameters:
EMS (∅30 × 20 mm); dipole moment: 16.6 Am2
EMS (80 × 80 × 105 mm); dipole moment: 780 Am2
IPM (∅15 × 25 mm); dipole moment: 5.20 Am2
Number of sensors: 72; range:±50 Gs
Sensor board size: 30 × 36 cm
Algorithm Parameters:
Weight for IPM: Q1 = 10−5I, Q2 = 10−5I
Number of particles: 50
Sampling parameters: η s = 0.01, η c = 90, s = 0.05

4.1. Localization under disturbances

Simulations show that the proposed method maintains higher accuracy and better continuity under parameter uncertainties and measurement noise compared with optimization-based methods and particle filtering.

In the simulation, the positions and headings of the IPM and EMSs are predetermined, and the theoretical measurements of the Hall sensors are calculated using a point-dipole model. Our method inversely solves for the IPM position and heading using the generated measurements as inputs. Both the optimization-based method and the particle filter require compensation for the magnetic field generated by the EMSs before estimating the IPM state from the residual field. The optimization-based method estimates the IPM position and moment by minimizing the difference between measured and theoretical magnetic fields. However, when the initial guess is far from the actual state, the optimization may converge to an inaccurate solution or fail to converge within the maximum number of iterations. The particle filter we use is similar to that in Taddese et al. (2018). Both use a random walk model as the process model, but the measurement model is different. Taddese et al. (2018) use an IMU to measure the roll and pitch, while a particle filter estimates its position and yaw. Since we do not use an IMU, we assume a known heading in simulations and apply the particle filter only for position estimation. This design is intended to compare the differences between the particle filter and our method in position estimation. We evaluate these methods under two types of disturbances. The first, parameter disturbance, simulates errors in EMS calibration by using inaccurate moments randomly sampled between 0.9 and 1 times the actual moment. The second, measurement disturbance, simulates sensor noise and external interference by generating measurements randomly sampled between 0.9 and 1.1 times the ideal values. All three methods estimate the IPM state based on these perturbed measurements.

In the first set of simulations, no disturbances are added to the measurements or parameters. The ideal trajectory in the xy-plane is a Lemniscate of Gerono, and the IPM trajectory is expressed by the equation: $x=0.11cos\left(\theta \right)$ , $y=0.11sin\left(\theta \right)$ , z = 0.16 − 0.05θ/360. The trajectory of the EMS is offset from the IPM trajectory by $\left[0.03,0.01,0.13\right]$ m, positioned above the IPM. The heading of the IPM is tangent to the trajectory in the direction of movement, while the heading of the EMS is opposite to that of the IPM. Figure 6(a) shows the results of three methods, and their mean absolute errors (MAE) of positions are given in Figure 6(b). The errors of our method, optimization-based method, and particle filter are 0.079, 0.001, and 0.808 mm, respectively. The maximum errors are 0.5, 0.1, and 1.6 mm, respectively. The optimization-based method achieves the highest accuracy without disturbances, followed by our method. The accuracy of the particle filter depends on the density of the sampled particles. Higher sampling density improves accuracy but will reduce computational efficiency.

Figure 6.

(a) and (b) are the localization results without disturbances in model parameters and measurements. (c)–(f) show the results when the EMS moment includes a 10% random noise. (a) The 3D localization performance of three approaches—our method, an optimization-based method, and particle filter. (b) The mean absolute localization error. When the model parameters and measurements are perfect, three methods have similar accuracy. (c) The localization performance of the three methods under inaccurate moment. (d) The mean absolute position error when the dipole moment is not accurate. (e) The heading error for the optimization-based method and our method. (f) The EMS moments used in the three algorithms are randomly sampled from 0.9 to 1.0 times the actual EMS moment. Our method does not rely on an accurate EMS moment and can rectify it in real time.

The second simulation adds disturbances to the EMS moment, and the trajectories of IPM and EMS are the same as in previous simulations. The results of localization are shown in Figure 6(c). The errors of position and heading are given in Figures 6(d) and (e). The mean position errors for the three methods are 0.089, 6.785, and 4.234 mm, respectively. The maximum errors reach 0.486, 17.21, and 7.427 mm, respectively. The mean heading errors for our and optimization-based methods are 0.014 and 2.893 degrees. The inaccurate moment is shown in Figure 6(f) by a blue line, and the red line is the rectified dipole moment through our method. When parameter disturbances exist, the accuracy of the optimization-based method decreases significantly. Our method achieves the highest accuracy, followed by the particle filter. Since our method does not require compensation for the magnetic field generated by the EMS, it does not rely on an accurate EMS calibration. Moreover, even though we use the inaccurate EMS moment for initialization, our method can rectify the EMS in real time, maintaining high IPM localization accuracy.

In the third simulation, disturbances are added to the measurements. The trajectories of EMS and IPM are Archimedean spiral. The IPM trajectory is expressed by: $x=\left(0.05+0.06\theta /360\right)cos\left(\theta \right)$ , $y=\left(0.05+0.06\theta /360\right)sin\left(\theta \right)$ , z = 0.16 − 0.05θ/360. The trajectory of the EMS is still offset from the IPM trajectory by $\left[0.03,0.01,0.13\right]$ m. The results are shown in Figure 7(a). The position errors for our method, optimization, and particle filter are 0.091, 6.848, and 2.825 mm, while the maximum errors are 0.512, 18.14, and 5.294 mm, respectively. The heading errors for our method and optimization are 0.010 and 1.963 degrees. Figure 7(d) shows the mean measurements of sensors after the disturbances are added. When disturbances are added to the measurements, our method still achieves the highest accuracy, followed by the particle filter. The accuracy of the optimization-based method decreases significantly under the influence of disturbances. Our method demonstrates more stable localization performance under disturbances and produces a more continuous trajectory.

Figure 7.

Localization results when measurements are inaccurate. (a) The localization performance of the three methods on an Archimedean spiral. (b) The mean absolute localization error. (c) The heading error for the optimization-based method and our method. (d) The average measurements from the sensor array and the disturbed measurements. The disturbances are generated by randomly sampling within ±10% of the current measurements. Our algorithm can achieve good localization accuracy and stability under disturbances.

4.2. Test of computational efficiency

This section tests the performance of our method in tracking fast-moving target and compares it with the optimization-based method.

In this experiment, a cylindrical magnet with a dipole moment of 16.6 Am² is mounted at the end of the manipulator. Its motion trajectory followed an ellipse with radii of 25 cm and 15 cm, with a height of 18 cm. The manipulator executed this trajectory at five different speeds, with motion periods of 1.5 s, 2 s, 3 s, 5 s, and 10 s, respectively. Our method estimates the magnet position based on the magnetic field measured by the sensor board.

As shown in Figure 8(a), our method maintains an approximate deviation from the desired path when tracking targets at different speeds. When the speed increases to one cycle per 1.5 s, a noticeable deviation appears, but the error quickly converges. We compared our method with the optimization-based method and summarized the position and heading errors at different speeds in Figure 8(b). As the motion speed increases from 10 s per cycle to 1.5 s per cycle, the position MAE of our method increases from 5.30 mm to 17.0 mm, while the optimization-based method increases from 4.8 mm to 25.8 mm. The angular MAE of the two methods increases from 0.87 degrees and 2.18 degrees to 2.67 degrees and 2.69 degrees, respectively. In this process, the maximum error of our method is 34.6 mm and 5.17 degrees, while the optimization-based method has a maximum error of 91.5 mm and 4.33 degrees. It can be observed that as the speed increases, the positioning accuracy of our method degrades more slowly. At higher target speeds, our method therefore achieves superior localization accuracy. This performance difference is primarily due to computational efficiency. As shown in Figure 8(c), when the motion period is 1.5 s, the optimization-based method exhibits noticeable tracking lag in the x and y directions. Solving the optimization problem requires an average of 37.5 ms, whereas our method requires only 1.31 ms on average. When tracking a slower-moving target, there is no significant difference between the two methods. However, as the speed increases, the higher computation time leads to more significant tracking errors. Meanwhile, the localization results of our method are also more densely distributed, allowing it to perceive sudden rapid movements more effectively. For example, in the scenario of magnetically actuated endoscopy, when the endoscope gets stuck in the colon, it may experience rapid slipping under the influence of the magnetic field, moving several centimeters within tens of milliseconds. Our method can quickly localize this process in real time.

Figure 8.

Tracking trajectory at different moving speeds of the magnet. (a) Desired trajectory and tracking results. (b) Position and heading errors of our method versus the optimization-based method. (c) Tracking results in the x and y directions. The optimization-based method shows a larger phase difference due to its higher time cost.

4.3. Localization at different sensor distances

We verified the stability of the method at different sensor–IPM distances and demonstrated its effectiveness in IPM localization under out-of-plane motion.

In the first experiment, two manipulators moved two magnets above the sensor board along a Lemniscate of Gerono trajectory. One magnet served as the EMS (16.6 Am²), and the other as the IPM (5.2 Am²) to be localized. We assessed the IPM localization performance at four heights: 10 cm, 12 cm, 16 cm, and 20 cm above the sensor board. We compared our method with the optimization-based approach, and the localization results are shown in row i of Figure 9. As the distance between the IPM and the sensors increases, the accuracy of both methods gradually declines. However, our method maintains higher accuracy and produces smoother localization results. The MAEs of both methods are shown in row ii. When the height is below 12 cm, both methods achieve similar accuracy, with MAEs within 2 mm for position and 3 degrees for heading. The maximum errors of our methods stay below 4 mm and 6 degrees, while the optimization-based method remains under 5 mm and 9 degrees. At 16 cm, both methods have a position MAE of 2.4 mm. However, the maximum errors for our method are 6.2 mm and 5.83 degrees, whereas those for the optimization-based method are 9.7 mm and 10.65 degrees. At 20 cm, the position MAEs are 4.70 mm for our method and 5.21 mm for the optimization-based method. The heading MAEs are 4.17 degrees and 4.72 degrees, respectively. Our maximum errors are 11.3 mm and 10.32 degrees, while those of the optimization-based method reach 27.3 mm and 28.7 degrees. Overall, our method maintains reliable localization performance up to an IPM height of 20 cm, with EMS heading estimation errors remaining within 6.3 degrees.

Figure 9.

(a)–(d) show the localization results of the IPM at different distances from the sensor board under the influence of an indeterminate EMS. Row i presents the actual trajectories of the EMS and IPM, together with the localization results from PIA and an optimization-based method. Row ii shows the mean absolute position and heading errors of both methods. Row iii displays the sensor reading ranges and magnetic field errors. (e) illustrates the result for out-of-plane motion, where the error increases with height and decreases effectively as the height is reduced.

Row iii of Figure 9 shows the estimated dipole moments and the magnetic field errors. The red area indicates the range of the magnetic field measurements. The red and black lines represent the percentage magnetic field errors for our method and the optimization-based method, respectively, defined as the mean absolute magnetic field error relative to the maximum sensor reading. As the IPM height increases, the measured magnetic field decreases from a maximum of 30 Gs to below 4 Gs, while the magnetic field error gradually increases. At a height of 20 cm, the percentage error reaches a maximum of 12.6% for our method, which is 2.1 times the maximum percentage error at 10 cm. The increase in magnetic field error leads to inaccurate localization. However, compared with the optimization-based method, our approach yields smaller magnetic field errors. As the magnet height increases, the dipole moments we estimate in real time change from 17.22 and 5.57 Am² to 16.54 and 5.32 Am², gradually approaching their nominal values. Treating the dipole moment as a variable parameter in the model allows it to adapt to real-time magnetic field changes, which partially compensates for model errors. The optimization-based method uses fixed dipole moments to separate the magnetic fields, and the resulting separation errors reduce the localization accuracy of the IPM. In contrast, our method estimates the dipole moments of both the EMS and IPM in real time during localization, avoiding magnetic field separation. This leads to smaller magnetic field errors and, therefore, better localization performance.

In the second experiment, the manipulator guided the IPM along a helical trajectory relative to the sensor array, with distances ranging from 8 cm at the closest point to 18 cm at the farthest. Out-of-plane motion tests demonstrate that the PIA stably localizes the IPM while it moves relative to the sensor array. In Figures 9(e–i), the IPM performs a helical motion perpendicular to the sensor plane, with minimum and maximum distances of 8 cm and 18 cm from the sensors, respectively. Figures 9(e–ii) shows that the maximum position and heading errors are 8.55 mm and 5.34 degrees, while the minimum errors at lower heights are 0.21 mm and 0.91 degrees. Figures 9(e–iii) shows that when the IPM is within 14 cm of the sensor array, the position error remains below 2 mm, and maximum errors increase with height. In Figures 9(e–iv), the percentage error represents the proportion of magnetic field error relative to the maximum measured field. As the IPM moves farther from the sensor array, field fluctuations caused by random noise increase, resulting in reduced localization accuracy. Although the error grows when the IPM moves away from the sensor array, our method enables rapid convergence and error reduction as the IPM approaches the sensors.

4.4. Localization under different indeterminate EMSs

This section tests our method under indeterminate EMSs through three experiments using permanent magnets, rapidly varying electromagnets, and hybrid magnetic sources. The results demonstrate that our method effectively localizes the IPM under different magnetic sources.

In the first experiment, as shown in Figure 10(a), the manipulator controls the motion of two permanent magnets, while the IPM (5.2 Am²) is manually moved above the sensor board. The hardware setup is similar to the previous experiment as shown in Figure 10(b). We estimate the position of the IPM using both our method and the optimization-based method, with the results shown in Figure 10(c). Both methods can localize the IPM in real time. However, in the optimization-based method, the estimated position occasionally deviates from the ground truth, shifting closer to the EMS position. The errors corresponding to the two methods are shown in Figures 10(d) and (e). The optimization-based method yields a position MAE of 4.3 mm and a heading MAE of 6.28 degrees, while our method achieves a position MAE of 3.2 mm and a heading MAE of 4.21 degrees.

Figure 10.

The localization results of IPM under the influence of two permanent magnets. (a) The experimental platform involves two manipulators controlling the motion of two EMSs. A sensor board is installed at the bottom to estimate the position of the IPM. Markers are attached to these magnets to measure their actual positions. (b) The hardware configuration. (c) The motion trajectories and estimated trajectories of the three magnets. (d) The position MAE of PIA and optimization-based method. (e) The heading MAE of our and optimization-based method. (f) The magnetic fields generated by the three magnets.

The blue regions indicate periods when the estimated values exhibit significant deviations. In these regions, the optimization-based method reaches maximum errors of 15.4 mm and 32.8 degrees. The maximum errors of our method are only 9.54 mm and 12.9 degrees. The theoretical magnetic fields generated by the three magnets are shown in Figure 10(f). As the error increases, the magnetic field generated by the IPM accounts for a smaller proportion of the total magnetic field. Combined with the disturbances caused by inaccurate dipole moment, this results in significant deviations. By considering historical estimates of the IPM position and real-time calibration of the dipole moment, our method maintains good stability despite interference from multiple EMSs.

In the second experiment, we evaluated the localization performance of PIA under rapidly varying electromagnetic fields, which the Navion generates. Navion is a robotic navigation system equipped with three coils, capable of generating magnetic fields in arbitrary directions over a large workspace (Boehler et al., 2022; Gervasoni et al., 2024). As shown in Figure 11(a), the Navion generates a magnetic field with a magnitude of 10 Gs, rotating around the x-axis at a frequency of 1 Hz. The sensor board is placed vertically approximately 0.5 m in front of the Navion. By applying a set of known currents to Navion, we can determine the positions of three coils in the sensor frame based on the measurements. The IPM, mounted on the end of a KUKA manipulator, followed a three-turn helical trajectory (6 cm radius) away from the sensor board, moving from 11 cm to 17 cm. An NDI camera measures the actual IPM position. The hardware setup between the sensors, manipulator, and computer is the same as in the previous experiment. A separate computer with Ubuntu system controls Navion and communicates with the host via TCP. Navion provides current feedback to the host in real time. In this experiment, the currents in the three coils vary sinusoidally, with a peak amplitude of about 1.8 A.

Figure 11.

The localization results of IPM under electromagnetic field. (a) The experimental platform, hardware setup and measured currents. The magnetic field is generated using the Navion, and the IPM is controlled by a manipulator to move along a helical trajectory. An NDI camera records the actual trajectory. (b) The actual IPM trajectory and the results estimated by PIA and optimization-based method. (c) The position and heading error of two methods. The error associated with the optimization-based method increases as the IPM moves away from the sensor, whereas our method demonstrates small error variation. (d) The average magnetic field on sensors generated by Navion and the average total magnetic field. (e) The magnetic field error obtained using our method and the norm of calibrated IPM moment.

The localization results of the PIA and optimization-based method are given in Figure 11(b), with the latter exhibiting significantly larger errors. Figure 11(c) shows the position and heading errors for three turns. As the number of turns increases, the IPM moves further from the sensor board. For the optimization-based method, the mean absolute error was 5.01 ± 11.70 mm in the first turn, increased to 6.66 ± 25.47 mm in the second, and reached 9.89 ± 29.73 mm in the final turn. In contrast, the mean absolute errors for PIA were 3.76 ± 7.31 mm, 4.55 ± 10.34 mm, and 4.24 ± 7.28 mm for the respective turns. The heading errors for two methods are 5.90 ± 33.83 and 1.95 ± 5.01 degrees, respectively. As the IPM moves away from the sensor, the optimization-based method shows a 97% increase in average error and an increase of 1.7 times in maximum error. In contrast, PIA exhibits minimal variation, maintaining an average error within 5 mm. Figure 11(d) shows the average magnetic field generated by the coils and the total field measured by the sensor. As the IPM moves farther away, the coil-generated field increasingly dominates. Inaccurate model parameters and current feedback delays hinder the optimization-based method from accurately isolating the magnetic field of the IPM, causing localization errors. In contrast, our method maintains stable performance despite high-frequency magnetic field variations. Figure 11(e) shows the real-time calibrated IPM moment obtained by PIA, with its magnitude deviating by up to 13.6% from the nominal moment. The red line represents the mean absolute error between the magnetic field computed from the calibrated moment using the dipole model and the measured field, with a maximum error of 0.19 Gs.

In the third experiment, we localize the position and moment of the IPM in a hybrid magnetic field generated by both permanent magnet and electromagnets. Figure 12(a) shows the configuration of the experimental platform. Six electromagnets are installed 0.14 m beneath the sensor board, with a 0.1 A current applied alternately to each for 0.5 s. The permanent magnet, initially positioned 0.45 m above the sensor board, performs six-DoF sinusoidal motion in accordance with the robot movement as shown in Figure 12(b). The real-time position of the permanent magnet is determined based on the manipulator kinematics, while the positions of the six electromagnets are obtained through calibration and treated as fixed values during the localization process. The IPM, equipped with markers, is manually moved above the sensor board, and the ground truth is recorded using an NDI camera. We compared our method with the optimization-based method, and the estimated trajectories are shown in Figure 12(c). The optimization-based method exhibits significant deviation from the actual trajectory in the hybrid magnetic field, whereas our method provides a stable estimation of the IPM position. Figure 12(d) presents the mean absolute errors of our method, with position and heading errors of 4.01 mm and 6.29 degrees, respectively, and maximum errors of 26.91 mm and 35.50 degrees. The optimization-based method is unstable because the fixed-moment point-dipole model struggles to segment the magnetic field of the IPM from the hybrid magnetic field. As shown in Figures 12(e) and (f), the EMS moment estimated by our method fluctuates around its nominal value, with a maximum deviation of up to 16%. When identical currents are applied, the moments of the electromagnets exhibit variability, resulting in magnetic field errors. The cumulative effect of these errors across multiple electromagnets causes the estimated position to deviate from the actual position. Our method demonstrates strong stability under hybrid magnetic fields, attributed to real-time moment calibration during the localization process.

Figure 12.

The localization results of IPM under hybrid magnetic field generated by a permanent magnet and six electromagnets. (a) The hybrid magnetic system consists of one permanent magnet and six electromagnets. The IPM moves above the sensor board, while the ground truth is collected using an NDI camera. (b) The motion of the manipulator end and the desired currents in electromagnets. Six electromagnets are supplied with a 0.1 A current in an alternating manner. (c) The motion trajectories of magnets and the localization results of the PIA and optimization-based method. (d) The position and heading error of our method. (e) The nominal and calibrated moment of permanent magnet. (f) Calibrated moments of six electromagnets.

4.5. Test in colon phantom

Our method enables simultaneous localization and actuation of a magnetic endoscope in a colon phantom. The experimental setup is illustrated in Figure 13(a). A magnetic endoscope equipped with an IPM and an onboard camera was actuated by an EMS mounted on the KUKA LBR Med 14. The endoscope was deployed within a colon phantom placed in a human abdominal model to simulate realistic conditions. The sensor array, consisting of 72 Hall sensors, was positioned beneath the abdominal model to measure magnetic fields across the entire abdominal region. The ideal IPM trajectory within the phantom is shown in Figure 13(b). During the experiment, the distance between the IPM and the sensor array varied from 17 cm to 20 cm. The endoscope heading followed the trajectory tangent, and the trajectory was divided into five segments, with lower velocities applied in Segments II and IV to enhance stability during turns. The manipulator controller runs on an Ubuntu system, with the localization and control algorithms implemented as two ROS2 nodes. The localization node receives magnetic field measurements and manipulator joint positions. It outputs the estimated position and magnetic moment of the IPM, as well as the EMS moment, at 100 Hz. A model predictive controller (MPC), detailed in Appendix B, is employed to control both the manipulator and the endoscope. The MPC operates at 40 Hz, taking the manipulator joint positions and the estimated IPM position and heading as inputs, and outputting joint commands to track the desired trajectory.

Figure 13.

Experimental results of simultaneous magnetic localization and actuation of a magnetic endoscope. (a) The experimental setup for the mock-up test. (b) The control diagram and desired trajectory. (c) Trajectories of the EMS and IPM using the optimization-based method and the PIA. The optimization-based method fails under the same desired trajectory. (d) Position errors of the two methods. (e) Heading error and the calibrated EMS moment. (f) The percentage error of the magnetic field and the range of measurement. (g) With the desired trajectory speed halved, the optimization-based method successfully passes through the colon phantom. (h) Position error using the optimization-based method. (i) The time cost of the two methods.

We perform localization using both the optimization-based approach and our proposed method. Figure 13(c) compares the trajectories obtained using the two localization methods. With the proposed method, the magnetic endoscope completed the entire trajectory. In contrast, the optimization-based method lost control after completing Segment i, leading to task failure. Figure 13(d) shows the position deviation from the ideal trajectory. Both methods exhibit triangular fluctuations, which are caused by stick–slip motion of the magnetically driven IPM, making smooth trajectory tracking difficult. The proposed method maintains a deviation below 2.5 cm throughout the experiment. However, the optimization-based method exhibits much larger fluctuations, with a maximum deviation of 4.5 cm, exceeding the phantom diameter. The resulting inaccurate localization degrades feedback quality and prevents the optimization solver from maintaining stable control. In Figure 13(e), the IPM heading error remains within 10° most of the time. Larger errors, up to 20°, occur in Segments ii and iv, where large rotations are required to change direction. During motion, the EMS dipole moment is calibrated online and fluctuates between 750 and 850 Am². Figure 13(f) presents the magnetic field measurement errors, with measured field magnitudes within 20 Gs. The real-time magnetic field percentage error is computed as the mean field error normalized by the maximum measured field. The proposed method maintains the error below 6.5%, whereas the optimization-based method shows larger fluctuations, reaching up to 9.6%. Such errors account for more than 20% of the maximum IPM field (black curve), severely degrading localization accuracy. The magnetic dipole moment estimated by the proposed method is adaptively adjusted using real-time field feedback, resulting in lower magnetic field errors. Figure 13(i) compares the computational efficiency of the two methods. The proposed method achieves an average computation time of 1.21 ms, whereas the optimization-based method requires 35.5 ms on average. Due to varying iteration numbers, the optimization-based method reaches a maximum computation time of 68 ms, while the proposed method remains consistently below 2.5 ms. As a result, the optimization-based approach cannot meet the 40 Hz control frequency requirement, further contributing to its unstable control performance. The optimization-based method successfully passes through the colon phantom after the desired trajectory speed is halved, as shown in Figure 13(g). The speeds of Segments ii to iv are reduced to allow gradual heading adjustment. A lower trajectory speed reduces the sensitivity of the controller to unstable localization and timing mismatch, preventing rapid error accumulation. However, as shown by the position error in Figure 13(h), when entering Segment iii, the localization becomes noisier at larger sensor distances, which induces oscillations in the EMS during actuation. The experimental results show that our method provides stable and fast localization under EMS actuation and remains robust to trajectory height variations.

5. Conclusion

This article proposes a method for fast localization of the position and heading of an internal permanent magnet (IPM) under the influence of external magnetic sources (EMSs). The method does not require prior knowledge of the dipole moments of the magnets. Through importance analysis, the proposed method can simultaneously localize the IPM and calibrate the dipole moments of both the IPM and EMSs with high efficiency. The proposed adaptive particle sampling method dynamically adjusts the sampling range, enabling efficient computation with sparse particles while maintaining stability. Simulation results demonstrate that, compared to the optimization-based algorithm, the proposed method achieves higher accuracy and better robustness to measurement noise and errors in EMS dipole moments. The proposed method achieves an average computation time of 1.2 ms on a sensor board with 72 sensors. This high computational efficiency enables effective tracking of rapidly moving magnets. We also discuss the difficulty encountered by existing methods in maintaining high accuracy at large distances between the magnet and the sensor board, and experimentally validate the effectiveness of our method under such conditions. For a small magnet (5.2 Am²) at a height of 20 cm, our method achieves an average localization accuracy of 4.70 mm and 4.17 degrees, and the maximum error is only 41% of that of the optimization-based method. Additionally, we experimentally proved the effectiveness of our method under permanent magnets, electromagnets, or the hybrid magnetic field. Compared to the optimization-based method, our method demonstrates higher accuracy and smoother results, particularly when the magnet is subjected to interference from multiple EMSs. Finally, we further validate our method in the simultaneous magnetic localization and actuation of the endoscope by conducting tests in a colon phantom.

The main limitation of our current method is that it still requires the positions of the EMSs. Simultaneous localization of multiple magnets without this information may result in ambiguities, leading to localization errors. Our future work will focus on addressing the challenge of multi-target localization. The particle filter faces difficulties in handling multiple targets, while optimization-based methods are limited by low computational efficiency. Once this issue is resolved, localization based on the magnetic field can be performed without relying on the calibration of EMSs. Moreover, this approach can be extended to a broader range of scenarios that require fast and simultaneous localization of multiple targets.

Supplemental material

Supplemental material

Supplemental material

Supplemental material

Supplemental material

Supplemental material

Supplemental material