Constructing High-Order Signed Distance Maps from Computed Tomography Data with Application to Bone Morphometry

An algorithm is presented for constructing high-order signed distance fields for two phase materials imaged with computed tomography. The signed distance field is high-order in that it is free of the quantization artifact associated with the distance transform of sampled signals. The narrowband is solved using a closest point algorithm extended for implicit embeddings that are not a signed distance field. The high-order fast sweeping algorithm is used to extend the narrowband to the remainder of the domain. The order of accuracy of the narrowband and extension methods are verified on ideal implicit surfaces. The method is applied to ten excised cubes of bovine trabecular bone. Localization of the surface, estimation of phase densities, and local morphometry is validated with these subjects. Since the embedding is high-order, gradients and thus curvatures can be accurately estimated locally in the image data.


I. INTRODUCTION
This paper develops a method of constructing signed distance fields from computed tomography datasets. Primarily, the method is concerned with two-phase materials such as trabecular bone, synthetic biomaterials, and soil samples. Being biphasic defines the interface between the two materials as an orientable surface, advantageous for modelling changes over time and measuring metrics of the surface such as volume, area, curvatures, and thickness.
Classically, a signed distance field would be constructed in two steps. First, an appropriate segmentation method would generate a binary image defining the object in space. Second, the exact signed distance transform would produce the distance map from the binary image [1]. However, the exact signed distance transform computes a quantized version of the true signed distance field [2]. This produces errors for geometric flow and curvature computation algorithms because the algorithm is only initialized first-order accurate, resulting in large errors in measurements and curvature computation. Furthermore, constructing high-order signed distance maps from binary sampled signals is a non-unique problem [3], permitting multiple embeddings that satisfy the same binary image.
The task presented here is to construct the signed distance field with a high-order of accuracy from data measured using computed tomography. The field will be constructed directly from grayscale data where gradients are well-defined. Once constructed, the densities of the respective phases can be estimated. Finally, it is demonstrated that morphometry can be performed accurately on these embeddings, a task that is too noisy when computed from the exact signed distance transform.

II. METHOD A. Problem Definition
The problem is to construct a signed distance field, φ, from a density-calibrated computed tomography image, ρ. Both will be represented on a sampled, rectilinear n-dimensional grid Ω ⊂ Z n . Let the domain Ω be split between the two phases, Ω 1 and Ω 2 such that Ω = Ω 1 ∪ Ω 2 . Associate to each phase a density, ρ 1 and ρ 2 . Γ is the interface between the two materials. We represent the interface implicitly as the zero level set of an embedding function u.
If the embedding is known, the density field can be reconstructed knowing the density of the two phases.
where θ is the Heaviside function and the convention inside-isnegative is used. This construction originates in the simulation of incompressible two-phase flow [4]. In many cases, we want the embedding function to satisfy an additional constraint, namely the Eikonal equation: where ∇ is the gradient operator, | · | is the vector norm, and F is a spatially varying speed function which yields the distance function when F = 1. Here, we used the symbol φ to distinguish the signed distance field from any general embedding u.
A signed distance field is highly desirable for many reasons. First, it can be used directly for many morphometric measures such as measuring areas, volumes, curvatures, and thicknesses [5], [6], [7]. Second, the signed distance field can represent arbitrary topology surfaces and evolve them without splitting and merging rules as would be needed for parametric meshes [8]. Third, during surface evolution problems, the Eikonal equation provides a way of rapidly redistancing the embedding during evolution [4], [9]. Constructing φ from ρ is the central challenge of this paper.

