Analysis and optimization of surgical electromagnetic tracking systems by using magnetic field gradients

Background: Electromagnetic tracking systems (EMTSs) are widely used in surgical navigation, allowing to improve the outcome of diagnosis and surgical interventions, by providing the surgeon with real-time position of surgical instruments during medical procedures. Objective: the main goal is to improve the limited range of current commercial systems, which strongly affects the freedom of movement of the medical team. Studies are currently being conducted to optimize the magnetic field generator (FG) configuration (both geometrical arrangements and electrical properties) since it affects tracking accuracy. Methods: In this paper, we discuss experimental data from an EMTS based on a developed 5-coils FG prototype, and we show the correlation between position tracking accuracy and the gradients of the magnetic field. Therefore, we optimize the configuration of the FG by employing two different metrics based on i) the maximization of the amplitude of the magnetic field as reported in literature, and ii) the maximization of its gradients. Results: The two optimized configurations are compared in terms of position tracking accuracy,


INTRODUCTION
Reducing invasiveness of medical measurement systems is driving the research toward the development of new systems and techniques to provide high-quality information with a minimallyinvasive approach, such as diagnostic techniques [1], miniaturized biosensors[2], techniques to monitor vital parameters [3]- [5], and tracking systems for surgical navigation [6]. Among these, surgical navigation systems are currently being paid increase attention for the purpose of increasing the tracking distance. Two main type of tracking systems are used in surgery: electromagnetic and optical systems [7].
Electromagnetic tracking systems (EMTS) use a small coil sensor, inserted inside the surgical instrument, to determine its position and orientation by measuring the amplitude of the magnetic field of known geometry generated by a Field Generator (FG) [8]. The small size of the sensor and the independence from a direct line of sight allow to overcome the limitations of optical systems [7], thus allowing the use of flexible instruments such as needles, catheters, and endoscopes [8], [9].
EMTS technology presents two main limitations: a) it has high sensitivity to EM interferences provided by electronic devices and to magnetic field distortions due to metal objects; b) current commercial systems are not able to provide accurate position estimations at a large distance from the signal source, due to the degradation of the magnetic field amplitude with the distance from the FG. Hence, the signal source must be placed too much near the operating volume (i.e., patient's table), thus hindering the medical staff during the operation. Currently, commercial EMTSs provide accurate estimation of surgical instrument position only for limited distance from the FG, generally not more than 0.5 m [8], [10]. It should be noted that a higher tracking distance can be achieved by using some systems, such as Polhemus Long Ranger trackers [11], but they employ larger FGs and bigger sensors, which are not suitable to track surgical instruments such as needles and endoscopes and cannot be used in many applications where a very small size is required.
The design of the FG (electrical and geometrical parameters, arrangements of transmitting coils, processing techniques) determines the distribution of the magnetic field in space and the size of the FG itself, and it affects tracking performances [12].
Often the transmitting coils must be placed inside a welldefined space due to practical needs, such as the configuration of the clinical environment, weight limitations, or application requirements. Moreover, the tracking volume is usually proportional to the size of the FG (i.e., the magnetic field intensity) [8]. Four types of generators can therefore be distinguished [8]: • Standard FGs: they are the most common type and are produced by different manufacturers. They are small and cover an operating volume with linear sizes of a few tens of centimeters. • Flat FGs: they are flat and thin, designed to be placed directly under the operating table where the patient lies. Shielding the back of the FG greatly reduces the sensitivity to interference from below, such as induced eddy currents in metal structures. • Mobile FGs: they are small FGs, which offer a small operating volume. Their advantage is that they can be positioned near the operating area, thus offering good accuracy despite the reduced range. • Long-range FGs: available only from Polhemus Inc., they can cover areas up to 2 m. Different manufacturers propose their own FGs, but they do not provide their design principle. Hence, ongoing research is focused on designing FG's parameters with the aim of increasing tracking volume and accuracy [13]- [15].
In [16] we proposed a solution to study the effect of system parameters and FG configuration on tracking accuracy, and we compared the performance of two FG configurations, whose transmitting coils were identical in their geometrical and electrical parameters, except for their number and their pose (i.e., position and orientation) in space. In particular, one configuration consisted of a flat FG with six coplanar transmitting coils, whereas the other one was a representation of the 5-coils FG prototype presented and characterized in [17]- [19], which was developed to overcome the state of the art of current commercial systems by increasing the tracking distance beyond 0.5 m from the FG. In [16] we did not optimize any FG configurations, but we only compared the performances of chosen configurations.
In [20] the configuration optimization of the fixed sensor array of an EMTS was proposed by employing a mobile transmitting coil, by running two-step evolutionary algorithm and using the position RMSE (Root Mean Square Error) as performance metric. The drawback is that the computational cost is large, and the performance metric is valid only for that specific reconstruction algorithm.
In [12], a simulator to evaluate tracking accuracy -according to certain performance metrics-by setting system parameters was developed. In this study a test volume of 500 mm × 500 mm × 500 mm, which is typical for medical applications was considered. The authors optimized the configuration of a FG composed of eight transmitting coils, for both a compact FG volume or a wide 3D volume, by defining and comparing two objective functions to i) maximize the voltage induced in the sensor (i.e., the magnetic field) or ii) minimize the positional RMSE of the sensor coil within the defined testing volume. The best performances were obtained when choosing the sensor position RMSE as objective function, which however depends on that specific used reconstruction algorithm.
In this paper we consider the same measurement volume as in [12], and we argue that the FG configuration can be optimized by maximizing the gradients of the induced voltage with respect to sensor position (as highlighted from preliminary experimental results), rather than maximizing the induced voltage as done in [12]; therefore, we follow an approach similar to [21], where the authors defined a metric based on the Fisher information matrix. Moreover, positional RMSE in [12] depended on the chosen reconstruction algorithm, whereas our proposed metric is valid regardless of the position reconstruction algorithms.
In [22] we have already discussed the utility of considering the magnetic field gradients to assess system performance, and we supposed position error was influence by gradients, whereas in this paper we perform simulations to employ information from gradients to analyze FG configurations. It should be noted that EMTSs that employ fixed sensor array and a moving emitting coil (instead of a fixed FG) are also present in literature. The approach that we develop in this paper is general, and it is valid for both cases, i.e., i) multiple fixed transmitting coils and one moving sensor, or ii) one moving transmitting coil and multiple fixed sensors.
This paper is structured as follows. In Section 2 we describe the model of the magnetic field obtained by approximating the transmitting coils as magnetic dipoles, which leads to define the sensor induced voltage and its spatial gradient. In Section 3 we discuss experimental data from an EMTS based on a novel 5coils FG prototype, and we show the correlation between position error and the gradients of the magnetic field. In Section 4 we define two metrics based on i) the amplitude of the magnetic field and ii) the gradients of the magnetic field, and we employ them to optimize the arrangement of the transmitting coils of an 8-coils FG. Then, the obtained FG configurations are compared in terms of position tracking accuracy, and the role of the gradient distribution is highlighted. Conclusions are drawn in Section 5.

