Passive Markers for Tracking Surgical Instruments in Real-Time 3-D Ultrasound Imaging

A family of passive echogenic markers is presented by which the position and orientation of a surgical instrument can be determined in a 3-D ultrasound volume, using simple image processing. Markers are attached near the distal end of the instrument so that they appear in the ultrasound volume along with the instrument tip. They are detected and measured within the ultrasound image, thus requiring no external tracking device. This approach facilitates imaging instruments and tissue simultaneously in ultrasound-guided interventions. Marker-based estimates of instrument pose can be used in augmented reality displays or for image-based servoing. Design principles for marker shapes are presented that ensure imaging system and measurement uniqueness constraints are met. An error analysis is included that can be used to guide marker design and which also establishes a lower bound on measurement uncertainty. Finally, examples of marker measurement and tracking algorithms are presented along with experimental validation of the concepts.


I. Introduction
ULTRASOUND imaging is a useful modality for guiding certain types of minimally invasive interventions due to its portability and safety.Given the larger field of view now provided by real-time 3-D imaging over the traditional 2-D imaging, more complex ultrasound-guided procedures may be possible [1].Of particular interest are procedures in which surgical instruments and tissue structures are imaged simultaneously to achieve precise manipulations.
In cardiac surgery, real-time 3-D imaging also enables off-pump beating-heart procedures, which eliminate the significant morbidity associated with cardio-pulmonary bypass, including delay of neural development in children, decline in adult cognitive performance, mechanical damage from vessel cannulation and stroke [2]- [4].In one example, real-time 3-D ultrasound imaging is being investigated for guiding repairs of atrial septal defects inside the beating-heart [5].In this case, surgeons insert a long, thin anchor driver through the Published approaches for addressing the instrument imaging problem include instrument modification, image processing techniques, and active tracking sensors.Instrument modification involves the application of coatings or surface modifications to reduce the specularity of reflections or to increase absorption [10], [11] [12].Image processing methods apply search techniques to locate an instrument in an image.Examples include use of the Hough transform to detect biopsy needles in 2-D ultrasound images [13], use of principal component analysis to also detect biopsy needles [14], and use of the radon transform to detect the axis of an instrument shaft in 3-D ultrasound volumes [15].
Active tracking uses sensors to detect instrument position relative to the ultrasound image.Methods, such as electromagnetic and optical tracking which track objects with respect to an inertial frame, require sensors for both the instrument and the ultrasound transducer to locate the instrument within the image.In comparison, they offer a trade-off between accuracy and accessibility.Electromagnetic tracking can track sensors within the body, but suffers from environmental noise and distortion due to metal objects in and near the tracking field.The latest generation of sensors allows needles to be tracked at the tip [16].Optical tracking provides highly accurate tracking data, but requires a line of sight between sensors and the tracking camera.As a result, all position tracking must be done outside the body.In this case, orientation tracking errors cause instrument tip position errors in proportion to the distance of the tip from the sensor.In addition to system measurement errors, tracking uncertainty arises from calibration errors.Various metrics on probe tracking calibration are used in the literature, including reconstruction accuracy, precision and navigation accuracy [17].Calibration accuracy of 0.6 mm has been achieved with optical tracking of the ultrasound transducer [18] in a phantom study, and 0.90-1.44mm with electromagnetic tracking [16] in a water-tank study.
Another active tracking method involves using an active ultrasound receiver mounted on a cardiac catheter to detect the ultrasound beam strength and direction to determine the catheter's position within the image volume directly [19].Mean accuracy between 0.22 ± 0.11 mm and 0.47 ± 0.47 mm in in vitro experiments, depending on the distance between the catheter and the ultrasound transducer, was achieved using this method.Multiple sensors would be required to determine instrument orientation, however.
In summary, these approaches all possess limitations for image-guided procedures such as intracardiac surgery.While optical tracking can be considered the gold standard for tracking outside the body, its line of sight constraint reduces its accuracy for applications inside the body.Electromagnetic tracking is suitable for intra-body tracking, but is subject to field distortion errors arising from metallic tools moving within the surgical field.Active ultrasound receivers yield acceptable accuracy but may be subject to artifact error when placed in a cardiac chamber with other instruments acting as acoustic reflectors.
This paper presents a novel method for tracking instruments in 3-D ultrasound images using passive instrument markers, which can be easily added to existing surgical instruments.Markers are constructed to possess two properties: 1) they appear clearly when imaged along with tissue regardless of instrument appearance, and 2) marker position and orientation can be determined from a single image using simple image processing.This approach does not require tracking of the ultrasound transducer.
Markers are similar to devices known as stereotactic frames, which have been studied extensively for imaging modalities such as CT and MRI [20], [21].A stereotactic frame consists of a shape that appears uniquely when imaged at various positions and orientations and is constructed of material easily seen in a particular imaging modality.Implantable markers have been utilized in many ultrasound based treatments, such as intraoperative ultrasound to preoperative MRI registration in neurosurgery [22] and prostate localization in abdominal radiotherapy [23].Typically, these fiducials are point markers used for three degree-of-freedom tissue localization by point-cloud registration.Thus, they are not appropriate for accurate six degree-of-freedom instrument tracking.One exception is the fluoroscopic fiducial marker of Jain et al. [24] for use in orthopedic surgery.In similarity with the markers proposed here, Jain's marker consists of a mathematically optimized set of geometric entities.
The paper is arranged as follows.The next section gives an overview of the instrument markers including the design constraints defining the family of possible marker shapes and measurement error analysis.Section III describes marker design techniques for minimizing measurement errors.Section IV then presents experimental validation of marker-based instrument tracking in real-time 3-D ultrasound.Section V presents a discussion of experimental results and illustrates the clinical value of marker-based tracking.

