Efficient computation of enclosed volume and surface area from the same triangulated surface representation

Computerized Medical Imaging and Graphics 35 (2011) 460–471

Contents lists available at ScienceDirect

Computerized Medical Imaging and Graphics journal homepage: www.elsevier.com/locate/compmedimag

Efficient computation of enclosed volume and surface area from the same triangulated surface representation Ingela Nyström a , George J. Grevera b,∗ , Bruce E. Hirsch c , Jayaram K. Udupa d a

Centre for Image Analysis, Uppsala University, Uppsala, Sweden Computer Science Department, St. Joseph’s University, Philadelphia, PA, USA Dept. of Neurobiology and Anatomy, College of Medicine, Drexel University, Philadelphia, PA, USA d Medical Image Processing Group, Dept. of Radiology, University of Pennsylvania, Philadelphia, PA, USA b c

a r t i c l e

i n f o

Article history: Received 9 June 2008 Received in revised form 20 October 2010 Accepted 8 November 2010 Keywords: Digital boundary Surface tracking Shell structure Marching Cubes Isosurface Shape parameters Boundary measurement Volume measurement Rendering

a b s t r a c t We demonstrate that the volume enclosed by triangulated surfaces can be computed efficiently in the same elegant way the volume enclosed by digital surfaces can be computed by digital surface integration. Although digital surfaces are effective and efficient for visualization and volume measurement, their drawback is that surface area measurements derived from them are inaccurate. On the other hand, triangulated surfaces give more accurate surface area measurements, but volume measurements and visualization are less efficient. Our data structure (called t-shell) for representing triangulated digital surfaces retains advantages and overcomes difficulties of both the digital and the triangulated surfaces. We create a lookup table with area and volume contributions for each of the 256 Marching Cubes configurations. When scanning the shell (e.g., while creating it), the surface area and volume are incrementally computed by using the lookup table and the current x co-ordinate, where the sign of the x component of the triangle normal indicates the sign of the volume contribution. We have computed surface area and volume for digitized mathematical phantoms, physical phantoms, and real objects. The experiments show that triangulated surface area is more accurate, triangulated volume follows digital volume closely, and that the values get closer to the true value with decreasing voxel size. © 2010 Elsevier Ltd. All rights reserved.

1. Introduction Currently, a number of imaging devices that generate volumetric (3D) images are available and are extensively used especially in medical imaging. There is usually an object of interest for which such an image is generated, for example, an organ, a tissue component, a tumor, a physical phantom, or a mathematical phantom. The 3D imaging operations typically performed on these images may be classified into four groups: preprocessing, visualization, manipulation, and analysis [1]. The preprocessing operations aim at improving the definition of and/or extracting the object of interest. The purpose of visualization is to assist humans in perceiving and comprehending the 3D object and its full form and shape. The manipulation operations allow humans to interactively change the objects, mimicking surgery. Finally, the analysis operations are intended to quantify the object, e.g., geometrically, morphologically, architecturally, and functionally. In this paper, we touch upon

∗ Corresponding author at: Computer Science Department, St. Joseph’s University, 5600 City Avenue, Philadelphia, PA 19131, USA. Tel.: +1 610 660 1535. E-mail addresses: [email protected] (I. Nyström), [email protected] (G.J. Grevera), [email protected] (B.E. Hirsch), [email protected] (J.K. Udupa). 0895-6111/$ – see front matter © 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.compmedimag.2010.11.001

some aspects of visualization, but concentrate mainly on analysis, particularly on estimating surface area and volume. A simple way to compute object volume is to count all voxels of the object while using a connected component labelling algorithm [2]. However, boundary tracking methods are more economical since they visit only elements in the vicinity of the boundary without having to visit the interior region. In fact, such methods can also label connected components [3]. Therefore, it makes sense to seek ways of estimating volume from boundary information only. Situations exist, such as when manipulating objects, where we are given only the boundary, not the scene data, and we need quantitative measurements also of modified objects which are also available only as boundary data. For representing objects by their boundary surfaces, the two classical approaches are digital surface tracking, as in the approach to boundary detection [4], and isosurface construction by polygonal surfaces, commonly by using triangles as in the Marching Cubes (abbreviated from now on by MC) algorithm [5,6]. Actually, there is a close relationship between the two approaches, since a digital surface can be transformed directly into a triangulated isosurface [7,8]. Digital boundaries and surfaces have many elegant geometrical and topological properties that lead to their efficient construction

I. Nyström et al. / Computerized Medical Imaging and Graphics 35 (2011) 460–471

(from image data) [9]; ultra fast rendering [10,11]; manipulation [10]; measurement [4]; and to elegant generalizations of the concepts to higher dimensional spaces [12]. The MC algorithm [5,6] was an attempt to bridge the gap between digital surface methods [13–15] (which predate the MC family of algorithms) and the traditional polygonal approaches established in computer graphics. The spirit of this bridging has been identified in [16], wherein it was recognized that the surfaces produced by MC algorithms have a digital embedding, and therefore, can be cast within the same efficient framework, namely that of a shell [9,17], which is used for efficiently representing, rendering, manipulating, and measuring digital surfaces. This recognition led to the idea of the t-shell [16] (short for triangulated shell) for dealing with triangulated surfaces. In spite of their many desirable properties, one drawback of digital surfaces has been the lack of a simple method to accurately estimate their area. The identity between the t-shell and the shell provides a natural means to overcome this problem. The volume enclosed by a digital surface, on the other hand, can be computed trivially by digital surface integration [4] while tracking the surface in a gray or binary image (or after its creation) without having to visit (and count) any voxel in its interior. Conversely, while area computation of triangulated surfaces is trivial, the computation of the volume enclosed by them is challenging. The purpose of this paper is to provide an effective solution to this dilemma. The basis of this solution is the recognition of the fact that the triangulated surfaces produced by the MC family of methods are indeed digital surfaces. Therefore, the ideas underlying digital surface integration methods can be adopted to the case of MC surfaces. Thus we will have efficient and accurate surface area and volume estimation from the same triangulated surface representation. From a survey of the literature, we observed that, because of the difficulty in computing the volume enclosed by triangulated surfaces, surface area is computed mostly from triangulated surfaces, whereas the volume is estimated by counting voxels – in effect, from digital surfaces. In situations where both measures are needed for the same physical object, this dichotomy amounts to using one geometric object (triangulated surface) for surface area estimation and a different geometric object (digital surface) for volume estimation, which implies that the two measures correspond to two different, albeit similar, objects. The purpose of this paper is to solve this “different objects” problem by representing a simple and efficient method for computing both surface area and enclosed volume from the same MC representation of an object. We note also that there exist a variety of approaches in the literature for estimating surface area more accurately than directly via digital surfaces. Some of these have attempted to estimate local surface area by examining the local structure of digital surfaces [18–20], and others have attempted to create better triangular approximations of the digital surface [21,22]. Their improved accuracy has been demonstrated mostly on mathematical and convex objects. This is mainly due to the difficulty in physically measuring the true surface area of real objects of complex shape. In some of the digital approaches [18–20], the issue of “different objects” mentioned above arises, and it is not clear as to how to estimate volume in those approaches so that the surface area and volume correspond to the same geometric object. Although a careful investigation of the real improvement in surface area accuracy achieved by these more sophisticated methods over the straightforward MC method on real objects of complex shape is certainly worthwhile, such a study does not exist. Perhaps owing to this gap, and owing to the simplicity of the MC method, the most commonly used surface area and volume estimation methods are via the MC triangulated surface and voxel counting, respectively. Our focus in this paper is not on the most accurate methods of estimating these measures, but on most efficiently estimating these entities by using the most commonly

