 Research
 Open Access
 Published:
Nonorthogonal onestep calibration method for robotized transcranial magnetic stimulation
BioMedical Engineering OnLinevolume 17, Article number: 137 (2018)
Abstract
Background
Robotized transcranial magnetic stimulation (TMS) combines the benefits of neuronavigation with automation and provides a precision brain stimulation method. Since the coil will normally remain unmounted between different clinical uses, hand/eye calibration and coil calibration are required before each experiment. Today, these two steps are still separate: hand/eye calibration is performed using methods proposed by Tsai/Lenz or Floris Ernst, and then the coil calibration is carried out based on the traditional TMS experimental step. The process is complex and timeconsuming, and traditional coil calibration using a handheld probe is susceptible to greater calibration error.
Methods
A novel onestep calibration method has been developed to confirm hand/eye and coil calibration results by formulating a matrix equation system and estimating its solution. Hand/eye calibration and coil calibration are performed to confirm the pose relationships of the marker/end effector ‘X’, probe/end effector ‘Y’, and robot/world ‘Z’. First, the coil is fixed on the end effector of the robot. During the onestep calibration process, a marker is mounted on the top of the coil and a calibration probe is fixed at the actual effective position of the coil. Next, the robot end effector is moved to a series of random positions ‘A’, the tracking data of marker ‘B’ and probe ‘C’ is obtained correspondingly. Then, a matrix equation system AX = ZB and AY = ZC can be acquired, and it is computed using a leastsquares approach. Finally, the calibration probe is removed after calibration, while the marker remains fixed to the coil during the TMS experiment. The methods were evaluated based on simulation data and on experimental data from an optical tracking device. We compared our methods with two classical methods: the QR24 method proposed by Floris Ernst and the handheld coil calibration method.
Results
The new methods outperform the QR24 method in the aspect of translational accuracy and performs similarly in the aspect of rotational accuracy, the total translational error decreased more than fifty percent. The new approach also outperforms traditional handheld coil calibration of navigated TMS systems, the total translational error decreased three to fourfold, and the rotational error decreased six to eightfold. Furthermore, the convergence speed is improved 16 to 27fold for the new algorithms.
Conclusion
These results suggest that the new method can be used for hand/eye and coil calibration of a robotized TMS system. Two complex steps can be simplified using a leastsquares approach.
Background
Transcranial magnetic stimulation (TMS) is a noninvasive and painless method for stimulating the cerebral cortex nerve [1,2,3]. Based on the principle of electromagnetic induction, an electric current is created on the cerebral cortex using a magnetic coil that is manually placed on top of the patient’s head. Recently, singlepulse and repetitive TMS have been used in clinical study for the therapy of mental disease [4,5,6,7,8,9]. However, TMS is still not widely promoted because its therapeutic effect changes between subjects [10, 11].
The variability is partly because of how accurate the stimulus process is carried out [12,13,14]. Initially, the magnetic stimulation coil was located over the subject’s head manually and imprecisely, without the help of any navigation system [10,11,12]. To solve this problem, neuralnavigation technology has been developed to position the coil on the subject’s head more accurately [13, 14]. Magnetic resonance imaging (MRI) of both the subject and the optical tracking device are used in the neuronavigation system, and navigation software can indicate the actual stimulation point of the coil on the subjects’ brain in real time. However, the stimulation coil, which usually weighs more than 2 kg (Fig. 1), is difficult for the operator to hold for more than 30 min in each procedure. Compared with handheld approach, a robotassisted coil positioning method is more stable and repeatable, which is beneficial to the clinical or academic application of TMS.
Recently, two types of robotic TMS systems have been reported and shown to improve stimulation accuracy [15,16,17,18]. First, a robotic TMS system has been developed using a general industrial robot [15,16,17]. A sixjointed industrial robot (Adept Viper s850) was used in a robotic TMS system designed by Lars et al. [17]. For this kind of robotic system, the coil was mounted to the robot’s end effector, a Polaris Spectra infrared tracking system was used for navigation. After calibration, the magnetic coil is placed quite precisely and directly over a selected target region by the robot. However, the safety of this kind of robotic system has been queried, because the system is equipped with actuators selected for highspeed motions like any other industrial robot [18]. Second, a dedicated robotic system for TMS has been designed based on mechanical architecture and a control strategy [18]. A sevenjointed dedicated robot is designed from the kinematic scheme of the mechanism, including the arm, the prismatic joint, and the wrist [19]. The design and control of the robotic system optimizes the safety of the procedure, and allows the force applied between the robot and the head to be controlled. Robotguided neuronavigated TMS has been used in some recent studies of psychiatric diseases [19,20,21].
Robotized TMS systems are designed and requested for highprecision stimulus, which is primarily determined by calibration [17,18,19,20,21,22]. Besides, in order to ensure the safety of the robotic TMS system, a marker is mounted on the top of the coil to detect the coil position in real time [17,18,19,20,21,22], but the actual effective position of the coil is at the bottom centre of the coil, as shown in Fig. 1. Thus, hand/eye calibration and coil calibration are the most important steps in the workflow of a robotized TMS, and these two steps directly determine the accuracy of the system. Coil calibration involves confirming the first unknown pose matrices between the actual effective positions of the coil and the marker. Traditionally, a handheld coil calibration method is used, as shown in Fig. 1. A marker is fixed on the top of the coil (Fig. 1A), then a probe (Fig. 1C) is placed under the coil and the pointer is aimed at the coil centre labels (Fig. 1B, D). Finally, the positions of the marker and the coil centre in the coordinates of the tracking device can be detected, and the first unknown pose matrices can be calculated. Hand/eye calibration involves accurately computing the other two unknown pose matrices, which are the endeffector/marker matrix and the robot/trackingdevice matrix [22, 23]. Traditional hand/eye calibration involves solving an AX = ZB form matrix equation, where the matrices ‘A’ and ‘B’ are known, and ‘X’ and ‘Z’ are unknown, as shown in Fig. 2D [22,23,24,25,26]. Typically, ‘A’ represents the transform from the robot’s base to the robot’s end effector, and ‘B’ represents the position and orientation data providing the full six degrees of freedom (DOF) of a marker obtained from a tracking device. ‘X’ represents the transforms from the robot’s end effector to the marker, and ‘Z’ represents the transform from the robot’s base to the tracking device.
The handheld coil calibration method typically suffers from larger calibration errors; a translational error of 3 mm is acceptable for handheld coil calibration [27,28,29,30]. In addition, the traditional hand/eye calibration method was reported separately by Shiu and Ahmad [23, 24] and Tsai and Lenz [25, 26]. Matrix algebra and the special properties of homogeneous matrices were used for determining the unknown matrices mentioned above. A review of these calibration algorithms was proposed by Wang et al. [31]. All these methods require that orthogonal homogeneous unknown matrices can be found. To overcome this limitation, the QR24 method, which uses three different variations on the basis of a naïve leastsquares solution of the equation system, was developed for simultaneous hand/eye calibration [22]. However, all the hand/eye calibration methods mentioned above are used for determining two unknown pose matrices. For a robotized navigated TMS system, there are three unknown pose matrices, which need to be solved. The leastsquares method can be applied to get unknown data and to minimize the error sum of squares between the gotten data and the actual data. For this reason, it is not unreasonable to assume that the simultaneous calibration of two steps in a naïve leastsquares solution of the equation system might result in improved accuracy and efficiency.
In order to simplify the two complex calibration steps and improve the accuracy, this report presents an innovative approach for simultaneous hand/eye and coil calibration. A matrix equation system AX = ZB and AY = ZC was acquired at different robot positions, as shown in Fig. 2D. A linear equation system was obtained based on the matrix equation system, and a leastsquares solution was calculated for determining ‘X’, ‘Y’, and ‘Z’. The feasibility and effectiveness of the method were demonstrated by comparing with two classical calibration methods: the QR24 method and the handheld coil calibration method.
Methods
Synchronous hand/eye and coil calibration
In this paper, we present a calibration method for robotized TMS with six DOF, which can be applied to any situation where there are three unknown pose matrices that need to be solved in the robotic system. As shown in Fig. 2A, B, typically, a coil is fixed to the end effector E of the robot R. A marker M with four nearinfrared reflector balls is attached to the top of the coil and, synchronously, a probe P is fixed at the bottom centre of the coil. The mechanical drafting of the probe is shown in Fig. 2C. The tracking device T is placed so that the coil can be adjusted to the centre of the view of T. The design of this new approach to synchronous hand/eye and coil calibration is based on the following relationships:
where the matrix ‘^{R}T_{E}’ which is obtained by forward kinematic represents the transform from the robot’s base to the robot’s end effector, and the matrices ‘^{T}T_{M}’ and ‘^{T}T_{P}’ which are obtained directly from the tracking system represent the tracking data providing full six DOF of the marker and the probe, respectively. These three pose matrices are known parameters. The matrices ‘^{E}T_{M}’ and ‘^{E}T_{P}’ represent the transforms from the robot’s end effector to the marker and the probe, respectively, and ‘^{R}T_{T}’ represents the transform from the robot’s base to the tracking device. These three pose matrices are unknown parameters. To obtain the unknown parameters in Eq. 1, the end effector of the robot is moved to n different positions to obtain Eq. 2.
Generally, the position of the robot’s end effector in Eq. 2 is chosen randomly in a sphere of radius r. The rotation angles of yaw, pitch, and roll for each position are also selected randomly between ± d degrees.
Let A_{i}= (^{R}T_{E})_{i}, X = ^{E}T_{M}, B_{i} = (^{T}T_{M})_{i}, Y = ^{E}T_{P}, C_{i} = (^{T}T_{P})_{i}, and Z = ^{R}T_{T}. Then we obtain Eq. 3:
Let
Thus, a linear system of equations can be derived from Eq. 3:
where \( D \in \mathbb{R}^{24n \times 36} \;{\text{and }}b \in \mathbb{R}^{24n}.\) To be more specific,
where
The general form of the homogenous transformation matrix ‘A _{ i} ^{−1} ’ can be presented as
where \( \text{R} [A_{i}^{  1} ]\in {\mathbb{R}}^{ 3\times 3} \) represents the rotation matrix of ‘A _{ i} ^{−1} ’ and \( T [A_{i}^{  1} ]\in {\mathbb{R}}^{ 3} \) represents the translation vector of ‘A _{ i} ^{−1} ’. ‘Zero_{m × n}’ is a zero matrix of m × n size, and ‘Eye_{k}’ is the identity matrix of k × k size. A leastsquares method with QRfactorization can be used to solve a linear system of equations such as Eq. 4 [22]. In the definition of linear algebra, QRfactorization is a decomposition of the matrix ‘D’ into a result ‘D = QR’, where ‘Q’ is an orthogonal matrix and ‘R’ is an upper triangular matrix. All elements of unknown matrix ‘X’, ‘Y’ and ‘Z’ can be determined from the solution vector w of Eq. 4 based on the decomposed matrix ‘Q’ and ‘R’. All elements of unknown matrices ‘X’, ‘Y’ and ‘Z’ can be determined from the solution vector w of Eq. 4. It is important to point out that this method gives the optimum solutions for ‘X’, ‘Y’ and ‘Z’; thus, they minimize Eq. 7
where *_{F} is the Frobenius norm. The method presented above is called the QR36 calibration method, because a linear system of equations with 36 unknown parameters can be solved by QRfactorization.
To use these matrices, orthonormalization is the final and most important step, because ‘X’, ‘Y’ and ‘Z’ are not orthogonally solved in the QR36 method. Without this step, the homogeneous coordinate transformations obtained from ‘X’, ‘Y’ and ‘Z’ would be incorrect. In this paper, the singular value decomposition (SVD) method was used for orthonormalization. For a position N obtained from a tracking device, for which we want to calculate the position in robotic coordinates, we compute ZN to obtain a nonorthogonal matrix; next, the SVD of ZN is calculated as UΣV^{T} = ZN. Finally, the orthonormalized ZN can be calculated with (ZN)^{⊥ }= UV^{T}.
Calibration errors
With the calibration matrices ‘X’, ‘Y’ and ‘Z’ calculated by the QR36 method, we can obtain the calibration error using Eq. 1:
where ‘E_{M}’ is defined as the calibration error matrix of the marker, and ‘E_{p}’ is defined as the calibration error matrix of the probe. To obtain the calibration quality, we define the translational error and rotational error based on Eqs. 9–13. Let O_{M}= E _{ M} ^{⊥} and O_{p}= E _{ p} ^{⊥} . Thus, the translational error of the system for QR36 method is defined as
where e_{t} is defined as the total translational error of the system, e_{trans}[O_{p}] is defined as the translational error of the probe, e_{trans}[O_{M}] is defined as the translational error of the marker.
Let (k_{M}, θ_{M}) and (k_{p}, θ_{p}) are the axisangle representations of R(O_{M}) and R(O_{p}), respectively. Thus, the total rotational error of the system for QR36 method is defined as
where θ_{t} is defined as the total rotational error of the system, θ_{p} is defined as the rotational error of the probe, θ_{M} is defined as the rotational error of the marker.
Finally, if the marker and probe are calibrated separately using the QR24 method [22], then the calibration matrix Z will not be identical for both optimal calibration procedures. Thus, another translational error e_{z} is defined for the system:
where ‘Z_{p}’ and ‘Z_{M}’ are the calibration results of the probe and the marker calculated using the QR24 method. Thus, the translational error of the system for QR24 method is defined as
where e_{t} is defined as the total translational error of the system, e_{trans}[O_{p}] is defined as the translational error of the probe, e_{trans}[O_{M}] is defined as the translational error of the marker.
Let (k_{z}, θ_{z}) be the axisangle representation of R(Z_{p} − Z_{M}); thus, the rotational error is defined as θ_{z}. Thus, the rotational error of the system for QR24 method is defined as
where θ_{t} is defined as the total rotational error of the system, θ_{p} is defined as the rotational error of the probe, θ_{M} is defined as the rotational error of the marker.
Data acquisition
Simulated and optical tracking data were obtained to estimate the accuracy and robustness of the QR36 method. To compare the quality of our method with the QR24 and handheld methods, the QR24 method and the QR36 method were performed on the same dataset. The results are presented in the following section.
Simulated data
A set of simulated data was generated to test the proposed method [22]. Random orthogonal matrices X, Y, and Z with same definition in Eq. 3 were created first. To generate realistic values, the elements of T[X] and T[Y] were chosen randomly from [− 50, 50] mm, and the elements of T[Z] were chosen randomly from ±[500, 2000] mm. Subsequently, five hundred totally random orthogonal robot positions matrix A_{i} were created and used to calculate the corresponding tracking matrices B_{i} and C_{i}. To get the simulated data of an optical tracking device, a marker and a probe with four markers was defined as
In order to distort the tracking matrices \( \tilde{B}_{i} \) and \( \tilde{C}_{i} , \) the marker and the probe were moved to the pose described by B_{i} and C_{i}, and according to Horn’s algorithm, the matrix was computed without scaling [26], that is,
Here,
where s_{i, j, k} was obtained from the Gaussian distribution with standard deviation equal to 0.05 mm.
A different method was used to distort the robot matrices: the rotational matrix was distorted by multiplication with a rotation matrix around a random axis a_{i}. The rotation angles θ_{i} was drawn from the Gaussian process with standard deviation equal to 0.05 mm. Gaussian noise was used to distort the translational vector:
where
s_{i, k} was obtained from the Gaussian distribution with standard deviation equal to 0.05 mm.
Data from optical tracking device
As shown in Fig. 2B, C, an optical tracking device (Polaris Spectra, Northern Digital Inc., Waterloo, Ontario, Canada) was mounted on tripods, and a coil model was fixed at the end effector of a 6axis robot (VS060, Denso Inc., Aichi, Japan). The marker was attached to the top of the coil, and the probe was fixed at the bottom centre of the coil. Then, the position of the end effector and the tracking device were adjusted so that the positions of the marker and probe were near the centre of the tracking device’s work volume (that is, x < 100 mm, Y < 100 mm, and Z = − 1500 ± 100 mm for both the marker and the probe), and the position of the end effector was defined as first position P_{0}. Next, the end effector was moved to 500 random positions around the first position P_{0}. The positions were confirmed by
where R_{x, y and z} (θ) is a rotational matrix around the x, y and zaxis, in which θ was randomly selected from − 10° to 10°, and T is a translation matrix with translation vector t that was randomly selected in the range from − 100 to 100 mm. 100 measurements (A, B, C) were averaged to reduce the errors from the robot and optical tracking device.
Results
Evaluation and implementation
The algorithms were evaluated on a standard personal computer (Intel core E7500 CPU with 4 GB of RAM, running 64bit Windows 10 OS). The algorithms were all implemented in MATLAB 2010a. The performance period was less than 0.1 s for the QR24 and QR36 algorithms.
Calibration errors
The QR36 algorithm was performed using poses P_{1},…,P_{n} for n = 5,…,250. The rotational and translational errors were computed for each remaining testing pose P_{251},…,P_{500} using Eqs. 9 and 10, and the average of all 250 testing poses was determined. The QR24 algorithm was calculated on the marker and the probe separately using poses P_{1},…,P_{n} for n = 5,…,250. The rotational and translational errors of the marker and the probe were computed for each remaining testing pose P_{251},…,P_{500} using Eqs. 12 and 13, and the average of all 250 testing poses was determined. Then, error E_{z} and θ_{z} were determined for each remaining testing pose P_{251},…,P_{500}, and the average of all 250 testing poses was determined.
Simulated results
Figure 3 presents the average errors for the simulated data using the QR24 method (top) and the QR36 method (bottom).
Clearly, using more than 20 poses hardly influences the translational and rotational errors of the QR36 method. However, using the QR24 method to calibrate the marker and probe separately could create two different calibration matrices Z. This can create a new calibration error e_{z}, which significantly influences the translational error of the QR24 algorithm. In the simulated data, the resulting error came very close to the error obtained using the correct matrices, and remained within 1%. The minimum total translational error decreased 59.54% for the simulated data. The minimum total translational error e_{t} of the QR24 method was 0.3419 mm, with an e_{z} of 0.1267 mm, while the minimum total error of the QR36 method was 0.2143 mm. In other words, without error e_{z}, the minimum total errors of the QR24 method and the QR36 method are approximately equal; the same result was obtained from the data for the 25th, median, and 75th errors, but not the maximum total error. Insufficient sample data produced an extremely large error e_{z}; therefore, our method is more useful for a small data set. Overall, both convergence speed and precision of the translational error were improved for the QR36 method compared with the QR24 method. However, although the QR36 method outperformed the QR24 method regarding the rotational error of θ_{t}, the effect of θ_{z} was very slight. The θ_{z} of the QR24 method only had a 0.78% influence on the minimum total rotational error.
Additionally, Table 1 presents the statistics of the errors for the simulated data shown in Fig. 3. The minimum, 25th, median, 75th, and maximum total translational error e_{t} and total rotational error θ_{t} and corresponding parameters are summarized in Table 1.
Experimental results
Figure 4 presents the average errors for the optical tracking device. We obtained similar results from the experimental data. The translational errors of the QR24 and QR36 methods were hardly changed by using more than 30 poses. In the experimental data, the resulting error became much larger than the error obtained from simulation. We also found that the error e_{z} had a greater influence on the translational error of the experimental data. The minimum total translational errors were decreased by 54.93%, as 0.9353 mm and 0.6037 mm for QR24 and QR36, respectively. e_{z} had a minimum total translational error of 0.3301 mm. The effect of θ_{z} on the rotational error of θ_{t} was very small, according to the experimental data. The maximum total translational errors were 8.3719 and 1.2384 mm for QR24 (P_{n}= 5) and QR36 (P_{n}= 6), respectively. e_{z} had a translational error of 5.7679 mm, indicating that fewer tracking data will lead to a larger calibration error e_{z}.
Additionally, Table 2 presents the statistics of the errors for the experimental data shown in Fig. 4. The minimum, 25th, median, 75th, and maximum total translational error e_{t} and total rotational error θ_{t} and corresponding parameters are summarized in Table 2.
Convergence properties. To acquire the convergence properties of the two algorithms, the simulated and experimental data were assessed again. It is clear from Figs. 3 and 4 that the QR36 method converged faster than the QR24 method. Additionally, it was possible to determine whether and how fast the algorithms converged to set values, which were presented in this paragraph. The results are shown in Table 3. Again, we can see that the convergence speed is improved 16 to 27fold for the QR36 algorithms, stabilizing at a translational error of less than 0.4 mm and 0.3 mm at P_{n}= 7 and 14, respectively. For the QR24 method, translational errors of < 0.4 and 0.3 mm were found at P_{n}= 118 and none, respectively. The translational errors of the experimental data followed a similar pattern: the convergence properties of the QR36 method were much better than those of the QR24 method. However, the rotational errors showed a different picture: the accuracy of less than some set values is achieved on both algorithms similarly.
Workspace sizes
Workspace sizes have a substantial influence on the calibration error in many medical robotic applications. To demonstrate the quality of our method in this regard, additional calibration experiments with the Polaris Spectra system and increased workspace sizes were performed. Five different workspace sizes were chosen for the experiments (100 mm, 10°; 200 mm, 10°; 300 mm, 10°; 300 mm, 15°; and 300 mm, 20°). In all cases, 500 samples were successfully recorded. The minimum errors for different workspace sizes are presented in Fig. 5. Interestingly, the total translational error of the QR36 method increased with workspace size. However, the total translational error of the QR24 method was influenced by error e_{z}. e_{z} did not increase linearly with increasing workspace size; however, it did create a smaller calibration error for a large workspace (e_{z}= 1.04 mm for a workspace of 300 mm, 20°). The total rotational error of the QR24 and QR36 methods all increased with increasing workspace size, but error e_{z} had a small influence.
Handheld coil calibration
To compare the quality of our method with the handheld coil calibration method, two operators with similar experience of TMS participated in this experiment. The QR36 method and handheld calibration method were performed by the two operators. The translational error e_{trans}[O_{p}] and rotational error θ_{p} of the probe, which were defined in Eqs. 9 and 10, were estimated to evaluate the calibration accuracy.
A set of five additional data acquisition experiments using the robotic TMS system in a workspace of 300 mm (10°) were performed by each operator. For the experimental procedure, refer to the experimental data acquisition paragraph in the Methods section. For each experiment, the probe was remounted on the coil by the operator, and 500 measurements (A_{i}, B_{i}, C_{i} _{i = 1,…,500}) were collected, with the end effector moving to 500 random positions around the first position P_{0}. Next, the robot was moved to the first position P_{0}, and the probe was unloaded from the coil and held at the bottom centre labels of the coil by the operator. An additional measurement (B_{a}, C_{a}) was acquired with an optical tracking device. The acquired 500 measurements and the additional measurement were defined as the dataset of each experiment.
The QR36 algorithm was performed using poses P_{1},…,P_{n} for n = 5,…,250 of each experiment. The rotational and translational errors of the probe were computed for each remaining testing pose P_{251},…,P_{500} using Eqs. 9 and 10, and the average of all 250 testing poses was determined. Then, the optimal translational and rotational calibration errors (MIN (e_{trans}[O_{p}])_{n} and MIN(θ_{p})_{n} for n = 5,…,250) were determined for each experiment. Finally, the mean value and standard deviation of the optimal calibration errors for five experiments per operator were determined. The results are presented in Fig. 6.
The traditional separate twostep calibration with QR24 and the handheld method were carried out for comparison with the QR36 method. Based on the dataset of marker (A_{i}, C_{i} _{i = 1,…,500}), the optimal calibration result (X, Z) was determined by the QR24 method with the minimal translational error. Subsequently, Y was calculated using Eq. 13 with the additional position pair (B_{a}, C_{a}) [24]. For each experiment, the translational error e_{trans}[O_{p}] and rotational error θ_{p} of the probe were computed for each testing pose (B_{i}, C_{i} _{i = 251,…,500}) with the calibration result (Y, Z), and the average of all 250 testing poses was determined as the handheld calibration error of the experiment. Finally, the mean values and standard deviations of calibration errors of five experiments for each operator were determined. The results are presented in Fig. 6. All the symbols (A, B, C, X, Y, Z) in this section are as defined in Eq. 3.
As shown in Fig. 6, both the mean values and the standard deviations of the rotational and translational errors significantly decreased when the QR36 method was used; the total translational error decreased three to fourfold, and the rotational error decreased six to eightfold. For the two operators, the minimum total translational errors for QR36 method were decreased to 0.7931 and 0.8361 mm, respectively, as compared with 2.4997 and 2.5354 mm for handheld method. The minimum total rotational errors were decreased to 0.5509° and 0.5624° using QR36 method, as compared with 4.2808° and 3.6202° using handheld method. The rotational error showed greater improvement, probably because there was still a positioning label on the coil for the position probe, whereas there is no method for reducing the perspective error of the operator. Moreover, because the probe was fixed on the coil mechanically and an optimization algorithm was used, no significant differences in either the translational or the rotational errors of the QR36 method were found for different operators. Finally, the QR36 method showed better robustness to choice of operator than the handheld coil calibration method.
Discussion
An innovative approach for synchronous hand/eye and coil calibration has been presented and evaluated. A leastsquares method that has been widely used in medical imaging and robotics was used to estimate the optimal calibration matrix [32, 33]. The approach has been validated with synthetic data and experimental data from an optical tracking device. The calibration effect of our method was compared with that of traditional methods such as QR24 and the handheld coil calibration method.
These results show that QR36 calibration method is advisable for use in the robotized TMS system. We have shown that both the QR24 method and our new QR36 method perform very well, with errors below 1 mm and 0.2°, using both simulated and experimental data. This shows that the two methods can both be used for typical calibration. However, whereas more than 200 different robot positions were required to achieve a calibration error below 1 mm for the QR24 method, we found that fewer than 20 positions were typically required by our method. This means that the QR36 method is especially suitable to handle cases where fewer tracking data are available.
Traditional handeye calibration method proposed by Shiu and Ahmad [23, 24] and Tsai and Lenz [25, 26] calculates the rotational and translational parts of the unknown matrices separately using matrix algebra. Li and Betsis applied a geometric approach and leastsquares solution for handeye calibration [34]. Dual quaternion approach was also used for handeye calibration [35]. All those methods expect that orthogonal homogeneous matrices can be found. Moreover, a robust realtime handeye calibration method was proposed by Lars and Floris present [17]. A marker is attached to the robot’s third link for realtime handeye calibration. However, this method is not as precise as the QR24 algorithm. The total translational errors were 0.88 mm and 1.36 mm for QR24 method and realtime calibration method, respectively [17]. A robot system will not be calibrated perfectly, so orthogonality is not necessary. It is accepted or requested to permit nonorthogonality of the matrices in our calibration method. Our results demonstrated that the calibration method used in this study is more accurate than the classical handeye calibration approaches [23, 25]. In terms of translational accuracy, our method also outperforms the QR24 method. Therefore, the calibration method proposed in this study is more suitable for robotized transcranial magnetic stimulation.
We should point out that the maximum optimal translation errors for five handheld experiments per operator can reach 3 mm, as shown in Fig. 6a. This result is acceptable for a handheld navigated TMS experiment. During the robotassisted TMS stimulation, the head motion is tracked using the marker on the subject’s head and compensated by the robot [18]. But, without robotic assistance and head motion compensation, the relative motion between the subject’s head and the handheld stimulation coil is greater than 3 mm during handheld TMS experiments [30]. Moreover, the calibration error of registration of the subject’s MRI and optical tracking device is also controlled at around 3 mm for handheld navigated TMS experiments [36]. However, for a robotic TMS system, which is designed for highprecision stimulation, the handheld coil calibration method is not suitable [37,38,39].
It is important to note that the new method has three limitations. It can only be applied when:

It is accepted or requested to permit nonorthogonality of the matrices in our calibration method;

The calibrated robotic TMS system is used in the same space where calibration was carried out;

There are three unknown matrices that need to be solved in the robotic system.
Conclusion
We have developed a new onestep calibration method to acquire three pose relationships from a navigated robotic system with a leastsquares approach. The new method can significantly improve the accuracy of the robotic TMS system. Besides, the convergence speed is improved for the new method, which means that our method is particularly suited to handle fewer tracking data. The capability of the new method has been demonstrated for synchronous calibration and determination of the pose relationship of marker/end effector, probe/end effector, and robot/world for a robotic TMS system, which can be used to perform precision TMS experiment. Finally, for many robotic applications, where three unknown matrices need to be solved, the method presented in this paper provides an alternative solution to the classical approaches. In the future, it will be interesting to quantitatively compare the stimulation effect of robotic and manual techniques. Investigation on patients will then be pursued to evaluate the medical benefits of robotic TMS system.
References
 1.
Barker AT, Jalinous R, Freeston IL. Noninvasive magnetic stimulation of human motor cortex. Lancet. 1985;1:1106.
 2.
Hallett M. Transcranial magnetic stimulation. Nature. 2000;406:147–50.
 3.
Hallett M. Transcranial magnetic stimulation: a primer. Neuron. 2007;55:187–99.
 4.