II. Passive Instrument Markers
As depicted in Fig. 1, markers consist of two parts, a cylindrical sleeve that can be fit over the shaft of a surgical instrument and ridges of constant height and width fixed to the outer surface of the sleeve.The ridges trace out prescribed paths on the sleeve's surface, which, when imaged, indicate the marker's position along and orientation about the cylinder's axis.The cylindrical shape allows the markers to have a low profile on the instrument shaft to fit through access ports used in minimally invasive surgery.The family of markers is characterized by a variable number of ridges and a variety of ridge paths; these are together referred to as the marker shape.It is assumed here that instruments always possess a straight cylindrical shaft, to which a marker can be attached.The instrument shaft is often visible in the ultrasound image and is used in conjunction with the marker shape.
As shown in Fig. 2, instrument tracking consists of determining the transformation , which relates the instrument's body coordinate frame to the ultrasound image coordinate frame.This transformation can be decomposed into three elements The transformation relates an intermediate coordinate frame, A, located on the instrument shaft, to the image coordinate frame.Frame A is defined such that z A corresponds to the shaft axis and points toward the instrument tip, while y A lies in the plane defined by z A and the apex of the ultrasound volume, as shown in the inset of Fig. 2.
The second transformation, , relates the marker's body frame M to the intermediate frame A. It is based on the two degrees of freedom provided by the marker.They are denoted θ and t, the rotation about and translation along the instrument shaft from frame A, respectively.
Finally, the constant transformation relates the instrument body frame (not shown) to the marker body frame.This is determined via a number of possible methods when the marker is attached to the instrument.One option is to apply the marker with the use of a jig, which places the marker and instrument in a precisely known relative position.Another option is to precisely measure the marker location with respect to the instrument using calipers or a calibrated camera.
Calculation of and is performed via two image processing steps.The first involves detection and segmentation of the instrument shaft.This returns the shaft axis vector z A as well as the linear equation in image coordinates of the shaft axis.Using this equation, a reference point on the shaft axis, r imA , corresponding to the origin of frame A in image coordinates, can be computed.Note that this point is not being tracked in the image, but instead its coordinates are computed by intersecting the line of the tracked shaft axis with the known constant boundary of the scanned ultrasound volume.
The second step image processing step involves detecting the marker on the instrument shaft and measuring the locations of the ridges along the top surface.The coordinates of the ridges, along z A , are collected into a vector l = [l 1 … l n ] T , where n is the number of marker ridges, as shown in Fig. 1.
As shown in Fig. 3, marker shape is described by a vector-valued function , which is related to the measurement l by (2) Here, t, as shown in Fig. 2, is the distance along z A from frame A to frame M. The components of f(θ) describe the coordinates of the marker ridges along the shaft axis versus rotation angle θ about the shaft axis.The continuity of the components of f(θ) is an important aspect of marker shape.Marker shapes can be categorized into three groups.
In all cases, df n /dθ is defined from above at discontinuities.It is assumed that the components of f(θ) are piecewise smooth functions.
Uncertainty arises in the measurement l from three sources: noise, marker manufacturing tolerance, and finite image resolution.Noise encompasses both electronic noise in the ultrasound transducer and speckle.Four physical aspects of ultrasound lead to finite resolution of objects: frequency, finite scan line spacing, transmit/receive beam-forming, and filtering after image formation.Uncertainty in l is also coupled to uncertainty in the shaft axis, z A in Fig. 2. Note that this error can be different for each element of l.
Uncertainty in the shaft axis is inversely proportional to both the amount of the shaft within the ultrasound image and the shaft's echogenicity.Echogenicity can be increased by extending the marker sleeve over the entire length of the shaft.Widening the field of view can also increase the amount of the shaft within the image, but may reduce resolution of the marker ridges.Since it is assumed that the instrument shaft is straight, measurement error can also arise if the shaft is bent or flexes due to mechanical loading.
A full treatment of material selection for instrument markers to achieve maximum visibility and minimum measurement error as a function of probe type and frequency and signal processing settings has not been performed and is beyond the scope of this paper.
Appropriate materials for the marker sleeve, however, should in general have much lower acoustic impedance than metal to avoid specular reflection and side-lobe artifacts and should reflect sound from a wide range of incident angles.Appropriate materials for the marker ridges are more variable and are likely to be more sensitive to probe type and imaging system settings.Large ridges may be constructed of soft materials, which scatter evenly, to increase their detail, while small ridges may be constructed of hard materials, which reflect strongly, to increase their contrast.
Fig. 4 shows an example marker shape, through which the constraints described in the next subsection can be illustrated.The shape is described by (3) where α = 1/π, β = 2, γ = 2, and τ = 8.As shown in the figure, marker length is 10 mm and the minimum spacing is 2 mm.

