HyBoDT: Hybrid Bounded Distance Transforms of Trimmed NURBS Models

Distance fields have been used extensively for a wide range of applications such as shape manipulation [1], shape and volume representation [2], constructive solid geometry (CSG) Boolean operations [3], graphics rendering [4], path planning [5,6], collision detection [7–9], morphing and sculpting of shapes [10,11], surface offset computations [12,13], etc. Computing and representing the distance field everywhere in 3D space might not be computationally tractable. However, one of the key properties of the distance fields is that it is C⁰ continuous. Hence, a discretized distance field or distance transform (DT), where the minimum distance is computed in a regularly spaced grid, is more common in applications using distance fields. The discrete representation of distance fields is analogous to a digital representation of a continuous field or a voxelized representation of a smooth object. Due to the continuous nature of distance fields, the distance field value at any point can be approximated using interpolation from the regular grid of values with bounded error.

Significant work has been done to compute the DF and DT of polygons, triangle soups, and parametric surfaces by using a multitude of methods such as fast marching methods [14], graphic rasterization methods [15], axis-aligned bounding box tree-based search [16], scanning based methods [17,18], and jump flooding [19]. However, most methods to compute the distance fields are computationally expensive, restricting their efficient use. We discuss several related works for computing distance transforms with and without any acceleration in Sec. 2. However, to the best of our knowledge, these works do not consider bounds on the computed distance values due to the CAD model representation. Specifically, the de facto representation for CAD models, trimmed non-uniform rational B-splines (NURBS), is not directly used for computing the DTs. These are usually tessellated into triangles, which are then used for computing the DT. In this work, we quantify the error in the DT caused by tessellating the original trimmed NURBS representation of a CAD model.

Computing the distance fields of CAD models made of trimmed NURBS surfaces is challenging, and most existing methods do not compute the distance fields directly. A standard approach to deal with such CAD models is to tessellate the trimmed surfaces into triangles and perform the distance computations. The major limitation of computing exact-DFs of trimmed NURBS through tessellation is the presence of cracks [20] between triangles of adjacent trimmed NURBS surfaces. These cracks or gaps in the tessellated model render the model non-manifold or not watertight. Methods based on ray-tracing usually fail for non-watertight models. Several previous works have a focus on creating a watertight tessellation of CAD models made of trimmed NURBS surfaces by tessellating the common edges or healing the cracks between adjacent surfaces using thin triangle strips [21,22]. We deal with gaps in our approach by first voxelizing the model and then computing the DT of the voxelized model. Our initial NURBS tessellation is computed to be smaller than the voxel size, which effectively deals with small gaps. In addition, we compute the bounds based on the resolution of the NURBS surface tessellation and use it to compute the bounds of the DT.

One popular approach to compute the DT is to use a marching algorithm starting from the boundary of the object. However, there are two main disadvantages of the marching approach. First, the marching process is inherently serial, which is slow in the absence of any acceleration approaches (such as bounding volume hierarchy). Second, the marching introduces an error in computing the distances in the absence of any correction (Please see the Appendix for a detailed explanation of the marching error). In this work, we perform the marching operation in parallel on the GPU. We devise a GPU parallel algorithm that accounts for the race condition during marching. In addition, we use a “hybrid” approach to marching, where we store the distance to the closest boundary and the index of the closest boundary voxel while marching. This hybrid approach overcomes the error due to marching. Finally, we also compute the error due to tessellating the NURBS surfaces and use it to calculate the distance bounds. Thus, our approach is both fast and provides theoretical bounds for the distances, making it practical to compute the DT of real-world CAD models represented using trimmed NURBS surfaces.

In this paper, we use graphics processing units (GPUs) to compute DT for CAD models consisting of parametric surfaces (specifically, trimmed NURBS surfaces) as shown in Fig. 1. We first explain the GPU accelerated computation of distance transforms using our hybrid marching wavefront method (see Sec. 3). We then analyze the distance bounds for a trimmed NURBS surface representation of the CAD model. The specific contributions of this paper include the following:

1. Hybrid GPU-accelerated distance transforms for trimmed NURBS objects using a marching wavefront algorithm that eliminates the marching error. The error in the final DT is only due to tessellation of the surface.

2. Theoretical bounds for the computation of distances from tessellation of trimmed NURBS surfaces, which can accurately bound the distance transform.

The rest of the paper is arranged as follows. We discuss some research works directly related to our approach in Sec. 2. We then explain our hybrid marching algorithm in Sec. 3. We show the calculation of the bounds on the DT based on the NURBS representation of the object in Sec. 4. Finally, in Sec. 5, we provide some experimental validation of the bounds and computational time for DTs of different objects.

The work done by Rosenfeld et al. [23] provided the preliminary idea of computing distances or performing distance transforms in a digitized picture. Initially, distance transforms were performed in two-dimensional image space [24–26]. Recently, with the advancement in computing power and the need to represent volume data, researchers are focused on performing distance transforms in a volumetric or three-dimensional grid [27]. To accomplish this task, several algorithms have been developed such as chamfer and vector distance methods [28–31], adaptive and complete methods [2,32], level sets and fast marching methods [14]. In the chamfer distance transforms, the minimum distance from the boundary of the object is computed from the center of a voxel to the minimum distance from its neighboring voxels, and then, this value is marched according to a marching scheme (such as wavefront or sweeping schemes). In vector distance transforms, a Euclidean distance from the voxel center to the boundary of the object is computed, and the minimum of those distances is considered for the transform. This results in increased accuracy of the computed distance field. Sethian [14] developed a fast marching method to compute distance transforms that use Eikonal equation solutions at the boundary for advancing fronts. Frisken et al. [2] used the adaptive subdivision of voxels employing an octree structure to develop an adaptive distance field representation. The method developed by Xu et al. [33] uses winding numbers to compute signed distance field for triangular meshes that also includes non-manifold meshes by offsetting the non-manifold surface distance transforms to fill small gaps and holes in the meshes. Further, several other works use tree-based traversal for computing distances quickly and efficiently [16,34]. In this paper, we implement a combination of chamfer and vector distance transforms, with a hybrid wavefront marching scheme that stores the indices of the closest boundary voxel rather than the distance itself to compute the minimum distances.

Additionally, there have been considerable efforts in accelerating the distance field computations using optimized algorithms and graphics hardware (GPUs) [35,36]. Sud et al. [15] presented a fast distance field computation algorithm that uses graphics hardware features such as culling and clamping to effectively compare the distance fields. Cao et al. [37] and Bastos and Celes [38] used GPUs to compute vector distance fields and adaptively sampled distance fields, respectively. However, accelerating chamfer and vector distance transforms is challenging on the GPU, since they involve sorting the voxels according to priority queues, which is difficult to implement in data-parallel GPU architectures. Recently, Zampirolli and Filipe [17], Man et al. [39], and Manduhu and Jones [40] perform exact distance transforms using GPU acceleration for 2D images. However, the algorithms described are limited to application on 2D images and are relatively difficult to implement. In this work, we accelerate the process by computing each march of the voxels from the object boundary, as parallel operations, using the GPU. Our algorithm naturally maps to the GPU architecture, which makes implementing our method relatively straight-forward.

To the best of our knowledge, there is little on computing the bounds for distance fields. Schneider et al. [41], Rong and Tan [19], and Cignoni and Sochor[42] compute approximate distance fields using GPU algorithms. They provide theoretical bounds on the computation of approximate distance transforms as compared to actual distance fields. However, they do not take into consideration the error due to CAD model representation. In this paper, we account for both the errors by using a hybrid approach to correct the error due to marching and quantifying the error due to the CAD model representation. To the best of our knowledge, this is the first work to compute distance bounds and computation of distance fields for trimmed NURBS representation of CAD models.

