Calibration of an orientation sensor for freehand 3D ultrasound and its use in a hybrid acquisition system

Background Freehand 3D ultrasound is a powerful imaging modality with many potential applications. However, its reliance on add-on position sensors, which can be expensive, obtrusive and difficult to calibrate, is a major drawback. Alternatively, freehand 3D ultrasound can be acquired without a position sensor using image-based techniques. Sensorless reconstructions exhibit good fine scale detail but are prone to tracking drift, resulting in large scale geometrical distortions. Method We investigate an alternative position sensor, the Xsens MT9-B, which is relatively unobtrusive but measures orientation only. We describe a straightforward approach to calibrating the sensor, and we measure the calibration precision (by repeated calibrations) and the orientation accuracy (using independent orientation measurements). We introduce algorithms that allow the MT9-B potentially to correct both linear and angular drift in sensorless reconstructions. Results The MT9-B can be calibrated to a precision of around 1°. Reconstruction accuracy is also around 1°. The MT9-B was able to eliminate angular drift in sensorless reconstructions, though it had little impact on linear drift. In comparison, six degree-of-freedom drift correction was shown to produce excellent reconstructions. Conclusion Gold standard freehand 3D ultrasound acquisition requires the synthesis of image-based techniques, for good fine scale detail, and position sensors, for good large scale geometrical accuracy. A hybrid system incorporating the MT9-B offers an attractive compromise between quality and ease of use. The position sensor is unobtrusive and the system is capable of faithful acquisition, with the one exception of linear drift in the elevational direction.


Background
3D ultrasound [1,2] is an emerging medical imaging modality with a wide range of potential applications [3]. The data can be acquired using dedicated 3D probes incorporating either a 2D array or an oscillating head which sweeps the B-scan plane over a fixed volume. The alternative, freehand approach involves the clinician manually sweeping a conventional probe over the target: by attaching a position sensor to the probe, each B-scan can be labelled with its position and orientation. The B-scans thus form a 3D data set which can be visualised and processed in a number of ways to extract clinically useful information. The freehand approach offers the advantages of arbitrary acquisition volumes, with translation as well as rotation of the scan head, low cost, and compatibility with the full palette of existing 2D probes. There are also disadvantages, including slow acquisition: freehand acquisition is not suitable for 4D scanning.
The 2D array and oscillating head approaches are the current focus of commercial activity. Perhaps the greatest barrier to more widespread uptake of freehand scanning is the add-on position sensor. Most position sensors for freehand 3D ultrasound fall into two categories: optical sensors, which employ two or more cameras to track targets attached to the probe, and magnetic sensors, which use a small receiver mounted on the probe to measure a spatially varying magnetic field generated by a fixed base station. Both types of sensor require careful calibration [4] and impose constraints on the scanning protocol. For optical position sensors, the user must maintain a clear line of sight between the probe and the cameras, and must be careful not to stray outside the cameras' field of view. Magnetic sensors also suffer from a limited operating region: furthermore, the immediate vicinity must be kept clear of metallic objects and stray magnetic fields. 3D reconstructions based on position sensor readings also suffer from fine scale jitter artefacts. The jitter arises through a combination of noisy sensor readings and misregistration caused by small probe pressure distortions of the anatomy. As such, it is not possible to completely remove the jitter simply by processing the sensor readings.
It is against this background that we evaluate an unconventional sensor that has received no attention from the 3D ultrasound community. The Xsens MT9-B [5] uses MEMS magnetometers, accelerometers and rate gyros to determine its orientation. MEMS devices are extremely compact and could easily be built into a probe assembly without inconveniencing the user. There is no line of sight requirement, no restricted operating region, just a moderate sensitivity to proximate ferromagnetic materials, although this need not be an issue in typical clinical practice with appropriate consideration of the scanning environment. The obvious limitation is that the MT9-B is a three degree-of-freedom device, measuring orientation but not position. How such a device might be useful for freehand 3D ultrasound acquisition will be explained shortly.
Freehand 3D ultrasound can also be acquired, without a position sensor, by deducing the probe's motion from the B-scan images themselves. Consider two neighbouring Bscans A and B in a freehand sequence. Any in-plane motion between A and B (translation in the axial and lat-eral directions, roll around the elevational axis) is readily determined using standard 2D image registration techniques [6,7]. Perhaps surprisingly, the out-of-plane motion components can also be estimated from the images [8][9][10]. This is because the focusing of the ultrasound beam is far from perfect. Consequently, the backscattered signal at any point in a B-scan is a function of the scatterers in a certain resolution cell around that point. The resolution cells are particularly elongated in the elevational direction and there is considerable spatial overlap between cells on neighbouring B-scans -see Figure 1. The echo signals in corresponding patches on A and B are therefore correlated, with the degree of correlation depending on the patches' elevational separation. The correlation between three (non-collinear) pairs of patches can therefore be used to infer the three patch separations and hence the out-of-plane separation, tilt and yaw of A relative to B.
In our own recent work on sensorless freehand 3D ultrasound, we have extended its capabilities to encompass arbitrary tissue types [11] and arbitrary probe motion [12,13]. Compared with position sensor-based alternatives, sensorless reconstructions exhibit superior fine scale Principle of elevational speckle decorrelation Figure 1 Principle of elevational speckle decorrelation. The inplane motion between scans A and B (translation in the x and y directions, roll around the plane normal) is readily determined using conventional 2D image registration techniques. This leaves three degrees of freedom: translation in the elevational direction, tilt (rotation about x) and yaw (rotation about y). Consider corresponding patches in scans A and B (the shaded ellipses). Because of the imperfect elevational focusing, the contents of the patches depend on scatterers within overlapping resolution cells (the hollow ellipsoids) and are therefore correlated. The correlation coefficient depends on the degree of overlap and hence the elevational separation. It follows that, given a suitable decorrelation curve, a measured correlation ρ 1 can be used to look up the corresponding separation d 1 . Repeating this process for three (or more) non-collinear patches determines the out-of-plane separation, tilt and yaw of A relative to B.  [13]. However, it is difficult to eliminate all sources of bias from the elevational offset estimates. Consequently, there is a cumulative drift error as the inter-frame displacements are concatenated to build up the overall 3D reconstruction. It is at this point that we return to our earlier discussion of the MT9-B orientation sensor. Clearly, this device could be used to correct any angular drift in sensorless reconstructions. Less obviously, it has the potential to ameliorate linear drift too, by filtering out incorrect elevational separation estimates that are inconsistent with the measured orientation.
In the following sections, we investigate the MT9-B with this drift-correction application in mind. We describe its calibration and measure the calibration precision and the reconstruction accuracy. Finally, we describe how orientation measurements can be used to reduce drift in sensorless reconstructions.

