Surface extraction with a three-dimensional freehand ultrasound system

Ultrasound in Med. & Biol., Vol. 30, No. 11, pp. 1461–1473, 2004 Copyright © 2004 World Federation for Ultrasound in Medicine & Biology Printed in the USA. All rights reserved 0301-5629/04/$–see front matter

doi:10.1016/j.ultrasmedbio.2004.08.020

● Original Contribution SURFACE EXTRACTION WITH A THREE-DIMENSIONAL FREEHAND ULTRASOUND SYSTEM WAYNE Y. ZHANG,* ROBERT N. ROHLING* and DINESH K. PAI† *Department of Electrical and Computer Engineering, University of British Columbia, Vancouver, British Columbia, Canada; and †Division of Computer and Information Sciences, Rutgers–the State University of New Jersey, Piscataway, New Jersey, USA (Received 10 March 2004; revised 2 August 2004; in final form 13 August 2004)

Abstract—This paper presents a system for acquiring three-dimensional ultrasound data and extracting surfaces of the examined structures. The acquisition is performed freehand with a PC-based two dimensional ultrasound machine and an optical tracker. The extraction of surfaces from ultrasound data are normally inhibited by speckle, shadowing and gaps in the acquired data. A new method is developed that extracts a surface directly from the irregularly spaced, noisy freehand ultrasound data. The freehand data are first interpolated with radial basis functions and then a mesh is formed along an isosurface of the functional interpolation. The ability of radial basis functions to smooth speckle and interpolate across gaps is demonstrated on a series of experiments with phantoms and human tissue in a water bath. The geometry of the extracted surface matches the external measurements with an average difference ranging from 0.8 to 2.9 mm. These differences are within the range of errors from calibration, resolution and landmark localization. The experiments also show the ability to create continuous and realistic surfaces from scans that require multiple sweeps over a structure. (E-mail: [email protected]) © 2004 World Federation for Ultrasound in Medicine & Biology. Key Words: Three dimensional ultrasound, Freehand, Surface Extraction, Radial basis functions.

angled portions of the boundary are poorly received. Figure 1 shows an example of this. Speckle, shadowing and other artifacts also obscure organ boundaries, as shown in Figure 2. Another problem arises from the nature of the 3D acquisition technique. The two most popular ways to acquire 3D ultrasound are with a freehand technique or with a specialized 3D probe. Freehand ultrasound uses a conventional 2D ultrasound probe combined with a position sensor to track the motion of the probe as it is moved across the body. By sweeping the probe from tip to tail of an organ, a set of 2D images and their corresponding positions are obtained. Figure 3 shows an example. The spacing among the images is determined by the operator’s hand, so it is nonuniform and typically much larger than the spacing of the pixels in the individual images. Specialized 3D probes are also hand-held but a volume of data is acquired automatically under the probe. The field-of-view of these specialized 3D probes is still limited by the footprint of the probe, so complete images of larger organs, such as the fetus, must be created by moving the probe to different viewpoints (e.g., front, back and sides) and combining the views

INTRODUCTION The objective of this study is to demonstrate a working three-dimensional (3D) ultrasound system for obtaining surface models of anatomical structures. The main application of current 3D ultrasound systems has been improved visualization of the anatomy. Another area is computer assisted surgery where 3D ultrasound is used increasingly as a tool for building models of organs (Fenster et al. 2002; Welch et al. 2000). The latest generation of real-time 3D ultrasound machines offer the advantages of speed, cost and portability over other modalities such as magnetic resonance (MR) and computed tomography (CT) imaging. The extraction of organ surfaces from 3D ultrasound is a bigger challenge than from MR or CT. One of the problems is the nature of ultrasound. The directionality of pulse-echo imaging makes it difficult to produce uniform boundaries in the image because echoes from

Address correspondence to: Dr. R. Rohling, Department of Electrical and Computer Engineering, University of British Columbia, 2356 Main Mall, Vancouver, B.C. V6T 1Z4, Canada. E-mail: rohling@ ece.ubc.ca 1461

1462

Ultrasound in Medicine and Biology

Fig. 1. Ultrasound image of a continuous rubber surface (toy duck) in a water bath. Two arrows indicate the gaps that exist between the wing and the torso. These sections of the surface are approximately parallel to the pulse-echo beam, so no reflections are received.

together. The end result with either 3D acquisition method is a noisy 3D data set with likely large gaps, irregular spacing, overlaps and ambiguities in the location of the organ boundary. If the main goal is visualization, then reconstructing the ultrasound data into a regular voxel array and applying ray-tracing techniques can produce realistic renderings of the anatomy. Many visualization techniques to overcome the difficulties of 3D ultrasound are available in the literature (Steen and Olstad 1994; Deng et al. 1996; Sakas and Walter 1995; Nelson and Elvins 1993; State et al. 1994; Fenster and Downey 2001). Applica-

Fig. 2. Ultrasound image of a human finger. Label 1 shows speckle and enhancement artifacts. Label 2 shows the skin surface. Label 3 shows the bone surface. Label 4 shows a vein and Label 5 shows a dark shadow interior to the bone where ultrasound cannot penetrate.

Volume 30, Number 11, 2004

Fig. 3. A freehand 3D ultrasound examination consists of a set of cross-sectional images whose nonregular spacing is determined by the movement of the operator’s hand. This is a single sweep of a human finger.

tions in computer assisted surgery usually require quantitative analysis, such as organ deformations, landmarks and organ volumes. For these applications, an explicit model of the organ surface is required. This paper focuses on the problem of extracting the surface of organs bounded by fluid. Examples include the fetus in amniotic fluid, the bladder, gall bladder and other fluid-filled organs. The surface extraction is therefore characterized by following bright boundaries that are bounded by darker fluid. A common method of generating an explicit surface is to use the reconstructed voxel array with surface extraction techniques used in MR and CT. By first applying a segmentation technique, such as thresholding, one of the popular marching methods can be used to extract a surface (Lorensen and Cline 1987). Early marching algorithms encountered problems with nonmanifold surfaces and with low quality surfaces constructed from many low-aspect-ratio triangles. The manifold surface problem has been addressed (e.g., Nielson and Hamann 1991). Attempts to improve the surface quality include the marching tetrahedra techniques (Treece et al. 1999b) and mesh simplification techniques (Cignoni et al. 1998). Even with new algorithms, the initial reconstruction step remains a problem because it must combine images from multiple viewpoints together while solving the problems of data overlap and gaps. Reconstruction normally reduces the resolution of the data and can introduce artifacts (Rohling et al. 1999). Setting the voxel size equal to the pixel size quickly overwhelms the memory of a personal computer and results in many unfilled voxels. Normally, the voxel size is chosen to be large enough to bridge gaps in the data, which results in many pixels falling into a single

