Estimating Volumes of Near-Spherical Molded Artifacts

The Food and Drug Administration (FDA) is conducting research on developing reference lung cancer lesions, called phantoms, to test computed tomography (CT) scanners and their software. FDA loaned two semi-spherical phantoms to the National Institute of Standards and Technology (NIST), called Green and Pink, and asked to have the phantoms’ volumes estimated. This report describes in detail both the metrology and computational methods used to estimate the phantoms’ volumes. Three sets of coordinate measuring machine (CMM) measured data were produced. One set of data involved reference surface data measurements of a known calibrated metal sphere. The other two sets were measurements of the two FDA phantoms at two densities, called the coarse set and the dense set. Two computational approaches were applied to the data. In the first approach spherical models were fit to the calibrated sphere data and to the phantom data. The second approach was to model the data points on the boundaries of the spheres with surface B-splines and then use the Divergence Theorem to estimate the volumes. Fitting a B-spline model to the calibrated sphere data was done as a reference check on the algorithm performance. It gave assurance that the volumes estimated for the phantoms would be meaningful. The results for the coarse and dense data sets tended to predict the volumes as expected and the results did show that the Green phantom was very near spherical. This was confirmed by both computational methods. The spherical model did not fit the Pink phantom as well and the B-spline approach provided a better estimate of the volume in that case.

The Food and Drug Administration (FDA) is conducting research on developing reference lung cancer lesions, called phantoms, to test computed tomography (CT) scanners and their software. FDA loaned two semi-spherical phantoms to the National Institute of Standards and Technology (NIST), called Green and Pink, and asked to have the phantoms' volumes estimated. This report describes in detail both the metrology and computational methods used to estimate the phantoms' volumes. Three sets of coordinate measuring machine (CMM) measured data were produced. One set of data involved reference surface data measurements of a known calibrated metal sphere. The other two sets were measurements of the two FDA phantoms at two densities, called the coarse set and the dense set. Two computational approaches were applied to the data. In the first approach spherical models were fit to the calibrated sphere data and to the phantom data. The second approach was to model the data points on the boundaries of the spheres with surface B-splines and then use the Divergence Theorem to estimate the volumes. Fitting a B-spline model to the calibrated sphere data was done as a reference check on the algorithm performance. It gave assurance that the volumes estimated for the phantoms would be meaningful. The results for the coarse and dense data sets tended to predict the volumes as expected and the results did show that the Green phantom was very near spherical. This was confirmed by both computational methods. The spherical model did not fit the Pink phantom as well and the B-spline approach provided a better estimate of the volume in that case. gradations to determine liquid displacement. However, in the case of these phantoms, the material used to manufacture these phantoms was porous and the phantoms would have to be coated. This, of course, would affect the "ground truth" volume estimate. As a result, the approach chosen for this study was based on a fundamental theorem in calculus, called the Divergence Theorem (see Taylor [2]), an analogue of Green's Theorem in two dimensional space. In the Divergence Theorem a volume integral is shown to be equal to a surface integral. Therefore, we surmised that, if a model of the surfaces of the phantoms could be developed, the Divergence Theorem would help to estimate their volumes. In order to develop a surface model we needed to obtain data about the surfaces of the phantoms. This was done using a coordinate measuring machine (CMM) in the Manufacturing Engineering Laboratory (MEL) at NIST. This machine produced a set of (x, y, z) points on the surface of each phantom. The data were then transformed to spherical coordinates and modeled using a set of basis functions, called B-splines. After fitting, the B-spline model was employed to generate a grid of values on the surface. These values were used to form surface triangles that were then used to compute the necessary surface integrals and finally the volume. The quality of the volume estimates depended on the surface grid sizes, as will be clear from the discussion below. The method of B-spline surface modeling is not new, in that it was suggested in the book by Dierckx [3]. However the application in the current case, that uses the Divergence Theorem, appears to be new. A reader can also consult the Dierckx references [4] and [5].
The paper is divided as follows. Section 2 describes how a volume can be computed using the Divergence Theorem. Section 3 describes the surface point generation experiments using the CMM. Section 4 introduces the two methods of data modeling used to estimate the phantoms' volumes. The first uses a spherical model. The second uses the B-splines as a basis for a least squares fit with regularization to model the phantoms' surface data. The fitted functions were then used to generate the grid data necessary to apply the Divergence Theorem. Section 5 presents the specific surface triangulation used to implement the Divergence Theorem. Section 6 describes the volumetric results from applying the spherical and B-spline models to the CMM data. A summary is given in Sec. 7.

Volume Estimation by the Divergence Theorem
In this section we state the Divergence Theorem and indicate how it can be used to estimate the volume of a polyhedron. This provides the motivation for the need to determine surface data points from the phantoms and then to build a surface model that is used to create data points at surface triangulation vertices. As the number of triangles was increased, the volume estimates were expected to tend to a fixed value. As will be seen, this expectation is confirmed below.

Divergence Theorem in 3-D and Volume Computation
Let Ω be a simply connected region in three dimensional space and Γ the surface boundary. Then (1) where F (x, y, z) = (F 1 (x, y, z), F 2 (x, y, z), F 3 (x, y, z)) is a differentiable vector field, (2) is the divergence of the vector field, n is an outwardpointing unit normal vector to Γ, and dσ is an infinitesimal element of the surface.

∫∫∫ ∫∫
Note that the left hand side is simply the volume of the region Ω. The 1/3 factor comes from the definition of F (x, y, z) so that ∇ → · F = 1. Now, if Γ is approximated by disjoint polyhedra (planar surface patches), Γ i , then (4) and (5) where nˆi is the unit outer normal to Γ i . Here we will model the surface patches by planar facets and, in particular, triangular facets. The plane for Γ i is given by nˆi ⋅ (x, y, z) = c i , where c i is a constant associated with each triangular facet. The sum of the integrals over the facets thus reduces to (6) This method of computing the volume of an object from a surface integral can be found in Schneider and Eberly [6].