A. Constraints on Marker Shape
Marker shape can be designed to satisfy the requirements of particular surgical interventions.Factors including instrument range of motion and measurement uncertainty in position and orientation are affected by the marker shape.Measurement uncertainty is affected by three qualities of marker shape: visibility of marker ridges, the choice for f(θ), and the uniqueness of solutions for θ and t.These requirements are represented in the following seven constraints governing the space of allowable marker shapes.In the following, it is assumed for simplicity that uncertainties in the components of l are independent and identically distributed, as N(0, δl).

1)
Constraint 1: The marker length h, shown in Fig. 4, must be small enough such that ridges are always contained within the image field-of-view (4)

2)
Constraint 2: Marker ridges must be distinguishable.Finite ultrasound resolution leads to a widening of the marker ridges.This can cause adjacent ridges to blend together if they are too close to each other.The minimum ridge spacing, s, shown in Fig. 4, must be large enough such that all ridges are distinguishable.This is defined as the shortest distance between adjacent ridges along the shaft axis (5) The minimum ridge spacing is determined physically by many factors including the axial, lateral, and elevation resolution of the ultrasound transducer, the marker ridge size and material, and the system transmit power.An explicit definition of a function for s based on all relevant parameters is, however, outside the scope of this paper.In the following sections, s is considered an independent variable for simplicity.

3)
Constraint 3: Each measurement l must correspond to a unique combination of θ and t.This constraint is expressed by ( 6) Combining ( 2) and ( 6) gives (7) which expresses the constraint in terms of f(θ).By ( 2), a marker pattern must possess at least two ridges (n ≥ 2) to admit a unique solution for θ and t.By ( 7), the curves describing these two ridges must differ.Equation ( 7) also implies (8) For marker shapes satisfying ( 7) and ( 8), solutions for θ and t can be found by the following procedure.First note that t can be expressed explicitly in terms of θ by (9) It follows that θ is given by (10) Note that if u T f(θ) = 0, then t can be solved independently.

4)
Constraint 4: Given bounded uncertainty in the measurement l, we want to achieve bounded uncertainty in θ.The following is the distance between any two points f(θ) and f(α) in the sense of ( 10) Given that ( 10) is a minimum norm solution, minima of (11) must be bounded from below.Therefore, marker shapes are further constrained such that (12) Since δl is the expected residual for (10), q should at least be greater than twice this uncertainty.It can be increased arbitrarily, however, as a safety margin.Note that this constraint does not apply to two-ridge marker shapes.
To better illustrate these constraints, an n-dimensional Cartesian space is introduced.In this space, is a 1-D curve, called the marker curve.Equations ( 4) and ( 5) define a set of planar constraints, forming a simplex bounded region within which they are satisfied.The constraint simplex in three dimensions is shown in Fig. 5. Two planes are described by (4); they are perpendicular to the f 1 and f 3 axes.The other two planes are defined by ( 5); they are parallel to f 1 = f 2 and f 2 = f 3 but are offset by s along the f 2 and f 3 axes, respectively.Also shown is the example marker curve defined by ( 3) and illustrated in Fig. 4. Note that it lies on a boundary of the constraint simplex.
The coordinates of the vertices, when translated by (0,-s,-2s), can be combined into a matrix (13) where each column represents a vertex.The same result would be obtained by arbitrary combination of h and s, such that h − 2s is constant.Thus, ( 4) and ( 5) together define an effective marker length .Extending this to n dimensions, the effective marker length becomes and maintains the same structure.The matrix A thus provides a convenient way to calculate the impact of ( 4) and ( 5) when n > 3.
Fig. 6 shows both the example marker curve and the vertices of the constraint simplex when projected on a plane orthogonal to the vector u.This projective space illustrates (6), which specifies that the projected path must not intersect.It also illustrates (10) where, given a measurement l, the point of f closest to l in the projective space is the solution for θ.Consequently, to prevent large errors in θ and t, (12) specifies that the projected path must have some minimal separation q.This is equivalent to the projected path having a width of q as outlined in the figure on either side of the projected path.Note that the figure is not shown to scale; to satisfy (11), the projected path (in 2-D) must have a radius of curvature greater than q/2.Projection of the constraint simplex vertices is calculated by pre-multiplying them by a matrix P, the rows of which describe an orthogonal basis for the space of vectors x satisfying 0 = x T u.The projected vertices for an n-ridge marker shape are therefore given by (14) In extending the projected marker curve concept to n dimensions, first note that (12) employs the two-norm, which is equivalent to the path having a spherical cross section of radius q/2.Since both projection on a subspace orthogonal to u and taking a cross section of a path are , the cross section of the projected marker curve has n − 2 dimensions.This can be seen in Fig. 6 where the projection of the example marker curve has a 1-D cross section.
Taking into account that larger diameter shafts can accomo-date higher curvatures, larger slopes, and more discontinuities without distortion, the following three constraints can be defined.For these constraints, parameters would likely be determined empirically.
Constraint 5: The component functions of f(θ) must have bounded curvatures (15) Constraint 6: The component functions of f(θ) must have bounded slopes For markers with more than two ridges, the seven constraints described above admit an infinite number of marker shapes.The design space can be further constrained by determining marker curves that possess minimum measurement uncertainty.This topic is addressed below.where f′ = df/dθ.Then, least squares solutions for δt and δθ are obtained by dropping γ 0 and taking the pseudo-inverse of (19).Finally, substituting ϕ(θ) = cos −1 (f′u/∥f′∥u), uncertainties are given by ( 20)

