%* 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: