Nicolas Ngan
email: ngann@esiee.fr

Eva Dokladalova
email: e.dokladalova@esiee.fr

Mohamed Akil
email: akilm@esiee.fr

Franc ¸ois Contou-Carrère
email: francois.contou-carrere@sagem.com

Dynamically adaptable architecture for real-time video processing

In this paper, we present a new adaptable ringbased architecture for video processing applications. The proposed architecture allows handling pipelined and parallel organization of computation for multiple video flows. A simplified version with four nodes has been implemented on an FPGA for a video application to show the adaptation mechanism between a pipelined and parallel structure.

I. INTRODUCTION

In addition to the quality of image sensors, the modern portable video equipments need a performant embedded computing system enabling more and more of advanced functionalities such as panoramic, picture-in-picture or image fusion applications. Those advanced video processing chains involve the management of multiple video flows depending on the spatial and temporal complexity of the application. Moreover, those complex applications have to be implemented in a variable but limited surface (depending of the portability of the system) and have to meet variable real-time constraint from 25 to 50 frames per second with an image size up to the 1080p HD resolution (1920x1080).

The last decades of research have shown that the performance and flexibility are determined by the balance between its programmability and reconfigurability levels (granularity) of computing systems [START_REF] Hartenstein | A decade of reconfigurable computing: a visionary retrospective[END_REF]. In parallel to the highly parallelized architectures, lots of linear architectures has proven that dedicated pipelined architectures [START_REF] Clienti | A system on chip dedicated to pipeline neighborhood processing for mathematical morphology[END_REF] [START_REF] Dokládalova | Carvision: SOC architecture for dynamic vision systems from image capture to high level image processing[END_REF] are the most performant and natural structure for video-based system. However, pipelined architecture depth is always limited by the surface allowed on a chip and the only flexibility remains on the programmability of the processing elements. Furthermore, the size of the pipeline has a direct impact on the latency and the only solutions are to either increase the processing elements efficiency or bypass some processing elements to shorten the datapath which results a loss of computing power.

Despite the numerous dedicated architectures for different image processing problems, there are few systems considering the efficient processing multiple video flows. In general, the processing chains are duplicated or video chains are merged at differents steps. This results a very specific system which excludes any major modification of the structure. Thus, in the context of image sensor evolution (resolution, pixel granularity, video acquisition speed) and as a consequence the video processing algorithms modifications, those rigid solutions could not be applied anymore except at a high development cost between each generation of equipment.

In order to deliver required performances and to offer the flexibility, the interconnections of computing ressources has to be chosen very carefully. In the last few years, ring-based architectures become very trendy in audio [START_REF] Creative | X-fi ring architecture[END_REF], graphics [START_REF] Ati | Radeon x1800 memory controller -whitepaper[END_REF][6], image such as [START_REF] Van Der Horst | A real-time stereo smartcam, using fpga, simd and vliw[END_REF] [8] [START_REF] Sassatelli | Systolic ring: A new approach for dynamical reconfigurable architectures[END_REF] and general purpose applications [START_REF] Collette | Symphonie massively parallel computer : Modelling and design[END_REF] [START_REF] Kistler | Cell multiprocessor communication network : Built for speed[END_REF]. These propositions allow to satisfy the performances but, unfortunately, most of those architectures are either unable to manage different parallel input flows or are simply not adapted for embedded systems. Thus, we propose to adapt the datapath for both pipeline and parallel execution according to the video application and the processing element.

In this paper, we present a new adaptable ring-based architecture for video processing applications. The proposed architecture allows handling pipelined and parallel organization of computation for multiple video flows. A simplified version with four nodes has been implemented on an FPGA for a video application to show the adaptation mechanism between a pipelined and parallel structure.

The paper is organized as follows. The Section II presents processing patterns in modern applications and resulting architecture specifications. The architecture proposition is described in Section III and its first FPGA implementation is proposed in Section IV.

II. SPECIFICATIONS FOR EFFICIENT ARCHITECTURE

DESIGN

From the discussion in the Introduction, we can conclude that design of the efficient computing pipeline is the most important aspect to keep performance but it has to be flexible enough to handle the multiple possible situations such as video loop-back, multiple inputs pipeline execution mode or multiple dependencies. In this paper, we call a pipeline a highly optimized datapath used to connect in a series computing ressources in order to process a given data flow. We say that it is full dedicated solution when the length of used pipeline is minimal and the obtained computing latency is minimal.

The analysis of numerous video processing chain has allowed to identify several repetitive data flow patterns, representing the most frequent pipeline implementations. Let T i the tasks to be executed, then the different implementations can be illustrated as follows:

Pipeline: it is the most performant and memory cost efficient method. It allows to store and process directly on the flow. Figure 1 illustrates this execution mode. The arrows represent the data dependencies between tasks. Multiple pipeline: some video applications need different input flows of the same nature to be combinated such as creating a panoramic with several images or incrust a picture in another (picture-in-picture technique). In this execution mode, multiple flows have to be processed and synchronised according to the tasks (Figure 2). Multiple pipeline with parallel tasks: this pattern adds new tasks in the video processing chain and those new tasks could be simply chained with the others (increasing the pipeline) or need to work in parallel of the main video processing flow. They can also create new data dependencies as shown in Figure 3. Frame loopbacks execution mode: the video processing chain might require to buffer several image frames. Let us consider that each task T i is implemented in a Processing Element P E i . The main idea of the adaptation is to morph dynamically between the different data flow patterns (Fig. 1234). The principle is also partially shown in Figure 5. For multi-tasking purposes, this adaptation should load or unload the pipeline from PEs dynamically and use the unloaded ones to realize parallel processing on different data flows.

Those adaptations can only be done physically by adding some additionnal dedicated hardwares (called DF Rs i.e. DataFlow Routers) to create the different connections allowed between PEs. We propose to interface the PEs with DataFlow Routers (DF Rs) which are interconnected on a traditional ring structure. All datapath modifications are realized inside the DataFlow Routers. We call DF R I/O a DF R node equipped with external input and output video flows. The ring topology has the advantage to be relatively simple to implement and each node has only and always two neigbours which means that routers and interconnections can be implemented at a limited cost [START_REF] Saldana | The routability of multiprocessor network topologies in fpgas[END_REF] compared to much more complex mesh structures. The structure can be used in a standard pipeline style computation. In addition, in the case of high area constraint, the ring topology facilitates the modifications of PEs (internal parallelization strategy or increasing the working frequency to keep the video rate).

III. ARCHITECTURE PRESENTATION

The proposed architecture is described in Figure 6. It is composed of DF R and DF R I/O nodes linked together forming a ring topology. Note that N DF R and N DF R I/O represent the number of each DF R nodes type. To be able to interface directly K video sources, we use K + 1 rings: K unidirectionnal data flow rings (main video channels) and one ring in the inverse direction to synchronize the flow and to communicate some information or processing elements states. A. Architecture description 1) DataFlow Router: A DataFlow Router (Figure 7) is in charge of redirecting all the video data flows required for its Fig. 7. DataFlow Router internal structure: K = 3 connected PE. Its main actions are switching input data flows from video channels, PE and memory. The router contains buffers to synchronize input flows, mutiplexers, configuration registers and pixel counters to track the status of the processing element video output. It receives a command from the global controller containing its routing configuration and potential PE modification information. The DataFlow Router I/O (DF R I/O ) has a similar structure but is dedicated to video inputs-outputs and contains line buffers for each video channel.

2) Frame Buffer Server: The Frame Buffer Server stores all the frames needed for the video processing chain and results from PEs. It is a shared memory composed by several buffer memory slots dedicated to each PEs and those slots could be exchanged to share frame results. A frame buffer server controller manages all the frame requests from the global controller for the PEs.

3) Global controller: The global controller contains the video processing chain sequence and all the router adaptation required. Its role is to configure all the routers in the architecture and make all the necessary memory request to the frame server to give the right video flow when needed. At the beginning, the controller configures every router for each video channel (parallel or pipeline mode) and the DF R I/O is configured for video loopbacks or not. It can request video flows from the frame buffer server which capture permanently the input videos. The controller can decide to unload the pipeline from one processing element to do some parallel tasks for the following frames. It dynamically configure the DataFlow Router for a specified video channel. The video channel bypasses the router and the processing element can work transparently on a video data flow from memory. The DF R i reports to the controller that the P E i task is done when the processed pixel counters for the processing element is reached. Let T COM the time for the controller to sets the DFR configuration registers, the processing element task and the quantity of pixels to process (a complete frame or frame part) by setting the DFR pixel counters. T DF R is the DFR reconfiguration time and T P E the processing element one. The total adaptation time T ad for a DF R i is the sum of those times [START_REF] Hartenstein | A decade of reconfigurable computing: a visionary retrospective[END_REF].

T ad (i) = T COM (i) + T DF R (i) + T P E (i) (1)
4) Processing Elements: Each PE can accept up to three different video flows to be combinated and it outputs the processed video flow. In this paper, PEs can be considered as black boxes and all the local memory managements is done inside them according to the input video flows from different sources.