III. Measurement Uncertainty
The terms ∊ t and ∊ θ are referred to as error factors.Geometrically, f′(θ) is the tangent vector to the marker curve at θ, and ϕ is the angle between f′ and u.Thus, 1/∊ θ (θ) is the magnitude of f(θ) when projected along u.
The error factors for the example shape of Fig. 4 are shown in Fig. 7.The peaks in ∊ θ coincide with rotation angles where f′ is small.

A. Marker Shape Design for Error Minimization
Since the design constraints presented allow a family of marker shapes, individual marker shapes can be selected based on their task appropriateness, i.e., the requirements of an interventional task may indicate appropriate values for the constraints described in Section II-A.For instance, the desired range of motion for an instrument would suggest a maximum allowable marker length h.Imaging depth and probe frequency also indicate appropriate values for q, s, p, d, and e.Furthermore, a grasping task would likely require high accuracy in orientation, while a puncture task would likely require high accuracy in position, indicating a tradeoff between ∊ t and ∊ θ .
A cost function for determining optimal marker shapes is as follows: (22) (23) (24) The infinity norm refers to the maximum value of the cost function over θ and the continuous functions ω t (θ) and ω θ (θ) allow the error factors to be weighted according to their importance for the task.An optimal marker shape minimizes J(f(θ)) while satisfying all design constraints, here taken as independent variables determined by the task.Because ∊ t ∝ 1/∥u∥, the complete problem is solved by first finding locally optimal candidate marker shapes for each allowable number of ridges (n = 2, 3, 4, …) and comparing these to find the minimum overall.
Solutions to (22) are difficult to obtain in closed form due to infinite dimensionality in the space of marker shapes and the sensitivity of solutions to shape constraints.The following subsections address three related subproblems: minimizing ( 23) and ( 24) independently for given projected marker curves and determining the lower bound on (24) for arbitrary projected marker curves.

