Jan Bartovsky
email: j.bartovsky@esiee.fr

Petr Dokládal
email: petr.dokladal@mines-paritech.fr

Eva Dokladalova
email: e.dokladalova@esiee.fr

Michel Bilodeau
email: michel.bilodeau@mines-paritech.fr

Mohamed Akil
email: m.akil@esiee.fr

Jan Bartovský

Eva Dokládalová

Real-Time Implementation of Morphological Filters with Polygonal Structuring Elements

Keywords: Mathematical Morphology, Hardware Implementation, Alternating Sequential Filter, Parallel Computation, Polygonal Structuring Element

HAL is

Introduction

Since its first introduction in the late 1960's, mathematical morphology has becom in the field of image processing a useful tool for analysis of the shape or the form of spatial structures [START_REF] Matheron | Random sets and integral geometry[END_REF][START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF][START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF]. Over time it has found its application as a widely-used image processing technique [START_REF] Coster | Image analysis and mathematical morphology for civil engineering materials[END_REF][START_REF]Mathematical Morphology: From Theory to Applications[END_REF].

Thanks to the recent technological development of sensors, the resolution of images increased to tens of megapixels. Certain morphological operations, e.g., top-hat transformation, ultimate openings, granulometry, alternating sequential filters (ASF) [START_REF] Serra | An overview of morphological filtering[END_REF], etc., on such large images require a large structuring element (SE), since its size should be proportional to the size of the image and its contents.

Existing hardware implementations either support rectangles using SE decomposition (fast computation, but angular anisotropic), or support arbitrarily-shaped SEs at the cost of significant performance decrease. Our work supports polygonal SEs at the performance rate of the rectangular SEs.

The paper is organized as follows: Section 2 makes a short survey of existing morphological algorithms and architectures. Section 3 outlines the basic aspects of morphological dilation and erosion, and show how to decompose the polygons into a set of lines. Section 4 describes the algorithm, and its use to decompose polygons while preserving the sequential access to data, minimal memory consumption and latency. Section 5 gives the functional implementation of the algorithm. The principal result, a parallel version using two levels of parallelism (temporal and spatial) is presented in Section 6. Finally, Section 7 presents experimental results obtained on an FPGA.

State of the art

First algorithms of morphological dilation were based on the definition (see Sec. 3), efficient for small structuring elements. The high computational complexity of large or arbitrarily-shaped SE motivated the search of decompositions, either by i) separation into lower dimensions, or ii) by decomposition into atomic shapes, using the Minkowski addition. This section reviews the literature on the most known algorithmic decompositions and most efficient implementations.

i) The separation into lower dimensions allows the decomposition of n-dimensional rectangular SE into 1-D segments. The simplest method to compute 1-D dilation is an exhaustive search for maximum in the scope of SE B according to the definition (Eq. 1). Clearly, this naive solution tends to need a large number of comparisons. The number of comparisons is considered as a metric of algorithm complexity, so the naive algorithm has complexity O(l) as it has to carry out l -1 comparisons for l pixel long SE. Such complexity suggests that naive algorithm is inefficient for any large SEs. Pecht [START_REF] Pecht | Speeding-up successive minkowski operations with bitplane computers[END_REF] proposed a method to decrease complexity based on logarithmic SE decomposition, thereby achieving O(⌈log 2 (l)⌉) complexity.

The first 1-D algorithm that reduced complexity to the constant is often referred to as HGW (it was published simultaneously in two papers: van Herk [START_REF] Van Herk | A fast algorithm for local minimum and maximum filters on rectangular and octagonal kernels[END_REF], and Gil and Werman [START_REF] Gil | Computing 2-D min, median, and max filters[END_REF]). The computation complexity is constant, i.e., of O [START_REF] Adams | Radial decomposition of discs and spheres[END_REF], which means the upper bound of computation time is independent of the SE size. The HGW algorithm uses two buffers, which are filled by the forward, and backward, propagation of local maxima, respectively. Both buffers are then merged into the result. The major drawback of this algorithm is the requirement of two data scans: forward and reverse (so-called causal and anti-causal), which need temporary input data storing.

Lemonnier and Klein [START_REF] Lemonnier | Fast dilation by large 1D structuring elements[END_REF] propose another O(1) algorithm that also identifies local extrema and propagates their values. Again, the limiting forward and reverse scans are needed for every non-causal SE. Lemire [START_REF] Lemire | Streaming maximum-minimum filter using no more than three comparisons per element[END_REF] proposes a fast stream-processing algorithm O(1) for causal line SEs. It replaces the line buffers of the previous algorithms by more dedicated memory structure-double-ended FIFO (queue). The author proposed that only locally monotonous signal is necessary for a given operation. The double-ended queue serves well for such a purpose; however, the algorithm works with causal SEs only.

This downside was solved later in Dokládal and Dokládalová [START_REF] Dokládal | Computationally efficient, onepass algorithm for morphological filters[END_REF] who proposed another queue-based algorithm (see Section 4 for further description of his algorithm). The advantages of these queue-based algorithms are low memory requirements, zero latency, and strictly sequential access to data, the paper presents also a comparative study of these features (Table 1).

ii) the decomposition into atomic shapes allows decomposition of large SE into a sequence of small SE, see Sec. 3.1 for more information on SE composition. Xu [START_REF] Xu | Decomposition of convex polygonal morphological structuring elements into neighborhood subsets[END_REF] claims that any 8-convex polygon (convex on 8-connectivity grid, hence 8-convex) is decomposable into a class of 13 nontrivial indecomposable convex polygonal SEs Q 1 -Q 13 shown in Fig. 1 (a). Normand [START_REF] Normand | Convex structuring element decomposition for single scan binary mathematical morphology[END_REF] reduces the class of shapes to only four 2-pixel SEs, see Fig. 1 (b), by allowing the union operator to take place in SE decomposition. For instance, Q 12 by Xu is obtained as