461

employed methods and by addressing the “different objects” problem. The main purpose of this paper is to demonstrate that, via the same single boundary representation scheme, namely t-shell, we can combine the strengths and overcome the shortcomings of both digital and triangulated surfaces and compute both surface area and volume trivially and efficiently, and still retain the digital setting. We describe the basic concepts of the two methods of surface representation and of t-shells in Section 2, and the new method of area and volume estimation in Section 3. Our experiments and results comparing the digital and triangulated methods are presented in Section 4, and Section 5 states our conclusions. A preliminary version of this paper was presented in [23]. 2. Object representation by boundary We refer to any gray-level volume (3D) image as a scene and represent it by a pair C = (C, f ), where C, called the scene domain, is a rectangular 3D array of cuboidal volume elements, called voxels, and f is a function that assigns a value f (v) to every voxel v ∈ C called scene intensity. The range of f is usually a set of integers. When this range is {0, 1}, we call C a binary scene. Scenes contain information about certain objects of interest, which are denned and delineated in the scene by a scene segmentation operation. The result sought very frequently is either a binary scene or a surface representation of the object. For what we wish to accomplish in this paper, it is the latter form that is relevant. In this section, we shall outline several common methods of representing objects by surfaces, wherein object information is derived from scenes. 2.1. Digital surface representations We will consider two methods of digital surface representation: (i) by a set of voxels [10], and (ii) by a set of oriented faces of the voxels [12,13]. In both cases, the boundary elements are all of identical size and shape, which makes efficient storage of the surface by using clever data structures (called shells) possible. Efficiency comes from the fact that the elements need not be stored explicitly in terms of their co-ordinates since they are all identical in size and shape, and have only a small number of distinct orientations. These surfaces satisfy properties such as closure, orientedness, and connectedness. In case (i), the digital surface is represented by a set of those voxels in the object that have a 6-adjacent neighbor (i.e., a face neighbor) in the background. The information associated with each boundary voxel consists of the x co-ordinate of the voxel (y and z coordinates are given implicitly by the data structure) and a neighbor code describing the configuration of its 6-neighborhood. Hence, it is known not only that a voxel is facing the background, but also whether it is facing the background in positive or negative (x, y, z) directions. This is useful information when computing the volume enclosed by the surface, as will be described in Section 3.2. In case (ii), the surface is represented by the set of oriented faces of voxels. The number of boundary elements in an object becomes larger in this case than in case (i). The faces are assigned an orientation, which means that a face with its normal vector pointing in the +x direction is distinguished from a face at the same location with the normal vector pointing in the −x direction. The shell data structure consists of a one-dimensional pointer array PA and a list LS of shell elements that comprise the object surface. PA contains as many elements as there are rows in the 3D scene that contains the object surface. Each pointer in PA points to the location in LS that corresponds to the first shell element in LS corresponding to that row. Each shell element is stored in LS as a 3-tuple xs , ncs , nvs , where xs is the x co-ordinate of the shell element s (voxel or voxel

462

I. Nyström et al. / Computerized Medical Imaging and Graphics 35 (2011) 460–471

Fig. 1. Three m-cubes of 2 × 2 × 2 voxels, where the voxels (actually the voxel centers) marked by 䊉 are inside the object (value 1) and the other voxels are outside (value 0). The surfaces for these configurations consist of one, three, and four triangles, respectively. (Fourth voxel marked 䊉 is not visible in (c).)

face), ncs is the neighbor code of s, and nvs is the surface normal vector associated with s. For the purpose of this paper, the first two elements in the 3-tuple are adequate. The third element is utilized for rendering the surface. See [10] for more details. 2.2. Polygonal surface representation A commonly used approach to polygonal surface representation is by a set of triangles, as in the MC algorithm [5,6]. We use the term m-cube (to indicate an individual Marching Cube, sometimes also called a cell in the literature) in a given scene domain to denote the cube bounded by the centers of the eight voxels in a 2 × 2 × 2 neighborhood of voxels in the scene domain. Each corner of the mcube corresponds to a voxel. The 6-, 18-, and 26-adjacency between two voxels can be described in terms of the m-cube as follows: connected by an edge of the m-cube, having a face of the m-cube in common, and having the whole m-cube in common, respectively. The possible number of configurations of the eight voxels forming an m-cube in a binary scene is 28 = 256. In the MC algorithm [6], each of these configurations yields zero to four triangles constituting the part of the surface passing through the m-cube corresponding to the configuration. See Fig. 1 for examples. The configurations are commonly grouped into symmetry and complementary cases, resulting in (14 or) 15 distinct cases out of the 256 possible cases [6,24]. There are two configurations with no triangles, which arise when the m-cube is situated completely inside or completely outside the object. For a given m-cube, the exact location of a vertex of a triangle along an m-cube edge is determined in one of two ways. Either it is considered to be the midpoint of the edge, or its location on the edge is determined by linear interpolation of the original scene intensities f(c) and f(d) at the voxels c and d constituting the endpoints of the m-cube edge. The midpoint solution is used in the rest of the paper for simplicity; we will come back to the other case in Section 5. The correct connection among intersection points for forming triangles is tricky for some configurations of voxels in the m-cube. This problem in the original MC algorithm was pointed out by Dürst

in 1988 [25]. If the same surfaces are used for such m-cube configurations and for their complementary cases (i.e., when object and background are inverted), closed surfaces are not produced. The problem requires careful consideration since any object with finite support must be enclosed by a closed surface to make topological sense. Further, the act of measuring the volume enclosed by a surface is meaningful only when the surface is closed. This problem has been examined by numerous investigators in the literature [24,26,27]. We have taken the approach of [24] where the problem is avoided by using more complex surfaces for these cases, as exemplified in Fig. 2.

2.3. t-Shell representation In previous shell data structures used for representing digital surfaces, the shell elements have been voxels or oriented faces of voxels [10,12], as briefly described in Section 2.1. In the t-shell [16], the shell elements are m-cubes, and the t-shell structure is very similar to the shell structure. The t-shell contains only the m-cubes on the boundary of the object, which are stored rowwise. This means that m-cubes completely inside or completely outside the object are not stored in the t-shell. The element xs corresponds to the x co-ordinate of the m-cubes’s front-upper-left voxel (other co-ordinates can be derived from this information), and ncs corresponds to the m-cube configuration code. In short, the data structure of t-shell is identical to that of shell; only the interpretation of the meaning of shell elements is different.

3. Surface area and enclosed volume We use notations A(u) and V(u) to denote area and volume functions, respectively, where u may be the original object O, the original surface S, a digital surface Sd , a triangulated surface St , or a triangle t. Sd and St represent different approximations of S. We are interested in the following entities for any object denned in a given scene:

Fig. 2. (a) A tricky m-cube configuration. (b) The complement of (a), which would result in surfaces that are not closed. (c) A set of five triangles that will produce a closed surface for this example. (Fifth voxel marked 䊉 is not visible.)

I. Nyström et al. / Computerized Medical Imaging and Graphics 35 (2011) 460–471

