- Email: [email protected]
Polygonization of volumetric skeletons with junctions
Computer-Aided Design 45 (2013) 822–828
Contents lists available at SciVerse ScienceDirect
Computer-Aided Design journal homepage: www.elsevier.com/locate/cad
Polygonization of volumetric skeletons with junctions Takashi Michikawa ∗ , Hiromasa Suzuki Research Center for Advanced Science and Technology, The University of Tokyo, Japan
article
info
Keywords: Volumetric data Medial axis Non-manifold meshes Reverse engineering
abstract This article presents a method for polygonizing skeleton structures computed from volumetric data, mainly of a thin-plate object. The main idea of our method is to simplify the topological structure of volumetric skeletons by sub-sampling. Given volumetric skeletons and distance fields of binarized volumes, we first apply sub-sampling of the skeleton voxels so that the junction and boundary voxels are preferentially sampled. By using sampled points as Voronoi sites, we compute Voronoi diagrams on the skeleton voxels. Polygonization can be computed by using connectivity of the diagram. This approach can build correct junctions where the skeleton branches in several directions, whereas conventional voxelbased methods tend to result in small cavities around them. This paper also demonstrates several results computed from CT-scanned engineering objects. © 2011 Elsevier Ltd. All rights reserved.
1. Introduction This paper outlines a method for converting volumetric skeleton structures or medial axes of thin-plate objects to polygonal meshes for use in digital engineering applications. Thin-plate engineering objects are made of metal plates, and are formed by stamping and welding; they are often found in cars and home electronics. Industrial companies measure such objects as volumetric images using X-ray CT scanners or nondestructive scanning devices (Fig. 1(a)). Since each voxel has a CT value that identifies its material, we can extract volumetric data with certain threshold values from CT images. Companies need polygonal meshes of these CT images for use in accelerating engineering processes such as comparison with CAD models and CAE simulation. It is well known that polygonal meshes can be obtained by isosurface contouring of volumetric data [1,2]. This is effective for solid objects, but can create closed surfaces with thin-plate objects (Fig. 1(b)). The center-line structure of the objects shown in Fig. 1(c) are preferable because common CAD systems represent thin-plate objects as open surfaces [3]. Thin-plate objects are usually welded to other objects for reinforcement. The medial surfaces of such objects become nonmanifold, and their welded parts become junctions. Since our focus here is on engineering application, these forms must be reconstructed correctly. One of the challenges involved is capturing non-manifold junctions from CT images correctly as shown in Fig. 1(c).
∗
Corresponding author. Tel.: +81 3 5452 5183. E-mail addresses: [email protected] (T. Michikawa), [email protected] (H. Suzuki). 0010-4485/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.cad.2011.06.003
The computation of medial surfaces from volumetric data has been extensively discussed in previous studies [4,3,5]. The basic approach first creates medial voxels from input volumetric data, and polygonal meshes are then created based on the resulting voxel connectivity. However, this approach creates small cavities around junctions even though no such cavities exist in the input data. This is due to the ambiguity of triangulation from voxel connectivity. Prohaska and Hege [4] introduced a look-up table for non-manifold topology polygonization. However, local voxels are still checked due to the existence of bad cases of cavities. Liu et al. proposed a robust algorithm by pruning the input models as cell complexes [6]. Their method can create topologically correct surfaces, however their polygonization are axis-aligned and it is difficult to represent higher-valence junctions. We present a method of medial surface reconstruction from volumetric data of thin-plate objects including junctions. Given medial voxels with their distance to the boundary, we first apply sub-sampling to these voxels so that junction and boundary voxels can be preferentially selected. This computation is performed by covering the medial voxels with spherical supports [7]. We then perform Voronoi diagram computation for medial voxels using the sampled points as Voronoi sites to build vertex connectivity, from which triangular meshes are then created. This approach offers various advantages in medial surface reconstruction from volumetric data. In particular, spherical supports in the sub-sampling phase remove unnecessary voxels around junctions. This contributes to the absorption of various types of artifacts arising in the junctions regions, including branches and small cavities. In addition, our method can be used for computing not only medial surfaces but also 1D skeletons and their coupling is also possible. This paper consists of five sections. We review related work in Section 2, the proposed algorithm is introduced in Section 3, we
T. Michikawa, H. Suzuki / Computer-Aided Design 45 (2013) 822–828
823
polygonization of non-manifold topology. However, local voxels are still checked due to the existence of bad cases of cavities. Liu et al. recently proposed a method for computing medial structures (wires and skeletons) by pruning the input models as cell complexes [6]. However, their polygonization depends on the topology of cell complex and it is usually insufficient for representing higher-valence junctions.
Fig. 1. CT images of welded aluminum and its polygonal meshes. The resolution of volumetric data (a) is 350 × 330 × 220 and its isosurface (b) and medial surface have 658,648 and 7539 faces (Computation Time: about 1 min.) respectively. Note that the junction edges (red) of the medial surfaces created using our method (c) are correctly reconstructed. Lower row shows corresponding illustrations. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
outline the experimental results and discuss them in Section 4, and Section 5 concludes the paper. 2. Related work Our issue is related to a medial axis whose formulation was first introduced in [8]; that is, the set of center points of a collection of maximal inscribed spheres. Many approaches involving MATs (medial axis transforms) have been studied over a period of four decades in the image processing and geometric modeling community. This section reviews several approaches related to our objective. Volume-based methods. Volume-based approaches compute medial voxels from binarized volumetric data. This produces a volumetric version of the medial axis. A typical strategy is topological thinning, or voxel pruning, which leaves the topology of the input data unchanged (such voxels are called ‘‘simple voxels’’) and preserves the ends of surface voxels. This process is iteratively applied until no more voxels can be removed. Sequential thinning removes simple voxels one by one. This approach preserves the topology of the original data, but the order of pruning may preserve excess simple voxels and cause many branches or bumps. Although several solutions to this have been discussed (e.g., sub-cycles [9], distance fields [10] and relaxed conditions [11]), the issue has not yet been completely resolved. Parallel thinning removes all simple voxels on the boundary of an object simultaneously. Tsao and Fu [12] formulated conditions for building smooth surfaces, but the topology of the input model may be lost with this technique. The resulting surfaces are relatively smooth, but it is difficult to preserve topology because many voxels are removed at once. There are also topologypreserving methods for parallel thinning [13,14], but these cause a deterioration in voxel quality. Prohaska and Hege [4] used the geodesic distance between two closest boundary voxels to determine pruning priority. This global metric is more stable for noisy data, and hardly any unnecessary branches are generated. Accordingly, it is used in a number of engineering applications [3,5,15]. Polygonal meshes can be computed from medial voxels. For instance, we can cast this issue to surface reconstruction problems (e.g., [16]). However, non-manifold surfaces cannot be handled. Another approach builds triangular meshes using voxel connectivity, but this method produces ambiguous cases of triangulation from such connectivity, and small cavities may be created. Prohaska and Hege [4] introduced a look-up table for the
Polygon-based methods. Alternative approaches involve computing medial axes from polygonal meshes. Although an exact solution can be computed by solving algebraic equations, this is not practical for large models, and most polygon-based methods deal with approximate solutions of MAT. One popular approach is the use of Voronoi diagrams [17–20]. The structures of medial axes computed using the above methods are usually complex, and may have small branch surfaces. For simplified medial axes, these are usually removed using certain criteria. One popular criterion is the separation angle of sheets [17,21]. Dey and Zhao [20] used additional conditions to improve enhancement of the quality of medial surfaces. Sud et al. [22] also improved the medial axis simplification method [21] so that the homotopy of the input data is preserved. Recently, the scale axis transform [23,24] algorithm was introduced to simplify medial surfaces by scaling the medial ball of the objects. However, they depend on definition of the medial axis. For instance, it is difficult to represent a higher-valence junctions because they may be represented as a set of lowervalence junctions. When these techniques are applied to our issue, polygonal meshes must be computed using isosurface extraction methods. This is straightforward, but the resulting meshes have several problems including numerical issues for complex models and artifacts caused by noise and blur in CT images. Accordingly, polygon-based approaches are not practical in achieving our objective. 3. Medial surface reconstruction from medial voxels The goal of this algorithm is to polygonize a skeleton structure with junctions of volumetric data. Fig. 2 shows an overview of our method. In this article, we suppose that medial voxels (from Fig. 2(a)–(d)) are already computed by conventional methods (e.g. [4,3,5,15]) (Fig. 3). We first compute binary volumes generated by thresholding volumetric data. The threshold value is a user-defined parameter. We next compute distance fields D = {di } from the binarized volumes (Fig. 3(c)). Once distance fields are computed, medial voxels M = {vi } are computed by evaluating the geodesic distance between the two closest boundary points neighboring to each voxel. Fig. 3(d) shows medial voxels computed from CT images in Fig. 3(a). The distance value at the medial voxels are drawn in color map in the figure. A method for converting medial voxels to polygons consists of three steps: Sub-sampling—distributes spherical supports on voxels iteratively so that they completely cover the input medial voxels (Fig. 2(e)–(g)). Voronoi diagram—computes a Voronoi diagram on voxels to build mesh edges (Fig. 2(h)). Non-manifold meshing—computes polygonal meshes by finding three vertices connected to each other and applies topological cleaning to remove excess non-manifold edges (Fig. 2(i)).
824
T. Michikawa, H. Suzuki / Computer-Aided Design 45 (2013) 822–828
a
b
c
a
b Sub-sampling
d
e
f
g
h
i
Fig. 4. A 2D comparison of meshing results with random distribution (a) and topology-dependent distribution (b). In each figure, the voxels shown in red denote sampled points, and dotted circles represent the spherical supports of voxels. The red lines show medial surfaces connected by a Voronoi diagram on medial voxels. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
(a) #component = 1. Fig. 2. An overview of our method. (a) CT images. (b) Binarized volumes (FG and BG voxels are drawn in black and white respectively). (c) Distance field of (b). (d) Medial voxels. Each voxel has distance to the closest boundary voxel of the object. (e)–(g) Sub-sampling of medial voxels. Neighboring voxels (gray) are removed. (h) Voronoi diagram on medial voxels using sampled points as Voronoi sites. (i) Polygonization.
a
b
c
d
Fig. 3. Medial voxel extraction from binary volumes. (a) Input CT images. (b) Binary volume. (c) Distance field of (b) (Cut-away view). (d) Medial voxels. Distance value is assigned to each voxel. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
3.1. Sub-sampling of medial voxels The first step of the algorithm involves medial voxel subsampling. Given a set of medial voxels M, we select one voxel vi from M and set it to vertex list V . Then, we remove neighboring voxels vj such that ∥vi − vj ∥ < α di from M, where vi denotes the position of vi , and α denotes the radius scaling factor. Note that a larger value for α usually makes robust triangulation at very thin parts of the object because the radii of such parts are usually small and may cause bad sampling of medial voxels. On the other hand, too large α also remove many voxels neighboring to vi and it causes ill shaped triangulation. Accordingly, a value of α = 2 usually makes good results. Iteration is performed until M is empty and then all the voxels are covered by spherical supports.
(b) #component = 2.
(c) #component = 3.
Fig. 5. End point classification. We count the number of components in the 5 × 5 × 5 masks with the exception to 3 × 3 × 3 mask (shaded region). If the number of components is 2, the voxel is not end point (b). Otherwise, the point is end point ((a) and (c)).
An important point here is how to find these vi voxels from M, because their distribution affects mesh topology directly. Fig. 4 illustrates a 2D example. Random distribution may not involve placement at junctions, and sampling points around junctions may be connected to each other, thereby creating small cavities as shown in Fig. 4(a). On the other hand, in Fig. 4(b), spheres are distributed at junction and boundary voxels first, and correct polygonization is obtained. We therefore introduce the topologydependent distribution as described below. Topology criterion control for sampling is implemented so that points are well distributed on junctions and boundary edges. This contributes to fastening low-dimensional topology features such as junctions and boundaries. We use it here as the primary priority. A similar idea can be found in the bubble mesh algorithm [25] for polygonal meshes with explicit topology, and we apply similar approach to volumetric data. The basic classification algorithm was studied in [26] (see Appendix for a brief summary). In addition, we also consider end points of junctions or boundary, because such points must be sampled first. After the classification, we extract wire structures consisting of boundary and junction voxels. We then count the number of components of the wire structure within 5 × 5 × 5 mask with the exception to 3 × 3 × 3 mask (Fig. 5). If the number of components is 2, the point is a normal boundary (Fig. 5(b)), otherwise it is an end point (Fig. 5(a) and (c)). We also consider the distance values of voxels as a secondary criterion. This is because medial voxels have relatively larger distances than neighboring voxels, and it is natural that such voxels should be selected first. Finally, the topology-dependent sub-sampling algorithm is summarized as below. 1. Find the vi with the top priority in M and copy it to V . 2. Remove vj from M such that ∥vi − vj ∥ < α di . 3. Iterate 1 and 2 until M becomes empty. Fig. 6 shows a result of topology classification of the medial voxels (Fig. 3(d)) and the sub-sampling result (Fig. 6(b)). Color of each sphere means a sequence of sub-sampling (red sphere is sampled first). We can see medial voxels on non-manifold features are preferentially selected.
T. Michikawa, H. Suzuki / Computer-Aided Design 45 (2013) 822–828
(a) Topology classification.
(b) Sampling result.
Fig. 6. Topology-dependent sub-sampling. (a) Topological classification of medial voxels. Surface, boundary, junction and end point voxels are drawn in gray, blue, red and green respectively. (b) Sub-sampling result. A center point of each sphere is sampled. Colors of each sphere shows sequence of sub-sampling (red spheres are selected first and blue spheres are selected later.) (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
825
Fig. 9. 2D illustration of neighboring voxels in surface voxels. Note that the planer graph can be defined between F (v) and B0 (v) (indicated by the bold red line). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
(a) Voronoi diagram.
(b) Polygonization.
Fig. 7. Voronoi diagram definition on medial voxels. Identically colored voxels belong to the same sites.
(a) Non-manifold surface.
(b) Manifold decomposition.
Fig. 8. Medial voxels can be separated into three simple sheets consisting of boundary and surface voxels.
3.2. Building a Voronoi diagram on medial voxels We next build mesh edges from the set of sampled vertices V . In our target objects of thin-plate structures, we can assume that medial voxels can be decomposed into a collection of subsets, each of which is an analog of a two-manifold surface (See below for details). This subset therefore inherits the properties of a twomanifold surface. In particular, it allows us to define Voronoi diagrams on this subset. Fig. 7 illustrates the formulation of Voronoi diagrams on medial voxels. Given sites V = {vi } ⊆ M, we define the region R(vi ) as follows: R(vi ) = x|d(vi , x) ≤ d(vj , x), ∀j ̸= i ,
(1)
where d(vi , x) denotes the shortest path length between vi and x in the voxel space. In practice, we compute it using the wavefront propagation of distance from the sites. Fig. 10(a) shows the result of the diagram for medial voxels M. We define an edge (vi , vj ) if there is a neighboring voxel pair (x, y) such that x ∈ R(vi ) and y ∈ R(vj ) in 26-adjacent. Validity of polygonization. Medial voxels in our method can be divided into surface, boundary and junction types by topology classification (Fig. 8(a)). They can be separated into a set of sheets consisting of boundary and surface voxels by junction voxels
(c) Polygonization result.
(d) Mesh cleaning result.
Fig. 10. Polygonization with a Voronoi diagram and mesh cleaning.
(Fig. 8(b)). We can also consider that the Voronoi diagram is defined on a surface composed of the faces of foreground voxels shared by background voxels. It should be noted here that a voxel is regarded as a cubic volume. Based on the definition of topological classification, each surface voxel v in a sheet separates its 26-neighboring voxels N (v) into two background components Bi (v) (i = 0, 1) using one foreground component F (v) (Fig. 9). When we focus on B0 (v) and F (v), the boundary surface between their voxels is topologically equivalent to a disk, because B0 (v) is 6-connected and its surface must be two-manifold. This indicates that every surface piece at a voxel on the sheet is topologically equivalent to a disk and thus the total surface of the sheet is topologically equivalent to a twomanifold surface. The Voronoi diagram is therefore defined on surface voxels and we can generate a polygonal mesh by taking the dual of the Voronoi diagram.
826
T. Michikawa, H. Suzuki / Computer-Aided Design 45 (2013) 822–828
(a) Input.
(b) Sub-sampling.
(c) Result.
Fig. 11. A result for side member model (Upper, 400 × 400 × 576 voxels, 7500 faces, 2.79 min.) and door panel model (Lower, 1370 × 880 × 500 voxels, 21,748 faces, 21.71 min.)
(a) Input.
(b) Sub-sampling.
(c) Result.
Fig. 12. A result for Knot model (Upper, Wire only, 512 × 512 × 512 voxel, 254 wire edges, 62.80 s) and coupled model (Lower, Wire and surfaces, 512 × 512 × 512 voxels, 2267 faces, 38 wire edges, 79.50 s.)
3.3. Non-manifold meshing from vertex connectivity
3.4. An extension to wire–skeleton structure extraction
In the final step, polygonal meshes are constructed from the vertices and edges computed in the previous steps (Fig. 10(a)). This computation is performed by finding triples (vi , vj , vk ) that are connected to each other (Fig. 10(b)). This approach creates non-manifold meshes with no small cavities caused by invalid sampling, but may also create overlapping triangles on surface voxels (Fig. 10(c)) for similar reasons related to the degeneration of Delaunay triangulation. Such triangles must be eliminated so that necessary ones around junctions are not removed. Based on medial voxel topology classification, we decompose all vertices into manifold vertices (on surfaces) and non-manifold vertices (on boundaries or junctions), and then apply manifold cleaning [7] to the manifold vertices only. Simply speaking, this method performs iteration to find a planar graph from the one-ring neighborhood of each vertex such that the number of vertices in the graph is maximized. Triangles excluded from the planar graph are removed from the mesh. Fig. 10(d) shows the results of cleaning. Note that some overlapping triangles may still be left if all vertices belong to junctions. Such issues are generally seen when junctions are very close.
An extension of our method to 1D skeleton extraction is straightforward. Moreover, the 1D skeleton or the wire structure can also be computed at once. Fig. 12 show results for a knot model and a volumetric data involving surface and wire structures. We can see both of results have correct topology. The procedure is almost same as the original described in this section. By Voronoi diagrams on medial voxels, we already have connectivity of wire edges. We just need to manage these edges so that they cannot be removed by surface cleaning. 4. Result and discussion The experimental results of meshing are shown in Figs. 11 and 12. Each figure shows input binary volumes, the sub-sampling result and the final meshing result. Figs. 1 and 13 are also obtained by our method. All examples were tested on a 3.16 GHz computer, and the computation times are provided in the caption of each figure. There is still room for improvement here, as our prototype is not yet optimized. The current bottleneck lies in Voronoi diagram creation and mesh cleaning. Since Voronoi diagrams compute using a multi-source Dijkstra’s algorithm, the computation cost is
T. Michikawa, H. Suzuki / Computer-Aided Design 45 (2013) 822–828
O(n log n). Note that although connecting spherical supports [7] may create polygons fast, the technique truncates topology information hidden in medial voxels, whereas the Voronoi-based method considers such topological information. The non-manifold mesh data structure used in cleaning is based on an IndexedFaceSet (a.k.a. Wavefront OBJ format), and takes a long time to search neighboring vertices. An advanced data structure for non-manifold meshes (e.g., [27]) is expected to improve the level of performance. Handling large objects is an important issue because CT images of engineering objects are usually large. Since our method manages medial voxels by octree, our method is efficient for memory usage and can handle large objects as shown in the lower row in Fig. 11. Using topology-dependent sampling, medial surfaces with junctions were reconstructed from CT-scanned engineering objects. In addition, sub-sampling with spherical supports contributes to the recovery of non-manifold junctions. Fig. 13 shows an efficient example with junctions connecting to multiple sheets. In the voxel space, it is hard to represent such complex junctions exactly, so they are represented by a set of simple junctions as an approximation (Fig. 13(b) and (c)). Since conventional methods create polygonal meshes from them, the meshing result also includes a set of simple junctions. Whereas, our method removes neighboring junction voxels during the sub-sampling phase using spherical supports. This is because correct junction voxels have larger distance values, and unnecessary ones are included in the spherical supports of correct junction points. Fig. 13(d) shows the outcome, and non-manifold edges are also shown in Fig. 13(e). Note that the center axis of the object is represented by single line sequence. On the other hand, our method failed to polygonize some special cases. First case is handling thick objects. Current implementation imposes medial voxel computation to thick objects (e.g. bolts in the lower row in Fig. 11(c)), although they should be polygonized as closed solid surfaces. Then, the radius of the sub-sampling sphere around solid objects becomes large and many necessary voxels around there are absorbed. As a result, the topology structure is broken. Other issue can be found in overlapping of thin-plate objects at the bottom part of Fig. 11(c). We first binarize input volumes to compute medial voxels and then overlapping region is considered as a single surface and the thickness at point on the region will be large. Our method contributes to the acceleration of digital engineering processes. One possible application is in conversion to parametric surfaces (a.k.a. reverse engineering). Fig. 14 shows the results of converting NURBS surfaces from a part of the side member model (Fig. 11) using a reverse engineering software. In this procedure, we decomposed the polygonal mesh into a set of manifold polygons and applied surface fitting to the manifold polygon. Once CAD models are converted, they can be used for a number of advanced applications such as CAE and quality assurance. 5. Conclusion and future work This paper has proposed a method for polygonization of skeleton structures of volumetric data. The technique applies subsampling on medial voxels with spherical supports depending on digital topology so that junctions and boundary voxels are preferentially sampled. Triangular meshes are computed by connecting neighboring triple vertices. Since sampling points are placed on singular voxels, junction edges are well reconstructed. In addition, this paper also showed our method can be used not only for polygonizing wire–skeleton structures but also for polygonizing different skeleton structures at once. For future work, we consider following two issues. First, our method does not guarantee that junction voxels are always connected, even though all examples in this paper are correctly
(a) Input volume.
827
(b) Medial voxel.
(d) Medial surface.
(c) Intersection of (b).
(e) Junction edges.
Fig. 13. A result for an artificial volume (256 × 256 × 256) with high-valence junctions. Note that the center axis can represent only single edge sequence in (d) and (e).
Fig. 14. CAD surface reconstruction from the bottom part of side member model (Fig. 11). Black curves, green points are boundary curves of each segments and control points of fitted surfaces respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
polygonized. For instance, sub-sampled voxels may propagate over pieces due to sparse density at junction voxels. We would also like to combine this method with solid surface reconstruction. We believe that this sub-sampling-based approach provides a unified framework for polygonization of different and mixed dimensional objects. Acknowledgments The authors would like to thank Koichi Matsuzaki for converting polygonal meshes to CAD data using his reverse engineering software. Door panel data (Fig. 11) and Knot data (Fig. 12) are courtesy of Hitachi and Robert G. Scharrein respectively. This research was partially supported by the Ministry of Education, Science, Sports and Culture, Grant-in-Aid for Scientific Research (A), (No. 22246018) and Grant-in-Aid for Young Scientists(B) (No. 22700091). Appendix. Topology classification of voxels The topology of each voxel can be estimated by counting the number of foreground and background components in neighboring voxels [26]. In this method, the following two entities are used:
• C ∗ : The number of 26-connected foreground voxel components in 26-neighbor voxels of v that are 26-adjacent to v . • C¯ : The number of 6-connected background voxel components in 18-neighbor voxels of v that are 6-adjacent to v with the exception of v itself. Using these entities, junction, boundary and surface voxels can be classified as follows:
• Junction voxel: C ∗ = 1C¯ > 2.
828
T. Michikawa, H. Suzuki / Computer-Aided Design 45 (2013) 822–828
(a) Junction.
(b) Boundary.
(c) Surface.
Fig. A.15. Examples of topology classification. Foreground and background voxels are indicated in black and white, respectively.
• Boundary voxel: C ∗ = 1C¯ = 1. • Surface voxel: C ∗ = 1C¯ = 2. Fig. A.15 shows examples of classification. References [1] Lorensen WE, Cline HE. Marching cubes: a high resolution 3D surface construction algorithm. In: SIGGRAPH’87: proceedings of the 14th annual conference on computer graphics and interactive techniques. New York (NY, USA): ACM; 1987. p. 163–9. [2] Ju T, Losasso F, Schaefer S, Warren J. Dual contouring of hermite data. In: SIGGRAPH’02: proceedings of the 29th annual conference on computer graphics and interactive techniques. New York (NY, USA): ACM; 2002. p. 339–46. [3] Fujimori T, Suzuki H, Kobayashi Y, Kase K. Contouring medial surface of thin plate structure using local marching cubes. In: International conference on shape modeling and applications. 2004. p. 297–306. [4] Prohaska S, Hege H-C. Fast visualization of plane-like structures in voxel data. In: Proceedings of the conference on visualization’02. Washington (DC, USA): IEEE Computer Society; 2002. p. 29–36. [5] Suzuki H, Fujimori T, Michikawa T, Miwata Y, Sadaoka N. Skeleton surface generation from volumetric models of thin plate structures for industrial applications. In: IMA conference on the mathematics of surfaces. LNCS, vol. 4647. 2007. p. 442–64. [6] Liu L, Chambers EW, Letscher D, Ju T. A simple and robust thinning algorithm on cell complexes. Computer Graphics Forum 2010;29(7):2253–60. [7] Ohtake Y, Belyaev A, Seidel H-P. An integrating approach to meshing scattered point data. In: Proceedings of the 2005 ACM symposium on solid and physical modeling. New York (NY, USA): ACM; 2005. p. 61–9. [8] Blum H. A transformation for extracting new descriptors of shape. In: WathenDunn W, editor. Models for the perception of speech and visual form. Cambridge: MIT Press; 1967. p. 362–80. [9] Palágyi K, Kuba A. A parallel 12-subiteration 3D thinning algorithm to extract medial lines. In: Proceedings of the 7th international conference on computer analysis of images and patterns. London (UK): Springer-Verlag; 1997. p. 400–7. [10] Toriwaki J, Mori K. Distance transformation and skeletonization of 3D pictures and their applications to medical images. In: Digital and image geometry. LNCS, vol. 2243. 2001. p. 412–28. [11] Ju T, Baker ML, Chiu W. Computing a family of skeletons of volumetric models for shape description. Computer Aided Design 2007;39(5):352–60.
[12] Tsao YF, Fu KS. A parallel thinning algorithm for 3-D pictures. Computer Graphics and Image Processing 1981;17(4):315–31. [13] Manzanera A, Bernard T, Preteux F, Longuet B. Medial faces from a concise 3D thinning algorithm. In: IEEE international conference on computer vision. Vol. 1. 1999. p. 337. [14] Borgefors G, Nystrom I, Baja GSD. Computing skeletons in three dimensions. Pattern Recognition 1999;32(7):1225–36. [15] Michikawa T, Nakazaki S, Suzuki H. Efficiend medial voxel extraction from large volumetric models. In: Proceeding of WSCG 2009. 2009. p. 169–76. [16] Dey TK, Li K, Ramos EA, Wenger R. Isotopic reconstruction of surfaces with boundaries. Computer Graphics Forum 2009;28(5):1371–82. [17] Attali D, Montanvert A. Computing and simplifying 2D and 3D continuous skeletons. Computer Vision and Image Understanding 1997;67(3): 261–273. [18] Etzion M, Rappoport A. Computing the Voronoi diagram of a 3-D polyhedron by separate computation of its symbolic and geometric parts. In: Proceedings of the fifth ACM symposium on solid modeling and applications. New York (NY, USA): ACM; 1999. p. 167–78. [19] Amenta N, Choi S, Kolluri RK. The power crust. In: SMA’01: proceedings of the sixth ACM symposium on solid modeling and applications. New York (NY, USA): ACM; 2001. p. 249–66. [20] Dey TK, Zhao W. Approximate medial axis as a Voronoi subcomplex. In: Proceedings of the seventh ACM symposium on solid modeling and applications. New York (NY, USA): ACM; 2002. p. 356–66. [21] Foskey M, Lin MC, Manocha D. Efficient computation of a simplified medial axis. In: Proceedings of the eighth ACM symposium on solid modeling and applications. New York (NY, USA): ACM; 2003. p. 96–107. [22] Sud A, Foskey M, Manocha D. Homotopy-preserving medial axis simplification. In: Proceedings of the 2005 ACM symposium on solid and physical modeling. New York (NY, USA): ACM; 2005. p. 39–50. [23] Giesen J, Miklos B, Pauly M, Wormser C. The scale axis transform. In: Proceedings of the 25th annual symposium on computational geometry. New York (NY, USA): ACM; 2009. p. 106–15. [24] Miklos B, Giesen J, Pauly M. Discrete scale axis representations for 3D geometry. ACM Transactions on Graphics 2010;29: 101:1–101:10. [25] Shimada K, Gossard DC. Bubble mesh: automated triangular meshing of nonmanifold geometry by sphere packing. In: Proceedings of the third ACM symposium on solid modeling and applications. New York (NY, USA): ACM; 1995. p. 409–19. [26] Malandain G, Bertrand G, Ayache N. Topological segmentation of discrete surfaces. International Journal of Computer Vision 1993;10(2):183–97. [27] Masuda H. Topological operators and Boolean operations for complex-based nonmanifold geometric models. Computer-Aided Design 1993;25(2):119–29.
Contents lists available at SciVerse ScienceDirect
Computer-Aided Design journal homepage: www.elsevier.com/locate/cad
Polygonization of volumetric skeletons with junctions Takashi Michikawa ∗ , Hiromasa Suzuki Research Center for Advanced Science and Technology, The University of Tokyo, Japan
article
info
Keywords: Volumetric data Medial axis Non-manifold meshes Reverse engineering
abstract This article presents a method for polygonizing skeleton structures computed from volumetric data, mainly of a thin-plate object. The main idea of our method is to simplify the topological structure of volumetric skeletons by sub-sampling. Given volumetric skeletons and distance fields of binarized volumes, we first apply sub-sampling of the skeleton voxels so that the junction and boundary voxels are preferentially sampled. By using sampled points as Voronoi sites, we compute Voronoi diagrams on the skeleton voxels. Polygonization can be computed by using connectivity of the diagram. This approach can build correct junctions where the skeleton branches in several directions, whereas conventional voxelbased methods tend to result in small cavities around them. This paper also demonstrates several results computed from CT-scanned engineering objects. © 2011 Elsevier Ltd. All rights reserved.
1. Introduction This paper outlines a method for converting volumetric skeleton structures or medial axes of thin-plate objects to polygonal meshes for use in digital engineering applications. Thin-plate engineering objects are made of metal plates, and are formed by stamping and welding; they are often found in cars and home electronics. Industrial companies measure such objects as volumetric images using X-ray CT scanners or nondestructive scanning devices (Fig. 1(a)). Since each voxel has a CT value that identifies its material, we can extract volumetric data with certain threshold values from CT images. Companies need polygonal meshes of these CT images for use in accelerating engineering processes such as comparison with CAD models and CAE simulation. It is well known that polygonal meshes can be obtained by isosurface contouring of volumetric data [1,2]. This is effective for solid objects, but can create closed surfaces with thin-plate objects (Fig. 1(b)). The center-line structure of the objects shown in Fig. 1(c) are preferable because common CAD systems represent thin-plate objects as open surfaces [3]. Thin-plate objects are usually welded to other objects for reinforcement. The medial surfaces of such objects become nonmanifold, and their welded parts become junctions. Since our focus here is on engineering application, these forms must be reconstructed correctly. One of the challenges involved is capturing non-manifold junctions from CT images correctly as shown in Fig. 1(c).
∗
Corresponding author. Tel.: +81 3 5452 5183. E-mail addresses: [email protected] (T. Michikawa), [email protected] (H. Suzuki). 0010-4485/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.cad.2011.06.003
The computation of medial surfaces from volumetric data has been extensively discussed in previous studies [4,3,5]. The basic approach first creates medial voxels from input volumetric data, and polygonal meshes are then created based on the resulting voxel connectivity. However, this approach creates small cavities around junctions even though no such cavities exist in the input data. This is due to the ambiguity of triangulation from voxel connectivity. Prohaska and Hege [4] introduced a look-up table for non-manifold topology polygonization. However, local voxels are still checked due to the existence of bad cases of cavities. Liu et al. proposed a robust algorithm by pruning the input models as cell complexes [6]. Their method can create topologically correct surfaces, however their polygonization are axis-aligned and it is difficult to represent higher-valence junctions. We present a method of medial surface reconstruction from volumetric data of thin-plate objects including junctions. Given medial voxels with their distance to the boundary, we first apply sub-sampling to these voxels so that junction and boundary voxels can be preferentially selected. This computation is performed by covering the medial voxels with spherical supports [7]. We then perform Voronoi diagram computation for medial voxels using the sampled points as Voronoi sites to build vertex connectivity, from which triangular meshes are then created. This approach offers various advantages in medial surface reconstruction from volumetric data. In particular, spherical supports in the sub-sampling phase remove unnecessary voxels around junctions. This contributes to the absorption of various types of artifacts arising in the junctions regions, including branches and small cavities. In addition, our method can be used for computing not only medial surfaces but also 1D skeletons and their coupling is also possible. This paper consists of five sections. We review related work in Section 2, the proposed algorithm is introduced in Section 3, we
T. Michikawa, H. Suzuki / Computer-Aided Design 45 (2013) 822–828
823
polygonization of non-manifold topology. However, local voxels are still checked due to the existence of bad cases of cavities. Liu et al. recently proposed a method for computing medial structures (wires and skeletons) by pruning the input models as cell complexes [6]. However, their polygonization depends on the topology of cell complex and it is usually insufficient for representing higher-valence junctions.
Fig. 1. CT images of welded aluminum and its polygonal meshes. The resolution of volumetric data (a) is 350 × 330 × 220 and its isosurface (b) and medial surface have 658,648 and 7539 faces (Computation Time: about 1 min.) respectively. Note that the junction edges (red) of the medial surfaces created using our method (c) are correctly reconstructed. Lower row shows corresponding illustrations. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
outline the experimental results and discuss them in Section 4, and Section 5 concludes the paper. 2. Related work Our issue is related to a medial axis whose formulation was first introduced in [8]; that is, the set of center points of a collection of maximal inscribed spheres. Many approaches involving MATs (medial axis transforms) have been studied over a period of four decades in the image processing and geometric modeling community. This section reviews several approaches related to our objective. Volume-based methods. Volume-based approaches compute medial voxels from binarized volumetric data. This produces a volumetric version of the medial axis. A typical strategy is topological thinning, or voxel pruning, which leaves the topology of the input data unchanged (such voxels are called ‘‘simple voxels’’) and preserves the ends of surface voxels. This process is iteratively applied until no more voxels can be removed. Sequential thinning removes simple voxels one by one. This approach preserves the topology of the original data, but the order of pruning may preserve excess simple voxels and cause many branches or bumps. Although several solutions to this have been discussed (e.g., sub-cycles [9], distance fields [10] and relaxed conditions [11]), the issue has not yet been completely resolved. Parallel thinning removes all simple voxels on the boundary of an object simultaneously. Tsao and Fu [12] formulated conditions for building smooth surfaces, but the topology of the input model may be lost with this technique. The resulting surfaces are relatively smooth, but it is difficult to preserve topology because many voxels are removed at once. There are also topologypreserving methods for parallel thinning [13,14], but these cause a deterioration in voxel quality. Prohaska and Hege [4] used the geodesic distance between two closest boundary voxels to determine pruning priority. This global metric is more stable for noisy data, and hardly any unnecessary branches are generated. Accordingly, it is used in a number of engineering applications [3,5,15]. Polygonal meshes can be computed from medial voxels. For instance, we can cast this issue to surface reconstruction problems (e.g., [16]). However, non-manifold surfaces cannot be handled. Another approach builds triangular meshes using voxel connectivity, but this method produces ambiguous cases of triangulation from such connectivity, and small cavities may be created. Prohaska and Hege [4] introduced a look-up table for the
Polygon-based methods. Alternative approaches involve computing medial axes from polygonal meshes. Although an exact solution can be computed by solving algebraic equations, this is not practical for large models, and most polygon-based methods deal with approximate solutions of MAT. One popular approach is the use of Voronoi diagrams [17–20]. The structures of medial axes computed using the above methods are usually complex, and may have small branch surfaces. For simplified medial axes, these are usually removed using certain criteria. One popular criterion is the separation angle of sheets [17,21]. Dey and Zhao [20] used additional conditions to improve enhancement of the quality of medial surfaces. Sud et al. [22] also improved the medial axis simplification method [21] so that the homotopy of the input data is preserved. Recently, the scale axis transform [23,24] algorithm was introduced to simplify medial surfaces by scaling the medial ball of the objects. However, they depend on definition of the medial axis. For instance, it is difficult to represent a higher-valence junctions because they may be represented as a set of lowervalence junctions. When these techniques are applied to our issue, polygonal meshes must be computed using isosurface extraction methods. This is straightforward, but the resulting meshes have several problems including numerical issues for complex models and artifacts caused by noise and blur in CT images. Accordingly, polygon-based approaches are not practical in achieving our objective. 3. Medial surface reconstruction from medial voxels The goal of this algorithm is to polygonize a skeleton structure with junctions of volumetric data. Fig. 2 shows an overview of our method. In this article, we suppose that medial voxels (from Fig. 2(a)–(d)) are already computed by conventional methods (e.g. [4,3,5,15]) (Fig. 3). We first compute binary volumes generated by thresholding volumetric data. The threshold value is a user-defined parameter. We next compute distance fields D = {di } from the binarized volumes (Fig. 3(c)). Once distance fields are computed, medial voxels M = {vi } are computed by evaluating the geodesic distance between the two closest boundary points neighboring to each voxel. Fig. 3(d) shows medial voxels computed from CT images in Fig. 3(a). The distance value at the medial voxels are drawn in color map in the figure. A method for converting medial voxels to polygons consists of three steps: Sub-sampling—distributes spherical supports on voxels iteratively so that they completely cover the input medial voxels (Fig. 2(e)–(g)). Voronoi diagram—computes a Voronoi diagram on voxels to build mesh edges (Fig. 2(h)). Non-manifold meshing—computes polygonal meshes by finding three vertices connected to each other and applies topological cleaning to remove excess non-manifold edges (Fig. 2(i)).
824
T. Michikawa, H. Suzuki / Computer-Aided Design 45 (2013) 822–828
a
b
c
a
b Sub-sampling
d
e
f
g
h
i
Fig. 4. A 2D comparison of meshing results with random distribution (a) and topology-dependent distribution (b). In each figure, the voxels shown in red denote sampled points, and dotted circles represent the spherical supports of voxels. The red lines show medial surfaces connected by a Voronoi diagram on medial voxels. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
(a) #component = 1. Fig. 2. An overview of our method. (a) CT images. (b) Binarized volumes (FG and BG voxels are drawn in black and white respectively). (c) Distance field of (b). (d) Medial voxels. Each voxel has distance to the closest boundary voxel of the object. (e)–(g) Sub-sampling of medial voxels. Neighboring voxels (gray) are removed. (h) Voronoi diagram on medial voxels using sampled points as Voronoi sites. (i) Polygonization.
a
b
c
d
Fig. 3. Medial voxel extraction from binary volumes. (a) Input CT images. (b) Binary volume. (c) Distance field of (b) (Cut-away view). (d) Medial voxels. Distance value is assigned to each voxel. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
3.1. Sub-sampling of medial voxels The first step of the algorithm involves medial voxel subsampling. Given a set of medial voxels M, we select one voxel vi from M and set it to vertex list V . Then, we remove neighboring voxels vj such that ∥vi − vj ∥ < α di from M, where vi denotes the position of vi , and α denotes the radius scaling factor. Note that a larger value for α usually makes robust triangulation at very thin parts of the object because the radii of such parts are usually small and may cause bad sampling of medial voxels. On the other hand, too large α also remove many voxels neighboring to vi and it causes ill shaped triangulation. Accordingly, a value of α = 2 usually makes good results. Iteration is performed until M is empty and then all the voxels are covered by spherical supports.
(b) #component = 2.
(c) #component = 3.
Fig. 5. End point classification. We count the number of components in the 5 × 5 × 5 masks with the exception to 3 × 3 × 3 mask (shaded region). If the number of components is 2, the voxel is not end point (b). Otherwise, the point is end point ((a) and (c)).
An important point here is how to find these vi voxels from M, because their distribution affects mesh topology directly. Fig. 4 illustrates a 2D example. Random distribution may not involve placement at junctions, and sampling points around junctions may be connected to each other, thereby creating small cavities as shown in Fig. 4(a). On the other hand, in Fig. 4(b), spheres are distributed at junction and boundary voxels first, and correct polygonization is obtained. We therefore introduce the topologydependent distribution as described below. Topology criterion control for sampling is implemented so that points are well distributed on junctions and boundary edges. This contributes to fastening low-dimensional topology features such as junctions and boundaries. We use it here as the primary priority. A similar idea can be found in the bubble mesh algorithm [25] for polygonal meshes with explicit topology, and we apply similar approach to volumetric data. The basic classification algorithm was studied in [26] (see Appendix for a brief summary). In addition, we also consider end points of junctions or boundary, because such points must be sampled first. After the classification, we extract wire structures consisting of boundary and junction voxels. We then count the number of components of the wire structure within 5 × 5 × 5 mask with the exception to 3 × 3 × 3 mask (Fig. 5). If the number of components is 2, the point is a normal boundary (Fig. 5(b)), otherwise it is an end point (Fig. 5(a) and (c)). We also consider the distance values of voxels as a secondary criterion. This is because medial voxels have relatively larger distances than neighboring voxels, and it is natural that such voxels should be selected first. Finally, the topology-dependent sub-sampling algorithm is summarized as below. 1. Find the vi with the top priority in M and copy it to V . 2. Remove vj from M such that ∥vi − vj ∥ < α di . 3. Iterate 1 and 2 until M becomes empty. Fig. 6 shows a result of topology classification of the medial voxels (Fig. 3(d)) and the sub-sampling result (Fig. 6(b)). Color of each sphere means a sequence of sub-sampling (red sphere is sampled first). We can see medial voxels on non-manifold features are preferentially selected.
T. Michikawa, H. Suzuki / Computer-Aided Design 45 (2013) 822–828
(a) Topology classification.
(b) Sampling result.
Fig. 6. Topology-dependent sub-sampling. (a) Topological classification of medial voxels. Surface, boundary, junction and end point voxels are drawn in gray, blue, red and green respectively. (b) Sub-sampling result. A center point of each sphere is sampled. Colors of each sphere shows sequence of sub-sampling (red spheres are selected first and blue spheres are selected later.) (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
825
Fig. 9. 2D illustration of neighboring voxels in surface voxels. Note that the planer graph can be defined between F (v) and B0 (v) (indicated by the bold red line). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
(a) Voronoi diagram.
(b) Polygonization.
Fig. 7. Voronoi diagram definition on medial voxels. Identically colored voxels belong to the same sites.
(a) Non-manifold surface.
(b) Manifold decomposition.
Fig. 8. Medial voxels can be separated into three simple sheets consisting of boundary and surface voxels.
3.2. Building a Voronoi diagram on medial voxels We next build mesh edges from the set of sampled vertices V . In our target objects of thin-plate structures, we can assume that medial voxels can be decomposed into a collection of subsets, each of which is an analog of a two-manifold surface (See below for details). This subset therefore inherits the properties of a twomanifold surface. In particular, it allows us to define Voronoi diagrams on this subset. Fig. 7 illustrates the formulation of Voronoi diagrams on medial voxels. Given sites V = {vi } ⊆ M, we define the region R(vi ) as follows: R(vi ) = x|d(vi , x) ≤ d(vj , x), ∀j ̸= i ,
(1)
where d(vi , x) denotes the shortest path length between vi and x in the voxel space. In practice, we compute it using the wavefront propagation of distance from the sites. Fig. 10(a) shows the result of the diagram for medial voxels M. We define an edge (vi , vj ) if there is a neighboring voxel pair (x, y) such that x ∈ R(vi ) and y ∈ R(vj ) in 26-adjacent. Validity of polygonization. Medial voxels in our method can be divided into surface, boundary and junction types by topology classification (Fig. 8(a)). They can be separated into a set of sheets consisting of boundary and surface voxels by junction voxels
(c) Polygonization result.
(d) Mesh cleaning result.
Fig. 10. Polygonization with a Voronoi diagram and mesh cleaning.
(Fig. 8(b)). We can also consider that the Voronoi diagram is defined on a surface composed of the faces of foreground voxels shared by background voxels. It should be noted here that a voxel is regarded as a cubic volume. Based on the definition of topological classification, each surface voxel v in a sheet separates its 26-neighboring voxels N (v) into two background components Bi (v) (i = 0, 1) using one foreground component F (v) (Fig. 9). When we focus on B0 (v) and F (v), the boundary surface between their voxels is topologically equivalent to a disk, because B0 (v) is 6-connected and its surface must be two-manifold. This indicates that every surface piece at a voxel on the sheet is topologically equivalent to a disk and thus the total surface of the sheet is topologically equivalent to a twomanifold surface. The Voronoi diagram is therefore defined on surface voxels and we can generate a polygonal mesh by taking the dual of the Voronoi diagram.
826
T. Michikawa, H. Suzuki / Computer-Aided Design 45 (2013) 822–828
(a) Input.
(b) Sub-sampling.
(c) Result.
Fig. 11. A result for side member model (Upper, 400 × 400 × 576 voxels, 7500 faces, 2.79 min.) and door panel model (Lower, 1370 × 880 × 500 voxels, 21,748 faces, 21.71 min.)
(a) Input.
(b) Sub-sampling.
(c) Result.
Fig. 12. A result for Knot model (Upper, Wire only, 512 × 512 × 512 voxel, 254 wire edges, 62.80 s) and coupled model (Lower, Wire and surfaces, 512 × 512 × 512 voxels, 2267 faces, 38 wire edges, 79.50 s.)
3.3. Non-manifold meshing from vertex connectivity
3.4. An extension to wire–skeleton structure extraction
In the final step, polygonal meshes are constructed from the vertices and edges computed in the previous steps (Fig. 10(a)). This computation is performed by finding triples (vi , vj , vk ) that are connected to each other (Fig. 10(b)). This approach creates non-manifold meshes with no small cavities caused by invalid sampling, but may also create overlapping triangles on surface voxels (Fig. 10(c)) for similar reasons related to the degeneration of Delaunay triangulation. Such triangles must be eliminated so that necessary ones around junctions are not removed. Based on medial voxel topology classification, we decompose all vertices into manifold vertices (on surfaces) and non-manifold vertices (on boundaries or junctions), and then apply manifold cleaning [7] to the manifold vertices only. Simply speaking, this method performs iteration to find a planar graph from the one-ring neighborhood of each vertex such that the number of vertices in the graph is maximized. Triangles excluded from the planar graph are removed from the mesh. Fig. 10(d) shows the results of cleaning. Note that some overlapping triangles may still be left if all vertices belong to junctions. Such issues are generally seen when junctions are very close.
An extension of our method to 1D skeleton extraction is straightforward. Moreover, the 1D skeleton or the wire structure can also be computed at once. Fig. 12 show results for a knot model and a volumetric data involving surface and wire structures. We can see both of results have correct topology. The procedure is almost same as the original described in this section. By Voronoi diagrams on medial voxels, we already have connectivity of wire edges. We just need to manage these edges so that they cannot be removed by surface cleaning. 4. Result and discussion The experimental results of meshing are shown in Figs. 11 and 12. Each figure shows input binary volumes, the sub-sampling result and the final meshing result. Figs. 1 and 13 are also obtained by our method. All examples were tested on a 3.16 GHz computer, and the computation times are provided in the caption of each figure. There is still room for improvement here, as our prototype is not yet optimized. The current bottleneck lies in Voronoi diagram creation and mesh cleaning. Since Voronoi diagrams compute using a multi-source Dijkstra’s algorithm, the computation cost is
T. Michikawa, H. Suzuki / Computer-Aided Design 45 (2013) 822–828
O(n log n). Note that although connecting spherical supports [7] may create polygons fast, the technique truncates topology information hidden in medial voxels, whereas the Voronoi-based method considers such topological information. The non-manifold mesh data structure used in cleaning is based on an IndexedFaceSet (a.k.a. Wavefront OBJ format), and takes a long time to search neighboring vertices. An advanced data structure for non-manifold meshes (e.g., [27]) is expected to improve the level of performance. Handling large objects is an important issue because CT images of engineering objects are usually large. Since our method manages medial voxels by octree, our method is efficient for memory usage and can handle large objects as shown in the lower row in Fig. 11. Using topology-dependent sampling, medial surfaces with junctions were reconstructed from CT-scanned engineering objects. In addition, sub-sampling with spherical supports contributes to the recovery of non-manifold junctions. Fig. 13 shows an efficient example with junctions connecting to multiple sheets. In the voxel space, it is hard to represent such complex junctions exactly, so they are represented by a set of simple junctions as an approximation (Fig. 13(b) and (c)). Since conventional methods create polygonal meshes from them, the meshing result also includes a set of simple junctions. Whereas, our method removes neighboring junction voxels during the sub-sampling phase using spherical supports. This is because correct junction voxels have larger distance values, and unnecessary ones are included in the spherical supports of correct junction points. Fig. 13(d) shows the outcome, and non-manifold edges are also shown in Fig. 13(e). Note that the center axis of the object is represented by single line sequence. On the other hand, our method failed to polygonize some special cases. First case is handling thick objects. Current implementation imposes medial voxel computation to thick objects (e.g. bolts in the lower row in Fig. 11(c)), although they should be polygonized as closed solid surfaces. Then, the radius of the sub-sampling sphere around solid objects becomes large and many necessary voxels around there are absorbed. As a result, the topology structure is broken. Other issue can be found in overlapping of thin-plate objects at the bottom part of Fig. 11(c). We first binarize input volumes to compute medial voxels and then overlapping region is considered as a single surface and the thickness at point on the region will be large. Our method contributes to the acceleration of digital engineering processes. One possible application is in conversion to parametric surfaces (a.k.a. reverse engineering). Fig. 14 shows the results of converting NURBS surfaces from a part of the side member model (Fig. 11) using a reverse engineering software. In this procedure, we decomposed the polygonal mesh into a set of manifold polygons and applied surface fitting to the manifold polygon. Once CAD models are converted, they can be used for a number of advanced applications such as CAE and quality assurance. 5. Conclusion and future work This paper has proposed a method for polygonization of skeleton structures of volumetric data. The technique applies subsampling on medial voxels with spherical supports depending on digital topology so that junctions and boundary voxels are preferentially sampled. Triangular meshes are computed by connecting neighboring triple vertices. Since sampling points are placed on singular voxels, junction edges are well reconstructed. In addition, this paper also showed our method can be used not only for polygonizing wire–skeleton structures but also for polygonizing different skeleton structures at once. For future work, we consider following two issues. First, our method does not guarantee that junction voxels are always connected, even though all examples in this paper are correctly
(a) Input volume.
827
(b) Medial voxel.
(d) Medial surface.
(c) Intersection of (b).
(e) Junction edges.
Fig. 13. A result for an artificial volume (256 × 256 × 256) with high-valence junctions. Note that the center axis can represent only single edge sequence in (d) and (e).
Fig. 14. CAD surface reconstruction from the bottom part of side member model (Fig. 11). Black curves, green points are boundary curves of each segments and control points of fitted surfaces respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
polygonized. For instance, sub-sampled voxels may propagate over pieces due to sparse density at junction voxels. We would also like to combine this method with solid surface reconstruction. We believe that this sub-sampling-based approach provides a unified framework for polygonization of different and mixed dimensional objects. Acknowledgments The authors would like to thank Koichi Matsuzaki for converting polygonal meshes to CAD data using his reverse engineering software. Door panel data (Fig. 11) and Knot data (Fig. 12) are courtesy of Hitachi and Robert G. Scharrein respectively. This research was partially supported by the Ministry of Education, Science, Sports and Culture, Grant-in-Aid for Scientific Research (A), (No. 22246018) and Grant-in-Aid for Young Scientists(B) (No. 22700091). Appendix. Topology classification of voxels The topology of each voxel can be estimated by counting the number of foreground and background components in neighboring voxels [26]. In this method, the following two entities are used:
• C ∗ : The number of 26-connected foreground voxel components in 26-neighbor voxels of v that are 26-adjacent to v . • C¯ : The number of 6-connected background voxel components in 18-neighbor voxels of v that are 6-adjacent to v with the exception of v itself. Using these entities, junction, boundary and surface voxels can be classified as follows:
• Junction voxel: C ∗ = 1C¯ > 2.
828
T. Michikawa, H. Suzuki / Computer-Aided Design 45 (2013) 822–828
(a) Junction.
(b) Boundary.
(c) Surface.
Fig. A.15. Examples of topology classification. Foreground and background voxels are indicated in black and white, respectively.
• Boundary voxel: C ∗ = 1C¯ = 1. • Surface voxel: C ∗ = 1C¯ = 2. Fig. A.15 shows examples of classification. References [1] Lorensen WE, Cline HE. Marching cubes: a high resolution 3D surface construction algorithm. In: SIGGRAPH’87: proceedings of the 14th annual conference on computer graphics and interactive techniques. New York (NY, USA): ACM; 1987. p. 163–9. [2] Ju T, Losasso F, Schaefer S, Warren J. Dual contouring of hermite data. In: SIGGRAPH’02: proceedings of the 29th annual conference on computer graphics and interactive techniques. New York (NY, USA): ACM; 2002. p. 339–46. [3] Fujimori T, Suzuki H, Kobayashi Y, Kase K. Contouring medial surface of thin plate structure using local marching cubes. In: International conference on shape modeling and applications. 2004. p. 297–306. [4] Prohaska S, Hege H-C. Fast visualization of plane-like structures in voxel data. In: Proceedings of the conference on visualization’02. Washington (DC, USA): IEEE Computer Society; 2002. p. 29–36. [5] Suzuki H, Fujimori T, Michikawa T, Miwata Y, Sadaoka N. Skeleton surface generation from volumetric models of thin plate structures for industrial applications. In: IMA conference on the mathematics of surfaces. LNCS, vol. 4647. 2007. p. 442–64. [6] Liu L, Chambers EW, Letscher D, Ju T. A simple and robust thinning algorithm on cell complexes. Computer Graphics Forum 2010;29(7):2253–60. [7] Ohtake Y, Belyaev A, Seidel H-P. An integrating approach to meshing scattered point data. In: Proceedings of the 2005 ACM symposium on solid and physical modeling. New York (NY, USA): ACM; 2005. p. 61–9. [8] Blum H. A transformation for extracting new descriptors of shape. In: WathenDunn W, editor. Models for the perception of speech and visual form. Cambridge: MIT Press; 1967. p. 362–80. [9] Palágyi K, Kuba A. A parallel 12-subiteration 3D thinning algorithm to extract medial lines. In: Proceedings of the 7th international conference on computer analysis of images and patterns. London (UK): Springer-Verlag; 1997. p. 400–7. [10] Toriwaki J, Mori K. Distance transformation and skeletonization of 3D pictures and their applications to medical images. In: Digital and image geometry. LNCS, vol. 2243. 2001. p. 412–28. [11] Ju T, Baker ML, Chiu W. Computing a family of skeletons of volumetric models for shape description. Computer Aided Design 2007;39(5):352–60.
[12] Tsao YF, Fu KS. A parallel thinning algorithm for 3-D pictures. Computer Graphics and Image Processing 1981;17(4):315–31. [13] Manzanera A, Bernard T, Preteux F, Longuet B. Medial faces from a concise 3D thinning algorithm. In: IEEE international conference on computer vision. Vol. 1. 1999. p. 337. [14] Borgefors G, Nystrom I, Baja GSD. Computing skeletons in three dimensions. Pattern Recognition 1999;32(7):1225–36. [15] Michikawa T, Nakazaki S, Suzuki H. Efficiend medial voxel extraction from large volumetric models. In: Proceeding of WSCG 2009. 2009. p. 169–76. [16] Dey TK, Li K, Ramos EA, Wenger R. Isotopic reconstruction of surfaces with boundaries. Computer Graphics Forum 2009;28(5):1371–82. [17] Attali D, Montanvert A. Computing and simplifying 2D and 3D continuous skeletons. Computer Vision and Image Understanding 1997;67(3): 261–273. [18] Etzion M, Rappoport A. Computing the Voronoi diagram of a 3-D polyhedron by separate computation of its symbolic and geometric parts. In: Proceedings of the fifth ACM symposium on solid modeling and applications. New York (NY, USA): ACM; 1999. p. 167–78. [19] Amenta N, Choi S, Kolluri RK. The power crust. In: SMA’01: proceedings of the sixth ACM symposium on solid modeling and applications. New York (NY, USA): ACM; 2001. p. 249–66. [20] Dey TK, Zhao W. Approximate medial axis as a Voronoi subcomplex. In: Proceedings of the seventh ACM symposium on solid modeling and applications. New York (NY, USA): ACM; 2002. p. 356–66. [21] Foskey M, Lin MC, Manocha D. Efficient computation of a simplified medial axis. In: Proceedings of the eighth ACM symposium on solid modeling and applications. New York (NY, USA): ACM; 2003. p. 96–107. [22] Sud A, Foskey M, Manocha D. Homotopy-preserving medial axis simplification. In: Proceedings of the 2005 ACM symposium on solid and physical modeling. New York (NY, USA): ACM; 2005. p. 39–50. [23] Giesen J, Miklos B, Pauly M, Wormser C. The scale axis transform. In: Proceedings of the 25th annual symposium on computational geometry. New York (NY, USA): ACM; 2009. p. 106–15. [24] Miklos B, Giesen J, Pauly M. Discrete scale axis representations for 3D geometry. ACM Transactions on Graphics 2010;29: 101:1–101:10. [25] Shimada K, Gossard DC. Bubble mesh: automated triangular meshing of nonmanifold geometry by sphere packing. In: Proceedings of the third ACM symposium on solid modeling and applications. New York (NY, USA): ACM; 1995. p. 409–19. [26] Malandain G, Bertrand G, Ayache N. Topological segmentation of discrete surfaces. International Journal of Computer Vision 1993;10(2):183–97. [27] Masuda H. Topological operators and Boolean operations for complex-based nonmanifold geometric models. Computer-Aided Design 1993;25(2):119–29.