Q 1 Q 2 Q 3 Q 4 Q 5 Q 6 Q 7 E 1 E 2 E 3 E 4 Q 8 Q 9 Q 10 Q 11 Q 12 Q 13 (a) Xu (b) Normand
(E 3 ⊕ E 3 ) ∪ E 4 by Normand.
The drawback of this approach is that large SEs need a long sequence of atomic operations, and may need iterate several times over the image.

For large, rotation-symmetric polygons, the problem has been solved in Soille et al. [START_REF] Soille | Recursive implementation of erosions and dilations along discrete lines at arbitrary angles[END_REF], who propose an approximation of polygons by a set of line SEs rotated by different angles. The complete dilation by a polygon still requires several iterations over the image. However, the number of iterations does not depend any longer on the size but rather on the shape. It needs three iterations for a hexagon, and four for an octogon. Each diagonal line is computed by the fast 1-D HGW algorithm oriented by the desired angle. The orientation of the SE is achieved through image partition into discrete lines (parallel, with no overlap), along which the image is processed.

More complex SEs can also be computed directly by dedicated algorithms, see Van Droogenbroeck and Buckley [START_REF] Van Droogenbroeck | Fast computation of morphological operations with arbitrary structuring elements[END_REF].

Hardware implementations

The hardware implementations of mathematical morphology are often called dataflow architectures in literature as they process an image in stream. They can be classified into three groups: (i) 3 × 3 neighborhood processors, (ii) partialresult reuse (PRR), and (iii) implementing efficient 1-D algorithm with O(1).

One of the first 3 × 3 neighborhood architectures was the texture analyzer [START_REF] Klein | The texture analyser[END_REF]. It was optimized for linear and rectangular SE by decomposition into line segments. In [START_REF] Klein | Pimm1, an image processing ASIC based on mathematical morphology[END_REF] the authors devised PIMM1 (Processeur Intégré de Morphologie Mathématique) ASIC that contains one numerical unit for gray-scale images and 8 binary units. However, the mechanism of computation was not communicated. Ruetz and Brodersen [START_REF] Ruetz | Architectures and design techniques for real-time image-processing IC's. Solid-State Circuits[END_REF] proposed another ASIC chip that separates the supported 3 × 3 SE into line segments as well. Then only 4 diadic comparisons are necessary to compute the 3 × 3 SE.

More recently, Velten and Kummert [START_REF] Velten | Implementation of a high-performance hardware architecture for binary morphological image processing operations[END_REF] propose another delay-line based architecture for binary images supporting arbitrarily shaped 3 × 3 SEs. The computation of dilation is realized by OR gates (topology was not communicated, probably a tree of diadic OR gates) achieving good performance, which was further improved by spatial parallelism.

Clienti et al. [START_REF] Ch | A system on chip dedicated to pipeline neighborhood processing for mathematical morphology[END_REF] proposes a highly parallel morphological System-on-Chip. It is a set of neighborhood processors PoC optimized for arbitrarily shaped 3×3 SE interconnected in a partially configurable pipeline. The dilation itself is carried out by a tree of diadic max operators, which is pipelined for better performance.

All previous architectures use the naive method to compute the morphological operations. In order to decrease a number of comparisons for large SEs, the PRR method (name proposed in [START_REF] Chien | Partial-result-reuse architecture and its design technique for morphological operations with flat structuring elements[END_REF]) takes a partial result of a morphological operation by some neighborhood B 1 in an early stage, delays it, and reuses it later in computation by some other neighborhood B 2 obtaining thus even other, larger B 3 .

One of the first PRR architectures for 1-D dilation was proposed in Pitas [START_REF] Pitas | Fast algorithms for running ordering and max/min calculation[END_REF] and improved in Coltuc and Pitas [START_REF] Coltuc | On fast running max-min filtering[END_REF]. The principle is based on so-called logarithmic SE decomposition. Even though it reduces a number of comparisons from naive l -1 to ⌈log 2 (l)⌉, it restricts a family of possible SE shapes to rectangles only.

The family of SE shapes has been enriched by Chien [START_REF] Chien | Partial-result-reuse architecture and its design technique for morphological operations with flat structuring elements[END_REF]. The authors presented more general concept of PRR that builds the desired SE by a set of distinct Xu elementary neighborhoods, see Fig. 1 (a). The decomposition process is computed by a dedicated algorithm. As a result, it supports arbitrary 8-convex polygon at the cost of additional comparisons. In [START_REF] Chien | Partial-result-reuse architecture and its design technique for morphological operations with flat structuring elements[END_REF] the PRR method was implemented as an ASIC chip supporting 5-diameter disk SE. Despite decent performance, the chip lacks possibility to control the shape of the SE.

A similar approach has been published by Déforges et al. [START_REF] Déforges | Fast recursive grayscale morphology operators: from the algorithm to the pipeline architecture[END_REF]. Based on Normand SE decomposition [START_REF] Normand | Convex structuring element decomposition for single scan binary mathematical morphology[END_REF] (a SE is decomposed into a number of causal 2-pixel SEs, which are applied in sequence or in parallel, see Fig. 1 (b)) and combined with a stream implementation, the authors propose a methodology for pipeline architecture design supporting arbitrary 8-convex SEs. The practical limitation comes from the need to create a long pipeline of atomic modules for large SE. Even thought it is not specified in the paper, such pipeline seems to be dedicated to the given SE shape and size.

Regarding the third group, implementations based on efficient algorithms, there are only two proposals in literature. Clienti et al. [START_REF] Ch | An efficient hardware architecture without line memories for morphological image processing[END_REF] implemented 1-D HGW and Lemonnier algorithms. In order to avoid a hardware-expensive reverse scan, they devised a ping-pong mirroring buffers that pro-vide reverse-scanned data with minimal latency and memory requirements. The major drawback of the two implementations is incapability of supporting vertical, or 2-D, SE along with horizontal scan data reading.

Prior to this paper, Bartovsky et al. [START_REF] Bartovský | Parallel implementation of sequential morphological filters[END_REF] reported an efficient parallel design based on the 1-D dilation algorithm [START_REF] Dokládal | Computationally efficient, onepass algorithm for morphological filters[END_REF]. However, it supports rectangular SEs only.

