Theoretical and Empirical Analysis of a Fast Algorithm for Extracting Polygons from Signed Distance Bounds

Recently there has been renewed interest in signed distance bound representations due to their unique properties for 3D shape modelling. This is especially the case for deep learning-based bounds. However, it is beneficial to work with polygons in most computer-graphics applications. Thus, in this paper we introduce and investigate an asymptotically fast method for transforming signed distance bounds into polygon meshes. This is achieved by combining the principles of sphere tracing (or ray marching) with traditional polygonization techniques, such as Marching Cubes. We provide theoretical and experimental evidence that this approach is of the $O(N^2\log N)$ computational complexity for a polygonization grid with $N^3$ cells. The algorithm is tested on both a set of primitive shapes as well as signed distance bounds generated from point clouds by machine learning (and represented as neural networks). Given its speed, implementation simplicity and portability, we argue that it could prove useful during the modelling stage as well as in shape compression for storage. The code is available here: https://github.com/nenadmarkus/gridhopping


Introduction
A Signed Distance Field (SDF) for some shape S is a function f S : R 3 → R such that f S (x, y, z) returns a geometric distance from the point (x, y, z) ∈ R 3 to S. If (x, y, z) is within S, then the returned distance is negative.In practice we often cannot obtain the exact distance to the shape, but work with distance bounds that underestimate the distance some of the time and never overestimate it.This enables us to efficiently represent a larger class of shapes.J. C. Hart [9] (see also [10]) uses the definition that f S is a signed distance bound (SDB) of S if and only if for all x ∈ R 3 we have where f −1 S (0) = {z : f S (z) = 0}.The earlier mentioned sign rule applies to the interior of S. It should be noted that the term SDF is often used in the literature even though SDB is more appropriate [9,10].
Signed distance bound representations show significant potential due to their rendering speed, storage efficiency, implementation simplicity and shape distribution modelling capabilities.They are powerful tools which are used in computer graphics, simulation, and medical imaging domains.Here are some typical applications: • shape compression -it may be more storageefficient to store a signed distance bound (e.g., represented by a small neural network [5]) than as a triangle mesh; • statistical shape representations for generative modelling -sampling from such a representation [3,18] can potentially help add variety to virtual worlds; • shape recovery from point clouds and noisy measurements [12,6]; • shape manipulation -such as in the paper by Hao et al. [8]; • computer-aided design (CAD).
In most graphics applications (e.g., classical 3D game engines), due to GPU hardware specifics and well-developed prior practices, it is beneficial to work with polygons even though some assets are encoded as SDBs.Therefore, converting a signed distance bound into a polygon mesh warrants investigation 1 .
The objective of this paper is to present and evaluate a method for transforming signed distance bounds into polygon meshes.The contributions are: 1. we describe a fast algorithm for polygonizing SDBs based on sphere tracing (ray marching), 2. we theoretically show that this algorithm is asymptotically faster than the obvious approaches (e.g., Marching cubes [14] over the whole grid), 3. we empirically confirm this in practice for signed distance bounds encoded with deep neural networks.
The next section provides context, basic definitions and related prior work.The description of our method and its analysis/comparisons follows after that.