B. Solver
We present a method of constructing φ from ρ in computed tomography datasets given the assumptions of Section II-A. An overview of the method is given in Figure 1.
1) Constructing an Implicit Surface: An implicit representation of the surface is constructed by shifting the image by a threshold. This causes the zero level set to be at the level of the threshold. Standard preprocessing can also be performed on the image such as smoothing to remove noise. A simple expression is used in this work: where G σ is a Gaussian kernel with standard deviation σ, * is convolution, and T is a threshold. The symbol ψ is used to denote that this is not a signed distance transform. Importantly, this is an implicit representation of the surface Γ that is subvoxel. This construction is inspired by a sub-pixel distance map method [10]. Of note, the threshold function can be made spatially varying, T = T (x), important to high resolution bone scans where cortical and spongy bone compartments are often segmented with different thresholds [11]. Similarly, alternative preprocessing filters can be used that match the specific application.
2) Solving the Narrowband: A modification of the closest point method [12] is used to construct the signed distance function in the narrowband. The concept is to take each point on the rectilinear grid, x, that is in the narrowband and find the closest point, y, on the implicit surface: where d(·, ·) is the Euclidean distance. The distance between these two points gives the embedding function: where the sign is determined from the sign of ψ(x). This is demonstrated graphically in Figure 2. The closest point method uses gradients of ψ to iteratively move a point in the narrowband onto the zero level set.
However, the closest point method originally used a signed distance field, which is not available in this case. The algorithm is modified for the embedding in Equation 4.
Consider a point in the narrowband. This point is moved in the negative direction of the gradient until it is sufficiently close to the zero level set. However, the embedding gradient points "uphill" (away from the medial axis) in both the positive and negative sections of the curve. As such, the direction is modified by multiplying the sign of the embedding so both interior and exterior points move towards the zero level set. The update equation is then: where λ is a step size, sgn(·) is the sign function, and n(y) is the normal at the point y. As is standard in image processing and level set methods, the normal is a function of the gradient: The process is repeated till the embedding at the point y is sufficiently close to zero. Equation 7 can be rewritten as a continuous partial differential equation: Visual representation of the problem. To construct the signed distance field φ, the distance to the implicit curve Γ is found. This is done first in the narrowband (gray) and then marched outwards. The surface Γ is shown as the gray level set. Vector length approaches zero as they get closer to the curve due to the regularized sign function.
where τ is a virtual time parameter. The closest point method is then a forward Euler approximation for the initial value integration problem where the point is integrated through a stationary vector field: This vector field is plotted in Figure 3. To solve Equation 9, an interpolator is needed such that the field ψ can be assessed between sample points. We construct the interpolator using basis splines: where β n is the n th order normalized basis spline, c(q) are the coefficients of the interpolant, andψ(x) is the interpolating function for ψ. Conveniently, once the coefficients have been constructed, derivatives can be evaluated quickly: This allows the computation of normals and magnitude gradients without having to compute finite differences in the image. B-splines have been well studied [13], [14], [15] and are particularly appropriate for this task. First, they have a finite support, making them rapid to evaluate. Second, the smoothness of the interpolating function is the order of the normalized basis spline, allowing a simple selection of the order of accuracy. The drawback is that the construction of the basis splines can take some time and requires additional memory to store the coefficients. However, the time cost of construction is made up by the speed of evaluation.
Combining the interpolator and vector field integrator gives the closest point algorithm presented in Algorithm 1. The integration time step is selected as λ = h, making the time step on the order of a voxel. Third order basis splines are used in the interpolant for accurate computation of normals. The integration stops when the point is sufficiently close to the implicit surface, given by h 3 . Since a signed distance function Algorithm 1 Closest Point Algorithm, CP . The vector field v is defined in Equation 10 and computed using the regularized sign equation and basis spline interpolator. Input: ψ, x, h Output: y Initialization 1: y = x 2:ψ = BuildBSpline(ψ) Process 3: while |ψ(y)| > h 3 |∇ψ(y)| do 4: y ← y − h sgn(ψ(y)) n(y) 5: end while 6: return y is not used, the stopping criterion is modified as described in Appendix A.
A numerical approximation is made for the sign function: where δ is a small constant preventing the denominator from being zero, selected to be the sample spacing, h. The regularization constant δ is weighted by the gradient magnitude for the same reasons presented in Appendix A. Such a regularization has been used previously [9]. Regularization of the sign function has the additional benefit of causing the rate of Equation 7 to approach zero as the point y n approaches the implicit curve, preventing oscillations around the zero crossing. For images with anisotropic spacing, h is selected as the minimum of the voxel edge lengths. An issue with Algorithm 1 is that it does not enforce collinearity with the surface [16] where the vector from x to y must be collinear with the surface normal at y. Even with high-order Runge-Kutta methods, integration through the vector field does not satisfy the collinearity property. In two dimensions, the candidate point y can be corrected by moving along the tangent [16]. However, in high dimensions, there is not a single tangent but a tangent plane. We extend the collinearity method to arbitrary dimensions here.
Let x be the candidate point and define p to be the current solution from CP . Define the vector q as pointing from p to x. Define T p M to be the tangent plane at the point p on the manifold. The candidate point is moved by a vector being the projection of q onto T p M .
This is shown graphically in Figure 4. The new candidate point is now the original plus this projection.
The projected vector z will be zero when it is collinear with the normal of the tangent plane, giving a natural stopping condition. At areas of large curvature, Equation 16 can cause oscillations in the candidate point. To avoid oscillations, the tangent plane vector is weighted: where β ∈ [0, 1] is the step size, taken to be 0.5 in this work. The magnitude of z can be understood as the size of the open ball at p n and enforcing the radius to be less than a power of the spacing enforces an order of accuracy. This produces slower convergence in flat areas where β = 1 would be exact, but greatly improves the convergence in high curvature areas.
The surfaces of interest have very few flat locations, making this trade-off advantageous. The order of accuracy is selected the same of CP , h 3 . This is termed the Collinear Closest Point Algorithm [16], CP ⊥ , and is described in Algorithm 2.
In rare cases, CP or CP ⊥ will not converge. This can occur when the initial point is at the medial axis of the image, the gradient of the edge is too flat causing noise in the normal estimation, or the noise in the density image is too large. A maximum iteration is set (100 in this work for both algorithms) and if the maximum iterations are reached, the point is solved with the method of Section II-B3. This causes a degradation of accuracy locally, but keeps the implementation stable.
z ← proj TpM ( yx) 5: y ← CP (ψ, y + 0.5z, h) 6: end while 7: return y 3) Extending off the Narrowband: After solving the narrowband, the signed distance field must be extended to the rest of the domain. This is done using the fast sweeping method of Zhao [17], [18]. The fast marching method [19], [20], [21] and fast iterative method [22] are also suitable algorithms for extending the narrowband, but the authors prefer the simplicity of the fast sweeping method.
As signed distance fields are non-smooth surfaces (technically, Lipschitz continuous), they cannot be better in general than first-order accurate at discontinuities [18], such as the medial axis where characteristics flow into each other. However, away from the medial axis, arbitrary order of accuracy can be achieved.
The sweeping method proceeds in three steps. First, the absolute value of the narrowband voxels are taken to remove the sign condition. This solves for the unsigned distance field from the surface. Next, an upwind quadratic solver is used to solve the quadratic equations at each voxel. Importantly, as the method is upwind, the solution is propagated from the narrowband to the rest of the domain. The solver is fast because alternating sweeping directions cascades the solution from the surface to the domain, Finally, the sign condition is reintroduced by multiplying the solved distance field by the sign of ψ. There are techniques to include the sign condition in the solver [3], but they are considerably slower due to the branching requirements and are considerably more difficult to implement.