We compute the distance transforms of manifolds (of dimension < = 2) embedded in 3D space or a curve embedded in 2D space. Without loss of generality, we assume that these are the boundary-representation (B-rep) of a solid model in 3D or the closed boundary of a shape in 2D. The indicator function value (0 or 1) is used to represent if any part of the object is contained in a voxel. Our marching algorithm to compute the distance transforms starts from the boundary voxels. Boundary voxels can formally be defined as follows:

The function neighbor(v) depends on the definition of neighborhood of a voxel in a voxel grid. In 3D, this could be defined as the 6 neighboring voxels that share a face with the current voxel, or the 26 neighboring voxels that share a face, edge, or vertex with the current voxel. We use all the 26 neighboring voxels for calculating the distance transforms in 3D. Similarly in 2D, we use the eight neighboring pixels. We define the distance transform function as follows:

In other words, the distance between each boundary voxel and the voxel v ∈ V is computed, and the minimum of these distances is assigned as the distance transform of the voxel v. This is done for each voxel in the grid. A naive implementation of this algorithm has a work complexity of $O(d.|VB|)$⁠, |V_B| refers to the size (cardinality) of the set of boundary voxels.

We use an hybrid marching wavefront scheme to compute the distance transform. Note that we are only computing the unsigned distance transforms here; hence, we do not use any more information relating to the inside-outside of the B-rep solid model as computed in Ref. [43] (thus giving us the ability to compute the unsigned distance fields for a voxel grid corresponding to non-manifold entities that need not represent a solid model).

The algorithm to compute the bounded distance transforms consists of two steps: (i) computing the distance values and the bounds for the boundary voxels and then (ii) propagating that distance value and bounds to other voxels. The minimum distance from the object boundary to the voxel center in case of a boundary voxel is computed based on the boundary representation. For a voxel-based representation, the bounds in distance to object boundary are $(−g/2,g/2)$⁠. The nominal distance value of the boundary voxels is considered zero, and hence, all the boundary voxels are assigned this value. However, we account for the error in computing the distance from the object boundary to the boundary voxel center for different representations. We cover the case of trimmed NURBS representation in Sec. 4. We have provided the bounds for a triangle soup representation in Appendix  B.

Once the distance for the boundary voxels is computed, we propagate the distances from the boundary voxels by one neighborhood of voxels (based on the neighborhood definition) in each pass (see Fig. 2). For example, in 2D, the marching can be performed in two directions, orthogonal and diagonal. Examples of cases where both directions of marching is possible are shown in Fig. 2. Although the distance obtained from the orthogonal marches might be temporarily assigned to a voxel, in some cases, the minimum distance might be in the diagonal direction, thereby incurring a temporary error in the distance fields. We correct this error by propagating the closest boundary voxel index instead of the distance value. At each pass, the field value computation marches by one neighborhood ring. The algorithm stops if the total sum of the distance field values do not change for two successive iterations or a maximum number of marches of $floor(2×max(nx,ny,nz))$ is performed.

One of the significant disadvantages of performing marching based approaches is the error accumulated because of performing orthogonal and diagonal marches (see Appendix  A for details). To circumvent this issue, instead of storing the minimum distance computed at every step, we store the argmin of the distances (i.e., the location of boundary voxels leading to the minimum distance). This approach allows us to correct the error accumulated due to the orthogonal and diagonal marches. In Sec. 4.2, we prove that this approach computes the correct DT.