IV. FPGA IMPLEMENTATION

We propose to implement a video application illustrated by Figure 8 combining several pipelined image filters. It inputs two video flows that are combinated in T 3 . We will then compare the performance of our solution compared to a full dedicated one. To implement in a three-nodes solution for instance, we divide the application in two parts: a pipeline part and a parallel part. Those two computing parts will be temporally multiplexed in the ring and our structure has to be dynamically adapted to realize the different datapaths. Obviously, making a loop with less PEs will have an impact in the global latency unless we can increase the frequency. Our goal is to show that with only three PEs, it is possible to implement the application above by one video loopback and a pipeline-parallel adaptation with our proposed structure.

A. Implementation

This design has been implemented on an Altera Stratix III EP3SL150 FPGA (150K logic elements) [START_REF] Altera | Stratix iii -user handbook, www.altera.com[END_REF]. To test our proposition, we prototype the architecture with three DF R nodes (N DF R = 3), one DF R I/O (N DF R I/O = 1) and three main video channels (K = 3) as illustrated in Figure 6 and7.

B. Surface and Performance results

Let S DF R and S DF R I/O the surface occupied by DF R and DF R I/O respectively. S ctrl is the surface of the controllers and the communication links. The surface overhead cost S cost compared to a full dedicated implementation is given in Equation [START_REF] Clienti | A system on chip dedicated to pipeline neighborhood processing for mathematical morphology[END_REF].

S cost = NDF R i=1 S DF R (i) + NDF R i/o j=1 S DF R i/o (j) + S ctrl (2)
Results in FPGA surface are shown in Table I Figure 9 shows the pixel outputs of each PEs for a burst of pixels a 1080p image (burst size = 1920). We can notice the processing time and the latency comparison between a full dedicated solution and our adaptable solution. We can see in this figure that Lat DF R is negligible or even hidden by PEs latencies. Note that our solution needs an additionnal buffering (during T BufferP E3 ) in DF R I/O to realize the data loopback. We measured a latency of 1998 clock cycles (it is related to burst size value) with three DFR nodes. However, by keeping the same structure and using six DFR nodes, we have similar performances in latency compared to the full dedicated one as shown in Figure 9 because Lat DF R is negligible. Note that V ideo in 0 from Figure 8, is directly stored in the frame buffer and injected when needed on the right node for T 3 contrary to a full dedicated structure with a fixed size buffer for synchronization.

From this implementation, we can conclude three things: First, for a given application, we could use the right number of nodes to get performance closed to a full dedicated structure and when needed, by keeping the same structure, our system can work with less number of nodes by making several video loopbacks thanks to the ring structure and the pipelineparallel adaptation mechanism. Thus, the remaining nodes can work on others tasks in parallel. Secondly, our system can work on different parallel video flows and the PEs can work transparently without knowing the source either from a PE or the frame server which is important for design reuse. Finally, depending on the PE complexity in time and surface, we show that the cost in surface for implementing our proposition is relatively interesting for gaining such flexibility on datapath and performance.

V. CONCLUSION

In this paper, we have presented an adaptable ring-based architecture proposition that could handle both pipeline and parallel processing in the same application for an input video flow. A simplified version with four nodes has been designed on an FPGA for a video application with two video flow inputs showing the parallel-pipeline adaptation mechanisms. Future works will focus on the study of scalability of the structure. Consumption aspects will also be analyzed.

Fig. 1 .

 1 Fig. 1. Basic Video Pipe from Task 0 to Task 3

Fig. 2 .

 2 Fig. 2. Video Pipe with multiple input data flows

Fig. 3 .

 3 Fig. 3. Video Pipe with a Task 4 and 5

Fig. 4 .

 4 Fig. 4. Video Pipe with frame buffer requests

Fig. 5 .

 5 Fig. 5. Datapath modification requirement examples

Fig. 6 .

 6 Fig. 6. Global architecture: N DF R = 3, N DF R I/O = 1 and K = 3

Fig. 8 .

 8 Fig. 8. Video application example

Fig. 9 .

 9 Fig. 9. comparison between architecture solutions

  and S cost represents 762 Logic Elements (LEs) for our implementation. Similar PEs of 1452 LEs are used and a global surface comparison is given in Table II. Let Lat DF R the latency

  We also use PEs which have the same latency (Lat P E (i) = 18 clock cycles). The global controller is dedicated to the structure with specific input and output ports and finite state machines. All the DFRs have a direct connection on the global controller (T COM (i) = 3) and the DFR modifications are done by multiplexers controlled by registers (T DF R (i) = 2 and Lat DF R = 2). Consequently, T ad = 6 clock cycles for our implementation.

		TABLE I
	FPGA SURFACE RESULTS FOR ROUTING AND CONTROL ELEMENTS
	(ALTERA STRATIX III EP3SL150)
		S DF R	S DF R I/O	S ctrl
	Surface (LEs)	122	181	215
	between two neighbours PEs across the routers and Lat P E (i)
	the latency for P E i . For our implementation, the PEs are
	dedicated structure that works on the flow and changing the
	functionality of a PE is simply done by switching coefficients
	(T P E (i) = 1).


Harold Phelippeau

Mohamed Akil

Breno Dias Rodrigues

Hugues Talbot

Stefan Bara

Bilat- Eral Bayer

H Phelippeau Ab

Bayer Bilateral denoising on TriMedia 3270

Keywords: data load, fetches, custom operations, data cache

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

In photographic digital imaging, sensor noise continues to be a major problem for the output image quality. Embedded devices impose power processing constraints and thus algorithmic limitations. The bilateral filter is an interesting solution because of its simplicity and of its efficient noise removing performances. We focus on the Bayer bilateral filter proposed in [START_REF] Phelippeau | Bayer bilateral denoising filter[END_REF], a sensor adapted version. The Bayer bilateral filter avoids the perturbation of the noise distribution introduced by the demosaicing step and allows an efficient noise removing. It also has a lower computational complexity due to a lower amount of data to process [START_REF] Phelippeau | Bayer bilateral denoising filter[END_REF]. In this paper we present the execution simulation and the code optimization of the Bayer bilateral filter on the TM3270 media-processor from the NXP Semiconductors TriMedia family. Our goal is to perform real time-processing at VGA resolution using the maximum processor clock ratio (350MHz). We first present a common photographic imaging processing chain. We then recall the bilateral filter properties and its formulation for the Bayer pattern. In section 5 we present the TM3270 architecture and its programming environment. In section 6, we present the techniques used for the code optimization. It includes common methods (such as LUT, loop unrolling, convenient data type) and custom TriMedia operations. In section 6.6, we present the optimization performances, which we then discuss. In section 7, we propose a new Bayer bilateral filter formulation dedicated for the TM3270 architecture. In section 7.2, we present our new optimization performances results. Finally, in section 8, we present the general conclusion and the perspectives of this work.

PHOTOGRAPHIC IMAGING PROCESSING CHAIN

As it is presented in [START_REF] Phelippeau | Bayer bilateral denoising filter[END_REF], a digital image is the result of three main steps, the optical image formation trough the lenses system, the conversion of the light into electrical signal and some image processing operations, including: demosaicing, white balancing, noise removing, gamma curve adjustment, etc. A complete documentation concerning computational photography is presented in [START_REF] Raskar | Computational photography[END_REF]. Figure 1 shows a block diagram of a common imaging system. Such a system includes the optical lenses, the color filter array (CFA) and the image sensor (CCD, CMOS). It also includes the main control systems, as automatic gain control (AGC), analog to digital converter (ADC), auto focus and auto exposure circuitry. The digital signal path is completed with color and digital image processing, and is finally sent to the baseband for storage, or to the interface for visualization.

BILATERAL FILTER PROPERTIES

The bilateral filter is a non-linear filter well suited for denoising applications. It exhibits demonstrated effectiveness properties and its formulation simplicity contributes to its popularity. The bilateral filter replaces a pixel value in an image with a weighted mean of its neighbors considering both their geometric closeness and photometric similarities [START_REF] Tomasi | Bilateral filtering for gray and color images[END_REF][START_REF] Aurich | Non-linear gaussian filters performing edge preserving diffusion[END_REF][START_REF] Smith | Susan a new approach to low level image processing[END_REF][START_REF] Paris | A fast approximation of the bilateral filter using a signal processing approach[END_REF]. The most popular version is the Gaussian bilateral filter, it is expressed as follows:

( ) ( ) 1 ( ) exp exp ( ) 2 2 ( ) x y u x u y v x u y C x h β ρ -- - - = ∑ , (1)
where β represents the sliding window, y is a set of 2-D pixel positions in the sliding window, and x is the 2-D position of the centered pixel in the sliding window. ( ) u x is the intensity of the pixel at the x position in the original image, ( ) v x is the estimated pixel at the x position, ρ and h are respectively the standard deviation of the Gaussian distribution of the geometrical and the intensity weight. The filter behavior depends on the parameters setting. Figure 2 shows an example of the bilateral filtering behavior using a varying h .