Basic definitions and related work
In our case, the polygonization volume is an axisaligned cube centered at the origin.The shape represented by its SDB is withing this cube.This cube is partitioned into a rectangular grid with N 3 cells by subividing each of its sides into N intervals of equal size.If we assume we are not dealing with fractals, only O(N 2 ) cells asymptotically contain the surface of our shape.Thus, the only triangles that need to be computed are passing through these cells.This puts the lower bound on the complexity of the triangulation algorithm.However, the challenge is to isolate just these O(N 2 ) cells.
The simplest solution, which leads to O(N 3 ) complexity, is to check each of the N 3 cells.This process is called enumeration [2] (Section 4.2.2 in the book) and for some purposes it is too slow.The original applications of the Marching cubes algorithm [14] applied this slow approach.However, these appli-cations usually dealt with real-world data obtained through measurments (e.g., CT scans), not well defined mathematical objects.
An improvement to this basic enumeration scheme is known as continuation [2] (Section 4.2.3 in the book).This idea of this approach is quite simple.
Starting from a single "seed" cell that intersects the surface of the shape, new cells are propagated across the surface until the entire surface is enclosed.The complexity of this algorithm is O(N 2 ) because only the surface cells are analyzed.Thus, we also call it surface tracing/crawling.Its main disadvantage is that it produces a single mesh component for each seed cell.I.e., we need a separate seed cell to polygonize all disjoint shape components.This can be problematic in practice.Another issue with this approach is caused by limited numerical precision of digital computers when representing seed cell parameters (e.g., their center point coordinates) -this complicates the positioning of the seed cell.Of course, this problems only show up for high grid resolutions, i.e., large values of N .
We also mention the unpublished work of Christopher Olah [17,16].This method aims at similar goals (transform an SDB into a triangle mesh), but uses a different subdivision strategy over the grid (recursive partitioning instead of sphere tracing) and no detailed complexity analysis was performed.It is not a formal academic work, but an implementation.
Besides the mentioned work that warrants a direct experimental comparison (enumeration and continuation [2], Marching cubes [14]) to our method, we also mention more broad work for additional context.There are a number of algorithms for speeding up isosurface extraction (as polygon meshes) from real-world data such as MRI and CT scans.E.g., based on octrees [24] or other more exotic data structures [4,13].However, these methods assume that data points on the grid are already known (sampled).Thus, direct application of these algorithms to polygonizing SDBs is not possible.
There is also a large body of work of generative models for geometry.In this context, we mention those that directly output polygons.Gao et al. [7] represent geometry as a signed distance field defined on a deformable tetrahedral grid and incorporate this into a learning framework.Their end result is a neural representation that can directly generate polygonal models from the training data distribution.In the PolyDiff paper [1], the authors introduce a mesh generator based on a denoising diffusion process [20,11].In the MeshGPT paper [19], the authors use a decoder-only transformer neural network to directly output the polygons of a mesh.All these works are learning-based, i.e., require a large dataset in order to be effective, are not easy to implement/use in a production system like a game engine and require significant computational resources.
The next section provides a detailed (mathematical, algorithmic) description of the gridhopping method to enable its analysis.

Method
Another way to speed up the triangulation process is to apply a variant of ray marching called Sphere Tracing (described by John C. Hart [10,9]).Note that sphere tracing is originally a 3D rendering scheme.However, we re-purpose it here for polygonization.The basic idea is to define a ray and move along its direction until you find the intersection with the shape or exit the rendering volume.In sphere tracing, the marching step is set to be equal to the (estimated) distance of the current point to the shape.This approach greatly speeds up the process of finding the intersection.However, unlike in the ray marching-based rendering of images, we do not stop the marching process at the surface.The marching along the ray in continued (starting at the next cell along direction of the ray) until the end of the polygonization volume is reached.Let us call this method gridhopping.This process enables us to efficiently isolate the O(N 2 ) cells that contain the surface of the shape: we present theoretical and experimental evidence that this method extracts a triangle mesh from a signed distance bound in O(N 2 log N ) steps for a grid with N 3 cells.Note that this is a significant speed improvement over the basic approach that analyses each of the N 3 cells, leading to O(N 3 ) complexity.
Without loss of generality, we assume that our polygonization volume is a unit cube centered at the origin.The grid resolution is specified by N : there are N 3 cubic cells in the grid, each with a volume equal to 1 N 3 .Each cell is assigned a triplet of integers (i, j, k) with i, j, k ∈ {0, 1, 2, . . ., N − 1}.The centroids of the cells are computed according to the following rules: A total of N 2 rays are cast in the +z direction from the plane z = −0.5 + 1 2N .Such rays have the following vector parameterization for λ ≥ 0: with ) T is the origin of ray R ij and d = (0, 0, 1) T is its direction.The (x i , y j ) pairs (N 2 of them) are computed according to equations (1).
We move along each ray using the ray marching (sphere tracing) [10,9] method.If the polygonization volume contains a shape S described by its signed distance bound f S , the following iteration describes this process: The iteration starts at r 0 = o ij and continues until |f S (r n )| is sufficiently small (indicating we are very close to the surface of S, by definition of f S ).In our case, we are only interested to move close enough to the surface to determine the (i, j, k) triplet determining the cell.Simple algebra shows that a cell possibly intersects the surface of S and we have to call a polygonization routine if the distance |f S | is less than or equal to The pseudocode in 3 contains the details, and Figure 1 provides a visual summary/illustration.
If the ray intersects the surface and we denote the closest intersection to r 0 with r * , then the above iteration converges to r * .This is because 2. on the ray between r 0 and r * , f S (r) = 0 only for r = r * ; 3. the iteration will never "overshoot" r * because f S is a signed distance bound.