Avenanti A, Bueti D, Galati G, et al. Transcranial magnetic stimulation highlights the sensorimotor side of empathy for pain. Nat Neurosci. 2005;8:955–60.
 5.
Fitzgerald PB, Brown TL, Marston NAU, et al. Transcranial magnetic stimulation in the treatment of depression: a doubleblind, placebocontrolled trial. Am J Psychiatry. 2003;160:835.
 6.
Hirayama A, Saitoh Y, Kishima H, et al. Reduction of intractable deafferentation pain by navigationguided repetitive transcranial magnetic stimulation of the primary motor cortex. Pain. 2006;122:22–7.
 7.
Hoffman RE, Cavus I. Slow transcranial magnetic stimulation, longterm depotentiation, and brain hyperexcitability disorders. Am J Psychiatry. 2002;159:1093.
 8.
Mcnamara B, Ray JL, Arthurs OJ, et al. Transcranial magnetic stimulation for depression and other psychiatric disorders. Psychol Med. 2001;31:1141–6.
 9.
PascualLeone A, Rubio B, Pallardó F, et al. Rapidrate transcranial magnetic stimulation of left dorsolateral prefrontal cortex in drugresistant depression. Lancet. 1996;348:233.
 10.
Herwig U, Padberg F, Unger J, et al. Transcranial magnetic stimulation in therapy studies: examination of the reliability of “standard” coil positioning by neuronavigation. Biol Psychiatry. 2001;50:58–61.
 11.