From the paragraphs above we can see that there are few hardware architectures capable of supporting polygonal SEs, and none of them is optimized for polygons. These architectures are usually suitable for small SEs but lose hardware resources efficiency for large SEs. In this paper we propose an architecture primarily dedicated to large polygonal SEs using an efficient algorithm.

Novelty

Most previous implementations can efficiently compute a single dilation or erosion with large or even arbitrarilyshaped SE using the decomposition into simpler, atomic operations. Provided the atomic operations concatenate, the computation can be efficiently implemented in a pipe. Obviously, more complex or larger shapes will require longer pipes. If the size is a priory unknown (e.g. function of a parameter) the complete pipe cannot be instantiated and one will need to iterate several times over the image, and store intermediate results in the memory.

Moreover, the morphological dilation is never (or rarely) used alone, but rather as a basic brick in filters. Also in a typical application consisting of several stages -for example: i) denoising, ii) object detection, and iii) measures -the dilation can even be used tens of times, with various SEs, from small (denoising) to large ones (object detection). Lastly, in geodesic operations (e.g. morphological reconstruction) the dilation can be used a variable number of times. Even though the resources of a long pipe are not excessive, such a realization lacks polyvalence and flexibility. It will only fit the targeted operator or application.

In this proposition, the atomic operation is the dilation by a large polygon (similarly to [START_REF] Bartovský | Parallel implementation of sequential morphological filters[END_REF] with rectangles). Such atomic operation situates at a higher level of abstraction which allowing simpler decompositions, and consequently shorter pipes.

This operator has been embedded in a programmable block (dilation/erosion and the SE size and shape). It also retains all beneficial features of the inherent algorithm [START_REF] Dokládal | Computationally efficient, onepass algorithm for morphological filters[END_REF], e.g. the sequential access to data, zero latency and low memory requirements. These advantages are especially beneficial in more challenging applications, such as ASF filters, that can be computed in one or a few image scans.

Basic Notions

Let δ B , ε B : Z 2 → R be a dilation and an erosion on grayscale images, parameterized by a structuring element B, as-sumed to be flat (i.e., B⊂Z 2 ) and translation-invariant, defined as

δ B (f ) = b∈B f b ; ε B (f ) = b∈ B f b (1)
The

ϕ B (f ) = ε B [δ B (f )] ; γ B (f ) = δ B [ε B (f )] (2)
Furthermore, the concatenation of the closing and opening forms sequential filters, e.g., ASF. The λ-order ASF (referred to as ASF λ ) is composed of the sequence of λ closings and λ openings with a progressively increasing SE. It starts with either the closing or opening

ASF λ = ϕ λ γ λ ϕ λ-1 γ λ-1 . . . ϕ 1 γ 1 (3)
ASF λ = γ λ ϕ λ γ λ-1 ϕ λ-1 . . . γ 1 ϕ 1 (4)
Then i.e. the ASF λ starting by closing, can be written as

ASF λ = δ B λ ε B λ ε B λ δ B λ δ B λ-1 ε B λ-1 . . . ε B1 δ B1 . (5)
The initial number of morphological operators 4λ can be reduced using associativity of dilations and erosions. Hence, every two consecutive dilations or erosions may be merged into one to obtain only 2λ + 1 operators, such as

ASF λ = δ B λ ε B λ ⊕B λ δ B λ ⊕B λ-1 . . . ε B1⊕B1 δ B1 . (6)

SE Decomposition

The separability of n-D morphological dilation into lower dimensions is a well-known property. The decomposed dilations are then applied in a sequence according to the following equation

δ R (f ) = δ H⊕V (f ) = δ H (δ V (f )) (7)
where R denotes a rectangle, H and V horizontal and vertical line segments, respectively. This decomposition applies, in general, to convex shapes. In order to suppress the angular anisotropy of rectangles (note the difference between a side length and a diagonal length), one prefers using circles. Regarding the implementation aspects, circles are often approximated by regular polygons (all sides have the same length) that are easily decomposable, originally described in [START_REF] Adams | Radial decomposition of discs and spheres[END_REF][START_REF] Xu | Decomposition of convex polygonal morphological structuring elements into neighborhood subsets[END_REF].

A 2n-top (n ∈ N) regular polygon SE P 2n can be decomposed into a set of n line SEs L αi [START_REF] Coltuc | On fast running max-min filtering[END_REF] oriented at angle α i , such as

P 2n = L α1 ⊕ • • • ⊕ L αn n times
α i = (i -1) 180 • n [ • ] ; i ∈ N, i ≤ n (9)
The length of all L αi is equal to the side of the desired polygon and can be computed from the circumcircle radius R as

L αi = 2R sin 180 • 2n (10)
For example, a hexagon can be obtained by three Hence, from ( 7) and ( 8) a 2-D dilation by a 2n-top polygon δ P2n of some function f : R 2 → R can be obtained by n consecutive 1-D dilations δ Lαi by line segments oriented by

L αi ori- ented in α i = { 0 • , 60 • , 120 • } on a 6-
α i δ P2n (f ) = δ Lα1 ( . . . δ Lαn n times (f )). ( 11
)
The aforementioned decomposition holds true for the unbounded support Z 2 . However, when using real images with a bounded support

D ⊂ Z 2 , D = [1..M ] × [1.
.N ], decomposition boundary effects appear if at least one L αi = {0 • , 90 • } is used. The cause is that the Minkowski addition of all decomposed line segments in [START_REF] Coltuc | On fast running max-min filtering[END_REF], which are cropped by image boundaries after every L αi of that concatenation, does not necessarily correspond to P 2n cropped by image boundaries just once as desired. It is expressed by the following expression where D∩ represents intersection with the image support

