%* Last edited: Sep  5 21:07 1996 (misha)
\chapter{Vision Algorithm Implementation}


This section describes the implementation of several vision and
communication algorithms on the Abacus architecture. 
%It will eventually
%contain images processed by the high-level Abacus simulator, as well as
%characterizations of the inner loops.


\section{Basic Vision Algorithms}


%We described a number of vision algorithms on the Silt architecture in
%\cite{Bolotski93}. The Abacus machine alters implementation
%details of only the arithmetic operations; the vision algorithms
%remained unchanged.

\paragraph{Edge Detection}

The widely used Marr-Hildreth edge detection algorithm consists of
smoothing the image with a Gaussian filter, computing the Laplacian of the
filtered image, and locating the zero crossings of the result
\cite{Horn86}.  Convolution with a Gaussian can be approximated by repeated
convolution with a triangular filter, with weights $\frac{1}{4}$,
$\frac{1}{2}$, and $\frac{1}{4}$, requiring only arithmetic shifts and
accumulates. Since the Gaussian filter is separable, it can be implemented
by two one-dimensional convolutions.  Further, since zero crossings are
independent of absolute magnitude, the expensive scaling operation of the
Laplacian kernel is not required.

%\begin{figure}[htbp]
%  \setlength{\unitlength}{1mm}
%  \begin{picture}(150,90)
%%  \put(50,10){\epsffile{\mypsdirectory zerox1.ps}} \put(65,5){(d)}
%%  \put(100,10){\epsfffile{\mypsdirectory zerox2.ps}} \put(115,5){(e)}
%%  \put(50,50){\epsffile{\mypsdirectory gauss1.ps}} \put(65,45){(b)}
%%  \put(100,50){\epsffile{\mypsdirectory gauss2.ps}} \put(115,45){(c)}
%%  \put(0,10){\epsffile{\mypsdirectory bob.ps}} \put(15,5){(a)}
%  \end{picture}
%\caption[Edge Detection Example]{%
%Edge detection with the $\nabla^2G$ filter. (a) original image; (b) smoothed with a Gaussian,
%2 iterations; (c) smoothed with a Gaussian, 4 iterations; (d) and (e) zero
%crossings of $\nabla^2G$.}
%\end{figure}


%% \mypsfig{chip_photo_bw.ps}{Chip Photograph}{chip}{Chip photograph.}

%The algorithm was implemented with 16-bit clusters to maintain precision
%when adding and dividing 8-bit intensity values. One iteration requires 103
%cycles.  The calculation of the zero crossings of convolution with a
%discrete $3 \times 3$ Laplacian kernel required approximately 1200 cycles.
%Total edge detection time for a Gaussian with a $\sigma$ of 2.0, requiring
%four iterations of the triangle filter, is approximately 1800 cycles. The
%performance metric is for a 128 by 128 array processing 256 by 256 pixel
%images.



\paragraph{Surface Reconstruction}

Surface reconstruction is the task of computing the surface shape (height
and slope) from a set of potentially sparse and noisy measurements.  Harris
\cite{Harris86} introduced the coupled depth/slope model for surface
reconstruction and developed an iterative solution technique suitable for a
\mesh-based massively parallel computer. The update equations in each
iteration are simple and local and consist of five or six additions and
subtractions, followed by a division.

The summation of multiple values at a pixel is computed efficiently by the
multi-operand addition algorithm.  Division is expensive on a RBP
architecture, as it is on a word-parallel bit-serial machine.  This
algorithm only requires division by a constant, which can be implemented by
a special purpose sequence of shifts and adds.


\paragraph{Optical Flow}

In correlation-based optical flow \cite{Little88} the original image is
displaced in a variety of directions. At each displacement the difference
between the two images is computed and summed over a window.  If the images
are processed to produce binary image features, the difference at a pixel
is the exclusive OR of the shifted and unshifted images in a window around
the pixel.  neighborhood.  Finally, each pixel chooses the displacement
with the smallest difference value in a winner take all step.  Images
manipulated by the algorithm may consist either of simple brightness values
or of more complex features such as edges.  This algorithm has been
implemented on the simulator for the binary feature case. The edge
detection algorithm discussed earlier could be used to obtain the features
(in this case edges) from a raw intensity image.