BAYER BILATERAL FILTER FORMULATION

In this section we made a rapid presentation of the Bayer bilateral filter proposed in [START_REF] Phelippeau | Bayer bilateral denoising filter[END_REF]. We present first the Bayer pattern mainly used in single sensor devices. We present then the adaptation of the bilateral filter to the Bayer pattern.

Bayer color filter array

A major defect of single sensors is that they can get only one signal information per photosite. To introduce color information a color filter is superimposed on the sensor. There exist different color filter schemes, the most widely used being the Bayer arrangement [START_REF] Bayer | Color Imaging Array[END_REF], it is shown on figure 3.

Bilateral filter formulation for the Bayer pattern

In [START_REF] Phelippeau | Bayer bilateral denoising filter[END_REF], the authors propose to apply the bilateral filter before the demosaicing step. Indeed, filtering the color filter array allows for better noise removing performance by avoiding the perturbation of the noise distribution induced by the demosaicing step. The Bayer bilateral filter consists of processing each color separately. Green, red and blue filtered values are respectively calculated using green, red and blue pixels. The Bayer bilateral filter is expressed as follows:

( ) ( ) 1 ' ' ( ) ( ')exp exp ( ) ' 2 2 ( ) c c c c x y u x u y c c c c v x c c u y C x h β δ ρ - - - - = - ∑ (2)
The symbols are the same as in equation ( 1), δ is a discrete Dirac function, c is the color filter characteristic of the current pixel, ' c is the color characteristic of a pixel used for the weighted mean, and c β is the size of the square windows. Due to the Bayer pattern arrangement, the size of c β vary considering the processed color. This way we keep the same number of pixels for the weighted mean calculation irrespective of the different color arrangements. Figure 4 shows the kernel windows c β used in the case of a 25 weighted mean pixels. c β is a 7x7 square window for the green pixels convolution, which becomes a 9x9 square window for the red and blue pixels convolution.

Computational complexity

The Bayer bilateral filter is a mean weighted filter, it has a linear computational complexity O(n), where n represents the number of processed pixels. Calculations performed for each pixel are listed on Table 1.

Table 1. Bayer bilateral filter algorithm computational complexity.

TM3270 ARCHITECTURE OVERVIEW

TriMedia TM3270 Processor

The TM3270 is a VLIW media-processor from the NXP Semiconductors TriMedia family. Complete characteristics concerning the TriMedia architecture can be found in [START_REF] Rathnam | An architectural overview of the programmable multimedia processor, tm-1[END_REF][START_REF] Halfhill | Philips power up for video[END_REF][START_REF] Hennessy | Computers Architecture: A Quantitative Approach[END_REF][START_REF] Van De Waerdt | The TM3270 media-processor[END_REF].

TriMedia Programming Environment

The TriMedia processors are supported by a TriMedia Compilation and Simulation System (TCS) that speeds up creation of multimedia applications entirely in the C and C++ programming languages, with an ANSI-compliant C/C++ compilation system. The TCS provides a suite of system software tools to compile and debug multimedia applications, analyze and optimize performances, and simulate execution on a TriMedia processor by a cycle-close machine-level simulator.

CODE OPTIMIZATION ON TM3270 ARCHITECTURE

In this section, we seek to optimize the Bayer bilateral filter for the TM3270 architecture using the TriMedia Compilation System. Various optimization methods supported by the TCS as well as techniques for exploiting the finegrain parallelism can be used. We use common code optimization techniques as Look Up Table, loop unrolling, and TriMedia custom operations to take advantage of the TriMedia architecture.

Naïve code version

The naïve code version is the version of the code without any optimizations. This version uses if and else conditions, complex mathematic operations (natural exponential) and none loop optimizations. The pseudo-code of the Bayer bilateral filter is shown on figure 5.

Look up table

The Bayer bilateral filter in the Gaussian version (equation ( 1)) uses natural exponential computations depending on pixels intensities and locations. The computational complexity can be decreased by replacing these functions by simple memory fetches using look up tables (LUT).

Location weights

Using a 9x9 square window for the red and blue pixels convolution and a 7x7 square window for the green pixels convolution, we deduce that there are 80 possible pixel positions, which lead to use a LUT of 80 float values for the position weights. To reduce the size of the LUT we can take advantage of the Bayer pattern symmetry properties. Indeed, the position of a pixel comparing to the centered pixel is calculated using the 2-d Euclidian distance, it induces no direction dependencies. If we consider an orthogonal axis (I,J) centered at the current pixel position of a 7x7 square window (figure 6), a green pixel can take the coordinates: 1, 2 or 3, on the I and J axis. We can deduce the number of 2-D coordinates permutations without repetition using the combinatory laws: 2

(3) 6 = . Excluding the coordinate (0,0), we find 5 unique 2-d Euclidian distances. They are presented using white dots on figure 6(a). Considering now the 9x9 window used for the red and blue pixels convolution, using the same reasoning we find also 5 unique 2-d Euclidian distances, they are represented using white dots on figure 6(b). It can be reduced to 3 possible cases, because 2 of them have already been counted in the green case. They are represented using black dots on figure 6(b). Finally, we concluded that all the location cases comparing to the centered pixel of the convolution kernel, can be stored using a LUT of only 8 float values.

LUT fusion

Using two LUT induces two fetches and one multiplication for the calculation of one pixel weight. We can decrease the number of operations by making a fusion between the location LUT and the intensity LUT. Using this idea, we create 8 LUT of 255 float values containing pre-multiplied location and intensity weight values. As a consequence, we decrease our computation const from two fetches and one multiplication to a single fetch. Interpreting these results, we see that by using all these optimizations it is still not possible to achieve real-time processing at VGA resolution. We seek to find why the algorithm formulation is not adapted for the TM3270 processor:

• The first problem appears by looking at the kernel convolution size. Indeed loading lines of 7 and 9 pixels of 1byte values is not adapted for an efficient use of the super_ld32 function (that allows loading 8-bytes per call).

Particularly, loading a line of 9 pixels of 1-byte values requires two super_ld32 function call. We estimate that using a 8x8 square convolution window would increase the yield of the super_ld32 function.

• The second problem is the utilization efficiency of the loaded values. Processing the pixels one by one reduce the utilization efficiency of the loaded values. Indeed because of the Bayer matrix arrangement 50% of the values loaded to process one pixel are not used. We conclude that we must take advantage of the Bayer arrangement and process red, green, and blue values simultaneously to increase the utilization efficiency of the loaded values.

In the following section, we seek to adapt the Bayer bilateral filter formulation considering the discussion above to increase the optimization performances.

NEW BAYER BILATERAL FILTER FORMULATION DEDICATED FOR TRIMEDIA 3270 7.1 New bayer bilateral filter formulation

Considering the conclusion of section 6, we seek to propose a new Bayer bilateral filtering formulation:

• To optimize the load function calls, we impose an 8x8 square kernel windows. This configuration enables the loading of a complete kernel line per load instruction, so all the pixels values of kernel can be charged with 8 super_ld32 instructions. Two 32-bits destination variables are used to store the pixels values loaded, in an unsigned char codification (4 pixel values are stored in each destination variable), allowing the four-operand instructions use. The pixel values loaded are ready to use by TriMedia special functions, like dspuquadabssub().

• To optimize the utilization efficiency of the loaded values we impose to process red, green and blue interpolation simultaneously. To do so, we can take advantage of the Bayer pattern arrangement. As we can see on figure 9(a), the Bayer pattern is the repetition of a 2x2 square element E that includes red, green, and blue pixel. We propose to process this kernel at each iteration. This way we process all the loaded pixels values from the window kernel by computing simultaneously the values of the element E. In addition, custom TriMedia instructions are used more effectively by avoiding data ranking.

The new convolution kernel can be represented as in figure 9(b), it as a 8x8 square size and a 2x2 centered element E.

G R G R G R G R G R B G B G B G B G B G B G B G B G B G R G R G R G R G R G B G B G B G B G R G R G R G R G R G B G B G B G B G R G R G R G R G R G B G B G B G B G R G R G R G R G R G B G B G B G B G E G R G R G R G R G R B G B G B G B G B G B G B G B G B G R G R G R G R G R G B G B G B G B G R G R G R G R G R G B G B G B G B G R G R G R G R G R G B G B G B G B G R G R G R G R G R G B G B G B G B G E (a) (b)

Experimental results

In this section we present the performances optimization of the new Bayer bilateral filter formulation. As in section 6, we use the TCS environment on the TM3270. As optimization measure we use the number of clock cycles, the % of improvement compared to the naïve code version, and the number of frames per second giving a VGA (640x480) image at the maximum clock ratio (350 MHz).

We used the new formulation approach with a convolution kernel line size up to 8 bytes (8x8 window), as discussed before. Loop unrolling was implicitly performed (we treat 4 pixels in a 2x2 central square at once) and explicitly by processing two central square element at one loop iteration. We use the LUTs system presented in section 6.

CONCLUSION

In this paper, we have simulated and optimized the Bayer bilateral filter execution on the TM3270 using the TCS tool.

Our goal was to reach real time processing using a VGA resolution image at the maximum processor clock ratio (350MHz). We have used common optimization method as LUT utilization, loop unrolling, convenient data type representation; we also used the custom TriMedia operations to take advantage of the TriMedia architecture. Despite using all this optimization tools and method it was still impossible to achieve real time video processing using the original formulation of the Bayer bilateral filter. The maximum speed was found to be 17.78 frames/sec. By analyzing the code and the custom TriMedia function call, we have found that the Bayer bilateral filter formulation does not allow an efficient utilization of the data load function (super_ld32). We have also found that is prevents an efficient utilization of the loaded data (we notice that 50% of the loaded data are not usef). To resolve these problems we have proposed a new Bayer bilateral filter formulation dedicated for the TM3270. It uses a 8x8 (or a 6x6, depending on the version) square convolution kernel , that results in an efficient use of the load function, and a 2x2 centered square element that results in a full utilization of the loaded data. Using this formulation it is finally possible to reach real-time processing with speeds of 35.27 frames/second for the 8x8 convolution kernel version, and 46.9 frames/second for the 6x6 convolution kernel. It improves the naïve version algorithm by 99.46% and 99.6% respectively. The real-time processing margin obtained with this two versions makes it possible to perform consecutive real-time processing.

Fig. 1 .

 1 Fig. 1. General block diagram of a digital camera.

Fig. 2 .

 2 Fig. 2. Gaussian bilateral filtering using a 5x5 square window as β and a variable h , from left to right

Fig. 3 .

 3 Fig. 3. Bayer color filter arrangement

Fig. 4 .

 4 Fig. 4. From (a) to (c), are presented the different convolution kernels used for the green, red and blue bilateral filtering

Fig. 5 .

 5 Fig. 5. Bayer bilateral filter naïve pseudo-code

Fig. 6 .

 6 Fig. 6. (a) shows with white dots the elementary positions used for the green interpolation, (b) shows with white and black dots the elementary positions used for the red and blue interpolation. Black dots mean common position between green and red or blue elementary position.

Fig. 9 .

 9 Fig. 9. Scheme of the new Bayer bilateral filter formulation

Fig. 10 .

 10 Fig. 10. Scheme of the new Bayer bilateral filter formulation

Table 2

 2

, gives an overview of the main TM3270 architectural features. It operates at up to 350MHz and can execute at up 5 instructions per clock cycle, tuned to match the performance demands of standard definition video processing. It supports SIMD multimedia and IEEE-754 floatingpoint operations. It uses 128 32-bits registers (unified register-file), a 128kB data cache (4-way set associative) and a 64kB instructions cache (8-way set associative), with 35 execution units and an extensive set of multimedia instructions for video and image processing.

  The first row of the table presents the results of the naïve version. We obtain 6026.1 cycles/pixel, which results in 0.19 frames/second. The second row presents the results obtained by using LUT, we can see the effectiveness of this code transformation. It allows having 92.6% of improvements comparing to the naïve version, giving 445.7 cycles/pixel and 2.56 frames/second. The third row presents the results of the loop unrolling optimization. It shows an important optimization gain with an improvement of 50.82% compared to the LUT version. It performs at 219.2 cycles/pixel and allows the CPU to process 5.2 frames/second. The passage to integer type values reduces the number of cycles/pixel to 155.5 or 7.33 frames/second, an improvement of 97.42% compared to the naïve version and 65.12% compared to the LUT version. Finally the TriMedia custom operations results in an efficient optimization, improving by 98.94% the naïve version and of 85.62% the LUT version. It makes it possible to process 17.78 frames/second and it reduces the number of cycles per pixel to 64.1.

Table 5 .

 5 Bayer bilateral filter algorithm optimization results on TriMedia TM3270 processor.

  [START_REF] Hennessy | Computers Architecture: A Quantitative Approach[END_REF] and we work with unsigned char types values. We also experiment another version of the algorithm using a convolution kernel with 6-bytes of line size (6x6 square window). The results of these experiments are shown on table 6 and on figure11.We can see that using our new algorithm formulation allows to improve the algorithm execution efficiency by 49.61% compared to the full optimized original version (corresponding to the TriMedia custom functions row in the table). It allows us to process 35.27 frames/second, thus achieving real time processing. Reducing the convolution kernel size to a 6x6 square improves the algorithm execution efficiency further to 61.93% compared to the full optimized original version. It allows us to process 46.9 frames/second, giving thus an important margin for real time processing.

Table 6 .

 6 Bayer bilateral filter algorithm optimization results on TM3270 processor: new formulation Fig. 11. Histogram of the Bayer bilateral filter algorithm optimization results on the TM3270 processor

		Frame/second (640x480) @350MHz
	50						46,9
	45
	40
						35,27
	35
	30
	25
	20				17,78
	15
	10		5,2	7,33
	5	2,56
	0,19
	0
	Naïve	Look Up	Loop	Integer	Trimedia	New	New
	version	Table	unrolling		custom	formulation-8x8	formulation-6x6
					functions	windows	windows

Computational complexity

Calculations performed for each pixel are listed on Table 3. Compared to the naïve code computational complexity (table 1), the LUT utilization allows an important number of simplifications. We note that the 50 exponential calculations are replaced by 25 fetches; we also note a gain of 75 multiplications, 50 additions, and 50 divisions that are directly due to the weights pre-calculation. Execution results of the code by using this method are presented on table 5.

Loop unrolling

Loop unrolling is done by replicated several iterations in a single sequence. In digital image processing it is done by computing several neighboring pixels in the same iteration. This technique increases significantly code execution by eliminating the iteration's jump instructions. It provides better processor pipeline utilization and operations scheduling.

By looking at equation ( 2), we remark that the Bayer bilateral filter computes different filtering treatments according to the pixel color. It induces control flow to select the best treatment (if conditions). We can apply loop unrolling to a Bayer matrix image by processing two or more consecutive pixels in an image line at each iteration. This can partially eliminate the conditional choices generated by the Bayer pattern because two different colors of a line are processed in the same loop iteration. Figure 7 illustrates a loop unrolling in the Bayer bilateral filter algorithm. We note that the control flow appears at each line in the unrolled version, while it appears at each pixel in the original naïve version. Experimental results of using this technique are given in section 6.6.

Convenient data type representation and custom TriMedia operations

Data type representation has an important impact in code optimization. It has influences on the number of cycles per operation, cache memory usage, and storage. To increase the code execution performances it is important to use data types as small as possible. For example, working with char values also reduces the memory required for the Look Up Tables from 8192 bytes of memory (8 LUT of 255 floating point values) to 2048 bytes (8LUT of 255 char values). Optimization results obtained by working with unsigned char values are listed on table 5. Furthermore, using the TM3270 multimedia processor, it is preferable to work with integer values to take advantage as much as possible of the TM3270 custom operations that use 4-bytes registers. Indeed, it is possible to represent 4, 2 or 1 pixel(s) in each 4-bytes registers depending of their data representation (8-, 16-or 32-bits). These operations allow the processor to compute up to four operations concurrently and reduce the total number of cycles of the code section. Among this set of custom operations the TriMedia library proposes absolute value calculation, maximum or minimum function that can be executed faster than usual C/C++ operations. Figure 8 shows an example of such an operation using the dspuquadabssub function used to calculate the rank in one of the 8 intensity weight LUT. This function takes in input two 4-bytes registers that contain four 1-bytes words (pixels). It returns the results of four absolute values differences in a 4-bytes register. We use the function super_ld32 to retrieve consecutive 8-bytes data from memory in two 4-bytes destination registers. This doubles the memory-processor bandwidth and allows retrieving an image window's convolution line using one instruction. Table 4 illustrates the number of processor cycles necessaries to perform four pixel treatment. The experimental results of custom TM3270 functions use are listed on table 5. 4. Processor cycles gain using custom TriMedia operations

Experimental Results

In this section we present the performances of the optimization steps presented above. We use the TCS environment on the TM3270. As optimization measure we use the number of clock cycles, the % of improvement compared to the naïve code version, and the number of frames per second giving a VGA (640x480) image at the maximum clock ratio (350 MHz). VGA is the common embedded devices video resolution and real-time processing is possible from 25 frames/second. Our goal is to obtain a real-time processing at VGA resolution. Optimization results are listed on Table 4.


Harold Phelippeau
email: phelipph@esiee.fr

Hugues Talbot
email: talboth@esiee.fr

Stefan Bara
email: stefan.bara@nxp.com

Mohamed Akil
email: akilm@esiee.fr

Shot-noise adaptive bilateral filter

Keywords: Digital imaging, image sensors, Poisson noise model, noise filtering, noise reduction, noise estimation

published or not. The documents may come   L'archive ouverte pluridisciplinaire

INTRODUCTION

The bilateral filter is a non linear filter used among others for denoising applications. In this paper we propose a method to estimate automatically this standard deviation using the variance and the mean of smooth image region.

BILATERAL FILTER CHARACTERISTICS

The bilateral filter consists to replace a pixel by a weighted mean of its neighbors considering both their geometric closeness and their photometric similarities. In this paper, we focus on the Gaussian bilateral filter because all the practical applications use this version. The Gaussian bilateral filter is defined as follows:

v(x) = 1 C(x) β e -|x-y| 2 ρ 2 e -|u(y)-u(x)| 2 h 2
.u(y).dy [START_REF] Tomasi | Bilateral filtering for gray and color images[END_REF] Where β represents the sliding window, y represents a set of 2-D pixel positions in the sliding window, x represents the centered pixel in the sliding window, u(x) is the intensity of the pixel at the x position in the original image, v(x) is the estimated pixel at the x position, ρ and h are respectively the standard deviation of the Gaussian distribution used for the geometrical weight calculation and the intensity weight calculation.

The parameters ρ and h define the bilateral filtering behavior. The geometrical parameter is chosen considering the size of the convolution kernel. The intensity parameter h has to be chosen considering the level of filtering needed for the application. Indeed, the more the standard deviation is high the more the filter is low-pass in the intensity space. An illustration of the h influence is shown on Fig. 1. In the case of noise removing the parameter h has to be chosen considering the noise level. The best h could be selected among a set of h by chosing the one that has the lowest mean squarre error measure between the denoised image and the original version. Repeating this task for each image filtering is not acceptable in term of time consumption and computational complexity. In order to avoid this task, a solution consists to calibrate the best h for each noise level. Using this calibration, estimating the best h amouts to estimate the noise level in the image. In the following section we present a new method to estimate shot noise level in a poisson processus acquired image.

BEST H PARAMETER ESTIMATION BASED ON IMAGE CONTENT

Sensor noise assumption

Many sources are cause of noise generation in CCD and CMOS sensors. These can be categorized in four main factors [START_REF] Costantini | Virtual sensor design[END_REF]: [START_REF] Tomasi | Bilateral filtering for gray and color images[END_REF] The photon shot noise -associated with a random Poisson process governing the number of incident photons reaching a photosite; (2) the Photon Response Non Uniformity -caused by small sensitivity differences between photosites;

(3) the dark current noise -produced by minority carriers thermically generated in the sensors wells, also associated with a random Poisson process; and (4) the readout noise -resulting from thermal noise cause by MOSFET amplifiers.

(a) h=0

(b) h=10 (c) h=20 (d) h=40
Fig. 1. Bilateral filtering using a 5x5 square window as β and a variable h

A sensor noise model complete with respect to these factors is proposed in [START_REF] Costantini | Virtual sensor design[END_REF]. According to this model, photon shot noise has, in standard illumination conditions, the most important influence on the output image.

Automatic gain control

In a digital imaging system, an automatic gain control (AGC) allows to capture images in varying light conditions by adjusting the average intensity of the output signal (see Fig. 2). By having access to the gain factor and knowing the sensor resolution, it is possible to deduce the number of photons that have reached the sensor and estimating then the noise level. In the following section we propose a method to estimate the applied gain by using images statistics and the proporties of a poisson distribution process.

Photon density from image statitics

In the presence of gain manipulation only, we expect the noise variance to be proportional to average pixel values in nearconstant regions, due to its Poisson characteristics. Let consider I the image before AGC, and I G the adjusted image after the AGC. Assuming G to be the gain factor, we can write:

I G = G × I; (2
From image I G we can write:

V ar(T -(I)) = E(T -(I G ) 2 ) -E(T -(I G )) 2 (4)
Using 2 we can write:

T -(I G ) = G × T -(I) (5)
E(T -(I)) = G × E(T -(I G )) (6)
Replacing the expressions 5 and 6 in 4 we obtain an estimation of the gain:

G = E(T -(I)) V ar(T -(I)) (7)
The image mean E(I) and variance V ar(I) can be simultaneously calculated at each point in I by using an efficient sliding-window algorithm [START_REF] Sun | Fast algorithm for local statistics calculation for n-dimensional images[END_REF]. We seek to establish a correlation between all E(I) and V ar(I) at each point where T - is non-zero. Using all the millions of pixels generally available in modern camera sensors is not necessary, so a random sampling of perhaps approximately 10 4 pixels was found sufficient. Since the mean and the variance are expected to be proportional, a robust linear, intercept-free, least-square correlation should be performed, the slope coefficient can then be interpreted as the gain factor transforming photon numbers to recorded pixel values.

Regions of near-constant values

In order to estimate the gain, we need to select areas T -(I) in the noisy image I that correspond to approximately constant values. These can be readily estimated from the image gradient

T -(I) = ε w (∇(G σ ⋆ I) > ϑ), ( 8
)
where I is the image, G is a Gaussian convolution of variance σ, ε w the morphological erosion on a window of size w × w. ϑ is an intensity parameter.

Best h from photon density

The filtering properties of the Gaussian bilateral filter varies with the set of its two parameters. Using a fixed sliding window, the filtering properties vary with the h parameter. Intuitively we can assume that a low h should be required in low noise level, wheras a high h is necessary in high-noise conditions. A calibration of the best h for each noise level is necessary to establish a correspondance between the gain factor estimation and the adapted h. In the following we use a photon shot noise model to do this calibration.

EXPERIMENTS AND RESULTS

In this section seek to validate the previous assumptions.

Photon shot noise simulation

To do so we simulate sensor aquisition via a regionalized cumulative spatial Poisson point process [START_REF] Shinozuka | Digital simulation of random processes and its applications[END_REF] using scene images as probability distribution functions. This simulates individual photons being recorded at each pixel location, as in a Monte-Carlo process. As scene images we used the Kodak PhotoCD database, extracted at the 512x768 resolution. These virtually noise-free images were scanned at the 3 full color samples per pixel from film original following the Kodak professional PhotoCD procedure. Fig. 3 illustrate this simulation.

(

Photon number estimation

We ran trials with a photon density per pixel between 20 and 320, from an extremely noisy to a clean image. We estimated the number of photons Φ with the following formula:

Φ(I) = nbpix(I)E(I) G = (x,y)∈I I(x, y) G , (9)
where nbpix(I) is the total number of pixels in I. The excellent correlation between simulated and estimated photons is shown on the log-log plot of Fig. 5

CONCLUSION

In this paper, we have proposed a methodology for digital image filtering in the presence of photon shot noise and a color filter array. For this, we have studied a subset of known filters that can be efficiently implemented in low-power CPUs and DSPs found in consumer digital cameras. We have shown that modeling noise in this context by a Gaussian distribution over an RGB image is not accurate, especially at high noise levels. We have also shown that when the noise standard deviation is sufficiently low, then it is better not to filter at all, even with the best filters, using fixed parameters. It is therefore important to adapt the level of filtering to the level of noise. Among the tested filters, the bilateral filter was shown to be the one with best performance. In addition, we showed how this filter can be optimised if the noise level is known. We have proposed a simple and efficient noise level detection algorithm. This algorithm was shown to exhibit good correlation between estimated and real noise level. By applying this

Fig. 2 .

 2 Fig. 2. Digital conversion of a photon flow by a camera sensor, the AGC adjusts the average image output intensity

)

  Considering a Poisson distribution process, we can assume the following relation where T -(I) represents regions of near constant values in the image I. E(T -(I)) = V ar(T -(I))

Fig. 3 .

 3 Fig. 3. Shot noise simulation using a sample representative image as a source and a per-pixel Poisson process. From (a) to (c) the mean number of photons per pixel over the whole image are 20, 80, 160. From (d) to (f) the corresponding images with a gain control.

  (R 2 = 0.996).

Fig. 4 .Fig. 5 .

 45 Fig. 4. Correlation between simulated and estimated photon numbers.

Fig. 6 .

 6 Fig. 6. Shot noise simulation using a sample representative image as a source and a per-pixel Poisson process. From (a) to (c) the mean number of photons per pixel over the whole image are 20, 80, 160. From (d) to (f) the corresponding images with a gain control.

noise level detection method to set the bilateral filter h parameter, we obtain a novel adaptive bilateral filter, which exhibits best-of-class performance at all noise levels.


Petr Dokladal
email: petr.dokladal@ensmp.fr

Eva Dokladalova
email: e.dokladalova@esiee.fr

GREY-SCALE 1-D DILATIONS WITH SPATIALLY-VARIANT STRUCTURING ELEMENTS IN LINEAR TIME

Spatially variant morphological operators can significantly improve filtering capabilities or object detection score of various applications. Whereas an effort has been made to define the theoretical background, the efficient implementation of adaptable algorithms remained far less considered.

In this paper, we present an efficient, one-scan, linear algorithm for 1-D grey-scale dilations/erosions with spatially variant structuring elements. The proposed algorithm processes data in stream, can work in place and produces results with minimal latency. The computing time is independent of the structuring element size.

INTRODUCTION

The concept of adaptive filters is well known in the image processing community since years [16,[START_REF] Saint-Marc | Adaptive smoothing: A general tool for early vision[END_REF][START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF]. The general principle consists in applying a filtering process which depends on internal, local properties of image [START_REF] Lerallut | Image filtering using morphological amoebas[END_REF], [START_REF] Braga-Neto | Filters by adaptive neighborhood structuring functions[END_REF] or on external, application dependent context, as e.g. in [16,[START_REF] Hedberg | Binary Morphology with Locally Adaptive Structuring Elements: Algorithm and Architecture[END_REF], or [START_REF] Roerdink | Mathematical morphology for structures without translation symmetry[END_REF] where the distance-to-camera or wide-angle deformations require adapting the structuring element.

In the domain of mathematical morphology, a considerable effort has been made to define the theoretical background. In [START_REF] Roerdink | Group morphology[END_REF], Roerding has assembled his work on spatially variant morphology in the general framework of Group Morphology. Charif-Chefchaouni and Schonfeld [START_REF] Charif-Chefchaouni | Spatiallyvariant mathematical morphology[END_REF], and also Bouaynaya et al [START_REF] Bouaynaya | Theoretical foundations of spatially-variant mathematical morphology -part i: Binary images[END_REF] and [START_REF] Bouaynaya | Theoretical foundations of spatially-variant mathematical morphology -part ii: Gray-level images[END_REF] examine a theory of spatially variant binary and functional mathematical morphology where the structuring elements can vary both in size and shape.

There is an evident lack of efficient implementations on gray-scale images needed for practical applications, e.g. Lerallut [START_REF] Lerallut | Image filtering using morphological amoebas[END_REF] states running time in 3-D images up to hours.

Efficient implementations

Whereas the efficient algorithms have beeen proposed for translation invariant structuring elements (SE) for both binary and grey-scale images, the efficient algorithms addressing the spatially variant SE are very limited.

Van Herk [START_REF] Van Herk | A fast algorithm for local minimum and maximum filters on rectangular and octagonal kernels[END_REF] proposes a 1-D fast algorithm, independent of the SE size, requiring three scans. The extension to translation invariant SE is not straightforward as the algorithm relies on a fixed division of the processed line by the size of SE.

Lemonnier and Klein [START_REF] Lemonnier | Fast dilation by large 1D structuring elements[END_REF] (also [START_REF] Lemonnier | Architecture électronique dédiée aux algorithmes rapides de segmentation basés sur la morphologie mathématique[END_REF]) propose a fast dilation with large 1D for symetric SE that completes in two scans of each line. The principle is the propagation of maxima according to the SE size. The first scan can result an error to be corrected by the second scan. The extension to adaptable SEs is either difficult or impossible.

Van Droogenbroeck and Talbot [START_REF] Van Droogenbroeck | Fast computation of morphological operations with arbitrary structuring elements[END_REF] propose an algorithm based on a histogram, updated as the SE slides over the image. The extension to spatially variant SE is direct. Computing the histogram however requires additional memory resources and brings in an additional computing complexity. Notice also that the histogram is not usable on floating-point data.

In the domain of adaptable SE, several efficient algorithms have been published for binary domain, see e.g. [START_REF] Cuisenaire | Locally adaptable mathematical morphology using distance transformations[END_REF]. It proposes a spatially variant binary morphology, using as kernels balls in various norms. Dilations are obtained by thresholding the distance function to objects. Local adaptability is achieved by thresholding at different levels. The computational complexity is the one of computing the distance. However, the only possible extension to the grey-scale domain is by decomposing the image into a set of binary images, bringing in an overwhelming increase of computational burden.

Hedberg et al. [START_REF] Hedberg | Binary Morphology with Locally Adaptive Structuring Elements: Algorithm and Architecture[END_REF] propose an efficient algorithm for spatially variant binary morphology. The algorithm accepts non centered resizable rectangles, uses little memory and its latency strictly equals the one imposed by data dependency.

In the gray-scale domain, Cheng [START_REF] Cheng | Adaptive morphological operators, fast algorithms and their applications[END_REF] proposes a fast algorithm for SE with adaptable shape and size. The algorithm proceeds by SE decomposition in smaller 2-D SE. Hence, the computing complexity heavily depends on the SE size and shape.

The algorithm presented in this paper is one-scan, linear algorithm for 1-D dilation with spatially variant, non centered SEs. The computing time is independent of the SE size. It's can be considered as an extension to gray-scale domain of the Hedberg algorithm.

The paper is organized as follows. Section 2 discusses the constraints applying to adaptable SEs. Then the algorithm is described in the Section 3. Finally Section 4 gives some application results. The text concludes with discussion of the properties and extension to higher dimensions.

ALGORITHMIC ISSUES

Consider a 1-D binary object X : R → {0, 1}, see Fig. 1a. Consider a collection of flat, one-sided structuring elements B defined for all x ∈ R called origin by the length β : R → R measured from x. Hence, every B is an interval in R delimited by (x, x + β (x)) ⊂ R. Notice that β can be positive or negative for SEs stretching right-or leftwards from the origin. 1b shows an example of dilation of X by a function B 1 , illustrated at i and j (the origin of B 1 (i) is given by the black dot above i). Notice that B 1 ( j) may be different from B 1 (i). The dilated object δ B 1 X(x) = 0, for x = i, j.

(a) (b) (c) 0 1 0 1 0 1 X B 1 (i) x x x B 2 (i) B 2 ( j) B 1 ( j) δ B X δ B X i j i j i j Figure 1: (a) A 1-D binary object X. (b) Dilation δ B X with a flexible SE B(x). (c) Example of a discontinuous result δ B X with B varying "too fast". The morphological dilation of X by B is defined δ B X = ∪ b∈B X + b. Fig.
Fig. 1c shows the same X dilated by a different function B 2 . The dilated object δ B 2 X(i) = 0 but δ B 2 X( j) = 1 since B 2 ( j) doesn't extend inside X as does B 1 ( j) in Fig. 1b. Dilating by B 2 has created an echo of X at i instead of dilating it.

SE size map B preventing this undesired behaviour needs to verify following restrictions :

∂ β ∂ x ≤ 1 (1)
For the rest of the paper, consider Z instead of R, and bidirectional structuring elements as intervals

B : Z → (a, b) ⊂ Z, a ≤ x ≤ b, encoded by distance from the origin x towards left and right L(x) = x -a and R(x) = b -x. Hence, ∀x ∈ Z we have B(x) = (L, R), with L, R ∈ Z + .
The condition Eq.1 becomes

∆L ∆x ≤ 1 and ∆R ∆x ≤ 1, (2)
These conditions allow a fast implementation of dilations, and appear several times below throughout the description of the algorithm.

ALGORITHM DESCRIPTION

Dilation and thresholding commute. This well known property is utilized for brute-force implementations of dilations of functions. A function is decomposed into a collection of binary cuts obtained by thresholding it at every level. Every cut is dilated and the function recomposed by superposition. The fact that not all binary dilations in this collection are necessary gave birth to the following algorithm. The algorithm consists of two main stages.

1) Anti-causal to causal conversion

Suppose input data f arriving in a stream, see Fig. 2b. Any operation utilizing a bi-directional structuring element B = (L, R) with R > 0 is an anticausal operation, Fig. 2a. Computing the result δ f B (x) depending of future samples of the signal will require storing f in memory. Utilising the property that dilation and translation commute

δ B-t f (x) ≡ δ B f (x + t)
any anticausal B = (L, R) has a causal equivalent B < = (L + R, 0), see Fig. 2c, yielding, after translation, an identical result

2) Propagating values

The algorithm reads in a stream the input data and the local structuring element size and writes in a stream the output, every sample is read/written once and only once. No jumps ahead or backwards are used. Notice that the reading and writing positions are not the same. Their difference is the algorithm latency.

See the result by definition Fig. 2f, dilating f by some B assigns δ B(x) f (x) the maximum found in the interval covered by B(x). Using a causal B < can be seen as propagating some values f (x) a number of positions towards right. The propagation distance, see Fig. 2d, is determined by finding the last B < (x + k), k > 0, that still includes x. Formally, the signal value f (x) is dilated at most until x + d(x), i.e. d(x) steps ahead:

d(x) = max k≥0 {k | x + k -L(x + k) ≤ x} (3)
Utilizing Eq. 2, it can be shown that Eq. 3 can be implemented efficiently by

d(x) = max k≥d(x-1) {k | x + k -L(x + k) ≤ x} (4)
a scheme which allows sequential reading on both signal and structuring element data. Algorithm Not all values can propagate as soon as they are read. This happens whenever a higher value is already being propagated. Smaller, more recent values need to be stored for possible later propagation. However, not all stored values are necessarily propagated. They can get masked by another incoming higher value.

Values stored for later propagation are memorized in a queue structure. This queue behaves like a FIFO, supporting inserting new elements with push and retrieving the oldest Refer to the pseudo code in Appendix. The propagation algorithm consits of following stages :

1) Read next signal value Fnext 2) Compute distance to propagate using scheme Eq. 4.

