The medical image data acquired using tomography imaging machines such as CT and MRI typically forms a regular distribution data field in 3D space, which is also referred to as volume data. Three-dimensional surface reconstruction from the volume data is usually performed for visualization, feature computing and computer-aided design. The cost of rendering and manipulating surface models is much less than that of volume models, and the surface rendering process can be highly accelerated with modern graphic hardware which is critical in real-time application. Therefore, surface model is a widely used tool for data visualization [1]. In addition, surface model is necessary for some modern industrial manufacture such as 3D printing.

Several algorithms have been proposed on surface reconstruction. Ekou [2] developed a contour-connecting method in which closed contours in each slice are first generated, then the whole surface is reconstructed by connecting the contours between adjacent slices. Yang [3] proposed an SMC method which connects the triangular facets in each cube through estimating state value of each vertex of the cube. Among all the surface reconstruction algorithms, the marching cubes (MC) algorithm proposed by Lorensen in 1987 is the most commonly used one due to its simplicity on understanding and implementation [4]. The algorithm proceeds through all elements (voxel) in the volume data, taking eight neighbor voxels at a time to define a cube, then computing the polygons which are the local surface in the cube. After processing all cubes, the polygons are merged to form a whole surface.

The original MC algorithm needs scan all cubes in the volume data, and much time is wasted on empty cubes where no polygon exists. This promotes some improved versions of MC algorithm aiming to accelerate the surface generating process by ignoring empty cubes. One of those is regional growth MC (RGMC) algorithm, which first selects a seed cube exactly intersects the target surface, then extract polygons in the cubes adjacent to the seed. These new cubes are queued as new seeds and the operations are repeated on these new cubes until the queue becomes empty. Since only cubes intersect the target surface would be considered during this process, the time is highly saved.

Parallel processing is another powerful way to speed up the extraction process by dividing the computation into several parts and process each part with a single core in CPU. This promotes the work in this paper. Among all the parallel computing frameworks, message passing interface (MPI) has become a de facto standard for communication among processes that model a parallel programming due to its openness and portability. In this work, we choose MPICH as the implementation of MPI for its high performance and easiness to use.

In this paper, we combine parallel computing and RGMC to further accelerate the process of surface extraction. The volume data is partitioned into several parts, and a single seed is initially selected manually for each part. Then RGMC is applied on each part in a single process started using MPICH. After extracting, the whole surface is smoothed to make a better visual effect.

2.1Marching Cubes

Marching cubes is a classical algorithm of surface extracting in 3D data visualization field. The algorithm assumes that the original data is a discrete three-dimensional regular data field defined on a volume data. When a cube intersects the surface, the algorithm first computes the intersecting vertices of the surface and cube edges using linear interpolation technique, then connects the vertices to form triangles which are the approximation of the surface in the cube. After traversing all the cubes, a whold surface is extracted.

Let\( f\left( {x,y,z} \right) \) be a scalar function defined on the 3D volume data\( V = \{ {\text{v}}_{\text{i,j,k}} \}_{i,j,k = 0,0,0}^{M,N,K} \), with\( v_{i,j,k} \) to be voxel at\( \left( {i,j,k} \right) \). The MC algorithm divides the space of the volume data into a series of small cubes whose eight points are the voxels on each pair of adjacent slices along z direction. Mathematically, the surface corresponding to a given isovalue\( c_{0} \) is the point set defined in the following way [5]:

In discrete case, the surface is represented as a set of polygons which are the intersection of the surface and the cubes. Thus the basic question is how to determine whether a single edge of any cube intersects the surface. A simple answer to this question is: when the function value at one end of the edge is greater than or equal to\( c_{0} \) (+), and the function value at the other end is less than\( c_{0} \) (−), the edge intersects the surface. Since each point of a cube has two states (±), there are altogether 256 cases for any cube which forms a case table. There is a single index of any cube in the case table and from the index we can easily determine the states of the points in the cube. For more details, please refer to [4]. After finding the index of a cube in the case table, the coordinates of the vertices on an edge can be computed by the following linear interpolation process. Let\( {\text{e}} \) be the edge of a cube and\( S \) be the surface, the intersection of\( S \) and\( {\text{e}} \) is:

where\( {\mathbf{X}} \) is the intersection vertex of\( {\text{e}} \) and\( S \),\( {\mathbf{P}}_{{\mathbf{1}}} \),\( {\mathbf{P}}_{{\mathbf{2}}} \) are the two endpoints of\( {\text{e}} \),\( V1 = f\left( {{\mathbf{P}}_{1} } \right) \) and\( V2 = f\left( {{\mathbf{P}}_{{\mathbf{2}}} } \right) \) are the values of the scalar function\( \text{f} \) at the two endpoints. After obtaining the intersection of\( S \) and all the edges in each cube, the polygons in that cube can be determined considering the index of the states of the cube [6].