D D ∩ (L α1 ⊕ . . . ⊕ L αn ) = D ∩ (L αn ⊕ . . . . . . ⊕ D ∩ (L α2 ⊕ D ∩ (L α1 ))). ( 12
)
The illustrative example of such boundary effects with a hexagonal SE is depicted in Fig. 3. We can see that the composition α 1 ⊕ α 2 is incomplete compared to the desired one in Fig. 2; a small part of the SE is missing. It holds true even for the entire hexagon, the composition α 1 ⊕ α 2 ⊕ α 3 is also incomplete. It is caused by the right boundary cropping not only the final P 2n , but also all intermediate results. The cropped values are later missing to form an appropriate polygon section. This issue is solved by adding a padding to the image. The section of P 2n contained inside the image support is then complete, the missing part of P 2n is located in the padded area. The added padding contains recessive values, i.e., values that do not affect the computation of a particular morphological operator (for f :D → V , ∧V for dilation, ∨V for erosion). The thickness of the padding is different in the horizontal and vertical direction and is determined by the size of oblique segments, particularly by the half of vertical and horizontal projection

0 0 0 N N M 0 M missing part missing part (a) α1 ⊕ α2 incomplete (b) α1 ⊕ α2 ⊕ α3 incomplete
B H = L αi cos (α 2 ) /2 [pixels] ( 13
) B V = L αi sin (α 2 ) /2 [pixels]. ( 14
)

Algorithm Description

This section explains the algorithmic principles involved in this paper. First, we expose the 1-D algorithm used for the dilation by line segments with arbitrary orientation. Next, we show how to combine these 1-D computations in order to obtain polygons running in stream. The issue of boundary handling is addressed afterwards.

1-D Dilation Algorithm

The implementation of (1) consists of searching the extremum of f within the scope of SE B

[δ B (f )](x) = max b∈B [f (x -b)] (15)
[ε B (f )](x) = min b∈B [f (x + b)] (16)
The algorithm used below is based on property [START_REF] Dokládal | Computationally efficient, onepass algorithm for morphological filters[END_REF] that for some B(x) (which contains its origin) the computation of the dilation δ B f (x) needs only those values of f (x i ) that can "be seen" from x when looking over the topographic profile of f , see Fig. 4. The values shadowed by the mountainsthat cannot be maxima -are immediately excluded from the computation. From the implementation point of view, assuming a sequential access to the input data f , the dilation δ B f (x) depends on points read after x. We say that B is non causal. One can transform a non causal SE to a causal one by utilizing the property that dilation commutes with translation

δ B+t f (x) = δ B f (x -t), ∀t ∈ D (17)
These two principles are used by Alg. 1. For each pixel of some input signal f : Z → R, the algorithm reads one pixel F = f (rp) at the so-called reading position rp and writes back one result pixel dF = δ B f (wp) at the current writing position wp, such as rp > wp. The wp coordinate coincides with the origin of the SE, rp conforms to the most recent input pixel of B. Indeed, the reading position rp is its right-hand side end, which conforms to the intuitive necessity of having read all the samples covered by the SE before computing the dilation.

Alg. 1 is to be called from an outer loop iterating over the writing position in δ B f , such as while wp<N. The writing position wp is to be incremented whenever Alg. 1 outputs a valid value. We give below the pseudocode, for detailed description see [START_REF] Bartovský | Parallel implementation of sequential morphological filters[END_REF].

Stream-Preserving Decomposition of Polygons

Alg. 1 can be used to compute the dilation by L αi segments in a stream. Its properties make it suitable for composing concatenated operators, namely the sequential access to input and output data, and minimal latency. Therefore, when

O((N + 2B H )(M + 2B V )) (18)
for an N × M image, and B H , B V padding sizes. Provided that of B H ≪ N and B V ≪ M , it reaches the complexity of rectangles O(N M ), see [START_REF] Bartovský | Parallel implementation of sequential morphological filters[END_REF].

Discretely Inclined 1-D Segments

The oblique segments included in a hexagon and an octagon, i.e., L αi , α i = 45, 60, 120, 135

col shift = (col + line/k i ) mod (N + 2B H ), (19)
where the inclination deviation from the vertical direction line/k i is called offset and changes only with a new image line.

Hardware Implementation

In this section, we present hardware implementation (called line unit LU) of Alg. 1 for dilation by L αi with emphasis on the inclined segment computation. Then, we chain several elementary LU units into a pipeline to form the polygonal processing unit PU. Finally we propose the parallel polygonal unit PPU.

1-D Algorithm Implementation

Alg. 1 along with the col shift addressing feature is seen as a simple Mealy finite state machine (FSM). This FSM controls all algorithm operations over the queue, rp and wp pointers etc. The state diagram (in Fig. 6) of the algorithm behavior consists especially of 2 main states {S1,S2} and one auxiliary state EOL. The basic operation of the direct algorithm implementation can be found in [START_REF] Bartovský | Stream implementation of serial morphological filters with approximated polygons[END_REF][START_REF] Bartovský | Pipeline architecture for compound morphological operators[END_REF]. The auxiliary state EOL is entered only at the end of every image line. Its main purpose is to update the offset value, to determine the shifted column addressing. The generation of the necessary inclination is extremely easy since it requires only elementary operations like incrementing, decrementing or stalling.

1-D Line Unit Architecture

The architecture of the LU unit capable of dilation by different line segments is shown in Fig. 7. The basic description of the preceding version supporting only horizontal and vertical orientation can be found in [START_REF] Bartovský | Pipeline architecture for compound morphological operators[END_REF]. Several modifications have been applied to the former version to allow inclined L αi . We have mainly added the Slope control unit that is highlighted in Fig. 7.

The LU comprises two parts: an FSM part and a memory part containing a collection of double-ended queues. FSM manages the whole computing procedure and temporarily stores values in the memory part. The memory instantiates one queue in the case of horizontal segment, N queues in the vertical case (N is the image width), or N + 2B H queues in the oblique case. Input and output ports are multiplexed; hence a multiplexor select signal can easily address one queue to work with. The shifted column address (col shift ) is used as the select signal.

The processing of one pixel proceeds as follows: In the beginning of S1, the last queued pixel is invoked by the Back() operation from the queue and fetched to Comparator 1 where it is compared with the current sample, and dequeued if necessary (code lines 1-2). The current pixel is simply extended with the value of position counter 1 and enqueued (line 3).