3) Compute the right and left edge (RE), (LE) of the structuring element at given writing position (wp). Note that before outputting a value (at position wp), all input data covered by its structuring elements need to be read: (cycle while 1), broken by one of breaks as soon as new value can be written.

4) Dequeue all stored smaller values that will never propagate, being masked by Fval.

5) Detect a downstep and store the corresponding triplet for the future propagation.

6) Compare the value currently being propagated with Fval. If Fval is higher, use Fval for propagation. Otherwise, dequeue all smaller values than Fval, then store Fval, its position x and distance d(x).

7) Detect the end of current propagation TT(2)+TT(3)<wp, meaning compare whether its initial position+distance is shorter that the current output position wp. If yes, find another value to propagate: erase too old values from the queue (no longer covered by current structuring element), and pop -if exists -next stored value higher that Fval, otherwise, use Fval for propagation. Memorize in TT its value, position and distance to propagate. Whenever the reading position reaches RE(wp), write at wp of the output dF the maximum of Fval and the currently propagated value T.

Algorithm Features

Input/Output data are read/written sequentially in a stream, every sample is read/written once and only once.

The algorithm latency (difference between reading and writing position) is at any moment strictly equal to the data dependency.

Data don't need to be allocated in the memory, i.e. there's no memory limitation to size of images.