Calibrating the MT9-B position sensor
The MT9-B was mounted on a 5-10 MHz linear array probe connected to a Dynamic Imaging Diasus [14] ultrasound machine. The depthf setting was 4 cm with a single focus at 2 cm. Analogue RF ultrasound signals were digitised after receive focusing and time-gain compensation, but before log-compression and envelope detection, using a Gage Compuscope CS14200 14-bit digitiser [15]. Sampling was at 66.67 MHz, synchronous with the ultrasound machine's internal clock: this synchronisation minimises phase jitter between vectors. The acquired vectors were filtered with a 5-10 MHz filter, then envelope-detected using the Hilbert transform. Each B-scan comprises 127 vectors, with 3818 samples per vector. The resolution is approximately 0.01 mm per sample in the axial direction and 0.3 mm per vector in the lateral direction.
Knowing the MT9-B's orientation is quite different to knowing the B-scan's orientation. To relate the two, we need to determine the three degree-of-freedom rotation matrix between the B-scan and MT9-B coordinate systems. This is a restricted version of the well-known spatial calibration problem for six degree-of-freedom position sensors. While there are many techniques for solving the six degree-of-freedom problem [16], most are not suitable for our restricted, orientation-only problem, since the rotation parameters are not adequately decoupled from the translation parameters.
One technique which does lend itself to the restricted problem is the plane-based calibration technique of [17]. This involves scanning a flat plane at the bottom of a water bath, producing a straight-line echo in the B-scan. For the full, six degree-of-freedom calibration problem, there are eleven unknowns: six for the desired rigid body transformation between the sensor and B-scan coordinate systems, three for the position and orientation of the scanned plane, and two image scale parameters. The scale parameters are determined separately using the ultrasound scanner's length measurement facility [18]. For the remaining nine unknowns, each image of the plane, coupled with the corresponding reading from the position sensor, provides two constraints. Provided the probe is moved in such a way as to exercise all degrees of freedom [19], a set of images and sensor readings provides sufficient constraints to solve for these position and orientation unknowns. [20] describes a three-stage nonlinear optimisation process for determining the solution. Significantly, the first stage solves for the rotation parameters using only the orientation readings from the position sensor. This is precisely what is required to calibrate the MT9-B sensor.
Since there are only five degrees of freedom at this stage (three for the orientation of the MT9-B with respect to the B-scan, and two for the orientation of the plane), a simpler set of probe motions provides the necessary constraints, as illustrated in Figure 2. With this one exception, we adopt the same calibration protocol as in [19], including the use of a "Cambridge phantom" to produce echoes from a flat plane without sensitivity to beam-width effects.
Given only orientation data, there is one ambiguity in the calibration solution. It is possible to rotate the B-scan 180° in-plane, this "mirror" solution remaining consistent with the MT9-B readings and the orientations of the echo lines in the B-scans. We overcome this problem by making use of the known geometry of the experimental setup. The coordinate system in which the MT9-B provides its readings is not arbitrary. The x direction is aligned with magnetic north and the z direction is upwards, defined by gravity. In addition, the calibration plane is roughly horizontal and viewed from above. The average B-scan axial direction should therefore have a negative z component in MT9-B coordinates. By checking whether this is the case, we can resolve the ambiguity and correct the calibration if necessary.