The normal vector at\( {\mathbf{X}} \) is also required in the surface rendering tasks, the normal vector can be computed at the same time when computing the coordinate of\( {\mathbf{X}} \). Central difference is used to calculate the normal vectors at each point of\( {\text{e}} \), and the normal vector of\( {\mathbf{X}} \) is computed by linear interpolating the normal vectors at\( {\mathbf{P}}_{{\mathbf{1}}} \) and\( {\mathbf{P}}_{{_{{\mathbf{2}}} }} \). For any voxel\( v_{i,j,k} = \left( {x_{i} ,y_{j} ,z_{k} } \right) \), the gradient of\( \text{f} \) at\( v_{i,j,k} \) is:

where\( \Delta x \),\( \Delta y \),\( \Delta z \) are the lengths of the edges of the cube respectively. Then the normal vector of\( S \) at\( v_{i,j,k} \) is:

After normalizing\( {\mathbf{g}} \) to a unit vector, the normal vector of\( S \) at\( {\mathbf{X}} \) can be computed using the following linear interpolation:

Where\( {\mathbf{N}}_{{\mathbf{1}}} = {\mathbf{g}}\left( {{\mathbf{P}}_{{\mathbf{1}}} } \right) \) and\( {\mathbf{N}}_{{\mathbf{2}}} {\mathbf{ = g}}\left( {{\mathbf{P}}_{{\mathbf{2}}} } \right) \) are the normal vectors of\( S \) at\( {\mathbf{P}}_{{\mathbf{1}}} \) and\( {\mathbf{P}}_{{\mathbf{2}}} \) respectively.

Marching cubes algorithm flow:

1. 1.

   Read the three-dimensional discrete data field and construct state table;

2. 2.

   Scan each two layers of the data to construct cubes whose eight points are taken from two adjacent layers;

3. 3.

   Scan all cubes. For each cube, the function value of each point of the cube is compared with the given isovalue of the surface, and the values are used to determine the index of current cube in case table;

4. 4.

   Compute the intersection of the cube edge and the surface by linear interpolation method using (2) and generate polygons in the cube;

5. 5.

   Compute the normal vectors of intersection vertices using (7);

6. 6.

   Combine polygons in all cubes into a whole surface mesh.

2.2Regional Growth MC Algorithm

The classic MC algorithm traverses all the cubes in the volume data, however, in general most cubes do not intersect the surface thus much time is wasted on empty cubes. To solve this problem, Kai [7] proposed a regional growth MC algorithm. First, a seed cube, which is a cube intersects the surface, is searched among all the cubes. The seed cube is immediately pushed into a queue and marked as visited. Then a loop process is applied on the queue. While the queue is not empty, the head cube will be popped out to compute the polygons inside it. Then its neighbors that have not been marked as visited and intersect the surface are pushed into the queue. This pop-compute-push process is repeated until the queue becomes empty. This way, traversing empty cubes is avoided and much time is saved. The workflow of RGMC is shown in Fig. 1.

Regional growth MC algorithm flow chart

Full size image

2.3MPI Programming

As a de facto standard, MPI has been widely accepted in high-performance computing systems [8,9] for it provides a language independent platform. It provides bindings with C, C++ and Fortran languages which makes it a feasible tool in programming practice. MPI has the advantages of good portability and high efficiency, and can be used in heterogeneous environments. There are many different free, efficient and practical versions of MPI [10].

In most parallel programming, there are six most commonly used basic interfaces of MPI, and we only need four of them [11,12].

In this section, we will explain the details of our method of surface extracting from volume data, combined RGMC and MPI. The process is based on VTK/ITK pipelines and can be divided into the following steps: read, reconstruct, post-process and output.

The flow chart of this article is shown in Fig. 2.

A pipeline for three-dimensional reconstruction

Full size image

3.1Read Data

Digital Imaging and Communications in Medicine (DICOM) is a standard for handling, storing, printing, and transmitting information in medical imaging. DICOM is widely used in radiotherapy, cardiovascular imaging and radiation diagnostic equipment (X-ray, CT, nuclear magnetic resonance, ultrasound, etc.). It has a large data dictionary, which contains almost all the medical environment under the common data, and can be a complete description of a variety of medical equipment, image format data and patient-related information.

Many software can read DICOM files, of which ITK and VTK are taken as the best combination for they can not only read and write DICOM files, but also provide a lot of useful medical image processing functionalities. That makes ITK and VTK a widely used set of tools for manipulating medical images and creating medical applications.

In this work, we use ITK to read DICOM data and use itk::ImageToVTKImageFilter to prepare data for VTK.

3.2Reconstruct

For processing the volume data in parallel, the volume data is divided into multiple parts and the RGMC procedure described in Sect. 2.2 is applied on each part separately, with a single process created and assigned to each part using MPICH.