Algorithm is able to run in place, i.e. input data can be overwritten by output data.

Computation Complexity and Memory requirements

The algorithm complexity is O(N), N being the amount of data (e.g. the number of pixels). Every sample is read and processed (or dequeued) only once. The complexity is constant with respect to the size of structuring element.

The data can be read and output in a stream, and don't even need to be stored in memory. Therefore, the only memory occupation is given by the maximal depth of the FIFO. The worst case memory occupation of the algorithm is obtained on a non-increasing interval, larger than the structuring element size used at this position.

EXPERIMENTAL RESULTS

Consider detail-preserving noise filtering or contour-aware image simplification. We illustrate the performance on the Manet's Le fifre painting, given by Fig. 3a. The contours to preserve (c) can be detected with any usual edge detector (we omit details here). Obviously, the contour preserving filtering must use spatially variant (SV) SE that do not cross these contours. Moreover, techniques based on mathematical morphology using SV SE yield in less iterations better results than with translation invariant (TI) SE.

After decomposition into 1-D, the present algorithm can be used as follows: (e) shows horizontal then vertical dilation by linear SV SE of size locally equal the distance (d) to contours. In this way it approximates1 2-D dilation with SV squares. The panel (f) gives approximation of closing, and (g) approximation of opening•closing (one stage ASF filter). Compare these results with TI ASF (h), obtained with