Assessing calibration precision and reconstruction accuracy
The calibration precision was assessed by repeating the calibration process 25 times and calculating the standard deviations of the three significant components of the solution (the orientation of the B-scan with respect to the MT9-B in terms of roll, yaw and tilt). At each repetition, the probe was removed and remounted in the Cambridge phantom.
Reconstruction accuracy was assessed by mounting the probe and sensor assembly on a mechanism that could be set accurately to known positions and orientations [21].
The mechanism allows rotation about two axes, and was set to eleven tilt values in the range -16.4° to +2.2° from the vertical, and seven yaw values in the range -5.5° to +5.5°, giving 77 orientations in total (the angular ranges were limited by the design of the mechanism). Since the probe was attached to the mechanism in such a manner as to align the B-scan approximately with two of the mechanism's principal axes, the 77 measured orientations were readily converted to B-scan orientations.
These B-scan orientations were also obtained from the calibrated MT9-B readings. The rotation between every pair of frames from the 77 was calculated, using each of the 25 MT9-B calibrations, giving a total of 73,150 inter-frame orientations. Each orientation was expressed in terms of rotations about the frame pair's average x, y and z axes, which we refer to as tilt, yaw and roll respectively. The MT9-B's accuracy was assessed in terms of the mean and standard deviation of the errors in these angles, using the corresponding mechanism readings as ground truth.
This experiment was performed with the probe mounted in two orientations on the mechanism. For the first experiment, the probe was aligned so that tilt and yaw of the mechanism was approximately equivalent to tilt and yaw of the B-scan. The probe was then remounted at approximately 90° to its original orientation, so that tilting the mechanism produced roll in the B-scan. This second mounting also brought the sensor closer to a ferromagnetic component of the mechanism, allowing some assessment of the MT9-B's magnetic sensitivity.

Sensorless drift correction using orientation measurements
While the MT9-B is unobtrusive, it measures only orientation and is therefore not suitable for use in a position sensor-based freehand 3D ultrasound system. However, it could be used to reduce drift errors in sensorless acquisition systems. Figure 3 illustrates the two predominant drift errors, namely stretch and tilt [13]. Speckle decorrelation occurs for reasons other than elevational separation, including probe rotation, physiological motion, probe pressure-induced deformation and electrical noise. The reduced correlation between B-scan patches is misinterpreted as increased elevational separation, leading to an erroneously stretched reconstruction in the elevational direction. Furthermore, there is a depth-dependent component of the bias, resulting in angular drift in the tilt direction, as illustrated in Figure 3.
In describing how the MT9-B might be used to ameliorate these drift errors, and assessing its efficacy in this role, it is informative to begin with six degree-of-freedom drift correction using a six degree-of-freedom position sensor. In the sensorless reconstruction, the position of each frame f is expressed as a homogeneous transformation matrix S f . We correct these positions by pre-multiplying S f by T f , where The correction C f , relative to the previous frame in the sequence, is most naturally expressed in B-scan coordinates (since the various sources of bias depend on B-scan depth). Multiplying by T f-1 takes account of all the adjustments made to earlier frames in the sequence.