Surfaces from 3D ultrasound ● W. Y. ZHANG et al.

voxel. There have been previous efforts to minimize the loss of resolution by representing the ultrasound data with smooth functionals (Rohling et al. 1999) or tetrahedron data structures (Roxborough and Nielson 2000). Performance remains a problem and is balanced against resolution to achieve a compromise. None of the ultrasound reconstruction papers have directly addressed the problem of surface extraction. Another class of methods to extract the surface is active surface models. Active surface models work on 3D data by expanding a “balloon” to match the boundary of the organ (Cohen 1991; Singh 1998). The properties of the model that constrain the allowable shapes are normally determined a priori. This constrains the topology of the surface that can be modeled. Another limitation of this technique is that the balloon must be initialized manually and the final surface depends on the initialization. It is difficult to fit active surface models to organs when there is ambiguity in the boundary data, the object’s shape is unknown and the surface is allowed to branch. Another method is to segment the surface contour in each of the 2D cross-sectional images and then to join the contours together to form a single mesh of triangles. This is a two step procedure. Segmentation can be performed manually (Prager et al. 2002) or by other means. Active contours are commonly used to semiautomate the procedure when the contours are clearly defined (Chen et al. 2001). Terzopoulos et al. (1988) initially proposed the active contour extraction algorithm, often called snakes. The snake is usually initialized manually as a closed contour near the organ boundary. The snake is then iteratively deformed to the nearest organ boundary by minimizing a function based on the image features and snake curvature. Typically, the weighting factors in the function are determined a priori. Choosing a large allowable snake curvature means that complex shapes can be fit but are more prone to leak out of the organ when the boundary has gaps. McInerney and Terzopoulos (1996) extended the geometric and topological adaptability of snakes, while retaining all the features of traditional snakes by superposing a simplicial grid over the image domain and using this grid iteratively to reparameterize the deforming snake model. But these methods are inherently 2D. This means that the 3D nature of the surface is not used and ambiguities in the surface boundary must be resolved with the individual 2D contours. For example, a thin tube is best described by a set of circular contours taken along the long axis of the tube, but contours derived along other axes may be more difficult to combine into a tube-like surface (Treece et al. 1999a). The second step of combining the contours into a single surface is done by making correspondences between neighboring contours while resolving branching

1463

surfaces. A low quality surface can easily result, so newer techniques, such as shape-based interpolation, have been developed to improve surface quality (Treece et al. 2000). The set of contour points extracted from the 2D images can also be combined and considered as a set of unorganized surface points in 3D. This avoids the explicit matching and correspondence problem of contour matching. Explicit methods of this type attempt to triangulate a surface directly from the points, such as by using Delaunay triangulation and Voronoi-based techniques (Boissonnat 1984). These often include an additional step to remove nonsurface triangles after the initial triangulation (Edelsbrunner and Mucke 1994; Bernardini et al. 1995; Amenta et al. 1998). The quality of the surface is sensitive to noise in the points and large gaps cannot be smoothly filled. Implicit methods define the surface from a functional representation. An example is zero sets (Hoppe et al. 1993), where a signed distance function in R3 is created from the points. The level set method is another example. The level set method was devised by Sethian (1996) as a simple and versatile method for computing and analyzing the motion of an interface in two or three dimensions. The original idea of the level set method is quite simple. A closed simplex (a curve in 2D and a surface in 3D) is the zero level set of a higher dimension function. Leventon et al. (2000) successfully applied the level set method in segmentation for 3D medical images. Zhao et al. (2001) constructed a weighted minimal surface-like model using differential geometry and partial differential equations and optimized its level set formulation. Radial basis functions (RBFs) can also be used to create a distance function where the surface is defined by an isovalue of the function. Radial basis functions are a mathematical tool for interpolation. They consist of placing a symmetric function (the basis function, such as a Gaussian function) at each data point then weighting the functions so that their sum passes through all of the data points. The sum of the basis functions can then be evaluated anywhere, including between the data points. RBFs have enjoyed a recent rise in popularity because they can smoothly interpolate across gaps. They can also approximate or interpolate the data depending on the level of noise. We propose a similar method using RBFs but work directly on the pixels from the set of 2D images instead of the extracted surface points from segmentation. An earlier study (Zhang et al. 2002) showed the feasibility of using RBFs on ultrasound data, but was limited to open surfaces, one main viewpoint and small gaps in the data. In the following sections, we will describe a method that attempts to

1464

Ultrasound in Medicine and Biology

Fig. 4. The freehand 3D ultrasound system consists of the Ultrasonix ultrasound machine, a standard 2D probe, a group of LEDs attached to the probe and a tracking device to measure the movement of the probe over an object.

● ● ● ● ● ●

combine multiple views of an object, fit a function to noisy irregularly spaced ultrasound pixel data, extract an explicit surface along an isovalue of the function, solve the problem of holes in the data by smooth interpolation across them, retain the high resolution details of the surface, and compute the results efficiently for practical applications.

The purpose of this paper is to extract the surfaces of several measurable objects using RBFs. These objects will first be measured externally, then scanned with a 3D freehand ultrasound system in a water bath. The goal is to demonstrate that the surfaces extracted with RBFs from the 3D ultrasound data correspond to the real objects. Comparisons will be made visually and by calculating the difference between the extracted surfaces and the external measurements.

Volume 30, Number 11, 2004