Area of a Planar Polyhedron in 3-D
In order to estimate the surface area of a facet, as used in (6), we will describe the process of computing the area of a planar polyhedron in space by use of Stokes' Theorem and then we will particularize it to a planar triangle in space. Here the assumption is made that there is surface data available at the planar triangle vertices. Stokes' Theorem states that if C is a piecewise smooth boundary curve, oriented positively, of a surface Γ, and if F is a differentiable vector field defined on Γ and nˆis a unit normal satisfying the right-hand rule relative to the boundary orientation, then (7) where the curl of F (x, y, z) = (F 1 (x, y, z) , F 2 (x, y, z), F 3 (x, y, z)) is (8) The factor 1/2 comes from the definition of F (x, y, z) We will now specialize the Stokes formula to find the area of a spatial triangle, γ, in terms of its three vertices oriented positively. Let the three vertices, in positive orientation, be specified by v 1 = (x 1 , y 1 , z 1 ), v 2 = (x 2 , y 2 , z 2 ), v 3 = (x 3 , y 3 , z 3 ). Parameterize the boundary of the triangle as follows. For t ∈ [0, 1], ). The area of the spatial triangle, in terms of the parameterized boundary, can be written as  Institute of Standards and Technology   151   1   1 1 (  , , ) ( , , ) , n y t z t z t y t dt z t x t x t z t dt x t y t y t x t dt By straightforward integration over each of the parameterized segments it is easy to show that (11) The normal vector to the oriented triangle can be com-

Surface Metrology of the Artifacts
In this section we discuss the experimental method used to obtain surface data for the three objects used in this study. As a test object for the modeling process, a well calibrated sphere was selected. This object, along with the FDA phantoms formed the study artifact set. Data points, in the form of (x, y, z) coordinates, were created by the probing action of the CMM. Figure 2 shows the CMM system used to measure the artifacts. The system is computer controlled and touches an object to be measured at programmed points in order to produce (x, y, z) values at the probed points. Figure 3 shows one of the phantoms, held by a device called a vacuum chuck, as it is being touched by the CMM probe. The probe itself can be programmed to approach an object at various angles. In the background of the figure is a high quality reference steel sphere. Before an object is probed the reference sphere is measured to determine the effective diameter of the CMM probe tip. This reference sphere is separate from the calibrated sphere used to test the modeling process. The difference between the known diameter of the reference steel sphere and the apparent diameter of the measured reference metal sphere gave the calibrated effective probe diameter.  Institute of Standards and Technology   152   1 2  2 1  2 3  3 2  3 1  1 3   1 2  2 1  2 3  3 2  3 1  1 3   1 2  2 1  2 3  3 2  3 1  1 n y z y z y z y z y z y z 1   3  1  3  1  3  1  2   1  3  1  2  1  3  1  2  1  3  1  2   1  3  1  2  1  3  1  2  1  3  1   1  2  1  3  1  2 puted as follows. Let (  ,  ,  ) and   The FDA phantoms were created by a molding process and the mold marks on both spheres were visible. These mold marks were used to align the phantoms with the coordinate system of the CMM. The marks were laid out like lines of latitude, leading to the use of a natural nomenclature of latitude and longitude like those of the Earth. For each phantom, the North Pole was chosen to be the one with darker, deeper, or more obvious mold marks.
For the purpose of understanding the measurement process we will describe the physical fixture positioning of the phantoms. They were held by a vacuum chuck (see Fig. 3) to minimize distortion and reduce the chance of damaging the spheres' surfaces. As Fig. 3 shows, the chuck has a shallow cone to hold the phantoms. The phantoms contacted the cone around a circle of latitude at about -45°. This was measured from the equatorial circle around the middle of the phantoms. For example, the point at the top of the phantoms would be at + 90°. The wall vacuum of the cone was strong enough to hold the spheres sufficiently that they did not move significantly, as shown by the repeatable results from run to run at the 1 μm level. The setup of the phantoms in the vacuum chuck was accomplished by eye alignment using the visible mold marks and minor imperfections as guides to the eye. The expanded uncertainty of alignment was estimated to be approximately ± 2° for the vertical angle (i.e., keeping the equator horizontal) and ± 4° for the azimuthal angle. For a presentation of the guidelines for how expanded uncertainty of parameters is computed see Taylor and Kuyatt [7]. Essentially, the guideline for estimating the expanded uncertainty involves applying a multiplier, k = 2, times the square root of the variance about a predicted parameter value (i.e., two times the standard error). For a full discussion of confidence intervals for parameter estimation see Draper and Smith [8].
Three sets of measurements were planned for each of the phantoms. In each, points were probed on only a hemisphere at a time and later the data sets were postprocessed to form a spherical data set. In the first set of measurements the North Pole was set as the top point. A set of points was programmed to probe the top hemisphere down to the equator. The phantom was then re-positioned in the chuck so that the South pole was the top point. The same hemispherical points were probed. The rotation was done in such a manner that the mold marks representing the longitudinal lines were kept as aligned as possible. Post-processing of the data associated the correct signs with the measured coordinates relative to the CMM coordinate system. The third measurement involved re-positioning the phantoms so that the equatorial circle was vertical and the North-South axis was horizontal. Again two hemisphere sets of probe points were measured. This position was not feasible for the Pink phantom in one of the experiments described below.
Visually, the pink sphere was noticeably out-ofround, in the shape of an oblate spheroid. A dial caliper gave diameter measurements given in Table 1. The uncertainty of caliper measurements on hard steel surfaces is about 0.1 mm, and is estimated to be about 0.3 mm on the sample spheres due to the potential that the contact force would distort the soft surfaces of the spheres. All estimated uncertainties were k = 2 expanded uncertainties.
Two probing experiments were performed on each of the phantoms and the calibrated sphere. They created what we will call a coarse data set and a dense data set. For the coarse data set the plan was to measure each phantom on the CMM three times in each position, with 61 coordinate points per hemisphere. For the green sphere, each measurement set consisted of three separate sets of points: North pole up, South pole up, and prime meridian/equator intersection up. The pink sphere could not be held sideways in the first experiment, as the out-of-roundness prevented an effective vacuum seal. Therefore, a measurement data set for this sphere had only the North Pole up and South Pole up data. The third data set in this case was a re-measure of the North Pole up position. Each data set consisted of 122 points. The plots in Fig. 4 and 5 show the radial deviation from a best-fit sphere for the full data sets. The figures show the radial residuals obtained by fitting sphere models to the Green and Pink data with an algorithm ordinarily used during sphere calibration work. In particular, they represent the residual errors between the distance from the fitted sphere center to the probed points and the fitted radius of the sphere model. The residuals, in the case of the Green phantom, range from approximately -0.1 mm to + 0.1 mm, whereas the residuals, in the case of the Pink phantom, range from approximately -0.57 mm to + 0.23 mm. This suggests the slight non-spherical nature of the Pink phantom. Figure 6 shows the typical distribution of the probe points on a sphere. The plot is a transparency so  that probe points on the opposite side of the sphere are visible. The calibrated sphere on the shaft, with a diameter of 19.05 mm, was also measured with 122 points, repeated five times. It was measured in one position, since it was permanently mounted on a support shaft. For the dense data set 181 points were taken on the phantoms and the calibrated sphere, with five repeats in each position. The positions were taken the same as those in the first experiment. That is, the alignment was selected with North pole up (position 1), South pole up (position 2), and prime meridian/equator intersection up (position 3). In this case it was possible to hold the Pink phantom in the sideways position 3. The calibrated sphere was also measured at 181 points with five repeats.

