A reliable surface reconstruction system in biomedicine

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 6 ( 2 0 0 7 ) 141–152

journal homepage: www.intl.elsevierhealth.com/journals/cmpb

A reliable surface reconstruction system in biomedicine Ying-Cheng Chen a , Yung-Chang Chen a,∗ , Ann-Shyn Chiang b , Kai-Sheng Hsieh c a b c

Department of Electrical Engineering, National Tsing Hua University, Hsinchu 300, Taiwan, ROC Institute of Biotechnology, National Tsing Hua University, Hsinchu 300, Taiwan, ROC Veterans General Hospital, Kaohsiung 813, Taiwan, ROC

a r t i c l e

i n f o

a b s t r a c t

Article history:

For common biomedical imaging facilities, such as CT, MRI, and confocal microscopy, the

Received 29 July 2005

acquired scans are sequential parallel sections. The object of interest in each section image

Received in revised form

can be extracted by segmentation procedure to form serial parallel planar contours. How

30 January 2007

to reconstruct a trustworthy surface from these contours is a crucial issue in biomedical

Accepted 30 January 2007

3D visualization. In this paper, we propose an automatic, fast, and reliable surface reconstruction system. An improved correspondence-determining algorithm is proposed in the

Keywords:

system to provide more reasonable contour-correspondences than the existing algorithms.

Surface reconstruction

It can handle more general input data, and does not produce wrong reconstruction results. A

Contour-correspondence

hybrid tiling algorithm is presented to tile the corresponding contours without the require-

Tiling

ment of a contour-matching procedure. It can also handle the branching problem without

Branching

any modification. For degenerate cases and branches, intermediate contours are introduced

Degenerate

by means of contour interpolation to enhance the reconstruction results. The surface area and volume are also calculated to facilitate the practical applications. © 2007 Elsevier Ireland Ltd. All rights reserved.

1.

Introduction

Computer technology advances very fast in the past two decades, which benefits biomedical imaging and facilitates its 3D visualization. We can reconstruct 3D objects and render them without altering traditional 2D imaging apparatus. 3D visualization techniques enable better understanding of the topology and shape of an object, and enable measurements of its geometrical characteristics. The extracted information is helpful in clinical diagnosis, surgical planning, and biological research. For common biomedical imaging facilities, such as CT, MRI, and confocal microscopy, the acquired scans are sequential parallel sections. Therefore, how to reconstruct a trustworthy surface from the sequential parallel 2D sections becomes a crucial issue in biomedical 3D visualization. Some previous methods solve this problem by estimating isometric volume data from original parallel section images. The isometric volume data are also composed of parallel

∗

sequential digital images, but the stack of images is more compact. Assuming original section images to be parallel to the xy-plane, the isometric volume data can be obtained by raster interpolation of adjacent images along z-axis. The object of interest in the volume data can be extracted by 2D or 3D segmentation procedures. For further refinement of the rugged boundary, the marching cubes algorithm [1] is usually adopted to construct a polyhedral model for representation of the object. This high resolution surface model comprises many tiny triangles that are about the same size as the voxels in the volume data. Algorithms of this category provide good reconstruction results, but are very time consuming. The other methods, including our approach, deal with this problem by extracting the object of interest in each original section image first. The object boundary in each image is assumed to be composed of several non-crossing simple closed contours. The object surface is reconstructed from these parallel planar contours directly. Meyers et al. [2] have described the

Corresponding author. Tel.: +886 3 573 1153; fax: +886 3 571 5971. E-mail address: [email protected] (Y.-C. Chen). 0169-2607/$ – see front matter © 2007 Elsevier Ireland Ltd. All rights reserved. doi:10.1016/j.cmpb.2007.01.011

142

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 6 ( 2 0 0 7 ) 141–152

main problems in the reconstruction process, including the correspondence, tiling, and branching problems. The simplest solution to the correspondence problem is to assign the correspondences manually [3], but it is not convenient. For automation, some methods deal with this problem by matching contour segments [4–6], others by evaluating overlaps of areas enclosed by the contours [7,8], and still others by constructing the areas of difference [9–11]. These methods either restrict the input data to the case without branching or handle specific data set only. They fail to provide reasonable results for some particular patterns of contour pairs. For the tiling problem, the greedy algorithm that tiles surfaces by local advancing rules is the most effortless method, but will probably introduce a self-intersecting surface while a contour and its correspondent in adjacent images are very different in shape, size, or orientation [3,12]. Some methods tile the surface by computing the optimal surface with predefined criteria [13,14]. These criterion-optimizing algorithms provide good reconstruction results but take a lot of computation time. Some hybrid approaches [4–6] are based on the two kinds of method mentioned above. A contourmatching procedure is applied to contours in adjacent images to obtain the similar contour segments. These contour segments are tiled by the greedy algorithm and the remainder by the criterion-optimizing algorithms. The hybrid methods provide reconstruction results almost as good and do not take as much time as the pure criterion-optimizing algorithms. Boissonnat [15,16] reconstructs object surfaces with a collection of tetrahedra defined from the Delaunay triangulation. The method solves the branching problem to some extent, but needs an extra correction scheme while three typical examples mentioned by Boissonnat occur. We propose an automatic system that provides a solution to the ill-posed problem of surface reconstruction from serial parallel planar contours. An improved correspondencedetermining algorithm is applied in the system to provide more reasonable contour-correspondences than the existing algorithms for some particular patterns of contour pairs. A hybrid tiling algorithm follows to tile contours without the requirement of a contour-matching procedure, and can also handle the branching problem without any modification. The reconstruction results of degenerate cases and branches can be further refined by means of contour interpolation. The system also calculates the surface area and volume of the object to facilitate practical applications. This paper is organized as follows. Section 2 presents an overview of our system, including system architecture and computation procedures. In Section 3, the proposed algorithms are described in detail, including contourcorrespondence analysis, surface tiling, rendering, and quantitative analysis. The implementation and experimental results are discussed in Section 4, and concluding remarks are made in Section 5.

