Computing shape description transforms on a multiresolution architecture

CVGIP:

IMAGE

UNDERSTANDING

Vol. 55, No. 3, May, pp. 287-295,

1992

Computing Shape Description Transforms on a Multiresolution Architecture L. CINQUE Dipartimento di Matemutica, Uniuersitri di Roma, Rome, Italy

C. GUERRA Dipartimento di Matematica, Universitd di Roma, Rome, Italy, and Department sf Computer Science, Purdue University, West Lafayette, Indiana 47907 AND S. LEVIALDI Dipartimento di Matematica, Universitd di Romd, Rome, Ituly Received June 19, 1990; acccptcd June 25, 1991

gree 4 and log II levels; the leaves are at level 0 of the tree and correspond to blocks that need no further subdivision and are therefore labeled BLACK or WHITE nodes. All nonleaf nodes are GRAY nodes and correspond to blocks consisting of both BLACK and WHITE pixels. The QMAT is the skeleton of the quadtree consisting of a subset of the BLACK nodes of the quadtree labeled with their associated distances from the background. The QMAT represents an image region as the union of squares (allowing overlap) having sides whose lengths are sums of powers of 2. Figures 3 and 4 show the QMAT for the image region of figure 1. The above representations have a number of applications in image processing, computer vision, and computer graphics as discussed in