Modeling Methodologies
In this section two forms of data modeling will be discussed. Since the phantoms seemed to be nearly spherical the natural tendency was first to consider fitting a spherical model to each phantom and estimating the volume of the fitted spheres. However, in order to develop a potentially more accurate volume estimation model, the surface data was also fit using tensor products of B-splines and the volumes estimated by the Divergence Theorem.

A Spherical Model
In this section we consider how close the data could be modeled by a spherical model for the data. The calibrated sphere, of course, could be modeled via a spherical model.
In particular, let c = (c 1 , c 2 , c 3 , c 4 ) and set (12) The unknown parameters c 1 , c 2 , c 3 represent the center of the sphere and c 4 is the radius. All of the data were measured in millimeters so that the parameters naturally have millimeter units. Define (13) where O, the objective function, is a measure of the residual for the fitted sphere defined by the coefficients c. Since the function O is a nonlinear implicit function of the parameters we needed to use a nonlinear minimization algorithm to find the best fit, i.e., to solve the problem There are a number of algorithms for fitting least squares models to data on geometric shapes (see Shakarji [9]). There are also various algorithms for minimizing general nonlinear functions, such as (13). All of the algorithms involve iterative minimization of some form. Many require computing derivatives of the objective function in order to generate search directions along which to identify a minimum. Others do not involve derivatives but may be somewhat slower in the minimization search. The algorithm selected here, because of the relatively few parameters involved, and the fact that derivatives are not required, is a form of polyhedron search method called the Nelder-Mead method (see Sauer [10]). The Newton method, requiring derivatives, was initially used to estimate the parameters, but the Nelder-Mead tended to produce the smallest value to (13).
The Nelder-Mead algorithm works iteratively through steps that involve reflections, vertex extensions, and multidimensional polyhedron contractions until the volume of the polyhedron becomes less than a prescribed tolerance. The polyhedron vertices then enclose the minimum. The median value of the polyhedron vertex values is taken as the minimum value.
The uncertainties of the estimated center and radius were computed using the methods proposed in Draper and Smith [8] for nonlinear regression. In particular, if c = (ĉ 1 ,ĉ 2 ,ĉ 3 ,ĉ 4 ), then define the n × 4 matrix with elements (15) The i-th row is given by (16) The standard error of ĉ i , s.e. (ĉ i ), is given by the square root of the i-th diagonal element of (Ẑ' Ẑ ) -1 s 2 where s = O (ĉ) / (n -4). The expanded uncertainty is given by 2s.e.(ĉ i ) as defined in Taylor and Kuyatt [20]. The uncertainty limits ofĉ i are ĉ i ± 2s.e. (ĉ i ). All units are in millimeters except volume, which is in cubic millimeters.

A B-spline Model
Given measured data points on the surfaces of the calibrated sphere and the two phantoms, surface models can be constructed using B-splines as basis functions. In this section we will define cubic B-splines and show the construction of tensor products of B-splines.

Cubic B-Splines in One Variable
Suppose that a function y = f (x), x ∈ [a, b], is known at the n points (x 1 , y 1 ) ,···, (x n , y n ), where a < x 1 < x 2 < ··· < x n < b, y q = f (x q ), q = 1,···, n. a and b are finite interval bounds. It is known that a polynomial of degree n -1, P (x), can be constructed to pass through these n points. In the case of highly accurate data points this polynomial can be constructed to interpolate these n points by, for example, Lagrange polynomial or Newton divided difference algorithms. But, for large n, it is also well known that polynomials of high degree can produce unwanted oscillations between the interpolated points. It is crucial then to approximate sets of data with as low degree polynomials as possible. Of course these polynomials may or may not interpolate the data points but may be made to come as close to them as possible. The ability to create highly flexible approximants from low-degree polynomials is a significant advantage of functions called splines.
A well known mathematical technique to construct complicated functions is to form linear combinations of simpler functions, called basis functions. In the current paper, data sets will be approximated by linear combinations of special spline functions, called B-splines, for Basis splines. It is known that any spline function can be represented in terms of B-splines. This particular basis has the advantage of leading to computational algorithms that are elegant, efficient, and stable. Only B-splines of order four will be considered here. They are cubic splines that are non-zero only over four adjacent intervals between knots (see Fig. 7). The notation for a B-spline is where in this work p = 4. To simplify notation let N i (x) = N 4, i (x). Then a B-spline of order four, or cubic B-spline, is a cubic spline with knots ξ i -4 , ξ i -3 , ξ i -2 , ξ i -1 , ξ i that is zero everywhere except in the range ξ i -4 < x < ξ i . N i (x) is defined uniquely except for a scaling factor and is conventionally taken to be positive throughout the range ξ i -4 < x < ξ i and has a single maximum value. Since N i (x) is a cubic spline it has continuous derivatives of order one and two at ξ i -4 and ξ i . These derivatives are zero at the endpoint knots.
There are computational advantages in using cubic B-splines. For any given x, all but four adjacent A least squares curve fitting problem to the data a < x 1 < x 2 < ··· < x n < b, y q = f (x q ), q = 1,···, n , becomes one of determining the coefficients c i as the least squares solution to the equations (18) These may be written in matrix notation as (19) where A is the n × r matrix whose element in column i of row q is N i (x q ) and c, y are column vectors with elements c i , y q , q = 1, ···, n , respectively. If the data points are arranged in increasing order of x, then the matrix A becomes a banded matrix with bandwidth four. For a more thorough discussion of B-splines and their computation see de Boor [11] and Cox [12].
There are two functions in the MATLAB Spline Toolbox that implement the evaluation of B-splines. A set of knots can be augmented at the ends by the function augknt and the B-splines and their derivatives can be evaluated by the function spcol.

Tensor Products of Cubic B-Splines in Two Variables
In this section the B-spline concept will be extended to two dimensions in order to fit two dimensional scattered data by a surface function. In this surfacefitting problem data points (x q , y q ), q = 1, 2, ···, n, and values at these points, z q = f (x q , y q ), q = 1, 2, ···, n, are given. The surface model used to fit these data points will involve sums of products of B-splines.
To introduce this model, define a rectangle, say R, by

Ac y ≈
does not restrict the data points to the Euclidean plane. They could, for example be angular coordinates, as will be seen later. The rectangle is subdivided by sets of where rˆ, sˆare indices, not necessarily equal and a, b, c, d are the bounds on the rectangle R. The knots are then extended by eight in each dimension as done in the one dimensional case. These knots divide the rectangle R into rectangular panels in the plane given by R i j , i = 1, 2, ···, rˆ, j = 1, 2, ···, sˆ. Then, a basis set of splines for this pairing of knots can be constructed as products of B-splines N i (x) N j (y). In fact, the surface spline model is given by (20) called a tensor product of splines, where r = rˆ+ 4 and s = sˆ+ 4. As in the one dimensional case, these tensor product splines have a number of advantages. First of all, the basis functions N i (x) N j (y) are each non-zero over a rectangle composed of sixteen panels in a 4 × 4 arrangement. In particular N i (x) N j (y) is non-zero only when ξ i-4 ≤ x ≤ ξ i and η j -4 ≤ y ≤ η j . Next, if (x q , y q ), q = 1, 2, ···, n and values at these points, z q = f (x q , y q ), q = 1, 2, ···, n are given, then the fitting problem can be formulated as finding the least squares solution of (21) Again, this can be written in a matrix form as (22) where A is now a matrix with n rows and rs columns, and z is a column vector with n rows. The elements in row q, q = 1, 2, ···, n, of A are formed as follows. We start with a fixed j, j = 1, 2, ···, s, and let i, i = 1, 2, ···, r vary for that j. Then, the element in column (j -1)r + i for row q is given by N i (x q ) N j (y q ) . The j is then incremented and the i 's varied again. The end resulting column values in A for row q would look like . For a full discussion of tensor products of B-splines see de Boor [11], Eberly [13], and Rogers and Adams [14].

An Issue in Computing Tensor Product B-Splines and the Relation to Least Squares
Unfortunately, for some choices of knots the resulting matrix A in (22) might be rank deficient and it would not be possible to use the normal equations to solve the least squares problem. This can potentially happen in the case of widely scattered data, where some of the knot panels do not contain data points. This problem can be solved, though, by using the matrix singular value decomposition.
Assume that a tensor product spline has been computed as described in Sec. 4.2.2. A least squares fit can be done to the scattered data as follows. Since the matrix A could be rank deficient, with potentially zero rows or columns, we cannot rely on the standard normal equations for determining the coefficients, but using the matrix singular value decomposition provides a convenient substitute (see Golub and Van Loan [15]). In the singular value decomposition of the matrix A, with the number of rows larger than the number of columns, the matrix is decomposed into the product of a column-orthogonal matrix U, a diagonal matrix S, with diagonal elements S i , and the transpose V' of a square orthogonal matrix, so that A = USV'. In order to solve the problem Ac = z we compute c = V(S + (U' z)), where (23) defines the generalized inverse of S and U' is the transpose of U. A tolerance is used to determine which of the small singular values in the decomposition should be considered zero.

A Knot Selection Algorithm
In a least squares data fitting process involving B-spline basis functions, the resulting model residuals are sensitive to the knot placement for the B-splines. The selection of B-spline knots in order to achieve as small a residual as possible during the least squares process is a nontrivial task. One would be extremely lucky to manually select a set of knots that could achieve a very small least squares residual. In this section we will describe a heuristic algorithm that, in practice, generates a set of knots that produces small , , is that thesum, at the -th iteration of Resid , i s ta- ken over all pair of points ( , ) such that .
k q y υ η + < least squares residuals. The strategy involves an iterative knot insertion algorithm. An initial set of knots is selected by the algorithm user and, at each iteration of the algorithm, new knots are inserted in the vicinity of the previous fit where the local residuals are the largest. The knots are not moved once they have been inserted. The iterative algorithm has a stopping criteria based on a statistical test.
First we will discuss the knot insertion algorithm. It is based on a strategy suggested by Dierckx [3]. At the beginning of each iteration the assumption is that there exists a current set of knots. In the first iteration of the algorithm these would be an initial set chosen by the user. The tensor product of the B-splines is computed at the data points, a least squares fit is made to the data, and the current absolute residuals of the fit are computed at each data point. The current knots divide the (x, y) plane into rectangles. Some rectangles have data points and others don't. Let the knots at the k-th iteration be labeled a < ξ 1 (k) < ξ 2 (k) < ··· < ξ r k (k) < b along the x-axis and c < η 1 (k) < η 2 (k) < ··· < η s k (k) < d along the y-axis. The (k + 1)-iteration begins by associating all of the data points with the knot panels in which they fall. Let the index pair ij indicate the R i j -th panel defined by ξ i ≤ x ≤ ξ i+1 , η j ≤ y ≤ η j +1 . Then suppose the data values (x 1 (ij ) , y 1 (ij ) ), ··· , (x ri j (ij) , y r ij (ij) ), fall into the R ij -th panel. Let F k (x, y) be the least squares B-spline function fit to the data at the k-th iteration. We next compute the residuals of the fit at all of the data points (24) From these residuals we form the sums of squares of the residuals that fall within the knot intervals. The sums are separated into the x direction and the y direction as follows.
What is meant here, for example in the case of We next find the maximums of these sums of squared residuals.
where u = i for some i = 1,2,···, r k and υ = j for some j = 1,2,···, s k . The next step is to add one knot at a time at each iteration. In particular, a knot is added in the x upon a weighted average of x or y data points in the columns or rows determined by the knot intervals with the maximum residual errors determined by (26). In particular, the positions are given by (27) We note that (28) so that (27) represents weighted averages of all of the The knots are then reindexed as necessary and F k+1 (x, y) is computed as the least squares B-spline function fit to the data at the (k + 1)-iteration based on the new knot set. The iterations continue until the stopping criteria is met.
The stopping criteria for the k-th iteration used in this algorithm is based on the use of the R 2 -statistic, called the coefficient of multiple determination (see Draper and Smith [8]). R 2 is the square of the correlation Resid ( , ), 1, , .
direction, if or added in the direction, if . The knot added is positioned based data pairs ( , ) satisfying and between the vector of observed data, z q , q = 1,···, n, and the least squares predicted data, F k (x q , y q ), q = 1,···, n. The statistic satisfies 0 ≤ R 2 ≤ 1. This statistic is often used as a measure of how well the regression equation explains the variation in the data. It is known that in building models based on adding terms in the regression equation, care must be taken in using this statistic. However, in the current algorithm, the statistic is used in a somewhat non-conventional manner. We use it as a measure of the benefit of adding more knots to the tensor product spline model. The knot selection algorithm is terminated when R 2 > 0.98.
Although the knot placement algorithm is heuristic, coupling it with the R 2 -statistic has shown good convergence in practice. It is reasonable to expect this, since knots are placed in intervals in which the fits at a previous iteration showed the largest local error. Placing a knot there allows extra flexibility for those areas.

Data Smoothing by Tikhonov Regularization
As noted in Sec. 4.2.3, the data distribution can lead to a rank deficient least squares matrix. Even though a fit can be computed using SVD, evaluations of the fitted function at some points can lead to unwanted oscillations. It is necessary to introduce a smoothing procedure by regularization. Regularization is a way of introducing extra information to the least squares objective function in order to control overfitting of the data. It is the overfitting of data that can lead to the unwanted oscillations. In the present work a penalty term is introduced to control the smoothness of the resulting fitted function. The objective function will then balance the data fitting with the smoothness of the fit. The second partial derivatives of the tensor product B-splines will be introduced as the smoothing operators. The objective function will take the form (29) λ is called a Tikhonov parameter.
In the objective function the data values are given by (x q , y q , z q ), q = 1, ···, n. The smoothing terms will be evaluated at a new set of points chosen so that every knot panel has an equal number of points assigned to the panel. In the current case there will be a total of t points throughout the knot panels given by (u p , υ p ), p = 1, ···, t .
To write this in matrix form we will define two new matrices B 1 and B 2 as follows. For m = 1,2,···, rs let (i m , j m ) be such that m = ( j m -1) r + i m . Then define the matrix elements where Q is (n + 2t) × (n + 2t) and R is (n + 2t) × rs. Q is orthogonal, so that Q T Q = I, and R is zero except in its
At this point we need to discuss the selection of the Tikhonov parameter, λ. A number of methods for choosing the parameter have been discussed in the literature (see [16,17,18,19,20]). For this work, though, we used a graphical method that led to the selection of a parameter in a few iterations. Before continuing to discuss the selection of the Tikhonov parameter for the current study we need to change the coordinates of the original probe points to spherical coordinates defined on a rectangle.

A Coordinate Transformation
Since our tensor product B-spline requires a rectangular grid, we convert our CMM data to spherical coordinates. We start with a given set of n points, (x i , y i , z i ), i = 1, 2, ···, n, for some n, on a surface. As measured, these points are given relative to the CMM origin. The first step is to center the data by using the center-of-data mass of the data points. This establishes an origin at the center-of-data mass point. It is done in order to establish a common reference point interior to the measured data points. It also simplifies writing vectors from the origin to the data points and allows the introduction of spherical coordinates. The centerof-data point is given by Since the data set is enclosed in a near-spherical bounded region, it is reasonable to identify the Euclidean data points with spherical coordinates. This will map the points on the sphere to a rectangular surface where the surface coordinates are designated by θ, φ and the height of a surface point is given by  R (φ , θ ). In order to use spherical coordinates to represent points on the boundary of a surface, we need to restrict our analysis to surfaces that are called star-shaped. These are surfaces in which a ray drawn from the center-of-data mass intersects the boundary in a unique point. Whereas, in the definition of the B-splines, we used the coordinates (x, y) we will now use (φ , θ ) and build B-splines in terms of these spherical coordinates. Euclidean coordinates will now refer to the measured data points. This switch in notation should not, it is hoped, cause too much confusion.
As an illustration, the results of the conversion of the coarse probe points for the Green phantom from Euclidean coordinates to φ, θ coordinates are shown in Fig. 9. We note the density of data points near the equator is higher than towards the poles at θ = 0 and θ = π. Unfortunately the distribution of data points was dictated by the software controlling the CMM. The lack of data points near the poles leads to a well known problem, called the Pole Problem in the literature (see Dierckx [3]). It creates rank deficient matrices during the least squares fitting process. Section 4.2.3 discussed how the singular value decomposition can be used to handle this problem.

Choosing a Tikhonov Parameter
Once the original probe data had been converted to spherical coordinates, the iterative selection of the Tikhonov parameter for the current volume estimation problem proceeded as follows. First, a set of knots was selected with λ = 0 (i.e., no regularization terms) in order to produce an R 2 > 0.98 as described in Sec. 4.3. We found that, for all of the data sets examined, only a very few extra knots were added beyond the initial set used to begin the knot selection process. This led to a rapid convergence to a set of final knots in all cases. These knots formed panels in the rectangle [0, π] × [0, 2π]. Next, nine points were uniformly chosen in each panel and the regularization terms formed. Numerical experimentation showed that the use of nine points in each panel provided sufficient extra data to the regularization terms in order to smooth the final fits. An initial Tikhonov parameter, λ, was chosen and the objective function (29) was minimized by the QR method discussed above. An initial grid of 40 θ and 80 φ values was generated in the rectangle [0, π] × [0, 2π]. This grid was finally fixed on for selecting a Tikhonov parameter since further numerical experimentation showed that the response surface of the radii for denser grids did not change the final range of the radii. The coefficients, c, computed to minimize the objective function (29), were then inserted into the linear form Ac, where A was the tensor product matrix formed from all of the 40 × 80 grid points. The end resulting radii at the grid values were then computed and the spherical volumes for those radii were computed. The maximum and minimum sphere volumes over the entire 40 × 80 grid were determined. To select the appropriate Tikhonov parameter, λ values over a range were used and the differences between the maximum and minimum volumes were plotted. The value of λ that produced the minimum difference was selected as the working λ. For the current volume estimation problem, the final λ was λ = 0.32.
To illustrate the effect of using regularization terms to smooth the least squares fit of the radii see Fig. 10.

Octant Number
Cartesian Coordinate Spherical Coordinate This figure shows the radii data for the case of λ = 0, or no regularization. Note the large radii oscillations at the boundaries (a range of approximately 10,000 mm). Now see Fig. 11 and note the narrow range (approximately 0.1 mm) of the radii over the entire grid. This clearly shows the effect of the regularization terms on the least squares fit. The data used for these plots was the coarse Green phantom data.

Volume Estimation by Surface Triangulation
We could now estimate the volume using the triangulation method described in Sec. 2 by dividing the phantom surfaces into triangular surface patches as follows. First of all we partitioned the phantom surfaces at grid points located at the colatitude angles θ 1 = 0 < θ 2 < ···< θ υ < θ υ + 1 = π from the north pole to the south pole and azimuthal angles φ 1 = 0 < φ 2 < ··· < φ h < φ h + 1 = 2π around the phantom surface, where υ stands for vertical and h for horizontal. Since these grid points did not necessarily fall at the measured data points, a radius, R(φ, θ ), was calculated at the grid points from the regularized fitted B-spline surface model. At the north and south poles the radius was taken as the median value, R north , of the grid point values R(φ i , 0), i = 1,···, h, for the north pole and R south , the median value of R(φ i , π), i = 1,···, h, . The spherical coordinates of all of the grid points were converted to Euclidean coordinates on the surface by (35) At this point we could apply the Divergence Theorem method of Sec. 2. The patches at the north poles were easily constructed to be triangular as part of the process of determining the contribution of each patch to the phantom volumes. In particular, at the north pole the designated point was (x 1 , y 1 , z 1 ) = (0, 0, R north ). We then iterated through the spherical coordinate points (φ i , θ 2 ), i = 1,···, h. At each of these angle pairs there was a value R(φ i , θ 2 ), i = 1,···, h. We generated the volume by adding up the contributions of each patch to the volume total. We did this by initializing a variable, vol, for the volume, to zero. We then started to generate the contribution from the first layer of patches at the north pole. As noted above, we set the north pole to (x 1 , y 1 , z 1 ) and then set ( , )cos( ). z R φ θ θ = These were set in a positive orientation.
In order to make this process more concrete, we have included a sample surface triangulation with υ = 5 colatitude angles and h = 8 azimuthal angles in Fig. 12. In this figure the Euclidean grid points have been indexed. The North Pole (x 1 , y 1 , z 1 ) is designated by the index 1. In the first step described above the points (x 2 , y 2 , z 2 ) and (x 3 , y 3 , z 3 ) are indexed in Fig. 12 by points 2 and 3 respectively. Next we computed the outward normal to the triangle patch by forming the vectors υ 1 = (x 2 , y 2 , z 2 ) -(x 1 , y 1 , z 1 ), υ 2 = (x 3 , y 3 , z 3 )-(x 1 , y 1 , z 1 ), for triangle 123 in Fig. 12 and then formed the normalized cross product (36) We then computed the contribution that this patch made to the volume as where (38) is set in order to compensate for the cyclical vertex indexing around the triangle patch. Formula (37) is a combination of Eq. (6), where (n · (x 1 , y 1 , z 1 )) is the coefficient c in (6), and the second factor is a compact form of Eq. (11). The factors 1/3 from (6) and 1/2 from (11) were combined as a multiple of 1/6 after all of the summations had been performed. x y x y We proceeded to the next patch in the North Pole layer. Again (x 1 , y 1 , z 1 ) was the North Pole, indexed by 1 in Fig. 12, and we then used the previous computation to get x 2 = x 3 , y 2 = y 3 , z 2 = z 3 and set In Fig. 12 the new point (x 3 , y 3 , z 3 ) is indexed by 4 in Fig. 12. The triangle of interest is now 134 in terms of indices. We computed the normalized cross product as for the first patch and then computed the contribution of the second patch to the volume using Eq. (37). We continued this process for θ 2 , φ i , i = 1,···, h. In Fig. 12 we would have proceeded with computing contributions to the volume by working through the indexed triangles 123, 134, 145, 156, 167, 178, 189, and 192.
Finally, at the South Pole there were h triangles to include in the volume calculation. Their vertices were identified as follows x 2 = 0, y 2 = 0, z 2 = -R south .
As an example, in Fig. 13 one of these triangles would be indexed by the points 23 26 24, where the South Pole is point number 26, although the South Pole index may be hard to see in the figure. The contribution of these triangles to the volume was computed as above. After all of the triangle contributions to the volume were computed, the final volume was taken to be vol = (1/6) vol. The factor 1/6 comes from the product of 1/3 from Sec. 2.1 and 1/2 from Sec. 2.2. The reader can refer back to these sections to see how the factors arose. The triangulation process can clearly be generalized to a denser surface triangulation.  Institute of Standards and Technology   166   3  3  2  2  3 ( , )sin( )cos( ), ( , )cos( ). z R φ θ θ =

Computational Results
Two computational processes were involved in generating the results for each object. The objects involved were the calibrated sphere and the two phantoms. First, sphere models were fit to an object. The second process involved three steps. The first step was to compute a good selection of knots for the tensor product spline function. The next step was to fit the tensor product spline function with the objective function modified by the regularization terms to the object. Finally, once the tensor spline function had been developed they were used to generate data at the vertices of the surface triangulation and to compute the volume using the Divergence Theorem.

Computational Results for the Spherical Model
Tables 3, 4, and 5 present the results obtained by fitting spherical models to the data from the calibrated sphere and the two phantoms. For the calibrated sphere there were five repeats for each of the coarse and dense data sets. Since the results of the point location repeats differed only at the micrometer level the data sets were combined by averaging to form two data sets representing the coarse and dense data sets for the calibrated sphere. Similarly, the point locations of the three position data sets for both the coarse and dense data sets for the phantoms differed only at the micrometer level. Therefore, by averaging of the three position data sets for coarse and dense data sets, two working data sets were formed for the Green and Pink phantoms. Since the output of the Nelder-Mead algorithm was a final polytope surrounding the minimizing parameter with the polytope, the final reported parameters were selected as the median value of the vertex values. These are the first four entries in the tables: Center x, Center y, Center z, and Radius. The units are millimeters. The fifth table entries are the spherical volumes based on the median radius value in cubic millimeters. The radius residuals were computed as the absolute value of  The sixth and seventh entries in the tables give the mean value and standard deviation of these residuals in millimeters. The eighth through the eleventh entries are the expanded uncertainties for the estimated Center x, Center y, Center z, and Radius, also in millimeters.
The results in these sections on spherical model fitting involved the non-linear model (39) fitting to the sparse data sets generated by the CMM. This process did not involve the B-spline algorithm. Therefore any differences between coarse and dense results were most likely the consequence of the point selection in the metrology of the artifacts. Table 3 gives the results of the fit of the sphere model to both the coarse (122 points) and the dense (181 points) data sets for the calibrated sphere. In both cases the average of the five repeat data sets was used for the fitting process. As can be seen, the two volume estimates differ by approximately 0.04 %. The radii estimates differ by about 0.01 %. There is a slightly larger difference when the spherical fit results are compared to the volume estimate based on the manufacturer estimated calibrated sphere diameter of 19.05 mm. This would lead to a volume estimate of 3619.8 mm 3 . Since the method of estimating the diameter by the manufacturer was proprietary we could not independently verify the measurement. We could only compare it against our sphere model fit results. The manufacturer measured radius differs from the computed radii in Table 3 by about 0.08 % whereas the volume estimate based on the measured radius differs from the volume estimates in Table 3 by about 0.2 %. If we expand the radii in Table 3 to a diameter we find the values 19.0651 millimeters in the coarse case and 19.0629 millimeters in the dense case. These differ by approximately 0.01 mm from the manufacture's measured diameter. It would seem that these differences were sufficiently small so that the difference in value estimates for radius and volume were not considered significant. Although the estimated centers are slightly different, all other values are of the same order of magnitude. Based on these results we can accept that the CMM is producing accurate position data for spherical metallic artifacts. Tables 4 and 5, however, begin to show the consequences involved with attempting to measure slightly non-spherical and non-metallic artifacts with the CMM.

Results for the Spherical Model Fit to the Phantoms With Coarse and Dense Data Sets
From Tables 4 and 5 it is clear that the Green phantom is more spherical than the Pink phantom as indicated by the residuals and the expanded uncertainties. This simply confirms the fact that the Pink phantom was more difficult to measure using the vacuum chuck due to its lack of sphericity. For example, the Mean Radial Residual for the Pink data is two orders of magnitude larger than for the Green data. The Standard Deviations for the Pink data are an order of magnitude greater, as are the Expanded Uncertainties. It is not clear how much the phantom material affected the results since it was difficult to set up a specific probe test for metallic versus phantom material. The artifacts would have to have been exactly the same size, positioned at exactly the right location, and probed at exactly the same coordinates to separate those factors from the material factor difference. There were no such comparable artifacts. We can make a guess, though, that there might be some effect due to probe force against the non-metallic material if we look at Table 3 and Table 4 for the Green phantom, the most spherical of the two phantoms. The Expanded Uncertainties for the sphere fit to the calibrated metallic sphere and the Expanded Uncertainties for the sphere fit to the Green phantom data differ by one to three orders of magnitude. Since the repeatability of the CMM measurements is at the 1 μm level, it is likely then that material difference had some significant affect on the difference in the uncertainties. It is a conjecture that this and the non-spherical shape of the Pink phantom account for a large part of the differences between the Pink phantom Expanded Uncertainties in Table 5, the Green phantom Expanded Uncertainties in Table 4, and the calibrated sphere Expanded Uncertainties in Table 3. For the Green phantom the radii estimates in Table 4 between the fits of the coarse and dense data sets is approximately 0.3 % and the volumes differ by approximately 0.8 %. Whereas, for the Pink phantom the radii estimates in Table 5 between the fits of the coarse and dense data sets are 0.4 % and the volumes differ by approximately 1.3 %. In both Tables 4 and 5 the expanded uncertainties for the dense data sets are approximately an order of magnitude better than the results for the coarse data sets. This suggests that the radius and volume estimates in the case of the dense data sets for both the Green and Pink phantoms are probably the more realistic values, at least in terms of a spherical model fit.

Spherical Fit Results to Pink Phantom
Properties Coarse Dense Table 1 shows that the diameter of the Green phantom is approximately 20 mm and that the diameter of the Pink phantom falls approximately between 20.0 mm and 18.4 mm. This would imply that the sphere model for the Green phantom should have a volume of about 4188.8 mm 3 .With the estimate of 0.3 mm uncertainty using the calipers, this would place the Green volume in the range of 4380.1 mm 3 and 4003.1 mm 3 . The estimated volumes in Table 4 are within this range for both the coarse and the dense data. The Pink phantom should have a volume between 4445.2 mm 3 and 3104.8 mm 3 . Table 5 for the Pink phantom shows that for the coarse and dense data the volume estimates fall from about 3862 mm 3 to about 3913 mm 3 . These volumes are within the expected range of the caliper measurements.

Computational Results for the B-Spline Model
In this section we will describe the method used to estimate volume uncertainties and discuss the results of using a B-spline surface model and the Divergence Theorem to estimate the phantom volumes. A nonparametric method to estimate the volume uncertainties was chosen since there did not exist any "ground truth" values for the Green and Pink FDA phantom volumes. All of the expanded uncertainty results are displayed in row two of Tables 6, 7, and 8.
For each data set we conducted twenty one estimates of the volumes and volume uncertainties by increasing the density of the surface triangulation. The twenty one were selected since beyond that point the computer time became large and the volume increment per case was minimal. In order to plot the volume results in terms of grid density we use the term grid size by which we mean the product of the number of θ values times the number of φ values for a particular grid density. As shown in Tables 6 through 8 we begin with a grid of ten θ and twenty φ values and compute the volume and uncertainty. The grid size in this 10 × 20 grid is then 200. The θ values were then incremented by ten and the φ values by twenty for each grid case until the last case of 210 θ values and 420 φ values, giving a grid size of 88200 (8.82 × 10 4 ). The volumes in Figs. 14 through 19 grow very rapidly and appear to approach a fixed value as the grid size increases. They simply reflect the volume data in the Tables 6 through 8.

Estimating Volume Uncertainties
The uncertainties were estimated by the nonparametric "bootstrap" method. It is a computer intensive technique for estimating uncertainties and it involves repeated Monte Carlo resampling from the spherical coordinates of the original measured data sets with radii values modified by the fitting residuals, refitting the model to estimate new volumes, and finally computing an uncertainty for the process from the set of computed volumes. For a full discussion of the bootstrap see Efron and Tibshirani [21] but we will give a brief description here of how we applied the bootstrap method in the current study.
The object of our application of the bootstrap was to develop volume uncertainty as a function of grid size. The process began with the conversion of the Euclidean data points to spherical coordinates. An automatic selection of knots was done. These steps were done once for a given data set. It was assumed that the modifications of the data made during the bootstrap would be small and not affect the knot selection. During the first pass of the bootstrap algorithm the radii of the spherical coordinates were fit by a tensor product of B-splines. The predicted values and residuals from this initial fit were called the master predicted values and master residuals respectively. The computed volume was then put in a list of volumes that was added to in subsequent passes of the algorithm. The algorithm then iterated two hundred times as follows. In the first iteration the master residuals were sampled randomly uniformly with replacement. The resampled residuals were then added to the master predicted radii and a new fit was performed. The new computed volume was added to the volume list and the master residuals sampled randomly uniformly with replacement for the next iteration. The process continued for all two hundred iterations. The standard deviation of the volumes in the volume list was then used as an estimate of the volume uncertainty and the average volume was taken as the reference volume for the chosen grid size. We found that in all cases the expanded uncertainties remained approximately constant and that is why only upper bounds for the expanded uncertainties are reported in the tables below.

Calibrated Sphere Volume Estimates and Uncertainties for the Coarse and Dense Data
As a test of the accuracy of the B-spline volume estimation procedure, we applied the method to the coarse and dense data sets for the calibrated sphere. If the volumes computed by the B-spline method in Table 6 are compared with the sphere fit results in Table 3, we see that the B-spline method and the sphere fit method results differ by approximately 0.01 % for both the coarse and dense data sets. This indicates that the B-spline approach is performing as adequately as the straightforward sphere fit model and thus can be relied upon to produce good results when applied to the phantom data.

Phantom Volume Estimates and Uncertainties for the Coarse Data
In this section we discuss briefly the results of computing the volumes and volume uncertainties for the two phantoms using the coarse data sets. As described earlier, twenty one surface grid cases were used and the results are shown in Tables 7 and 8. Figures 16 and 18 show the trends of the volume estimates as the grid sizes increase. The estimates rise rapidly, and both tend to what appears to be stabilized values. The figures graphically display the values in the tables.
The spherical fit to the Green phantom data in Table 4 indicates an estimated volume of 4331.3 mm 3 . Table 7 indicates an estimated volume of 4331.84 mm 3 for the largest grid size of 210 θ by 420 φ. This is approximately a 0.01 % difference suggesting that the Green phantom is very nearly spherical.
The spherical fit to the Pink phantom data in Table 5 indicates an estimated volume of 3862.6 mm 3 . Table 8 indicates an estimated volume of 3854.68 mm 3 for the largest grid size. The difference is approximately 0.21 %, suggesting the slight non-spherical shape of the Pink phantom.
From Tables 7 and 8, the uncertainties are shown bounded on the coarse surface grid by 12 mm 3 for the Green phantom and 28 mm 3 for the Pink phantom. This larger uncertainty for the Pink phantom might be due to the non-spherical nature of the Pink phantom. As indicated here, the uncertainties for the Pink phantom are approximately twice those of the Green phantom.

Volume Estimates and Uncertainties for the Dense Data
In this section we report the results of volume estimates and the extended uncertainties for the dense data sets for both phantom surfaces. The first thing of interest is that the stabilized volume estimates for the Green dense data in Table 7 differ from the volume estimate for the Green coarse data by about 0.82 %. This is in the negative direction in that the volume estimates for the Green dense data were smaller than for the Green coarse data. In the case of the volume estimates for the Pink dense data in Table 8, the differences with the volume estimate for the Pink coarse data is approximately 1.33 %. This is in the positive direction in the sense that the volume estimate for the dense data is slightly larger than the volume estimate for the coarse data. These results are consistent with the results of the spherical model fit to the Green and Pink phantoms in Tables 4 and 5. Another thing to notice is that the uncertainties are lower for the dense data sets as shown in Tables 7 through 8 compared to those estimated from the coarse data, although the uncertainty for the Pink phantom is about twice that for the Green phantom. These uncertainties are consistent with the uncertainties computed during the fits of the sphere models in Tables 4 and 5. It is not clear why the results for the dense data differ slightly in the directions they do from the results for the coarse data. A complete analysis of these differences is beyond the scope of this paper, but would be appropriate for a detailed study related to an analysis of the fitting algorithms and their sensitivity to the surface data distribution.

Summary
The B-spline surface model joined with the Divergence Theorem seems to be a viable approach for estimating volumes of the near-spherical molded phantoms, but the volume results seem to depend strongly on the distribution of the data points on the surface. In the case of sparse data the problem of extending a surface fitted to the data to a grid on the surface leads potentially to unwanted oscillations in the neighborhood of some of the sparse data points. In the current algorithm this problem was dealt with by adding regularization terms to the objective function. These provided a balance between the fit and the surface smoothness in order to control overfitting of the data. The final regularization parameter was λ = 0.32 for all data sets.
The results obtained from both the coarse and dense data distributions appear to be consistent with the results from the spherical model fits. As shown in Table 7 the volume estimates for the Green phantom are very close to the estimates obtained by the spherical model in Table 4 suggesting that the Green sphere is essentially spherical in shape. The volume estimates in Table 5 for the Pink sphere are higher by about 0.19 % for the coarse data and 0.15 % for the dense data than the B-spline model estimates in Table 8 suggesting that the Pink phantom might be close to spherical, although it does show some non-spherical tendencies. The uncertainties in the Pink phantom case suggest though that the Pink phantom may indeed have a slightly non-spherical shape. This is also suggested by the difficulty with holding the Pink phantom in the vacuum chuck during one of the surface data probing experiments.
In summary we conclude that the combination of using the CMM to obtain surface data and the B-spline/Divergence Theorem volume estimation method can produce useful results but that there are some limitations. First, the CMM is primarily used for probing manufactured metallic artifacts and its use in probing non-metallic artifacts such as the FDA phantoms can likely lead to some larger uncertainties in volume estimation. Second, the distribution of probed points can lead to results that indicate the models used are sensitive to the point locations. Third, the methods employed may not provide useful volume estimation for more complex artifacts that surely would arise in developing simulated lung cancer phantoms. The current artifacts were near-spherical and even these did not provide fully expected results based on data distribution for the Pink phantom. Further research on the affect of surface data distribution is required. Fourth, the Pole Problem is definitely a limitation that required regularizing of the objective function. Finally, the number of data points obtained on the surfaces was limited by the relative size of the artifact and the CMM probe size. code by Prof. Timothy Sauer, George Mason University. Finally, we gratefully acknowledge Dr. Nicholas Petrick of FDA for the loan of the Green and Pink phantoms used in this study.

Disclaimer
Certain commercial software products are identified in this paper in order to adequately specify the computational procedures. Such identification does not imply recommendation or endoresement by the National Institute of Standards and Technology nor does it imply that the software products identified are necessarily the best available for the purpose.