Extension to higher dimensions

Although 2-D dilation by a TI rectangles can be implemented by decomposition into two 1-D dilations (horizontal and vertical), separability of SV SE is not straightforward. For the present algorithm, restrictions on the SE map have to be introduced during its generation. Exact formulation of these restrictions is under investigation and will be published separately.

Benchmarks

Experiments confirm linear execution time with respect to size of data and constant time with respect to mean size of SE, see Fig. 4, obtained on a 2GHz Intel Core 2 CPU, with 800MHz DDR2 RAM.

CONCLUSION

This paper presents fast, 1-D, spatially variant, grey-scale dilation, obtained in one scan of data. Its complexity is linear with respect to size of data, and constant with size of structuring elements. Its latency is the strict minimum imposed by the operator dependency on input data. The memory requirements are only bounded by the maximal size of structuring element. Input/output data being read/written sequentially, the implementation is independent of size of data. This is extremely interesting for filtering large data with cascade filters, e.g. ASF, under severe time or memory constraints. These interesting properties make this algorithm useful not only for spatially variant morphological filtering, but also with translation-invariant structuring elements.

In general, this algorithm represents a next step in the effort towards efficient implementation of grey-scale morphology with spatially variant arbitrary-shaped structuring elements.

Figure 2 :

 2 Figure 2: 1-D dilation algorithm: (a) input two-sided structuring element (circle denotes origin), (b) input signal, (c) causal structuring element, (d) distance to propagate the input signal samples, (e) algorithm result, (f) result by definition.