C. Component Estimation
Now that a non-parametric representation of the interface is available, we would like to measure the parameters ρ 1 and ρ 2 . We follow here the method of Chan and Vese [6]: where ρ 1 is the density of the phase inside the embedding and ρ 2 is the density of the phase outside the embedding. In the context of trabecular bone, ρ 1 is the tissue mineral density and ρ 2 the marrow mineral density. A nice feature of this method of estimating the phase densities is that the reconstructed density image (Equation 2) will have an average density near the original image, allowing for the computation of average density from the embedded image and the phase densities. Of note, components are estimated from the smoothed density image, G σ * ρ, not the implicit surface, ψ, or the original density image, ρ. This keeps the phase density estimation and surface localization consistent with any filtering. Numerical approximations are needed for the integration and Heaviside function. Numerical integration is performed with Simpson's rule, particularly important for surfaces with few samples across their thickness. A regularized version of the Heaviside function is used: where 2 is the finite support of the numerical Heaviside θ . A finite support is advantageous for structures with few samples across their thickness as it avoids extending into the medial axis.

D. Narrowband Size
The voxels in the narrowband must be determined. Since the initial embedding is not a signed distance transform, the standard absolute value threshold cannot be used: where |·| is the absolute value function and γ is the narrowband size. Instead, the set of voxels in the narrowband is determined by thresholding the embedding below zero and performing morphological operators: where ⊕ is morphological dilation, is morphological erosion, and S is a suitable structuring element. This is considerably slower than the thresholding approach but selects a technically correct narrowband. The cross structuring element is used in this work and the radius is set to half the stencil size minus one.
This makes the narrowband large enough that a finite difference stencil placed outside the narrowband does not extend across the narrowband. The image J then has the value one for all voxels in the narrowband and zero for all elements outside the narrowband.

