Versatile robotic probe calibration for position tracking in ultrasound imaging

Within the field of ultrasound-guided procedures, there are a number of methods for ultrasound probe calibration. While these methods are usually developed for a specific probe, they are in principle easily adapted to other probes. In practice, however, the adaptation often proves tedious and this is impractical in a research setting, where new probes are tested regularly. Therefore, we developed a method which can be applied to a large variety of probes without adaptation. The method used a robot arm to move a plastic sphere submerged in water through the ultrasound image plane, providing a slow and precise movement. The sphere was then segmented from the recorded ultrasound images using a MATLAB programme and the calibration matrix was computed based on this segmentation in combination with tracking information. The method was tested on three very different probes demonstrating both great versatility and high accuracy.


Introduction
When using ultrasound in image-guided therapy, it is essential to know the position of the ultrasound images in space; i.e. where in space the objects that appear in the images are located. This is necessary both to create three-dimensional volumes from the images and subsequently to navigate on those volumes. To determine its position, the ultrasound probe is usually equipped with a position sensor, whose position is measured in real time. The most common types of position sensors are optical sensors consisting of multiple infrared light sources which are tracked by infrared cameras and electromagnetic sensors consisting of small coils whose positions can be determined by setting up a controlled, varying magnetic field and measuring the voltages induced in the coils. Since the position sensor is fixed to the probe, the position of the ultrasound images relative to the sensor is constant and can be found in a process referred to as probe calibration. By combining the result of this calibration with the real-time measurements of the sensor, the position of each ultrasound image can be determined.
Substantial research has been done in the field of probe calibration in the last two decades and this has resulted in a number of fast, automatic and accurate calibration methods. Thorough reviews of this work are given in Mercier et al (2005) and Hsu et al (2009). While these methods are usually developed for a specific probe, they can in principle easily be adapted to other probes. In practice, however, this is not the case. Ultrasound probes are becoming more and more specialized and tailored to an increasing number of applications. The differences between probes are large with respect to properties such as shape, field of view, resolution, contrast and noise. As a result, the adaptation of calibration methods often proves tedious, requiring modifications to several central components, especially phantoms and image processing algorithms. A side-looking probe may, for instance, not be able to get to the surface of a phantom made for an ordinary forward-looking probe; a probe with a small field of view may not be able to image all of the wires in a typical wire phantom; and a highresolution probe imaging a bead phantom may produce quite different reverberation artefacts than an ordinary probe.
This tedious adaptation is impractical in a research setting, where ultrasound guidance is regularly applied to new medical fields, all of which have their own more or less specialized ultrasound equipment. Particularly for initial trials and small-scale feasibility studies, the need for a customized phantom is an obstacle. Therefore, we have developed a method which is not particularly fast, nor completely automatic, but which is accurate and robust and, most importantly, can be applied to a large variety of probes without adaptation.