B. Minimizing J t and J θ
Minimizing J t emphasizes marker shapes where ∥u∥ sin(ϕ(θ)) is maximized.It is achieved by two steps: 1) Maximizing the number of ridges n, thereby maximizing ∥u∥ and 2) selecting marker shapes where ϕ = ±π/2, thereby maximizing sin(ϕ(θ)).These markers have paths which lie entirely in an (n − 1)-dimensional subspace orthogonal to u.
In minimizing J t , design constraints allow different results for discontinuous marker shapes than for semi-continuous and continuous shapes.It is possible to minimize J t and J θ independently with discontinuous marker shapes.For other shapes, however, the cost functions are coupled.For example, in the case where n = 3, Fig. 8 shows both the design constraints and two planes orthogonal to u.The area of any plane within the bounding simplex is less than the area of the projected bounding simplex, as shown in the figure inset.Since semi-continuous and continuous marker shapes have no discontinuities in the range θ ∈ 0 … 2π, their projected marker curves must fit within this smaller area to minimize J t .
Discontinuous patterns, on the other hand, consist of several path segments.The segments can be arranged such that each lies on a different plane orthogonal to u, at varying offsets from the origin.Thus, minimum values for both J t and J θ can always be achieved with discontinuous marker patterns because the minimizations can be decoupled.J t is minimized for arbitrary marker shapes f(θ) with n ridges by conversion to a marker shape satisfying sin(ϕ(θ)) = 1 by (25) where P is the projection matrix shown in (14).The new marker shape, however, is not guaranteed to satisfy all design constraints.Note that (25) also alters ∊ θ .This is illustrated using the example marker shape of Fig. 4. The modified marker shape is shown in Fig. 9(a) and the modified marker curve is shown in Fig. 9(b).While the original curve is aligned with the f 1 = 0 constraint, the new curve lies in a plane orthogonal to u.
Note that the new marker shape, however, has an increased marker length.
Minimization of J θ emphasizes marker shapes where ∥f′(ϕ(θ))∥ is maximized.Note that (26) is the arc length of the projected marker curve.Thus, a marker shape f(θ) must be found with maximum projected path arc length.While there are an infinite number of shapes that satisfy this requirement, J θ is minimized by selecting the shape with appropriate ∊ θ , given ω θ , as follows: (27) This result is shown as follows.Note that minimizing J θ is equivalent to maximizing
Minimization of J θ is illustrated using the example marker of Fig. 4. Given the weighting function ω θ shown in Fig. 10(a), substitution of the θ*(θ) shown in Fig. 10(b) gives the modified marker shape shown in Fig. 11 and a constant value of J θ .Note that the sharp corners in the example modified marker shape are caused by regions of θ*(θ) where the slope is close to zero and ∥d 2 f* / dθ 2 ∥ ≫ ∥d 2 f / dθ 2 ∥.
To achieve the required marker length, a discontinuity is manually introduced in f 3 at θ = π.
The new discontinuous marker shape is of the required length, but it does not meet all design constraints.For instance, the discontinuity reduces q, as described in (12).Thus, although discontinuous shapes can decouple J t and J θ , they must still be designed to meet all constraints.

C. Example Marker Shapes
Four additional marker shape examples are presented to further illustrate the range of possible marker shapes and to illustrate marker design methods.

2)
Example 2: Continuous Shape-Fig.13(a) shows a three-ridge continuous marker shape described by (33) Parameter values are α = 1.10 mm and β = 3.90 mm.The marker length is 10 mm and the minimum spacing is 2 mm.The shape has constant error factors of and ∊ θ = 0.74.The bounding simplex, path, and projected path are shown in Fig. 13(b).The marker curve lies entirely in a plane perpendicular to u, so ∊ t is minimized.

3)
Examples 3 and 4: Discontinuous Shapes-Since cost functions are decoupled for this type of marker shape, error factors can be minimized more easily by working directly with the projected curve.The marker shape can then be obtained from the projected path by suitably embedding it in f-space such that it meets all the design constraints.Marker shapes are derived from the projected paths by (34) where denotes the projected path, denotes the resulting path, and P is the projection matrix (14).Fig. 14(a) shows an example path design for a three-ridge marker using multiple parallel lines spaced at intervals of q = 1.42 mm.Fig. 15(a) shows a similar path design for a fourridge marker.In this case, lines are spaced horizontally and vertically at intervals of q = 0.71 mm.This value of q was chosen to ensure that the projected path length has no more than three discontinuities.The error factors of these designs are close to the minimum possible values.Note that aligning the parallel lines with an edge of the bounding region maximizes the path length.
Example patterns for three-and four-ridge markers based on the parallel line path designs are shown in Fig. 14(b) and Fig. 15(b), respectively.For these designs, marker length is 10 mm and minimum ridge spacing is 2 mm.The error factors for these markers are as follows: and ∊ θ = 0.64 for the three-ridge marker, and and ∊ θ = 0.79 for the fourridge marker.