The S2 invokes the oldest queued pair {F, rp} by the Front() operation. This read pixel is a correct result if the set of output conditions (code line 6) is fulfilled. The deleting of an outdated value is managed by comparing the stored position value with the current one in Comparator 2.

The purpose of Slope control is to select the corresponding queue memory which is currently used by Alg. 1. The queues are addressed by the Col shift counter, which is incremented with every pixel of the input image and reset at the end of the image line. The initial reset value of the col shift counter is offset (see Section 4.3). The offset is updated at the end of every image line (state EOL); its value is incremented or decremented either every line or every other line according to k i .

Polygon Unit Architecture

The LU units described above can be arranged in a sequence to form a 2-D Polygon Unit (PU). The overall architecture of the PU unit (see Fig. 8) is composed of three differentpurpose parts: computation part, controller, and padding part.

The computation part mainly contains four LUs for distinct L αi orientations. There are the horizontal unit (α 1 = 0 • ), the first inclined unit (α 2 = 45 • or 60 • ), the second inclined unit (α 3 = 135 • or 120 • ), and the vertical unit (α 4 = 90 • ) connected in a simple pipeline; the output of each unit is read by the successive unit which processes the image by further L αi . The computation part is able to operate either with a hexagon or octagon SE. In the case of the hexagon SE, the vertical unit is bypassed.

Note that the output of every computation unit is an intermediate result image, which can be brought out for another purpose, e.g., a multi-scale analysis descriptor. Then the dilation by line, rectangle, and octagon SEs (all centered) can be obtained during a single image scan (considering units reordering). Only the Remove padding block is to be copied several times for each output data stream.

According to the boundary effects mentioned in Section 4, the inclined units need padding to extend the original image before the processing. The padding is removed after the last 1-D unit. It is carried out by a pair of dual padding blocks at the beginning and the end of the computation part.

The controller ensures the correct global system behavior. It accepts the SE diameter and the shape select signal, then it determines the particular SE sizes for every LU and padding from them, and initiates the computation. The entire set of parameters, i.e., the image width and height, all SE features (size and shape), and the morphological function select, is run-time programmable at the beginning of the frame. These parameters are run-time programmable within the upper bound specified during the synthesis.

For instance, consider the SE size upper bound a 91×91 bounding box, and octagon-capable architecture. This means, the architecture has four LU; each LU supporting at most l=31 pixels segments. During the operation -at the beginning of a frame -the SE can programmed to either of the following: a line up to 31 pixel long, a rectangle up to 31 × 31 pixels, or an octagon up to 91 × 91.

To enable processing a uniform input stream, one needs to handle unequal processing rates of LUs. It is caused by variable algorithm latency to compute a dilation for one pixel. Therefore, the balancing FIFO memories are inserted in front of each 1-D unit, and to the input and output ports.

The depth of input and output FIFOs depends on the timing of input data stream (possibility of stalling, synchronization, etc.).

Memory Requirements

The most significant memory demand is made by the set of queues. Although the algorithm works with separated queues, the queues within each LU are merged into a single dual-port memory, mapped side by side in a linear memory space. Every queue has a related pair of front and back pointers which must be retained throughout the entire computation process in the pointer memory. This approach leads to more efficient implementation. The LUs have the following memory requirements (considering N×M image including padding, L αi with bounding boxes W x ×H x , and bpp bits per pixel):

M hor =W H (bpp + ⌈log 2 (W H -1)⌉) [bits] ( 20
) M ver =N ((H V -1)(bpp + ⌈log 2 (H V -1)⌉) + 2⌈log 2 (H V -1)⌉) [bits] (21)
M slope =(N + W S )((H S -1)(bpp + ⌈log 2 (H S -1)⌉) + 2⌈log 2 (H S -1)⌉) [bits] (22)
Example: Consider a dilation of 8-bit, SVGA image (i.e., 800×600=N ×M ) by a hexagon with radius 41 px. Such a SE is decomposed into horizontal SE 21 px wide, and 2 slope SE each 11 px wide and 19 px tall (hexagon SE bounding box is 41×37 px).

The computation memory (the queues) requires (20-22)

M hor =21(8 + 5) = 273 [bits] M slope =(811 + 11)((19 -1)(8 + 5) + 10) = =200 ′ 568 [bits]
resulting in total consumption of M all = M hor + 2M slope ∼ = 392 kbits for the 2-D dilation by hexagon. This is far below the mere size of the image itself M image = 800 × 600 × 8bpp ∼ = 3.66 Mbits which does not need to be stored at any moment.

Parallel Implementation

This section describes the Parallel Polygon Unit (PPU) that aims at increasing the computational performance while maintaining as much as possible the beneficial streaming properties of the proposed algorithm.

Partition of the Image

The parallelism is obtained by use of concurrently working units that simultaneously process different parts of the image (spatial parallelism). The number of instantiated units defines the parallelism degree (PD). Since the processing runs in stream, we propose a solution that transforms the input stream into a set of PD streams in a way to minimize the waiting-for-data periods of all units. For the sake of clarity, we use PD=2 in the description hereafter. A similar method has proven to be useful in [START_REF] Bartovský | Parallel implementation of sequential morphological filters[END_REF].

The partition of the input image is twofold, see Fig. 9: an interleaved line-by-line partition for the horizontal α 1 units, and vertical stripes for the vertical and inclined α 2 , α 3 , α 4 units. The final image partition of 2-D image is the intersection of both. Intuitively, the streams have to be transformed from one type to the other between α 1 -α 2 , and α 4 -output in the PU. The transformation is done by simple circular stream switching when a partition edge is encountered. With the beginning of the image, it starts with the H1•V1 segment on the first line. When the end of this segment is reached, the streams are switched so that segments H1•V2 (1 st line) and H2•V1 (2 nd line) are processed at the same time. Later, it proceeds to segments H2•V2 (2 nd line) and H1•V1 (3 rd line) and so forth. In general PD segments located on a backward diagonal run simultaneously throughout the image (note that the streams are mutually delayed by N/PD pixels).