A. Correctness on Ideal Surfaces
A series of experiments are performed on ideal surfaces to validate the correctness of the algorithms. The signed distance field of a sphere and torus are known analytically: where r is the radius of the sphere, a is the radius of the torus, and c is the distance to the center of the torus. This allows us to construct a ground truth φ, create a density image ρ, and assess the ability of the algorithm to recover φ from ρ. Additionally, two intersecting spheres are constructed, producing a continuous but non-smooth implicit surface. The surface parameters are r = 5 for the sphere, a = 2, c = 3 for the torus, and r = 5 for the double spheres with separation distance 3 √ 3. Images were instantiated on a 20 × 20 × 20 grid (nominal physical units).
We are particularly concerned with the order of accuracy of the method. Decreasing the grid size will always increase the accuracy of the method, but with appropriate methods the accuracy improves to some power of the sampled spacing. This is important for measuring morphometric quantities such as curvature that require a certain amount of smoothness in the embedding.
In general, the order of accuracy, m, is estimated by: where φ is the ground truth, φ h is the solution at spacing h, φ h/2 is the solution at spacing h/2, and · p is an p norm.
Since the distance transform is only Lipschitz continuous at the medial axis, different norms are used to approximate measures of the global (everywhere) and local (away from medial axis) order of accuracy. Because the signed distance map is Lipschitz continuous, the method can only be firstorder accurate when measured with the ∞ norm but can achieve arbitrary order of accuracy when measured with the 1 norm [18]. 1) Order of Accuracy in the Narrowband: The order of accuracy for solving the narrowband is measured. The analytic signed distance transform of the sphere, torus, and doublesphere are instantiated on a grid with a prescribed spacing. Using Equation 2, phase densities ρ 1 = 100, ρ 2 = 0 (in arbitrary density units), and Heaviside smoothness constant = 2, the density image is instantiated. The narrowband is solved and the 1 and ∞ norm measured on the narrowband voxels. The grid size is divided by two and repeated, allowing an estimation of the order of accuracy of the narrowband initialization. No smoothing is performed such to preserve the zero crossing location and the threshold is set to the average of the phase densities: One can see this threshold is correct from Equations 2 and 4 noting that θ(0) = 0.5.
2) Order of Accuracy in Sweeping: The order of accuracy of the signed fast sweeping method is measured. The analytic signed distance transform of the sphere, torus, and doublesphere are instantiated at the narrowband with a prescribed spacing. The remainder of the domain is solved from the narrowband and the 1 and ∞ norms are measured excluding narrowband voxels. The grid size is divided by two and the process repeated, allowing an estimation of the order of accuracy of the extension method.
3) Visualizing Result: The ideal signed distance transform of a double sphere is sampled at h = 0.25 and mapped into density as in Section III-A1. From the density image, the signed distance transform is solved using the narrowband and subsequent fast sweeping method. Different level sets of the ideal and solved distance transform are visualized using the Marching Cubes algorithms [23] to demonstrate the visual smoothness and correctness of the method.