Theoretical background
To fully understand what a probe calibration does, it is helpful to know the various coordinate systems that are involved. These are illustrated in figure 1. In the following, M r ← s denotes a 4-by-4 transformation matrix that, through multiplication, converts the coordinate vector p s of a given point in coordinate system s into the coordinate vector p r of the same point in coordinate system r, i.e.
First, a position sensor is usually attached to the patient or the operating table to act as a fixed reference and this sensor defines the reference coordinate system r. This means that the positions of all other objects, such as the instruments, images and the patient him/herself, are given in this system. Then, the position sensor attached to the ultrasound probe has its own coordinate system s. The tracking system continuously measures the positions of these two first coordinate systems and calculates the spatial relationship between them, i.e. the rigid transformation M r ← s converting position sensor coordinates into reference coordinates. Finally, the two-dimensional ultrasound images are defined in a third coordinate system i. The spatial relationship between this system and the one defined by the sensor on the probe is fixed and the goal of the probe calibration is to find this relationship in the form of a rigid transformation M s ← i . The basic principle of most probe calibration methods is to image an object whose appearance in the ultrasound images makes it easy to accurately measure its position within the image. This is referred to as an imaging phantom and it usually includes one or more imaging targets, which are features that are easily identified in the images. To enable the ultrasound imaging, the phantom must be built in or submersed in an acoustic coupling medium such as water. It is also equipped with a reference sensor r and the position p r of the imaging target relative to this reference is measured, usually using some kind of tracked pointer. As the ultrasound images of the phantom are recorded, the positions of both the phantom and the ultrasound probe are captured by the tracking system and saved together with the images. The imaging target's position p i in the coordinate system i is also extracted from the ultrasound images and the calibration matrix can then be found by minimizing the distance between the target positions in the two coordinate systems, i.e. as A lot of different phantoms have been proposed, but most of them belong to one of the following five groups (Mercier et al 2005): Point target phantoms: Here, a single, small object serves as an imaging target, e.g. a small bead (State et al 1994, Lindseth et al 2003, Brendel et al 2004, Wang et al 2011, a wire cross (Detmer et al 1994, Meairs et al 2000, Yaniv et al 2011 or the tip of a stylus   (Khamene and Sauer 2005, Hsu 2007, Lang et al 2012. These phantoms are simple and thus easy to make, but they only provide one datapoint per ultrasound image and must therefore be imaged many times for each calibration. The small size of the bead or wire also means that its appearance in the ultrasound image varies a lot between probes; what appears as a focused and distinct spot in one image may appear blurred and noisy when imaged with a different system. A target that is suitable for one probe may therefore not be suitable for another. Multiple point targets phantoms: These are similar to the point target phantoms, except they include several small objects that are to be imaged either successively (Trobaugh et al 1994, Sauer et al 2001 or simultaneously (Lindseth et al 2003, Leotta 2004). This limits the number of recordings required for each calibration. In the last case, however, there is the additional problem of properly aligning the ultrasound plane with the targets in order to image them all at the same time. Z-fiducial phantoms: These phantoms include one or more Z-fiducials (Comeau et al 2000, Lindseth et al 2003, Chen et al 2006, each of which consists of three thin wires stretched between the walls of the phantom forming a Z on the axial plane. When imaged from above, the three wires appear as three bright points in the ultrasound image and by measuring the relative distance between these points, the line of intersection between the ultrasound plane and the wires can be determined. With this approach, a large number of independent datapoints can be collected in one recording. However, again the appearance of the wire in the ultrasound image varies a lot between probes. In addition, the geometry of the Z-fiducials must be adapted to the size of the image plane. 2D shape alignment phantoms: These phantoms include a membrane with an irregular, jagged edge (Sato et al 1998, Langø 2000. The ultrasound plane is aligned with the membrane so that this edge appears as a jagged line in the ultrasound image and the corners of this line can then be located in the image. As with the multiple point targets phantoms, the proper alignment of the ultrasound plane with the membrane is a challenge. Wall phantoms: Here, one or more plane surfaces, such as walls (Prager et al 1998) or membranes (Langø 2000, Baumann et al 2006, are imaged, producing bright lines in the ultrasound image. The most basic versions just use the bottom of a water tank as an imaging target. The lines are easily identified in the ultrasound images and the technique thus lends itself to automatic image processing. However, the visibility of the imaged surface in the images is very dependent on the angle between the probe and the surface. Reverberation artefacts caused by multiple reflections of the sound is also a challenge. The references given here are only meant as examples. More exhaustive references and even more phantom variants can be found in Mercier et al (2005) and Hsu et al (2009).