Figure 3 :

 3 Figure 3: Image simplification application

Figure 4 :

 4 Figure 4: Execution time (a) vs. size of data, (b) vs. mean SE size (image 800x800 pixels)

unlike TI, unconstrained SV rectangular SE are not separable into 1-D.

Appendix: One-Scan Dilation Algorithm pseudo code

 ------------------------The following symbols denote '←' affectation, '=' equal test and '!' logical negation.

------------------------


Jean Cousty
email: jean.cousty@inria.fr

Laurent Najman

Jean Serra

Raising in watershed lattices

Keywords: Watershed, border kernels, lattice, hierarchies, data merging, mathematical morphology

come

INTRODUCTION

The watershed segmentation is a popular tool in image processing [START_REF] Beucher | Morphological segmentation[END_REF][START_REF] Vincent | Watersheds in digital spaces: An efficient algorithm based on immersion simulations[END_REF][START_REF] Bertrand | On topological watersheds[END_REF]. It is often combined with connected filters [START_REF] Salembier | Flat zones filtering, connected operators and filters by reconstruction[END_REF][START_REF] Breen | Attribute openings, thinnings and granulometries[END_REF] which simplify the considered function (also called map in this paper), and lead to watershed partitions with larger classes. However, it is not only a matter of "partition growing": a connected filter by flat zones increases the partition by flat zones of the considered function, but it does not necessarily increase its watershed partition. Indeed, the crest lines (thus, the watershed) of the map obtained after the filtering step may be "shifted" compared to the ones of the original map (see Fig. 1, where F is the initial map, F 3 a watershed of F , F 4 a filtering of F 3 , and F 5 a watershed of F 4 ).

On the other hand, there is an ambiguity in the very definition of a watershed. In some cases, we want to effectively build the crest lines (divide lines). In some other cases, we are interested in the partition of the space into its catchment basins and, in this case, we do not know on which side lies the crest lines.

Recent works by J. Cousty et al. [START_REF] Cousty | Watersheds, minimum spanning forest and the drop of water principle[END_REF][START_REF] Cousty | On watershed cuts and thinnings[END_REF] allow us to make precise these notions at least in the discrete cases of finite edge-weighted graphs. Starting from the space E made of the edges of a graph rather than of its vertices, an original idea consists of introducing the border thinnings on the edge maps. Let us briefly recall these notions. Hereafter, the workspace is a finite family of edges E (i.e., E is a set made of pairs of points), whose extremities define the vertex space E * . Any element X ⊆ E induces a family of vertices X * in E * . The class made of the maps F : E → K (where K is any finite subset of Z) is denoted by F . A map F ∈ F weights the edges of E.

F 2 = θ v (F 1 ) F 3 = θ(F )
Given a map F ∈ F, to define the border thinnings, it is convenient to also consider F * : E ⋆ → K which maps each element x in E * to the minimal value of an edge in E that contains x, i.e., F * (x) = min{F (u) | u ∈ E and x ∈ u}. Let u = {x, y} ∈ E, we say that

• u is a separating edge (for F ) if F (u) > max(F * (x), F * (y));
• u is a border edge (for F ) if F (u) = max(F * (x), F * (y)) and F (u) > min(F * (x), F * (y)); and that

• u is an inner edge (for F ) if F * (x) = F * (y) = F (u).
The border thinnings are the idempotent applications θ acting on F and generated by compositions of the elementary operator θ u , with u ∈ E:

(θ u F )(u) = min x∈u {F * (x)} if u is a border edge for F, (θ u F )(u) = F (u) if u is not a border edge for F, (θ u F )(v) = F (v) if v = u.
The map θ(F ) obtained by a border thinning is called a border kernel (of F ). These notions are illustrated in Fig 1 by the maps F, F 1 , F 2 and F 3 . Any edge of a border kernel is either an inner edge (in this case, it belongs to a minimum) or a separating edge (in this case it is not in a minimum but its extremities are both in a minimum). Therefore, any border kernel induces a (connected) partition of E * . Each class of the partition is a set of vertices induced by all edges in a single minimum. It satisfies the strong following properties.

Theorem 1 ( [6, 7]

) 1 If H is a border kernel of F then: 1/ the set S of all edges in E whose extremities are in two distinct minima of H is a watershed of F and furthermore for any u ∈ S, H(u) = F (u); and 2/ the union of all minima of H is a minimum spanning forest for F relative to the minima of F .

This theorem invites us to study the structure of the set of all border kernels. Indeed, building a hierarchy of watersheds means that we are able to define an order relation over this set.

Could we go further and build a lattice? It would then allows us to construct pyramids of watersheds and to combine watersheds stemming from several sources. How could we furthermore describe basic operations (e.g., dilations, closings) on this lattice and study their properties? We begin in this paper to investigate such structures and operators.

LATTICE OF BORDER KERNELS