The individual computations for each voxel can be computed independently in parallel in our marching algorithm. We map the voxel grid to the 3D block and grid structure in NVIDIA CUDA [44]. The GPU kernel computes the distances for each i, j, and k voxel. The kernel is called once for each pass with multiple passes for the marching. Several neighboring voxels might contribute to the minimum distance of a specific voxel. Still, not all the contributing neighboring voxels have the correct minimum distance value computed yet, which causes a race condition shown in Fig. 3. To overcome this race condition, we make use of a ping-pong approach where, we store two copies of distance fields at each pass, the current and previous distance fields (or DTM-A and DTM-B in Fig. 3). We use DTM-A as the reference for performing all the distance comparisons and assign the computed distance values in DTM-B (hence, avoiding the race condition). We then synchronize the threads after each marching pass and then swap the previous and current pointers to the distance fields. In each kernel call, the distance values march by one layer of voxels, and we check if a particular neighboring voxel distance value has changed by the march and recompute the minimum distance using the newly computed distances for the neighbors and update the distance value for the voxel in DTM-B. After each march, we check for the change in the sum of the field values from the previous march and stop if the change is smaller than a user-defined tolerance. The complete algorithm is shown in Algorithm 1.

Another important aspect of the implementation is the memory complexity of the proposed algorithm. To avoid the race conditions, we store two copies of the distance field matrix. If the voxel grid resolution is (n_x, n_y, n_z), then the total number of values that need to be stored in GPU memory is 2n_xn_yn_z. If the maximum dimension of the voxel grid is n = max(n_x, n_y, n_z), the memory complexity of Algorithm 1 is $O(n3)$⁠.

Since a NURBS surface is a mapping from a two-dimensional parametric space to a three-dimensional Euclidean space, any tessellation (triangular, quadrilateral, etc.) applied to the NURBS surface in the parametric space has a direct equivalent tessellated structure in the three-dimensional Euclidean space. Performing the tessellation in the parametric space is an efficient method for spline surfaces. It eliminates the complexity of the three-dimensional triangulation algorithms while ensuring that the shape of the surface is maintained during the tessellation operation [45]. There are several ways to tessellate trimmed surfaces directly in the 3-dimensional Euclidean space [46]; however, this approach would disconnect the NURBS surface from the parametrization, adding another layer of complexity. Although it is also possible to apply quadrilateral tessellation for surfaces, we use triangular tessellation for better compatibility with triangle soup representation. Figure 4 illustrates our tessellation algorithm for trimmed and non-trimmed surfaces.

Before the triangular tessellation, we evaluate the surface using a relatively small evaluation interval delta, δ = 0.01 in both parametric dimensions to generate a large number of surface points to triangulate. We use δ value to change the number of triangles generated. In addition, we use the vertex spacing, δ_v, to control the size of the triangles. Vertex spacing can be considered an integer multiplier to the δ value, which increases the jump distance between the surface points; therefore, making the edges of the triangles longer. Using two separate variables for triangulation allows us more flexibility and control over the resolution of the final shape. Nonetheless, using a smaller δ value with δ_v = 1 would allow us to generate more triangles. The tessellation algorithm considers four parametrically adjacent points at a time, which forms a quadrilateral structure in a counter-clockwise direction and creates triangles depending on the intersection of the quadrilateral structure with the trim curve (indicated with green dashed line).

When there are no trim curves attached to the surface, the algorithm continuously triangulates the surface using the user-specified δ and δ_v values. If there are any trim curves, the algorithm first determines which points are inside or outside the trim curve with the point-in-polygon algorithm implemented using a winding number approach [47]. The algorithm then uses a ray intersection test to determine which points intersect with the quadrilateral structure and the trim curve. After obtaining the outside and intersection points, a closed polygon is constructed and triangulated. The barycentric coordinate of each new triangle is also checked with the point-in-polygon algorithm to ensure that all triangles are outside the trimming region. This last step is only required if the trim curve is rational, and the weights are different from 1.

Our tessellation algorithm provides a mechanism to control and change the tessellation density and the number of triangles. We directly generate triangles using a user-specified specified δ and δ_v variables in one pass. Our trimming and tessellation algorithms are easily parallelizable since each set of four points can be processed in parallel. Once the tessellation of the trimmed NURBS surface is obtained, computing the distance transforms follows the procedure mentioned in Sec. 3.