Materials and methods
Our method uses a point target phantom with a plastic sphere in a water tank as an imaging target. This is similar to the approach taken by Sauer et al (2001). The setup includes an optical tracking system (Polaris Spectra, NDI, Waterloo, ON, Canada) that measures the positions of small, retroreflective plastic spheres with a diameter of 11.5 mm. We designed a calibration arm made mainly from plastic incorporating seven such spheres: one sphere was mounted at the tip of the arm to function as an imaging target, while the other six were mounted in a particular pattern on a plastic plate at the other end of the arm to function as a tracking reference. The arm is shown in figure 2(a). The position p r of the centre of the imaging target relative to the reference was measured using the tracking system. The arm was then attached to a six-axis robot arm (UR5, Universal Robots, Odense, Denmark) mounted next to a water tank on a bench and positioned so that the end of the arm with the imaging target reached into the water. The robot had a repeatability of ± 0.1 mm. The ultrasound probe that was to be calibrated, which was also equipped with a position sensor, was positioned directly above the imaging target and rigidly attached to the bench. Finally, a computer running an in-house navigation system (CustusX, SINTEF, Trondheim, Norway) (Askeland et al 2011) was connected both to the ultrasound scanner and to the tracking system. The complete setup is shown in figure 2(c).
The robot now moved the imaging target to a starting position just outside the ultrasound image plane. It then moved it slowly, at a speed of 1 mm s −1 , first straight through the image plane and then in the opposite direction until it was back at the starting position again. This is illustrated in figure 3. While the target was moving, the ultrasound images produced by the  ultrasound scanner were recorded by the navigation system. For each image k, the system also recorded the position ← M r s k of the ultrasound probe relative to the reference sensor on the calibration arm. The imaging target was then moved to a new starting position. This process was repeated nine times so that the imaging target was passed through the ultrasound image plane a total of 18 times at nine different positions. The distances between the positions were chosen so that they were evenly distributed throughout the plane.
This process was repeated nine times so that the imaging target was passed through the ultrasound image plane a total of 18 times at nine different positions evenly distributed throughout the plane.
The next step was to find the position of the imaging target within each of these 18 recordings. In practice, this meant that for recording j we needed to find the image k j that passed through the centre of the sphere and then to determine the image coordinates p i kj of the sphere centre within this image. This is again illustrated in figure 3. To achieve this, the recorded images were read into the software MATLAB (Mathworks, Natick, MA, USA) and displayed. As shown in figure 4, the images changed as the imaging target moved through the image plane. However, since the target was spherical and moved with a constant speed and direction, the appearance of the images was symmetric around the centre of the sphere. By drawing a rectangular box around the sphere, summing the intensities of all the pixels within the box for   each ultrasound image and plotting this sum against the image number, the result was therefore a graph that was also symmetric. This can be seen in the bottom panel of figure 4. Using this graph for support, it was easy to flip through the displayed images and find the one going through the centre of the imaging target. The resulting image k j only showed the circular surface of the plastic sphere, but knowing its diameter and the pixel size of the ultrasound image, a circle of the same size as the sphere was drawn in the image. By zooming in, this circle could be moved manually so that its circumference corresponded with the surface of the imaged sphere and its centre thus corresponded with the centre of the sphere, as shown in figure 5. This was repeated for each of the 18 recordings of the imaging target. Now, the image k j showing the centre of the imaging target and the image coordinates p i kj of the centre within this image had been found for each of the 18 recordings, i.e. for j = 1, …, 18. The transformations ← M r s kj representing the position and orientation of the ultrasound probe corresponding to each of these images were then extracted from the navigation system.
The positions ← M r s kj of the ultrasound probe corresponding to each of these images were then extracted from the navigation system. Combining this information with the previously measured coordinates p r of the imaging target relative to the reference on the calibration arm, its position relative to the position sensor on the ultrasound probe was found as: The result was two sets of coordinates for each of the 18 selected images, both of which described the position of the imaging target: one in the image coordinate system and one in the probe's coordinate system. The calibration matrix was then given as: and this was calculated using a closed-form method given by Horn (1987).