Lisanby SH, Kinnunen LH, Crupain MJ. Applications of TMS to therapy in psychiatry. J Clin Neurophysiol. 2002;19:344–60.
 12.
Sparing R, Buelte D, Meister IG, et al. Transcranial magnetic stimulation and the challenge of coil placement: a comparison of conventional and stereotaxic neuronavigational strategies. Hum Brain Mapp. 2008;29:82.
 13.
Herwig U, Schönfeldtlecuona C, Wunderlich AP, et al. The navigation of transcranial magnetic stimulation. Psychiatry Res Neuroimaging. 2001;108:123–31.
 14.
Neggers SFW, Langerak TR, Schutter DJLG, et al. A stereotactic method for imageguided transcranial magnetic stimulation validated with fMRI and motorevoked potentials. Neuroimage. 2004;21:1805.
 15.
Kantelhardt SR, Fadini T, Finke M, et al. Robotassisted imageguided transcranial magnetic stimulation for somatotopic mapping of the motor cortex: a clinical pilot study. Acta Neurochir (Wien). 2010;152:333–43.
 16.
Lancaster JL, Narayana S, Wenzel D, et al. Evaluation of an imageguided, robotically positioned transcranial magnetic stimulation system. Hum Brain Mapp. 2004;22:329.
 17.
