Pattern Recognition Letters 23 (2002) 523–531 www.elsevier.com/locate/patrec
A three-dimensional holes closing algorithm Zouina Aktouf, Gilles Bertrand, Laurent Perroton
*
ESIEE Cit e Descartes, A2SI Laboratory, ESIEE 2, Bd Blaise Pascal – B.P 99, Noisy-Le-Grand Cedex 93162, France Received 11 July 2000; received in revised form 3 July 2001
Abstract Contrary to the 2D case, a 3D hole is not a subset of the 3D space. It is therefore not possible to use connected components search algorithms for detecting and suppressing 3D holes. In this paper, we propose an algorithm for suppressing 3D holes. It is based on properties of the previously introduced notion of topological number. Our algorithm is linear in time and it allows to control the ‘‘size’’ of the holes which are closed. As far as we know, this is the first 3D holes closing algorithm. 2002 Elsevier Science B.V. All rights reserved. Keywords: Digital topology; Topological numbers; Topological hull; Holes; Distance transform
1. Introduction In 3D biomedical images, the reconstructed objects which result from the segmentation process have sometimes unwanted holes. These holes can be considered as noise, and may not be avoided by segmentation. However, not all of the holes are considered as noise. In many applications, the noisy holes can be characterized by their size: larger holes are actually features of the object, whereas a large amount of small holes are irrelevant. In this paper, we present an algorithm for suppressing 3D holes which original implementation uses the distance transform to control the holes that are suppressed. We will call this operation closing the holes. Fig. 1 il-
*
Corresponding author. E-mail address:
[email protected] (L. Perroton).
lustrates an example of our algorithm on a reconstructed skull from tomography. Only half of the skull is actually represented so that the hole can be clearly seen. The image has been segmented by a simple threshold which created a big hole in the middle of the shape of the skull which can be seen on the left image. The right image is the result of applying our algorithm to close the hole. In our approach, we consider the notion of holes from a topological point of view. It is important to make the difference between holes, cavities and concavities. Concavities are concave shapes of the contour of the object, as shown in Fig. 2(a). Cavities are hollows in an object, or more formally, bounded connected components of the background. The notion of hole is not very simple to visualize in 3D. For example, a torus has one hole and no cavity, whereas a hollow torus with one cavity has two holes. More formally, the presence of a hole in an object is detected
0167-8655/02/$ - see front matter 2002 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 7 - 8 6 5 5 ( 0 1 ) 0 0 1 5 2 - 0
524
Z. Aktouf et al. / Pattern Recognition Letters 23 (2002) 523–531
Fig. 1. Half of a skull from CT with a hole and after closing the hole.
Fig. 2. (a) A 2D object with two holes. (b) A torus with one hole (detected by one path p). (c) The hole of the torus has been closed.
whenever there is a closed path which cannot be transformed into a single point by a sequence of elementary local deformations inside the object, see Fig. 2(b). The closing of holes in 3D images is much more difficult than in 2D. Indeed, a hole in 2D image is a cavity and thus a region, that is a bounded welldefined set of points, whereas it is not necessarily the case in 3D. Fig. 2(c) illustrates the closing of the hole of a torus. The closing of the hole is actually performed by building a surface such as a kind of patch which obstructs the hole. The article is an extension of the paper (Aktouf et al., 1996) and is organized as follows. In the next section, we will review some basic definitions related to discrete topology of 3D. Then, we will introduce the notion of topological hull
which will be used to define our algorithm from a purely topological point of view. In the following section, we will show how distance transform can be used in the topological algorithm previously introduced to enhance its result by adding a geometrical constraint and obtaining a centered topological hull. It also closes irrelevant holes. The algorithm and the detail of the implementation involving hierarchical lists will be described in Appendix A.
2. Basic topological notions We denote E ¼ Z3 . A point x 2 E is defined by ðx1 ; x2 ; x3 Þ with xi 2 Z. We consider the three neighborhoods:
Z. Aktouf et al. / Pattern Recognition Letters 23 (2002) 523–531
n N26 ðxÞ ¼ x0 2 E; Max½jx1 x01 j; jx2 x02 j; o jx3 x03 j 6 1 ; n N26 ðxÞ ¼ x0 2 E; Max½jx1 x01 j; jx2 x02 j; o jx3 x03 j 6 1 ; n N6 ðxÞ ¼ x0 2 E; jx1 x01 j þ jx2 x02 j o þ jx3 x03 j 6 1 ; n N18 ðxÞ ¼ x0 2 E; jx1 x01 j þ jx2 x02 j o þ jx3 x03 j 6 2 \ N26 ðxÞ: We denote N6 ðxÞ ¼ N6 ðxÞ n fxg, N26 ðxÞ ¼ N26 ðxÞn fxg, N18 ðxÞ ¼ N18 ðxÞ n fxg. Two points x and y are said to be n-adjacent (n ¼ 6, 18, 26) if y 2 Nn ðxÞ. An n-path p is a sequence (possibly empty) of points x0 ; . . . ; xk , with xi n-adjacent to xi1 for i ¼ 1; . . . ; k. If p is non-empty, p is a path from x0 to xk and the length of p is equal to k. If x0 ¼ xk , p is closed. An object X E is said to be n-connected if, for any two points of X, there is an n-path in X between these two points. The equivalence classes relative to this relation are the n-connected components of X. Note that, if X is finite, the infinite connected component of X is called the background, the other connected components of X are called the cavities. The set composed of all the n-connected components of X is denoted Cn ðX Þ. The set of all n-connected components of X n-adjacent to a point x is denoted Cnx ðX Þ. Note that Cn ðX Þ and Cnx ðX Þ are sets of subsets of X and not sets of points. Furthermore, if S is a set, we will denote #S the number of its elements. As in 2D, if we use an n-connectivity for X we have to use another n-connectivity for X , i.e., the 6-connectivity for X is associated to the 18- or the 26-connectivity for X (and vice versa). This is necessary for having a correspondence between the topology of X and the topology of X . Furthermore, it is sometimes necessary to distinguish the 6-connectivity associated with the 18-connectivity and the 6-connectivity associated with the 26-
525
connectivity. In order to handle this distinction, a 6þ -notion will indicate a 6-notion associated with the 18-connectivity. So we can have ðn; nÞ ¼ ð6; 26Þ, ð26; 6Þ, ð6þ ; 18Þ or ð18; 6þ Þ. 2.1. Notion of hole The notion of a hole is not simple to define. The presence of a hole in X is detected whenever there is a closed path in X that cannot be deformed in X to a single point. For example a hollow ball has one cavity and no hole, a solid torus has no cavity and one hole, a hollow torus has one cavity and two holes. For a more formal definition of holes based on digital geometry, the interested reader will refer to Kong (1989). In the remaining of the paper, p, and a; . . . ; z (with possible subscripts) will denote, respectively, paths and points. The concatenation of the path p ¼ xp1 y by the path p0 ¼ yp01 z is the path p p0 ¼ xp1 yp01 z. We give below a more formal definition of deformation which principle is that c0 will be considered as an elementary deformation of c if c and c0 are the same except in a small neighborhood which depends on the connectivity used. Let X E and p 2 X be a point, called the base point. Let c and c0 be two closed m-paths in X. We say that c0 is an elementary n-deformation of c, which is denoted c c0 , if c and c0 are of the form c ¼ pp1 x1 p x2 p2 p, c0 ¼ pp1 x1 p0 x2 p2 p and: • for n ¼ 6, we have m ¼ 6 and x1 p x2 , x1 p0 x2 are included in a unit square (a 2 2 square); • for n ¼ 6þ , 18, 26, we have, respectively, m ¼ 6, 18, 26 and x1 p x2 , x1 p0 x2 are included in a unit cube (a 2 2 2 cube). We say that c0 is an n-deformation of c or c ’ c0 if there is a sequence of closed n-paths c0 ; . . . ; ck such that c ¼ c0 , c0 ¼ ck and ci1 ci for i ¼ 1; . . . ; k. We say that a set X E has no hole if every path included in X may be n-deformed to a single point. 2.2. Topological numbers We introduce now some definitions which allow to extract some topological characteristics of a point.
526
Z. Aktouf et al. / Pattern Recognition Letters 23 (2002) 523–531
Let X E and x 2 E. The geodesic n-neighborhood of x inside X of order k is the set Nnk ðx; X Þ defined recursively by Nn1 ðx; X Þ ¼ Nn ðxÞ \ X and Nnk ðx; X Þ ¼ [fNn ðyÞ\ N26 ðxÞ \ X ; y 2 Nnk1 ðx; X Þg. k In other words Nn ðx; X Þ is the set composed of all points y of N26 ðxÞ \ X such that there exists an npath p from x to y of length less than or equal to k, all points of p, except possibly x, belonging to N26 ðxÞ \ X . The geodesic neighborhoods Gn ðx; X Þ are defined by G6 ðx; X Þ ¼ N62 ðx; X Þ, G6þ ðx; X Þ ¼ 2 N63 ðx; X Þ, G18 ðx; X Þ ¼ N18 ðx; X Þ and G26 ðx; X Þ ¼ 1 N26 ðx; X Þ. We can now give a definition of topological numbers (Bertrand and Malandain, 1994; Bertrand, 1994): Definition 1 (Topological numbers). Let X E and x 2 E. The topological numbers Tn ðx; X Þ are defined by Tn ðx; X Þ ¼ #Cn ½Gn ðx; X Þ . If we use the n-connectivity for X and the nconnectivity for X , Tn ðx; X Þ and Tn ðx; X Þ are considered for describing the topological characteristics of the point x. In particular, the topological numbers allow to detect if a point is simple. Intuitively, a point is simple if its removal does not ‘‘change the topology of the image’’, see Kong (1997) for a precise definition. We have (Bertrand and Malandain, 1994; Bertrand, 1994): Let X E and x 2 X : x is an n-simple point () Tn ðx; X Þ ¼ 1 and Tn ðx; X Þ ¼ 1. The topological numbers allow to detect other types of points. A point such that Tn ðx; X Þ ¼ 0 is an isolated point. If Tn ðx; X Þ ¼ 0, we have an interior point and border points are characterized by Tn ðx; X Þ 6¼ 0. Let us consider a point x such that Tn ðx; X Þ P 2. If we delete x from X, we locally disconnect X. We can say that such a point is a 1D isthmus. In the same manner, a point such that Tn ðx; X Þ P 2 may be called a 2D isthmus since its deletion leads to locally merge connected components of X . Point such that Tn ðx; X Þ ¼ 2 and Tn ðx; X Þ ¼ 1 is a simple 1D isthmus. We also have a simple 2D isthmus if Tn ðx; X Þ ¼ 2 and Tn ðx; X Þ ¼ 1.
3. Topological hull We introduce the notion of topological hull: Let X be a subset of E. We denote SðX Þ the class composed of all the subsets of E containing X. Definition 2 (Topological hull). Let X E and Y 2 SðX Þ. We say that Y is a topological hull of X , if Y has no holes and no cavities and if, 8y 2 Y n X , Y n fyg has a hole or a cavity. The following theorem gives a local characterization of the class of sets which are topological hulls relatively to the class of sets which have no holes and no cavities: Theorem 3. Let X E and Y 2 SðX Þ. Suppose that Y has no cavities and no holes. Then Y is a topological hull of X if and only if for each y of Y n X , y is an interior point or a 2D isthmus for Y , i.e., if and only if Tn ðy; Y Þ 6¼ 1. The interested reader will find the formal proof of this theorem in (Aktouf et al., 1996). Let R be a binary relation on a set and let U be an element of a set. We will denote RðU Þ ¼ fV =U RV g; we define recursively the binary relation Rk by V 2 Rk ðU Þ if 9W such that V 2 RðW Þ and W 2 Rk1 ðU Þ with R1 ¼ R. We define R1 ðU Þ ¼ fV ; such that 9k; V 2 Rk ðU Þ and RðV Þ ¼ ;g. The following definition introduces a binary relation H on SðX Þ used to formalize the process to build a topological hull: Definition 4 (Binary relation H on SðX Þ). Let X E. We define the binary relation H on SðX Þ: 8U ; V 2 SðX Þ; V 2 HðU Þ if 9u 2 U n X such that V ¼ U n fug and Tn ðu; U Þ ¼ 1. The following corollary is a direct consequence of the proof of Theorem 3: Corollary 5. Let X ; Y E and B 2 SðX Þ such that B has no cavities and no holes. If Y 2 H1 ðBÞ, then Y is a topological hull of X . The previous definition and corollary give a method for extracting a topological hull of a set
Z. Aktouf et al. / Pattern Recognition Letters 23 (2002) 523–531
527
Fig. 3. (b) and (d) are, respectively, topological hulls of (a) and (c).
X. We compute a bounding box B which has no cavities and no holes and which contains X. We iteratively delete the points of B which do not belong to X and which are not interior point and not 2D isthmus. We repeat this procedure until idempotence. Assuming sequential computational model, the following result may be easily proved: Proposition 6. A topological hull of an object X may be obtained in linear time with linear space with respect to the size of its bounding box B. In Fig. 3(b), an example of topological hull is given. Another example is given Fig. 3(d). It should be noted that a given object have several topological hulls and that the added surface produced by the method previously described may be not centered.
4. W-topological hull The most important drawback that we would want to avoid is to obtain a topologicall hull which has the tendency to ‘‘go away’’ from the object, as illustrated in Fig. 3(d). In the following, we are going to show how the distance transform can be introduced in the computation of the topological hull to control the process to obtain a centered topological hull. First, we introduce the notion of topological hull which is directed by a function. Intuitively, the function is used to induce an order in the processing of the points.
Definition 7 (W-topological hull). Let X E and let W be a function from E to N. We define the binary relation K on SðX Þ: 8U ; V 2 SðX Þ, V 2 KðU Þ if 9 u 2 U n X with Tn ðu; U Þ ¼ 1 such that: (i) V ¼ U n fug and (ii) WðuÞ ¼ maxfWðyÞ, 8y 2 U n X with Tn ðy; U Þ ¼ 1g. Let B 2 SðX Þ be a set with no holes and no cavities. If Y 2 K1 ðBÞ, we say that Y is a W-topological hull of X. Proposition 8. Let B E be a set with no holes and no cavities and let N ¼ #B. A W-topological hull of a set X included in B may be obtained in OðN log N Þ time complexity with OðN Þ space complexity provided the knowledge of WðbÞ for all b 2 B. Proof. Let X B. We suppose that we know WðxÞ, 8b 2 B. We consider a procedure which uses a list L, an element of L consists in the coordinates of some point. We will evaluate the complexity in terms of operations performed on the individual points in the list L. At any step of the procedure, the points of L are sorted according to their increasing W-values. We initialize L with all points belonging to the boundary of B. This may be done in OðN log N Þ time with linear space by using classical sorting algorithms. We also initialize Y ¼ B. The procedure consists in extracting the last point of the list; we check if this point satisfies the conditions for deletion. If yes, the point is removed from Y and the neighbors y of b which belong to N26 ðbÞ \ Y \ X are inserted in L according to their W-values. N26 ðbÞ is looked up as the deletion of b may affect any point within
528
Z. Aktouf et al. / Pattern Recognition Letters 23 (2002) 523–531
N26 ðbÞ \ Y \ X . This is independant from the connectivity of the object/background which only affects the topological numbers computation. The insertion of a point in L may be done in Oðlog N Þ time complexity (Cormen et al., 1990). The number of times a given point is inserted in L is bounded by a constant (26). The complexity may be improved if we have more information about W. Let us define a restricted class of functions W. This case is of interest, as it includes the classical discrete distance transforms used in the literature that we will apply later in the article. Definition 9 (Neighborhood transformation). Let W be a function from E to N. We will say that W is a neighborhood transformation if 9K 2 N=8x 2 E; 8y 2
N26 ðxÞ;
jWðxÞ WðyÞj < K:
Proposition 10. Let B E be a set with no holes and no cavities and N ¼ #B. Let W be a neighborhood transformation. A W-topological hull of a set X included in B may be obtained in OðN Þ time complexity with OðN Þ space complexity. Proof. As W is a neighborhood transformation, let K 2 N be the lowest value such that 8x 2 E; 8y 2 N26 ðxÞ; jWðxÞ WðyÞj < K. Then let N ¼ #B; M ¼ N K, and W0 ¼ minx2B fWðxÞg. We set W0 ¼ W W0 . Then, we have 8x 2 B; W0 ðxÞ < M. Note that OðMÞ ¼ OðN Þ. The procedure uses a table LðiÞ, i ¼ 1; . . . ; M, of M lists. At any step of the procedure, LðiÞ is a list of coordinates of points x such that W0 ðxÞ ¼ i. We initialize Y ¼ B and the table L is initialized with all border points of B. Let imax be the highest index i such that LðiÞ is not empty. The procedure consists in extracting the first point x of Lðimax Þ. If x was the last point of Lðimax Þ, we set imax :¼ imax 1. Then, we check if this point satisfies the conditions for deletion. If yes, x is deleted from Y and the neighbors y of x which belong to N26 ðxÞ \ Y \ X are inserted in lists according to their W-values. If a neighbor y has a value W0 ðyÞ > imax we set imax ¼ W0 ðyÞ. The procedure is repeated until all the lists LðiÞ are empty.
As the lists are indexed, all the list operations (insertion, extraction) are in constant time. The number of times a given point is inserted in L is bounded by a constant, that is its number of neighbors (26). After Lðimax Þ has been processed, Lðimax 1Þ may be empty, producing a possible sequence of successive empty lists. Nevertheless, it may be seen that each time a point is deleted, it induces at most K such empty lists, hence bounding the global number of empty lists. Thus, the algorithm has a linear time and linear space complexity. It should be noted that the W-topological hull of a given object may depend on the order in which the points of the lists will be processed. 5. Application of the distance transform Let d be a distance over E. The distance transform (Bergefors, 1984) is a function which takes as input a binary image I containing two classes of points: a subset X of object points and a subset X of background points. This process produces another image D in which the value of each point is its distance to the closest object points. The value of an object point is 0. The image D can be computed from the input image by starting from the object points and iteratively cumulating the elementary distances from the immediate neighbors: • for each x 2 X , we set DðxÞ ¼ 0; • we iterate until stability the following: 8x 2 X ; DðxÞ ¼ minfDðyÞ þ dðx; yÞ; 8y 2 CðxÞg; where CðxÞ is the successor function of the graph of a neighborhood relationship. In 3D, we will typically use C ¼ N6 or N26 and we will call D6 the distance in which all the points in N6 ðxÞ are considered distant of 1 from the point x, and we can define similarly D26 in which the points in N26 ðxÞ are at distance 1 from x. The choice of the distance values of the immediate neighbors allows to get closer to the Euclidean distance while maintaining acceptable computational cost by the use of integer arithmetic; for example this is the case of the so called chamfer distance (Chassery and Montanvert, 1991). We consider the W-topological hull as previously defined with W ¼ D. It corresponds to set a
Z. Aktouf et al. / Pattern Recognition Letters 23 (2002) 523–531
constraint in order to restrict the choice of the possible points to build the D-topological hull. The points are sorted according to their decreasing distance from the object. Then, the previously described algorithm to construct the D-topological hull is applied. The furthest points from the object are going to be processed first, and then, successively, the closest points. When a point is processed, it is eliminated if it does not satisfy the surviving criteria. In this case, its neighborhood is scanned and the points which are not already inserted in the list to be processed are then inserted. By ordering the points according to their distance from the object and processing them in the decreasing order, we are going to eliminate first the furthest points. The fundamental geometrical feature of this algorithm is that the first points which are going to survive to the elimination process will be the closest to the object. The reader can see on Fig. 3(b) that the resulting topological hull is very well centered within the original object. Let us notice that the distance transforms based on the classical discrete distances that we previously mentioned are neighborhood transformations (Definition 9) leading to a linear complexity algorithm (Proposition 10).
6. An -closing algorithm An interesting feature of the D-topological hull is that we have the ability to select the holes which
529
are going to be closed. Indeed, as we previously mentioned, the generic holes closing process is going to process all of the points which have been sorted according to the decreasing order of their distance value. Let us suppose now that we do not start from the highest distance value, but from an arbitrary given threshold . Then, all the points which are further from the object than will not be modified, and thus only selected smaller holes will be closed. We call the result of this modified version of the algorithm an -closing. The detailed algorithm is given in Appendix A. The Fig. 4 gives an example of application of our -closing algorithm to a biomedical image of a vertebra. Figs. 4(a)–(c) illustrate how the threshold can be used to keep relevant holes and eliminate the holes which correspond to noise.
7. Conclusion A 3D holes closing algorithm have been presented. It is based on some properties of the topological numbers and on the principle of the topological hull. One of its novel qualities is to allow the closing of the 3D-holes of different sizes, especially to remove smaller ‘‘noisy’’ ones while preserving those which are meaningful. Further, our algorithm is linear in time complexity with respect to the size of the bounding box of the considered object, on a sequential computational model.
Fig. 4. Application of the 3D holes -closing algorithm to an image of a vertebra. Case of the 26-connectivity for the object: (a) initial image (7 slices); (b) ¼ 2; (c) ¼ 3.
530
Z. Aktouf et al. / Pattern Recognition Letters 23 (2002) 523–531
Fig. 5. The 3D-holes closing algorithm.
Appendix A. An algorithm based on hierarchical lists Our algorithm is implemented by using a hierarchical list (Meyer, 1991; Vincent, 1990). A hierarchical list is an array of ordered lists. Each individual list is processed as a regular FIFO list. However, the lists themselves are totally ordered according to a certain criteria. We will call this criteria a priority level because it refers to the stage when a given list is going to be processed. In the case of our algorithm, the elements of the lists will be the points to be considered in the shrinking operation and the priority level of the list will be the distance value of the points in the list. For example, at the initial step, all the points with the same distance transform value belong to the same list. The input of the algorithm are a distance transform image W, a distance threshold , the object X and the bounding box B (which may be a solid cube), and the result is a centered topological hull Y of the initial image object in which the holes
smaller than are closed. The n-connectivity is used for the object and the n-connectivity for the complementary. In addition to that, we will use a binary image E which associates a boolean flag (TRUE or FALSE) to each point of the image to let us know if the corresponding point already belongs to the hierarchical list. We also use a function FIRST(L) to extract the first element of a list L. The algorithm is represented in Fig. 5.
References Aktouf, Z., Bertrand, G., Perroton, L., 1996. A 3d holes closing algorithm. In: 6th Discrete Geometry for Computer Imagery. Lecture Notes in Computer Science, Vol. 1176. Springer, Berlin, pp. 36–47. Bergefors, G., 1984. Distance transforms in arbitrary dimension. Comput. Vision Graphics Image Process. 27, 321–345. Bertrand, G., 1994. Simple points, topological numbers and geodesic neighborhoods in cubic grids. Pattern Recognition Lett. 15, 1003–1011.
Z. Aktouf et al. / Pattern Recognition Letters 23 (2002) 523–531 Bertrand, G., Malandain, G., 1994. A new characterization of 3d simple points. Pattern Recognition Lett. 15, 169– 175. Chassery, J.M., Montanvert, A., 1991. Geometrie discrete en analyse d’images. HERMES. Cormen, T.H., Leiserson, C.E., Rivest, R.L., 1990. Introduction to Algorithms. McGraw-Hill, New York. Kong, T.Y., 1989. A digital fundamental group. Computer Graphics 13, 159–166.
531
Kong, T.Y., 1997. Topology-preserving deletion of 1’s from 2-, 3-, 4-dimensional binary images. In: 7th Discrete Geometry for Computer Imagery. Lecture Notes in Computer Science, Vol. 1347. Springer, Berlin. Meyer, F., 1991. Un algorithme optimal de ligne de partage des eaux. In: RFIA, Vol. 2. AFCETpp. 847–857. Vincent, L., 1990. Algorithmes morphologiques a base de files d’attente et de lacets. Extension aux graphes, Ph.D. Thesis, Ecole Nationale Superieure des Mines de Paris.ty