MODELING THE MAGNETIC FIELD AND ITS SPATIAL GRADIENTS
If the wavelength of the generated magnetic fields is greater than the considered tracking volume, the magnetic field produced by the FG can be calculated by treating the transmitting coils as magnetic dipoles. Therefore, the magnetic moment produced by the i-th transmitting coil can be described as a function of coil parameters and excitation [18], [23], [24]: where ̂, is the versor orthogonal to the surface , of the i-th coil, and , , , and denote the coil radius, the number of turns and the RMS value of the sinewave excitation current, respectively. The same size and current are considered for all transmitting coils in the FG design, but slight differences in real quantity values have been observed. To account for coils' properties mismatch, these parameters are denoted by the subscript .
The RMS magnetic field generated by the i-th transmitting coil in a generic point = [ , , ] T is [18], [23], [24]: where = | |, with = − , represents the vector distance between and the center , of the i-th transmitting coil, and ̂, is its associated versor.
By considering a homogeneous magnetic flux on the surface of the sensing coil , voltage induced by the i-th transmitting coil may be calculated as: where is the number of coil sensor turns, is the excitation current (different for each transmitting coil) and ̂ is the versor orthogonal to the sensor surface (Figure 1), given by We can rewrite (3) in matrix notation: is the × matrix of sensor sensitivities at each frequency, with = 2 π , and is the ×3 matrix of the generated magnetic fields. It is convenient to define = . The gradient of the voltage induced by the i-th transmitting coil with respect to the position of the sensor (i.e., its spatial gradient) can be expressed by [21]:

ANALYSIS OF MAGNETIC FIELD GRADIENTS
In [25] experimental tests have been carried out to evaluate position RMS error of an EMTS based on a novel 5-coils FG prototype based on FDM (Frequency Division Multiplexing), which was designed to reduce mutual inductances between transmitting coils, still maintaining a high induced voltage in the test volume, and in compliance with IEEE standards for safety levels [26], [27]. In this section we refer to those experimental data and we carry out new analyses to show the correlation between position error and the gradients of the magnetic field.