This paper describes parallel algorithms for computing the medial axis transform and the quadtree medial axis transform on a pyramidal architecture; it also presents pyramidal algorithms to derive some geometric properties from the transformed image. For an n x n image, the algorithm for the quadtree medial axis takes O(log n) time on an n x n based pyramid, improving on a previous O(log* n) algorithm; all the other algorithms achieve 0( rP) time. 0 1992 Academic Press, Inc. 1. INTRODUCTION There are a number of image representations based on the concept of distance. Such representations involve various types of maximal block decomposition for an image region, where a block is some geometric primitive element. The medial axis transform (MAT) represents an image region R as the union of maximal blocks (squares or discs) that it contains [ 1, 91. The maximal blocks can be of any size and at any position and may overlap. Other methods, also based on distances, try to use the notion of a maximal block in a more systematic way by constraining the sizes and the locations of the primitive elements. The quadtree medial axis transform (QMAT) is one such representation [ 10, 1I]. The quadtree version of a binary image is formed by repeatedly decomposing the image into quadrants, subquadrants, etc., until each block contains a single value (either BLACK or WHITE, corresponding to the image region or the background, respectively), as illustrated for Fig. 1 in Fig. 2. For an n X n image the quadtree is a tree where each node has de-

[9, 101. In this paper we present parallel algorithms for computing the QMAT and the MAT on a pyramidal architecture; we also present pyramidal algorithms to derive properties such as the perimeter and the area of an image region represented by the MAT. For an n x II image the algorithm for the computation of the QMAT takes O(log n) time on a pyramid having at the base a mesh of size n x II. This result improves on a previous O(log2 n) time algorithm by Fan and Li [S]. Our algorithm is also simple since it only uses elementary operations. All other algorithms based on the MAT which are presented in this paper achieve O( n”*) time and use a divide-and-conquer technique. A number of parallel algorithms have been proposed in the literature for the determination of the MAT and the computation of geometric properties. Systolic algorithms 287

All

1049.Y660/92 55.00 Copyright 0 1992 by Academic Press. Inc. rights of reproduction in any form reerved.

288

CINQUE, GUERRA, AND LEVIALDI

FIG. 3. Image QMAT.

FIG. 1. Sample binary image.

The geometric structure of the pyramid naturally matches many types of concurrent computations in image processing; the hierarchical nature of the QMAT and MAT lends itself to efficient pyramidal implementations. 2.

COMPUTING

THE QMAT ON A PYRAMIDAL

for the MAT extraction using both serial and local operaCOMPUTER tors are discussed in [6]. A PRAM algorithm for the area Before presenting the algorithm, we introduce some and perimeter determination is presented in [4] achieving @log n log log n) time with O(n) processing elements. An notation and give a more precise definition of the QMAT. optimal O(n) time algorithm for the above two problems The distance transform DIST for a quadtree is defined as a function that gives for each BLACK node in the quadon an n x 12mesh is given in [ 121.An O(n) time algorithm using dynamic programming has also been proposed on tree the chessboard distance from the center of the correan n x n mesh or on a linear array of n processing ele- sponding block to the nearest point which is on the edge between the BLACK and WHITE elements. Thus the ments [4]. We will consider as a computational model a pyramidal width of a block is twice the associated DIST value, as computer [2]; it can be regarded as a combination of a may be seen on Fig. 4. A block has 8 neighbors, 4 neightapering array of mesh computers and a tree machine bors are adjacent along a boundary having a N, S, E, or architecture (for a survey on multiprocessor architec- W direction and 4 are adjacent along a corner pointing tures for computer vision see [3]). A pyramid of size n x NW, NE, SW, or SE. The QMAT is a new quadtree whose nodes are a subn (n a power of 2) has at the base a mesh of IZ x n processing elements (PEs) and a total of 4 n2-QPEs. The set of the nodes of the original quadtree. More precisely, base is considered to be at level 0; the levels above are the QMAT is defined as follows: let B = {b;} be the set of numbered from 1 to log n, where the log n level repre- the BLACK nodes in the original quadtree and let R(bi) sents the apex of the pyramid. Every PE at level i is be the square centered about block bi with side 2 DIST connected (via bidirectional unit-time communication (bi). The QMAT is the subset T = {t,} of B with the links) to its 9 neighbors: 4 neighboring PEs on the mesh at associated DIST values such that the following properthe same level i, 4 children at level i - 1, and a parent at ties hold: level i + 1. It is assumed that each PE has a fixed number of words of size O(log n), and all the arithmetical and Boolean operations and the information exchanges between adjacent PEs take unit time.

FIG. 2. Image quadtree.

FIG. 4. Block decomposition of QMAT with associated DIST values.

289

SHAPE DESCRIPTION TRANSFORMS

(1) the union of the squares R( tj) covers the entire image region; (2) the elements of T are the blocks with largest DIST values, that is, for any tj there is no element bi in B and not in T such that R( tj) is a subset of R( bi): this property is known as the subsumption property; (3) a block in B that does not belong to T is subsumed by a single element of T.

carried out during one bottom-up process; such phase is called pruning phase and takes an additional O(log n) time, leading to a total O(log n) time. We now describe phases 1 and 2 of step i that are similar to those found in [5]; we then present phase 3 of step i, and, finally explain phase 4 (the pruning phase).

Property (3) is used for the purpose of simplifying the algorithm for the derivation of the QMAT; however, it does not yield a minimal block decomposition, as shown in [lo]. We present an algorithm for computing the QMAT of an n x n image on an n x n pyramidal architecture (see Section 1). We assume that the n x n image is stored at the base of the pyramid, one pixel for each PE. A PE at level i of the pyramid is the apex of a subpyramid whose base is a block of size 2’ x 2’: thus it represents a node in the QMAT that corresponds to an image block of size 2’ x 2’. Note that neighboring nodes in the QMAT that are adjacent along a corner do not correspond to neighboring PEs of the pyramid. The algorithm consists of log IZ + 1 computation steps. Step i, 0 I i 9 log n, consists of three phases which are outlined here and described in more detail below:

At level 0 the input to phase 1 is the binary image (as seen on Fig. 1) and phase 1 simply labels the 1 pixels BLACK and the 0 pixels WHITE. At level i the inputs for phase 1 to a PE are the outputs of its four children from the level below of the same phase. The output of a PE, at level i, is either a common labeling of a single color (BLACK or WHITE) if all of its children have the same color, or a list of colors of its children that corresponds to a GRAY node in the quadtree. The action of this phase is depicted in Fig. 5. Not all BLACK nodes found at step i will be preserved in the final QMAT; in fact, BLACK nodes whose fathers are BLACK have to be removed from the QMAT. This will be done in the final pruning phase of the algorithm.

phase 1 (labeling phase) labels the image blocks of size 2’ x 2’ as BLACK, WHITE, or GRAY; phase 2 (distance computation phase) computes for each BLACK node at level i of the quadtree the associated DIST value; phase 3 (subsuming phase) determines whether a BLACK node at level i can subsume any BLACK node at level j < i. Phases 1 and 2 take O(1) time at each step and are performed in a bottom-up fashion along the pyramid in log n + 1 steps: at step i, the two phases are carried out by the PEs at level i. Last, phase 3 of step i (i > 0) checks whether a BLACK node at level i of the quadtree can subsume any node at a lower level. This is done by a top-down process along the pyramid, originated at level i, that terminates, in the worst case, at level 0: thus phase 3 takes O(i) time at step i. Phase 3 is repeated for each level i, in a bottomup fashion along the pyramid for a total of log n steps. The top-down processes, originated by phase 3 at successive levels, may be pipelined along the pyramid allowing different subsuming processes to be concurrently executed; thus an O(log n) time is obtained for computing the log 12+ 1 steps of the algorithm. During the above three phases, no pruning of the quadtree is done, that is, subsumable nodes are marked so but are not removed from the quadtree. Finally, at the end of the three phases of all the log n + 1 steps, phase 4 is

2.1 Labeling the Block (Phase 1)

2.2.

Computation of the DIST Function (Phase 2)

Distances will be computed at each step for any BLACK node that corresponds to a region whose extent is already known at that step; otherwise a temporary labeling will be performed with a value V larger than the image size. We first consider phase 2 at the base of the pyramid: the inputs to the PEs are outputs of phase 1 at level 0. Two cases may arise, one in which a BLACK node (corresponding to a single pixel) is embedded in eight neighboring BLACK nodes and another in which one neighboring WHITE node is present. For the first case the temporary value V is assigned to the node. For the second case (a BLACK node that has a neighboring WHITE node), phase 2 assigns a DIST value of 0.5, based on the

level i

level i-l FIG. 5. BLACK, WHITE, and GRAY labeling

290

CINQUE, GUERRA, AND LEVIALDI

assumption that the size of a block at the base is 1. Corner nodes in a quadtree do not correspond to neighboring PEs; thus, one routing step is needed for the PEs to collect the necessary information. At level i > 0, climbing along the pyramid, the inputs to a PE are the outputs of phases 1 and 2 at the previous level i - 1. Phase 2 performs two different actions: (1) examines whether a BLACK node p has any of its children having DIST < V; if so, it will assign to p the value DIST(p) = minh {DIST(ph)} + 2jm2, whereph, h = 1,. . . , 4, denotes a child of p, and 2i-2 is the block size at the level below, that is, one fourth of that at level i; otherwise it will assign value V to p. (2) tests if any BLACK node p with V DIST value has any of its neighboring nodes q labeled GRAY. The children of q that are neighbors of p (see Fig. 5) must all be BLACK (otherwise a DIST value less than V would have been assigned to one of the children of p at the previous step, and therefore to p itself). The algorithm will assign to P

DIST(p) = min {DIST(p), minh {DIST(qh)} + 2’-I + 2’-2). where the minimum is taken over all neighboring GRAY nodes q of p and qh , h = 1,2, . . . , 12, are the neighboring children of q. A correctness proof (based on exhaustion) of this phase is contained in [5]. 2.3. Determination (Phase 3)

of Subsumable

BLACK Nodes

Phase 3 determines whether a BLACK node in the quadtree can subsume any BLACK node below and, if so, marks this latter WHITE. Intuitively, let q be a node at level i and R(q) the black region at the base with side 2DIST( q); if R( q) contains a subregion R(p) corresponding to a BLACK node p at a lower level, we will say that p is subsumed by q. It is useful, at this point, to give a more precise definition of subsumption. Node p at level j is said to be subsumed by q at level i (j < i) if D(p, q) = DIST(p), where D(p, q) = DIST(q) - 2’-’ - 2-l. A node r subsumable by q may not be adjacent to q; this happens, for example, if r is subsumable by a node p adjacent to q and lower than q. In this case, it is D(p, q) = DIST(p) > 3-l. The node r does not belong to the QMAT and therefore has to be removed; however, this can only happen if the test of subsumption for r is carried out before the test forp, that is, beforep is removed from

FIG. 6. Corner block (in SE quadrant) can be subsumed by a node but not by its father.

the quadtree. Since our algorithm examines the nodes according to increasing size of the corresponding blocks, it will handle this situation correctly. Note that the test for subsumption at each level has to be done only for those BLACK nodes whose father is a GRAY node. A BLACK node identified at level i can disappear in the next step if its father is a BLACK node. For such node no test has to be done since it may subsume a node that cannot be subsumed by its father (an example can be seen in Fig. 6). A BLACK node belongs to the final quadtree structure if not all of its brothers are BLACK. This test can be done without moving data one level up along the pyramid, because all the required information is available at the same level. We are now ready to describe phase 3 of step i of the algorithm. Every PE at level i containing a BLACK node q whose father is a GRAY node determines whether it can subsume any adjacent BLACK node which is a descendent of a neighboring GRAY node p at the same level. It does so by sending a message to the neighboring GRAY node that initiates a top-down process. This message contains the following information: the level i at which the message is generated, the value DIST(q), and the direction (N, E, W, S, NE, NW, SE, SW) along which node q is adjacent to p. The GRAY node p sends the message down to its four children. This operation is repeated recursively for each GRAY node that (a) receives the message and (b) is adjacent to q, until a BLACK or WHITE node is encountered. If a BLACK node k at levelj receives the message then a test is done to determine whether k is adjacent to q and subsumable by q. If so, the BLACK node k is changed into a WHITE node. The top-down process stops when all BLACK or WHITE nodes are encountered which, in the worst case, occurs at level 0 of the pyramid. A GRAY node at level i may receive up to eight messages from neighboring BLACK nodes; top-down pro-

SHAPE DESCRIPTION TRANSFORMS

cesses are originated for all messages one after the other. Property 3 (see beginning of Section 2) guarantees the correctness of the algorithm since the order in which the test for subsumption is carried out for nodes at the same level is irrelevant. Any PE at a lower level of the pyramid also receives, at most, eight messages. Thus the total time for the top-down process initiated at level i is O(i) in the worst case. It has to be noted, fortunately, that the top-down processes can occur concurrently with the bottom-up process; that is, as soon as phase 2 at level i is done, the top-down process of phase 3 can start concurrently with phase 1 at the next level. This is a crucial observation that allows us to obtain a total of O(log n) time for the algorithm. Top-down propagation originated from successive levels can be pipelined along the pyramid without conflict or congestion of data. A message originated at level i + 1 reaches a lower levelj at least one unit of time later than any message originated at level i. At some point the PEs belonging to O(log n) different levels may be active simultaneously. The correctness of this phase is ensured by the fact that the tests for subsumption for nodes at different levels are performed on nodes corresponding to blocks of increasing size. At phase 3 no pruning of the quadtree is done; that is, BLACK nodes are marked WHITE as a result of a subsumption but are not removed from the quadtree. Therefore, there may be WHITE nodes whose sons are all WHITE. The pruning of the quadtree is postponed until the entire quadtree structure has been constructed and all the subsumptions have been discovered. At this point, phase 4 is performed. It consists of a bottom-up pass along the pyramid to remove all WHITE (BLACK) nodes with WHITE (BLACK) father. This is a simple test which can be done in constant time at each level of the pyramid, leading to an additional O(log n) time. 3.

COMPUTING

THE MAT

In this section, we present a parallel pyramidal algorithm for computing the MAT from an n x n binary image. Given an image region R, the MAT can be defined as the set of centers and radii of the maximal disks that are contained in R. This definition, first introduced by Blum [l], has been extended in several ways by different authors according to the application areas and the features to be used. One common definition, the one used in this paper, is that the MAT represents a region of a digital image as the union of the maximal upright black squares that are contained in that region. These squares may overlap and the number of them required to cover a specific region may be large. One of the most important consequences is that the computation of geometrical properties from the MAT is expensive, so that the MAT may be

291

not an efficient representation scheme. One way to reduce this time and speed up the computation is to perform parallel computation. Before presenting the pyramidal algorithm for the computation of the MAT, we briefly review the sequential algorithm [9]. For each pixel of the binary image I placed at the coordinates (x, y) in the plane, the algorithm determines the maximal upright square consisting of black pixels whose top-left corner is at (x, y). Let MAT(x, y) be the length of the side of such square. The algorithm initializes the entry MAT(x, y) to the value of the pixel 2(x, y), that is, either 0 or 1. It also creates a list of all black pixels in the image with their associated MAT values. It then performs a number of scans over this list: at each new scan it increments the entry MAT(x, y) if the values MAT(x + 1, y + l), MAT(x + 1, y) and MAT(x, y + 1) are all the same; otherwise, it removes the element from the list, since its MAT value is already the final one. The algorithm will stop when the list is empty. The running time of the sequential algorithm is then O(kn?), where k is the size of the largest maximal square in the image (n in the worst case). Our parallel algorithm is based on a bottom-up divideand-conquer technique that consists of O(log n) steps. During the first step, the image is divided into square subimages of constant size and maximal squares are found within each such subimage. During the next steps the partial results from four adjacent square subimages are combined until the entire image is assembled. The combining operation is based on the following simple observation. Consider four adjacent subimages (quadrants) that are combined into a single subimage M at a given step. The MAT value at (x, y) in M differs from the value computed at the previous step only if the maximal square at (x, y) extends over the border of the adjacent quadrants of M. Referring to Fig. 7, let (x, y) be the leftmost pixel along row x whose corresponding maximal square extends over the border and let dy (dx) be the

FIG. 7. A

squareextendingover the border of a quadrant.

292

CINQUE, GUERRA, AND LEVIALDI

distance of (x, y) from the right-vertical (bottom-horizontal) border. The MAT value at (x, y) has to be updated according to the formula MAT(x, y) = MAT(x, y) + min {MAT(x, y + dy), MAT(x + dx, y), MAT(x + dx, y + dy)}.

(1)

Note that for a pixel to the right of (x, y) along the row x in the same quadrant, the above formula has to be used with the same value of MAT(x, y + dy). For each row x, the algorithm determines the value MAT(x, y) of the leftmost point only. These are the only values that are needed at the next step; the correct values of all the other pixels will be determined later at the end of the bottom-up process. Similarly, along each column, only the MAT value of the topmost pixel whose corresponding square overlaps the border is needed for the next step. Thus the number of values needed at each step is proportional to the boundary size of the subimage considered. We now describe, in more detail, all the steps of the divide-and-conquer strategy. We assume from now on that n is a power of 4. Step 1 is performed by the PEs belonging to the base of the pyramid. The input image is partitioned into subimages of size 8 x 8 and, correspondingly, the base mesh is partitioned into blocks of PEs of size 8 X 8. The computation in each block proceeds in a systolic fashion [4] starting from the bottom right PE and moving antidiagonal-by-antidiagonal until the top left PE is reached (see Fig. 8). All the elements on an antidiagonal are processed simultaneously. At the start, the bottom right PE assigns the MAT value 1 or 0 to the pixel stored in it, depending on whether it is black or white. When a new antidiagonal is processed, the MAT value of a black pixel (x, y) will be computed from the MAT values of the

FIG. 8. Antidiagonal-by-antidiagonal pyramid.

processing at the base of the

FIG. 9. Bold cells contain leftmost points.

neighboring pixels already processed according to the following formula: MAT(x, Y> = 1 + min(MAT(x + 1, y), MAT(x, y + l), MAT(x + 1, y + 1). At the end of step 1, the leftmost (topmost) element in each row of an 8 x 8 subimage is identified. This is done by a simple process that originates at each right-vertical (bottom-horizontal) border PE of the block and stops when a pixel with a white neighbor is reached. A marker is used to indicate a row (column) that has a 0 value on the border. Finally, the values !, (0 5 x 2 15) of the y coordinate of the leftmost elements and their associated MAT values that have been obtained are moved within the block and stored one every fourth PE, as illustrated in Fig. 9. At each successive step the size of the subimage being analyzed is doubled, so that the size of a subimage analyzed at step i, 1 < i 5 log II - 2, is 22+1 x 22+i. The computation at step i is performed by the PEs belonging to the level Li/2_1 of the pyramid. The mesh at this level is partitioned into square blocks of PEs of size

Each block analyzes the subimage of size 22+’ x 2?+’ whose pixels are stored in a subset of the PEs at the base. Such a subset consists of all the PEs of the bases of the subpyramids having apexes inside the block. Step i uses results from step i - 1; such results are stored either at the same level Li/2-I of the pyramid or one

293

SHAPE DESCRIPTION TRANSFORMS

level below, depending on whether i is odd or even. In the former case, the values are already distributed among the PEs of the block one every fourth PE; in the latter case, they have to be moved one level up. In either case, every PE will contain at most one value. For each row of a left quadrant of M, step i updates the MAT value of the leftmost pixel according to (1). This computation can be performed in each submesh by using standard routing operations that take time of the order of the side of the submesh. A leftmost element in a quadrant (determined at step i - I) is not necessarily a leftmost element in M. This is the case, for instance, of a leftmost element in the top-left quadrant if the leftmost element of the same row in the top-right quadrant is not on the border between the two quadrants. The identification of leftmost elements in M is done at the end of step i and also takes time

r0

-1 0 0 0

100

0

0

-1

0

1

0

0

m

0

0

0

-1

0

0

1

-r

1

0

1

0

-2

0

0 0

0

0

-1

0

1

0

0

N

0

0

0

0

0

00

”

10000 3

4

5

67

0

0

0

0

12

-1

(b) r-

0

Y,

111110

0

0

v1

‘122000

-9

1

2

3

1

I?

0

0

0

1

N

0

1

1

1

-

0111110 3

4

12

0 0

1 1 I

5

1

0

1

0

1

0

67

Cd)

Crucial to this strategy is the reduction of the data to be computed at each stage. This number is proportional to the side of the corresponding subimage, which guarantees that no more than one value has to be stored at each PE. Let us analyze the time complexity of the algorithm. Step 1 takes O(1) time. The time for step i is proportional to the block size, that is O(22”ri’2’); then the total time of all the log IZ - 2 steps is given by 2

log /I/2

C

22trd21

i=l

which equals 0 ( H”~), At the end of the log n - 2 computational steps the MAT values will be available at level (log n)/2 of the pyramid, that is, a mesh of size n1’2x H”~. Remark that exactly 2n values are stored at (log n)/2 level corresponding to the leftmost and topmost pixel of each column and row respectively. From that level a top-down process originates to assign final values to all the other pixels in the image. This last process uses a standard pyramidal operation described in detail in [S], and referred to as funnel read in that paper. The funnel read operation obtains the final values for all the pixels in O(n1’2) time. 4.

AREA AND PERIMETER

FROM THE MAT

Given a set of rectangles represented by the coordinates of their corners, we compute the perimeter and the area of a region covered by the union of these squares. We assume that the corners of these squares are mapped into pixels of an n x n image as follows: Let R be a square whose corners (starting from the top-left corner)

FIG. 10. (a) Four rectangles Rl, R2, R3 and R4; (b) representation of the four rectangles on the mesh; (c) node values at the end of the vertical sweep; (d) horizontal and vertical edges of the contour marked after the horizontal sweep.

in the clockwise direction have coordinates (i, j), (k, j), (k, I), and (i, I) (with i < k and 1 < j). Then pixel I( i, j) (upper-left corner) has a value - 1, Z(k, j) (upper-right corner) has a value + 1, Z(k, I) (lower-right corner) has a value - 1, and I( i, j) (lower-left corner) has a value + 1 (Fig. 10(a) and (b)). Using this representation of a square, we assume that if (x, y) is the corner vertex of more than one square, then Z(x, y) has the algebraic sum of the values assigned by these squares. Like the mesh algorithms presented in [ 121,our algorithms first compute for each pixel a value Vxy that represents the number of squares it belongs to. From the Vxy values the area and perimeter are easily computed as will be described later. We first review the mesh algorithm to compute the Vxy values [12]. It consists of a vertical sweep of the image followed by a horizontal sweep (Fig. IO(c) and (d)). Vertical Sweep in an n X n Image. For x + 1 to 2n DO (in parallel) tx g 0 For y + 1 to 2n DO Processor PExy adds tx to its value Vxy; it sets tx to this new value and sends it to adjacent PE on the next row. The horizontal sweep computation is similar to the previous one except that data move horizontally from left to right.

294

CINQUE, GUERRA, AND LEVIALDI

The parallel algorithm that computes the Vxy values is based on a divide-and-conquer technique that achieves a total O(n”‘) time on a pyramidal architecture. The divide-and-conquer technique consists of O(log n) steps. In the first step the image is partitioned into adjacent square subimages of constant size, say 8 x 8, and in each subimage the partial values Vxy are computed. The first step consists, like the mesh algorithm, of a vertical sweep followed by a horizontal sweep inside each subimage. The computation of step 1 is carried out by the base mesh of the pyramid that is partitioned into adjacent blocks of size 8 x 8. The vertical sweep line, starting at the bottom boundary of a block, reaches the top boundary of the same block, where the obtained values are stored. We refer to such values as boundary values. This step takes constant time. The generic step i adds values from the four adjacent subimages Mh (h = 1, . . . , 4) of size 22-i x 22-’ analyzed at the previous step i - 1 to obtain new boundary values for the combined subimage. Precisely, if Mh is a subimage and Mk is its neighbor in the northern direction, the boundary values of Mh and Mk corresponding to the same column are added. Their sum is assigned to the top boundary of the combined subimage. Similarly, if M, is the subimage that is a neighbor of M,, in the eastern direction, their boundary values on the same row are added. The data movements along the pyramid necessary to carry out the above operations are similar to the ones described in the previous section and will not be repeated here. Since at each step the number of operations is proportional to the boundary size of the subimages being analyzed, the total time of the algorithm is O( ~l’~). Once the Vxy values are available, the computation of the area of a region is trivial since it requires that the number of pixels that have nonzero Vxy values be counted. This can be done on a pyramid in O(log n) time. Similarly in O(log n) time we can determine the perimeter by identifying and counting the border points of the region defined by the nonzero Vxy values. 5. MULTIPLE OBJECTS

In the previous sections we have implicitly assumed that a single object (black region) was present in the input image. In this section we show that the previous algorithms can be extended to the case of multiple objects. For this purpose, we use a technique first introduced by Miller and Stout in [7]. The interesting result is that the same time bounds as for single objects are obtained for the new algorithms. Since the condition of dealing with a single object was never imposed on the computational process, the algorithms of Section 3 for computing the MAT and the algorithm of Section 4 for determining the Vxy values can also

handle the case of multiple objects with no modifications. Both of them take O(n’12) time. Thus we only need to be concerned with the algorithms for determining the areas and perimeters of all the image regions from the Vxy values. Counting the number of pixels in a single region is done in O(log n) on a pyramid with a single bottom-up process. If we have many regions then there might be congestion toward the apex of the pyramid if many data try to move upwards: at the same time this can be avoided if the data move only halfway up the pyramid. One way of solving the problem is to set up a divideand-conquer strategy similar to that in Section 3. It consists of dividing the image into subimages, computing for each region the area within the subimage, and then combining results from adjacent subimages. For each row of the subimage, there can be at most one region extending over the border of the subimage, so that for a k x k subimage, we have a reduction of the data from k x k to O(k) (the boundary size of the subimage). As described in Section 3, this reduction allows O(n”*) time complexity to be obtained. 6. CONCLUSIONS

The shape of binary images may be quantitatively described by means of special transforms such as skeletons (S), medial axis transforms (MAT), and recently, based on quadtrees, the quadtree medial axis transform (QMAT). Much work has been performed in trying to speed up such computations, for instance by using parallel architectures like meshes, PRAMS, and pyramids. This paper presents an O(ni”) time algorithm for extracting the MAT and an O(log n) time algorithm for deriving the QMAT on a pyramidal computer. Basically, a bottom-up divide-and-conquer strategy has been used in the algorithm for obtaining the MAT. In order to extract the QMAT, the pyramidal architecture is exploited in a double concurrent pipeline mode: sending and computing data by overlapping bottom-up and topdown processing to obtain the labeled QMAT nodes. This algorithm is optimal in the asymptotic sense and, moreover, simple to implement. The area and perimeter from the MAT can also be obtained by an extension of the divide-and-conquer algorithm and the shape of multiple objects may also be extracted without excessively burdening the computation. We finally remark that the area and the perimeter of a region from the QMAT representation can be easily obtained in O(log n) time by adding contributions from all the BLACK nodes of the QMAT. REFERENCES 1. H. Blum, A transformation for extracting new descriptors of shape, in Models for the Perception of Speech and Visual Form, pp. 362370, MIT Press, Cambridge, 1967.

SHAPE DESCRIPTION TRANSFORMS 2. V. Cantoni and S. Levialdi (eds.), Pyramidal Systemsfor Computer Vision, NATO ASI, Series F, Vol. 25, Springer-Verlag, Berlin, 1986. 3. V. Cantoni and S. Levialdi, Multiprocessor computing for images, Proc. IEEE 76(8), 1988, pp. 959-969. 4. C. Chandran and D. Mount, “Optimal shared memory parallel algorithms and the medial axis transform, CAR-TR411, CS-TR-2165, University of Maryland, 1988. 5. N. P. Fan and C. C. Li, Computing quadtree medial axis transform by a multi-layered pyramid of LISP-processor arrays, in Proceedings, Conference Computer Vision and Pattern Recognition, 1988, pp. 628-634. 6. C. Guerra, Systolic algorithms for local operations on images, IEEE Trans. Comput. C-35, 1986, 73-77.

295

7. R. Miller and Q. F. Stout, Simulating essential pyramids in Pr? ceedings, Computer Vision and Pattern Recognition, 1988, pp. 912-917. 8. R. Miller and Q. F. Stout, Data movement techniques for the pyramid computer, SIAM J. Comput. 16(l), 1987, 38-60. 9. A. Rosenfeld and J. L. Pfaltz, Sequential operations on digital picture processing, J. Assoc. Comput. Mach., 1966, 471-494. 10. H. Samet, The quadtree and related hierarchical data structures, ACM Comput. Surv. 16(2), 1984, 187-268. 11. H. Samet, A quadtree medial axis transform, Comm. ACM 26(9), 1983, 680-693. 12. A. Y. Wu, S. K. Bhaskar, and A. Rosenfeld, Parallel computation

of geometric properties from the medial axis transform. Comprrt. Vision Graphics Image Process. 41, 1988, 323-332.