z xy
Figure 1: A 2D illustration of the gridhopping method: grid cells are dashed gray squares, the shape with a SDF is in red, the z axis and the xy plane are in black.The dashed green line (parallel to the z axis) represents one of the rays along which sphere tracing is applied.Since the SDF of the red shape is known, we can compute the maximum step size along the ray: the green dots are locations obtained from these calculations.The Marching cubes polygonization is first applied to the light blue cell since it is the first one along the ray that satisfies the condition related to equation 4.

Theoretical analysis of computational complexity
We analyze the asymptotic number of steps required by the method from previous section to polygonize a shape defined through its signed distance bound.
For non-fractal shapes, there are at most O(N 2 ) cells that contain polygons.The challenge is to isolate these cells in a fast manner.The trivial way is to check all N 3 cells.This may be too slow for some applications when high resolution (large N ) is required.Our claim is that the algorithm from the previous section is faster than that: its complexity is O(N 2 log N ).
We provide evidence for this in the following steps: 1. provide a proof for polygonizing planes; 2. provide a proof for polygonizing axis-aligned boxes; 3. argue that any non-fractal shape can be approximated as a union of boxes.
These steps are explained in the following three subsections.

Polygonizing planes
A plane is a flat, two-dimensional surface that extends infinitely far.Of course, we are interested in polygonizing only the part that intersects with the polygonization volume.
The exact signed distance from a point r to a plane P is given by the following equation [23]: where r P is some point lying on P and n P is P 's normal vector such that n T P • n P = ||n P || 2 2 = 1.We analyze three different cases: two cases of axisaligned planes and one case for a plane in general position.The first case is when the plane and the rays are perpendicular.In our case, since the rays are cast in the +z direction, this corresponds to the plane z = C for some constant C. The second case is when the plane and the rays are parallel (plane specified by x = C or y = C).The third case is the plane in a general position.
Case #1: plane P and rays are perpendicular.First, notice that the orientation of the plane does not matter since the ray marching always uses the absolute value of the computed distance bound.Thus, we have two sub-cases: approaching the plane and escaping the vicinity of its surface.The approaching phase is performed in a single step for each ray since D P provides the exact distance estimate.Hence, its complexity is O(1), independent of N , the position of the plane along the z axis and the starting point.Since there are N 2 rays, the complexity of the approaching phase is O(N 2 ).Escaping the plane's surface requires more work and the following analysis holds for each of the N 2 rays that need to be cast.Let r 0 denote the starting point in the vicinity of P 's surface.Note that D P (r 0 ) is about 1  N in size immediately after the polygonization routines for P 's cells have been invoked (at most two in this case, only one containing polygons).Analyzing the iteration (3), it is easy to see that D P (r 1 ) = 2•D P (r 0 ) and in general the following holds: i.e., the method escapes the surface in steps of exponentially increasing size.If D P (r 0 ) is about 1 N , then the number of steps required to exit the polygonization volume is O(log N ).Given that there are N 2 such rays, the complexity of the escaping phase is O(N 2 log N ).The approaching and escaping phase are performed sequentially.Thus, the complexity of polygonizing a plane in this scenario is O(N 2 log N ).
Case #2: plane P and rays are parallel.First, notice that if the distance between a ray and P is equal to k N in this case, then the method exits the polygonization volume in approximately N k steps.There are N rays marching through cells that contain P .The algorithm takes approximately N steps along each of these rays before terminating.Next, notice that there are N rays above and N rays below P , parallel to P and of distance approximately k N to P for some integer k < N .These rays require about N/k steps before exiting the polygonization volume.Thus, a conservative estimate for the number of steps S for all N 2 rays is where H N is N th partial sum of the harmonic series [22]: The number H N is about as large as log N .The reason for this comes from the comparison of H N and the integral Case #3: the general case.Due to easier exposition and without loss of generality, we assume that the plane passes through origin (i.e., r P = 0).Combining this assumption with equations (3) and ( 5), we get the following iteration for the z coordinate: where (n x , n y , n z ) T is the unit normal of the plane and (x 0 , y 0 , z 0 ) T is the origin of the ray.Let us denote with z * the intersection of the ray and the plane: Without loss of generality, we assume that n z < 0. There are two sub-cases: (1) the method approaches the plane along the ray and (2) the method moves away from the plane along the ray.In the first subcase, we have n x x 0 + n y y 0 + n z z n > 0. In this scenario, it is easy to see that Since 1 + n z is between 0 and 1, we have that the number of iterations n has to be about O(log N ) so that D P (r n ) becomes less than √ 6 2N (at which point the polygonization routine is invoked and we can move to the other side of the plane).In the second sub-case, we have n x x 0 + n y y 0 + n z z n < 0. Now the following holds: iterations along the ray are needed to exit the polygonization volume.Given that there are N2 rays in total, the complexity of case #3 is O(N 2 log N ).