We provide the following three results: (i) the correctness of the marching (Theorem 1), (ii) the guarantee for computing the DT in $⌊(2max(nx,ny,nz))⌋$ passes (Lemma 1), and (iii) the bounds for the distance computation from trimmed NURBS representation of the CAD geometry (Theorem 3).

We extend Theorem 2 to obtain a bound on the distance transforms for trimmed NURBS representation. For this work, we use key results from the theorems in Krishnamurthy et al. [48]. This work theoretically computes the bounds for distance computation between a point and the trimmed NURBS surface tessellation. We could formally state the theorem as follows.

In this section, we provide some anecdotal results and comparisons for the proposed algorithms and bounds in the above sections. We implemented the DT algorithm in python3.6 [49]. The backend for the GPU-accelerated code is written in c++ using the Pybind11 API [50] and CUDA toolkit [44] for GPU acceleration. We also use numba [51] for performing a just-in-time compilation of the CPU version of the same code to perform a fair comparison between the CPU and GPU performance. The Trimesh library² is used for the read/write operation of triangle files and subdivision of triangular meshes for performing the surface voxelization. For NURBS surfaces, we implemented the tessellation algorithm using c++ on the CPU. We make use of the OpenCASCADE library to read NURBS models from STEP files. We used the GPU-accelerated voxelization code developed by Young and Krishnamurthy [52] to compute the initial boundary voxelization and bounds computation.

To the best knowledge of the authors, there are no publicly available source codes with CPU or GPU implementation for distance fields computation of trimmed NURBS surfaces. Therefore, before we demonstrate the performance of our proposed approach on trimmed NURBS surfaces, we first show our results for voxel-based representation (in 2D pixels and 3D voxels) and triangle soup representation and perform comparisons with publicly available source codes. Then, we show the performance of our approach for trimmed NURBS surfaces.

Examples of the 2D distance transforms are shown in Fig. 6. They also show the maximum and minimum bounds of the distance transforms for the two binary images. The error bound shown in the figure is the difference between the maximum bound DT and the minimum bound DT. It can be seen the maximum values of the bounds in the DTs occur in the neighborhood of the boundary voxels.

The 3D distance transforms (unsigned and signed) for the tessellated model of an Armadillo and Bunny are shown on the left of Figs. 7 and 8, respectively. The unsigned DTs are cross-sectional views, and the signed distance fields are volume renderings. The actual distance transform is the minimum distance to the triangle soup computed at each voxel center using a brute force approach. This is computed by using a KDTree-based representation of the triangles implemented in the Trimesh³ library to calculate the distances efficiently on the CPU. The distance of each voxel center to the triangles present inside the voxel is computed using the CPU, while the marching to compute the DTs is performed on the GPU. The rightmost images in Fig. 8 shows the minimum and maximum bounds of the DTs. It can be seen that the minimum bound DT has a considerable gap (lower values) at the boundary of the object, while this gap is not present in the maximum bound DT. This is in accordance with the min-max bounds computed.

The DTs for trimmed NURBS representations of the Ducky, Hammer, and Scooby are shown in Figs. 9 and 7, respectively. The NURBS models are tessellated using a δ of 0.05 and a voxel grid resolution of 128³. The minimum and maximum bounds of the computed distance transforms are shown in Fig. 10 for the ducky and hammer models and in Fig. 7 for the Scooby model.

The CPU and GPU timings for the computation of the DTs for voxels, tessellated 3D models, and NURBS representations are shown in Table 1. The timings are an average of three runs for all the models. We also provide a comparison with an implementation of distance transform computation using Scipy [53].⁴ We observe that our GPU implementation performs better than CPU implementation and the Scipy implementation at higher resolutions. At lower resolutions, the speedup is smaller but still substantially faster. The smaller speedup at lower resolutions is due to the overhead involved in transferring the CPU data to GPU and performing the computations on the GPU. The GPU implementation performs well for even very high resolutions. The computational time scales almost linearly with the resolution as the resolution increases. We also report the GPU memory usage in Table 1 ⁵ the GPU memory for even very high-resolution voxel grids is approximately 250 MB, allowing our algorithm to run efficiently on current-generation consumer GPUs.

