Efficient Image Processing Algorithms on the Scan Line Array Processor.
AbstractWe develop efficient algorithms for low and intermediate level image processing on the scan line array processor, a SIMD machine consisting of a linear array of cells that processes images in a scan line fashion. For low level processing, we present algorithms for block DFT, block DCT, convolution, template matching, shrinking, and expanding which run in real-time. By real-time, we mean that, if the required processing is based on neighborhoods of size mm, then the output lines are generated at a rate of O(m) operations per line and a latency of O(m) scan lines, which is the best that can be achieved on this model. We also develop an algorithm for median filtering which runs in almost real-time at a cost of O(m${\rm log}$m) time per scan line and a latency of ${\bf \lfloor^{\underline m}_{\,\, 2}\rfloor}$ scan lines. For intermediate level processing, we present optimal algorithms for translation, histogram computation, scaling, and rotation. We also develop efficient algorithms for labelling the connected components and determining the convex hulls of multiple figures which run in O(n${\rm log}$n) and O(n${\rm log}$2n) time, respectively. The latter algorithms are significantly simpler and easier to implement than those already reported in the literature for linear arrays.
Introduction
1.1 Motivation
Real-time image processing and understanding has long been regarded as a particularly
demanding problem for computer implementation, both because of the computational complexity
and because of the large I/O bandwidth required by most of the tasks involved.
Consider, as an example, the I/O bandwidth required to perform real-time HDTV simula-
tion. Such a task typically involves the handling of 1K by 1K frames at the rate of 60 frames
per second and results in a bandwidth requirement of approximately 500 Mbytes per second
for a progressively scanned image. Not surprisingly, such problems generally lie well beyond
the capacity of existing sequential processors. Consequently, a great deal of effort has been
devoted to developing parallel architectures and algorithms for real-time image processing.
The simplest category of the proposed architectures is the two-dimensional array, or mesh.
Examples of this class include the MPP [1], the CLIPP series [2], the MasPar [3], the DAP
[4], and the GAPP [5]. The general intent behind these SIMD (Single Instruction Stream,
Single Data Stream) machines is that the dimensions of the mesh should match those of
the input image, and that the pixels should be assigned to processors so as to maintain
the spatial relationships of the image. As a consequence, these machines can perform the
local window operations typical of low level image processing with extreme efficiency. Of
course, there is a considerable cost usually associated with such a large number of processors.
Furthermore, while nearest-neighbor links might make local inter-processor communication
quite fast, communication between two processors on either end of an n \Theta n array will require
n) time.
To reduce the cost of global communication while retaining the advantages of the mesh,
a second class of architectures has been proposed. This category improves on the mesh by
linking it to a sequence of progressively smaller meshes, each mesh having dimensions that
are one half those of its predecessor. As a result of this pyramidal structure, the global
communication costs on an n \Theta n array can be reduced from \Theta(n) to \Theta(log n) steps. Thus,
this pyramidal architecture is more appropriate for intermediate level tasks which require
the global exchange of information. Examples of this pyramidal machine include the PAPIA
[6] and the GAM [7] systems. Of course, when compared to the mesh, there is an increased
cost associated with the pyramid due to its increased complexity.
In general, high level image processing requires the simultaneous execution of several
independent tasks, which is more appropriate for a MIMD (Multiple Instruction Stream,
Multiple Data Stream) machine. But the disadvantage with using a purely MIMD machine
is that, while they can simulate SIMD operation, they do so only with a considerable overhead
for control and synchronization. Hence, they can not be expected to perform the data
parallel computations typical of low level image processing with the same efficiency displayed
by truly SIMD machines. Recognizing the relative advantages of SIMD and MIMD oper-
ations, a number of machines have been proposed which generally consist of an elaborate
combination of reconfigurable SIMD and MIMD modules. The basic idea behind creating
these hybrids is to enable the programmer to utilize the most efficient architecture for whatever
particular problem is presented. Examples of this class of architectures include the
Image Understanding Architecture (IUA) [8], the Associative String Processor (ASP) [9],
and NETRA [10]. While this strategy for architecture design probably offers the best hope
for achieving optimal performance across the spectrum of processing tasks, it will also be
the most expensive.
Fortunately, not all applications require all levels of image processing. For those situations
which are limited to low and intermediate level tasks, either of the SIMD architectures
already mentioned might well be sufficient. However, there is one particularly simple example
that avoids most of the cost constraints mentioned so far and that can still achieve optimal
performance for low and intermediate level operations. This architecture has been proposed
and in some cases implemented under a variety of names, including the Scan Line Array
Processor (SLAP) [11] [12], the Princeton Engine [13], and the Sarnoff Engine [14]. From
here on, this architecture will be referred to as the SLAP. The basic topology of this SIMD
machine is a linear array of processors, in which the number of processors corresponds to
the number of pixels in each row of the image. Incoming data is loaded line by line into
the processor array, with a distinct processor receiving each image column. This of course
results in a very high I/O bandwidth. The memory in this architecture is entirely distributed
among the individual processors, and each processor can communicate directly only with its
two immediate neighbors, whenever they exist. The most obvious disadvantage with this
interconnection scheme is the very small communication bandwidth, which would seem to
make global communication very inefficient. However, there are a number of compensations
when compared with the alternatives. Anything else would involve either a more elaborate
interconnection network or a shared memory, with an accompanying increase in cost and
complexity. Further, scaling the SLAP simply requires adding more processors to the end of
the existing array, and unlike bus based systems there are no obvious technical limitations
on how many such processors can be linked together without degrading the performance.
Hence, it seems reasonable to suspect that this linear array architecture might offer one of
the lowest cost options for achieving optimal low and intermediate level image processing.
However, for this to be true, it is necessary to argue not just that the technology is relatively
cheap, but that efficient algorithms can be written to run on this architecture.
Surprisingly, the design of algorithms which utilize the linear array architecture has received
only modest attention to date in the literature. Specifically, algorithms have been
proposed for sorting [15], matrix multiplication [16], and the Hough transform [17]. Algorithms
have also been proposed to solve a number of other graph [18] and geometric problems
[19], but only by assuming that the image data is already partitioned among the processors
according to the shuffled row-major distribution. The advantage of this shuffled row-major
distribution is that it significantly reduces the cost of global communication between the
different contiguous regions of the image. The disadvantage is that it is rather time consuming
to implement. To appreciate why, recall that the image is loaded line by line into the
SLAP. Hence, each pixel in a line is initially assigned to a processor of like index. To achieve
the shuffled row-major distribution by interprocessor broadcasting would then
operations per processor for an n \Theta n image, which would make it unsuitable for real-time
applications. By contrast, the initial distribution requires only a single load operation per
processor for each scan line. Accordingly, we assume this straightforward data distribution
as a starting point for developing efficient algorithms for a selection of low and intermediate
level image processing tasks.
1.2 Computational Model
The computational model used in this paper can be defined as follows. The SLAP architecture
is a linear array of n processors, where each image to be processed is of size n \Theta n.
Each processor is indexed in ascending order from 0 to (n \Gamma 1). After each scan line is
received, each processor P k latches in that value with the corresponding column index k.
Each processor P k is connected by a bidirectional communication link to processors P (k\Gamma1)
and P (k+1) , whenever they exist. Processors P (k\Gamma1) and P (k+1) are referred to as the left and
right neighbors of processor P k , respectively. The bandwidth of each communication link
is assumed to be a word, where a word is defined to be O(log n) bits, and each processor
can communicate a constant number of words to its immediate neighbors in a unit time.
The processors operate together in a SIMD fashion, and each processor is a general purpose
sequential processor with the capacity for conditional command execution and local address
generation. In a unit of time, each processor can compute O(1) basic logic or arithmetic
operations. Finally, each processor has associated with it a random access memory which
can hold O(n) words. In a unit of time, each processor can access its local memory to read
or write a word.
1.3 Problem Formulation
In this paper, the problem of developing algorithms for low level image processing is treated
separately from that of intermediate level tasks. In the case of low level operations, it is
assumed that we need to achieve the smallest possible latency consistent with an optimal
running time. To this end, we require that the processing of a scan line begin immediately
after it is received, and that the corresponding output line be generated in an amount of
time independent of the image size. More formally, suppose that the required processing is
based on neighborhoods of size m \Theta m. Then, we require that the output lines be generated
at a rate of O(m) operations per line after an initial delay of O(m) lines. If an algorithm
achieves this level of performance, then it runs in real time.
In the case of intermediate level operations, the entire image must be available for examination
before an output line can be generated. Hence, we assume at the outset that the
entire image is already stored in the local memories of the processors in the same way it
is received, with each column stored at the processor of like index. As already noted, this
straightforward input distribution requires only O(n) time. Further, we require that any
subsequent preprocessing be included in our time estimates. For these intermediate level
operations, we call the algorithm optimal either if it runs in O(n) time or it can be shown
that no faster algorithm is possible.
1.4 Results
Our results can be summarized as follows. For low level operations, we develop real-time algorithms
for block DFT, block DCT, convolution, template matching, shrinking, and expanding
[20]. We also develop an algorithm for median filtering which runs in almost real-time at a
cost of O(m log m) time per scan line and a latency of b mc scan lines [20].
For intermediate level operations, we develop optimal algorithms for translation, histogram
computation, scaling, and translation [20], although for the sake of brevity, we will
omit discussion of the latter two algorithms. We also develop algorithms for labelling the
connected components and for determining the convex hulls of multiple components which
run in O(n log n) time and O(n log 2 n), respectively [20]. The complexities of these last two
algorithms compare favorably with those of the existing algorithms, which respectively require
O(n) and O(n log n) time but which also assume the shuffled row-major distribution
[19]. In addition, our algorithms are significantly simpler and easier to implement and are
thus expected to perform better in practice.
Processing Operations
2.1 Median Filtering
Median filtering is the process of replacing each pixel of a given image by the median of the
pixels contained in a window centered around that pixel. This filtering operation is useful in
removing isolated lines or pixels while preserving spatial resolution. More specifically, let X
be an input image of size n \Theta n and let the window be of size m \Theta m, where m is assumed
to be odd for convenience. The filtered image Y is defined by:
(1)
1). The straightforward approach to solving this problem on the SLAP
is to compute the median independently for each output pixel. Specifically, to produce row
j of the output image, each processor P k first accumulates the m 2 pixel values necessary
to compute Y [j; k]. The median of these m 2 values is then computed by the best known
sequential algorithm in \Theta(m 2 ) time. When the next scan line is received, each processor P k
updates its current set of pixel values by replacing the least recent m values by the current
located to either side of processor P k . Clearly, this simple method yields the
desired latency of O(m) scan lines, but it also
operations per scan line.
There is, however, a more efficient procedure that allows the median filtering of the image
to be performed at the rate of O(m log m) operations per scan line with the same desired
latency. This procedure is based on the observation that the window surrounding pixel [j; k]
and the window surrounding pixel [(j differ only by 2m pixels. This immediately
suggests that we should store the m 2 values contained in the window about pixel [j; k] in
an appropriate data structure that allows us to (1) efficiently update the data structure to
reflect the window about pixel [(j quickly find the median of the values
stored in the data structure. One such data structure is the order-statistic tree [21]. It allows
us to dynamically (1) delete an element, (2) insert an element, or (3) locate the p th smallest
element (for any integer p) in time proportional to the logarithm of the tree size. Hence, by
using this data structure, we can produce an algorithm which will perform median filtering
in O(m log m) operations per scan line.
Assume for our real-time implementation of this operation that the input is the j th row
of the input image X. Initially, each processor P k broadcasts the value of its input pixel
X[j; k] to the set of (m \Gamma 1) processors f(k \Gamma
windows include X[j; k]. Since this essentially involves a left shift of the input scan line
m)c positions followed by a right shift of b( m)c positions, it clearly can be accomplished
in O(m) time. Next, each processor P k deletes the m oldest values in its order-statistic tree
and then inserts both the value of X[j; k] and the values received from its neighbors.
As already noted, this can be accomplished in O(m log m) time. Finally, in O(log m) time,
each processor P k determines the median of its updated order-statistic tree and then outputs
this result as the value of pixel [(j \Gamma b( m)c); k] of the output image. Thus, we have shown
the following theorem:
Theorem 1: Given a SLAP with n processors, the median filtering operation on an n \Theta n
image with a window of size m \Theta m can be performed in almost real-time at a cost of
O(m log m) operations per scan line and a latency of b( m)c scan lines.
2.2 Block 2D-DFT and 2D-DCT
A block transformation of a given type on an image involves partitioning the image into
non-overlapping blocks and then applying the desired transformation on each of the resulting
blocks. Two of the most widely used transformations are the two dimensional discrete Fourier
transform (2D-DFT) and the two dimensional discrete cosine transform (2D-DCT). In this
section, we present an algorithm which can implement either transformation on the m \Theta m
blocks of an image matrix in real-time. Since the two problems are essentially analogous, we
restrict our discussion to the task of computing the block 2D-DFT.
Let the n \Theta n input image X be partitioned into non-overlapping m \Theta m blocks, where
we assume for convenience that m divides n evenly. If we denote one such block as X (a;b) ,
g, then the 2D-DFT of X (a;b) is defined as:
1). The straightforward approach to apply this
transformation on the SLAP would be to independently evaluate the contribution of each
input line to the computation of each coefficient, which would
line. However, because the 2D-DFT is a separable transform , we can rewrite its definition
as follows:
Notice that, for a given row g of the input block and a given column d of the output block,
there is a single value of Y 0
(a;b) [d; g] which can be computed with O(m) operations. Once
a;b [g; d] has been calculated, we can then evaluate the contribution of the g th input row to
the d th output column with O(m) operations. Thus, by taking advantage of its separability,
we can compute the block DFT in real-time.
Assume for our real-time implementation of this operation that the input is the j th row of
the image X. Initially, each processor P k broadcasts the value of its input pixel X[j; k] to the
set of (m \Gamma 1) processors f(b(k=m)c \Theta m); ((b(k=m)c \Theta m)+ 1); :::; ((b(k=m)c \Theta m)+m \Gamma 1)g
which share the same block, which requires no more than O(m) time. Following this, each
processor P k computes the (k mod m) th coefficient of the 1D-DFT of the (j mod m) th row
of its block. Specifically, each processor P k evaluates the expression:
Clearly, this involves no more than O(m) operations. Next, each processor P k uses this
result to update the column of partially computed 2D-DFT coefficients held in the array
at its local memory. More precisely, for each value of c, where 0 - c -
performs the following operation:
Again, this clearly involves no more than O(m) operations. When the last row of the
block is processed, the values held in the array Partial \Gamma DFT k will be the (k mod m) th
column of the of the 2D-DFT of the block just processed. Finally, each processor P k outputs
Y ((b(j=m)c\Gamma1);b(k=m)c) [(j mod m); (k mod m)] as pixel [(j \Gamma m); k] of the output image. Hence,
we have the following theorem:
Theorem 2: Given a SLAP with n processors, the block DFT or DCT of an image of
size n \Theta n with each block of size m \Theta m can be computed in real-time at a cost of O(m)
operations per scan line and a latency of m scan lines.
2.3 Convolution and Template Matching
Convolution and template matching are fundamental image processing operations that are
computationally demanding. In this section, we present a method that performs the convolution
of an image of size n \Theta n with a kernel of size m \Theta m in real-time - that is, at the rate
of O(m) operations per scan line. The same method can be extended to perform template
matching within the same time bounds.
Given an image X of size n \Theta n and a kernel W of size m \Theta m, the convolution of X with
W is the image Y of size (n +m \Gamma 1) \Theta (n +m \Gamma 1) defined by:
Here we are assuming implicitly that
equal to zero whenever (j \Gamma r) or s) is not in the interval [0; (n \Gamma 1)].
A straightforward computation of the convolution would require operations and
therefore would require a minimum of \Omega\Gamma operations per scan line on the SLAP. Suppose,
instead, that we employ the overlap-and-add strategy [22]. Specifically, we partition the input
image X into non-overlapping m \Theta m blocks referred to as X 0 and indexed by (a; b). We
then convolve each block X 0
(a;b)
with the m \Theta m kernel W to obtain the
(a;b)
. If we then let a
is easy to verify that the value Y [j; k] defined by (6) can now be obtained by simply adding
the four entries Y 0
Up until now, nothing apparent has been gained computationally by redefining our matrix
convolution problem as we have in terms of block convolution. However, we can now "zero
pad" each block of X according to the following rule:
(a;b)
(a;b)
1). We can also "zero pad" our kernel W in an analogous manner to
obtain W 00 . Having done this, we are now able to obtain the linear convolution of X 0
(a;b)
and
W by performing the circular convolution of their "zero-padded" counterparts, X 00
(a;b)
and
W 00 . The advantage of circular convolution lies in the fact that the circular convolution of two
matrices can be obtained simply by finding their respective DFTs, multiplying the identically
indexed elements of these two DFTs, and finally finding the IDFT of the resulting product
matrix. This can be expressed more concisely as follows. Let
J denote circular convolution
and
N denote element-wise multiplication. Then:
(a;b)
O
We have already demonstrated in the previous section that we can compute the block DFT
of an n \Theta n image with m \Theta m blocks in O(m) operations per scan line. Hence, we would
also expect to be able to obtain the DFT of our 2m \Theta 2m "zero-padded" blocks in O(m)
time as well. Further, with the resulting DFT of these blocks distributed two columns per
processor, it is clear that the element-wise multiplication of the DFT of X 00
(a;b)
with the DFT
of W would require no more than O(m) operations per processor as well. Finally, since the
computation of the IDFT is essentially analogous to the computation of the DFT, we would
expect that a modification of our algorithm for finding the block DFT would also yield the
IDFT of (DFT (X 00
(a;b)
operations per input line. Of course, we would
need to complete our computation of (DFT (X 00
(a;b)
before we could begin our
computation of its IDFT. Hence, we would actually overlap our computation of (DFT (X 00
(a;b)
not with the computation of the IDFT of (DFT (X 00
(a;b)
rather with the computation of the IDFT of (DFT (X 00
N DFT (W 00 )). Thus, we see
that by converting our problem of convolution to one of block convolution and then taking
advantage of the properties of Fourier transforms, we can reduce the cost of convolution
from O(m 2 ) operations per scan line to O(m) operations per scan line. For a more detailed
presentation of this algorithm, see [23].
Hence, we have the following theorem:
Theorem 3: Given a SLAP with n processors, the convolution of an image of size n \Theta n
with a kernel of size m \Theta m can be computed in real-time at a cost of O(m) operations per
scan line and a latency of 2m scan lines. Similarly, the template matching problem can be
solved at the same rate for a template of size m \Theta m.
2.4 Shrinking and Expanding
Given a positive integer m, the m-step shrinking of an n \Theta n image X is the n \Theta n image
defined recursively as follows:
1). Similarly, the m-step expansion of an n \Theta n image X is the n \Theta n
image E (m) defined by replacing the minimum by the maximum in the above definition. It
is easy to verify that the m-step shrinking and the m-step expansion simply involve the
replacement of each pixel in X by the respective minimum or maximum of the (2m
pixels contained in the (2m centered about that pixel. Because of
the essential similarity between the shrinking and expansion problems, we will only discuss
an algorithm for the former.
Assume for our real-time algorithm that each processor P k holds in its local memory
a queue, referred to as Queue k , which holds the previous (2m input to that
processor. Once a new pixel X[j; k] is received at processor P k , Queue k is updated and the
minimum min k of the values in Queue k is computed. Next, each processor P k broadcasts
min k to the set of processors whose windows encompass
column k. When this is completed, each processor has (2m which represent the
minima of each of the columns contained in the window centered at that processor. Each
processor P k then computes the smallest of these (2m which it then outputs as
pixel of the output image. Hence, we have the following theorem:
Theorem 4: Given a SLAP with n processors, the m-step shrinking or expansion of an
image of size n \Theta n can be computed in real-time at a cost of O(m) operations per scan line
and a latency of m scan lines.
3 Intermediate Level Image Processing Operations
3.1 Convex Hull
Given a set S of points distinguished by a common label, the convex hull is defined as the
smallest convex polygon which contains all the points of the set. The extreme points of S
are defined as the corners of this smallest polygon. In this section, we consider an n \Theta n
input image X in which each pixel can have any one of O(n) possible labels. The pixels
which share a particular label are said to belong to a set, even though they may be spatially
unconnected to one another. For each of these O(n) sets, we wish to determine the extreme
points of its convex hull.
To do this, we present an algorithm which can compute these O(n) convex hulls in
O(n log 2 n) time. To simplify our presentation, we concentrate only on the problem of
determining the upper hulls, since the task of finding the lower hulls is entirely analogous.
Our divide-and-conquer algorithm first divides the input image into two subimages, one
consisting of the leftmost ( n) columns and the other consisting of the rightmost ( n) columns.
For each of the O(n) labels, the upper hulls of the two subimages are recursively computed
in parallel, after which a novel strategy is used to merge the two upper hulls. The merging
procedure consists of first concatenating the sequence of extreme points which define the
right upper hull to the sequence of extreme points which define the left upper hull. Having
done this, we then remove any subsequences of points fp such that all these
points lie on or below the line connecting the adjacent points p l and p r . In our algorithm,
this is performed by making O(log n) left-to-right and right-to-left sweeps across the image
for each label, each time eliminating a fraction of these non-extreme points. Extensive
pipelining is used to insure that in O(n) time a sweep can be completed for every one of
the O(n) labeled sets. In spite of the simplicity of our algorithm, its analysis requires a
somewhat tricky geometric argument to establish that O(log n) sweeps are sufficient. Hence,
we first present the algorithm and then provide a proof of its correctness.
3.1.1 Algorithm
We begin our algorithm with two preprocessing steps, each of which reduces the complexity
of our task. In the first preprocessing step, each processor P k sequentially examines the
correspondingly indexed column of the input image X in its local memory. For each label it
encounters, it retains the uppermost occurrence of that label as a potential candidate point
for the upper hull of that set. In the second preprocessing step, for each label, we make a
left-to-right sweep across the processor array. As we pass each processor P k , we identify, for
the candidate point (if any) at that processor, the candidate point from processors 0 through
k which has the maximum row index. Similarly, for each label, we also make a right-to-left
sweep across the processor array. Again, as we pass each processor P k , we identify for the
candidate point (if any) at that processor the candidate point from processors k through
which has the maximum row index. By pipelining the sweeps for each of the O(n)
labels, the whole process can be completed in O(n) time. Each processor P k then examines
each of its candidate points to see if it lies below the line connecting its respective left
and right maxima. If it does, then that point can be eliminated as a possible candidate
point. For each label, the sequence of candidate points left after these two preprocessing
steps actually consists of two characteristic subsequences, one of which may be empty. The
first subsequence is monotonically increasing, and the other subsequence is monotonically
decreasing. The significance of this property will become apparent when we discuss our
proof.
The algorithm for computing the extreme points of the upper hulls proceeds as follows.
Divide the n \Theta n input image X into two subimages X 1 and X 2 , where X 1 consists of columns
consists of columns ( n) through (n \Gamma 1). For each of the O(n)
labels, we recursively compute in parallel UH 1 and UH 2 , which are respectively the upper
hulls of the sets of points with that label found in X 1 and X 2 . When this computation
is completed, UH 1 and UH 2 are represented for a given label by the remaining candidate
points at columns 0 through ( n\Gamma 1) and columns ( n) through (n \Gamma 1), respectively.
To merge UH 1 and UH 2 and thereby obtain the completed upper hull UH, we proceed
as follows. First, we compute for each candidate point p i the coordinates of its immediate
successor S(p i ). For each label, this can be accomplished by making a single right-to-left pass
across the processor array, carrying along the coordinates of the most recently encountered
candidate point with that label. Again, by pipelining these sweeps for each of the O(n)
labels, the whole process can be completed in O(n) time. Next, for O(log n) repetitions,
we perform the following computation. For each label, we make a left-to-right pass across
the processor array, again carrying along the coordinates of the most recently encountered
candidate point with that label. Hence, when we arrive at a candidate point p i , we have
available the coordinates of its current predecessor P (p i ). Denote the angle formed by
the points P (p i opens away from the interior of the convex hull as
At each candidate point p i , we check to see if P (p i
This corresponds to asking whether p i lies on or below the line segment connecting P (p i )
to S(p i ). If it does, then we eliminate this point from consideration and we carry along the
coordinates of P (p i ) as the new value of the predecessor of S(p i ). When we complete this
left-to-right pass, we then begin an analogous right-to-left pass back across the processor
array. Now, however, we have to be concerned that the deletions of the previous left-to-right
pass may have left some successor values outdated. Hence, we bring along to each candidate
point p i the coordinates of its immediate successor. Pipelining again insures that we can
complete these two passes for all O(n) labels in O(n) time. And, after O(log n) repetitions,
the candidate points remaining for each set will be the extreme points of the upper hull of
that set.
3.1.2 Proof
We have made the claim that O(n log 2 n) time is sufficient to determine the extreme points
of the convex hulls of O(n) arbitrary sets. The basis of this claim is our contention that
O(log n) passes across the processor array are sufficient to merge any two arbitrary upper
hulls. We will first prove this latter claim, and then we will demonstrate that the overall
running time follows.
The proof of the upper limit on the running time of our merging algorithm is based on
a geometric argument, and therefore it is most clearly presented by a diagram such as that
shown in Figure 1. Figure 1 illustrates the process of merging two arbitrary upper hulls,
defined by the subsequence of extreme points fp 1
and UH 2 is defined by the subsequence of extreme points fp (j+1) g. The order
of the points in each subsequence corresponds to the relative order of their column indices.
We assume that UH 1 spans some portion of columns 0 through ( n\Gamma 1) and UH 2 spans some
portion of columns ( n) through (n \Gamma 1). As such, each upper hull can be defined by up to
n) extreme points. According to the way we have drawn UH 1 and UH 2 in Figure 1, the
upper hull UH resulting from merging UH 1 and UH 2 will be the line segment
For the sake of clarity, we have allowed a number of distortions in Figure 1. First, we
have represented UH 1 and UH 2 as continuous curves rather than as concatenations of line
segments. Second, we have allowed the angle p (j opening away from the interior
of the convex hull to be less than 90 ffi , when in fact our preprocessing makes this impossible.
Finally, we have drawn UH 1 and UH 2 so that merging will eliminate the entire subsequence
g. Yet, this subsequence is composed of two smaller subsequences
one of which is shown as being monotonically decreasing
and the other of which is shown as being monotonically increasing. To understand
the difficulty with this, recall that our preprocessing procedure insures that the original se-
Figure

1: Merging of Two Upper Hulls
quence of candidate points is monotonically increasing up to one or two maximum points
and then monotonically decreasing thereafter. Clearly, these maximum point(s) must also
be extreme points of the completed convex hull, implying that only subsequences of candidate
points which lie either to the left or right can qualify for elimination. Therefore,
any subsequence of non-extreme candidate points which need to be eliminated must be either
monotonically increasing or decreasing, and so too with g. However, a
situation such as that pictured in Figure 1 could still arise if we allow ourselves to rotate
the image so that the line segment p 1 p k is parallel to the horizontal axis of the image grid.
We permit all these changes described because they enhance the clarity of our presentation
without affecting the validity of our proof.
Consider a arbitrary left-to-right pass across the image as described in our algorithm.
This sweep extends a tangent to UH 2 from the second-to- last point remaining in UH 1 ,
thereby eliminating all the intervening points which fall below this line. Assume that an
arbitrary number of passes have been made across the processor array, each time eliminating
points from the end of UH 1 , and that now pL is the last point in the remaining sequence. The
next left-to-right pass extends the line segment P (p L )p T , as shown in Figure 1. Following
this, the succeeding right-to-left pass draws the tangent S(p T )p A . Two facts concerning
the point pA must hold. First, pA must lie to the left of the line segment P (p L )G, the
perpendicular to p 1 passing through point P (p L ). To see why, recall that any subsequence
of candidate points which requires elimination must be either monotonically increasing or
monotonically decreasing. As a consequence, the angle formed by any three points p u , p v ,
and pw such that u must be greater than 90 ffi , which in turn means that the
perpendicular to p 1 must divide the subsequence into exactly two
parts. Thus, every point remaining from UH1 aside from pL and P (p L ) must lie to the left
of the line segment P (p L )G, and so too with pA . We also require that the point pA lie on
or above the line segment p T C, which is the line segment parallel to passing through
. This is guaranteed by the fact that the angle formed by any three consecutive points on
the upper hull must be greater than 180 ffi . Hence, if p u , p v , and pw are any three consecutive
points which determine UH 1 , the slope of p v pw defined with respect to P (p L )G must be
steeper than p u p v . Consequently, the slope of S(p T )p A defined with respect to P (p L )G must
be steeper than must lie above p T C.
Following the right-to-left pass that has drawn the tangent S(p T )p A , the succeeding left-
to-right pass draws the tangent P (p A )p B . Notice that for the same reason that point pA
must lie on or above the line segment must lie on or above the line segment
pA H, which is the line segment parallel to passing through pA . Moreover, point pB
must lie to the right of the line passing through points P (p L ) and p T , since the angle formed
by any three consecutive points on the upper hull must be greater than 180 ffi . Finally, if we
define point F to be the point on the line passing through points P (p L ) and p T which is the
same distance from p T as P (p L ), then we can divide the region to the right of this line into
two halves as indicated in Figure 1. This enables us to make the following two observations:
(1) If pB lies in region 1, then the projection of P (p A )p B on p 1 p k must be greater than or
equal to the projection of P (p L )F . Since by definition the projection of P (p L )F is exactly
twice the projection of P (p L )p T , it follows that the projection of the second tangent P (p A )p B
must be at least twice the projection of the previous tangent P (p L )p T .
(2) Since point pA must lie above p T C and to the left of P (p L )G, it follows that the slope
of P (p A )p B defined with respect to P (p L )G must be less than the slope of DpB . If pB lies
in region 2, then the slope of DpB must be less than the slope of DF . Since by definition
the slope of DF is exactly half the slope of P (p L )F , it follows that the slope of the second
tangent P (p A )p B must be less than half the slope of the previous tangent P (p L )p T .
Hence, we can conclude that the tangent drawn by each left-to-right pass is either half the
slope or double the projection of the tangent drawn by the previous left-to-right pass.
Because the image is digitalized, it can be easily shown that the smallest possible nonzero
projection of a tangent on either p 1 p k or its perpendicular is \Omega\Gamma 1
Of course, the maximum
possible projection of p u p v on either p 1 p k or its perpendicular is obviously O(n). It follows
from this that the maximum and minimum possible slopes of p u defined with respect to
are O(n 2 )
repectively. And, therefore, since the tangent drawn by each
left-to-right pass is either half the slope or double the projection of the tangent drawn by
the previous left-to-right pass, it follows that O(log n) such passes are sufficient.
To compute the upper bound on the overall running time for determining the O(n)
possible convex hulls, we note that the solution of the merging problem for two n \Theta psubimages (2 - p ! n) still requires \Theta(n log n)time. Hence, the upper bound on the running
time for our divide-and-conquer algorithm is governed by the following recurrence:
O(n log n) if
whose solution is O(n log 2 n). Hence, we have the shown the following theorem:
Theorem 5: Assume that we have a SLAP with n processors and an input image of size
n \Theta n distributed one column per processor. If we also assume that each pixel is labelled
with one of O(n) possible labels, then we can determine the extreme points of the convex
hulls of all these labeled sets in O(n log 2 n) time.
We can also use this algorithm to solve another problem involving convex hulls. Consider
an n \Theta n image in which each pixel can belong to any one of O(n 2 ) sets distinguished by a
common label. If we require that all the pixels which belong to a particular set must form a
connected component, then with appropriate preprocessing we can use our algorithm to find
the extreme points of the convex hulls of all O(n 2 ) sets in O(n log 2 n) time. To justify this
claim, consider an arbitrary vertical line that divides the image into two parts. Because of
our requirement that the members of a particular set all belong to a connected component,
no more than n sets can have members on both sides of this arbitrary line. Hence, when our
algorithm calls for us to make a pass across the processor array for a given set, we obviously
only need to make a pass across this arbitrary line if this set is one of those n sets. Therefore,
as we move from one processor to another, the coordinates of only O(n) candidate points
and their labels need to be carried along, and hence O(n log 2 n) time suffices for this case as
well. Hence, we also have the following theorem:
Theorem Assume that we have a SLAP with n processors and an input image of size
n \Theta n distributed one column per processor. If we also assume that each pixel is labelled with
one of O(n 2 ) possible labels with the restriction that all those pixels which share a particular
label must form a connected component, then we can determine the extreme points of the
convex hulls of all these labeled sets in O(n log 2 n) time.
3.2 Connected Components
Given a binary image X, two pixels are called neighbors if
In other words, two pixels are said to be neighbors if they are adjacent
horizontally, vertically, or diagonally. Two pixels p (j 1 are then said to be
connected if there exists a sequence of pixels p (j h ;k h sharing some property such
that p (j (h\Gamma1) ;k (h\Gamma1) ) is a neighbor of p (j h ;k h ) . In this section, we present an algorithm which
labels the connected components of an n \Theta n input image on the SLAP in O(n log n) time.
Our algorithm for finding the connected components employs the divide-and-conquer
strategy. Assume that each processor P k holds the k th column of the input image. We divide
the n \Theta n input image X into two subimages X 1 and X 2 , where X 1 consists of columns 0
through ( n\Gamma 1) and X 2 consists of columns ( n) through (n \Gamma 1) Recursively and in parallel,
we identify and label the connected components of the two subimages X 1 and X 2 . When
this is completed, we identify and merge those connected components which span columns
n). This can be easily done using existing sequential algorithms in O(n) time
by either processor P ( n\Gamma1) or processor P ( n) . We then broadcast any label changes across
the two subimages. Since there can be no more than O(n) connected components spanning
these two columns and hence no more than O(n) label changes, this too can be accomplished
by pipelining in O(n) time. Clearly, then, the whole merging operation can be completed in
O(n) time.
To compute the upper bound on the running time for labeling the connected components,
we note that the solution of the merging problem for two n \Theta psubimages (2 -
requires \Theta(n) time. Hence, the upper bound on the running time for our divide-and-conquer
algorithm is governed by the following recurrence:
whose solution is O(n log n). Hence, we have the following theorem:
Theorem 7: Given a SLAP with n processors and an input image of size n\Thetan distributed one
column per processor, the connected components of this image can be labeled in O(n log n)
time.
3.3 Translation
Translation is the process of mapping each pixel at location [j; k] in the n \Theta n input image
X to a new position [j in the n \Theta n output image Y such that:
1)g. We assume for simplicity that a and b are integers
and that jaj; jbj - (n=2), since the cases of when either jaj ? (n=2) or jbj ? (n=2) can be
reformulated so that jaj; jbj - (n=2).
The algorithm proceeds as follows. Initially, we shift each row jbj positions to the right
correctly relocates all but jbj columns. Following this, we
wrap around the remaining jbj columns to the other end of the processor array by pipelining
their column-by-column broadcast. The task of translating the pixels within a particular
column by jaj positions is easily accommodated by insuring that as each pixel reaches the
correct processor, it is stored at the properly shifted memory location.
It is straightforward to verify that the upper bound on the running time of this algorithm
is O(njbj). For comparison, a lower bound can be derived by noting that the problem requires
us to move (n \Gamma jbj) columns of pixels a distance of jbj processors and jbj columns a distance
of (n \Gamma jbj) processors. Since each processor can perform a constant number transfers per
unit time and there are n such processors, it follows that the lower bound is \Omega\Gamma njbj). Hence,
we have the following theorem:
Theorem 8: Given a SLAP with n processors and an input image of size n \Theta n distributed
one column per processor, the translation of this image by a distance a in the vertical
direction and a distance b in the horizontal direction can be completed in \Theta(njbj) time. This
bound is optimal.
3.4 Histogram Computation
The histogram of an n \Theta n input image X is a computation which determines the relative
frequency of occurrence of the various possible pixel values in the image. To develop an
algorithm for this, we assume for simplicity that the possible pixel intensities are constrained
to integer values ranging from 0 up to (m \Gamma 1), where nothing is assumed about the relative
sizes of m and n.
The algorithm proceeds as follows. First, each processor P k computes the histogram of
the k th column of the input image. To do this, each processor simply examines each of these
n pixels and increments the entry in an array called whose index corresponds to the pixel
value. Next, we compute the image matrix histogram by adding together the identically
indexed values in each of the local histograms H k . This simply requires m right-to-left
passes across the processor array, which will then leave the completed histogram values in
the array H 0 kept at processor P 0 . If we pipeline these m passes, then the whole procedure
can be completed in O(m n) time, which is clearly optimal. Hence, we have the following
theorem:
Theorem 9: We are given a SLAP with n processors and an input image of size n \Theta
distributed one column per processor. If we assume that the possible pixel values are
constrained to integer values in the interval [0; (m \Gamma 1)], then the histogram of this image
can be computed in O(n +m) time.
In this paper we have presented efficient algorithms for the SLAP for a variety of low and
intermediate level image processing tasks. For the low level operations, our algorithms run in
real-time, which is clearly the best that can be achieved on this model. For the intermediate
level operations, our algorithms are either optimal or compare favorably with the existing
algorithms. Moreover, our algorithms achieve their performance without assuming that the
input image is already partitioned according to the shuffled row-major distribution. Taken
together, our results suggest that the SLAP is a promising architecture for real-time image
processing.



--R

"Design of a Massively Parallel Processor"
"The CLIP7A Image Processor"
"The Design of the MasPar MP-1: A Cost Effective Massively Parallel Computer,"
Distributed Processor Array"
"The Geometric Arithmetic Parallel Processor"
"PAPIA"
"The GAM Pyramid"
"The Image Understanding Architecture"
"The Asp: A Cost-Effective Parallel Microcomputer"
"NETRA: An Architecture for a Large Scale Multiprocessor Vision System"
"Real-Time Image Processing on Scan Line Array Pro- cessors"
"Scan Line Array Processors for Image Computations"
"The Princeton Engine: A Real-Time Video System Simulator"
"The Sarnoff Engine: A Massively Parallel Computer for High Definition System Simulation"
"Optimal Sorting Algorithms for Parallel Computers"
"Modular Matrix Multiplication on a Linear Ar- ray"
"Computing the Hough Transform on a Scan-Line Array Processor"
"Optimal Graph Algorithms on a Fixed-Size Linear Array"
"Optimal Geometric Algorithms for Digitized Images on Fixed-Size Linear Arrays and Scan-Line Arrays"
Fundamentals of Digital Image Processing
Introduction to Algorithms

"Efficient Image Processing Algorithms on the Scan Line Array Proces- sor"
--TR

--CTR
Ronald I. Greenberg, Finding connected components on a scan line array processor, Proceedings of the seventh annual ACM symposium on Parallel algorithms and architectures, p.195-202, June 24-26, 1995, Santa Barbara, California, United States
Francesco Gregoretti , Roberto Passerone , Leonardo Maria Reyneri , Claudio Sanso, A High Speed VLSI Architecture for Handwriting Recognition, Journal of VLSI Signal Processing Systems, v.28 n.3, p.259-278, July 2001
Shorin Kyo , Shin'ichiro Okazaki , Tamio Arai, An Integrated Memory Array Processor Architecture for Embedded Image Recognition Systems, ACM SIGARCH Computer Architecture News, v.33 n.2, p.134-145, May 2005
parallel algorithms;image processing;scan line array processors;video procesor;linear array;SIMD algorithms