- Email: [email protected]
Following boundaries of discrete binary objects in space
Pattern Recognition, Vol. 23, No. 6, pp. 547-552, 1990 Printed in Great Britain
0031-3203/90 $3.00 + .00 Pergamon Press pie © 1990 Pattern Recognition Society
FOLLOWING BOUNDARIES OF DISCRETE BINARY OBJECTS IN SPACE ALAN BRYANT* and JACK BRYANT~f~ *Robotics Laboratory, Department of Mechanical Engineering and ~'Center for Approximation Theory, Department of Mathematics, Texas A&M University, College Station, TX ??843-3368, U.S.A.
(Received 26 July 1988; in revisedform 26 April 1989; receivedfor publication 1 June 1989) Abstract--Two algorithms for finding and following the boundaries of a set of discrete objects in ndimensional space are described; they are compared in three dimensions. In three dimensions, the new method is found to be about three times as fast as the straightforward approach. It operates on the boundary (surface) of the object rather than the whole object (i.e. the set of volume elements). The boundary, in turn, is found very rapidly using logical operations on the bit patterns by which the discrete object is stored. Applications to three dimensional shape recognition and display problems are discussed. Boundary following
Discrete objects
Boundary detection
INTRODUCTION In connection with a recent investigation,(1) we developed efficient algorithms to find the boundary of discrete simply connected planar binary objects. We have refined the method and can now locate and identify all the shapes in a 128 x 512 image in real t i m e - - t h a t is, as fast as the binary image can be obtained (computer: MIPS 2000). The key to the method is a fast boundary finding and extraction process and a fast classifier based on novel shape measures. In this paper we describe a generalization of the boundary finding and extraction portion to n-dimensions. Gordon and Udupa (2) have recently addressed this problem in three dimensions. Their paper contains a literature review concentrating on the problem of visualization; indeed, almost all of their 55 references are not about surface tracking at all. Their methods lead to a saving of about 55% over the simple approach. We have been unable to find previous work in which an approach comparable to ours has been taken to either the problem of finding the boundary as an n-dimensional object or to collecting the boundary points in a physically onedimensional array. Two methods for extracting the 3-D boundary (which in space might more naturally be called the surface) are compared and their relative efficiency is tested. Portable F O R T R A N code, although not given here, is available from the authors for the 3-dimensional case. Although the objects considered are presumed to exist in n-dimensional space, they are known to the computer as a binary array; that is, an array with each member being either 0 ("off") or 1 ("on"). Such :~To whom correspondence should be addressed. 547
binary objects arrise naturally in several types of medical image forming systems and in X-ray crystallography; they can also be generated artificially, as are the objects on which the methods in this paper were tested. We assume rectangular organization for the n-dimensional array, and that the point in the binary array corresponding to a cell in real or continuous space is on if the cell is contained in the real object and off if the cell is disjoint from the object. In our study of the planar case we found any sampling subject to these two restrictions adequate to represent real objects provided the mesh was fine enough; a similar result holds in higher dimensions for smooth objects with connected boundaries.
Digital topology Let Z = {... - 1, 0, 1, 2...} denote the set of integers. A discrete object is a finite subset of Z n. A discrete neighbor concept is determined by the number d of dimensions in which the coordinates of two points differ by at most one. In n-dimensions, a point has 2n nearest (d = I) neighbors and 2n z nearest or next nearest (d < 2) neighbors. For example, in the plane a point has four neighbors which differ by one in a single dimension and eight which differ in at m o s t two dimensions. In a three dimensional grid the.' corresponding numbers are 6 and 18. Reasons f o r concentrating on these choices of d are revealed when ~' one examines the connectivity of the boundary. (See* Rosenfeld (3'4) for more information on discrete anal-, ogs of topological concepts.)
Definition. A discrete set C is said to be d-connected if, when C is separated into non-void disjoint subsets A and B, then some point of A is a _< d-neighbor of some point in B.
548
ALANBRYANTand JACK BRYANT
Equivalently, C is d-connected if and only if each pair of points in C can be joined in C by a < dpath lying in C. The concept of connectedness and boundary thus depend on the value of d; one usually uses d = 1-neighbors for connectedness when considering solid objects. This selection minimizes computations and fits with some topological concepts of boundary and component.
Definition. The (internal) boundary of a discrete object is the set of points in the object having at least one 1-neighbor outside the object. This definition furnishes a straightforward method for finding the boundary of an n-dimensional object as is described in Algorithm 1 (Fig. 1). If p is a 1-neighbor of both ql and q2 and the three points ql-P-q2 are not on a line segment then ql is a 2-neighbor of q2. For this reason, if a point p is in the boundary then 2-neighbors of p are candidates to be boundary points. Each 2-component (i.e. maximal 2-connected subset of the boundary) C can be found by gathering the set B of all boundary points which join to any boundary point by a 2-path lying in B. This technique is described in Algorithm 2 (shown in Fig. 2). One feature which distinguishes space from the plane is the manner in which boundaries of domains are organized. In the plane, a smooth simply connected shape has as its boundary a simple closed curve; parameterization by arc length yields a one dimensional functional description of the boundary. The problems which stem from inadequate sampling are more difficult than might at first be thought, but can be overcome. In the discrete case we can extract an ordered "chain code" of the boundary. Nothing of the sort works in space; it is essential to note that the boundary of a discrete shape in Z", n _> 3, cannot be collected in a spatially natural order on a serial machine. When analyzing either the boundary or the solid shape itself one cannot process the points as fast as they are obtained; the best one can do is to store the location of unprocessed points in a queue for later processing. (The same procedure is sometimes unnecessarily followed in the plane.) This phenomenon is not an artifact of any particular method or algorithm. Experimentally, one finds (in Z a for shapes containing 106 points) that the process of collecting boundary points in a list gets several thousand points behind before it begins to catch up. Either the queue can be organized as a circular buffer or else a "garbage collector" can be used to control the size of the queue. We find the latter approach to be slightly faster; it is also safer with a clear indication of overflow. We adopt this strategy in the Algorithm 1 (Fig. 1). ALGORITHMS
According to Knuth, tSI an algorithm is "a preciselydefined sequence of rules telling how to produce
specified output information from given input information in a finite number of steps." He continues, "A particular representation of an algorithm is called a program, just as we use the word 'data' to stand for a particular representation of'information'." (See also Knuth. t6)) Thus an algorithm can be explained by giving a program written in a commonly used language; another presentation is to give Pidgin A L G O L pseudo-code. Neither extreme is attractive to the authors. The obvious advantage of the first method is that a reader fluent in that particular language will have little additional work to do in verifying our assertations. On the other hand, pseudo-code is more compact, and if one can obtain the requisite precision in the definition of the "sequence of rules" then the general reader may be better served since translation into any language should proceed routinely.* Our way out of this dilemma is to describe the methods as far as is possible in English, using pseudo-code to show the logical and loop structure, and to offer interested readers FORTRAN77 programs on request.
Practical considerations We assume the objects being studied are embedded in an array as bit maps, with one of the dimensions having the bits packed into w-bit computer words. It will develop that in most serial machines the proper value of w is the number of bits in the computer's data bus. Thus one of the dimensions r (for row), which we assume is the first (in storage order), is a multiple of the word size in bits. The number of computer words required to form a row is denoted by s (so that r --- s.w). In Z ~, for example, the two other dimensions are (say) c (columns) and p (planes). The array of computer words will be called A; it is indexed (in 3-D) as A[I,j,k] with 1 _< I _< s, 1 _
Algorithm 1 The boundary we locate is the internal 1-boundary defined above. If the objects have connected boundaries, then each object will be associated with one 2connected subset. The definition suggests that one locate the boundary by finding those points in the object which have a 1-neighbor not in the object. The problem with this approach is that the majority of points in the object may not be boundary points, leading to many needless tests. Two arrays the same * Of course, it is clear that many algorithms are not meant to be implemented by the reader, only considered.
Boundaries of discrete binary objects
Input: n-D binary array A Output: List of coordinates of boundary points bp, size boundary_size Coordinate queue: coord_que, size queue_size Pointers:
current, initially 1 last, initially 1 boundary.point, initially 0
Boolean variable: switch begin Copy array A to binary array B Locate a point p in B which is on Turn off the point p in B while last >_ current do switch := t r u e for each 1-neighbor q of p do if q is on in A t h e n last := last + 1 if last > queueJize t h e n if current = 1 t h e n Indicate shape too large; e x i t else Move the contents of the queue to the beginning last :-- last - current current
:= 1
endif Store the coordinates of q at location last in coord_que Turn off the point q in B endif elseif switch and the point at q in A is off t h e n switch := false boundary_point := boundary_point + 1 if boundary_point > boundary.size t h e n Indicate boundary too large; exit endif Store coordinates of p in bp at location boundary_point endif endfor current := current + 1 Fetch new point p from coord_que pointed to by current endwhile end Fig. 1. Algorithm 1--collecting the boundary from the solid.
549
550
ALAN BRYANTand JACK BRYANT
Input: n-dimensionalbinary array A Outputs:
List of coordinates of boundary points bp, size boundary_size n-dimensional binary array C containing the boundary
Pointers:
current, initially 1 last, initially 1
begin Using logical operations, form the array B of boundary points in A Copy the array B to C if the n-D boundary is needed as an array Locate a point p in B which is on Turn off the point p in B while last > current do for each 2-neighbor q of p do if q is on in B t h e n last := last + 1 if last > boundary_size t h e n indicate shape too large; exit endif
Store the coordinates of q at location last in.bp Turn off the point q in B endif e n d for
current :-- cur"rent + 1 Fetch new point p from bp pointed to by current endwhile end Fig. 2. Algorithm 2--collecting the boundary from the boundary as a discrete object. size as the original array will be needed to store the boundary being built and the places one has visited. (The need for the boundary as an array in Z ~ will be explained.) The original binary array should not be destroyed in the process, for it must be present to check whether a newly investigated point is a boundary point. However, only one new array is needed if the boundary itself is not needed as a spatially organized object. In fact algorithm 1 furnishes a list of coordinates of boundary points but not an n-dimensional binary arr.,y of the boundary. That is, the boundary is not organized naturally in a spatial sense. Using the definition, we find and collect the boundary at the same time. Algorithm 2
the logical and of O and the shifted set, then a left "end point" of O will be off in S, but points which have left neighbors will be on in S. On taking the exclusive or of S and O, one obtains the set of all points in O with no left neighbor in O--that is, the set of right-hand boundary points. The same process with left shift replaced by right shift yields the set of left-hand boundary points. The process works in any number of dimensions; indeed, the case n = 1 is the most troublesome, for in higher dimensions we process as many points as there are bits in the data bus in two logical steps (perdimension).* In the row direction, however, one must add special code for the low- and high-order bits shifted right and left. An exception to this occurs when the data bus is wide enough to hold one entire dimension.
A more elegant and much faster program for locating the boundary exploits the logical operations in Table 1 (logical and, shifts, and exclusive or) to achieve modest parallelism on a serial machine. To visualize how the process works in Z" it is enough to look at n = 1: if a set O of 1-D "pixels" is shifted to the right one position and a new set S is defined as
* More precisely, two high-level language statements per point must be executed. As usual, at the machine instruction level,data management loads, stores and indexing operations are needed, in addition to loop opening and closing operations, amounting slightly more than to eight machine instructions per word per dimension on a VAX; in this case, each word is 32 bits.
Boundaries of discrete binary objects
551
Table 1. Bit manipulation functions Name set reset point and eor
Ishift rshift
Action Sets a bit in an array Clears a bit in an array Tests a bit in an array (boolean) Bit-by-bit logical and of two words Bit-by-bit logical exclusive or Bit-by-bit logical left shift of a word Bit-by-bit logical right shift
Argument(s) (array, voxel) (array, voxel) (array, voxel) (word, word) (word, word) (word) (word)
Usage example set(A, i, j, k) reset (A, i, j, k) if point (A, i, j, k) then
W:= W:= W: = W:=
and(Y,Z) cot(Y,Z) Ishift(Y) rshift(Y)
linear relationship.) The break even point is p = 1.85; a solid with that small an interior is probably not being adequately sampled. The program which finds the boundary as an array is, as expected, almost instantaneous. For example, all the boundaries in a 128 × 128 × 128 volume are found in 2.4 s, or about 1.14/~s per voxel. This time is completely independent of the number of shapes in the region. One can imagine how the analogous planar boundary extraction method extracts boundary in Comparison of the methods in three dimensions real time when implemented on a faster computer. At first glance, the direct search of Algorithm 1 may An example shape had 153,635 points including 29,206 seem more efficient since only six neighbors need to boundary points: the program based on Algorithm 1 be examined; however, there are usually many more required over three C P U minutes, while the program points in the solid than in the boundary. In addition, using logical operations and Algorithm 2 used less two arrays must be managed in Algorithm 1 (the~ than a minute. (The time ratio predicted by the formula is 3.1031.) queue and the boundary point list); the result is that the more indirect two stage search is faster, at least in three dimensions, and furnishes the boundary as a Potential applications
After the boundaries are marked by the fast procedure just outlined the task of "collecting" them remains. In space, this requires looking at 18 neighboring points in our search for all boundary in one 2component of the set of all boundaries. The algorithm is structured like the procedure in Algorithm 1, with pointers to the current and last stored point and erosion of the existing boundary as the collection proceeds.
three dimensional image as well. The two approaches studied are not equivalent. Algorithm 1 will extract the entire boundary of a 1-component of the object; Algorithm 2 will extract a 2-component of the boundary of the object. For a l-connected object with a 2connected boundary the list of boundary elements furnished will be the same, and it is for such objects that we make these comparisons. Test data were generated to evaluate and compare the methods; the comparisons were made on a VAXstation 2000 under the VMS operating system. Various shapes were tested. We found the only variable which influenced the relative performance was the ratio of the number of points in the complete shape to the number in the boundary. Denote this ratio by p, the time required by the method of Algorithm 1 by ~t and the time required by the two-tiered search of finding the boundary using logical operations and collecting using Algorithm 2 by ft. We find a linear relation between the ratios:* = 0.6159p - 0.1368. (We tested methods on over one hundred 3-D shapes; the equation is the linear regression obtained using least squares. The correlation coefficient of 0.98 with sample size over one hundred indicates a significant * While the exact value of the parameters will no doubt depend on the machine, a linear relation should be expected to hold on any computer.
One application, for which data are potentially available, is the identification of a molecule from thresholded density 3-D images obtained from X-ray crystallography. Our approach is to develop 3-D "shape measures" analogous to the 2-D measures which are so successful. I1~ Part of the problem is the need to calculate integrals of the form s F(n, r)da, where F is a function of vectors r and n, r is the vector from the centroid of the object to a boundary element, n is the outward unit normal, and d a is the element of surface area. Many such integrals are naturally rotation- and translation-invariant, and with suitable normalization are magnification-invariant. The estimation of such an integral is difficult because the object is not continuous even though the underlying shape being represented may be C ®. Three related questions must be answered before the integral can be estimated numerically: - - H o w can one estimate the unit normal or, equivalently, the tangent plane? - - G i v e n the unit normal at each point, how can one assign a consistent "outward" direction? - - E v e n a continuous closed surface is often difficult to parameterize so as to compute da; what can one do to get the discrete element of surface area? The need for the boundary (or surface) as a 3-D object is now apparent: it is not enough to have the boundary as a list of points. One also needs the spatial
552
ALANBRYANTand JACK BRYANT
neighbors of a boundary point to estimate the tangent plane. One of us (J.B.) jointly with C. Krumvieda has recently solved the problem of estimating such integrals by using least-squares tangent plane approximation. A typical integral is
('L N
F
IIr - rol?l(r - r o ) ' n l ' d o "+q
(1)
where p and q are integers such that p > 0 and p + q > 0, IIr - r011 is the distance from the surface to the centroid, and n is the unit normal. In tests of the method, 36 integrals of this general form are computed; the time taken (including the use of Algorithm 2) varies depending on the size of the boundary and the ratio of the number of boundary points to volume points. For example, a sphere of nominal radius 16 required 8 C P U seconds to find the surface and volume and 18 C P U seconds to find the 36 integrals. The corresponding numbers for a sphere of radius 64 are 27 and 309 C P U seconds. The most complex shape we have processed to date took only 23 seconds for the surface-volume calculation but 389 CPU seconds for calculating the integrals. This amounts to about 10CPU seconds per integral on a VAXstation 2000, a machine slightly slower than the VAX11/780 (a popular machine used for benchmark comparisons of computers). (The computer performed floating point operations in software.) A report on this work is available from the authors. SUMMARY In this paper we have presented two algorithms
which extract boundary information from three dimensional discrete binary objects. The objects are assumed to be connected and have connected boundaries. One of the methods begins by finding a point in the object, and then simply examines every point in the object, looking for a neighbor not in the object; if the search succeeds the point is added to the boundary list. The other algorithm begins by finding the boundary of the object using logical operations on the bit patterns which encode the object. This step is a natural generalization of work of the authors in the plane, has been produced as a binary n-dimensional image, boundary points are extracted and placed in the boundary list. Although more neighbors must be checked to traverse the boundary image than the object, there are fewer boundary points in general. Tests in three dimensions of the method revealed a surprising gain in efficiency over the straightforward approach.
REFERENCES
1. A. Bryant and J. Bryant, Recognizing shapes in planar binary images. Pattern Recognition 22, 155-164 (1989). 2. D. Gordon and J. Udupa, Fast surface tracking in three-dimensional binary images, Comput. Vision Graphics Image Process. 45, 196-214 (1989). 3. A. Rosenfeld, Connectivity in digital pictures, J. Assoc. Comput. Mach. 17, 146-160 (1970). 4. A. Rosenfeld, Digital topology, Am. Math. Men. 86, 621630 (1979). 5. D. E. Knuth, Computer science and its relation to mathematics, Am. Math. Men. 81, 323-343 (1974). 6. D. E. Knuth, Algorithm and program; information and data, Communs ACM 9, 654 (1966).
About the Author--ALAN BRYANTis currently associated with the Robotics Laboratory, a division of the Mechanical Engineering Department at Texas A&M University. He has been studying robot vision problems, in particular the problem of recognizing planar shapes, for five years. Mr Bryant is a member of the Pattern Recognition Society. About the Author--JACK BRYANTis a mathematician at Texas A&M University. His main interest is the machine analysis of complex multi-dimensional natural images. He became interested in the problem of recognizing planar shapes only recently. His Ph.D. in Mathematics from Rice University, where he specialized in Fourier Analysis, was conferred in 1965. Dr Bryant is a member of the Pattern Recognition Society and is an associate editor of the journal Pattern Recognition.
0031-3203/90 $3.00 + .00 Pergamon Press pie © 1990 Pattern Recognition Society
FOLLOWING BOUNDARIES OF DISCRETE BINARY OBJECTS IN SPACE ALAN BRYANT* and JACK BRYANT~f~ *Robotics Laboratory, Department of Mechanical Engineering and ~'Center for Approximation Theory, Department of Mathematics, Texas A&M University, College Station, TX ??843-3368, U.S.A.
(Received 26 July 1988; in revisedform 26 April 1989; receivedfor publication 1 June 1989) Abstract--Two algorithms for finding and following the boundaries of a set of discrete objects in ndimensional space are described; they are compared in three dimensions. In three dimensions, the new method is found to be about three times as fast as the straightforward approach. It operates on the boundary (surface) of the object rather than the whole object (i.e. the set of volume elements). The boundary, in turn, is found very rapidly using logical operations on the bit patterns by which the discrete object is stored. Applications to three dimensional shape recognition and display problems are discussed. Boundary following
Discrete objects
Boundary detection
INTRODUCTION In connection with a recent investigation,(1) we developed efficient algorithms to find the boundary of discrete simply connected planar binary objects. We have refined the method and can now locate and identify all the shapes in a 128 x 512 image in real t i m e - - t h a t is, as fast as the binary image can be obtained (computer: MIPS 2000). The key to the method is a fast boundary finding and extraction process and a fast classifier based on novel shape measures. In this paper we describe a generalization of the boundary finding and extraction portion to n-dimensions. Gordon and Udupa (2) have recently addressed this problem in three dimensions. Their paper contains a literature review concentrating on the problem of visualization; indeed, almost all of their 55 references are not about surface tracking at all. Their methods lead to a saving of about 55% over the simple approach. We have been unable to find previous work in which an approach comparable to ours has been taken to either the problem of finding the boundary as an n-dimensional object or to collecting the boundary points in a physically onedimensional array. Two methods for extracting the 3-D boundary (which in space might more naturally be called the surface) are compared and their relative efficiency is tested. Portable F O R T R A N code, although not given here, is available from the authors for the 3-dimensional case. Although the objects considered are presumed to exist in n-dimensional space, they are known to the computer as a binary array; that is, an array with each member being either 0 ("off") or 1 ("on"). Such :~To whom correspondence should be addressed. 547
binary objects arrise naturally in several types of medical image forming systems and in X-ray crystallography; they can also be generated artificially, as are the objects on which the methods in this paper were tested. We assume rectangular organization for the n-dimensional array, and that the point in the binary array corresponding to a cell in real or continuous space is on if the cell is contained in the real object and off if the cell is disjoint from the object. In our study of the planar case we found any sampling subject to these two restrictions adequate to represent real objects provided the mesh was fine enough; a similar result holds in higher dimensions for smooth objects with connected boundaries.
Digital topology Let Z = {... - 1, 0, 1, 2...} denote the set of integers. A discrete object is a finite subset of Z n. A discrete neighbor concept is determined by the number d of dimensions in which the coordinates of two points differ by at most one. In n-dimensions, a point has 2n nearest (d = I) neighbors and 2n z nearest or next nearest (d < 2) neighbors. For example, in the plane a point has four neighbors which differ by one in a single dimension and eight which differ in at m o s t two dimensions. In a three dimensional grid the.' corresponding numbers are 6 and 18. Reasons f o r concentrating on these choices of d are revealed when ~' one examines the connectivity of the boundary. (See* Rosenfeld (3'4) for more information on discrete anal-, ogs of topological concepts.)
Definition. A discrete set C is said to be d-connected if, when C is separated into non-void disjoint subsets A and B, then some point of A is a _< d-neighbor of some point in B.
548
ALANBRYANTand JACK BRYANT
Equivalently, C is d-connected if and only if each pair of points in C can be joined in C by a < dpath lying in C. The concept of connectedness and boundary thus depend on the value of d; one usually uses d = 1-neighbors for connectedness when considering solid objects. This selection minimizes computations and fits with some topological concepts of boundary and component.
Definition. The (internal) boundary of a discrete object is the set of points in the object having at least one 1-neighbor outside the object. This definition furnishes a straightforward method for finding the boundary of an n-dimensional object as is described in Algorithm 1 (Fig. 1). If p is a 1-neighbor of both ql and q2 and the three points ql-P-q2 are not on a line segment then ql is a 2-neighbor of q2. For this reason, if a point p is in the boundary then 2-neighbors of p are candidates to be boundary points. Each 2-component (i.e. maximal 2-connected subset of the boundary) C can be found by gathering the set B of all boundary points which join to any boundary point by a 2-path lying in B. This technique is described in Algorithm 2 (shown in Fig. 2). One feature which distinguishes space from the plane is the manner in which boundaries of domains are organized. In the plane, a smooth simply connected shape has as its boundary a simple closed curve; parameterization by arc length yields a one dimensional functional description of the boundary. The problems which stem from inadequate sampling are more difficult than might at first be thought, but can be overcome. In the discrete case we can extract an ordered "chain code" of the boundary. Nothing of the sort works in space; it is essential to note that the boundary of a discrete shape in Z", n _> 3, cannot be collected in a spatially natural order on a serial machine. When analyzing either the boundary or the solid shape itself one cannot process the points as fast as they are obtained; the best one can do is to store the location of unprocessed points in a queue for later processing. (The same procedure is sometimes unnecessarily followed in the plane.) This phenomenon is not an artifact of any particular method or algorithm. Experimentally, one finds (in Z a for shapes containing 106 points) that the process of collecting boundary points in a list gets several thousand points behind before it begins to catch up. Either the queue can be organized as a circular buffer or else a "garbage collector" can be used to control the size of the queue. We find the latter approach to be slightly faster; it is also safer with a clear indication of overflow. We adopt this strategy in the Algorithm 1 (Fig. 1). ALGORITHMS
According to Knuth, tSI an algorithm is "a preciselydefined sequence of rules telling how to produce
specified output information from given input information in a finite number of steps." He continues, "A particular representation of an algorithm is called a program, just as we use the word 'data' to stand for a particular representation of'information'." (See also Knuth. t6)) Thus an algorithm can be explained by giving a program written in a commonly used language; another presentation is to give Pidgin A L G O L pseudo-code. Neither extreme is attractive to the authors. The obvious advantage of the first method is that a reader fluent in that particular language will have little additional work to do in verifying our assertations. On the other hand, pseudo-code is more compact, and if one can obtain the requisite precision in the definition of the "sequence of rules" then the general reader may be better served since translation into any language should proceed routinely.* Our way out of this dilemma is to describe the methods as far as is possible in English, using pseudo-code to show the logical and loop structure, and to offer interested readers FORTRAN77 programs on request.
Practical considerations We assume the objects being studied are embedded in an array as bit maps, with one of the dimensions having the bits packed into w-bit computer words. It will develop that in most serial machines the proper value of w is the number of bits in the computer's data bus. Thus one of the dimensions r (for row), which we assume is the first (in storage order), is a multiple of the word size in bits. The number of computer words required to form a row is denoted by s (so that r --- s.w). In Z ~, for example, the two other dimensions are (say) c (columns) and p (planes). The array of computer words will be called A; it is indexed (in 3-D) as A[I,j,k] with 1 _< I _< s, 1 _
Algorithm 1 The boundary we locate is the internal 1-boundary defined above. If the objects have connected boundaries, then each object will be associated with one 2connected subset. The definition suggests that one locate the boundary by finding those points in the object which have a 1-neighbor not in the object. The problem with this approach is that the majority of points in the object may not be boundary points, leading to many needless tests. Two arrays the same * Of course, it is clear that many algorithms are not meant to be implemented by the reader, only considered.
Boundaries of discrete binary objects
Input: n-D binary array A Output: List of coordinates of boundary points bp, size boundary_size Coordinate queue: coord_que, size queue_size Pointers:
current, initially 1 last, initially 1 boundary.point, initially 0
Boolean variable: switch begin Copy array A to binary array B Locate a point p in B which is on Turn off the point p in B while last >_ current do switch := t r u e for each 1-neighbor q of p do if q is on in A t h e n last := last + 1 if last > queueJize t h e n if current = 1 t h e n Indicate shape too large; e x i t else Move the contents of the queue to the beginning last :-- last - current current
:= 1
endif Store the coordinates of q at location last in coord_que Turn off the point q in B endif elseif switch and the point at q in A is off t h e n switch := false boundary_point := boundary_point + 1 if boundary_point > boundary.size t h e n Indicate boundary too large; exit endif Store coordinates of p in bp at location boundary_point endif endfor current := current + 1 Fetch new point p from coord_que pointed to by current endwhile end Fig. 1. Algorithm 1--collecting the boundary from the solid.
549
550
ALAN BRYANTand JACK BRYANT
Input: n-dimensionalbinary array A Outputs:
List of coordinates of boundary points bp, size boundary_size n-dimensional binary array C containing the boundary
Pointers:
current, initially 1 last, initially 1
begin Using logical operations, form the array B of boundary points in A Copy the array B to C if the n-D boundary is needed as an array Locate a point p in B which is on Turn off the point p in B while last > current do for each 2-neighbor q of p do if q is on in B t h e n last := last + 1 if last > boundary_size t h e n indicate shape too large; exit endif
Store the coordinates of q at location last in.bp Turn off the point q in B endif e n d for
current :-- cur"rent + 1 Fetch new point p from bp pointed to by current endwhile end Fig. 2. Algorithm 2--collecting the boundary from the boundary as a discrete object. size as the original array will be needed to store the boundary being built and the places one has visited. (The need for the boundary as an array in Z ~ will be explained.) The original binary array should not be destroyed in the process, for it must be present to check whether a newly investigated point is a boundary point. However, only one new array is needed if the boundary itself is not needed as a spatially organized object. In fact algorithm 1 furnishes a list of coordinates of boundary points but not an n-dimensional binary arr.,y of the boundary. That is, the boundary is not organized naturally in a spatial sense. Using the definition, we find and collect the boundary at the same time. Algorithm 2
the logical and of O and the shifted set, then a left "end point" of O will be off in S, but points which have left neighbors will be on in S. On taking the exclusive or of S and O, one obtains the set of all points in O with no left neighbor in O--that is, the set of right-hand boundary points. The same process with left shift replaced by right shift yields the set of left-hand boundary points. The process works in any number of dimensions; indeed, the case n = 1 is the most troublesome, for in higher dimensions we process as many points as there are bits in the data bus in two logical steps (perdimension).* In the row direction, however, one must add special code for the low- and high-order bits shifted right and left. An exception to this occurs when the data bus is wide enough to hold one entire dimension.
A more elegant and much faster program for locating the boundary exploits the logical operations in Table 1 (logical and, shifts, and exclusive or) to achieve modest parallelism on a serial machine. To visualize how the process works in Z" it is enough to look at n = 1: if a set O of 1-D "pixels" is shifted to the right one position and a new set S is defined as
* More precisely, two high-level language statements per point must be executed. As usual, at the machine instruction level,data management loads, stores and indexing operations are needed, in addition to loop opening and closing operations, amounting slightly more than to eight machine instructions per word per dimension on a VAX; in this case, each word is 32 bits.
Boundaries of discrete binary objects
551
Table 1. Bit manipulation functions Name set reset point and eor
Ishift rshift
Action Sets a bit in an array Clears a bit in an array Tests a bit in an array (boolean) Bit-by-bit logical and of two words Bit-by-bit logical exclusive or Bit-by-bit logical left shift of a word Bit-by-bit logical right shift
Argument(s) (array, voxel) (array, voxel) (array, voxel) (word, word) (word, word) (word) (word)
Usage example set(A, i, j, k) reset (A, i, j, k) if point (A, i, j, k) then
W:= W:= W: = W:=
and(Y,Z) cot(Y,Z) Ishift(Y) rshift(Y)
linear relationship.) The break even point is p = 1.85; a solid with that small an interior is probably not being adequately sampled. The program which finds the boundary as an array is, as expected, almost instantaneous. For example, all the boundaries in a 128 × 128 × 128 volume are found in 2.4 s, or about 1.14/~s per voxel. This time is completely independent of the number of shapes in the region. One can imagine how the analogous planar boundary extraction method extracts boundary in Comparison of the methods in three dimensions real time when implemented on a faster computer. At first glance, the direct search of Algorithm 1 may An example shape had 153,635 points including 29,206 seem more efficient since only six neighbors need to boundary points: the program based on Algorithm 1 be examined; however, there are usually many more required over three C P U minutes, while the program points in the solid than in the boundary. In addition, using logical operations and Algorithm 2 used less two arrays must be managed in Algorithm 1 (the~ than a minute. (The time ratio predicted by the formula is 3.1031.) queue and the boundary point list); the result is that the more indirect two stage search is faster, at least in three dimensions, and furnishes the boundary as a Potential applications
After the boundaries are marked by the fast procedure just outlined the task of "collecting" them remains. In space, this requires looking at 18 neighboring points in our search for all boundary in one 2component of the set of all boundaries. The algorithm is structured like the procedure in Algorithm 1, with pointers to the current and last stored point and erosion of the existing boundary as the collection proceeds.
three dimensional image as well. The two approaches studied are not equivalent. Algorithm 1 will extract the entire boundary of a 1-component of the object; Algorithm 2 will extract a 2-component of the boundary of the object. For a l-connected object with a 2connected boundary the list of boundary elements furnished will be the same, and it is for such objects that we make these comparisons. Test data were generated to evaluate and compare the methods; the comparisons were made on a VAXstation 2000 under the VMS operating system. Various shapes were tested. We found the only variable which influenced the relative performance was the ratio of the number of points in the complete shape to the number in the boundary. Denote this ratio by p, the time required by the method of Algorithm 1 by ~t and the time required by the two-tiered search of finding the boundary using logical operations and collecting using Algorithm 2 by ft. We find a linear relation between the ratios:* = 0.6159p - 0.1368. (We tested methods on over one hundred 3-D shapes; the equation is the linear regression obtained using least squares. The correlation coefficient of 0.98 with sample size over one hundred indicates a significant * While the exact value of the parameters will no doubt depend on the machine, a linear relation should be expected to hold on any computer.
One application, for which data are potentially available, is the identification of a molecule from thresholded density 3-D images obtained from X-ray crystallography. Our approach is to develop 3-D "shape measures" analogous to the 2-D measures which are so successful. I1~ Part of the problem is the need to calculate integrals of the form s F(n, r)da, where F is a function of vectors r and n, r is the vector from the centroid of the object to a boundary element, n is the outward unit normal, and d a is the element of surface area. Many such integrals are naturally rotation- and translation-invariant, and with suitable normalization are magnification-invariant. The estimation of such an integral is difficult because the object is not continuous even though the underlying shape being represented may be C ®. Three related questions must be answered before the integral can be estimated numerically: - - H o w can one estimate the unit normal or, equivalently, the tangent plane? - - G i v e n the unit normal at each point, how can one assign a consistent "outward" direction? - - E v e n a continuous closed surface is often difficult to parameterize so as to compute da; what can one do to get the discrete element of surface area? The need for the boundary (or surface) as a 3-D object is now apparent: it is not enough to have the boundary as a list of points. One also needs the spatial
552
ALANBRYANTand JACK BRYANT
neighbors of a boundary point to estimate the tangent plane. One of us (J.B.) jointly with C. Krumvieda has recently solved the problem of estimating such integrals by using least-squares tangent plane approximation. A typical integral is
('L N
F
IIr - rol?l(r - r o ) ' n l ' d o "+q
(1)
where p and q are integers such that p > 0 and p + q > 0, IIr - r011 is the distance from the surface to the centroid, and n is the unit normal. In tests of the method, 36 integrals of this general form are computed; the time taken (including the use of Algorithm 2) varies depending on the size of the boundary and the ratio of the number of boundary points to volume points. For example, a sphere of nominal radius 16 required 8 C P U seconds to find the surface and volume and 18 C P U seconds to find the 36 integrals. The corresponding numbers for a sphere of radius 64 are 27 and 309 C P U seconds. The most complex shape we have processed to date took only 23 seconds for the surface-volume calculation but 389 CPU seconds for calculating the integrals. This amounts to about 10CPU seconds per integral on a VAXstation 2000, a machine slightly slower than the VAX11/780 (a popular machine used for benchmark comparisons of computers). (The computer performed floating point operations in software.) A report on this work is available from the authors. SUMMARY In this paper we have presented two algorithms
which extract boundary information from three dimensional discrete binary objects. The objects are assumed to be connected and have connected boundaries. One of the methods begins by finding a point in the object, and then simply examines every point in the object, looking for a neighbor not in the object; if the search succeeds the point is added to the boundary list. The other algorithm begins by finding the boundary of the object using logical operations on the bit patterns which encode the object. This step is a natural generalization of work of the authors in the plane, has been produced as a binary n-dimensional image, boundary points are extracted and placed in the boundary list. Although more neighbors must be checked to traverse the boundary image than the object, there are fewer boundary points in general. Tests in three dimensions of the method revealed a surprising gain in efficiency over the straightforward approach.
REFERENCES
1. A. Bryant and J. Bryant, Recognizing shapes in planar binary images. Pattern Recognition 22, 155-164 (1989). 2. D. Gordon and J. Udupa, Fast surface tracking in three-dimensional binary images, Comput. Vision Graphics Image Process. 45, 196-214 (1989). 3. A. Rosenfeld, Connectivity in digital pictures, J. Assoc. Comput. Mach. 17, 146-160 (1970). 4. A. Rosenfeld, Digital topology, Am. Math. Men. 86, 621630 (1979). 5. D. E. Knuth, Computer science and its relation to mathematics, Am. Math. Men. 81, 323-343 (1974). 6. D. E. Knuth, Algorithm and program; information and data, Communs ACM 9, 654 (1966).
About the Author--ALAN BRYANTis currently associated with the Robotics Laboratory, a division of the Mechanical Engineering Department at Texas A&M University. He has been studying robot vision problems, in particular the problem of recognizing planar shapes, for five years. Mr Bryant is a member of the Pattern Recognition Society. About the Author--JACK BRYANTis a mathematician at Texas A&M University. His main interest is the machine analysis of complex multi-dimensional natural images. He became interested in the problem of recognizing planar shapes only recently. His Ph.D. in Mathematics from Rice University, where he specialized in Fourier Analysis, was conferred in 1965. Dr Bryant is a member of the Pattern Recognition Society and is an associate editor of the journal Pattern Recognition.