Even if two distinct series {θ up } can lead to two distinct limit products θ1 and θ2 , we always have θ2 θ1 (F ) = θ1 (F ), since, by definition, there is no border edge for θ1 (F ). Similarly, θ1 θ2 (F ) = θ2 (F ). It means that the invariance domain of the thinnings θi is the same for all thinnings θi . It is the set of all border kernels, or said differently the maps in F for which there is no border edge. For the sake of simplicity, in this paper, we will only consider the border kernels whose minima are all of altitude 0 (the minimal value of K) and we will denote this set by A.

It is convenient to provide the space A with the order by minima, or min order. Let F and G ∈ A. Denote by M(F ) and M(G) the sets of all edges lying in the minima of F and G respectively, and denote by S(F ) and S(G) the sets of their separating edges (M(F )∪ S(F ) = E). Map F is said smaller than map G for the min order, written F G, if:

1/ M(F ) ⊆ M(G), or equivalently S(F ) ⊇ S(G)
2/ F (u) ≤ G(u) for any edge u ∈ S(G). 1 The reader can refer to [START_REF] Cousty | Watersheds, minimum spanning forest and the drop of water principle[END_REF] for the precise definitions of a watershed (cut) and of a relative minimum spanning forest, as considered in Theorem 1.

We call hierarchy (of border kernels) any sequence

(F 0 , . . . , F n ) of maps in A such that F 0 • • • F n .
For instance, F 3 (Fig. 1) is smaller than F 7 (Fig. 2) but it is not smaller than F 5 (Fig. 1) since S(F 3 ) S(F 5 ).

Theorem 2 (Border kernel sup-semi lattice)

The min order generates on A a sup-semi lattice, denoted , whose greatest element is the constant border kernel F max for which the weight of each edge equals to 0 (the minimal value of K).

For any two elements F 1 and F 2 in A there exists a smallest upper-bound F = F 1 F 2 , or supremum defined by:

F (u) = max{F 1 (u), F 2 (u)} if u ∈ S(F 1 ) ∩ S(F 2 ) F (u) = 0 otherwise.
Furthermore, the partition of

E * induced by F 1 F 2 is equal to the partition D = D 1 ∨ D 2 , where D 1 (resp. D 2 )
is the partition induced by F 1 (resp. F 2 ) and where ∨ denotes the supremum of partitions [START_REF] Serra | A lattice approach to image segmentation[END_REF]. This definition of a supremum in A is illustrated Fig. 3. On the other hand, two kernels F 1 and F 2 in A do not necessarily admit an infimum since the intersection of the classes induced by F 1 and F 2 can lead to classes reduced to a single vertex and there is no set of edges that induces such a class. However, if we set, by hypothesis, a smallest border kernel F 0 (the zero element of a hierarchy), then the family of all border kernels greater than F 0 is a complete lattice (for the min order). Indeed is still the supremum for this family and, since F 1 F 0 , and F 2 F 0 , the set of the border kernels less than both F 1 and F 2 is not empty and this family admits a greater elements F 1 F 2 .

Corollary 3 Let F 0 ∈ A. The family A 0 of the elements in A greater than F 0 is a complete lattice whose supremum is the one of A and whose infimum is defined by:

F (u) = min{F 1 (u), F 2 (u)} if u ∈ S(F 1 ) ∩ S(F 2 ) or u ∈ M(F 1 ) ∩ M(F 2 ) F (u) = max{F 1 (u), F 2 (u)} otherwise.

RAISINGS

Which operations can we build for acting on the sup-semi lattice of border-kernels A? The existence of a supremum orients us towards dilations, and, since this supremum extends the minima, we take as the basic operation the elementary raising ρ u defined below. Fix an edge u that parametrizes the operation ρ u and define, for any F ∈ A, the elementary raising by edges by:

ρ u F (u) = 0 if u is adjacent to two distinct minima of F ρ u F (u) = F (u) if u is adjacent to one minimum of F ρ u F (v) = F (v) if v = u Theorem 4 Let u ∈ E.
The raising ρ u acts on A and is both a closing and a dilation with respect to the min order.

In other words, ρ u satisfies the five following properties:

i) ρ u (F ) ∈ A; ii) ρ u (F 1 F 2 ) = ρ u (F 1 ) ρ u (F 2 ); iii) F 1 F 2 implies ρ u (F 1 ) ρ u (F 2 ); iv) F ρ u (F ); and v) ρ u ρ u (F ) = ρ u (F ).
Furthermore, any product ρk obtained by successive compositions of a series {ρ uj | j ∈ [1, k]} (i.e., ρk = ρ u k . . . ρ u1 ) are still dilations and closings. Thus, by iv), the successive terms { ρj } generate a hierarchy of border kernels.

FLOODINGS

The flooding operation η, which we consider in this section, is a connected operator and it constitutes an interesting way for obtaining a raising. F. Meyer and L. Najman [START_REF] Meyer | Morphologie Mathématique: approches déterministes, chapter Segmentation, arbre de poids minimum et hiérarchies[END_REF] define it as any extensive operator acting on F such that:

[(ηF )(u) > max[(ηF )(v), v adjacent to u] (1) ⇒ (ηF )(u) = F (u).
(

) 2
The previous property is rather a characteristic property than a direct definition. In fact, we can show that any flooding η can be obtained as a composition product of elementary floodings η u associated to any edge u ∈ E and defined by:

(η u F )(v) = F (v) + 1 if u and v belong to a same minimum (η u F )(v) = F (v) otherwise.
The floodings (see for instance F 4 = ηw (F 3 ) in Fig. 1), that often give nice results, as connected filters in F , do not allow for producing a hierarchy of watersheds, even if we apply them to a map in F ∈ A (see the counter-example F 5 in Figure 1). To obtain a hierarchy from floodings, we have to complete them with a class of constrained border thinnings (more details will be provided in an extended version [START_REF] Cousty | Raising and watershed lattices[END_REF]).

The maps F 3 , F 4 (Fig. 1), F 6 and F 7 (Fig. 2) illustrates the four steps of the composition of a flooding ηw with a constrained border thinning ζ, whose product is indeed a raising and thus allows for producing a hierarchy of border kernels.

Hence, to produce a hierarchy of border kernels based on floodings and constrained border thinnings, one only needs a sequence of edges to parametrize successive floodings. To this end, given an initial border kernel, one may select one edge per minimum and order these edges thanks to attributes relative to the minima (area, dynamics, . . . , [START_REF] Breen | Attribute openings, thinnings and granulometries[END_REF][START_REF] Vachier | Extraction de caractristiques, segmentation d'images et morphologie mathématique[END_REF][START_REF] Bertrand | On the dynamics[END_REF]). Then, we can construct a hierarchy of border kernels, hence a hierarchy of watersheds. These watersheds can be stacked to build a new map so that each edge is weighted by the number of watersheds it belongs to. Such a map is called a saliency map [START_REF] Najman | Geodesic saliency of watershed contours and hierarchical segmentation[END_REF]. Note that any saliency map is a border kernel.

The saliency maps S 1 and S 2 obtained from Fig. 4 thanks to dynamics and surface attributes are depicted in Fig. 5. The saliency S 1 correctly discriminates the significant contours but it also strongly delineates many small highly-contrasted regions which correspond to noise. On the other hand, S 2 does not discriminates these noisy regions, but it divides some large homogeneous zones of the image into several parts. How could we combine the advantages of these two hierarchies? The framework settled in this paper precisely provides an answer: the infimum of S 1 and S 2 (S 1 S 2 ) is depicted in Fig. 5. In the second row of Fig. 5, we also show the segmentations into 100 regions obtained from S 1 , S 2 and S 1 S 2 .

CONCLUSION

In this paper, we associated a semi lattice structure to the family of border kernels. Based on this, we showed how to merge elements of this family and proposed a generic morphological operation to build hierarchies of border kernels. We outlined how to obtain interesting hierarchies based on the flooding connected-operator and how to merge these hierarchies.

The floodings only deal with scalar functions (so that we can define minima). However, raisings are adapted to data fusion, and can handle simultaneously several minima since they are dilations. Future works will focus on this last point.

F 4 5 Fig. 1 .

 451 Fig. 1. Watershed and connected filtering in edge-weighted graphs. The map F 3 is a border kernel of F and the dashed edges a watershed of F . The map F 4 is obtained from F 3 thanks to a connected filter (a flooding see Sec. 4). The map F 5 is a border kernel of F 4 and the dashed edges is a watershed of F 4 . Minima are depicted in bold.

7 Fig. 2 .

 72 Fig. 2. Map F 6 is obtained from F 4 by a first step of thinning constrained by the kernel F 3 and F 7 is a border kernel of F 4 (obtained by a second step of thinning) constrained by F 3 .

2 G 1 G 2 Fig. 3 .

 2123 Fig. 3. Illustration of the supremum for the min order.

Fig. 4 .

 4 Fig. 4. A grayscale image.

S 1 S 2 S 1 S 2 Fig. 5 .

 12125 Fig. 5. Saliency maps obtained from Fig. 4 (first row) and associated segmentations into 100 regions (second row) [see text].