\begin{Figure}{optflow}{Optical  Flow Example}{Optical flow calculation with maximum 
displacement of 2 pixels and a 5 \by 5 summation region. (a) original
image.  (b) displaced image.  The top shape moved up by one pixel and
partially off the image; the right shape moved down by one pixel and the
lower shape moved three units down and one to the left. }
  \setlength{\unitlength}{1mm}
  \begin{picture}(150,50)
  \put(20,10){\fbox{\epsffile{\mypsdirectory flow1.ps}}}
  \put(100,10){\fbox{\epsffile{\mypsdirectory flow2.ps}}}
  \put(35,0){(a)}
  \put(105,0){(b)}
  \end{picture}
\end{Figure}


\mypsfigy{2in}{optvec.ps}{Optical Flow Field}{optvec}{%
Optical flow field resulting from the raw data in \fig{optflow}. Some
confusion is caused by the movement of the top shape off the image, and by
the movement of the lower shape by a distance exceeding the maximum
displacement layer.}



\section{The DARPA IU Benchmark}

This section describes the subset of the DARPA Image Understanding
Benchmark that falls in the domain of early vision.

\paragraph{Connected Components: Broadcast.}

A common problem in image analysis involves labeling the connected regions
of constant pixel intensity, where unique regions are constituted from a
definition of either 4-connectedness (horizontally or vertically adjacent),
or 8-connectedness (horizontally, vertically, or diagonally adjacent).
This operation is discussed for binary images in the literature
(\cite{Levialdi72}, \cite{Leighton92}, \cite{Ziavras93},
\cite{Choudhary90}, \cite{Cypher90}) but has been extended for grayscale
values in the following three algorithms.  The first approach described
here is the most straightforward but requires the most time, the second
alternative is  faster but requires considerable amounts of memory,
and the third is a  modification of the second that compresses memory
requirements drastically at the cost of increased computation. 

In the simplest algorithm, each pixel of the intensity image is labeled uniquely
from its PE row and column.  This intensity and its corresponding label (of
$O(\log N)$ bits) are broadcast to each of its 8-connected neighbors. At
each PE, the neighbors of matching intensity are determined, and the
minimum of the their corresponding labels and the current label becomes the
new label.  This process continues until no more labels in the \mesh\ are
updated (checked by a global compare).  The number of broadcasting
operations required is proportional to the largest ``intrinsic diameter''
of all connected components in the image, defined as the maximal shortest
connected path between any two pixels in the region \cite{Leighton92}.  For
spirals and other high-curvature images, this intrinsic diameter can be as
high as $O(N^{2})$ in the worst case, resulting in $O(N^{2}\log N)$ bit
operations.

\paragraph{Connected Components: Shrinking.}
Levialdi's region shrinking operation \cite{Levialdi72} provides a
iterative method to directionally compress each region down to a single
pixel and then remove it entirely without fragmenting or fusing separate
regions.  If the results of each operation are saved away, the operation
can be reversed to generate a unique label when a region consists of only
one pixel, and that label can be transmitted to all possible neighbors in
the direction of expansion so that they can make a decision should they
become connected in the next stage.  The third algorithm modifies the basic
one by storing only a subset of the shrunk images and reconstructing on the
fly by repeating the shrinking operation.

%%See \fig{reg-shrink-alg}.

\paragraph{K-curvature Tracking and Corner Detection}

The connected components map is processed to produce $K$-curvature values
for those pixels on the component borders, which are then smoothed with a
Gaussian filter to eliminate multiple peaks near corners.  Pixels with
smoothed curvature values exceeding a threshold value (peaks) are
intersected with the zero-crossings of the first derivative of smoothed
curvature to extract candidate corners in the image.  Border pixels are
defined to be any pixel adjacent (N,E,W,S) to a pixel belonging to another
component.  $K$-curvature is defined at each border pixel as the interior
angle between the two lines passing through the current pixel and those $K$
border pixels away in either direction along the region's border.  See
\fig{k-curv-def}.

\mypsfig{K-curv.ps}{K-curvature
  definition}{k-curv-def}{K-curvature definition. In this figure, $\theta$
  is the K-curvature for $K$=3. Note that all pixels shown are edge
  pixels.}


\paragraph{Median Filter}

The median filter is a common image processing operation.  Unlike linear
filters, each pixel replaced not by a linear combination of its neighbors
but rather by the median value.  The operation is effective for removing
high frequency ``speckle'' noise without degrading the rest of the image.