Error propagation based on magnetic field gradients
The standard deviation of repeated measurements of the sensor position = [ , , ] T is calculated to evaluate the position repeatability error. We assume that uncorrelated noise affects the voltage components = [ 1 , … , ] T (n=5 for the examined case) induced in the sensor. The diagonal covariance matrix is: Let g = [g 1 , … , g n ] T : Χ ⫃ R 3 → R n be the function relating and : Expression (10) can be linearly approximated in the neighborhood of a generic point , by means of the Jacobian matrix ( ) ∈ ×3 , whose rows are the gradients ( ), defined as: with = 1, … , , which represents the slope of the function ( ) around along each direction. Expression (11) is a general approximation of (8), and it is independent on the chosen magnetic field model, hence can be used for complex systems of when the field model has not been provided or developed.
As discussed in [22], [25], the variance T of repeated measurements of sensor position can be estimated by: where + ( ) ∈ 3× is the Moore-Penrose pseudoinverse of the Jacobian matrix, ∘2 is a Nx1 column vector containing the element of the diagonal of in (9), ∘ is the Hadamard power operator and ∘ 2 denotes element-wise square. The Moore-Penrose pseudoinverse in (12) provides a least squares solution to the linearized expression of (10). For details, the reader should refer to [22], [25], [28].
Equation (10) is experimentally sampled by moving the sensor on a grid of evenly spaced points , = 1, … , , and measuring voltages . The RMS error, also known as mean radial spherical error in the context of 3D localization, has been determined for each grid point as = √ 2 + 2 + 2 .
The partial derivatives are experimentally estimated by the It should be noticed that this approach is applicable to any position reconstruction algorithm that can be linearized for error propagation. As shown in [22], the same approach can be used to evaluate the effects of residual drift in voltage measurements.

Analysis of the pseudoinverse of the Jacobian matrix
The Jacobian matrix introduced in Section 3.1 gives important information about what to expect in terms of accuracy in different regions of the tracking volume.
In this section we evaluate the elements of the pseudoinverse + of the Jacobian matrix in a volume of 400 mm x 500 mm x 400 mm, for different orientations of the magnetic sensor. For each test, the sensor has been placed by means of the industrial robot in a regular grid of = · · points , = 1, … , , where = 5, = 6, = 5 are the number of points along x-, y-and z-axis, respectively, with a step of 100 mm in each direction, hence each plane normal to z-axis contains · = 30 points, for a total of =150 points. The grid has been scanned by increasing the x-coordinate, then the y-coordinate is increased, and finally the z-coordinate, considering the reference system of Figure 1.
In Figure 2, Figure 3, and Figure 4 the element of + are shown for the sensor oriented along x-, y-, and z-axis, respectively, thus considering , , and components of the magnetic field. The point order has been changed for Figure 3, by increasing first the y-coordinate, then the x-coordinate, in order to highlight the presence of symmetry when switching xand y-coordinate. It can be noted that the elements of + generally increase when the distance from the FG increases, according to magnetic field model. The elements of + related to are much lower than those related to and . By comparing Figure 2 and Figure 3, it can be noted that when the sensor is aligned along x-axis, the elements related to the ycomponent are higher whereas the x-component is lower; viceversa when the sensor is aligned along y-axis; when the sensor is aligned along z-axis all elements assume low values. Moreover, it can be noted the presence of peaks in the elements of + .
It should be said that we conducted a related analysis in [22] by investigating the 3ꞏn elements of + , but in that case only one orientation of the sensor was considered instead of multiple orientations, and the sensor was moved along a single linear trajectory instead of a large volume as done in this case, hence we could not draw any considerations in that case.