T T S C S
S C S  where t is a value between 0 and 1, defining how far through the sequence the frame is. The cubic formulation ensures that the drift correction is smoothly varying along the sequence, and cannot over-fit to noise in the position sensor readings. We also allow an initial offset at the first frame, defined by an additional six parameters, independent of the drift correction in (3). The overall effect can be thought of as a rigid registration of the sensorless sequence to the position sensor readings (the additional six parameters), combined with the non-rigid deformation in (3).
The coefficients of Equation (3)   ρ decrease the angle between each frame pair, and this alone would move the end frame significantly closer to its correct position.
In addition, there is another way in which orientation measurements might be used to ameliorate linear, stretch drift. The orientation data could support a filter to remove biased elevational distance estimates. Figure 4 shows a typical frame pair and the resulting pattern of distance estimates in the axial direction, as determined by speckle decorrelation. In this example, the bias in the distance estimates is larger towards the bottom of the frame. This could be caused by a lack of signal in this region of the image, producing lower correlations and hence overestimated distance estimates. Fitting a B-scan plane (using least squares) to these distances results in a tilt bias and a small overestimate in the elevational centre offset, as shown in Figure 4(a).
The proposed filter is based on the Random Sample Consensus algorithm [23]. We consider each elevational distance estimate individually and fit a B-scan plane to it using the orientation measured by the MT9-B. We then count how many of the other distance estimates support the hypothesis that this plane is correct. The plane that attracts the largest support provides the set of "correct" distance estimates and any others are then ignored. In the example in Figure 4(b), the distance estimates highlighted would be ignored, resulting in the much better fit shown in the figure. The filter has one tunable parameter: the threshold specifying how close a distance estimate must be to the fitted plane to count as a "supporter". This is set to one standard deviation of the estimated uncertainty in the distances, determined using the method described in [13].
Drift correction was assessed by freehand scanning of a speckle phantom (from the Department of Medical Physics, University of Wisconsin, Madison, Wisconsin, USA) and animal tissue in vitro. The sensorless reconstructions were obtained using the technique described in [13] and ground truth was available from a Northern Digital Polaris [24] optical position sensor tracking a Traxtal AdapTrax [25] target attached to the probe. Temporal calibration was performed according to the technique in [19]. For reasons of comparison and algorithm verification, the sensorless reconstructions were first corrected using both position and orientation data from the Polaris, then using just the orientation data, in the manner described above. Note that the MT9-B was not used in these experiments, with good reason: the Polaris provides an essential ground truth against which to assess the orientation-only drift correction. However, we also present a qualitative in vivo scan of a human calf muscle corrected using the MT9-B. Table 1 shows precision statistics for the MT9-B calibration process. The in-plane roll, and the out-of-plane yaw and tilt angles, are expressed in the B-scan reference frame.

MT9-B calibration precision
The results indicate significant dependence on the accuracy of detecting the lines (the echoes of the calibration plane) in the B-scan images. For example, the tilt calibra-   Figure 2. Unfortunately, when the plane is imaged in this manner it is difficult to obtain a clear echo and measure its orientation reliably. It is therefore not surprising that the worst precision is in the tilt angle. Nevertheless, the calibration precision is comparable to previous results for six degree-offreedom sensors [19]. Calibration imprecision is one factor affecting the reconstruction accuracy determined in the next section. Table 2 summarises the results of the two accuracy experiments. The tabulated statistics relate to the accuracy of measuring angles between frame pairs and therefore give an indication of the accuracy for drift correction. It is immediately apparent from these results that the yaw accuracy is worse than the roll and tilt accuracy, and is particularly bad in the roll-yaw experiment. This difference is not unexpected, since the roll and tilt measurements depend predominantly on the MT9-B's accelerometer measurements of the gravitational vertical, while the yaw relies more on the magnetometer. In the roll-yaw experiment, the sensor was mounted closer to the field-distorting ferromagnetic material, causing larger errors in the yaw angle. Figure 5 shows some representative examples of the orientations measured using the MT9-B, compared to the actual orientations set by the mechanism. For each of the three rotation axes, the graphs are for a sequence of frames selected from the grid of orientations so that there is no rotation about the other two axes. The most notable feature of these graphs is the yaw drift in the first column of the figure. This is due to the magnetic field distortion caused by the small ferromagnetic parts of the controlling mechanism. The better yaw result in the second column is from the tilt-yaw experiment, where the sensor was mounted further away from the ferromagnetic part. In comparison, the roll and tilt graphs, which do not depend significantly on the MT9-B's magnetometer readings, show no sign of drift.