B. Trabecular Bone
10 previously imaged bovine trabecular bone cubes were used to verify the correctness of the algorithm [24]. Bone cubes were sawed to 10 mm in edge length and imaged at a nominal resolution of 20 µm. The bones were embedded with Equation 4 using a standard deviation σ = 20 µm and threshold T = 400 mg HA/cc. The signed distance transform was computed with the proposed method.
To increase the execution speed, the first-order fast sweeping method was used for extending the narrowband. A high accuracy is only needed at the narrowband for visualization, phase density estimation, and morphometry. The models here had an average of 143 million voxels, requiring roughly 20 minutes to compute the narrowband and roughly 40 minutes per iteration in the fast sweeping method on a consumer grade MacBook Pro (32 GB random access memory, 2.3 GHz processor with 8 cores). The fast sweeping method can be accelerated by only solving to a certain distance, but a highorder cannot be parallelized due to the lack of a monotonicity condition [18]. It should be noted that the narrowband method is a so-called "embarrassingly" parallelizable algorithm, where every narrowband voxel could be solved on a separate thread.
1) Embedding and Surface Localization: A visual comparison is made between the density field, the implicit surface, and the solved signed distance transform to verify the location of the level set coincides between images. The exact signed distance transform is also computed from the implicit surface to demonstrate the quantization artifact [25]. Meshes are extracted using the Marching Cubes algorithms [23] at the appropriate level set for the given field.
2) Phase Densities: Next, density estimation is validated. The phase components are estimated using Equations 18 and 19. Once estimated, the density image is reconstructed using Equation 2 and compared visually. The bone mineral density (vBMD [mg HA/cc]) between the original density image and the constructed density image are compared using regression and Bland-Altman analysis [26]. The regularization parameter was set to = 2h.
χ =K 2π (35) where V is the object volume, A is the object surface area, H is the mean curvature, K is the Gaussian curvature,H,K are the total curvatures, H , K are the surface average curvature, χ is the Euler-Poincaré characteristic, θ is the Heaviside function, and δ is the Dirac delta function. The curvatures are estimated using fourth-order central differences [7] and integration is performed using Simpson's rule from Section II-C. Furthermore, the numerical estimation of θ in Equation 20 is used. Taking the derivative of Equation 20 gives the numerical approximation to the Dirac delta function: From these definitions, traditional bone morphometrics [27] can be derived: can be computed simply as the ratio of BV to TV. The link between mean curvature and SMI was first described by Jinnai et al. [31] and extended to TBPf by Stauber and Müller [32]. The proposed method is compared to the corresponding standard methods for each measure. The original density image is blurred with the same Gaussian filter and thresheld to binarize the volume. It is standard in bone morphometry to perform a largest connected component filter to force one bone component and one marrow component in the image. In fact, connected component filtering is required to meet the assumptions on Betti numbers in connectivity density estimation [30]. As a result, the standard methods are applied to the binary image after largest connected component filtering. However, connected component filter cannot be performed in the proposed method in a straight forward manner. As a result, Euler-Poincaré characteristic and not connectivity density is used for validation. Additionally, since SMI and TBPf are both related to mean curvature, only SMI is used for comparison. BV/TV, BS, SMI, and χ are compared between the proposed and standard methods using regression and Bland-Altman plots. Standard analysis was performed with Image Processing Language (IPL v5.42, SCANCO Medical AG, Brüttisellen, Switzerland). Finally, curvature is visualized across the surface. The Marching Cubes algorithm [23] is run