2.

System overview

As mentioned in Section 1, the system takes sequential parallel planar contours as input to reconstruct the polyhedral model of an object, and then performs 3D rendering and

quantitative analysis for the polyhedral model. Fig. 1 shows the system flowchart. The proposed correspondence analysis algorithm is applied to contours in adjacent slices. If there is a one-to-one correspondence between two contours, they are tiled with a triangle strip by the proposed tiling algorithm directly. Otherwise, if the correspondence among contours is not one-to-one, i.e. the object has branches, an extra contour interpolation algorithm is applied before the tiling process to create intermediate contours. By tiling these additional contours, a better visualization result can be achieved. After going through all slices, a preliminary reconstruction result is accomplished. It is further refined by an algorithm of degenerate contour tiling. The reconstructed surface is smoother at the extremities of the object after the refinement. The final reconstruction result is represented by a polyhedral model. We render the polyhedral model to obtain 3D images of the object, and perform a quantitative analysis to obtain the surface area and volume of the object. All the procedures mentioned above will be described and illustrated in the next section.

3.

Methods

3.1.

Contour-correspondence analysis

3.1.1.

Initialization

A tree structure is constructed to describe the relationship among contours in each slice image of the image stack. The procedure is illustrated in Fig. 2 concisely. A contour with index k on the slice i is denoted as Ci,k . Ci,2 and Ci,3 lie in the inner part of Ci,1 , so Ci,2 and Ci,3 will be the children of node Ci,1 . Ci,4 lies in the inner part of Ci,2 , so Ci,4 will be the child of node Ci,2 . Here two types of contour, external and internal, are defined. A contour is external (internal resp.) if its immediate interior (exterior resp.) is the object’s material. In fact, the type of contour is different between parent and child nodes. Since Ci,1 is an external contour, its children, Ci,2 and Ci,3 , are internal contours. As a child of node Ci,2 , Ci,4 is an external contour. Once the tree has been built, it can facilitate the following procedures.

3.1.2.

Improved correspondence-determining algorithm

In general, the automatic correspondence-determining algorithms can be categorized into two main groups: overlapbased approaches [7–11] and matching-based approaches [4–6]. Neither of them is adequate to provide acceptable results in some particular cases. For the overlap-based approaches, two contours are considered connected if they overlap more than a predefined threshold. Due to the limitation of computing power in former days, the overlaps were evaluated by using 2D bounding boxes [8]. It may introduce unreasonable connections in some cases of concave contours, as shown in Fig. 3(a). Although Cj,1 and Cj+1,1 should not connect to each other, their 2D bounding boxes overlap significantly. Wang and Aggarwal [7] evaluate overlaps of the areas that are enclosed by the contours themselves. This approach cannot handle complex branching problems owing to the absence of an appropriate data structure that can be used to describe the relationship among contours. Another overlap-based method determines contour-correspondences by construct-

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 6 ( 2 0 0 7 ) 141–152

143

Fig. 1 – System flowchart.

ing the areas of difference [9–11]. The approach may fail to reconstruct the surface of a vessel with varying diameters. An area is considered degenerate if its corresponding area of difference is bounded by contours coming from one image only. As the example shown in Fig. 3(b), the two rings in slice i and i + 1 are considered degenerate, which makes the vessel cut into two separate pieces. The matching-based approaches solve the correspondence problem intuitively, but may yield ambiguous matching between the external and the internal contours. As shown in Fig. 4(a), Cj+1,1 and Cj+1,2 should connect

to Cj,1 and Cj,2 , respectively, but Cj+1,1 is most similar to Cj,2 . The external contour, Cj+1,1 , is wrongly connected to the internal contour, Cj,2 . A restriction is imposed on the types of connection [6] to avoid this problem. Contours can only be connected to others of the same type. The adoption of this restriction overcomes the ambiguous matching problems indeed, but new problems arise when contours must be connected to others of different types, as the example shown in Fig. 4(b). Ci+1,1 must be connected to both Ci,1 and Ci,2 . Fig. 4(c) shows the reasonable connection result. When the

144

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 6 ( 2 0 0 7 ) 141–152

Fig. 5 – Connection between contours of (a) the same type and (b) different types.

The first step of the algorithm is to compute the overlap ratio of contours in adjacent images. Instead of using 2D bounding boxes, we compute the ratio by using areas enclosed by the original contours to avoid the problem introduced by concave contours. After vertical projection, the overlap ratio  of two contours, Ci,k and Ci+1,l , is defined as the ratio of the overlapped area to the smaller area: (Ci,k , Ci+1,l ) = Fig. 2 – A tree structure that describes the relationship among contours in a slice image.

above restriction is applied, Ci+1,1 cannot be connected to Ci,2 but only to Ci,1 . The central hole is sealed after surface tiling as shown in Fig. 4(d), which is an unreasonable reconstruction result. An improved correspondence-determining algorithm is proposed to solve all the problems mentioned above.