\paragraph{Gradient Magnitude}

By computing the magnitude of a discrete intensity gradient in the image
and thresholding the result, strong direction-independent edges in the map
can be located.  Every pixel intensity in a 3x3 neighborhood about the
current PE is multiplied by a weight according to the Sobel X and Y masks
and summed to form the Sobel X and Y magnitudes.  Since the weights are
either 0, 1, or 2, the multiplications are converted to shifts.

The gradient magnitude is the square root of the sum of the squares of 
the X and Y magnitudes.  Results greater than a threshold value are 
flagged to create a boolean edges map.  The operation is constant in
space and time with respect to $N$.

\paragraph{Hough Transform}

The Hough transform partitions a binary image into discrete bands one pixel
thick and oriented at a certain angle, and sums the values within the band
to yield a set of projections that can be scanned to locate strong edges in
the original image at that angle \cite{Cypher89}.  Usually transforms are
computed for many angles, so that the aggregate data can provide insights
into the image properties.

The implemented algorithm partitions the image into bands of
constant line offset $\rho$ for a given angle $\theta$ according to the
equation: $\{(x,y): x \cos \theta + y \sin \theta = \rho\}$ where $(x,y)$
are PE coordinates and $\theta$ is assumed to be in the range $\pi/2 \leq
\theta < 3\pi/4$ (the other angles can be accommodated by pre-rotation of
the image).  A ``band total'' variable visits all PEs in its assigned band
by shifting east, and north if necessary (a result of the angle range is
that at most two pixels in the same column belong to the same band).  The
PEs in the first column are visited first, and then the variable travels
eastward across the columns.  Since there are many angles to be projected,
they are pipelined one column at a time, yielding $P$ projections on an
$N$x$N$ image in $O(N+P)$ time.

%%\input{hough_fig}

\section{Sorting and Routing Algorithms}

Sorting algorithms can be converted into routing algorithms by associating
data with the sort key and ensuring that the combined packet is moved
together.  This allows packets to be routed to a destination.  Two types of
sorting algorithms were examined: those with a fixed execution time, and
those with a data-dependent execution time.

\paragraph{ShearSort and RevSort}

%RevSort, like most of the oblivious mesh sorting algorithms depends on the
%compare-exchange primitive, which performs a two-element sort. Keys are
%compared, and the packet corresponding to the larger key is placed in the
%first data location, while the other packet is placed in the second
%location.

ShearSort is a $\sqrt{N}(\log{N} + 1)$ step sorting algorithm
\cite{Leighton92}. It consists of sorting rows in alternating opposite
directions starting at the left edge of the array, then sorting the columns
from top to bottom.  ShearSort relies on the TranspositionSort
\cite{Quinn87} procedure to get things sorted.  TranspositionSort is the
basic linear array sort.  It works similar to the uniprocessor BubbleSort:
comparing then exchanging (where necessary) alternating pairs of items.

An improvement of ShearSort is called RevSort \cite{Schnorr86}.  It is
identical to ShearSort, except CyclicSort is used to sort rows instead of
Transposition Sorts.  This algorithm finishes in $\sqrt{N}(\log{\log{N}} +
3)$ steps and the complexity of the each step is equivalent to each step of
ShearSort.  CyclicSort is a descendant of the TranspositionSort, and
performs as fast as its ancestor.  The smallest item is sorted to a chosen
processor instead of the leftmost processor.  



%%\paragraph{KielSort}

\paragraph{Mesh Greedy Routing Algorithm}

A very different packet-based routing algorithm was developed by Herbordt
\cite{Herbordt94}.  Each PE emulates two communication channels, one
vertical and one horizontal. Packets are moved one step through the X
channel until the correct X coordinate for the packet is reached. At that
point the packet is moved to the Y channel and proceeds vertically.  The X
and Y routing steps are interleaved, so that one of each occurs during each
iteration.  Notice that packets can be blocked from switching to the Y
channel if that section of the channel is full. When this occurs, packets
behind the blocked PE in the X direction are also blocked. The blocking
information propagates backward along the X channel one step at a time.
Since each PE contains space for two packets, collisions will not overwrite
data, and the blocking information can afford to propagate only one step
per cycle.