Position error evaluation through gradients
In [25] we evaluated the position RMS error of the 5-coils FG prototype, by applying the procedure shown in Section 3.1 in the same volume and by employing the same pseudoinverse values analyzed in Section 3.2. Hence, in this section we refer to those experimental data and we relate the obtained results to the elements of the pseudoinverse of the Jacobian matrix, which is influenced by the arrangement of the transmitting coils of the FG. Figure 5 and Figure 6 show the position error versus the point index, separately for each z-plane, for better visualization. The point order has been changed for Figure 6, by increasing first the y-coordinate, then the x-coordinate, in order to highlight the presence of symmetry when switching x-and y-coordinate.
The error generally increases when increasing the zcoordinate, in accordance with the magnetic field model and with previous results. When the sensor is aligned along the x-axis, we found that the error of the y-component is higher (not shown  here, see [25]), whereas the error of the x-component is lower; vice-versa when the sensor is aligned along y-axis. The error of the z-component (not shown here) is very low in both cases, about one order of magnitude lower. This is in accordance with previous results obtained in Section 3.2.
By considering the elements of + shown in Section 3.2, we can affirm that the increased error obtained for orientation along x-and y-axis is due to higher values of the elements of + (i.e., lower gradients of the magnetic field), rather than to higher voltage noise: in fact, this latter decreases with the distance from the F, whereas the position error increases.
In accordance with the results of Section 3.2, Figure 5 and Figure 6 highlight some regions where the position RMS error is lower than others, as well as the presence of peaks in correspondence of x = 0 mm for Figure 5 (y = ±50 mm for Figure 6) when the sensor is aligned along x-axis (y-axis for Figure 6); these peaks can be here related to peaks in the elements of + . Except for those peaks, position error is always ≤ 0.8 mm, which is still slightly higher than the values that we obtained in [22]: there, in fact, the sensor was more oriented toward the FG, and not in the xy-plane, thus presenting errors similar to the case of the sensor aligned along z-axis.

OPTIMIZATION OF FIELD GENERATOR CONFIGURATION
Following the observations from previous Section 3, in this section we optimize the configuration of the FG by employing two different metrics based on i) the maximization of the induced voltage (i.e., the magnetic field) as in [12], and ii) the maximization of the gradients of the magnetic field, and we evaluate the tracking accuracy associated with the optimized configurations.

Objective functions for optimization procedure
The FG configuration is optimized by minimizing two different cost functions, namely 1 and 2 : • 1 is the same proposed in [12]. It is defined as follows: where is the number of points of the volume considered for the optimization procedure, and is the voltage induced by the i-th transmitting coil when the sensor is in the j-th grid point. In [12] a constant noise level was considered in the whole test volume, and minimizing 1 equals maximizing the measured induced voltage in the test volume.
• 2 is based on the computation of the spatial gradients expressed by (8). We follow an approach similar to [21], where the authors optimized the allocation of measurements for an EMTS based on TDM (time-division multiplexing) by using a cost function based on the Fisher Information Matrix to select a subset of a large number (at least 1089) of fixed coplanar sensors, with dipole moments all oriented along z-axis. They considered a constant 2 in the whole test volume. However, they did not evaluate the tracking performance related to their configurations. Therefore, we define the following cost function 2 : where: is the Fisher Information Matrix (FIM), ∇ ( ) is the spatial gradient of , defined by expression (8), and 2 is the variance of the voltage noise at the j-th point and related to the i-th transmitting coil. Under the assumption that the noise terms are independent, it follows from the Cramér-Rao inequality that maximizing the FIM (in some sense) minimizes the attainable covariance of the estimated position, providing a lower bound [21], [29]: Since it is generally not possible to find an optimal FIM, the real-valued function 2 can be chosen instead [21], [30]. Minimizing 2 is related to the approach followed in Section 3.1 by using expression (12): in fact both the covariance (12) and the inverse of the FIM can be interpreted as the inverse of the Hessian of the weighted least-squares criterion [30].
Note that, in contrast to [21], here the model of the voltage noise is composed of two terms, i.e., 2 = 2 + 2 ( ), where 2 depends on the DAQ device and it is assumed to be constant in the whole test volume, whereas 2 depends on the transmitting coil's frequency and pose, on the noise of the excitation currents, and on the sensor position. The noise values for the simulations have been set according to experimental tests conducted on the 5-coils FG prototype, with an excitation current of 1 A rms for each transmitting coil, thus obtaining noise values acq = 20 nV and = 0.07 mA, as in [19].
MATLAB fmincon function is used to minimize the cost functions thus optimizing the FG configuration, by setting the boundaries of the transmitting coils positions within a box volume of 300 mm × 300 mm × 70 mm. A planar grid of 500 mm × 500 mm is considered at a height from the FG of z = 650 mm, and with a step of 50 mm along x-and y-axis, for a total of 121 positions. At each position, the sensor coil is aligned along the x-, y-, and z-axis, respectively, therefore = 121 × 3 = 363 sensor poses are considered.

