rupMC: a ray-unit parallel marching cubes algorithm on CPU/GPU heterogeneous architectures

ABSTRACT

The marching cubes (MC) algorithm is widely used for extracting isosurfaces from volume data and 3D visualizations because of its effectiveness and robustness but require extensive memory and computing time for large-scale applications. Additionally, MC isosurfaces lack topologic information, making them difficult to use in some geologic applications. To overcome these limitations, this study proposes an enhanced MC using CPU/GPU heterogeneous architecture called the ray-unit parallel MC (rupMC) algorithm. First, ray units form the basic voxel to determine how the surface intersects to reduce repeated computations and enhance efficiency. Then, rupMC uses multiple computing processes and threads on a CPU/GPU heterogeneous architecture to process points concurrently. Finally, the unique surface intersection indices are preserved to compose the surface triangles, and the topological surface information is directly embedded in the triangle compositions. Experiments on five stratum datasets of varying sizes demonstrated that, rupMC achieved approximately dozens of times faster than other serial MC and 4 times faster than a parallel DMC. rupMC demonstrated high scalability and adaptability to various CPUs/GPUs and datasets of various sizes. rupMC has remarkable capabilities for efficiently and feasibly extracting precise surface intersections and triangles, making it well-suited for large-scale and high-density applications.

KEYWORDS:

• Open-source Message Passing Interface (OpenMPI)
• Compute Unified Device Architecture (CUDA)
• CPU/GPU heterogeneous architecture
• Iso-surface extraction
• marching cubes (MC)

1. Introduction

Reconstructing a surface from a series of geological sampling points is a well-explored problem in geological modeling and computer graphics (Zhang, Zhong, and Wang 2022), such as the visualization of geological models (Li, Zhong, and Wang 2021) and the simulation of rock microstructures (Zhou and Xiao 2017). In geological modeling, an implicit modeling method is employed to interpolate geological sampling points and derive an implicit function that represents the optimal geological surface. To effectively visualize the interpolated implicit function, it is crucial to use the implicit surface reconstruction method to extract the target iso-surface of the interpolated implicit function.

Numerous implicit surface reconstruction methods are used to extract the isosurface, which can be divided into two categories (Chin et al. 2021). In the first category of surface rendering methods, a set of pertinent surfaces is first extracted from a given dataset. Subsequently, these surfaces are polygonized and rendered for display using conventional graphics hardware. This type of method includes the marching cubes (MC) (Lorensen and Cline 1987; Wang et al. 2021; Zhang, Zhong, and Wang 2022), marching tetrahedra (MT) (d’Otreppe, Boman, and Ponthot 2012; T. Zhao et al. 2021), and Delaunay triangulation methods (Jamin et al. 2015; Peiró, Sherwin, and Giordana 2008). In the second category of volume rendering methods, the volume is sampled directly with a set of rays originating from the eye viewpoint. One of the most popular volume-rendering methods is the ray-casting algorithm (Ling and Qian 2011).

As a classical surface-rendering method, the MC algorithm has been widely used owing to its effectiveness and robustness (Lorensen 2020). The MC algorithm triangulates sampling points by evaluating their corresponding functional values. First, regularly spaced evaluation points are sampled with a specified voxel size in the modeling space. Subsequently, based on the eight vertices of each cube and a predefined lookup table, the sampling points are triangulated to approximate the target isosurface. However, the MC algorithm also has shortcomings, such as the massive memory required, large time cost, and lack of topologic information. Numerous studies have been conducted to address these problems.

Montani, Scateni, and Scopigno (1994) alternated the midpoint of the edge to the surface intersection calculated by MC, which is able to reduce floating-point operations and rapidly extract surface triangles. Du et al. (2012) marked 15 basic intersection topologies that indicated which adjacent cubes needed to be visited to avoid traversing the entire data area and improve computing efficiency. Wang et al. (2021) proposed an MC algorithm based on edge growth to extract the 3D contours of parts connected to the seed. Zhang, Zhong, and Wang (2022) proposed an improved signed MC method that consists of rough and exact evaluation steps. The basic strategy of this method is to filter the invalid voxels that do not intersect with the target iso-surface in a low-cost manner, and to evaluate the valid voxels that intersect with the target iso-surface accurately. In these studies, the computational performance was improved by reducing the number of computations. However, the search for and precise calculation of surface intersections are overlaid on these methods, making it difficult to maintain the accuracy of the surface intersections.