There are two basic MPI parallel programming modes: peer-to-peer and master-slave. We use peer-to-peer mode because each process is responsible for the processing of different parts of the same data space.

MPI will start multiple processes and call the following MAIN procedure in each process. The number of processes and the id of current process could be extracted from the arguments passed to the MAIN procedure, based on MPI mechanism. Letimwidth andiheight be the width and height of the volume data along x-y direction, the reconstruct process is shown in the following algorithm:

In VTK, the classic MC algorithm is implemented in class VtkMarchingCubes. Since there is no RGMC implementation in VTK, we must extent VTK to realize this functionality. We wrote a class named vtkRGMarchingCubes, in which most code comes from VtkMarchingCubes except that some functions are modified to implement the parallel RGMC algorithm.

3.3Post-Process and Output

Before outputting the result surface, some further operations are needed for achieving a better visual effect. If the CT slices along the z-axis are not fine enough, the surface obtained will produce terraces (see Fig. 3(a)).

Meshes before and after smoothing. The original surface has terrace effect (a), and a surface smoothing process is applied on the surface using vtkSmoothPolyDataFilter (b).

Full size image

In order to solve this problem, we use the vtkSmoothPolyDataFilter to smooth the result data, which is a filter that adjusts point coordinates using Laplacian operation. We took the default settings of vtkSmoothPolyDataFilter, of which the iteration number is 20 and relaxation factor is 0.01. The smoothing result is shown in Fig. 3(b).

Since the meshes in different part are generated separately, there might be topology inconsistency along the common boundary of adjacent meshes (see Fig. 4(a)). This topological problem can be corrected using triangular refinement techniques (see Fig. 4(b)). For two parallel regions, look at their demarcation line (z-axis at k) for isolated nodes (e.g. point\( D_{1} \)). Then connect\( D_{1} \) to a point that is not connected, which belongs to a triangle that has a common edge on another region.

Topological inconsistency along common edges of meshes (a) and this problem can be solved using triangular refinement techniques

Full size image

4.1Programming Environment

• CPU: Intel(R) Core(TM)2 Quad CPU Q9400 @2.66 GHz 2.67 GHz

• RAM: 5.00 GB

• Operating System: Windows7 64bit

• Application software: ITK VTK

• Development environment: Microsoft Visual Studio 2013

4.2Results and Discussion

Three cases (skull, femur, forearm bone) were tested. The comparison of running time using MC, RGMC and our method 2 processes is shown in Table 1. As can be seen from the table that our method overcomes the other two significantly.

Table 2 lists the experimental results under different number of processes.

Figure 5(a) (b) (c) is the result of the RGMC algorithm, Fig. 5(d) (e) (f) is the result of our method with 2 processes. The reconstruction results of each process is colored in red and gray respectively.

The reconstruction result of RGMC and our method (Color figure online)

Full size image

4.3Evaluation

In order to assess the similarity between our reconstructed results and the RGMC reconstruction results, we used the hausdorff distance.

The hausdorff distance is a measure of the degree of similarity between two sets of point sets, which is a definition of the distance between two sets of points:

Assuming that there are two sets of sets\( A = \{ {\text{a}}_{ 1} ,\ldots , {\text{a}}_{\text{p}} \} \),\( {\text{B}} = \{ {\text{b}}_{ 1} ,\ldots {\text{b}}_{\text{q}} \} \), then the hausdorff distance between the two points is defined as:

where:\( {\text{h(A,B)}} = \mathop {\hbox{max} }\limits_{a \in A} \mathop {\hbox{min} }\limits_{b \in B} ||a - b|| \),\( {\text{h(B,A)}} = \mathop {\hbox{max} }\limits_{b \in B} \mathop {\hbox{min} }\limits_{a \in A} ||b - a|| \),\( || \bullet || \) is the distance between the A and B points (e.g., L2 or Euclidean distance),\( {\text{h(A,B),h(B,A)}} \) are called the one-way hausdorff distance from A set to B set and from set B to set A, respectively. That is,\( {\text{h(A,B)}} \) actually first sorts the distance\( ||a_{i} - b_{j} || \) between each point\( a_{i} \) in point set A to the point\( b_{j} \) closest to this point\( a_{i} \), and then take the maximum value in the distance as the value of\( {\text{h(A,B)}} \).\( {\text{h(B,A)}} \) is equally available. From the formula (1), the bidirectional hausdorff distance H (A, B) is the larger of the two-way distance h (A, B) and h (B, A), which measures the interval between two points maximum degree of mismatch.

Here we evaluate the results of these three datasets separately are 0.79511, 0.59136, 0.601142.

A parallel regional growth MC algorithm is presented which significantly overcomes the classic and regional growth MC algorithms in time efficiency. The scalar field is partitioned into several parts and the RGMC is applied for each part in different processes using MPI. Experiment illustrates that about 50% of time will be saved compared with classical MC algorithm and RGMC algorithms.