Evaluation of position tracking accuracy
To compare different FG configurations, we define a test volume consisting of a grid of points: a volume of 500 mm × 500 mm × 600 mm is considered (from z = 50 mm to z = 650 mm), with a step of 50 mm along x-and y-axis, and 100 mm along z-axis, for a total of 847 positions of the sensor coil. At each position, the sensor coil is aligned along the x-, y-, and z-axis, respectively, therefore = 847 × 3 = 2541 sensor poses are considered. For each point, we evaluate the position error: where , , denote the actual sensor position at the j-th grid point, whereas ̂, ̂, ̂ denote the estimated position by applying a reconstruction algorithm based on the dipole model approximation, as we have also done in [16], [18], [19]. In a real scenario, the position reconstruction algorithm will use the sensor pose estimation of the previous point as starting guess for the estimation of the following one. To simulate a real scenario, a position and orientation perturbation are added to the starting guess pose: for each Cartesian coordinate, values form Uniform distribution (−√3 mm, √3 mm) are added, whereas for polar angle and the azimuthal angle in spherical coordinates we consider (−1 °, 1 °). The Levenberg-Marquardt algorithm is used to solve the position reconstruction problem. For the details, see [19].

Results and discussions
Three different FG configurations are compared, as shown in Figure 7: two of them are composed of eight transmitting coils, whose arrangement in space has been optimized by employing F 1 and F 2 , whereas the other one represents the 5-coils FG prototype developed in [17] and it has been designed to minimize the mutual inductances still maintaining a high induced voltage in the test volume. As reported in [12], eight transmitter coils are commonly applied for real systems; this number of coils has been found elsewhere to be the minimum number of coplanar coils required to continuously track a sensor coil [13]. Hence, eight transmitting coils are considered in this work for a direct comparison with [12]; all transmitting coils are identical in their electrical and geometrical parameters, and are powered with sinusoidal currents of 1 A rms.
In contrast to [12], where the authors considered the same frequency of 1 kHz for all transmitting coils, here we consider different frequencies for each coil, with a frequency gap of 1 kHz according to the values used for first FG prototype, as shown in Table 2; the optimization procedure will take into account this difference. The magnetic field produced by the 8-coils FGs is in compliance with IEEE standards for safety levels [26], [27] beyond 10 cm from the FGs (i.e., z ≥ 10 cm in the FG reference system). Table 1 shows the position error metrics related to the three different FG configurations. It can be noted that the optimized configurations lead to more accurate position tracking with respect to the configuration with five coils. Moreover, optimizing with respect to 2 leads to better results than optimizing with  [17]; b) an 8-coils FG optimized by applying 1 as done in [12]; c) an 8-coils FG optimized by applying 2 . Same scale is applied to all subfigures.
respect to 1 . Note also that both configurations fall within the defined volume, and configuration (b) is more compact than configuration (c), thus resulting in a smaller size of the FG, which is an aspect to consider during system design. Figure 8 shows the value of the mean induced voltage (i.e., 1/ 1 , for better visualization) and of the cost function 2 for each grid point. Only points with index > 363 are considered, i.e., z ≥ 350 mm, since for z ≤ 250 mm the position error does not present significant difference between the different configurations. As expected, the maximum mean induced voltage is achieved for configuration (b), nevertheless, the best position tracking accuracy is obtained with configuration (c), which presents the lowest value of cost function 2 . Configuration (a) and (b) presents similar values of 2 , but the mean induced voltage of configuration (a) is much lower than configuration (b): this is probably responsible for the reduced accuracy of the 5-coils FG. It should be noted that both the mean induced voltage and the cost function 2 assume an oscillating trend, which resemble the behavior of the pseudoinverse of the Jacobian matrix shown in Section 3, thus it is in accordance with previous results.
Overall, the obtained position tracking accuracy is within the requirements of many surgical applications [8] and can be furtherly enhanced by improving the position reconstruction algorithm.

CONCLUSIONS
Commercial EMTSs for surgical navigation present a limited tracking volume, mainly due to the reduced sensitivity of the small magnetic sensors with the distance from the FG. Studies are currently being conducted to optimize the FG configuration of EMTSs for surgical navigation, since it affects tracking accuracy, hence an optimized configuration can increase tracking accuracy in regions far from the FG. In this paper, we discussed the influence of magnetic field gradients on tracking accuracy. Therefore, we optimized the configuration of the FG by employing a metric based on the maximization of the magnetic field gradients, obtaining better performance than the case where the amplitude of the magnetic field is maximized (as done in literature [12]). Overall, position tracking accuracy is within the requirements of many surgical applications. The spatial gradient of the magnetic field produced by the FG assumes an important role in the design of the FG itself, and it paves the way toward the exploitation of new metrics to achieve higher accuracy in regions far from the FG, thus increasing the tracking volume which is the main limitation of current commercial systems. Shim, Continuous glucose monitoring using a microneedle array sensor coupled with a wireless signal transmitter, Sens. Actuators B Chem. 281 (2019), pp. 14-21.