D. Lower Bound on ∥∊ θ ∥ ∞
Since a projected marker curve has a finite cross-sectional width q, it traces out finite content [(n − 1)-dimensional volume] along its length.Since the marker curve must not intersect itself by ( 12), its length is bounded by the content of the projected constraint simplex.The content of the projected constraint simplex is given by (35) where A p describes the vertices of the projected bounding simplex, as given by ( 6).
Following the derivation given in the appendix, the content swept out by the projected marker curve can be approximated by (36) where Ω is the projected path arc length and Γ is the gamma function.If n = 2, the projected path has no thickness.In this case, .Equating ( 35) and ( 36) and solving for 2πΩ gives an initial lower bound on |∊ θ | ∞ .
It is possible, however, for the projected marker curve to lie along the border of the projected bounding simplex, so some of the content swept out by the projected path lies outside the bounding simplex.To account for this, an approximate correction to the effective marker length, , was determined empirically.In doing so, it was assumed that three-and four-ridge discontinuous marker shapes with projected paths consisting of parallel line segments spaced at 2q have maximum projected path length.
Although not theoretically motivated, this correction increases the content of the projected bounding simplex in proportion to both q and n, which in part determine the content of a projected path.Fig. 16 shows a plot of min(∥∊ θ ∥ ∞ for arbitrary markers with 2, 3, 4, and 5 ridges, versus marker length with s = 2 mm and q = 0.71 mm.

IV. Instrument Tracking Experiments
Two sets of experiments were performed in a water tank to evaluate tracking accuracy using an instrument marker designed for intracardiac applications.In the first set of experiments, the instrument was the only object in the water tank.These experiments provided a best case evaluation of marker performance.To provide a more realistic testing scenario, a second set of experiments was also performed in which a porcine heart was placed in the water tank and the marker was positioned inside the heart.
For all experiments, ultrasound images were collected using a Philips Sonos 7500 ×4 transducer, which generates 3-D image volumes with a 4 MHz working frequency.The system was modified such that scan converted volumes could be obtained in digital form via a 1 gigabit per second Ethernet link at 20 Hz.The transducer was positioned above the water tank, with the crystal array just submerged.Standard image settings were used, 50% overall gain, 0% time gain, high density scan line spacing, and 5 cm image depth.The resulting voxel spacing was (0.301 mm, 0.386 mm, 0.314 mm) in x, y, and z, respectively, as defined in Fig. 2.

A. Tracking Algorithm
Instrument tracking was accomplished via three algorithms: instrument shaft segmentation, marker detection, and Kalman filtering.An Euler angle formulation of instrument marker orientation was used.In this case, the system state vector has 12 elements Instrument shaft segmentation determines the transformation in (1).It is achieved by first randomly subsampling voxels with intensity higher than a given threshold and collecting their coordinates.Principal component analysis (PCA) is then performed on multiple subregions using the subsampled points.To speed processing, the subregions considered can be reduced based on the predicted instrument coordinates.Adjoining regions having high eccentricity and similar alignment correspond to the instrument shaft.Points from these regions are collected and the RANSAC algorithm [25] is applied to remove outliers.Finally, PCA is again used on the collected points to determine the instrument shaft axis.
Marker detection makes use of the instrument shaft axis to determine the transformation in (1).It is achieved by extracting a 2-D image from the image volume in the plane shown in Fig. 2 inset.Fig. 17 shows a sample image with marker ridge locations highlighted.The image is first smoothed by a Gaussian kernel.A surface contour is then extracted by detecting points of high gradient in each column of pixels.A ridge detection kernel is also scanned laterally along the image midline, and its correlation with equally sized subregions of the image is recorded.Peaks in this ridge-likelihood contour are identified as candidate ridge locations.All three-ridge combinations of candidate locations are then collected.Combinations with excessively large or small ridge spacings are eliminated.The rest are weighted based on both their height in the surface contour and their distance from predicted ridge locations from the Kalman filter.The ridge combination with the highest weight is taken to be the instrument marker.

B. Marker Calibration
Due to manufacturing tolerances and small nonlinear image distortions, physical markers can appear slightly differently than the nominal marker designs.This leads to systematic measurement errors in t and θ, which must be corrected via an initial calibration step.
Calibration involves rotating the marker about its axis through 2π several times while continuously recording the ridge locations.This is repeated at a number of positions and shaft orientations within the image volume.For the experiment, pitch angles of 30°, 45°, and 60° were tested at five points within the image volume.Fig. 18 depicts the results from all positions and orientations for the marker of Fig. 17.
Two methods may be applied to compensate for this systematic error.The first consists of fitting curves f*(θ) to the measured ridge locations and subtracting variations in t from them.These curves are then used in place of the nominal marker shape in solving for t and θ.The second method involves fitting curves to the errors in t and θ.These curves are then subtracted from measured values based on the nominal marker shape.Fits were achieved via least squares with sine-cosine basis functions with periods of (π, π/2, π/4 …).The first method was used in this experiment.Since corrections are only defined for 15 image location (x, z) and pitch angle ϕ combinations, interpolation was performed for arbitrary combinations via a distance weighted average (38) where d m = ∥ϕ − ϕ m ∥ + ∥(x, z) − (x, z) m ∥ and c is a matrix of coefficients (α, β) defining f*(θ).

C. Water Tank Experimental Results
In these experiments, the instrument was the only object within the water tank.As shown in Fig. 17, the instrument consisted of a 15-cm-long, 3-mm-diameter stainless steel tube, covered with a thin layer of soft plastic (polyolefin heat shrink insulation).The marker shape was etched into the plastic covering and 0.5-mm-diameter hollow plastic tubes (electrical wire insulation) were glued into the grooves.A continuous marker pattern was used based on example shape 2 in Section III-C, as shown in the figure inset, with h = 10 mm and s = 2 mm.The error factors for this marker shape are and and ∊ θ = 0.74.The figure also shows a cross section of the instrument from an ultrasound volume; the line indicates the estimated instrument shaft axis and the dots indicate the estimated locations of the marker ridges.
The instrument was connected to a robot arm providing ±0.05 mm repeatability at the proximal end of the instrument shaft.The instrument mount was machined precisely and connected with the robot wrist via a keyway such that the instrument was aligned with the rotation axis of the most distal motor of the robot.Robot joint angles were controlled using a proportional-integral control law and trapezoidal velocity profiles.The instrument was commanded through a specified trajectory such that the marker remained inside the image volume.Yaw angle was held constant to ensure consistent instrument visibility given the modest thickness of the ultrasound image volume.The position and orientation of the instrument marker in robot coordinates was calculated by (39) where is the position and orientation of the robot wrist relative to the base, as calculated via the joint encoders, and is the position and orientation of the marker with respect to the robot wrist, which was constant.The quantity was initially estimated based on the instrument mount dimensions and caliper measurements of marker location and rotation.
Given that the instrument marker position is measured via both the robot and the ultrasound transducer, the experimental system forms a closed kinematic chain, which is defined by combining ( 39) with (1) (40) where , , and relates the robot base frame to the ultrasound image frame.The components of the experimental system are shown in Fig. 19.
The position and orientation of the instrument marker in the ultrasound image volume, was recorded as it was moved throughout the trajectory.These measurements were compared to those from the robot via a modification of (40) (41) where E is the tracking error.Since the probe was rigidly mounted with respect to the robot base, the transformation was constant.It was determined by a batch calculation over the complete tracking data set such that the position components of E were minimized in the least squares sense.Spatially, this is equivalent to aligning the marker positions measured in the ultrasound image volume to those measured via the robot such that their gross orientation is equal.This assumes no static measurement error.Between images 400 and 500, as shown in Fig. 20, marker velocity was too great (about 15 mm/s) for the filter to maintain tracking.The excessive error triggered the innovation threshold and the filter was reset when tracking did not converge.Note that this velocity exceeds the maximum linear velocity of 10 mm/s that is typically encountered in laparoscopic surgery [26].Velocities employed in intracardiac surgery are even lower due to workspace constraints and tissue sensitivity.
The water tank accuracy achieved, shown in the following table, surpasses many other available forms of instrument tracking.To evaluate marker tracking in a more realistic testing scenario, ex vivo experiments are described below.

D. Ex Vivo Cardiac Experimental Results
To characterize marker performance in a more realistic testing situation, five trials were conducted in which a porcine heart was suspended in a water tank and the marker was inserted inside the left atrium.The marker for these experiments has the same geometry and size as that of Fig. 17, however, it was created as an integrated unit with the instrument shaft from a photopolymer using a 3-D printer (Objet, Inc.) and is shown in Fig. 22.In the experimental system shown in Fig. 23, the instrument location is fixed inside the heart and the ultrasound probe was moved incrementally on a linear stage.Stage location was varied manually over a 1 cm range-of-motion and 3-D ultrasound volumes were recorded for each stage position.
In this experiment, the instrument tracking algorithm required an accurate initial detection for the first frame owing to instrument imaging artifacts and the low contrast between the instrument and surrounding tissue.To ensure accuracy, target detection in the initial frame was performed by reducing the ultra-sound transducer power to −16.2 dB.This had the effect of enhancing the contrast of the instrument at the expense of tissue visualization.All subsequent frames were processed at normal transducer power 0 dB, which provides visualization of both the instrument and the cardiac tissue.
During tracking, the search space was confined to a volume-of-interest centered at the previously estimated location.Using Kalman filtering, the subsequent detected marker position was introduced as a new measurement for updating the innovation filtering process.Since the real-time 3-D ultrasound volumetric images are updated at 28 Hz and typical instrument motions in intra-body surgeries are performed slowly [26], it is typicallly sufficient to confine the search space to ±5 voxels relative to the previous detected target.Fig. 24 plots the spatial path of the instrument tip as computed from probe displacement and as estimated from the marker.Comparing the ex vivo tracking error with that of the water tank experiments, average tracking error is 2.4 times larger at 1.7 mm.Standard deviations for the individual coordinate directions are reported in Table I.

V. Discussion
The experimental results indicate that marker-based instrument tracking is a promising method for improved image-based instrument navigation.Applications include augmented reality surgical displays and image-based servoing of robotic instruments.As demonstrated by the experimental results, the accuracy of the proposed approach depends on the quality of the ultrasound imaging.This contrasts with optical and electromagnetic tracking methods that maintain a fixed accuracy regardless of the target size or image resolution.The proposed markers are not subject, however, to errors arising from instrument flexure (optical techniques) and field distortion (electromagnetic tracking) and, furthermore, as high-frequency, high-resolution probes are used to image small targets, the accuracy of marker tracking will increase with the use of such probes.
Since the proposed markers are image based, they eliminate both the labor and the error involved in registering the tracking coordinate frame to the image frame.In addition, they are completely passive, which allows them to be implemented easily, without prior integration with the imaging system.Finally, using the design methods presented, marker shapes can be tailored for particular navigation tasks.
As suggested by the need for marker shape calibration, further study is required to realize the full potential of marker-based instrument tracking.Calibration adjusts for inaccuracies in marker manufacturing, which can be improved through more accurate construction techniques.It also adjusts for the shifts in ridge location produced by the real acoustic response of ridges, which are likely dependent on properties of the imaging system including the transducer beam spacing, focal depth, and point spread function.
In addition, the instrument size and orientation, ridge size and shape (curvature, slope and discontinuities), and the materials used for marker construction would all contribute.A complete study of the acoustic response of markers would produce a function, with which a theoretical marker shape may be transformed into an actual marker shape, as it appears in the ultrasound image.The theoretical analysis presented in this paper can be used to analyze the properties of shapes transformed by this function.Additional study may also be directed toward shape optimization, specifically to find continuous or semi-continuous shapes that minimize ( 22).This would involve determining a parameterization for projected marker curves, such that all possible solutions may be explored.Such a result would also enable more precise calculation of the lower bound on ∊ θ .Furthermore, the image processing methods demonstrated here were designed to illustrate the potential of the marker technology.Alternate methods will likely improve tracking accuracy and bandwidth.For example, a GPU-based shaft segmentation algorithm has been demonstrated to achieve real-time performance during in vivo cardiac experiments [27].Instrument marker attached to a laparoscopic grasping instrument.Coordinate frames used in marker-based instrument tracking.Example semi-continuous marker shape.Plot of the bounding simplex defined by marker design constraints, including the example marker curve.Plot showing the projected path of the example marker and the projected bounding simplex.The magnitude of q is exaggerated in the figure for illustration.The , axes of the plot are defined on the plane orthogonal to the vector u.Two planes orthogonal to u intersecting the bounding simplex for n = 3.The inset shows the projected constraint simplex and the projected area of intersection between the planes and the constraint simplex.Example marker shape from Fig. 4 after reorientation so that f′ T u = 0.The marker shape in   Example marker shape from Fig. 4 reoriented so that f′ T u = 0 and reparameterized to minimize J θ .Example two-ridge marker shape (a) and the corresponding marker curve (b).Example four-ridge discontinuous marker shape, (a), and the corresponding projected marker curve, (b).Experimental system for marker tracking inside an ex vivo porcine heart.Actual and estimated instrument path corresponding to linear probe motion.Units are in millimeters.
y, and z describe the marker position, ϕ and ψ describe the orientation of the instrument shaft axis, and θ describes the roll angle.Primed variables indicate derivatives with respect to time.The singularities associated with the Euler angle formulation are not problematic because they are outside the instrument's possible range of motion within the image.The Kalman filter is based on a constant velocity process model.It is used for trajectory smoothing and for added robustness in the shaft segmentation and marker detection algorithms.An innovation threshold is also implemented, whereby the filter is reinitialized automatically when the difference between the predicted and measured marker coordinates becomes too large.
Finally, the result for was used to refine the estimate of marker rotation about the shaft axis in such that measurement error in θ had zero mean.Measurements in which the innovation filter detected a large jump in position or angle were not considered in the calculation of tracking error.Since this determination is made on the basis of the innovation filter, this omission represents a realistic correction.Tracking results of water-bath experiment are shown in Figs.20 and 21.Solid lines indicate measurements within the ultrasound image volume; dotted lines indicate measurements from the robot.Segmentation algorithms operated at approximately 3 volumes per second when implemented in Matlab on a 3 GHz desktop.Algorithms were neither compiled nor optimized for speed; volume rates closer to that provided by the ultrasound system would be achieved through compilation and optimization.

Fig. 3 .
Fig. 3. Ridge location along shaft axis, f i , versus sleeve angle, θ.The plot in (a) corresponds to the 3-D marker in (b).

Fig. 7 .
Fig. 7.Error factors for the example marker shape depicted in Fig.4.
(a) corresponds to the marker curve in (b).

Fig. 17 .
Fig. 17.Instrument and marker used in tracking experiment with the marker shape shown in the inset (a) and a cross-sectional image of the marker extracted from an ultrasound image volume (b) The dots indicate the estimated ridge locations and the line indicates the estimated shaft axis.

Fig. 18 .
Fig. 18.Actual marker shape obtained by rotating the instrument of Fig.17in place about θ at the five image locations indicated.

Fig. 19 .
Fig. 19.Components of the experimental system.M: Marker's body frame; R: Robot base frame; im: ultrasound image frame; A: intermediate frame defining the plane of apex of ultrasound volume and instrument shaft.

Fig. 20 .
Fig. 20.Instrument tracking result for x, y, and z.The dotted lines indicate the actual instrument trajectory as measured by robot joint encoders.

Fig. 21 .
Fig. 21.Instrument tracking result for θ, ϕ, and ψ.The dotted lines indicate the actual instrument trajectory as measured by robot joint encoders.

Fig. 22 .
Fig. 22. Instrument and marker used in ex vivo tracking experiments.The marker shape and size are identical to Fig.17. .