area(Ci,k ) ∩ area(Ci+1,l ) min(area(Ci,k ), area(Ci+1,l ))

(1)

where area(·) is the function of the area enclosed by a contour, and can be implemented by the filling algorithms in [17]. The range of  is the interval from 0 to 1. Two contours of the same type tend to link to each other, if  of the two contours is higher. As the example shown in Fig. 5(a), Cj+1,1 and Cj+1,2 should be connected to Cj,1 and Cj,2 , respectively. However, contours of different types tend to link to each other while  of the two contours is lower. Fig. 5(b) illustrates the connection between contours of different types.

Fig. 3 – The particular cases encountered in the overlap-based correspondence-determining algorithms. The overlaps are evaluated (a) by using bounding boxes and (b) by constructing the areas of difference.

Fig. 4 – (a) and (b) are two particular cases encountered in the matching-based correspondence-determining algorithms. (c) Reasonable and (d) unreasonable reconstruction results of (b).

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 6 ( 2 0 0 7 ) 141–152

145

Although (Ci,2 , Ci+1,1 ) is close to zero, Ci+1,1 must be connected to Ci,2 besides Ci,1 . We define two thresholds, T1 and T2 , for determining the correspondences between contours. If Ci,k and Ci+1,j are of the same type: Correspondence (Ci,k , Ci+1,l )



=

connected,

T1 ≤ (Ci,k , Ci+1,l )

unconnected,

(Ci,k , Ci+1,l ) < T1

(2)

If Ci,k and Ci+1,l are of different types: Correspondence (Ci,k , Ci+1,l )



=

connected,

(Ci,k , Ci+1,l ) ≤ T2

unconnected, T2 < (Ci,k , Ci+1,l )

(3)

T1 and T2 are set based on our prior knowledge of the object to be reconstructed. A higher T1 implies a tendency to split the reconstructed object into separate pieces, and a lower T2 implies a tendency to seal holes in the reconstructed object. In fact, we may encounter numerous obstacles if we determine contour-correspondences by using Eq. (3) solely. Fig. 6 gives an example of the wrong connection between contours of different types. (Ci,3 , Ci+1,l ) is zero, but Ci,3 should not be connected to Ci+1,1 . For solving this problem, we consider the correspondence problem level-by-level with the tree structure in Section 3.1.1. Contour-correspondences at the child level are constrained by those at the parent level. Fig. 7(a) illustrates the restriction concisely. If we have known that Ci,1 is connected to Ci+1,2 , only Ci+1,2 , Ci+1,5 , Ci+1,6 , and Ci+1,7 can be the candidates for the corresponding contours of Ci,2 and Ci,3 . A practical example is shown in Fig. 7(b). Only Ci+1,1 is considered a candidate for the corresponding contour of Ci,3 , and only Ci+1,2 , Ci+1,3 , and Ci+1,4 are considered candidates for the corresponding contours of Ci,4 .

Fig. 7 – A restriction on acceptable connections is introduced to solve the problem in Fig. 6. (a) Illustration of the restriction. (b) A practical example.

gles will also make rendering and related applications difficult. Polygonal approximation of these contours is a viable alternative. A smaller (larger resp.) sampling step size implies a polygonal contour that is more (less resp.) similar to the original contour, but the polygon comprises more (fewer resp.) vertices and edges.

3.2.2. 3.2.

Surface tiling

3.2.1.

Polygonal contour approximation

The input contours for our system are extracted from the original section images. Each contour is assumed to be simple closed and is represented by a cyclic sequence of points. Every point is an eight-neighbor of its successive points in the cyclic sequence. The number of points on a contour varies with the image resolution. If the reconstruction procedure is applied to these contours directly, it will take much computation time. Furthermore, the considerable number of reconstructed trian-

Tiling triangles

Once two polygonal contours in adjacent images are determined to be connected, a triangle strip should be created by a tiling process to join the two contours. For the purpose of explanation, the lower of the two contours is denoted by C1 and the upper by C2 . C1 and C2 are defined by two sequences of distinct vertices, P1 , P2 , . . ., Pm and Q1 , Q2 , . . ., Qn , respectively. When the contours are viewed from the top, vertices of the external (internal resp.) contours are indexed in clockwise (counterclockwise resp.) order. The vertices are also indexed in a circular manner for convenience, i.e. Pi+m ≡ Pi and Qi+n ≡ Qi . The proposed tiling algorithm consists of two passes. In the first pass, the vertex pair, assumed to be (Pi , Qj ), of the shortest distance is selected as the initial span. Pi Qj is regarded as an edge of the first reconstructed triangle. The triangle can be obtained by including an adjacent vertex, Qj+1 or Pi+1 . If the criterion in Eq. (4) is satisfied, the strategy for advancing the vertex is based on the greedy algorithm that chooses the shorter span edge: min(D(Pi , Qj+1 ), D(Pi+1 , Qj )) < T3

Fig. 6 – Wrong connection between contours of different types.

(4)

where D(·, ·) is the function of the Euclidean distance between two vertices. T3 is a predefined threshold that depends on

146

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 6 ( 2 0 0 7 ) 141–152