• true surface area, A(S) (when known), • area of the surface of the object approximated by a digital surface, A(Sd ), • area of the surface of the object approximated by a triangulated surface, A(St ), • true enclosed volume, V(S) (when known), • volume enclosed by a digital surface approximation of the object surface, V(Sd ), • volume enclosed by a triangulated surface approximation of the object surface, V(St ). We will describe the underlying idea in the 2D case first, where the perimeter of and the area enclosed by a given boundary, which is represented by MC triangulation, are computed in an incremental fashion by scanning the boundary of the object represented by its t-shell. Understanding the 2D case helps when we later extend the concepts to the 3D case. The intuitive idea extends easily, but as is usually the case when going from two dimensions to three dimensions, it is not a pure generalization; each m-cube holds contributions from one or more triangles, unlike in the 2D case, where each configuration consists of only one contribution. We shall discuss surface area first, and then the enclosed volume. 3.1. Surface area Surface area estimation has been studied earlier [28,29], but not with the approach to simultaneously compute the volume enclosed by the same boundary representation. Our primary aim is to use the same t-shell boundary representation giving accurate results for both measures simultaneously and with minimal computation.

√ and weight 2 for diagonal steps are not optimal when measuring lengths of line segments [30]. In the digital case, another estimate of the perimeter is to count the number of pixel edges between the object and the background. This results in an overestimate. This can be improved by finding a better approximation to the boundary. One such improvement, via the MC approximation, is illustrated in Fig. 3. An m-square is the 2D analog of the m-cube, i.e., a square bounded by the centers of the four pixels in a 2 × 2 neighborhood. In this representation, the contribution to the boundary length may be precomputed and stored in a lookup table for the 24 different MC configurations. This idea will now be extended to the 3D case. 3.1.2. Extending to three dimensions In the purely digital case, where the digital surface Sd may be represented by a digital shell SD consisting of the faces of the voxels at the boundary, the number of faces gives an estimate of the surface area, A(Sd ). Similar to the 2D case, this is an overestimate. If the surface is represented by a set St of triangles, a better surface area estimate, A(St ), is obtained efficiently by table lookup, as described below. Let ST be a t-shell that represents a 3D boundary surface of an object in a given scene. Let c denote an m-cube (a t-shell element) of ST and let (c) be the set of triangles in the MC configuration ncc of c. Our idea is to create a lookup table that stores the contributions of the different MC configurations to surface area by summing the areas of the included triangles, and then to use this table to compute the area of the surface represented by ST . The surface area estimate, A(St ), is incrementally computed for each m-cube of the object: A(St ) =



A(t) =

C ∈ ST t ∈ (c)

3.1.1. Underlying idea in two dimensions In the 2D case, the equivalent of surface area is the perimeter of an object, i.e., the boundary length. In the digital case, one approach is to count how many horizontal, vertical, and diagonal steps are taken between pixel centers when scanning the boundary. Each type of step can then be assigned a weight. We wish to point out that the weight 1 for horizontal and vertical steps

463



A(c).

(1)

c ∈ ST

Treating the configurations one by one has been said to be “possible but tedious and error-prone” [6]. We, however, found it more convenient to create a lookup table (see Appendix A, Table III) with entries for each of the 256 configurations (this is done only once anyway). Hence, we do not need to take rotation, symmetry, or complementary configurations into account when scanning the tshell. Recall that m-cubes completely inside or completely outside

Fig. 3. (a) Examples of configurations of Marching Cubes, i.e., m-squares, in two dimensions. (b) Pixels inside the object are gray-colored. The boundary pixels are marked with 䊉. Only the m-squares indicated by dashed lines need to be stored in the t-shell. The boundary defined by the m-squares is delineated with solid lines.

464

I. Nyström et al. / Computerized Medical Imaging and Graphics 35 (2011) 460–471

the object are not stored in the t-shell; they will contribute neither to the area nor to the volume measurements. This is the beauty of closed (Jordan) surfaces. We note that the surface area computation may be further simplified by counting the number of appearances of each of the 256 configurations, and at the end multiplying the entries in this histogram by the surface area contributions by the individual configurations and adding up the result. However, this approach is not possible for the computation of the enclosed volume, where the current x co-ordinate is needed in the computations, as described below. 3.2. Enclosed volume One of the earliest works on estimation of volume enclosed by a given surface denned in a digital scene is based on a polyhedral approximation of the surface [31]. In this approach, the surface is utilized to decompose the slices of the object into polyhedra for which the volume can be computed analytically [31]. The slice-byslice nature of their approach was shown to be unnecessary in [4], where a true 3D digital approach was taken, leading to a tremendous simplification of the approach as well as of the computations. This approach used a digital version of the Gauss’ theorem: the integral of the normal component of any vector over any closed surface equals the integral of the divergence of the vector over the volume enclosed. In [32–35], this approach is utilized to compute volumes enclosed by surfaces that are represented as a set of polygons. In these approaches, the volume calculations involve three summations along the x, y, and z directions with quite complex weighting of the terms that reflect the orientation of the object. The volume is calculated by summing the vector product of the centroid, area, and normal of all triangles in the mesh surface. We observe that summations along three directions are not necessary. It is sufficient to sum along only one direction, exactly as in the purely digital surface case [4], and this is greatly facilitated by the t-shell data structure, as we describe below. In [36], an approach similar to ours is used involving accumulation of contributions from each facet, but their procedure is more complex due to splitting of facets into subfacets, as well as due to an integration over the facets. We demonstrate not only that summations along one direction are sufficient but also that all computations can be reduced to table lookup operations. 3.2.1. Underlying idea in two dimensions In the 2D case, the goal is to compute the area of the object from the given boundary. In the digital case, the area is accumulated as a sum of the x co-ordinates of the boundary elements (oriented pixel edges), where the contribution is −x if the boundary element is facing the background in the negative direction, +x if the boundary element is facing the background in the positive direction, and there is no contribution if the boundary element is facing the background in the ±y direction [4]. We noted in Section 3.1 that the perimeter estimation is more accurate if it is done from the m-squares than done from pixels. For consistency, we recommend that area estimation should also be done for the same representation in a quantitative analysis. In the triangulated case, in principle, the area contribution for a certain row of the object (with only two boundary elements per row) is simply the x co-ordinate of (the center of) the edge farthest to the right subtracted by the x co-ordinate of (the center of) the edge farthest to the left with some offset ı to adjust for the contribution of the m-square; see the 2D example in Fig. 3. The contributed area for the different m-square configurations can be precomputed and stored in a lookup table, as described below. The edge passing through the m-square is projected onto the y axis. The length of this projection is multiplied by the x co-ordinate

at the midpoint of the edge to determine the area contribution of the edge. The sign of the x component of the normal vector for the edge is sufficient to determine whether it is a positive or a negative contribution to the area. If the x component of the normal vector is positive (Fig. 3(a) top row), the incremental area contribution is positive. If this component is negative (Fig. 3(a) middle row), the incremental area contribution is negative. If the x component of the normal vector is equal to 0, i.e., the boundary element is parallel to the x axis (Fig. 3(a) bottom row), then such a configuration will not contribute to the area (but contributes only to the perimeter). From Fig. 3(b), it may appear that m-squares (2) and (3) should contribute to the area, but they will not. Their areas are accounted for by m-squares (1) and (4), which have negative and positive normal components in the x direction, respectively. 3.2.2. Extending to three dimensions In the purely digital case, where the surface is represented by the faces of the voxels at the boundary, the volume V(Sd ) enclosed by the digital surface Sd can be computed in a way similar to the 2D case by accumulating, while scanning the boundary [4], the −x and +x co-ordinate values of only those voxel faces that are orthogonal to the x axis. For the triangulated surface representation, we compute lookup table entries for each of the 256 m-cube configurations. Hence, both incremental area and incremental volume for each configuration will be known for this surface representation. We do not distinguish among rotational, symmetrical, and complementary cases. Actually, we take advantage of the fact that some of the configurations do not contribute to the volume, while a rotation of the same configuration may. This is specific to the scanning direction we have chosen, namely x (as opposed to y or z). Hence, besides being simple, it is a necessity to store all configurations instead of storing a reduced number of them. There are only 256 of them anyway. Let St be a set of triangles representing any closed surface obtained by using an MC type algorithm, and let ST be a t-shell representation of St . Let c be any m-cube in St of configuration ncc with a set of triangles (0), and let cx be the x co-ordinate of c. The signed volume contribution V(c) of c to the volume V(St ) enclosed by St can be thought of as consisting of two components (see Fig. 4) – contribution from outside of c and from inside of c: V(c) = ıo (c) + ıi (c).