Lars R, Floris E, Alexer S, et al. Robust realtime robotworld calibration for robotized transcranial magnetic stimulation. Int J Med Robot. 2011;7:414.
 18.
Zorn L, Renaud P, Bayle B, et al. Design and evaluation of a robotic system for transcranial magnetic stimulation. IEEE Trans Biomed Eng. 2012;59:805.
 19.
Yi X, Bicker R. Design of a robotic transcranial magnetic stimulation system. Robotics automation and mechatronics. IEEE, 2018. p. 78–83.
 20.
Hartmann A, Kalis R, Matthäus L, et al. P 226. Robot guided positioning of a magnetic coil for rTMS over the cortex in patients with cerebral lesions. Clin Neurophysiol. 2013;124(10):e173–4.
 21.
Quesada C, et al. Robotguided neuronavigated repetitive transcranial magnetic stimulation (rTMS) in central neuropathic pain. An update of longterm followup. Arch Phys Med Rehabilit. 2018. https://doi.org/10.1016/j.apmr.2018.04.013.
 22.
Floris E, Lars R, Lars M, et al. Nonorthogonal tool/end effector and robot/world calibration. Int J Med Robot. 2012;8:407–20.
 23.
Shiu Y, Ahmad S. Finding the mounting position of a sensor by solving a homogeneous transform equation of the form AX = XB. In: IEEE international conference on robotics and automation proceedings. 1987. p. 1666–71.
 24.