the sampling step size of contours. Pi Qj Qj+1 is selected to be tiled if D(Pi , Qj+1 ) < D(Pi+1 , Qj ). Otherwise, Pi Qj Qi+1 is selected. The new generated edge, Pi Qj+1 or Pi+1 Qj , becomes the current span. The procedure continues recursively until the criterion in Eq. (4) is not satisfied. The terminal span, assumed to be Pk Ql , is defined as a forward partition marker (FPM). The next available span, say, Pk+d∗ Ql+e∗ is sought along the sequences of vertices. d* and e* are defined as follows: ∗

e = min(e|D(Pk+d , Ql+e ) < T3 , d ≥ 0, e ≥ 0), d∗ = min(d|D(Pk+d , Ql+e∗ ) < T3 , d ≥ 0)

(5)

Pk+d∗ Ql+e∗ is defined as a backward partition marker (BPM), and becomes the current span. The first pass is completed after the current span returns to the initial span. In the second pass, a FPM is picked arbitrarily as the starting partition marker. A BPM is sought starting from the FPM along C1 in the index-increasing direction, and the visited vertices are recorded in order. After a BPM is found, the next FPM is sought along C2 in the index-decreasing direction, and the visited vertices are also recorded in order. This procedure continues recursively until it goes back to the starting FPM. A 3D closed vertex chain can be obtained eventually. We apply triangulation to the 3D polygon by using heuristic methods. The objective function F is defined as follows: F=



[wA(i ) + (1 − w)L(i )]

(6)

i ∈ S

where A(i ) indicates the area of triangle i , and L(i ) indicates the perimeter of triangle i . S denotes the set of triangles obtained via the triangulation. The minimization of F is achieved by using the dynamic programming methods [18]. The influence of different objective functions is investigated thoroughly in [4]. The objective function, F, is adopted to reduces the probability of self-intersection and thin triangles. After the optimal triangulation is applied to the 3D polygon, the involved FPMs and BPMs are removed. An FPM is picked arbitrarily among the remainder as the new starting partition marker. The procedure continues recursively until no parti-

tion marker remains. The tiling process is completed after the second pass is finished. The distance measure function used in Eqs. (4) and  (5) is defined as D(P, Q) = D((x1 , y1 , z1 ), (x2 , y2 , z2 )) = 2

2

2

2

2

(x1 − x2 ) + (y1 − y2 ) + (z1 − z2 ) . Since C1 and C2 are parallel to the xy-plane, (z1 − z2 )2 in D(P, Q) should be constant for any P and Q. A projected Euclidean distance measure function is defined as Dp (P, Q) = Dp ((x1 , y1 , z1 ), (x2 , y2 , z2 )) =



(x1 − x2 ) + (y1 − y2 ) . Since the square root function is monotonic, D(Pi , Qj ) > D(Pm , Qn ) can imply Dp (Pi , Qj ) > Dp (Pm , Qn ), and vice versa. Consequently, Dp (P, Q) can be substituted for D(P, Q) in Eqs. (4) and (5) to accelerate the computation without affecting the tiling results.

3.2.3.

Degenerate cases and contour interpolation

The degenerate cases occur when a contour cannot find its corresponding counterparts in both of the adjacent images. Fig. 8(a) gives an example. Two external contours, Ci−1,1 and Ci,1 , are suggested to be connected to each other, but Ci,1 cannot find any corresponding contours in slice i + 1. One can only surmise that the region enclosed by Ci−1,1 and Ci,1 vanishes somewhere between slice i and slice i + 1. The simplest assumption is that the region vanishes just at slice i. As shown in Fig. 8(b), Ci,1 is tiled as a regular 3D polygon in Section 3.2.2. The tiling result is not smooth, and looks like being truncated unexpectedly. Another assumption is that the region shrinks gradually and vanishes at slice i + 1. Several intermediate contours are created by a shape-based interpolation algorithm to facilitate the tiling process. The shape-based interpolation algorithm is first introduced in [19] and then in [20,21] with further studies. For the interpolation of two contours, two signed distance maps are calculated with distance transformation first. The signed distance map records the distance from each pixel to its nearest contour pixel. The value is set to positive (negative resp.), if the pixel is inside (outside resp.) the contour. By linear interpolation of the two distance maps, the intermediate contours can be extracted from the interpolated distance maps with zero-crossing detection. A pixel in the image plane is denoted as x, and DMi (x) indicates the distance map of slice i. DMi (x) and DMi+1 (x) are both required for

Fig. 8 – Illustration of different assumptions about the object in the degenerate cases. (a) The input contours. (b) The object vanishes at slice i. (c) The object vanishes at slice i + 1. (d) An unexpected artificial spike will be produced when the area enclosed by Ci,1 is relatively small. (e) Linear prediction of the location of a reasonable vanishing point. (f) The reconstruction result of the proposed method.

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 6 ( 2 0 0 7 ) 141–152

147

the interpolation process. DMi (x) can be computed from Ci,1 , but DMi+1 (x) cannot be obtained directly. Since the region is assumed to shrink gradually and vanish at slice i + 1, DMi+1 (x) is defined as follows: DMi+1 (x) = DMi (x) − ı

(7)

where ı is defined as max(DMi (x)). Fig. 8(c) gives the tiling x