Adaption to other tracking systems
Our navigation system interfaces not only with the optical tracking system, but also with an electromagnetic tracking system (Aurora, NDI, Waterloo, ON, Canada). This is used in settings where either it is hard to achieve a clear line of sight, or the optical position sensors are too bulky to be integrated properly with the instruments in question. This is typically the case for small or flexible instruments. To enable the calibration of ultrasound probes equipped with such electromagnetic position sensors, we mounted an electromagnetic reference sensor on the calibration arm close to the imaging target. This can be seen in figure 2(b). Since this tracking system could not measure the position of the target directly, we instead found the spatial relationship between the optical and the electromagnetic reference sensors. This was done by rigidly attaching the calibration arm to a cubic plastic box with each side measuring 25 cm. On each of the four side walls there were drilled four small holes 16 cm apart. The position of each of these 16 holes was then measured two times: first using an optically tracked pointer and then using an electromagnetically tracked pointer. Finally, the rigid transformation minimizing the distance between these two point sets was found, again using the method by Horn (1987). The position of the imaging target relative to the electromagnetic reference sensor was then found by applying the resulting transformation to the optically measured position.

Evaluation
To evaluate the calibration method, we used the ultrasound scanner SonixMDP (Analogic, Boston, MA, USA) and chose three ultrasound probes that differed a lot with respect to both shape, size and image resolution. The probes are listed in table 1 and shown in figure 6. One of them was equipped with an optical position sensor and the other two were equipped with electromagnetic position sensors. Each probe was calibrated five times, producing altogether 15 different calibrations. Since each calibration was based on 18 recordings of the imaging target, a total of 270 recordings were made.
To provide a measure of the overall accuracy that was independent of the calibration setup, a precisely engineered accuracy phantom was also used. This phantom consisted of a water tank with two nylon wires crossing each other at its centre. The tank was equipped with a tracking reference and the position p r of the wire cross relative to the reference was measured using a mechanical stylus with an accuracy of 0.01 mm. The tank is shown in figure 7(a). The wire cross was imaged four times for each calibration, resulting in 60 accuracy recordings altogether.
Based on these calibrations and recordings, we calculated four different quality measures (Lindseth et al 2003): the leave-one-out cross-validation error (LooCvE), the point reconstruction accuracy (PRAc), the calibration reproducibility (CR) and the three-dimensional navigation accuracy (3D NAc).
As described previously, each calibration was based on 18 recordings of the imaging target, from which the position of the imaging target both in the image coordinate system and in the probe's coordinate system was found. To determine the LooCvE, we left out the jth recording and calculated a new calibration ← ¬ M s i j based on the remaining 17. The resulting calibration was applied to the image coordinates p i kj of the sphere extracted from the chosen recording, transforming them to the probe's coordinate system. A partial error was then calculated as the euclidean distance between these transformed coordinates and the coordinates p s kj of the sphere measured by the tracking system. This process was repeated for all 18 recordings and the LooCvE was calculated as the average taken over the resulting 18 partial errors, i.e. as: The PRAc is similar to the LooCvE. The difference is that rather than using the same recordings for evaluation that are used to calculate the calibration, we used a separate set of recordings; since we had produced a total of five calibrations for each probe, we used the 72 recordings originally used to produce the other four calibrations. The given calibration was applied to the image coordinates of the sphere extracted from each of these 72 recordings and again the partial error was calculated as the euclidean distance between these transformed coordinates and the sphere coordinates measured by the tracking system. The PRAc for the mth calibration was then found as the average of these 72 partial errors, i.e. as:  The CR is a measure of the precision of the calibration method. Each of the five calibrations were applied to a virtual image point p i virtual , which in this case was chosen to be the lower right-hand corner of the recorded images in accordance with Lindseth et al (2003). This resulted in five points in the probe's coordinate system. The CR for the mth calibration was then calculated as the mean of the euclidean distances from the point transformed by this calibration to the points transformed by the other four calibrations, i.e. as: The last quality measure, the 3D NAc, is based on the ultrasound recordings of the accuracy phantom and it is thus the only one that is independent of the calibration setup. For each recording, the navigation system created a three-dimensional ultrasound volume based on the given probe calibration, the recorded images and the corresponding tracking data. This was done using the reconstruction algorithm Pixel Nearest Neighbour . The two wires were then segmented from the ultrasound volume using a fast,  automatic method for centre line extraction (Smistad 2012). Finally, these centre lines were registered to a model of the wire cross that was based on the mechanical measurements of the wires. This was done using the Iterative Closest Point algorithm (Besl and McKay 1992). The translational part t r of the resulting registration was used as a measure of the distance between the wire cross in the ultrasound volume and the mechanically measured wire cross. The 3D NAc was then found as the mean of this distance taken over all four volumes, i.e. as: where t r l is the translational part of the registration matrix based on the lth recorded ultrasound volume. Figure 7(b) shows both the mechanically measured wire cross and the ultrasound volume for one of the 60 accuracy recordings.
The 3D NAc is in fact a measure of the overall accuracy of the navigation system when navigating on a reconstructed three-dimensional ultrasound volume. Thus, it includes multiple error sources in addition to the probe calibration, such as sensor attachment repeatability, position sensor tracking, synchronization between position data and images and reconstruction algorithm (Lindseth et al 2002). This should therefore be regarded as an upper bound on the accuracy of the probe calibration.