With improvements in computer devices, in addition to improving the computing efficiency based on the algorithm itself, some enhanced algorithms based on improved hardware have been proposed (Yang et al. 2022; Zhou et al. 2018; Zhou, He, and Qiu 2017). Wei et al. (2017) parallelized the extraction process of the regional growth MC algorithm using a message-passing interface on a Central Processing Unit (CPU). Gao, Tang, and Wu (2015) presented a novel Graphics Processing Unit-based (GPU) parallel MC algorithm that created a GPU-based uniform grid structure to manage point clouds. D’Agostino and Seinstra (2015) presented a parallel implementation of the MC algorithm that can exploit multiple GPU devices. Chin et al. (2021) proposed a GPU-accelerated enhanced MC 33 algorithm and Grosso and Zint (2022) used a half edge data structure in a GPU-based parallel dual-MC approach. In these studies, computational efficiency was improved on a CPU or GPU. However, the computations of the basic cube units in MC are independent, making them suitable for parallelization on CPUs and GPUs. The CPU/GPU heterogeneous architecture can make full use of multiple CPU processes and GPU threads to further improve the computing efficiency.

Few methods have been presented for extracting the topologic relationships of an isosurface from the topologic information of the surface triangles. However, it is important to directly show the topological connections between the surface triangles when applying them to geological modeling. For example, when reconstructing implicit surfaces for geological models, repairing the raw structural topological relationships of geological bodies is important (Li, Zhong, and Wang 2021).

This study proposes an enhanced MC algorithm on a CPU/GPU heterogeneous architecture called ray-unit parallel MC (rupMC), which addresses the three aforementioned problems. First, a ray unit is used for the marching computations, which consists of a point and three rays along the three directions of the data composition. The lookup index for the ray unit is still created based on the state of the eight vertices of a cube and serves as a pointer in the precalculated table that provides all edge intersections for a given cube configuration. Edge intersections are preserved and calculated only once in the corresponding ray unit. Second, the algorithm is parallelized on a CPU/GPU heterogeneous architecture. The Input/Output (I/O) process is parallelized by multiple CPU processes. Once a sub-dataset is input, it is transmitted to another CPU process and computed on a GPU. Simultaneously, the next sub-dataset is input. Such pipelined task distribution saves time by reading large amounts of data. The computations are also parallelized using multiple GPU threads to further improve computational efficiency. Furthermore, an adaptive domain decomposition technique is implemented to enable the efficient processing of datasets of any size across various CPU/GPU environments with varying computing capacities. Third, based on the ray unit, the unique indices of the intersections are preserved to compose the surface triangles, and the topological information of the surface is directly embedded in the triangle compositions.

The experiments show that rupMC reduces the number of repeated final products and generates the same surface triangles with a unique intersection index. rupMC significantly reduces computation time and maintains the topological information of the surface triangles. The proposed rupMC demonstrates high scalability and adaptability to various CPUs/GPUs and datasets of various sizes. rupMC has remarkable capabilities for efficiently and feasibly extracting precise surface intersections and triangles, making it well-suited for large-scale and high-density applications.

2. Method

2.1. Introduction of the marching cubes algorithm

To improve our understanding of rupMC, we present an introduction to MC and its limitations. The MC algorithm comprises three main steps (for more details, please refer to Lorensen and Cline 1987).

1. A cube is created using eight points from two adjacent slices of the sampling points at different heights (as shown in Figure 1). A state of one is assigned to a cube vertex if the corresponding data value at that vertex exceeds (or equals) the value of the isosurface, indicating that these vertices are inside (or on) the surface. In contrast, cube vertices with data values below the surface are assigned states of zero and are outside the surface. We assume that the surface intersects the cube edges, where one vertex has a state of one (inside the surface) and the other has a state of zero (outside the surface). A lookup index is created for the cube based on the states of the vertices.