ital Ltd., Waterloo, ON, Canada) are attached to an 8 MHz linear array probe, as shown in Figure 4. The positions of the LEDs are tracked with an accuracy of 0.1 mm (Rohling et al. 1995) and converted into position and orientation of the probe. To convert into position and orientation of the image plane, another transformation is needed. This last transformation is determined by calibration before each examination. Calibration is performed here using the z-wire calibration phantom (Comeau et al. 2000; Pagoulatos et al. 2001). The phantom consists of a number of thin nylon wires placed in holes accurately drilled by a CNC milling machine. A closed loop is formed with the tracking system so that the calibration transformation remains as an unknown in the equations for the loop. Figure 5 shows L the loop where U T is the unknown transformation (in homogeneous matrix notation) between the coordinate system of the ultrasound image (U) and the LEDs on the probe (L). O L T is the measured location of the LED coordinate system with respect to the Optotrak (O) and O P T is the location of the z-wire phantom (P) with respect to the Optotrak. POT is measured by a pointing device and remains fixed, since the box is stationary with respect to P the Optotrak. U T is measured by locating the dots in the ultrasound image that correspond to the z-wires in the phantom. So, by observing the locations of the wires in an ultrasound image, the transformation equations can be written for the loop: T · ULT ⫽ POT · UPT.

O L

(1)

L Although it is possible to solve for U T from a single image of the z-wires, we use multiple images and solve