(2)

These components can be further broken up into contributions from the individual triangles in (c): V(c) =



(ıo (t) + ıi (t)) =

t ∈ (c)



(A(pt )cx + ıi (t)),

(3)

t ∈ (c)

where A(pt ) is the area of the projection pt of triangle t onto the y–z plane. The total enclosed volume V(St ) is then given by



V(St ) =

(A(pt )cx + ıi (t)).

(4)

c ∈ ST t ∈ (c)

We note that both A(pt ) and ıi (t) are signed entities whose sign is determined by the sign of the x-component of the normal vector assigned to t. Observing that A(pt ) and ıi (t) are both completely independent of the individual m-cube and depend only on the MC configuration, (4) can be simplified to

V(St ) =

 c ∈ ST

=



c ∈ ST

⎛⎡ ⎝⎣



⎤

A(pt )⎦ cx +

t ∈ (c)

(ıA (c)cx + ıi (c)),

 t ∈ (c)

⎞

ıi (t)⎠ (5)

I. Nyström et al. / Computerized Medical Imaging and Graphics 35 (2011) 460–471

465

Fig. 4. For any triangle t (in a triangulated closed surface ST ) such as EFG with a non-zero x-component of its normal vector, we think of two components to its volume contribution: contribution ıo (t) from outside of the MC c, in this case, the volume of the prism E F G HGF; and contribution ıi (t) from inside, that is, the volume of the tetrahedron HGFE.

where ıA (denoting the sum in (5)) represents the sum of the signed areas of the projections (onto the y–z plane) of the triangles in (c). Since both ıA and ıi (c) depend only on the MC configuration, we precompute these entities and store them in a table (see Appendix A for details on the formulas and Table III), along with the area contributions A to the surface area (see (1)), for each of the 256 MC configurations. To compute A(St ) and V(St ), we simply use the t-shell ST , the table, and Eqs. (1) and (5). The property that the sign and the magnitude of the normal vector describes the contribution is exactly the same as for the digital surface approach [4]. The sign of the x component of the normal vector is sufficient to determine whether it is a positive or a negative contribution to the volume. If this component is positive, the incremental volume contribution is positive for that triangle, and negative if the component is negative. If this component is 0 (i.e., the triangle is parallel to the x axis), then the triangle will not contribute to the volume (but contributes only to the surface area). The volume will be accounted for by other m-cubes that have triangles with non-zero normal component in the x direction. A rather complex example is the configuration with four voxels inside and four voxels outside the object, as illustrated in Fig. 1(c). In this particular configuration, two of the four triangles yield a positive volume contribution and the other two triangles yield a negative volume contribution, due to the direction of their normal vectors. Another complexity is that the projected triangle areas partially overlap. However, this does not matter since the net volume and the net area contribution for a configuration are only a summation of the respective contributions for each triangle. In this specific case, the net volume contribution is actually 0, while for configurations resulting from rotations of this same configuration of voxels, it can be easily verified that this is not true.

produced are closed, even when scene resolution is not isotropic. See Appendix A for details. 4. Experiments and results We have estimated surface area and volume for both digital surfaces as well as triangulated surfaces for mathematically denned synthetic objects, physical phantoms, and real objects. This section describes our experiments and presents the results. 4.1. Mathematical phantoms Our first object is a ball. The digitization of a ball centered at (x0 , y0 , z0 ) ∈ R3 of radius r ∈ R is generated by the following equation:



f (x, y, z) =

2

2

2

1, if (x − x0 ) + (y − y0 ) + (z − z0 ) ≤ r 2 , 0, otherwise,

where (x, y, z) ∈ Z3

.

3.3. Lookup table We have computed a lookup table (Table 3) listing surface area and volume contributions for each of the 256 m-cube configurations specifically for the case where the triangle vertices are midway on the m-cube edges. See Appendix A, Table III, where the area of the triangles (A(c)), the area projected onto the y–z plane (ıA (c)), and the inside volume contribution (ıi (c)) for each configuration are listed for cubic voxels of size l × l × l. The lookup table is universal and can be used to compute area and volume of any surfaces created by the MC family of algorithms when the surfaces so

Fig. 5. Surface area and volume estimated by different methods divided by the true surface area and volume for 5310 digitized balls of increasing radius. The true surface area and volume were determined for the continuous balls from their known radii. The radius is expressed in voxel units, and therefore, it also directly indicates resolution. It can be noted that balls with radius less than 10 voxels give unreliable measurements.

466

I. Nyström et al. / Computerized Medical Imaging and Graphics 35 (2011) 460–471

Fig. 6. An MRI slice of a dry talus bone in a contrast medium (left), the same slice after segmentation (middle), and a surface rendering of the bone (right).

We have studied the performance of surface area and volume estimation methods for 5310 digitized balls with varying radii. In this work, we chose 177 different radii in the range from radius 0.45 to 79.65 voxels. For each size, 30 balls were generated by using (6) with randomized alignment in the digitization grid. Thereafter, surface area and enclosed volume were estimated by using the digital and the triangulated approaches. The corresponding mean values (over 30 data sets) are presented for each size in Fig. 5. It can be noted that balls with radius less than 10 voxels give unreliable measurements, cautioning us to be careful when measuring features of digitized objects that are small. The digital volume, V(Sd ), and the triangulated volume, V(St ), both coincide well with the true volume, V(St ) always being slightly smaller than V(Sd ), due to the “cutting of corners” in convex objects. The digital area, A(Sd ), is not shown in this plot, because of its large overestimate (close to 50%). This overestimate was also investigated recently in [37]. The triangulated area, A(St ), on the other hand, is a much more accurate estimate. The plot shows surface area convergence to an overestimate of 8.8% for large balls. Ways to improve the surface area estimates by choosing different triangulations and assigning optimized local weights have been reported [29]. We point out that our proposed efficient methods of surface area and volume computa-

tion will work fine with all such surfaces provided the methods conform to the MC family. The radius expressed in terms of the number of voxels has the same meaning as resolution. Imagine that the ball is 10 cm in diameter and is digitized at 10 different resolutions. If the voxel size is in the interval 0.2–2.0 mm in x, y, and z directions. Then the resolution varies from 250 voxels down to 25 voxels. For these balls, A(S) = 4r 2 = 314.16 cm2 , and V(S) = (4/3)r 3 = 523.60 cm3 . 4.2. Physical phantoms The physical phantoms we have used in this work are five dry bone specimens of the talus – one of the bones in the area of the human hindfoot. Their surfaces were sealed with methacrylate and varnish, so water cannot seep in. The true bone volumes were then determined experimentally from water displacement measurements. The water displacement box is filled with water until it flows from the spout coming out from it. The bones were then put inside one by one, and the runoff was collected and measured. Before scanning via MRI, the talus bones were suspended in the center of plastic boxes. Thereafter, the boxes were filled with the contrast medium OmniscanTM (a gadodiamide solution) diluted