MT9-B reconstruction accuracy
These results demonstrate that the MT9-B is capable of high accuracy orientation measurement, provided it is kept clear of magnetic field distortions. While this requires careful consideration of the scanning environment, it is not difficult to achieve. For the small amounts of ferromagnetic material present in our case, a few centimetres clearance was sufficient to remove the yaw drift and obtain an accuracy within 1°. Figure 6 shows representative results for ten sweeps of in vitro animal tissue (a joint of beef) and an additional sweep of a speckle phantom with spherical inclusions. The reconstructions are shown in terms of their accumulated elevational length (left column) and their accumulated tilt (right column). Tables 3 and 4 show numerical results for all eleven data sets. Note the intentional nonmonotonic scanning motion in the second beef graph and the phantom data set. Data sets 1-10 were essentially linear sweeps in the elevational direction, with data sets 6-10 including a deliberate change of direction part way through the sweep. The phantom data set also featured deliberate probe rotation. Each plot in Figure 6 comprises four traces: the Polaris reconstruction, the sensorless reconstruction, and the drift-corrected sensorless reconstructions using six degree-of-freedom and orientationonly correction.

Sensorless drift correction
It is evident from the figure and tables that the six degreeof-freedom corrections achieve the goal of producing a geometrically accurate reconstruction. The orientationonly correction achieves the expected elimination of angular drift, but the filtering strategy has had little effect on the linear drift. From this, we can infer that any axial variation in elevational offset bias (this variation being the cause of tilt drift) is evenly distributed about the centre of the B-scan, so correcting the orientation does not affect the centre offset. Consequently, a small bias will remain in any measurements taken from the reconstructed data. Figures 7 and 8 show reslices though beef data set 1 and the phantom data set respectively, for the various recon-    beef data set 1 beef data set 7 phantom data set struction and drift correction methods. Although the sensorless reconstructions exhibit no qualitatively obvious length error, there is a significant tilt error which is particularly obvious in the longer phantom sequence. Also, the fine scale accuracy of the sensorless reconstruction is noticeably superior to that of the position sensor reconstruction, particularly in the beef data, which is more susceptible to probe pressure jitter. Both of the drift correction methods are able to significantly reduce the tilt error, without over-fitting to position sensor noise and thereby sacrificing the fine scale accuracy. It is evident from these reslices that tilt bias is a more significant problem than length bias. The limited orientation correction therefore provides a worthwhile improvement, despite the fact that it is unable to correct length bias. Figure 9 shows reslices through the in vivo calf muscle data set, corrected using the MT9-B. In this case, there are no six degree-of-freedom sensor readings for comparison. However, the data was recorded using an approximately linear   scanning protocol, so we can be fairly confident that the drift corrected result using the MT9-B is a more accurate representation of the scanning subject than that obtained using purely sensorless techniques.

Conclusion
Gold standard freehand 3D ultrasound acquisition, with accurate geometrical dimensions along with good fine scale detail, requires image-based reconstruction with six degree-of-freedom drift correction using a position sensor. Suitable position sensors, however, impose restrictions on the scanning protocol that make them somewhat obtrusive. A compromise is offered by the MT9-B orientation sensor. It is relatively unobtrusive and capable of correcting angular drift in sensorless reconstructions to within around 1°, although linear drift remains a problem. The MT9-B is easily calibrated using a variation of the established plane-based calibration method.