METHODS Acquisition system The freehand 3D ultrasound system is based on a PC-based ultrasound machine (Ultrasonix 500RP, Ultrasonix Medical Corporation, Burnaby, BC, Canada) and a programming interface. The combination is called PUPIL (programmable ultrasound platform and interface library). PUPIL has an open architecture that allows easy access to the image formation pipeline and control of the acquisition parameters from software (Rohling et al. 2003). Five LED sensors from a trinocular tracking device (Optotrak 3020, Northern Dig-

Fig. 5. The four coordinate systems used for calibration are the fixed Optotrak system (O), a system centered on the LEDs (L), a system at the top-left corner of the ultrasound image (U) and a system for the z-wire phantom (P).

Surfaces from 3D ultrasound ● W. Y. ZHANG et al.

an over-constrained set of equations, to minimize the effect of random error on the measurements. To test the reliability of the calibration, the z-wire phantom is replaced with a single point target that is imaged from 15 different points of view. This is similar to the calibration test performed by (Prager et al. 1998). The result of imaging the point target is a cloud of points with standard deviations of 0.39 mm, 0.37 mm and 0.58 mm in the x, y and z directions, respectively. The maximum errors (the maximum minus the minimum axes of the cloud) are 1.68 mm, 1.49 mm and 1.57 mm for x, y and z. This is approximately the same accuracy as that reported by (Prager et al. 1998) with a different position tracker, suggesting that the accuracy of calibration is determined mainly by the resolution of the ultrasound machine and not by the accuracy of the optical tracker used here. The calibrated system is now used to scan a number of objects. The main purpose of the experiments is to validate the surface extraction method with known objects. The objects are a toy duck, a human hand, a toy fish and a simulated fluid-filled cyst. The objects are scanned in a water bath to provide acoustic coupling. Each object is scanned from multiple points of view, because the objects are either too large for a single view or some surfaces (such as the backsides) are not visible from one viewpoint. Once scanned, the data are stored as a set of 2D images. Each image (i ⫽ 1,. . .,N) is labeled with its L corresponding position in space (O L T · UT). The surface extraction algorithm contains up to three steps. The first step is to interpolate the data using RBFs. The second step is to extract the surface along the isovalues of the function. The output is a mesh of triangles. The third step is optional: any holes in the mesh are repaired. From a complete mesh representation of the object, the surface can be viewed. We can also calculate the volume of a closed surface (e.g., the simulated cyst). In this paper, we wish to examine the ability of the new method to combine several views of an object, fill holes and reproduce realistic renderings of the surface, despite the presence of noise in the data. Each of the steps is now explained in more detail. Radial basis function interpolation Many interpolation techniques are available for irregularly spaced 3D data, but RBFs offer a unique combination of features. They are a compact form of representing sparse, nonuniform data. They can either interpolate or approximate the data points. They can be evaluated anywhere to extract a mesh at the desired resolution. Local gradients, higher derivatives and surface normals can be easily calculated from the smooth and continuous RBFs. Isovalues can be used to extract

1465

manifold surfaces and the surfaces can be easily smoothed (Carr et al. 2001). Ultrasound acquisition can be considered as a sampling of tissue structure with a finite point spread function equal in size to the resolution of the pulse-echo system. Since ultrasound is normally used for imaging soft tissue structures, it is reasonable to assume that the sampled surface of the structures can be modeled by a single smooth continuous function. The proposed RBFs can model such surfaces without any restrictions on topology. More formally, consider a set of 3D points at locations xi, i ⫽ 1,. . .,N, where each point has an intensity value pi. The basic idea of RBF interpolation is to find a spline S(x) that fits the data points as smoothly as possible. With this constraint, we find the S(x) that minimizes the following: N

兺 ⱍp ⫺ S共x 兲ⱍ ⫹ ␻I共s兲. 2

i

i

(2)

i⫽1

The first term is the deviation of the spline from the data points, which represents the interpolation error. The second term, I(s), is a smoothness function of surface s, normally chosen to minimize one or more partial derivatives of the function. The parameter ␻ is the weight which determines the relative contribution of these two components. The general solution can be expressed as follows: N

S共x兲 ⫽ t共x兲 ⫹

兺 ␭ R共ⱍx, x ⱍ兲 i

i

(3)

i⫽1

where |x, xi| is the Euclidean distance, t(x) is a polynomial called the trend function, R(x, xi) is a radial basis function, xi are the centers of the RBF and ␭i are the weights of each RBF. The form of the RBF depends on the choice of I(s). The biharmonic spline R(r) ⫽ r (where r ⫽ |x, xi| is the distance between the point of interest x and the i’th basis function located at xi) means that the basis function is zero at x ⫽ xi and increases linearly away from xi. The biharmonic spline basis function minimizes the integral of the second derivative. This is a reasonable measure of smoothness and is used here. Given the form of the solution from eqn 3 and a set of points pi, we write the set of linear equations for interpolating the points. Using the biharmonic spline in 3D means that the trend function is t(x) ⫽ c1 ⫹ c2x ⫹ c3y ⫹ c4z. ␭i and t(x) are the unknowns. To maintain a finite I(s), the following side conditions must also be satisfied: N

兺 ␭ t共x 兲 ⫽ 0. i

i⫽1

i

(4)

1466

Ultrasound in Medicine and Biology

The set of equations and the side conditions can be written in matrix form as follows:

冉 冊冉 冊 冉 冊 ␭

A T T

T

0

c

⫽

p 0

(5)

where

ⱍ ⱍ

Ai,j ⫽ R共 xi,x j 兲, i, j ⫽ 1, ..., N, Ti,j ⫽ t j共xi兲, i ⫽ 1, ..., N, j ⫽ 1, ..., l, p is the vector of intensity of the pixels pi, {t1,. . .,tl} is a basis for polynomials and ␭ is a vector of the coefficients ␭i. The vector p consists of all pixels on all ultrasound images used in the interpolation. This means they are spatially distributed in the manner of the freehand acquisition: they are regularly ordered in columns and rows within an image plane, but the planes are not regularly located in space. The planes may overlap or be spaced far apart, depending on the particular examination. It is assumed that the points in p follow a real surface, otherwise the equations become ill-conditioned (e.g., for a random cloud of points) and a good surface is not obtained. The vector c holds the coefficients (c1, c2, c3, c4) of the trend function t(x). T is a matrix where the ith row contains (1,xi, yi, zi) and l, the degree of this polynomial, equals 4. The weight factor ␻ of eqn 2, can be included as follows:

冉

A ⫹ ␻I T TT

0

冊冉 冊 冉 冊 ␭ c

⫽

p 0

(6)

A nonzero value of ␻ means the RBF passes close to, but not exactly, through the data points. The result is an approximating function that can be used to smooth speckle and other noise. In this paper, the term interpolation is used for either exact or approximating functions and specific cases of approximation are mentioned explicitly. Since solving these equations has complexity O(N3), the speed will be prohibitively long for the large number of points in 3D ultrasound. Instead, fast evaluation of RBFs is performed using the fast multiple method (Greengard and Rokhlin 1987). A particular adaptation (Beatson et al. 2001) of the fast multipole method was implemented elsewhere and tested on laser range finder data (Carr et al. 2001). This version has O(N log N) complexity and is used here. The complexity is reduced by recognizing that data points far from a point of interest have little effect on the local interpolation. The interpolation is computed more efficiently by simplifying the effects of distant points to a small, but finite, accuracy. Details of the approach are given by (Beatson and Greengard 1997).

Volume 30, Number 11, 2004

Looking again at the biharmonic spline R(r) ⫽ r, we see a noncompactly supported basis function. These basis functions will be shown later to be well suited to interpolation across gaps in the data. The main limitation is still speed. To further improve performance, the number of RBF nodes can be reduced in a process called center reduction (Carr et al. 2001). This simply reduces the number of nodes while still maintaining the interpolation to a desired accuracy. The accuracy is chosen to be much smaller than the surface features. A final tool to reduce the size of the problem is a masking step at the preprocessing stage. This tool (described by (Zhang et al. 2002) and kept here without modification) simply masks data points that have values physically far from the expected boundary. By manually applying the mask before computing the RBFs, the surfaces appear to be identical to those computed without masking, but the memory requirements are lower. The masking also makes it easier to visualize the raw data points. Surface extraction The goal of surface extraction is to produce a high quality mesh of triangles along an isovalue of the RBFs. A low quality mesh contains a large number of triangles with poor aspect ratio. Since conventional marching algorithms are susceptible to this problem, a marching tetrahedra variant (Carr et al. 2001) is used with a mesh optimization step from Treece et al. (1999b). Marching methods are normally designed for volume data where samples consist of regularly spaced voxels. It would be inefficient here to sample the RBF at high resolution throughout the volume, so a surfacefollowing algorithm is used instead. The idea is to use a seed point on the surface and emanate wavefronts of facets that follow the surface until they meet another wavefront or reach the bounding box. The seed point is chosen to be an RBF center that is located on the surface. If an RBF center is chosen that does not lie directly on the surface, the nearest point on the surface is calculated. These calculations are readily performed, because the value and gradient of the RBFs are easily computed anywhere. The seeds are chosen automatically and only a small number of seeds are needed to ensure that the wavefronts cover all surfaces. Mesh optimization is performed on-the-fly as the wavefronts advance. A mathematical and graphical description of this process is given in Carr et al. (2001). Mesh repair The extracted surface may still contain gaps where the ultrasound data insufficiently delineated the boundary (e.g., from shadowing). The main principle of mesh repair is to use the shape and continuity of the rest of the surface to fill the gaps. Again, the method by Carr et al.

Surfaces from 3D ultrasound ● W. Y. ZHANG et al.

Fig. 6. An edge view of the surface is shown. (a) The surface points and the surface normals. (b) The off-surface points.

(2001) is used. The main principle is to create a new set of data points that represent a signed distance-to-surface field. In other words, the original ultrasound data points are discarded after the first surface is created. A new set of points is created by first sampling along the first surface to get surface points and then generating addi-

1467

tional off-surface points. The surface points are given a value of zero. Each off-surface point is given a value that is the signed distance of the point to nearest spot on the surface. Points outward of the surface are given positive values; points inward are given negative values. To generate the off-surface points, the surface normals are used: see Figure 6. Given a point p on the surface, with normal direction n, off-surface points po (outward) and pi (inward) are generated by po ⫽ p ⫹ ␣n and pi ⫽ p ⫺ ␤n, where ␣ and ␤ are the projection distances of the off-surface points. Constant values for ␣ and ␤ can cause problems at thin sections and sharp corners, so the distance and density of the points are validated before the points are used. A simple method to validate the correct choice of the projection distance is to find the nearest point q on the

Fig. 7. Examples of ultrasound images. (a) An image of the duck phantom. Notice the noise and gaps in the boundary. (b) An image of the fish phantom. (c) An image of the hand near the middle finger knuckle. (d) An image of the simulated cyst.

1468

Ultrasound in Medicine and Biology

Volume 30, Number 11, 2004

Once a validated set of points representing distanceto-surface are obtained, they are used for fitting RBFs and isosurfacing in the same manner as explained for ultrasound data. The only difference is that the isosurface is extracted along an isovalue of zero. Because the erroneous portions of the ultrasound data (e.g., from shadowing) are no longer used, the RBFs can interpolate more easily across the gaps. The mesh repair algorithm can only fill small holes in the surface where there is sufficient neighboring data points to be able to estimate the surface across the gap. For organs that fall outside the field-of-view, or are obscured by a large amount of shadow, or are not closed, it is not possible to extract a closed surface. In those cases, only the visible portion of the surface will be extracted.

Fig. 8. Duck phantom. (a) A photograph of the real surface. (b) The point cloud. Notice the examination is performed with two sweeps with some large gaps between the sweeps. The first sweep covers the head from left to right starting from the beak. The second sweep covers the torso from right to left starting from the tail. (c) The surface extracted from using approximately 99,000 points. (d) Another examination of the duck but the surface is extracted from only nine thousand points. (e) A final examination where both the interior and exterior surfaces of the rubber phantom are extracted. A portion has been cut away to show the separation between the two surfaces.

surface for each off-surface point poff (created from onsurface point p). If q is different than p, then the projection distance is reduced until the nearest point is p.

Validation The extracted surfaces are normally visualized directly using the visualization toolkit (Schroeder et al. 2003). Since the surfaces of the objects can be examined directly, the dimensions of the extracted surface can be compared with external measurements. For the duck, fish and human hands, the external measurements are made with calipers. The geometry that is measured is chosen according the ability to locate well-defined landmarks. For these surfaces, landmarks are chosen as high curvature points that are clearly identified in both the real object and the extracted mesh surface. In addition to the external caliper measurements, the human hand was also measured using a laser range finding system (FastSCAN Cobra, Polhemus Inc., Colchester, Vermont). The laser range-finder was used because calipers require contact with the skin which can deform. The laser range-finder does not require contact; so measurements give an indication of the variability introduced in the measurements from contact measurements. The calipers have an accuracy of 0.3 mm and the laser range-finder has an accuracy of approximately 1

Table 1. Computation times for the examinations. All computations were performed by a 1 GHz Pentium III (Intel Corp., Santa Clara, CA, USA), except the fish, which used a 2 GHz Pentium IV. Experiment Duck Fish Hand-closed Hand-open Gel pocket

Figure

No. points

Fit time (s)

Surface time (s)

8(c) 8(d) 8(e) 9(d,e) 10(c) 10(d) 10(e) 11(c) 11(d) 11(e) 12(a,b,c)

98654 9907 99087 123874 99122 99122 99122 99174 99714 99714 241545

1195 58 1047 937 1236 1636 1803 1122 1288 1287

317 70 391 25 464 450 449 676 252 59 8420 (multistage)

Surfaces from 3D ultrasound ● W. Y. ZHANG et al.

Fig. 9. Fish phantom. (a) A photograph of the real surface. (b) and (c) The point cloud. Notice that the examination is performed with four sweeps with large gaps where the front sweeps meet the back sweeps. (d) The mesh of triangles representing the surface. (e) The extracted surface.

mm. The total accuracy of the external measurements is a combination of the measurement tool, landmark localization and tissue deformation. If the surface is closed, then the volume of the object can be calculated. We calculate the volume of the simulated cyst from the extracted surface mesh using an

Fig. 10. Human hand. (a) A photograph of the real hand. The distances measured for validation are labeled. (b) The point cloud. Notice that the examination is performed with three sweeps with some gaps and overlap among the sweeps. (c) The surface extracted without smoothing. (d) The surface extracted with moderate smoothing (␻ ⫽ 3). (e) The surface extracted with heavy smoothing (␻ ⫽ 9).

1469

Fig. 11. Human hand. (a) A photograph of the real hand. The distances measured for validation are labeled. (b) The point cloud. Notice that the examination is performed with five sweeps with some gaps and overlap among the sweeps. (c) The surface extracted with a resolution of 0.1 mm. The number of triangles in the surface is 245,420. (d) The surface extracted at 1.0 mm. The number of triangles in the surface is 63,200. (e) The surface extracted at 2.5 mm. The number of triangles in the surface is 10,776.

algorithm based on the discrete form of the divergence function (Alyassin et al. 1994). This volume is then compared with the real volume measured by water displacement in a graduated cylinder.

Fig. 12. Simulated cyst. (a) The first surface extracted from the ultrasound data. Some gaps remain where the boundary is obscured in the ultrasound data. (b) The mesh of triangles indicating good triangulation, but with several holes. (c) A fully enclosed surface created from the first surface after interpolating over the holes. Notice that the surface detail is almost exactly retained from the surface before mesh repair.

1470

Ultrasound in Medicine and Biology

Volume 30, Number 11, 2004

Table 2. Comparison of measurements made of the duck phantom. The first measurements are from the constructed surface mesh. The second measurements are from calipers Measurement

First point (x,y,z) (mm)

Second point (x,y,z) (mm)

Distance (mm)

Calipers (mm)

Difference (mm)

Bill to back of head Head height Bill to tail Bill to wingtip

(⫺138.0,1828.1,⫺96.9) (⫺133.7,1869.1,⫺78.7) (⫺138.0,1828.1,⫺96.9) (⫺138.0,1828.1,⫺96.9)

(⫺141.1,1899.7,⫺95.9) (⫺128.3,1874.4,⫺123.6) (⫺141.9,1939.2,⫺109.4) (⫺128.2,1937.4,⫺125.2)

71.7 45.5 111.9 113.3

75.0 49.5 111.5 115.1 Average:

3.3 4.0 0.4 1.8 2.4

RESULTS Examples of the ultrasound images used to construct the surfaces are shown in Figure 7. The surfaces are depicted as bright curves. The inherent speckle and shadows are apparent. The duck phantom was scanned with two sweeps, one for the head and another for the torso. The purpose of this examination is to demonstrate the ability to join sweeps together when substantial gaps exist between them, as shown in Figure 8. The surface is interpolated across the gaps with approximately the same shape as the real surface. The bumpy nature of the surface is a result of the speckle in ultrasound, yet many details remain. The gaps among the slices in the point data are visible as vertical creases. The creases represent misalignment among the cross-sectional slices and is likely due to errors in tracking and calibration. A second scan of the duck phantom was performed with fewer slices and fewer retained data points. The purpose was to demonstrate whether the surface is still realistic when derived from much fewer points. Although some detail is lost, the surface does not contain any obvious errors. The time to fit the RBFs and extract the surface is substantially reduced with fewer data points. Table 1 lists the computation times for the experiments. A final scan of the duck was done to extract both the inner and outer surfaces. The toy duck is made from a constant thickness material, so the inner and outer surfaces are similar. The last rendering in Figure 8 shows a cross-section with both surfaces together. The fish phantom examination demonstrates the ability to combine sweeps from opposite viewpoints. Four sweeps are used to cover both the front and back of

the left and right portions of the tail. Joining the front and back views is the more difficult task because of a large gap. Despite the gap, Figure 9 shows a surface that combines the views seamlessly while retaining a high level of surface detail. The mesh also shows the quality of the triangulation; few low aspect triangles are created. The human hand was scanned to show the ability to extract surfaces of living tissue. The first examination with the fingers parallel was performed with three sweeps and is shown in Figure 10. This examination also demonstrates the effect of smoothing on the surface quality. Figure 10 shows that increasing the level of smoothing reduces the magnitude of the horizontal creases. Some influence of speckle is also reduced. The larger features, such as the crevices between fingers, are retained. The computing time for fitting the RBFs also increases. The second examination of the hand with fingers splayed is shown in Figure 11. Five sweeps are needed (one for each finger) but they are joined together with a small crease where a large gap exists. Horizontal creases from the misalignment of the slices are still visible at high resolution. This examination also demonstrates the effect of evaluating the surface with different resolutions. As the resolution decreases (sampling spacing increases), the small details of the hand are lost but the overall shape is retained. The time to create the surface also decreases (see Table 1), but the fitting time remains the same. The final examination is used to demonstrate the ability to form closed surfaces. The simulated cyst is scanned with a single sweep, but some dark shadows produce gaps in the surface, as shown in Figure 12. By using RBFs a second time, the mesh is repaired. Once an

Table 3. Comparison of measurements made of the fish phantom. The first measurements are from the constructed surface mesh. The second measurements are from calipers Measurement

First point (x,y,z) (mm)

Second point (x,y,z) (mm)

Distance (mm)

Calipers (mm)

Difference (mm)

Fin tip to waist Fin tip to notch Fin thickness (middle)

(⫺283.9,2209.8,33.3) (⫺284.0,2210.3,33.6) (⫺256.3,2220.2,39.2)

(⫺230.7,2219.7,34.3) (⫺257.5,2243.0,34.3) (⫺257.3,2217.8,30.25)

54.1 42.1 9.3

56.0 42.2 8.9 Average:

1.9 0.1 0.4 0.8

Surfaces from 3D ultrasound ● W. Y. ZHANG et al.

1471

Table 4. Comparison of measurements made of the closed human hand. The first measure-ments are from the constructed surface mesh. The second measurements are from calipers. Each of the measurements is labelled in Figure 10 Measurement (label)

First point (x,y,z) (mm)

Second point (x,y,z) (mm)

Distance (mm)

Calipers (mm)

Difference (mm)

Thumb to index (a) Index to middle (b) Middle to ring (c) Ring to pinkie (d) Top palm width (e) Mid palm width (f)

(⫺316.6,1756.2,⫺115.8) (⫺317.1,1693.8,⫺143.3) (⫺320.8,1688.2,⫺174.5) (⫺322.9,1699.0,⫺184.9) (⫺315.6,1771.3,⫺122.0) (⫺306.9,1809.9,⫺111.8)

(⫺321.5,1692.9,⫺127.4) (⫺319.4,1684.4,⫺145.6) (⫺317.5,1696.1,⫺172.1) (⫺317.3,1730.3,⫺182.7) (⫺322.1,1785.1,⫺188.9) (⫺317.6,1823.1,⫺188.8)

64.5 9.9 8.9 31.9 68.6 78.9

60.5 16.0 9.7 29.0 70.3 76.7 Average:

4.0 6.1 0.8 2.9 1.7 2.1 2.9

enclosed mesh is created, the bounding volume is computed. Using an isovalue of 35, a volume of 9.3 mL is computed. The actual volume is measured by submerging the cyst in a water-filled graduated cylinder and observing a displacement of 10.5 mL ⫾ 2 mL. The relatively large variation arises from the finite thickness of the cyst material, especially where it is partly folded together. By varying the isovalue from 25 to 45, the computed volume changes from 8.3 to 10.6 mL. This variation shows that the computed volume has some uncertainty related to the thickness of the cyst membrane, the resolution of the ultrasound machine and the mean intensity of the boundary. Nevertheless, there is an approximate agreement between the computed and known volumes of the object. The results of the surface validation measurements are listed in Tables 2, 3, 4 and 5. The measurements cover a variety of lengths and directions. The difference between the caliper measurements and the extracted surface arises from several sources of error, including the surface extraction, 3D ultrasound calibration, ultrasound resolution, landmark localization, tissue deformation and caliper accuracy. Several of these errors are known to be relatively large. Calibration errors range from 0.39 to 0.58 mm.

The axial resolution is approximated as 0.4 mm (the length of two cycles of the pulse). The lateral and elevational resolutions depend on focusing and depth, but are approximately 1 to 2 mm, respectively. Landmark localization error varies according to the object. For the duck phantom, the average difference between the caliper measurements and the extracted surface is 2.4 mm. The landmark localization for the duck phantom ranges from 1 to 2 mm. For the fish phantom, the average difference between the caliper measurements and the extracted surface is 0.8 mm. The landmark localization for the fish phantom was approximately 1 mm, since the features were more clearly defined. For the closed hand, the average difference was 2.9 mm and the closed hand was 2.1 mm. Tissue deformation and landmark localization were larger for the hand and ranged from 3 to 4 mm. The laser range finder produced an average difference of 2.5 mm, which is not significantly different than the 2.1 mm error from calipers. DISCUSSION In general, the surfaces are realistic and the algorithm handles gaps and multiple views without generating any obvious errors in the surface. The main drawback

Table 5. Comparison of measurements made of the open human hand. The first measure-ments are from the constructed surface mesh. The second measurements are from calipers. The third measurements are from laser surface scanning. Differences are measured with respect to the first mesh measurements. Each of the measurements is labelled in Figure 11 Measurement (label)

Distance (x,y,z) (mm)

Calipers (x,y,z) (mm)

Difference (mm)

Laser (mm)

Difference (mm)

Finger 1: index (g) Finger 2 (h) Finger 3 (i) Finger 4 (j) Thumb width (k) Finger 1: index width (l) Finger 2 width (m) Finger 3 width (n) Finger 4 width (o)

71.9 89.7 89.4 62.3 16.0 16.6 16.6 15.7 13.0

74.3 85.9 84.8 60.5 15.4 18.1 18.0 17.1 14.5 Average:

2.4 3.8 4.6 1.8 0.5 1.5 1.4 1.4 1.5 2.1

76.9 86.8 85.4 65.6 21.4 16.5 15.9 14.7 13.1 Average:

5.0 2.9 4.0 3.3 5.5 0.1 0.6 1.0 0.1 2.5

1472

Ultrasound in Medicine and Biology

is speed. Despite the use of center reduction and masking at the preprocessing step, the computing times are still too high for clinical use. Increasing the resolution and the level of smoothing increases the computation time further. Reducing the number of data points reduces the overall time while retaining a realistic surface. So it is reasonable to consider more elaborate schemes for reducing the number of data points. If the clinical goal is to estimate volume, then objects of these complexities can be represented with as few as 10,000 points without a significant loss in accuracy. Yet if the goal is accurate surface detail, then at least 100,000 points are required. An alternative solution to the problem of speed is simply to incorporate hardware acceleration of the calculations. Since the Ultrasonix machine already utilizes large field programmable gate arrays (FPGAs) it may be feasible to implement RBF fitting, evaluation and surface extraction on a single inexpensive FPGA. When computing times drop to the order of a few seconds, then clinical use is feasible. One of the other advantages of this method is the ability to extract surfaces without needing to adjust a number of parameters a priori. The main parameters to adjust are the fitting and evaluation accuracies, the level of smoothing, the sampling resolution, the amount of center reduction and the isovalue to extract the surface. Out of these parameters, only the isovalue should be under the control of the user. The others can be adjusted for a particular ultrasound system and remain constant. The geometry of the extracted surface matches the external measurements with an average difference ranging from 0.8 to 2.9 mm. These differences are within the range of errors from calibration, resolution and landmark localization. This gives confidence in the fidelity of the extracted surfaces. This paper focused on the problem of extracting bright organ surfaces bounded by darker fluid. RBFs were therefore applied directly to the 3D ultrasound data. For other organs, it may be possible to use RBFs if a suitable image processing step is incorporated that converts the 3D data to the type of problem described here. For example, an image processing step could be used to convert different textures into different intensities, such that a particular organ boundary is depicted at a different intensity than neighboring tissue with a different texture. In the future, it should be possible to incorporate more sophisticated smoothing techniques to further decrease the effect of speckle. Currently, isotropic smoothing is performed while fitting the RBFs. But it is possible to extract a smoother surface by simply changing the basis functions while evaluating the RBFs (Carr et al. 2003). Moreover, the basis functions can be tuned so that the smoothing is nonisotropic and nonuniform. This would be an advantage with freehand data since more

Volume 30, Number 11, 2004

smoothing should be done to cover the gaps and misalignment among the images that create creases in the images. Future work will also determine the ability to extract objects with very complex surfaces, such as malignant lesions. The ability to extract extremely complicated shapes has not been shown here, but more complicated shapes have been extracted with RBFs using other imaging modalities (Carr et al. 2001). One of the limitations with ultrasound is resolution; the object surface must be discernable in the ultrasound data to be extracted with RBFs. Ultrasound resolution, noise, calibration and other errors will create a limit on the complexity of the shape that can be extracted. Work is also needed to investigate the choice of parameters for particular clinical applications. Surfaces with reasonable fidelity were obtained with a range of different parameter values for the objects in this paper, but other organs may require more careful tuning. CONCLUSIONS This paper has demonstrated a new 3D ultrasound system with the ability to extract surfaces using radial basis functions. The interpolation properties of the functions allows multiple views of an object to be combined while accommodating overlap and large gaps in the data. This method produces realistic surfaces with a high level of detail. No reconstruction to a regular voxel array is required, so there is no inherent loss of resolution. Data smoothing is possible and can be controlled with a single parameter. The main limitation is the computational demand. With current personal computers, the method requires several minutes to extract a surface at high resolution. By combining a hardware implementation with more advanced data selection techniques to reduce the number of data points, surface extraction should be possible on the order of seconds. With improved speed, clinical studies can proceed with the goal of obtaining multiple surfaces of human organs. Acknowledgments—This work is supported by the Natural Sciences and Engineering Research Council of Canada and the British Columbia Advanced Systems Institute. The infrastructure of the author’s (Rohling) laboratory was developed under a grant from the Canadian Foundation for Innovation. The authors thank Qi Wei for her assistance with the experiments and for performing the laser scanning. The authors also thank the staff at Ultrasonix Medical Corporation for their assistance with the low-level access to the ultrasound machine.

REFERENCES Alyassin AM, Lancaster JL, Downs JH, Fox PT. Evaluation of new algorithms for the interactive measurement of surface area and volume. Med Phys 1994;21(6):741–752. Amenta N, Bern M, Kamvysselis M. A new voronoi-based surface reconstruction algorithm. Proc SIGGRAPH 1998:415– 421.

Surfaces from 3D ultrasound ● W. Y. ZHANG et al. Beatson RK, Greengard L. A short course on fast multipole methods. In: Ainsworth M, Levesley J, Light W, Marletta M, eds. Wavelets, multilevel methods and elliptic PDEs. Oxford, UK: Oxford University Press, 1997:1–37. Beatson RK, Cherrie JB, Ragozin DL. Fast evaluation of radial basis functions: methods for four-dimensional polyharmonic splines. SIAM J Math Anal 2001;32(6):1272–1310. Bernardini F, Mittleman J, Rushmeier H, Silva C, Taubin G. The ball-pivoting algorithm for surface reconstruction. IEEE Trans Vis Comp Graph 1995;5(4):349 –359. Boissonnat JD. Geometric structures for three-dimensional shape representation. ACM Transact Graph 1984:266 –286. Carr JC, Beatson RK, Cherrie JB. Reconstruction and representation of 3D objects with radial basis functions. Proc SIGGRAPH 2001:67–76. Carr JC, Beatson RK, McCallum BC. Smooth surface reconstruction from noisy range data. Proc ACM GRAPHITE 2003:119 –126. Chen CM, Lu HHS, Hsiao AT. A dual snake model of high penetrability for ultrasound image boundary extraction. Ultrasound Med Biol 2001;27(12):1651–1665. Cignoni P, Montani C, Scopigno R. A comparison of mesh simplification algorithms. Computers and Graphics 1998;22(1):37–54. Cohen L. On active contour models and balloons. Computer Vision Graphics and Image Processing. Image Understanding 1991;53(2): 211–218. Comeau RM, Sadikot AF, Fenster A, Peters TM. Intraoperative ultrasound for guidance and tissue shift correction in image-guided neurosurgery. Med Phys 2000;27(4):787– 800. Deng J, Gardener JE, Rodeck CH, Lees WR. Fetal echocardiography in three and four dimensions. Ultrasound Med Biol 1996;22:979 –986. Edelsbrunner H, Mucke EP. Three-dimensional alpha shape. ACM Transactions on Graphics 1994;13:43–72. Fenster A, Downey DB, Cardinal HN. Three-dimensional ultrasound imaging. Phys Med Biol 2001;46:R67–99. Fenster A, Surry K, Smith W, Gill J, Downey D. 3D ultrasound imaging: application in image guided therapy and biopsy. Computers and Graphics 2002;26(4):557–568. Greengard L, Rokhlin V. A fast algorithm for partical simulations. J Comput Phys 1987;73:325–348. Hoppe H, DeRose T, Duchamp T. Surface reconstruction from unorganized points. Proc SIGGRAPH 1993:19 –26. Lorensen WE, Cline HE. Marching Cubes: A high resolution 3-D surface reconstruction algorithm. Comput Graph 1987;21(4):163– 169. Leventon ME Faugeras O, Grimson WEL, Wells WM, III. Level set segmentation with intensity and curvature priors. Proceedings of the Workshop on Mathematical Methods in Biomedical Image Analysis, Hilton Head, SC, USA. June, 2000;4 –11. McInerney T, Terzopoulos D. Deformable models in medical image analysis: A survey. Med Image Anal 1996;1(2):91–108. Nelson T, Elvins T. Visualization of 3D ultrasound data, IEEE Comp Graph Appl 1993;13(6):50 –57. Nielson G, Hamann B. The asymptotic decider: Resolving the ambiguity of marching cubes. Proc IEEE Visualization 1991:83–91.

1473

Pagoulatos N, Haynor DR, Kim Y. A fast calibration method for 3D tracking of ultrasound images using a spatial localizer. Ultrasound Med Biol 2001;27(9):1219 –1229. Prager RW, Rohling RN, Gee AH, Berman L. Rapid calibration for 3D free-hand ultrasound. Ultrasound Med Biol 1998;24(6):855– 869. Prager R, Gee A, Treece G, Berman L. Freehand 3D ultrasound without voxels: Volume measurement and visualization using the Stradx system. Ultrasonics 2002;40(1/8):109 –115. Rohling R, Munger P, Hollerbach J, Peters T. Comparison of relative accuracy between a mechanical and an optical position tracker for image-guided neurosurgery. J Image Guided Surg 1995;1(1):30 –34. Rohling R, Gee A, Berman L. A comparison of freehand 3D ultrasound reconstruction techniques. Med Image Anal 1999;3(4):339 –359. Rohling R, Fung W, Lajevardi P. PUPIL: Programmable ultrasound platform and interface library. Proc MICCAI 2003. Lecture Notes in Computer Science 2879:424 – 431. Roxborough T, Nielson GM. Tetrahedron based, least squares, progressive volume models with application to freehand ultrasound data. Proc Visualization 2000:93–100. Sakas G, Walter S. Extracting surfaces from fuzzy 3D-ultrasound data. Proc SIGGRAPH 1995:465– 474. Schroeder W, Martin B, Lorenson B. The visualization toolkit: An objectoriented approach to 3D graphics, 3rd ed. New York: Kitware Inc., 2003. Sethian JA. Level set methods evolving interfaces in geometry, fluid mechanics, computer vision, and materials science. Cambridge: University Press, 1996. Singh A, Goldgof DB, Terzopoulos D, Eds. In: Deformable models in medical image analysis. Los Angeles: IEEE Computer Society Press, 1998. State A, Chen DT, Tector C. Case study: Observing a volume-rendered fetus within a pregnant patient. Proc IEEE Visualization 1994:364 – 368. Steen E, Olstad B. Volume rendering of 3D medical ultrasound data using direct feature mapping. IEEE Trans Med Imaging 1994; 13(3):517–525. Terzopoulos D, Witkin A, Kass M. constraints on deformable models: Recovering 3D shape and nonrigid motion. Aritific Intel 1988; 36(1):91–123. Treece GM, Prager RW, Gee AH, Berman L. Fast surface and volume estimation from non-parallel cross-sections for freehand 3-D ultrasound. Med Image Anal 1999a;3(2):141–173. Treece GM, Prager RW, Gee AH. Regularised marching tetrahedra: improved iso-surface extraction. Computers and Graphics 1999b; 23(4):583–598. Treece GM, Prager RW, Gee AH, Berman L. Surface interpolation from sparse cross-sections using region correspondence. IEEE Trans Med Imaging 2000;19(11):1106 –1114. Welch JN, Johnson JA, Bax MR, Badr R, Shahidi R. A real-time freehand 3D ultrasound system for image-guided surgery. Proc IEEE Ultrasonics Symposium 2000:1601–1604. Zhang Y, Rohling R, Pai DK. Direct surface extraction from 3D freehand ultrasound images. Proc IEEE Visualization 2002:45–52. Zhao HK, Osher S, Fedkiw R. Fast surface reconstruction using the level set method. IEEE Workshop on Variational and Level Set Methods in Computer Vision. (VLSM 2001) July 2001.