Table 1 Volumes and surface areas derived from MR images of dry talus bones in a contrast medium. The true volumes are determined by measuring water displacement. Volumes are in milliliters. Areas are in cm2 . The two entries for the measured entities correspond to the two different orientations used for the bones while scanning. Talus

Surface area in cm2

Enclosed volume in ml Anisotropic

Isotropic

Anisotropic

Resolution in slice Isotropic

V(S)

V(Sd )

V(St )

V(Sd )

V(St )

A(Sd )

A(St )

A(Sd )

A(St )

No. 1

37.3

38.4 38.5

38.4 38.5

38.3 38.4

38.3 38.4

111.1 108.9

89.3 88.5

109.5 107.5

81.4 81.5

256 × 256

No. 2

29.1

29.9 29.6

29.9 29.6

29.9 29.5

29.9 29.5

93.6 93.3

79.2 79.4

92.7 92.2

69.3 68.3

512 × 512

29.1 29.1

29.0 29.1

29.0 29.0

29.0 29.0

93.0 91.2

74.3 74.0

92.0 89.6

67.3 66.7

256 × 256

25.2 25.2

25.2 25.2

25.1 25.1

25.1 25.1

83.1 84.0

70.7 72.0

82.1 83.2

62.1 61.6

512 × 512

24.7 24.7

24.7 24.7

24.6 24.6

24.6 24.6

81.5 81.5

65.3 65.8

80.5 80.5

60.3 59.9

256 × 256

23.6 23.6

23.6 23.6

23.5 23.5

23.5 23.5

81.4 82.5

69.2 70.3

80.5 81.5

61.4 61.2

512 × 512

23.5 23.3

23.5 23.2

23.4 23.2

23.4 23.2

81.1 82.2

65.5 66.2

79.8 80.9

60.3 60.0

256 × 256

35.2 35.1

35.2 35.0

35.1 35.0

35.1 35.0

106.0 109.6

85.8 88.5

104.1 108.1

78.6 79.4

256 × 256

No. 3

No. 4

No. 5

25.1

23.6

35.2

I. Nyström et al. / Computerized Medical Imaging and Graphics 35 (2011) 460–471

1:1000 in water. The phantoms were scanned twice (in different orientations) at resolution 256 × 256, yielding a voxel size of 0.31 mm × 0.31 mm × 1.20 mm. Three of the bones were also scanned twice (in different orientations) at resolution 512 × 512, yielding a voxel size of 0.16 mm × 0.16 mm × 1.00 mm. The original scenes have approximately 50 slices. The bones were thereafter segmented by using the interactive live-wire method [38]. A typical slice in the middle of a bone from its scene is shown before segmentation in Fig. 6, left, and after segmentation in Fig. 6, middle. Here, we consider the bone to be

467

a solid object; its porosity is not taken into account. A rendition of the same bone, created by using digital shell rendering, is displayed in Fig. 6, right. Digital and triangulated measurements of the segmented bones are presented in Table 1 both for the original anisotropic data and interpolated data, where a shape-based approach [39] has been used to obtain scenes with cubic voxels. It can be noted that the estimated volumes are slightly larger than the experimentally measured volumes, still being close to these values. The interpolated data seem to follow the original data well. The triangulated volumes are in general identical to the digital vol-

Fig. 7. Three-dimensional renderings of objects of different shapes and sizes segmented from CT and MRI data sets. Quantitative measurements are presented in Table 2.

468

I. Nyström et al. / Computerized Medical Imaging and Graphics 35 (2011) 460–471

Table 2 Real medical objects imaged at different resolutions. Volumes are in cm3 . Areas are in cm2 . Data sets

CT head, skin CT head, skull CT head, skull MR head, skin CT child skull CT dry skull CT knee MRA vessels

Enclosed volume in cm3

Surface area in cm2

V(Sd )

V(St )

A(Sd )

A(St )

2440.3 405.3 148.1 2600.8 313.9 556.4 92.3 86.8

2440.3 404.4 147.8 2607.9 313.8 555.8 91.4 82.5

2513.5 3214.0 1275.0 5210.5 2798.8 2949.1 1077.1 783.1

2094.6 2555.7 991.5 4038.3 2133.5 2316.0 1089.4 535.7

umes. The triangulated area measurements are never close to the digital area measurements, always smaller, but this is expected as a summation of voxel faces results in an overestimate. 4.3. Real objects We have tested our methods on real objects via clinical CT and MRI data sets obtained for various regions of the human anatomy. Surfaces were created from segmented scene data obtained by applying a simple gray-level threshold specifically tailored to the respective data set. We chose eight objects of different sizes and shapes, as described in Table 2. The triangulated volume estimates follow the digital volume estimates closely, while the triangulated area estimates are not as close, but this deviation can be expected, as described above. Fig. 7 shows renditions of some of the objects obtained by using t-shell rendering [16]. 5. Concluding remarks In this paper, we have presented an efficient method to simultaneously compute enclosed volume and surface area of objects represented by their 3D boundary surface in the form of a set of triangles. It is known that triangulated surface representations give a more accurate area estimate for the original object surface than digital surface representations. We have demonstrated that volume enclosed by triangulated surfaces can be computed efficiently in the same elegant way the volume enclosed by digital surfaces is computed by digital surface integration. We have shown that our triangulated shell data structure retains advantages and overcomes difficulties of both the digital and the triangulated approaches. We have also presented a lookup table that can be used for efficiently estimating the area and volume of closed surfaces produced by any method in the MC algorithm family. Note also that Eq. (12) in Appendix A shows that the central idea of this paper can be utilized to compute the volume enclosed by any closed triangulated or polygonized surface. Such surfaces may be produced, for example, by modified MC algorithms. If these algorithms produce triangles embedded in a cuboidal grid, then our table lookup method can be employed; otherwise, Eq. (12) will still yield an efficient method. Because of the way the surface is scanned and the normal components are utilized, not only solid objects can be handled, but also object cavities are dealt with correctly. Cavities are identified by the fact that surfaces bounding them yield a negative volume in our method. It is also possible to obtain measurements from several objects in an image in one scan. For this purpose, the described method is very efficient compared to connected component labelling, which is more costly than surface tracking [3]. In this work, the vertices of the surface triangles have been positioned midway between adjacent voxels. By employing a more sophisticated interpolation technique to determine the vertices, the resulting triangles may represent the true surface more

Image size in no. of voxels

Voxels size in mm

512 × 512 × 97 512 × 512 × 97 512 × 512 × 32 256 × 256 × 53 351 × 465 × 145 193 × 242 × 68 171 × 193 × 69 256 × 256 × 128

0.49 × 0.49 × 1.50 0.49 × 0.49 × 1.50 0.49 × 0.49 × 1.50 0.86 × 0.86 × 3.00 0.41 × 0.41 × 1.00 0.80 × 0.80 × 3.00 0.68 × 0.68 × 1.00 1.09 × 1.09 × 2.20