2. Each cube consists of eight vertices, with two states: inside and outside. Consequently, there are only 256 (2⁸) possible ways for a surface to intersect a cube. These 256 cases are generated using complementary and rotational symmetry from the 14 basic patterns (as shown in Figure 2). By enumerating these cases, a lookup table is created to identify surface–edge intersections, given the labeling of the cube vertices. The lookup index serves as a pointer in the table that provides all the edge intersections for a particular cube configuration.

3. The lookup index is used to determine which edge the surface intersects, the locations of the intersections are calculated by linear interpolation, and surface triangles are created. Once the algorithm determines the way the isosurface intersects the cube, it moves (or marches) to the next cube.

Figure 1. Composition of a cube.

Figure 2. Basic patterns of a surface intersecting a cube.

The first issue with MC is associated with the cubic unit used for basic computing. Owing to the regular arrangement of the cubes, each cube edge is shared by four adjacent cubes (Figure 3). Using cubes as the basic marching unit for the calculation results in four repeated computations and preservation of the surface intersections, wasting computing time and memory. Thus, a new unit for marching computing is necessary to reduce redundant computations and results.

Figure 3. Cube edge shared by four adjacent cubes.

The second issue with MC is its high I/O and computational intensity. The MC algorithm produces an accurate surface when high-density sampling point data are used but requires large reading and writing times. In particular, for some real-time visualization applications, the extensive I/O time is a significant limitation. In addition, the computations for each marching unit are similar and independent of the computations for the other cubes. Therefore, a different approach is necessary to improve the I/O and computational efficiency.

The last issue with MC is the preservation of surface intersections and triangular compositions. MC is used not only for data visualization but also for some specific spatial analyses, such as the reconstruction of raw structural topological relationships between geological bodies, which require topologic details. Therefore, an approach is required that preserves the topological information of the isosurface triangles.

2.2. rupMC: a ray-unit parallel marching cubes algorithm

To address the three aforementioned issues of MC, this study proposes a ray unit parallel MC algorithm called rupMC. As shown in Figure 4, the three key improvements in rupMC are as follows.

1. A ray unit, which consists of a pixel and three rays along the three data composition directions, is used for the marching computations. The lookup index for the ray unit is still created based on the state of the eight cube vertices and serves as a pointer in the predefined table that provides all edge intersections for a given cube configuration. Edge intersections are preserved and calculated only once in the corresponding ray units.

2. The algorithm is parallelized using a CPU/GPU heterogeneous architecture. The I/O process is parallelized by multiple CPU processes, and the computations are parallelized on GPUs to further improve the computational efficiency. Furthermore, an adaptive domain decomposition technique is implemented to enable the efficient processing of datasets of any size across various CPU/GPU environments with varying computing capacities.

3. Based on the ray unit, the unique indices of the intersections are preserved to consist of surface triangles, and the topological information of the surface is directly embedded in the triangle compositions.

Figure 4. rupMC flowchart.

2.2.1. Ray-unit computing units

A cube is the basic computing unit in the MC algorithm and is created from eight points, four each from two adjacent slices (as shown in Figure 1). Owing to the regular arrangement of the cubes, each cube edge is shared by four adjacent cubes (Figure 2). Thus, using cubes as the marching unit to identify and calculate edge intersections will require four repeated computations and preservations, wasting computing time and memory space. Thus, a new unit for marching computing is necessary to reduce redundant computations and results.

By simplifying the cube unit, the ray unit becomes sufficient to represent a three-dimensional dataset without overlapping edges (as shown in Figure 5). The ray unit consists of a point and three rays along the three directions of the data composition. Thus, each point has a corresponding ray unit. The lookup index for the ray unit is still created based on the state of the eight cube vertices and serves as a pointer in the predefined table that provides all edge intersections for a particular cube configuration.

Figure 5. Ray-unit of rupMC.

In rupMC, two sets of marks are created: one for recording the count of surface triangles in the cube, and the other to mark whether there are surface intersections on each edge of the ray unit. These marks are used for the calculation and preservation of the surface intersections (details are described in Section 2.2.3).

2.2.2. Parallelization on CPU/GPU architectures