Polygonizing rectangular boxes
A rectangular box can be obtained by intersecting six axis-aligned planes.Let d 1 , d 2 , . . ., d 6 be the distances from point (x, y, z) to each of these planes.Then the distance to the box is bounded by Since polygonizing each of the box sides takes O(N 2 log N ) steps, this is also the total complexity of polygonizing a box.
This can also be justified by the fact that the max operation partitions the polygonization volume into several regions.In each of these regions only the distance to one particular plane is relevant (largest d i ).All the rays passing through this region require at most O(N 2 log N ) steps before exiting the region.The conclusion about the total complexity follows from the fact that the number of such regions is finite.
As noted, Equation ( 12) bounds the distance to the box.The exact distance function can be constructed and this leads to more efficient marching in practice (a constant speed-up, not in the asymptotic sense) 2 .

Polygonizing other shapes
A shape can be approximated by K axis-aligned boxes.See Figure 2 for an illustration of this process.Of course, we can improve the quality of approximation by increasing K.It is important to note Figure 2: One possible approximation of a shape (red) as a union of axis-aligned boxes (blue).Note that the approximation would be much more efficient if non-axis-aligned planes are used on the boundary.
that K does not depend on grid resolution N .The efficiency of approximation can be increased by using non-axis-aligned planes at the boundary of the shape.This process is not unlike the use of triangle meshes in modern computer graphics.
Approximating the shape as a union of K boxes keeps the O(N 2 log N ) polygonization complexity.This is because the union of K boxes (and, in general, shapes) can be obtained by applying the min operation to combine all the individual distance bounds.The number of steps the method has to make in this case is asymptotically no worse than polygonizing each box on its own.Thus, the total number of steps scales as O(N 2 log N ) since the grid resolution N does not depend on K.

Experimental analysis
In this section we experimentally compare the three methods3 : gridhopping (ghop), enumeration (enum) and continuation (cont).The goal is to show that the enum method is asymptotically slower than cont and ghop.
The polygonization of a cell is obtained with the Marching cubes algorithm.To achieve this, we implement all methods and run them on different scenes for varying grid resolutions.It is important to note that all implementations produce exactly the same meshes when polygonizing signed distance bounds.

