In the following, we shall briefly review methods for estimating curvature in volumes, and then see howκmis computed for distance field volumes.
The common assumption is that the voxels correspond to samples of a scalar field f :R3→R, and that the surface whose curvature we are trying to estimate is an iso–surface{x|f(x) =τ}. Furthermore,f should be at least twice differentiable and the gradient∇f must exist at all points on theτiso–surface. In the context of the Level–Set Method, we are interested in the case wheref(·) = Φ(·, t).
7.5 Estimating Mean Curvature 149
The older of the two previous methods [141, 142, 115], consists of fitting para-metric patches to a neighbourhood around a given surface point. The position and normal of the surface point is used to initialize a local orthonormal coor-dinate system where one axis is parallel to the normal and the two other axes are perpendicular. A simple parametric patch z=h(x, y) is then fitted (using local coordinates) to the cloud of neighbouring points using e.g a Kalman filter [115]. The fitted function is the simple second order function
h(x, y) =ax2+bxy+cy2 (7.24) and once its coefficients a, b, and c are determined, the principal, mean and Gaußian curvatures can readily be computed using formulae from [32].
The other approach (and the one that is usually employed in the context of the LSM) is to estimate the first and second order partial derivatives of the scalar field (e.g. fx, fxy . . . ). From these derivatives the various curvatures may be computed directly [123, 116, 180]. (7.25) may be found in [37] or [156]5. A more detailed description of a similar method is given by Monga et al in [116]. In this paper, the authors wish to find the principal curvatures at previously detected surface points in acquired volumetric data sets. It is shown that the normal curvature of the iso–surface at a given pointpin the direction of a tangent vectorwis
κw= wTHw
||∇f|| (7.26)
where H is the Hessian matrix (i.e. the matrix of second order derivatives) of f evaluated atp. Monga et al. estimate the Hessian matrix using the so–called DeRiche filters to estimate the partial derivatives. Subsequently, (7.26) is used in conjunction with a Lagrangian Multipliers technique to find the principal cur-vatures, i.e. the normal curvatures in the principal directions. Unsurprisingly, this method turns out to be more efficient than the scheme based on fitting [116]. In the case of distance field volumes, the underlying scalar field is really a distance field which simplifies matters.
5In Sethian’s book [156] the factor 12 is missing. Another erroneous version of the formula is found in [180, 178]
Erich Hartmann has observed that for 3D distance fields which he calls3D Hesse Normalforms, the normal curvature in the direction of a tangent vectorw is
κw=wTHw (7.27)
where H is now the Hessian of a distance fieldV. (7.27) follows directly from (7.26), because the gradient is unit length in the case of distance fields. Further-more, two of the eigenvalues of the Hessian are the principal curvaturesκminand κmax, and the third eigenvalue is zero6 [77].
By definition, the mean curvature is κm = κmin+κ2 max. To compute the mean curvature, we could approximate the Hessian and find the eigenvalues, but since the last eigenvalue is zero,κmin+κmax is equal to the trace of H (the sum of eigenvalues being equal to the trace of a matrix). Hence, we can use the simpler formula [77]:
κm= Σ3i=1Hi,i
2 (7.28)
In order to apply this result to DFVs, we need a discrete operator for approxi-mating the Hessian. Fortunately, gradients are stored in the volume, and since the gradient is really a vector of first order partial derivatives we can estimate the Hessian simply by applying a central differences filter to the gradients:
H1,iestH2,iestH3,iestT
This method is a bit unorthodox, and should be compared to previous meth-ods. A well–known and more conventional way to compute second order partial derivatives that was used in [180] is the following scheme:
∂2f
If (7.31) is used to compute the second order derivatives along the diagonal of the Hessian, we get another estimate of the mean curvature.
κm=
Σ3i=1G[p+vi] +G[p−vi] −6G[p]
2 (7.33)
6Since a distance field changes linearly when moving in a directionnperpendicular to the surface, the gradient does not change in that direction, andHn=0
7.5 Estimating Mean Curvature 151
In other words, we have two schemes: One that employs the usual methods for computing second order derivative (i.e. (7.31)) and one that is based on the fact that the gradients already stored in the volume can be used in the computation of second order partial derivatives.
Of course, it is possible to use the standard method (7.25) with either of the two methods for computing the second order partial derivatives. This yields a total of four different ways to estimate curvature. In the following, the methods that rely on the fact that we are dealing with a sampled distance field will be called DF methods. The usual methods based on (7.25) will be called NDF (non distance field) methods. Methods that use (7.31) in the computation of the second order partial derivatives will be called CD (central differences) methods.
The methods that compute the second order partial derivatives using gradient information will be called GCD (gradient central differences) methods.
Hence (7.30) is the DF–GCD method and (7.33) is the DF–CD method. The two last methods based on (7.25) are NDF–GCD and NDF–CD.
7.5.1 Testing Curvature Filters
It is important to know how these four methods compare with respect to speed and precision. To test the latter, we need to compute the curvature for a voxel model where the exact curvature is known. The sphere and the ellipsoid were selected.
For the sphere, the experiment is as follows: For 12 radii in the range from 2.5 vu to 30vu
• Voxelize a sphere of given radius.
• Estimate curvatures at 100000 random surface locations.
• Print out the maximum and average errors.
• Select new radius (i.e. loop back).
The error is always computed as the percentage of the numerical error with respect to the true curvature, i.e
errκ= 100|κm−κm|
κm (7.34)
Ellipsoid Sphere
DF NDF DF NDF
CD 26.6 74.3 CD 12.4 57.5
GCD 28.5 40.5 GCD 15.1 27.1 Table 7.1: Timings (in seconds) of curvature filters.
The experiment with the ellipsoid is almost identical, the only difference being that instead of a sphere of a given radius, r, an ellipsoid with principal axes [r2r3r] is used. Notice that for the same choice ofr, the ellipsoid has a broad range of mean curvatures. In addition, the highest mean curvature for any given ris much higher than for the sphere. This means that we should expect greater errors for the ellipsoid than for the sphere.
Both the sphere and ellipsoid experiments were run using each of the four meth-ods for computing the curvature and in each case the mean and max errors were recorded for each choice of radius. The mean and max errors for the sphere ex-periments are plotted in the graphs in Figure 7.5 and the results for the ellipsoid experiments are shown in Figure 7.6. Timings are summarized in Table 7.1. The timings are in “wall clock” seconds measured on the Athlon platform and the times are the best out of three runs.
Clearly, the DF methods are faster than the NDF methods. The tests indicate that the DF methods are between 1.4 times and 4.6 times faster than the NDF methods. This can be attributed to the fact thatpow must be called because of the denominator in (7.25) which is raised to the power of 3/2. Because the mean curvature is used for curvature–based smoothing in interactive volume sculpting, the DF schemes are clearly preferable.
Precision is also an issue and by looking at the plots, it seems that DF–CD is the overall most precise scheme. However, interactive tests showed that the DF–CD scheme is, unfortunately, less stable than the DF–GCD scheme. In practice, smoothing with a DF–CD smoothing tool added some noise to the surface. Moreover, Figure 7.5 indicates that the DF–CD error actually increases slightly for very small curvatures. The superior stability of the DF–GCD scheme can probably be explained by the fact that the gradients themselves have been computed using central differences. Thus, theeffectivestencil used in computing the curvature with the DF–GCD scheme is much larger than with the DF–CD scheme.
Since stability and speed are the most important issues, the DF–GCD scheme has been selected. This is also justified by the fact that although the error is vast
7.5 Estimating Mean Curvature 153
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
0 5 10 15 20 25 30
radius
Spheres, mean curvature error in percentages of true value Mean error, DF-GCD
Mean error, DF-CD Mean error, NDF-GCD Mean error, NDF-CD
0 1 2 3 4 5 6 7 8
0 5 10 15 20 25 30
radius
Spheres, max curvature error in percentages of true value
Max error, DF-GCD Max error, DF-CD Max error, NDF-GCD Max error, NDF-CD
Figure 7.5: Curvature errors for sphere test
0 5 10 15 20 25 30 35
0 5 10 15 20 25 30
size
Ellipsoids, mean curvature error in percentages of true value Mean error, DF-GCD
Mean error, DF-CD Mean error NDF-GCD Mean error, NDF-CD
0 20 40 60 80 100 120 140
0 5 10 15 20 25 30
size
Ellipsoids, max curvature error in percentages of true value Max error, DF-GCD
Max error, DF-CD Max error, NDF-GCD Max error, NDF-CD
Figure 7.6: Curvature errors for ellipsoid test