result under the assumption. The visualization result in Fig. 8(c) is obviously better than that in Fig. 8(b). If the area enclosed by Ci,1 is relatively small, unexpected artificial spikes will be produced as shown in Fig. 8(d). Consequently, it is inappropriate to assume that the region vanishes at slice i or slice i + 1 simply. If the surface of the object is smooth in the vicinity of slice i, Ci,1 can be assumed to be similar to Ci−1,1 . As illustrated in Fig. 8(e), the location of a reasonable vanishing point is obtained from the following equation: Ai−1 Ai

 D + D 2 v z

(8)

Dv

where Ai−1 and Ai denote the areas enclosed by Ci−1,1 and Ci,1 , respectively. Dz denotes the distance between slice i and slice i + 1. Dv denotes the distance from the vanishing point to slice i. Since the object shrinks gradually and vanishes at a point between slice i and slice i + 1, a constraint must be placed on Dv . Dv is derived from Eq. (8), and we have:

Dv =

⎧  A  ⎪ ⎨ i−1 ⎪ ⎩

Ai

−1 −1

Dz ,

if

  A

i−1

Ai



−1

>1

(9)

otherwise

Dz ,

Analogously to Eq. (7), the corresponding DMi+1 (x) for the region to vanish at the point of a distance Dv from slice i is defined as follows: DMi+1 (x) = DMi (x) −

Dz ı Dv

(10)

Substituting Eq. (9) into Eq. (10) yields:

⎧ ⎨ DMi+1 (x) =

⎩

  DMi (x) − DMi (x) − ı,

Ai−1 Ai



−1

  ı,

if

Ai−1 Ai



−1

>1

otherwise (11)

For the example in Fig. 8(d), the tiling result of the intermediate contours obtained according to Eq. (11) is shown in Fig. 8(f). Compared with Fig. 8(d) and (f) gives a better visualization result. A chamfering process is used in the original shape-based interpolation algorithm to provide an approximation of the Euclidean distance map. For better interpolation results, a fast linear time algorithm [22] is adopted to compute the exact Euclidean distance map. The degenerate cases happen at the top and bottom of the object or terminals of branches only. In most situations, the surface refinement procedure does not take much additional computation time.

Fig. 9 – Illustration of the tiling process for branching objects.

The tiling of two complicated contours in adjacent slices cannot always provide valid triangulation. If the two contours differ much from each other in shape, self-intersecting surfaces may be introduced during the tiling process [11]. The addition of intermediate contours is suggested to alleviate the problem. The shape-based interpolation algorithm mentioned above can handle the problem. With the information of contours in both slices, the respective distance maps can be computed more directly than those in degenerate cases.

3.2.4.

Branching problems

The proposed tiling algorithm can handle the branching problem without any modification. For explanation, an example of simple branching is given in Fig. 9. After the contourcorrespondence analysis of the three external contours, C1,1 is assigned to connect to C2,1 and C2,2 . For C1,1 and C2,1 , the first pass of the tiling process tiles the similar contour segments and finds two partition markers, FPM1 and BPM1 . For C1,1 and C2,2 , the same procedure can find four partition markers, FPM2 , BPM2 , FPM3 , and BPM3 . In the second pass, two 3D closed vertex chains are obtained from the six partition markers. The first vertex chain starts from FPM1 to BPM3 via C1,1 in the index-increasing direction, and goes to FPM3 via C2,2 in the index-decreasing direction, and goes to BPM1 via C1,1 in the index-increasing direction, and goes back to FPM1 via C2,1 in the index-decreasing direction. The second vertex chain is obtained from the remaining partition markers, FPM2 and BPM2 . The two vertex chains are tiled by minimizing F in Eq. (6) to complete the whole tiling process. Unlike in the case of one-to-one connections, the vertex chain, such as the first vertex chain in the example, may comprise contour segments from more than two contours. Nevertheless, whether the object branches or not makes no difference to the proposed tiling algorithm. The tiling procedure for multibranching is analogous with that for simple branching, and is omitted here. Usually, an extra contour interpolation algorithm is applied before the tiling process to create intermediate contours. By tiling these additional contours, a better visualization result can be achieved. The improvement is illustrated by an example in Section 4.

3.3.

Rendering and quantitative analysis

3.3.1.

Three-dimensional rendering

The rendering of polyhedral models is composed of several procedures, including perspective projections, visible-surface determination, illumination, and shading [23]. If the process

148

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 6 ( 2 0 0 7 ) 141–152

is performed by pure software methods, real-time rendering of high resolution polyhedral models cannot be achieved. A widely used graphics API standard, OpenGLTM (Open Graphics Library, a trademark of Silicon Graphics Inc.), is adopted to facilitate the rendering process. OpenGLTM provides an interface between hardware and software, which enables the user to utilize the graphics hardware to perform the rendering efficiently. In some medical applications, polyhedral models and medical volume data set must be rendered simultaneously. The volume data can only be rendered by software methods [24–26], because OpenGLTM does not support the rendering of volume data. Meyer et al. [27] propose a pipelined hybrid rendering algorithm to render both polyhedral models and volume data simultaneously.

3.3.2.

Surface area and volume computation

Besides 3D rendering, quantitative analysis of the object is useful in practical applications. Since the reconstructed object is represented by a polyhedral model, the surface area and volume of the object can be approximated with the polyhedral model. The surface area of the polyhedral model is the sum of the areas of all triangular patches. The origin and the three vertices of each triangular patch can form a tetrahedron. For example, PQR and the origin O can form the tetrahedron PQRO. Its signed volume, V(PQR), is obtained by computing − → −→ − → the mixed product of the three position vectors, [OP OQ OR], and we have: V(PQR) =