accurately and the estimates may become more accurate. The complexity of the method would increase, however, since precomputed lookup tables may not be feasible to the same extent. Note that, by discretizing the possible locations of the triangle vertices along the m-cube edges, instead of using real precision for each triangle vertex, the lookup table idea still becomes feasible, although the size of the table will increase. Our experiments with this idea indicate that the resulting improvements in the quality of the surface are negligible when the number of vertex positions considered is more than 8. Therefore our implementation in 3DVIEWNIX [40] (see www.mipg.upenn.edu) considers discretization of vertex locations by using 3 bits. By using even the hard (midpoint) triangles, an improvement is achieved compared to using the voxel faces, as was shown in Fig. 5. The need to generate the Marching Cubes configuration code while generating digital surfaces leads to future work. A study of the relationship among digital surfaces, the Marching Cubes representation, and t-shell is worth pursuing. Because of the efficiency of the t-shell representation, we do not find it necessary to perform triangle decimation which is widely pursued by methods dealing with triangulated surfaces. In fact, decimation will be counterproductive since the elegance of the digital nature of the triangulated surface will be lost. (This phenomenon is also observed in t-shell rendering in [16], which has demonstrated a considerable gain in speed of rendering triangulated surfaces in software (via t-shells) compared to rendering via hardware.) If decimation is done, our method can still compute surface area A(St ) and enclosed volume V(St ), although at the cost of foregoing the lookup table efficiency. We believe that our ideas from projective geometry can be utilized to compute the volume enclosed by surfaces approximated by any planar polygonal representation. This needs to be verified in a future study. Acknowledgments The authors are grateful to Dewey Odhner at MIPG for help in understanding the “secrets of bits” of 3DVIEWNIX and to Joakim Lindblad at Center for Image Analysis for help with the mathematical phantoms. This work was initiated during the first author’s visit to MIPG in 2001. The support of DHHS grant EB004395 is acknowledged. Appendix A. Given three vertices P1 , P2 , and P3 of any triangle t, the following entities are needed in our method. The centroid G for t, where its x component is given by Gx =

P1x + P2x + P3x . 3

(7)

Table 3 Lookup table for surface area and volume computation for the 256 different Marching Cubes configurations (ncc ) in the t-shell data structure. the net area of the triangles (A(c)), the net area projected onto the y–z-plane (ıA (c)), and the net contributing volume (ıi (c)) for each configuration are listed for cubic voxels, where the voxel size is assumed to be 1 × 1 × 1. A(c)

ıA (c)

ıi (c)

ncc

A(c)

ıA (c)

ıi (c)

ncc

A(c)

ıA (c)

ıi (c)

ncc

A(c)

ıA (c)

ıi (c)

ncc

A(c)

ıA (c)

ıi (c)

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50

0.0000 0.2165 0.2165 0.7071 0.2165 0.4330 0.7071 1.1495 0.2165 0.7071 0.4330 1.1495 0.7071 1.1495 1.1495 1.0000 0.2165 0.7071 0.4330 1.1495 0.4330 0.9236 0.9236 1.5731 0.4330 1.1495 0.6495 1.2990 0.9236 1.5731 1.3660 1.1495 0.2165 0.4330 0.7071 1.1495 0.4330 0.6495 1.1495 1.2990 0.4330 0.9236 0.9236 1.5731 0.9236 1.3660 1.5731 1.1495 0.7071 1.1495 1.1495

0.000 0.125 −0.125 0.000 −0.125 0.000 −0.500 −0.375 0.125 0.500 0.000 0.375 0.000 0.375 −0.375 0.000 0.125 0.500 0.000 0.375 0.000 0.375 −0.375 0.000 0.250 0.875 0.125 0.750 0.125 0.750 −0.250 0.375 −0.125 0.000 −0.500 −0.375 −0.250 −0.125 −0.875 −0.750 0.000 0.375 −0.375 0.000 −0.125 0.250 −0.750 −0.375 0.000 0.375 −0.375

0.0000 0.0208 −0.1042 0.0000 −0.1042 −0.0833 −0.3750 −0.1458 0.0208 0.1250 −0.0833 0.2292 0.0000 0.2292 −0.1458 0.0000 0.0208 0.1250 −0.0833 0.2292 −0.0833 0.0208 −0.3542 0.0000 0.0417 0.3542 −0.0625 0.3750 0.0208 0.3750 −0.1250 0.1458 −0.1042 −0.0833 −0.3750 −0.1458 −0.2083 −0.1875 −0.5208 −0.3750 −0.0833 0.0208 −0.3542 0.0000 −0.1042 0.1250 −0.3750 −0.2292 0.0000 0.2292 −0.1458

51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101

1.0000 0.9236 1.3660 1.5731 1.1495 0.9236 1.5731 1.3660 1.1495 1.4142 1.9656 1.9656 0.7071 0.2165 0.4330 0.4330 0.9236 0.7071 0.9236 1.1495 1.5731 0.4330 0.9236 0.6495 1.3660 1.1495 1.5731 1.2990 1.1495 0.4330 0.9236 0.6495 1.3660 0.9236 1.4142 1.3660 1.9656 0.6495 1.3660 0.8660 1.5155 1.3660 1.9656 1.5155 1.2990 0.7071 0.9236 1.1495 1.5731 1.1495 1.3660

0.000 −0.125 0.250 −0.750 −0.375 0.125 0.750 −0.250 0.375 0.000 0.625 −0.625 0.000 −0.125 0.000 −0.250 −0.125 −0.500 −0.375 −0.875 −0.750 0.000 0.375 −0.125 0.250 −0.375 0.000 −0.750 −0.375 0.000 0.375 −0.125 0.250 −0.375 0.000 −0.750 −0.375 0.125 0.750 0.000 0.625 −0.250 0.375 −0.625 0.000 −0.500 −0.375 −0.875 −0.750 −0.875 −0.750

0.0000 −0.1042 0.1250 −0.3750 −0.2292 0.0208 0.3750 −0.1250 0.1458 0.0000 0.2708 −0.3542 0.0000 −0.1042 −0.0833 −0.2083 −0.1042 −0.3750 −0.3542 −0.5208 −0.3750 −0.0833 0.0208 −0.1875 0.1250 −0.1458 0.0000 −0.3750 −0.2292 −0.0833 0.0208 −0.1875 0.1250 −0.3542 −0.2500 −0.5000 −0.3542 −0.0625 0.2500 −0.1667 0.2708 −0.1250 0.0208 −0.3542 −0.0833 −0.3750 −0.3542 −0.5208 −0.3750 −0.5208 −0.5000

102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152

1.0000 1.1495 0.9236 1.4142 1.3660 1.9656 1.5731 1.9656 1.1495 0.7071 1.1495 1.5731 1.2990 1.1495 1.5731 1.9656 1.1495 0.7071 1.3660 1.9656 1.5155 1.2990 1.9656 0.4330 1.2990 0.2165 0.2165 0.4330 0.4330 0.9236 0.4330 0.6495 0.9236 1.3660 0.7071 1.1495 0.9236 1.5731 1.1495 1.2990 1.5731 1.1495 0.7071 1.1495 0.9236 1.5731 0.9236 1.3660 1.4142 1.9656 1.1495

−1.000 −0.875 −0.375 0.000 −0.750 −0.375 −0.750 −0.375 −0.875 −0.500 −0.375 0.000 −0.750 −0.375 −0.750 −0.375 −0.875 −0.500 −0.250 0.375 −0.625 0.000 −0.625 0.000 −0.750 −0.125 0.125 0.250 0.000 0.125 0.000 0.125 −0.375 −0.250 0.500 0.875 0.375 0.750 0.375 0.750 0.000 0.375 0.500 0.875 0.375 0.750 0.375 0.750 0.000 0.375 0.875