Primitives and simple shapes
We use four different scenes in this batch of experiments.The first scene contains 7 basic primitives: sphere, cube, cone, cylinder, torus, hexagonal prism and capsule.All these primitives have simple and efficient signed distance bounds.The second scene is a surface of genus 2 given by the implicit equation 2y( The third scene contains a knot with an explicit signed ditance bound.The fourth scene contains the Sierpinski tetrahedron.These scenes are visualized in Figure 3. Figure 4 shows the times needed to polygonize the scenes with the slow and the fast algorithm.Note that the axes in the graphs have logarithmic scale.We can see that the measured times for both methods appear as lines.This is expected since the computed theoretical complexities are polynomial (∼ N 3 and ∼ N 2 ).However, the fast method becomes significantly faster for large N , i.e., asymptotically.This aligns with the predictions from Section 4.

Experiments on learned SDBs
There has been considerable interest in learningbased methods for signed distance bound modelling.And this is lately especially true in the area of deep learning [3,18,5,8,6,26,27,15,21] 4 .
Let us constrict our analysis to a set of methods that use a neural network (NN) to approximate a SDB of a shape: NN θ (x, y, z) ≈ d S (x, y, z), where (x, y, z) ∈ R 3 is a point in 3D space and d S (x, y, z) denotes its Euclidean distance to the surface of the shape S. The parameters of the N N are denoted with θ.These parameters have to be tuned to make the approximation as best as possible.
Our hypothesis is that gridhopping is a faster algorithm for transforming such neural SDBs to triangle meshes than by using the basic enumeration technique.This should especially be true for higher grid resolutions.
To experimentally test our claim, we take six shapes from the Thingi10K dataset [28] (see Figure 5).These shapes are represented as triangle meshes and we need their SDBs.We first sample a dataset of point-distance pairs around each shape5 : An example can be seen in Figure 6.Next, we use the Adam stochastic gradient descent technique to tune θ (from Equation 13) so that NN θ (x i , y i , z i ) ≈ d i .The NN architecture is very simple: eight layers, each with embedding dimension equal to 64, ReLU is set as the activation function.To speed up training, we also use batch normalization in all layers except the last one.This results in about 30, 000 network weights (120kB of memory).By modern standards, this is a tiny network, but it is capable of accurately representing the shapes used in our experiments.Note that our methodology trivially extends to other NN approaches for SDB modelling.We decided to go with this one to keep things as simple as possible.
Recall that if certain criteria are not met (e.g., Lipschitz continuity), sphere tracing might "miss" parts  of the surface of the shape and we get an incomplete mesh with gridhopping [9].In practice, this means that the NN overestiamtes the surface distance for some points (e.g., due to the imperfection of the learning algorithm).To circumvent these problems, we shrink the signed distance approximation by factor λ: SDB(x, y, z) = NN θ (x,y,z) λ .This procedure slows down gridhopping and has no effects on the speed of the complete enumeration algorithm.We empirically found that λ = 1.5 works well.This (or smaller) value results in perfect mesh reconstructions for all shapes, i.e., all methods produce identical meshes.
We compare the times needed to polygonize the prepared SDBs for the complete enumeration algorithm and gridhopping.The grid resolution N varies from 64 to 512.First, all computations are performed on Intel(R) Core(TM) i7-9750H CPU @ 2.60GHz.Batching the NN evaluations in chunks of size ≈ 10, 000 leads to significant speed improvements for both methods.The results can be seen in Figure 7.Even though the enumeration method is sometimes faster for small grid resolutions, we can clearly observe that gridhopping has a significant asymptotic advantage.The resuts are in accordance with the theoretically predicted ones 4: ∼ N 3 vs.∼ N 2 computational complexity.
We also performed experiments when computations are performed on an Nvidia GeForce RTX 2060 Mobile GPU.Batching the NN computations in chunks of size ≈ 100, 000 seems best and we report timings in this setting, Figure 8.The theoretical complexity results hold in this setting as well.Also, GPU computations are significantly faster than CPU ones.However, gridhopping is somewhat slower for smaller grid resolutions.We attribute this fact to constant overhead needed to prepare the gridhopping run.

Conclusion
In this work we described and evaluated gridhopping -a method for transforming shapes defined by signed distance bounds into triangle meshes.Its main advantages are simplicity, ease of implementation and processing speed.We confirmed this in theory and experiments by comparing the method to the enumeration and continuation variants of the Marching cubes algorithm [14] which are the standard approaches for polygonizing signed distance bounds in practical applications.We used both simple shapes defined by short equations/code and more complicated ones defined through neural networks.We also showed that these conclusions hold when the computations are performed on the CPU and on the GPU, which is a further confirmation of the practicality of the method.
The main limitations of gridhopping are the Lipschitz continuity requirement, that it works only with abstract, mathematically defined, shapes and not real-world data like point clouds, and the fact that it generates too many triangles for some simple shapes.Note that the last limitation also holds for Marching cubes as both methods produce identical triangle meshes.Also, within the context of this paper, we would like to mention that the Marching cubes-based Figure 7: CPU experiments: elapsed times plotted against grid resolution for the six different shapes from Figure 5.The legend for all graphs is plotted in the top left one.Note that all axes are logarithmic.Figure 8: Same setting as in Figure 7, except the computations are done on the GPU.

Figure 3 :
Figure 3: The four scenes used in our experiments from Section 5.1.

Figure 4 :
Figure 4: Elapsed times plotted against grid resolution for the four different scenes from Figure 3.The legend for all graphs is plotted in the top left one.Note that all axes are logarithmic.

Figure 5 :
Figure 5: Shapes used in experiments from Section 5.2.

Figure 6 :
Figure 6: An example point cloud generated from a 3D model of a frog.