We compare our distance transforms timings with actual computation of the distance fields for triangle soups in Tables 2 and 3. All the computations were performed using an Intel Xeon 2.40 GHz CPU with 384GB RAM and an NVIDIA TITAN RTX GPU with 24GB GPU RAM. We perform actual distance field computations using a CPU implementation provided⁶ in Tables 2 and 3.⁷ For a fair comparison, we only compare the timings of the CPU implementation of DT. While the comparison here is performed only on CPU, we observe a good speed-up for more complex models like Bunny and Armadillo. However, with simpler models such as the Cessna, the time taken by libigl [43] is lower than our DT. Note that for these comparisons, we compute the DT on the CPU. However, our marching algorithm is essentially devised, keeping the GPU acceleration in mind. Whereas the other methods do not use a marching-based algorithm for performing the distance computations.

In this paper, we have developed an hybrid marching method to compute the bounded distance transform for trimmed NURBS models. Our algorithm eliminates the error due to marching by storing the propagating the closest boundary voxel instead of the distance values. We compute the theoretical bounds for the distance transforms obtained using our algorithm. We empirically show that the computed distance transforms are accurate in replicating the accurate distance fields within our computed bounds. We have a GPU implementation and provide running time comparisons between the CPU and GPU implementations.

Bounded distance transforms provide a better object representation, especially for applications such as object recognition and collision detection. Future work will involve improving our GPU implementation and experimenting with volumetric data structures such as Octrees to reduce memory requirements of storing the distance transforms.

This work was partly supported by the National Science Foundation under Grant Nos. CMMI:1644441 and OAC-1750865.

There are no conflicts of interest.

The authors attest that all data for this study are included in the paper. Data provided by a third party listed in Acknowledgment.

Directly propagating the minimum distance value in a marching algorithm leads to an accumulation of the error.

In this section, we details on the distance computations for objects represented using triangle soups. Since we only compute the unsigned distance fields, there are no specific requirements for the object representation to be two-manifold. If the geometry is manifold, we can compute the signed distance field by postprocessing the distance transforms using the inside-outside information from methods such as Barill et al. [43] and Jacobson et al. [54]. For computing the unsigned DTs, we first need to identify the boundary voxels (i.e., the voxels that intersect with the triangle soup).

We use a subdivision-based method for identifying the boundary voxels that intersect with the triangle soup. We subdivide the triangle soup until the maximum edge length is less than the resolution of the voxel grid. Thus, identifying the boundary voxel grid is simplified to a point membership classification of a point (vertices of the triangle) in an axis-aligned bounding-box (AABB). Then, we compute the distance for all the points in the voxel and check for the maximum and minimum distances. Similar to Sec. 4.2, we compute the maximum and minimum bounds of the distance fields.

The algorithm remains the same as the Sec. 4.2. The only difference is computing the bounds of the boundary voxels using triangle soups. Once the precomputed distance field D₀ is obtained, we perform the marching algorithm, which computes DTs for all the other voxels. Using the unsigned distance field, we compute the signed distance field by assigning the sign to the unsigned distance value based on the inside-outside test on manifold geometries.

Now, we compute bounds for DTs using the triangle soup representation. Theorem 2 can be extended to this case with ε_t = [−d_max, d_max] instead of $ϵv∈[−(g3/2),g3/2]$⁠, where ε_t represents the bound in computing the minimum distance from the triangle soup and d_max represents the distance from the voxel center to the object represented using triangle soup. d_max is maximum $g3/2$⁠, when the nearest vertex of the triangle is at the corner of the voxel. However, in general, the d_max would be the distance to the closest vertices of the triangle (distance to closest vertices of the triangle ensures that the gaps in closing the object do not affect the minimum distance and hence bounds the distance loosely). With this, the bound in computing the distance field from a triangle to the boundary voxel is ε_t ∈ [−d_max, d_max].