Shiu YC, Ahmad S. Calibration of wristmounted robotic sensors by solving homogeneous transform equations of the form AX = XB. IEEE Trans Robot Autom. 1989;5:16–29.
 25.
Tsai RY, Lenz RK. A new technique for fully autonomous and efficient 3D robotics hand/eye calibration. In: International symposium on robotics research. 1988. p. 287–97.
 26.
Tsai RY, Lenz RK. A new technique for fully autonomous and efficient 3D robotics hand/eye calibration. IEEE Trans Robot Autom. 1989;5:345–58.
 27.
Matthaus L. A robotic assistance system for transcranial magnetic stimulation and its application to motor cortex mapping. Ph.D. thesis, Universität zu Lübeck; 2008.
 28.
Richter L. Robotized transcranial magnetic stimulation. Dodrecht: Springer; 2013. p. 1794.
 29.
Horn BKP. Closedform solution of absolute orientation using unit quaternions. J Opt Soc Am A. 1987;4(4):629–42.
 30.
Richter L, Trillenberg P, Schweikard A, et al. Stimulus intensity for hand held and robotic transcranial magnetic stimulation. Brain Stimul. 2013;6(3):315–21.
 31.
Wang CC. Extrinsic calibration of a vision sensor mounted on a robot. IEEE Trans Robot Autom. 2002;8(2):161–75.
 32.