Results
A total of 15 different calibrations were performed. Each calibration took approximately 60 min, out of which setting up the equipment, acquiring data and processing data took around 20 min each. The time spent on setting up the equipment is of course reduced when multiple calibrations are performed at the same time.
The quality measures for the three different probes that were calibrated are shown in table 2. Both the accuracy and precision are good with PRAc below 1.07 mm and CR below 0.89 mm for all 15 calibrations. These results are further supported by the 3D NAc, which shows that the overall accuracy of the system is below 1.45 mm for all the probes and calibrations and as low as 1.09 mm for the C5-2 probe with optical tracking.

Discussion
The motivation for this work was the need for a calibration method which could be applied to any kind of ultrasound probe, regardless of shape, field of view, resolution, contrast or noise level, without any need for adaptation. This is an important requirement in a research setting where new and specialized probes are tested on a regular basis. As described in the introduction, a multitude of calibration methods and phantoms have already been presented, but after having tried a number of them during the last 15 years, we have still not found one that fulfills this requirement. The method presented here is therefore not very novel, but it is a practical solution adapted to our needs. We chose a point-based method mainly due to its simplicity: a point target is easy to get at with the ultrasound probe regardless of its size and shape and since it is spherically symmetric it can be imaged from any angle. This means that the angle of the ultrasound image plane does not have to be aligned with the calibration arm, which greatly simplifies the setup. We found that the imaging target should be relatively large compared to the resolution of the scanner. While small bead-like targets tend to be smeared out in the ultrasound images, the larger plastic spheres appear as a well-defined, semicircular shape, which is easily segmented from the images. The accuracy of the manual segmentation that we propose is hard to determine, as the ground truth is not known. However, the semicircular shape is usually only between one and two millimetres thick and in these cases the segmentation can probably be assumed to have submillimetre accuracy.
The use of a robot for moving the imaging target has both advantages and disadvantages. It facilitates the automation of the calibration process and it performs the movements very accurately and repeatably. Slow, steady motions of the target produces a better image quality, which again makes the subsequent segmentation easier and more accurate. The main advantage, however, is that it solves the problem of aligning the centre of the target with the ultrasound plane, which is one of the major problems of point-based methods (Hsu 2007). Since the motion of the imaging target is performed at a constant speed and follows a linear trajectory, the acquired ultrasound images are completely symmetric around the centre of the target. It is therefore straightforward to identify the image corresponding to this centre. Moreover, due to the low speed of the motion, the error made by missing the centre with a few frames is negligible.
One problem when operating with a research system in combination with a large number of different ultrasound scanners is the time lag between the images collected from the scanner and the position data collected from the tracking system. If the lag is constant, which it often is, this can be measured and compensated for. However, depending on the hardware being used, the lag may vary, making accurate compensation difficult. The calibration method presented here goes a long way towards eliminating this problem during calibration. This is because for each position where the imaging target was passed through the ultrasound image plane, it was moved both back and forth with a constant speed along the same linear trajectory. Assuming the time lag of the system remained constant during this movement, which lasted only 40 s, the position errors introduced by the lag would be exactly the opposite for the two movements and thus cancel each other out.
It is challenging to measure the exact accuracy of a probe calibration. The LooCvE and the PRAc are both based on data from the same phantom and setup that are used for the calibration itself. They can therefore provide information about the consistency of the collected data, but systematic errors, e.g. caused by an inaccurate characterization of the calibration phantom (in our case this is the calibration arm), will not be detected. The 3D NAc, on the other hand, is measured using a separate phantom and it will therefore also reveal systematic errors in the calibration method. However, as previously mentioned, this is a measure of the overall accuracy of the navigation system, which includes several other error sources. A high 3D-NAc (i.e. high value, not high accuracy) does therefore not necessarily mean that the calibration accuracy is poor, but a low 3D-NAc (i.e. low number, not low accuracy) means that the the calibration accuracy is good.
We achieved a 3D NAc below 1.45 mm for all the probes, which is sufficient for most clinical uses. This also means that the calibration accuracy is sufficient. The results for PRAc were similar to those presented by Lindseth et al (2003) (see their table 6) and slightly better than those presented by Hsu et al (2009) (see their table 1). Lindseth et al reported considerably higher maximum values, but this may be due to the fact that their results were based on 15 calibrations for each probe, which is three times as many as ours. The CR numbers were similar for all three studies. It is, however, important to note that the variation between different probes and different tracking systems is often larger than the variation between different calibration methods when it comes to accuracy. While both of the cited studies used fairly standard probes and optical tracking systems, two of the probes used in this study (the L13-7 and the LAP9-4) had a shape which would make them difficult to calibrate with most other calibration methods. For these probes we also used an electromagnetic tracking system which is much more vulnerable to disturbances from the surroundings than the optical systems. The fact that we also achieved comparable results with these probes is a testimony to its high performance.
The main disadvantage of this method is that the robot is both expensive and spacedemanding. However, industrial robots are becoming both smaller, cheaper and more easily programmable, with prices starting around €15 000, or about half of that of a traditional industrial robot. This may still be somewhat expensive if probe calibration is the only task that is to be performed, but a robotic arm like the UR5 is a very flexible tool which can be used to automize a wide range of laboratory tasks. Together, it may justify the investment. One should also note that the robot may be replaced by any mechanical device capable of providing a slow, linear movement at a relatively constant speed. This would, however, reduce the flexibility of the method and also the potential of automation.
Another disadvantage is that both data acquisition and data processing are relatively time consuming compared to other methods. However, we found that our method provides a reasonable compromise between time, flexibility and reliability. It also has a great potential for further automation. By integrating the robot with the navigation system, a simple program could perform the entire data acquisition process without any manual interaction. The data processing, on the other hand, is harder to automate as different probes can produce very different images of the same object. In our case, the image of the plastic sphere can vary from a thin line with a clear, semicircular shape, to a fuzzy dot which is somewhat flatter on one side than on the other. Still, tools could be made to aid the segmentation of the sphere from the images, e.g. by automatic symmetry detection and this could further speed up the process.

Conclusion
The proposed probe calibration method can be used to calibrate a wide range of different probes without any adaptation and with high accuracy and repeatability. It is thus especially suitable in a research setting where new and specialized probes are tested on a regular basis. Although the method involves some manual steps today, it has the potential to be made fully automatic.