A. Correctness on Ideal Surfaces
1) Order of Accuracy in the Narrowband: Measured order of accuracy is reported in Table I. For the sphere and torus, the narrowband has an order of accuracy of 4 in both the 1 and ∞ norms. This is an order higher than expected from analysis, a result seen in similar studies [21]. Since the narrowband is smooth, the ∞ and 1 norms have the same order of accuracy.
For sphere-sphere intersection, the 1 order of accuracy is 2 while the ∞ order of accuracy is one. The degenerate 1 order of accuracy comes from the implicit surface itself not being smooth, so the order of accuracy of the method cannot be achieved. Similarly, a degenerate 1 order of accuracy is seen, reflecting the presence of a non-smooth solution in the narrowband. The ∞ order of accuracy is worse than one at higher resolutions. This is caused by points in the narrowband that are non-smooth being mapped to a wrong but collinear point on the surface. At the sphere-sphere intersection, nonsmooth locations are mapped onto one of the two spheres and not the intersection, while still satisfying the collinear requirement.
2) Order of Accuracy in Sweeping: Measured order of accuracy in sweeping is reported in Table II. For the sphere and torus, an order of accuracy of 3 is seen in the 1 norm with an order of accuracy of 1 in the ∞ norm. The reduced ∞ corresponds to a shock at their center. Across all grid sizes, the number of iterations was less than 40 The non-smooth sphere-sphere intersection demonstrated an order of accuracy near 2 for the 1 norm and an order of accuracy of 1 in the ∞ norm. Errors were seen in areas of rarefaction, corresponding to the lens shape of the intersecting spheres, degrading the 1 norm order of accuracy. Outside the area of rarefaction the result was as expected. This result is well known in the sweeping method [18]. The number of iterations required for convergence was larger for the nonsmooth surface.
3) Visualizing Result: The ideal and solved double sphere surface is rendered at four level sets in Figure 6. The zero level set surfaces correspond very well with no visual indication of the quantization artifact. In the positive domain, the level set   in rarefaction. Finally, the dumbbell separates into two teardrop shapes, some errors can be seen in the roundness of the solution.
B. Trabecular Bone 1) Embedding and Surface Localization: The surface contour of ρ, ψ, and φ for a single slice are overlaid in Figure 7. Gaussian blurring reduces noise and smooths the surface while slightly modifying the level set location. However, no further modification to the level set is seen in the solved signed distance transform. The exact signed distance transform produces large errors as expected. The meshed level sets are displayed in Figure 8. Noise has been removed by blurring and the solved distance transform agrees with the filtered image. Importantly, the exact signed distance transform applied to the binarized image exhibits artifacts due to quantization of the metric.
2) Phase Densities: Measured phase densities are reported in Table III. The average bone phase density is 646.6 ± 34.0 mg HA/cc (3734.8 ± 192.0 HU) and the average marrow phase density is 3.8 ± 26.5 mg HA/cc (100.9 ± 149.8 HU), reported mean ± standard deviation. Average densities are separated on either side of the 400 mg HA/cc threshold. Negative effective densities are seen in the marrow space because the calibration model does not account for densities below water. However, the marrow space Hounsfield units are in appropriate ranges between fat (−100 HU) and general soft tissue (100 HU). The Hounsfield units of the bone phase are very large, but are reasonable given the high resolution of the imager. When weighting the bone and marrow components by their volume fraction (Section IV-B3), the image Hounsfield units range from 535.5 HU to 1575.9 HU, consistent with clinical computed tomography where partial voluming would average the components together. Consistency in the phase densities is expected since these were the same species at the same location and there were no known diseases affecting mineralization or marrow composition.
The filtered image and constructed image are shown in Figure 9. Noise has been removed due to replacing the phases with their average values. Similarly, because a regularized Heaviside function is used and the high-order signed distance transform has gradients, the constructed image has gradients at the edge between the phases. However, the constructed image has a thicker appearance than the original image with gradients not matching local hypo-intensities in the bone phase. This suggests the regularization parameter could be included in optimization. It should be noted that the zero-crossings of ψ, φ, and the constructed image correspond (Figure 7), so the object is not actually thicker. The object just appears thicker because of discrepancies in threshold, local intensities, and the regularization parameter.
A regression and Bland-Altman plot between the volumetric bone minearal density of the original (ψ) and constructed image are given in Figure 10. A small proportional bias is seen in the Bland-Altman analysis (R 2 = 66.1%) caused by selecting the same density threshold independent of the object densities. Both the bias and proportional bias is tiny (1.41%−0.69%) in comparison to the volumetric bone mineral density. The component estimation method of Equations 18 and 19 has the property that the constructed volume will always agree with the original image on average density. As a result, vBMD errors originate in the ground truth data not matching the assumptions of a biphasic material, either because of noise or that the phases do not exhibit a single density.
3) Morphometry: Regression and Bland-Altman plots are displayed in Figure 11. Agreement in bone volume fraction (BV/TV) is very high suggesting very few voxels were removed by component labelling. A slight proportional bias is seen, consistent with the proportional bias in vBMD ( Figure 9). Agreement in bone surface area (BS) is high (R 2 = 100.0%) with a small bias of 2.03 [−29.35, 33.42] (95% LoA). The proposed method underestimated surface area, consistent with a sub-voxel localization of the surface that is smaller than the binarized volume. Excellent agreement (R 2 = 100.0%) and small bias (0.02 [−0.06, 0.10]) is seen in structure model index (SMI). The standard method estimates SMI using a finite sized dilation of the surface on the order of the sample spacing. This appears to result in an overestimation of SMI relative to the infinitesimal definition of the proposed method. Finally, Euler-Poincaré characteristic (χ) shows reduced agreement (R 2 = 96.3%) with a modest bias of −1857.78 [−3241.42, −474.1265]. The standard method requires component labelling to estimated χ, removing many small components which would add to the total Euler-Poincaré characteristic. As a result, the proposed method overestimated χ relative to the standard method.
Mean and Gaussian curvature are displayed across the surface of a single sample in Figure 12. Rendered mean curvature is smooth across the surface, showing negative and positive curvature across rods and plates. The mean curvature histogram shows a spike around zero where the image border clips the bone flat. Rendered Gaussian curvature is less smooth, owing to the distribution being fat tailed across the surface. The fat tailed distribution is confirmed in the histogram. The smoothness of both the mean and Gaussian curvature agrees with the smoothness of ψ in Figure 8b.