Processing the partition segments separately introduces undesired border effects on each partition edge. A common solution -similar to padding at image borders -is to introduce an overlap. Contrary to the padding that adds recessive values, the overlap extends a partition by a portion of the neighboring partition. The width of the overlap depends on the size of the SE, and is equal to the width of the horizontal padding B H . Intuitively, the overlap introduces redundant computation, and slightly degrades the performance and minimal latency.

Parallel architecture

At this point, all the previously mentioned principles are brought together to form the Parallel Polygon Unit (PPU). The PPU (see Fig. 10) is scalable with respect to PD, the number of parallel streams it can process at the time. Each stream needs one pipeline of four LUs (α i , i = 1..4, just like the PU), two add overlap blocks in front of inclined LUs, two remove overlap blocks behind inclined LUs, add padding at the front end, and remove padding at the back end. The PPU also contains a pair of switches to transform the streams from one type to the other, and a controller (ommited in Fig. 10).

Figure 11 shows the introduction of the overlap in the course of the i-th image line. As we know, this line is split into two streams. The streams are labelled I1, I2 before the addition and O1, O2 after (refer also to Fig. 10). The entire I1 stream plus B H pixels of I2 form O1 output stream with overlap, and last B H pixels of I1 and the whole I2 stream form O2.

During the overlap sections, either I1 or I2 stream is mapped to both output streams at the same time. This data duplication does not temporarily allow for parallel processing of both streams and may result in stalling of either stream. However, the effect of the overlap is negligible as long as B H ≪ N .

Two important properties are to be noted: (i) input and output streams are mutually delayed by N/PD (ensured by stream switching); (ii) several PPUs can be chained into a pipe. The schematic of some application, e.g., ASF, may look like in Fig. 12. At the front end there is an input buffer transforming the input stream (which is PD-times faster than each of PD processing streams) into PD processing streams Hi (i = 1..PD). The transformation only needs i-th image line to be stored in {i mod PD}-th line buffer. In this manner, the processing streams are properly delayed by N/PD pixels. The output buffer transforms PD processing streams into one fast stream in the opposite way. One can place as many PPUs as desired between these two buffers in a pipeline or other topology.

The PPU involves the following limitations on the programmability: the image size is set before synthesis, the padding sizes B H , B V are computed for the maximal SE, specified before the synthesis. The reason is that handling the varying SE and image sizes would introduce an unreasonable hardware overhead of image partition, padding, and overlap features. The SE features (size and shape) remain fully programmable as in the non-parallel case, see Section 5.3.

Experimental results

Hereafter we discuss the results of the proposed implementation. First, we discuss the results of a single 2-D PU and PPU unit, followed by their performance in an application, an ASF filter. We conclude by comparising them with other architectures.

The proposed PU and PPU architecture has been implemented in VHDL and targeted to the Xilinx FPGA Virtex-6 device (XC6VLX240T). The ultimate specification conforms to the following: 8-bit gray-scale images of size up to 1080p (1920×1080 px), height of L αi up to 31 px (thus, hexagon SE up to 61 px, and octagon SE up to 91 px), and support of uniform stream processing. Notice that all three previous factors affect the memory requirements that (in contrast to the PC), have significant influence on the clock frequency. Our specification implies the clock frequency of 100 MHz.

The timing with respect to (shortly w.r.t.) the image size and the size of the SE (Table 2 and3) have been evaluated on a natural photo image. We report several measures. 1) The pixel rate gives the average number of clock ticks to process one pixel. It is given by the overall number of clock ticks divided by the image size. One can see that the rate is almost constant w.r.t. both the size of the image and the size of the SE. The slight variation is caused by the size of the SE which affects the size of the padding and overlap, increasing the number of effectively processed pixels, see [START_REF]Mathematical Morphology: From Theory to Applications[END_REF].

2) Latency is expressed in a number of image lines. Note that it is strictly half the SE size. This is a further irreducible factor corresponding to the dependency of the output on the input. This corresponds to the half-height of the SE that needs time to have read enough data to compute the dilation.

3) The last measure is the throughput in terms of the number of frames per second (FPS). The ultimate result we obtain is 105 fps for the 1080p resolution, allowing the 100Hz 1080 FullHD TV standard to be processed in real time.

The speed-up PPU vs PU measures the acceleration obtained from the parallelization. The difference from the ideal upper limit (PD=6) is due to the overlap. With increasing image size the acceleration converges towards 6 because of the SE size (and consequently the overlap) becomes negligible with regard to the image size.

Table 4 outlines efficiency of the scalability (that is the parallelism degree PD) in terms of the FPS and speed-up. One can see that the real speed-up is somewhat lower than the PD. The difference is due to two factors: (i) the overlap, which demands redundant computation, and (ii) the stream switching that needs inter-stream synchronization which may introduce wait cycles.

Table 5 reveals the cost of parallelization on FPGA resources in terms of registers, LUTs, and BRAMs of the PPU and the pair of input and output buffer as shown in Fig. 12.

Alternating Sequential Filter

The ASF filter is an essential method of morphological filtering, see example Fig. 13. Since the ASF is applied as a sequence of alternating dilations and erosions with a changing SE, it can be advantageously implemented by chaining instances of the proposed architecture into a pipeline structure. The output of each operator is immediately processed by a subsequent operator to achieve the following beneficial properties: (i) all the operators are being applied in parallel (temporal parallelism), (ii) the image is filtered with minimal latency inferred by the Minkowski addition of all SEs. Table 7 illustrates the performance of ASF λ in terms of the experimentally achieved FPS and the inferred latency. Note that the frame rate of the whole ASF with respect to the order λ since larger SE implies larger padding and overlap. However, the performance of the filters is comparable with the rate of a single unit in Table 2. The implementation and performance comparison of our architecture with the others is outlined in Table 6. At first, we take into account single 2-D units only. Clienti [START_REF] Ch | A system on chip dedicated to pipeline neighborhood processing for mathematical morphology[END_REF] yields a high throughput for an elementary SE 3×3. The Chien [START_REF] Chien | Partial-result-reuse architecture and its design technique for morphological operations with flat structuring elements[END_REF] ASIC chip achieves a reasonable throughput with a small 5×5 diamond SE. Both architectures use homothecy to obtain larger SEs. On the other hand, one Déforges [START_REF] Déforges | Fast recursive grayscale morphology operators: from the algorithm to the pipeline architecture[END_REF] unit supports various 8-convex SEs in one scan. As mentioned in the state of the art, the programmability of the modules, namely the possibility to control the SE shape after the synthesis is not clear. Smaller throughput is probably caused by usage of a less powerfull device than the rest of implementations.