1 − 1 − → −→ − → → −→ − → [OP OQ OR] = (OP × OQ) · OR 6 6

(12)

− → denote the normal vector of PQR, and is defined as (PQ × Let n − → : PR). The sign of V(PQR) depends on the direction vector of n − → · OP) sgn(V(PQR)) = sgn(n

(13)

Since each triangular patch obtained from the proposed tiling algorithm has an outward normal vector, the volume of the polyhedral model is the sum of signed volume of these tetrahedra.

4. Implementation and experimental results The performance of the proposed system is measured on a desktop computer with the AMD AthlonTM 64 X2 DualCore 3600+ CPU and 1 GB of main memory. We use Microsoft Visual C++ 6.0 to implement the proposed system under the operating system of Microsoft Windows XP Professional. Four parameters in the system need to be set: T1 , T2 , T3 , and w. T1 and T2 are set based on our prior knowledge of the object to be reconstructed. A higher T1 implies a tendency to split the reconstructed object into separate pieces, and a lower T2 implies a tendency to seal holes in the reconstructed object. For the examples in the implementation section, we set T1 to 0.3 and T2 to 0.7. T3 decides the percentage of triangles that are obtained from the greedy algorithm in the hybrid tiling algorithm. Increasing T3 will raise the percentage, and accelerate the tiling process. When the contours to be tiled are very different in shape, size, or orientation, unacceptable tiling results or self-intersecting surfaces may be introduced during the greedy tiling process. Consequently, there is a trade-off between efficiency and accuracy. When the distance functions D(·, ·) in Eqs. (4) and (5) are replaced by the projected distance function Dp (·, ·), T3 is set to



2

64 + (sampling step size of contours) empirically. We also empirically set w to about 0.8, which provides the best tiling result for most situations. Fig. 10 shows a synthetic example of the degenerate case. The input contours are shown in Fig. 10(a) and (b). Two external contours, Ci−1,1 and Ci,1 , are suggested to connect to each other, but Ci,1 cannot find any corresponding contours in slice i + 1. If the object is assumed to vanish just at slice i, the reconstruction result is shown in Fig. 10(c). The reconstructed polyhedral model looks like being truncated unexpectedly and therefore does not provide a smooth visualization result. If the object is assumed to shrink gradually and vanish at slice i + 1, the reconstruction result is shown in Fig. 10(d). It also produces an unexpected artificial spike. If the vanishing point is obtained from Eq. (9), the reconstruction result is shown in Fig. 10(e). Obviously, it gives the best visualization result.

Fig. 10 – A synthetic example of the degenerate case. (a) and (b) are the input contours. (c) and (d) are the reconstruction results under the assumptions that the object vanishes at slice i and slice i + 1, respectively. (e) The reconstruction result of the proposed method.

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 6 ( 2 0 0 7 ) 141–152

149

Fig. 11 – A synthetic branching example. (a) and (b) are the input contours. (c) Tiling the branch without intermediate contours. (d) Intermediate contours after contour interpolation. (e) Tiling the branch with the intermediate contours.

Some methods reconstruct smooth surfaces by fitting parametric surfaces to Ci−1,1 and Ci,1 , but they are complicated and very time consuming. Fig. 11 shows a synthetic branching example. The input contours are shown in Fig. 11(a) and (b). If we tile the branch without intermediate contours, the reconstruction result is shown in Fig. 11(c). The reconstruction result can be refined by means of contour interpolation. The sequence of intermediate contours is shown in Fig. 11(d) and (e) gives the result of the refinement. The branching point is moved to a position between slice i and slice i + 1, which gives a better visualization result. A practical example is presented in Fig. 12. The original dataset consists of 36 confocal microscopic images. The object that we try to reconstruct is the mushroom body of a cock-

roach. Contours of some slices are shown in Fig. 12(a). The parts enclosed by the dashed line are similar to the example given in Fig. 4. For a reasonable reconstruction result, we must keep the hole instead of sealing it. The reconstructed polyhedral model is shown in Fig. 12(b). It contains 4978 vertices and consists of 9956 triangular patches. As we can see, the calyxes of the mushroom body are well reconstructed, and the holes are not sealed. Fig. 13(a) shows the reconstruction result of another real data set obtained from confocal microscopy. The reconstructed Drosophila cortex also gives an accurate and good visualization result. A synthetic dataset is used to evaluate the effect of polyhedral representation in quantitative analysis. The synthetic object is defined by the inequalities: x2 + y2 ≤ ((z2 /100) + 10)2 and 0 ≤ z ≤ 100. After sampling, the polyhedral model is reconstructed and shown

Fig. 12 – (a) Contours of some slices and (b) the reconstructed polyhedral model of the mushroom body of a cockroach.

150

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 6 ( 2 0 0 7 ) 141–152

Fig. 13 – Reconstruction results of (a) real data and (b) synthetic data.

to branches usually enclose larger areas. The correlation coefficient between the computation time for and the number of branches is 0.71, which indicates a moderately strong positive relationship. The computation time for degenerate cases is less than 5% of the total computation time. The computation time for branches is about 20–45% of the total computation time. In real-time applications, the interpolation procedure for branches can be omitted to accelerate the rendering. Even if the whole procedures are performed, most of the datasets still can be reconstructed fast within 1 s. Table 2 shows a comparison of the computation time of the proposed tiling method and the typical method of Barequet and Sharir [5]. Since the proposed approach does not require a contour-matching procedure, it is faster than the typical method. The proposed method can save about 40% of the computation time for tiling. Two practical applications are based on the proposed system. The screen shots are shown in Fig. 14. The first application is for the 3D visualization of CT scans. It reconstructs several objects and displays them together in the same frame. As shown in Fig. 14(a), it is the reconstruction result of the human brain. The boundary of the brain is displayed in gray. The region of internal bleeding is displayed in red, and the region of ventricle dilation in yellow. The second application is for the quantitative analysis of confocal microscopic scans. As shown in Fig. 14(b), it is the reconstruction result

in Fig. 13(b). It comprises 8118 vertices and 16,232 triangular patches. The actual surface area and volume of the object are 8.572 × 104 unit2 and 8.692 × 105 unit3 , respectively. The surface area and volume of the polyhedral model are 8.555 × 104 unit2 and 8.663 × 105 unit3 , respectively. The percentage errors in surface area and volume are 0.19% and 0.33%, respectively. They are not significant for most practical applications. Table 1 summarizes the performance of the proposed system. Ten datasets of the confocal microscopic scans are tested. Each dataset is obtained from different Drosophila brains, and the object that we try to reconstruct is the mushroom body of the Drosophila. The image size of the datasets is 240 × 480. The correlation coefficient between the total computation time and the number of vertices is 0.91, which indicates a fairly strong positive relationship. Consistent with intuition, the total computation time highly depends on the size and the resolution of the reconstructed polyhedral model. The correlation coefficient between the computation time for and the number of degenerate cases is 0.25, which indicates that no significant relationship exists. This is because contours related to degenerate cases usually enclose small areas. The computation time for degenerate cases depends more on the size of the enclosed areas than on the number of degenerate cases. As compared with degenerate cases, contours related

Table 1 – Performance of the proposed system Data set No. 1

No. 2

No. 3

No. 4

No. 5

No. 6

No. 7

No. 8

No. 9

Number of slices Number of degenerate cases (I) Number of branches (II) Total computation time (ms)

109 8 4 690

107 10 6 595

117 10 6 624

114 10 6 696

104 10 6 710

108 12 8 1,004

113 10 6 858

113 10 6 895

112 11 7 722

110 10 6 622

Computation time for (I) (ms) (%) Computation time for (II) (ms) (%)

23 3.33 182 26.4

26 4.37 179 30.1

27 4.33 176 28.2

23 3.30 220 31.6

25 3.52 232 32.7

28 3.26 299 34.8

42 4.69 244 27.3

36 4.99 233 32.3

29 4.66 202 32.5

Number of vertices Number of patches

9,045 18,082

7,951 15,894

8,253 16,498

8,935 17,862

8,784 17,560

26 2.59 431 42.9 10,045 20,082

10,061 20,114

10,093 20,178

8,798 17,588

No. 10

7,436 14,864

151

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 6 ( 2 0 0 7 ) 141–152

Table 2 – Computation time of the proposed tiling method and the typical method of Barequet and Sharir [5] Data set

The proposed method (I) (ms) Barequet and Sharir (II) (ms) (I)/(II)

No. 1

No. 2

No. 3

No. 4

No. 5

No. 6

No. 7

No. 8

No. 9

No. 10

235 422

187 328

188 344

249 437

282 485

344 625

359 609

437 672

250 422

173 313

0.56

0.57

0.55

0.57

0.58

0.55

0.59

0.65

0.59

0.55

Fig. 14 – Two practical applications based on the proposed system. (a) A system for 3D visualization of CT scans. (b) A system for the quantitative analysis of confocal microscopic scans. (For interpretation of the references to color in the text, the reader is referred to the web version of the article.)

of the mushroom body of the Drosophila. The surface area of the mushroom body is 4.18 × 104 ␮m2 , and the volume is 2.69 × 105 ␮m3 .

5.

Conclusions

Different reconstruction results can be obtained from a single data set of parallel planar contours. It is impossible to reconstruct the object perfectly in the case of undersampling. For solving this ill-posed problem, well-posedness must be restored by imposing restrictions on the feasible solutions. For example, the reconstruction result is usually expected to be smooth and consistent with our prior knowledge. In this paper, we develop an automatic surface reconstruction system that provides a solution to the ill-posed problem. It brings improvements in accuracy, computation time and smoothness. The proposed correspondence-determining algorithm provides more reasonable contour-correspondences than the existing algorithms for some particular patterns of contour pairs, thus increasing the accuracy of reconstruction. The proposed tiling algorithm tiles the corresponding contours without the requirement of a contour-matching procedure, thus saving more computation time than ordinary hybrid approaches. The adoption of contour interpolation yields smoother reconstruction results in degenerate and branching cases. The system also calculates the surface area and volume of the object to facilitate practical applications. Since the polyhedral model is reconstructed from polygonal contours of a fixed step size, triangular patches distribute evenly on the surface of the polyhedral model. In our future works, a polygon reduction algorithm will be applied to reduce

the number of triangular patches on smooth areas. That will improve the efficiency of the follow-up applications. We also plan to enhance the proposed correspondence-determining algorithm. If the original scanned images are available, the texture information of the object can be utilized for contourcorrespondence analysis.

Acknowledgement This work was supported by the Brain Research Center of the University System of Taiwan and the National Center for Highperformance Computing.

references

[1] G.M. Nielson, On marching cubes, IEEE Trans. Vis. Comput. Graph. 9 (3) (2003) 283–292. [2] D. Meyers, S. Skinner, K. Sloan, Surfaces from contours, ACM Trans. Graph. 11 (3) (1992) 228–258. [3] D. Weinstein, Scanline surfacing: building separating surfaces from planar contours, in: Proceedings of the 11th IEEE Visualization 2000 Conference, 2000, pp. 283–289, 568. [4] G. Barequet, D. Shapiro, A. Tal, History consideration in reconstructing polyhedral surfaces from parallel slices, in: Proceedings of the IEEE Visualization ’96 Conference, 1996, pp. 149–156, 479. [5] G. Barequet, M. Sharir, Piecewise-linear interpolation between polygonal slices, Comput. Vis. Image Und. 63 (2) (1996) 251–272. [6] M.Y. Shih, D.C. Tseng, Versatile surface model reconstruction from serial planar contours, in: Proceedings

152

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15]

[16]

[17]

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 6 ( 2 0 0 7 ) 141–152

of the 22nd Annual EMBS International Conference, vol. 4, 2000, pp. 2910–2913. Y.F. Wang, J.K. Aggarwal, Surface reconstruction and representation of 3D scenes, Pattern Recogn. 19 (3) (1986) 197–207. A. Ekoule, F. Peyrin, C. Odet, A triangulation algorithm from arbitrary shaped multiple planar contours, ACM Trans. Graph. 10 (2) (1991) 182–199. J.F. Guo, Y.L. Cai, Y.P. Wang, Morphology-based interpolation for 3D medical image reconstruction, Comput. Med. Imaging Graph. 19 (3) (1995) 267–279. R. Klein, A. Schilling, W. Straßer, Reconstruction and simplification of surfaces from contours, Graph. Mod. 62 (6) (2000) 429–443. J.-M. Oliva, M. Perrin, S. Coquillart, 3D reconstruction of complex polyhedral shapes from contours using a simplified generalized Vorono¨ı diagram, Comput. Graph. Forum 15 (3) (1996) C397–C408. C. Gitlin, J. O’Rourke, V. Subramanian, On reconstructing polyhedra from parallel slices, Int. J. Comput. Geom. Appl. 6 (1) (1996) 103–122. Y. Shinagawa, T.L. Kunii, The homotopy model: a generalized model for smooth surface generation from cross sectional data, Visual Comput. 7 (2–3) (1991) 72–86. E. Welzl, B. Wolfers, Surface reconstruction between simple polygons via angle criteria, J. Symb. Comput. 17 (4) (1994) 351–369. J.D. Boissonnat, Shape reconstruction from planar cross-sections, Comput. Vis. Graph. Image Process. 44 (1) (1988) 1–29. J.D. Boissonnat, B. Geiger, Three-dimensional reconstruction of complex shapes based on the Delaunay triangulation, INRIA, Research Report no. 1697, 1992. J.D. Foley, A. van Dam, S.K. Feiner, J.F. Hughes, Advanced geometric and raster algorithms: filling algorithms, in:

[18]

[19]

[20]

[21]

[22]

[23]

[24] [25]

[26]

[27]

Computer Graphics: Principles and Practice, 2nd ed., Addison-Wesley, Massachusetts, 1996, pp. 979–986. T.H. Cormen, C.E. Leiserson, R.L. Rivest, C. Stein, Dynamic programming, in: Introduction to Algorithms, 2nd ed., MIT Press, Massachusetts, 2001, pp. 323–369. S.P. Raya, J.K. Udupa, Shape-based interpolation of multidimensional objects, IEEE Trans. Med. Imaging 9 (1) (1990) 32–42. G.T. Herman, J. Zheng, C.A. Bucholtz, Shape-based interpolation, IEEE Comput. Graph. Appl. 12 (3) (1992) 69–79. D. Cohen-Or, D. Levin, A. Solomovici, Three-dimensional distance field metamorphosis, ACM Trans. Graph. 17 (2) (1998) 116–141. C.R. Maurer Jr., R. Qi, V. Raghavan, A linear time algorithm for computing exact Euclidean distance transforms of binary images in arbitrary dimensions, IEEE Trans. Pattern Anal. Mach. Intell. 25 (2) (2003) 265–270. J.D. Foley, A. van Dam, S.K. Feiner, J.F. Hughes, The quest for visual realism: rendering techniques for shaded images, in: Computer Graphics: Principles and Practice, 2nd ed., Addison-Wesley, Reading, Massachusetts, 1996, pp. 612–615. M. Levoy, Efficient ray tracing of volume data, ACM Trans. Graph. 9 (3) (1990) 245–261. S. Parker, M. Parker, Y. Livnat, P.P. Sloan, C. Hansen, P. Shirley, Interactive ray tracing for volume visualization, IEEE Trans. Vis. Comput. Graph. 5 (3) (1999) 238–250. ¨ H. Hauser, L. Mroz, G.I. Bischi, M.E. Groller, Two-level volume rendering, IEEE Trans. Vis. Comput. Graph. 7 (3) (2001) 242–252. ¨ J. Meyer, S. Gelder, K. Kretschmer, K. Silkenbaumer, H. Hagen, Interactive visualization of hybrid medical data sets, in: Proceedings of the Fifth International Conference in Central Europe on Computer Graphics and Visualization (WSCG’97), vol. 2, 1997, pp. 371–380.