−0.5000 −0.3542 −0.3542 −0.2500 −0.5000 −0.3542 −0.3750 −0.3542 −0.3542 −0.1250 −0.1458 0.0000 −0.3750 −0.2292 −0.3750 −0.3542 −0.3542 −0.1250 −0.1250 0.0208 −0.3542 −0.0833 −0.3542 0.0833 −0.2083 −0.0208 0.0208 0.0417 −0.0833 0.0208 −0.0833 −0.0625 −0.3542 −0.1250 0.1250 0.3542 0.0208 0.3750 0.2292 0.3750 0.0000 0.1458 0.1250 0.3542 0.0208 0.3750 0.0208 0.2500 −0.2500 0.0208 0.3542

153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203

1.0000 1.3660 1.1495 1.5731 1.1495 1.9656 0.7071 0.4330 0.6495 0.9236 1.3660 0.6495 0.8660 1.3660 1.5155 0.9236 1.3660 1.4142 1.9656 1.3660 1.5155 1.9656 1.2990 1.1495 1.2990 1.5731 1.1495 1.3660 1.5155 1.9656 1.2990 1.5731 1.1495 1.9656 0.7071 1.9656 1.2990 0.4330 0.2165 0.7071 0.9236 0.9236 1.4142 1.1495 1.3660 1.5731 1.9656 1.1495 1.5731 1.3660 1.9656

1.000 0.750 0.875 0.750 0.875 0.375 0.500 0.000 0.125 −0.375 −0.250 −0.125 0.000 −0.750 −0.625 0.375 0.750 0.000 0.375 0.250 0.625 −0.375 0.000 0.375 0.750 0.000 0.375 0.250 0.625 −0.375 0.000 0.750 0.875 0.375 0.500 0.625 0.750 0.000 0.125 0.000 0.125 −0.125 0.000 −0.375 −0.250 −0.750 −0.625 0.375 0.750 0.250 0.625

0.5000 0.2500 0.5208 0.3750 0.5208 0.0208 0.3750 −0.0833 −0.0625 −0.3542 −0.1250 −0.1875 −0.1667 −0.5000 −0.3542 0.0208 0.2500 −0.2500 0.0208 0.1250 0.2708 −0.3542 −0.0833 0.2292 0.3750 0.0000 0.1458 0.1250 0.2708 −0.3542 −0.0833 0.3750 0.5208 0.0208 0.3750 0.2708 0.5417 0.0833 0.1042 0.0000 0.0208 −0.1042 0.0000 −0.1458 −0.1250 −0.3750 −0.3542 0.2292 0.3750 0.1250 0.2708

204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255

1.0000 1.1495 1.1495 0.7071 1.1495 1.5731 1.3660 1.9656 1.5731 1.9656 1.9656 0.4330 1.2990 1.1495 1.5155 1.2990 1.1495 0.7071 1.2990 0.2165 1.1495 1.3660 1.5731 1.9656 1.2990 1.5155 1.1495 1.2990 1.5731 1.9656 1.9656 0.4330 1.1495 1.2990 0.7071 0.2165 1.0000 1.1495 1.1495 0.7071 1.1495 1.2990 0.7071 0.2165 1.1495 0.7071 1.2990 0.2165 0.7071 0.2165 0.2165 0.0000

0.000 0.375 −0.375 0.000 0.375 0.750 0.250 0.625 0.000 0.375 −0.375 0.000 0.750 0.875 0.625 0.750 0.375 0.500 0.000 0.125 −0.375 −0.250 −0.750 −0.625 −0.750 −0.625 −0.875 −0.750 0.000 0.375 −0.375 0.000 −0.375 0.000 −0.500 −0.125 0.000 0.375 −0.375 0.000 −0.375 0.000 −0.500 −0.125 0.375 0.500 0.000 0.125 0.000 0.125 −0.125 0.000

0.0000 0.1458 −0.2292 0.0000 0.2292 0.3750 0.1250 0.2708 0.0000 0.0208 −0.3542 0.0833 0.3750 0.5208 0.2708 0.5417 0.1458 0.3750 −0.0833 0.1042 −0.1458 −0.1250 −0.3750 −0.3542 −0.3750 −0.3542 −0.3542 −0.2083 0.0000 0.0208 −0.3542 0.0833 −0.2292 −0.0833 −0.1250 −0.0208 0.0000 0.1458 −0.2292 0.0000 −0.2292 −0.0833 −0.1250 −0.0208 0.1458 0.3750 −0.0833 0.1042 0.0000 0.1042 −0.0208 0.0000

I. Nyström et al. / Computerized Medical Imaging and Graphics 35 (2011) 460–471

ncc

469

470

I. Nyström et al. / Computerized Medical Imaging and Graphics 35 (2011) 460–471

The normal vector for t is −−→ −−→ N = P1 P2 × P1 P3 .

(8)

The area of t is given by A(t) =

1 −−→ −−→ |P1 P2 × P1 P3 |, 2

(9)

where |·| denotes the length (the magnitude) of the vector. Note that the sign of the area will become positive or negative depending on whether the vertices are numbered clockwise or counterclockwise. The signed area A(pt ) of the projection of the triangle onto the y–z plane is A(pt ) = A(t)

Nx 1 = Nx , |N| 2

(10)

where the last expression certainly simplifies the computations. The two components of the volume contributed to V(St ) by a triangle in St are given by (the variables are as in Eq. (3)) ıo (t) = A(pt )cx , . ıi (t) = Gx A(pt ) − ıo (t)

(11)

We note that, even when St is any triangulated closed surface (not necessarily produced by an algorithm in the MC family), the following expression can be used to compute the enclosed volume: V(St ) =



Gx A(pt )

(12)

t ∈ St

Note that (12) is valid even when St is any polyhedron. The above equations are valid if the voxels are of size 1 × 1 × 1, i.e., the m-cubes are real cubes. To incorporate the proper resolution and handle anisotropic data, a scale vector s = (sx , sy , sz ) is needed. sx and sy represent the pixel size and sz denotes the slice spacing, i.e., the voxel size is sx × sy × sz . The normal vector N is modified by scaling the components of the vectors in (8): sx (P1x − P2x ), sy (P1y − P2y ), and so on. Also, when computing the volume under the triangle (12), the scale in the x direction must be used to modify the x component of the centroid: V(t) = sx Gx A(p). Eqs. (7)–(10) are general for any triangle. We have computed a lookup table (Table 3) listing area and volume contributions for each m-cube configuration specifically for the case where the triangle vertices are midway on the m-cube edges. Eqs. (9)–(11) are utilized for each triangle in the specific configuration and the results are summed which yields the three net contributions stored in the lookup table. The projected area ıA (c) times the current x coordinate plus the volume inside ıi (c) under all triangles gives the volume contribution for an m-cube. The x co-ordinate should be scaled with the voxel size in the x direction when computing the contribution for an m-cube configuration. The lookup table can be used to compute area and volume of any surfaces created by the MC family of algorithms when the surfaces so produced are closed. Note also that 24 among the 256 MC configurations contribute zero volume! References [1] Udupa JK, Herman GT, editors. 3D imaging in medicine. 2nd ed. Boca Raton, FL: CRC Press LLC; 2000 [ch. 1]. [2] Thurfjell L, Bengtsson E, Nordin B. A new three-dimensional connected components labeling algorithm with simultaneous object feature extraction capability. CVGIP: Graphical Models and Image Processing 1992;54(July):357–64. [3] Udupa JK, Ajjanagadde VG. Boundary and object labelling in threedimensional images. Computer Vision, Graphics, and Image Processing 1990;51(September):355–69. [4] Udupa JK. Determination of 3-D shape parameters from boundary information. Computer Graphics and Image Processing 1981;17:52–9. [5] Wyvill G, McPheeters C, Wyvill B. Data structures for soft objects. The Visual Computer 1986;2(7):227–34. [6] Lorensen WE, Cline HE. Marching Cubes: a high resolution 3D surface construction algorithm. In: Proceedings of the 14th ACM SIGGRAPH on Computer Graphics, vol. 21, no. 4. 1987. p. 163–9.