V. DISCUSSION
A method is presented for constructing high-order signed distance transforms of two phase materials imaged with computed tomography. The narrowband is solved by integrating voxel positions through a stationary velocity field. The narrowband is extended to the remainder of the domain using a fast sweeping method. The phase densities are estimated from the original density image and computed signed distance field.  (13c) At shocks, a point in the domain has two or more closest points, consistent with Blum's grassfire burn [33]. All other points have a one-to-one and onto mapping to the implicit curve for a given level set. Rarefaction and shocks should be understood as regions in the related problem of simulating hyperbolic conservative laws.
The utility of the method is demonstrated by measuring bone morphometrics directly from the embedding. The novelty of the method is to solve the narrowband from an embedding that is not the signed distance transform. This is done by integrating narrowband voxels through a stationary velocity field until they are both sufficiently close to the zero level set and sufficiently collinear with the normal. The modification of the convergence criterion for an image without a gradient magnitude of one was sufficient to obtain the prescribed order of accuracy. The computation of a distance map from a grayscale image was inspiried by, but varies considerably from, the sub-pixel distance map method of Kimmel et al. [10] where a redistancing algorithm is used to construct the high-order signed distance map (called sub-pixel in that work). Originally, the narrowband method was inspired by closest point methods [12], [16] but extended for arbitrary dimensions and embeddings that are not already a signed distance field. The same nomenclature for the CP and CP ⊥ algorithms are used. The narrowband algorithm is most closely related to the redistancing algorithm of Chopp [21] with a different interpolator and modified convergence criterion.
The foundational assumption of the method was the construction of an embedding function that had the proper zero crossing. A smoothing filter and threshold standard to bone analysis was modified to construct the implicit surface, relying on density calibration. In general, any implicit representation will suffice as long as the zero crossing corresponds to the surface, the embedding is monotone in the narrowband, and the embedding has well-defined gradients. However, many image analysis tasks cannot make this assumption. Integration with alternative definitions of edges [34], [35], [36], [37] would increase the utility of the technique. Such methods can even be directly integrated into active contour segmentation algorithms which have a penalty term forcing the embedding to be a signed distance field [38].
The limitation of the method is that the surface must be smooth, coming from smoothness assumptions in the interpolator and velocity field. This is a good assumption for the exact task of embedding biphasic materials imaged with computed tomography where filtering has smoothed the surface. However, the accuracy degraded at continuous but non-smooth corners of the surface. As such, we do not recommend our method for redistancing [4], [9], [21].
Establishing the embedding for non-smooth surfaces is very challenging. An attempt to alleviate the issue could be to force the divergence of the vector field to zero. This is the case in cortical thickness mapping where Laplace's equation provides a one-to-one mapping between two surfaces [39]. As a result, the integration paths from one surface to another cannot cross. However, the current problem does not permit a one-to-one mapping between level sets of the embedding because there can be multiple points on the surface that are closest (in the shock regions of advection) or multiple points on the same level set that map to a single point on the surface ( Figure 13).
An improvement to the method would be to estimate the regularization parameter, , from the image data. One approach would be to estimate in conjunction with ρ 1 and ρ 2 based on a non-linear least squares. It is expected that the magnitude of the regularization parameter corresponds to the size of the point spread function of the imager. Such a result would improve the accuracy of the gradients in the constructed image. Additionally, the method could be fully automated by including a threshold estimation step. Otsu's method [40] is particularly appropriate for biphasic materials, and would give a distribution instead of point estimate for phase densities.