Lets take as an application example of a compound morphological operator, consider ASF 6 = δ 13×13 ε 25×25 . . . ε 5×5 δ 3×3 that consists of 13 morphology operations. One Clienti's system instantiates 16 elementary 3×3 processing units. Hence, it will require 6 image scans (the entire image must be stored in the memory). Chien also uses the homothecy, therefore, as many as 45 scans are to be done. In the case of Déforges, neither FPGA occupation with respect to the the size of SE nor possibility of using multiple instances in a single chip was communicated. We consider that only one unit fits the FPGA, so 13 image scans are needed. Of course, if several units fit the FPGA surface, it will reduce proportionally the number of scans. This is true for all streaming architectures.

Obviously, between two consecutive scans the data are read/written from/into the memory that degrades performance and significantly increases latency to orders of several image scans. The dense memory traffic might also overwhelm the data bus.

From the estimated performance results for the ASF 6 in Table 6 we observe that the high use of homothecy tends to increase the number of necessary image scans. Indeed, all Clienti, Chien, and Déforges (a) whose solutions are efficient for small SE sizes and short concatenations become more or less penalized for longer concatenations; their performance drop down with the higher numbers of necessary image scans. On the other hand, Déforges (b) and our work, which does not need more than one image scan, attain the high performance for ASF 6 comparable to the performance of a single 2-D unit.

These features allow a temporal-parallel execution of all atomic operators that is essential to obtain the real-time performance for high demanding applications. In addition, the low memory requirements facilitate embedding several instances of proposed computation units into a single FPGA circuit.

Conclusions

It is widely known that the processing data in stream allows to reduce latency, memory consumption and increases the system throughput. Until recently, computing morphological dilations or erosions in stream was only possible for small, limited neighborhoods [START_REF] Chien | Partial-result-reuse architecture and its design technique for morphological operations with flat structuring elements[END_REF][START_REF] Ch | A system on chip dedicated to pipeline neighborhood processing for mathematical morphology[END_REF][START_REF] Déforges | Fast recursive grayscale morphology operators: from the algorithm to the pipeline architecture[END_REF], or large rectangles [START_REF] Bartovský | Parallel implementation of sequential morphological filters[END_REF]. Dilations by large polygons were computed iteratively, by using the homothecy. This required an external memory for intermediate data, limited the flexibility, and drastically increased system latency. This paper opens the possibility of stream execution to morphological dilation with large polygons. Although the decomposition of polygons into the Minkowski addition of inclined lines has been known for years [START_REF] Adams | Radial decomposition of discs and spheres[END_REF], we bring several suggestions that-combined together-allow the execution in stream.

We show how to implement dilation by inclined linear segments with sequential access to input and output data. We show how to handle border effects, and recall (since this is less known) that it requires large padding. Furthermore, we show how to partition an image to introduce efficient spatial parallelism while maintaining sequential access to data at all levels. This avoids increasing the system clock by dividing a fast data stream into several slower streams to process at a slower rate. We show how to efficiently handle the border effects on the partition.

The proposed polygon decomposition uses sequential access to both input and output data. This allows for temporal parallelism, where in concatenations like . . . δεδ . . . all these operators run simultaneously on the time-delayed data. We attain a very low (nearly optimal) latency, which has beneficial impacts on memory consumption. No external memory is used even for large SEs and large images. All these aspects brought together allow for a considerable data throughput for sequential morphological filters. We have implemented a programmable IP block, usable in industrial systems running under heavy timing constraints satisfying up to the 100Hz 1080p FullHD TV requirements.

Fig. 1

 1 Fig. 1 Classes of elementary SEs. Any 8-convex polygon SE is decomposable into either: (a) Xu class, or (b) Normand class while using union along with ⊕.

  hat denotes the transposition of the SE, equal to the set reflection B = {x | -x ∈ B}, and f b denotes the translation of the function f by some scalar b. The SE B is equipped with an origin x ∈ B. The basic concatenation of the dilation and erosion forms other morphological operators. The closing and opening on gray-scale images, ϕ B , γ B : Z 2 → R, parameterized by a structuring element B, are defined as

Fig. 2

 2 Fig. 2 Polygon SE composition of line SEs. (a) hexagon is composed of 3 segments, (b) octagon is composed of 4 segments. ⊕ operator stands for the Minkowski addition; αi stands for Lα i .

Fig. 3

 3 Fig.[START_REF] Bartovský | Stream implementation of serial morphological filters with approximated polygons[END_REF] Polygon SE composition without padding. The desired SEs presented in Fig.2are incomplete, a small triangle is missing.

Fig. 4

 4 Fig.[START_REF] Bartovský | Pipeline architecture for compound morphological operators[END_REF] Computing the dilation δBf (x): Values in valleys shadowed by mountains when looking from x over the topographic relief of f are useless.

Algorithm 1 : 2 Q 5 Q

 125 dF ← 1D DIL (rp, wp, F, SE1, SE2, N) Input: F -input signal sample f (rp); rp, wp -reading/writing position; SE1, SE2 -SE size towards left and right; Nlength of the signal Result: dF -dilated signal sample δBf (wp) Data: Q -Queue (first in, first out) 1 while Q.back()[1] ≤ F do .dequeue() ; // Dequeue useless valuess 3 Q.push({F , rp}) ; // Enqueue the current sample 4 if wp -SE1 > Q.front()[2] then .pop() ; // Delete too old value 6 if rp = min (N , wp + SE2) emptythe input image is read in a horizontal raster scan mode, i.e., line by line, and every line from the left to the right, the output of Alg. 1 instance conforms to the very same scan order, delayed by some latency defined by the distance between reading rp and writing wp positions. It allows a direct connection of several Alg. 1 instances in a sequence without any need of coupling elements. The resulting 2-D SE is then obtained with minimal latency, that is as soon as all necessary data have been read. The example of decomposition of a hexagon into three L αi is depicted in Fig.5. The image is sequentially read by horizontal L 0 • at the reading position of the polygon (a). The result of the horizontal segment is immediately provided as an input to the first oblique L 60 • at (b) so that the reading position of L 60 • coincides with the writing position of L 0 • . By the very same rule, the result of the L 60 • is brought as input data to the second oblique L 120 • at (c), the writing position of which is the writing position of the complete polygon (d). The total latency is then defined by distance between the reading (a) and the writing position (d) of the polygon.