[7] Lachaud J-O, Montanvert A. Digital surfaces as a basis for building iso-surfaces. In: Proc. of 5th IEEE Int. Conference on Image Processing (ICIP’98), vol. 2. 1998. p. 977–81. [8] Lachaud J-O, Montanvert A. Continuous analogs of digital boundaries: a topological approach to iso-surfaces. Graphical Models 2000;62:129–64. [9] Gordon D, Udupa JK. Fast surface tracking in three-dimensional binary images. Computer Vision, Graphics, and Image Processing 1989;45(February):196–214. [10] Udupa JK, Odhner D. Fast visualization, manipulation, and analysis of binary volumetric objects. IEEE Computer Graphics and Applications 1991;11(November):53–62. [11] Grevera GJ, Udupa JK, Odhner D. An order of magnitude faster isosurface rendering in software on a PC than using dedicated general purpose rendering hardware. IEEE Transactions on Visualization and Computer Graphics 2000;6(October):335–45. [12] Udupa JK. Multidimensional digital boundaries. Graphical Models and Image Processing 1994;56(July):311–23. [13] Artzy E, Frieder G, Herman G. The theory, design, implementation and evaluation of a three-dimensional surface detection algorithm. Computer Graphics and Image Processing 1981;15:1–24. [14] Udupa JK, Srihari SN, Herman GT. Boundary detection in multidimensions. IEEE Transactions on Pattern Analysis and Machine Intelligence 1982;4:41–50. [15] Morgenthaler DG, Rosenfeld A. Surfaces in three-dimensional digital images. Information and Control 1981;51(3):227–47. [16] Grevera GJ, Udupa JK, Odhner D. T-shell rendering. In: Mun SK, editor. Medical imaging 2001: visualization, display, and image-guided procedures, Proc. SPIE 4319, p. 413–25; 2001. [17] Udupa JK, Odhner D. Shell rendering. IEEE Computer Graphics and Applications 1993;13(November):58–67. [18] Mullikin J, Verbeek P. Surface area estimation of digitized planes. Bioimaging 1993;1(1):6–16. [19] Lindblad J. Surface area estimation of digitized planes using weighted local configurations. In: Nystrom I, Sanniti di Baja G, Svensson S, editors. Proceedings of Discrete Geometry for Computer Imagery (DGCI 2003), vol. 2886 of Lecture Notes in Computer Science, Springer-Verlag, p. 348–57; 2003. [20] Windreich G, Kiryati N, Lohmann G. Voxel-based surface area estimation: from theory to practice. Pattern Recognition 2003;36(11):2531–41. [21] Gerstner T, Pajarola R. Topology preserving and controlled topology simplifying multiresolution isosurface extraction. In: Proceedings of the Visualization. 2000. p. 259–66. [22] Wood Z, Hoppe H, Desbrun M, Schroder P. Removing excess topology from isosurfaces. ACM Transactions on Graphics 2004;23(April):190–208. [23] Nystrom I, Udupa JK, Grevera GJ, Hirsch BE. Area of and volume enclosed by digital and triangulated surfaces. In: Mun SK, editor. Medical imaging 2002: visualization, image-guided procedures, and display. Proc. SPIE 4681, p. 669–80; 2002. [24] Van Gelder A, Wilhelms J. Topological considerations in isosurface generation. ACM Transactions on Graphics 1994;13(4):337–75. [25] Dürst MJ. Letters: additional reference to “Marching Cubes”. In: Proceedings of ACM SIGGRAPH on Computer Graphics, vol. 22, no. 2. 1988. p. 72–3. [26] Nielsen GM, Hamman B. The asymptotic decider: Resolving the ambiguity in Marching Cubes. In: Proceedings of IEEE visualization’91. IEEE Computer Society Press; 1991. p. 83–90. [27] Kenmochi Y, Kotani K, Imiya A. Marching Cubes method with connectivity. In: Proc. of the IEEE international conference on image processing (ICIP’99). 1999. p. 361–5. [28] Kenmochi Y, Klette R. Surface area estimation for digitized regular solids. In: Latecki LJ, Melter RA, Mount DM, Wu AY, editors. Vision geometry IX, Proc. SPIE 4117, p. 100–11; 2000. [29] Lindblad J, Nystrm I. Surface area estimation of digitized 3D objects using local computations. In: Braquelaire A, Lachaud J-O, Vialard A, editors. Proceedings of Discrete Geometry for Computer Imagery (DGCI 2002), vol. 2301 of Lecture Notes in Computer Science. Springer-Verlag; 2002. p. 267–78. [30] Groen FCA, Verbeek PW. Freeman code probabilities of object boundary quantized contours. Computer Graphics and Image Processing 1978;7: 391–402. [31] Cook LT, Cook PN, Lee KR, Batnitzky S, Wong BYS, Fritz SL, et al. An algorithm for volume estimation based on polyhedral approximation. IEEE Transactions on Biomedical Engineering 1980;27(September):493–500. [32] Eberly D, Lancaster J, Alyassin A. On gray scale image measurements: II. Surface area and volume. CVGIP: Graphical Models and Image Processing 1991;53(November):550–62. [33] Lancaster JL, Eberly D, Alyassin A, Downs III JH, Fox PT. A geometric model for measurement of surface distance, surface area, and volume from tomographic images. Medical Physics 1992;19(March):419–31. [34] Alyassin AM, Lancaster JL, Downs III JH, Fox PT. Evaluation of new algorithms for the interactive measurement of surface area and volume. Medical Physics 1994;21(June):741–52. [35] Hughes SW, D’Arcy TJ, Maxwell DJ, Saunders JE, Ruff CF, Chiu WSC, et al. Application of a new discreet form of Gauss’ theorem for measuring volume. Physics in Medicine and Biology 1996;41(September):1809–21. [36] Battle XL, Cunningham GS, Hanson KM. 3D tomographic reconstruction using geometrical models. In: Hanson KM, editor. Medical Imaging 1997: Image Processing. Proc. SPIE 3034, p. 346–57; 1997. [37] Rajon DA, Patton PW, Shah AP, Watchman CJ, Bolch WE. Surface area overestimation within three dimensional digital images and its consequence for skeletal dosimetry. Medical Physics 2002;29(May):682–93.

I. Nyström et al. / Computerized Medical Imaging and Graphics 35 (2011) 460–471 [38] Falco AX, Udupa JK, Samarasekera S, Sharma S, Hirsch BE, Lotufo RA. Usersteered image segmentation paradigms: live wire and live lane. Graphical Models and Image Processing 1998;60(May):233–60. [39] Grevera GJ, Udupa JK. Shape-based interpolation of multidimensional greylevel images. IEEE Transactions on Medical Imaging 1996;15(6):881–92.

471

[40] Udupa JK, Odhner D, Samarasekera S, Goncalves R, Iyer K, Venugopal K, et al. 3dviewnix: an open, transportable, multidimensional, multimodality, multiparametric imaging software system. In: Medical Imaging vol. 1994. Proc SPIE 2164, p. 58–73; 1994.