VI. CONCLUSION
A method is presented for computing high-order signed distance transforms of biphasic materials imaged with computed tomography. Phase densities and general morphometry can be estimated from these embeddings. The algorithm is highorder, free of the quantization artifact associated with distance transforms of sampled binary images.

APPENDIX APPROXIMATING DISTANCE FROM THE SURFACE
It is desirable to be able to select the convergence criterion based on a distance from the surface. This allows measuring the convergence of CP in physical units, permitting one to prescribe an order of accuracy.
The geometric picture is given in Figure 14. A Taylor series approximation to the implicit contour is made at the zero crossing: where δx is a directional vector. The direction is selected as the normal vector of a prescribed length: where γ is the length. Second, we know ψ(x) = 0 by definition of the zero crossing. As such, we get the following approximation for the distance from the contour: We modify the convergence criterion to be: This exact expression avoids division by zero in flat areas of ψ. Furthermore, it prevents early stopping away from the interface since |ψ| is greater than zero but |∇ψ| approaches zero near the medial axis. Additionally, the expression |∇ψ| is already available from the computation of the normal in Algorithm 1, avoiding any additional computation. Finally, the error in the approximation decreases as the point approaches the surface at a rate O(γ 2 ). If ψ is a signed distance function, Equation 46 becomes the same convergence criterion as before since the gradient magnitude is unity. In general, |∇ψ| can be seen as a normalizing factor bringing the gradient magnitude close to one such that the distance travelled in the embedding is the same as the distance travelled in the domain.