Fig. 5

 5 Fig. 5 Stream concatenation of three Lα i into hexagonal SE P ; rp/wp -reading/writing position.

  S1 state manages the dequeuing loop and pushing of a new value, code lines 1-3; the S2 state handles the deletion of outdated values and returns the result, code lines 4-9. (); return (Q.front()[1] );

Fig. 6 Fig. 7

 67 Fig. 6 State diagram of Alg. 1. Conditions of state transitions are typed in bold, output actions are located in gray rectangles.

Fig. 8

 8 Fig.[START_REF] Coltuc | On fast running max-min filtering[END_REF] Overall architecture of the polygonal PU unit. It contains one LU for each δ L α i of (11), control, and padding units.

Fig. 9

 9 Fig. 9 Example of image partition for PD=2: line by line for horizontal orientation; vertical stripes for non-horizontal orientation.

Fig. 11 Fig. 12

 1112 Fig. 11 Addition of overlap on one image line

Fig. 10

 10 Fig. 10 Overall architecture of the parallel polygonal unit PPU for PD=2. The controller and balancing FIFOs are ommited.

Fig. 13

 13 Fig.[START_REF] Klein | Pimm1, an image processing ASIC based on mathematical morphology[END_REF] Example of ASF filtering. A zoom into original and the ASF 3 filtered "Mountain" image.

Table 1

 1 Comparison of fast 1-D dilation algorithms.

	Algorithm	SE type Compar.	Alg. Data	Working
			per pixel	lat.	mem. mem.
	Naive 1-D	User	l -1	0	N	0
	HGW	Sym	3 -4/l	l	N	2l
	Lemire	Causal	3	0	0	2l
	Lemmonier Sym	NC (O(1)) N	N	N
	Dokladal	User	3	0	0	2l
	Sym = symmetric SE; User = User-defined SE; l = length of a 1-D SE;
	N = line size; G = number of gray levels; NC = not communicated.

  connected grid, and an octagon by four L αi ,α i = { 0 • , 45 • , 90 • , 135

• } using an 8-connected grid, see Fig.

2

.

  • , need appropriate addressing to determine the pixels to process. Note that all inclinations verify α i ≥ 45 • , and the coefficients k i verify k i = tan α i ≥ 1. If we use 8-connectivity for k 45 • ,135 • = ±1, and 6-connectivity for k 60 • ,120 • = ±2, we can very easily generate the pixel adressing -for every inclination -by only modifying the original column index col by an additive constant line/k i such as

Table 2

 2 Timing of PU and PPU w.r.t. image size (SE size = 51 px, PD=6).

	Image Size	VGA SVGA XGA 1080p
	Pixel Rate (PU) [clk/px]	2.61	2.53	2.53	2.44
	Latency [image line]	25	25	25	25
	FPS (PU) [frame/s]	125	82	50	19
	FPS (PPU) [frame/s]	599	406	257	105
	Speed-up PPU vs. PU [-]	4.79	4.94	5.1	5.34

Table 3

 3 Timing of PU w.r.t. SE size (SVGA image)

	SE size [px]	21	31	41	51	61
	Rate [clk/px]	2.42 2.46 2.49 2.53 2.58
	FPS [frame/s]	85	84	83	82	81
	Latency [image line]	10	15	20	25	30

Table 4

 4 Speed-up of PPU w.r.t. PD (SVGA image, SE size = 31 px).

	Parallelism degree PD	2	3	4	5	6
	FPS [frame/s]	162	234	306	376	441
	Speed-up [-]	1.92 2.77 3.62 4.44 5.22

Table 5

 5 FPGA resources w.r.t. PD (SVGA image, SE size = 91 px).

	PD	1	2	3	4	5	6
	Registers (P)PU	787 1,644 2,469 3,215 4,019 4,850
	LUTs (P)PU	2,656 4,831 7,330 9,301 11,540 14,221
	Block RAM (P)PU 39	39	59	42	53	63
	Registers buf	0	251 355 466	590	671
	LUTs buf	0	1,296 1,929 2,545 3,158 3,748

Table 6

 6 Comparison of several FPGA and ASIC architectures concerning morphological dilation and erosion. N , M stand for the image width and height of respective architectures.

			Processing unit		Hardware System	Application Example ASF 6
		Technology	Supported	Throughput fmax Number of Supported Image	FPS	Latency
			SE	[Mpx/s]	[MHz]	units	image	scans [frame/s]	[px]
	Clienti [6]	FPGA	arbitrary 3×3	403	100	16	1024×1024	6	66.7	5N M + 84N
	Chien [5]	ASIC	disc 5×5	190	200	1	720×480	45	12.2	44N M + 84N
	Déforges [10]	FPGA	arbitrary 8-convex	50	50	1	512×512	13	14.7	12N M + 84N
	This paper	FPGA	regular polygon	195	100	13	1024×1024	1	185	84N

Table 7

 7 Timing of ASF λ ; SVGA image size, PD=6.

	Order of ASF λ	1	2	3	4	5	6
	Number of δ, ε	3	5	7	9	11	13
	Size of max. SE [px]	5	9	13	17	21	25
	FPS by PUs [frame/s]	88	87	86	86	86	85
	FPS by PPUs [frame/s] 491 483 466 440 415 387
	Latency [image line]	4	12	24	40	60	84
	7.2 Architecture Comparison