The key difference from the user's point of view is that this algorithm is
not guaranteed to finish in $O(\sqrt{N})$ iterations. In fact, the worst
case performance is $O(N)$, where $N$ is the number of PEs.  According to
Herbordt, the worst case is very unlikely to arise in either a completely
random, or typical routing conditions.  In fact, the algorithm generally
does not need to be iterated more than 2.5 times the optimal case of
$2\sqrt{N} - 2$. The paper demonstrates that this property holds for a
number of common permutations, such as transpose, bit shuffle, and
rotations.



%\section{Other Algorithms}

%\paragraph{Text Retrieval}

%With the advent of the Web, and ever-increasing amount of documents, the
%need for efficient document retrieval becomes more relevant. SIMD machines
%are very a natural fit to this task.  Several algorithms for this
%application have been developed for machines such as the CM-2.  An Abacus-1
%system has a capacity of 1 GB\footnote{A 256-chip Abacus-2 system has a
%  memory capacity of 32 MB.}. Performance numbers will be limited by memory
%  bandwidth to 4 bytes per cycle, or 500 MB/sec/chip, for a character
%  comparison rate of over 100 GB/sec.

%\paragraph{Seismic Migration}

%Discretized solutions of wave equations have many practical applications.
%All communication is local in space (within two grid distance) and in time
%(only three previous time points are necessary). The time points can be
%stored in memory, and dumped to off-chip memory as they recede into the
%past. A single update requires 11 move operations (some over distance two),
%10 summations, a multiply and a divide, both by compile-time
%constants. 

\section{Performance Summary}


In the vision algorithms table, computation times are shown for a 128 by
128 array processing 128 by 128 pixel images, without virtualization.
Processing larger images would not scale up linearly, since communication
time would be reduced, but data I/O time would be increased.
Memory-intensive algorithms would therefore degrade much faster.


%images.

%\begin{center}
%\small
%\begin{tabular}{|l|rr|rr|rr|}
%\hline
%\hdr{Operation} & \mcol{2}{|c|}{8-bit} & \mcol{2}{|c|}{16-bit} &
%\mcol{2}{|c|}{32-bit} \\
%%\cline{2-7} & \hdr{Cycles} & \hdr{GOPS} & \hdr{Cycles} & \hdr{GOPS} &
%%n\hdr{Cycles} & \hdr{GOPS} \\
%\hline 
%Add             & 4 & 4.0 & 4 & 2.0 &  5 & 0.7  \\
%Shift           & 2 & 8.0 & 2 & 4.0 & 2 & 2.0  \\
%Accumulate      & 3 & 5.2 & 3 & 2.6 & 3 & 1.3  \\
%Move            & 3 & 5.2 & 4 & 2.0 & 6 &  0.6 \\
%Compare         & 6 & 2.6 & 11 & 0.6 & 12 & 0.2\\
%%Multiply (16 $\times$ 16) & & & && 180 & 0.03\\
%\hline
%\end{tabular}
%\end{center}



\begin{Table}{visionperf}{Vision Algorithm Performance Summary}{Vision
    Algorithm Performance Summary}
\begin{tabular}{|l|r|r|r|}
\hline
\hdr{Algorithm} & \hdr{Cycles}  & \hdr{Time}  & \hdr{Memory}\\
                &  \hdr{1000s}  & \hdr{($\mu$sec)}  & \hdr{16-bit regs}\\
\hline
Edge Detection $\sigma=2.0$             & 0.45  & 3.6 & 3\\
Optical Flow, $\Delta=5$, 5 \by 5 region&  10  & 80  & 8 \\
Surface Reconstruction (1 iteration)    &  0.38 & 3   & 6 \\
Connected components, $D=256$           & 42.5 & 340  & 9 \\
Hough transform, $P=90$, multiplies     & 113  & 904  & 22 \\
Hough transform, $P=90$, differencing   & 25.6 & 204  & 22 \\
Hough transform, $P=90$, precompute     & 10.4 & 83  & 67 \\
K-curvature, $K=4$                      & 8.2  & 66  & 46 \\
Gradient Magnitude                      & 0.5  &  4  & 9 \\
\hline
Shear Permute, 1 data item             & 167  & 1340 & 15 \\
Rev Permute, 1 data item               &  78  & 624 & 15 \\
MGRA (average case),  1 data item       & 45  & 360 & 18\\
\hline
\end{tabular}
\end{Table}


% Local Variables: 
% mode: latex
% TeX-master: "main"
% End: 