Lu H, Li X, Hsiao IT, et al. Analytical noise treatment for lowdose CT projection data by penalized weighted leastsquare smoothing in the KL domain. Chin J Med Phys. 2002;4682(4):146–52.
 33.
Wang T, Nakamoto K, Zhang H, et al. Reweighted anisotropic total variation minimization for limitedangle CT reconstruction. IEEE Trans Nuclear Sci. 2017;64(10):2742–60.
 34.
Li M, Betsis D. Headeye calibration. In: Proceedings of the 5th international conference on computer vision (ICCV’95). 1995. p. 40–5.
 35.
Daniilidis K. Handeye calibration using dual quaternions. Int J Robot Res. 1999;18(3):286–98.
 36.
Herwig U, Schönfeldtlecuona C, Wunderlich AP, et al. The navigation of transcranial magnetic stimulation. Psychiatry Res Neuroimaging. 2001;108(2):123–31.
 37.
Kantelhardt SR, Fadini T, Finke M, et al. Robotassisted imageguided transcranial magnetic stimulation for somatotopic mapping of the motor cortex: a clinical pilot study. Acta Neurochir. 2010;152(2):333.
 38.
Wan Z, Wan N. Forcecontrolled transcranial magnetic stimulation (TMS) robotic system. In: 2010 8th world congress on intelligent control and automation (WCICA). IEEE. 2012. p. 377781.
 39.