MC calculates the edge intersections and surface triangles using the raw point values and marks generated in each ray unit. The computations of each ray unit are independent of those of the other units, which implies that the computational processes can be parallelized. rupMC takes advantage of the parallel computing capability of the CPU and GPU and uses multiple computational processes and threads to improve computing efficiency. The I/O processes of rupMC are parallelized based on the Message Passing Interface (OpenMPI), which is a high-performance message passing library for communication in CPUs (https://www.open-mpi.org/). The computing processes of rupMC are parallelized by the Compute Unified Device Architecture (CUDA), a parallel computing programming model for general computing on GPUs (https://developer.nvidia.com/about-cuda).

The parallel-computing framework of rupMC on a CPU/GPU heterogeneous architecture is shown in Figure 4. First, the parameters of the data are input into the algorithm by CPU Process 0. rupMC divides the entire dataset into multiple sub-datasets adaptively based on the number of CPU processes such that each sub-dataset can be accommodated and processed by each process. Each subdata point is then input and transmitted from Process 0 to the corresponding computing processes.

Second, each process receives a sub-dataset from Process 0 and adaptively divides the entire sub-dataset into multiple subdomains based on the available memory space of the GPU. Subsequently, the data of each subdomain are transmitted from the CPU memory to the GPU memory for processing. Each GPU thread is responsible for the data input and all computations of a ray unit. Thus, multiple GPU threads generate lookup indices and marks for multiple points simultaneously. Subsequently, the edge intersections and surface triangle constructions are calculated using GPU threads. Next, the results are transmitted from the GPU memory back to the CPU memory. Subsequently, all loops on the GPU are completed and the sub-results are transmitted from the corresponding processes to Process 1.

Finally, Process 1 receives the sub-results and writes them into the final files.

2.2.3. Index of intersections

The MC calculates the surface intersections and triangular compositions in each cube whose edges are shared by the four adjacent cubes. Thus, the locations of the intersections are repeated, and it is difficult to show the topological relationships between the surface triangles.

rupMC substitutes the ray unit for a cube and reduces the number of redundant computations. Whether intersections exist on the edges of a ray unit, and the number of surface triangles are marked. Thus, using the corresponding marks, rupMC only calculates the intersections at the edge of the corresponding ray unit, eliminating repeated computations and results. Surface triangles are composed of unique intersection indices; therefore, it is easy to determine the topological connections between surface triangles. For a detailed example, please refer to Figure 9 of our experiments (Section 3.2).

3. Experiments

3.1. Datasets and test environment

A set of stratum data (Figure 6) was used in the experiments to assess the performance of rupMC, which provided 102 × 151 × 702 soil moisture content values at a 15 m × 15 m × 1 m spatial resolution. To assess the accuracy and computational performance of rupMC better, the original data were interpolated into four datasets of different sizes (204 × 302 × 1404, 306 × 453 × 2106, 408 × 604 × 2808, and 510 × 755 × 3510).

Figure 6. Original stratum data, which provides 102 × 151 × 702 values at 15 m × 15 m × 1 m spatial resolution.

We compared rupMC to four other methods: the VTK-implemented MC, the VTK-implemented Flying Edges method (FE, Schroeder, Maynard, and Geveci 2015), C++-implemented MC, and a parallel Dual MC (DMC, Grosso and Zint 2022) method. MC is typically implemented using the Visualization Toolkit (VTK), which is open-source software for image processing and visualization (https://vtk.org). We reimplemented MC as a serial program using C++ to better evaluate the computational performance. The FE method is a high-performance iso-contouring algorithm for structured data, and is available in VTK. The parallel dual-MC approach generates meshes with better triangular and quadrilateral quality, which was implemented by CUDA on a GPU.

The experiments were conducted on a computer; the hardware and software environments are listed in Table 1. The VTK-implemented MC, VTK-implemented FE, C++-implemented MC, DMC, and rupMC were tested on Computer 1.

Table 1. Hardware and software environments on the computer.

Computer		Computer 1
CPU	Type	Intel Xeon Bronze 3204 CPU @ 1.90GHz
	Processors	12 × 6
	Memory	64 GB
GPU	Type	NVIDIA GeoForce RTX 2080 Ti
	Counts	2
	Kernels	4352
	Memory	11 GB
MPI	Version	10.1.2
CUDA	Version	11.0

3.2. Overall accuracy and visualization

The stratum datasets were extracted using the VTK-implemented MC, VTK-implemented FE, C++-implemented MC, DMC, and rupMC. The visualizations of the isosurface extracted by these methods were compared to assess the accuracy. As shown in Figure 7, the visualizations of the VTK-implemented MC, VTK-implemented FE, C++-implemented MC, and rupMC were the same. DMC used the MC algorithm to extract the isosurface and simplified these triangles, leading to fewer triangles but a rougher surface. The isosurfaces extracted by DMC were closed; thus, only the outermost isosurface was visible.

Figure 7. Visualizations of iso-surface extracted by five methods: (1) VTK-implemented MC; (2) VTK-implemented FE; (3) C++-implemented MC; (4) DMC; (5) rupMC.

The counts of the surface intersections (PTS in Figure 8) and triangles (TRA in Figure 8) obtained using these methods are shown in Figure 8. VTK is designed mainly for computer graphics and visualization. Thus, the MC algorithm in VTK used a simplified method to generate surface triangles for rapid visualization (Avila 2010). Therefore, the number of surface intersections and triangles generated by the VTK-implemented MC were lower, and their accuracies were lower than those generated by the MC. Compared with the C++-implemented MC, rupMC created the same surface triangles but with less intersection preservation, mainly because of the reduction in redundant computations. When the lookup index serves as a pointer in the table and gives all edge intersections for a given cube configuration, there are intersections on the edges of the ray unit, and the counts of the surface triangles are marked. Thus, rupMC calculated only the intersections on the edge of the corresponding ray unit, thereby reducing the number of repeated results. The results show that with these intersections and triangle marks, rupMC effectively reduced the number of repeated final products. The DMC used the face-coloring method to simplify quads that contain repeated intersections; thus, the number of surface intersections was low. However, the half-edge data structure in the DMC must be recomputed after the simplification step, which required many global memory accesses on the GPU and was time-consuming (Grosso and Zint 2022).

Figure 8. Counts of surface intersections (PTS) and triangles (TRA) by five methods.

The surface triangle compositions and point indices generated by the C++-implemented MC and rupMC are shown in Figure 9. Compared to the C++-implemented MC, rupMC generated intersections with a unique index, leading to a reduction in repeated intersections. Owing to the unique index of each intersection, it was easy to find a common edge between two surface triangles and the corresponding topologically adjacent relationships.

Figure 9. Surface triangle compositions and intersections generated by C++-implemented MC and rupMC. (1) Surface triangle compositions with repeated point indices generated by C++-implemented MC; (2) Surface triangle compositions with point index and corresponding locations generated by C++-implemented MC; (3) Surface triangle compositions with unique point indices generated by rupMC; (4) Surface triangle compositions with point index and corresponding locations generated by rupMC. The topologic adjacent relationships between surface triangles are shown.

3.3. Computational performance of rupMC

The overall running times of VTK-implemented MC, VTK-implemented FE, C++-implemented MC, DMC, and rupMC are listed in Table 2. The results showed that for all five datasets, rupMC on a CPU/GPU heterogeneous architecture achieved speedups of approximately 5 and 4 over the VTK-implemented MC and VTK-implemented FE, over 9 compared with the C++-implemented MC on a CPU, and over 2 compared with DMC on a GPU, indicating that rupMC significantly improved computational performance.

Table 2. Overall running time and speedups by VTK-implemented MC, VTK-implemented FE, C++-implemented MC, DMC and rupMC.

Size of dataset		Time (s)					Speedup
		VTK_MC	VTK_FE	C++_MC	DMC	rupMC	VTK_MC /rupMC	VTK_FE /rupMC	C++_MC /rupMC	DMC /rupMC
102 × 151 × 702	36.40MB	224	223	623	87	76	2.95	2.93	8.20	1.14
204 × 302 × 1404	294.00MB	1143	1096	2603	473	302	3.78	3.63	8.62	1.57
306 × 453 × 2106	1.00GB	3231	2952	6622	1495	729	4.43	4.05	9.08	2.05
408 × 604 × 2808	2.34GB	6935	5650	12815	4118	1392	4.98	4.06	9.21	2.96
510 × 755 × 3510	4.55GB	12708	10300	23812	6415	2397	5.30	4.30	9.93	2.68

The speedup was consistent across all five datasets, suggesting that rupMC is scalable for varying data volumes. In addition to the use of multiple processes on the CPU and multithreading on the GPU, the reduction of redundant calculations, along with the adaptive domain decomposition technique, aided in improving the computational performance. Moreover, adaptive domain decomposition enables rupMC to manage datasets of any size and to operate on various types of CPU/GPUs.

The computation times without the I/O process and the speedups of the five methods are listed in Table 3. The speedups of rupMC parallelized the computing processes by the multithread operation of the GPU, which significantly improved computational performance.

Table 3. Time and speedups of computing processes by five methods.

Size of dataset		Computing time (s)					Computing speedup
		VTK_MC	VTK_FE	C++_MC	DMC	rupMC	VTK_MC /rupMC	VTK_FE /rupMC	C++_MC /rupMC	DMC /rupMC
102 × 151 × 702	36.40MB	22	5	125	7	1	22.00	5.00	125.00	6.83
204 × 302 × 1404	294.00MB	138	29	478	11	2	69.00	14.50	239.00	5.60
306 × 453 × 2106	1.00GB	389	89	1129	24	5	77.80	17.80	225.80	4.85
408 × 604 × 2808	2.34GB	884	184	2353	51	20	44.20	9.20	117.65	2.57
510 × 755 × 3510	4.55GB	1646	350	4199	85	41	40.15	8.54	102.41	2.08

4. Conclusion and future work

This study proposes an enhanced MC (rupMC) on a CPU/GPU heterogeneous architecture to extract surface intersections and triangles, which has three key improvements. First, a ray unit is used in rupMC as the basic voxel for MC, consisting of a point and three rays along the three directions of the data composition. Second, rupMC uses multiple CPU processes for the I/O of data and multiple GPU threads for computing. This pipelined task distribution on a CPU/GPU heterogeneous architecture makes full use of the devices to process points concurrently. It automatically partitions the spatial domain of raw data into multiple subdomains for various GPUs with different computing capacities. Third, based on the ray unit, the unique indices of the intersections are preserved to compose the surface triangles, and the topological information of the surface is directly embedded in the triangular compositions.

In the experiments, the visualizations of the final products using rupMC were the same as those using VTK-implemented MC, VTK-implemented FE, C ++-implemented MC, and DMC. Compared with the C++-implemented MC, rupMC generated fewer intersections with a unique index and created the same surface triangles. The results showed that rupMC effectively reduced the number of repeated final products and made it easier to reveal the topological connections between surface triangles.

In addition, rupMC achieved approximately 50, 160, 10, and 4 speedups compared to the VTK-implemented MC, VTK-implemented FE, C++-implemented MC, and DMC, respectively, indicating significant improvements in computational performance. Through parallelization, rupMC significantly improved the computing performance and achieved hundreds of speedups for all five datasets over the C++-implemented MC. By using multi-process and multi-threading on CPU/GPU heterogeneous architectures and changing the basic computing unit, rupMC reduced redundant computations, and achieved approximately two speed-ups over the parallel DMC. In summary, the technique demonstrated its high scalability and adaptability for varying CPUs/GPUs and datasets of various sizes. rupMC has remarkable capabilities for efficiently and feasibly extracting precise surface intersections and triangles, making it well-suited for large-scale and high-density applications. The source code of rupMC is publicly available at https://github.com/HPSCIL/rupMC.

In our future work, rupMC will be adopted for applications such as surface reconstruction from LIDAR point clouds and real-time visualization of geological models.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

Data supporting the findings of this study are available in [rupMC] at [https://github.com/HPSCIL/rupMC].

Additional information

Funding

This study was supported by the National Natural Science Foundation of China (grant number 42171466).