Richter L, Bruder R, Schweikard A. Handassisted positioning and contact pressure control for motion compensated robotized transcranial magnetic stimulation. Int J Comput Assist Radiol Surg. 2012;7(6):845.
Authors’ contributions
HW and JNJ conceived and designed the study. YL and XW performed the experiments. HW wrote the paper. HW and TY reviewed and edited the manuscript. All authors read and approved the final manuscript.
Acknowledgements
H.W. acknowledges the grant support from PUMC Youth Fund and the Fundamental Research Funds for the Central Universities (2017320025). T.Y. acknowledges grant support from CAMS Initiative for Innovative Medicine (CAMSI2M 2016I2M1004) and the Natural Science Foundation of China (No. 81772003).
Competing interests
The authors declare that they have no competing interests.
Availability of data and materials
Not applicable.
Consent for publication
At the time of their initial briefing, all study participants were informed of the likelihood that the data would be part of a publication.
Ethics approval and consent to participate
Written informed consent was obtained from the subject for the publication of this report and any accompanying, in accordance with the regulations of the local ethics committee (Committee on Health Research Ethics, Chinese Academy of Medical Sciences).
Funding
CAMS Initiative for Innovative Medicine (CAMSI2M 2016I2M1004). Natural Science Foundation of China NSFC (No. 81772003). PUMC Youth Fund and the Fundamental Research Funds for the Central Universities (2017320025).
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Author information
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 Transcranial magnetic stimulation
 TMS
 Medical robots and systems
 Hand/eye and coil calibration
 Leastsquares approach