Real Time Image Processing manuscript No. (will be inserted by the editor)

Gilles Perrot^{1} · Stéphane Domas^{1} · Raphaël Couturier^{1} · Nicolas Bertaux^{2}

Fast GPU-based denoising filter using isoline levels

Received: date / Revised: date

Abstract In this study, we propose to address the issue of image denoising by means of a GPU-based filter, able to achieve high-speed processing by taking advantage of the parallel computation capabilities of modern GPUs. Our approach is based on the level sets theory first introduced by Matheron (1975) in 1975 but little implemented because of its high computation costs. What we actually do is try to guess the best isoline shapes inside the noisy image. At first, our method involved the polyline modelling of isolines; then we found an optimization heuristics which very closely fits the capabilities of GPUs. So far, though our proposed hybrid PI-PD filter has not achieved the best denoising levels, it is nonetheless able to process a 512x512 image in about 11 ms.

1 Introduction

Denoising has been a much studied research issue since electronic transmission was first used. The wide range of applications that involve denoising makes it uneasy to propose a universal filtering method. Among them, digital image processing is a major field of interest as the number of digital devices able to take pictures or make movies is growing fast and shooting is rarely done in optimal conditions. Moreover, the increase in pixel density of the CCD or CMOS sensors used to measure light intensity leads to higher noise effects and imposes high output flow rates to the various processing algorithms.

In addition, it is difficult to quantify the quality of an image processing algorithm, as visual perception is subject to high variation from one human to another.

1

FEMTO-ST institute Rue Engel Gros, 90000 Belfort, France. forename.name@univ-fcomte.fr

2

Institut Fresnel, CNRS, Aix-Marseille Université, Ecole Centrale Marseille, Campus de Saint-Jérôme, 13013 Marseille, France. nicolas.bertaux@ec-marseille.fr

So far, the advent of GPUs has brought high speedups to a lot of algorithms, and many researchers and developpers have successfully adressed the issue of implementing existing algorithms on such devices. For example in McGuire (2008),Chen et al (2009) and Sanchez and Rodriguez (2012), authors managed to design quite fast median filters. Bilateral filtering has also been successfully proposed in Yang et al (2009). Still, most high quality algorithms, like NL-means Coll et al (2011) or BM3D Dabov et al (2009) make use of non-local similarities and/or frequency domain transforms. However, speedups achieved by their current GPU implementations, though quite sigificant (as shown for example with NL-means in Palhano Xavier De Fontes et al (2010)), do not come near those achieved by local methods such as gaussian, median or neighborhood filters, as they have not originally been designed against GPU architecture. In order to fully benefit from the capabilities of GPUs, it is important that the approach to designing algorithms be more hardware-oriented, keeping in mind, from the very beginning, the intrinsic constraints of the device which is actually going to run those algorithms. Consequently, this often results in unusual options and even apparently sub-optimal solutions, but the considerable speed benefits obtained would possibly make it at least a good compromise or even the only current way to realtime high-definition image processing.

2 Contribution

As early as 1975 Matheron (1975), it was found that, under the conditions mentioned in section 5, an image can be decomposed into a set of level lines. Accordingly, reallife images fulfill the above conditions and since then, with the increase of computing capabilities, researchers have succeded in implementing such level-lines based algorithms as in Caselles et al (1997) and Monasse and Guichard (1999). A few years ago, in Bertaux et al (2004), authors proposed an original method which significantly reduces speckle noise inside coherent images, using the level lines in the image to constrain the minimization process. Those level lines are actually iso-gray-level lines, which are called isolines. In Bertaux et al (2004), isolines consist in neighborhoods of polyline shapes determined by maximum likelihood optimization. This method proved not only to bring good enhancement but also to preserve edges between regions. Nevertheless, the costs in computation time, though not prohibitive, did not allow real-time image processing; as an example, the authors of Bertaux et al (2004) managed to process an almost 2Mpixel image within a minute on an old PIII-1GHz.

Our work started by designing a set of GPU implementations with various optimization heuristics, in order to find out which tracks could be followed towards minimizing loss in quality and preserve admissible execution times. Those algorithms have been tested with reference images taken from den (2007) for which various processing results have been published. Some of the more interesting ones are listed and compared in Buades et al (2005). Statistical observations (to be detailed below) made on the output images produced by the method proposed in Bertaux et al (2004), led us to propose a very fast and simple parallel denoising method which gives good results in terms of average gray-level error, but also avoids the blurring of edges.

On the basis of the BM3D timings listed in Dabov et al (2009) and with our own measurements, our proposed GPU-based filter runs around 350 times faster and thus is able to process high definition images at over 16fps. It also achieves good denoising quality.

3 Plan

In the following, section 4 briefly focuses on recent Nvidia GPU characteristics. Section 5 will introduce the theory and notations used to define isolines. Then, in section 6, we will describe the two isoline based models that led to the final hybrid model, while section 7 details the parallel implementation of the proposed algorithm. Finally, we present our results in section 8 before drawing our conclusions and outlining our future works in section 9.

4 NVidia’s GPU architecture

GPUs are multi-core, multi-threaded processors, optimized for highly parallel computation. Their design focuses on a SIMT model (Single Instruction Multiple Threads) that devotes more transistors to data processing rather than data-caching and flow control (see CUD (2010) for more details). For example, a C2070 card features 6 GBytes global memory and a total of 448 cores bundled in several Streaming Multiprocessors (SM). An amount of shared memory, much faster than global memory, is avalaible on each SM (up to 48 KB for a C20xx card)

Writing efficient code for such architectures is not obvious, as re-serialization must be avoided as much as possible. Thus, code design requires one pays attention to a number of points, among which:

– the CUDA model organizes threads by a) thread blocks in which synchronization is possible, b) a grid of blocks with no possible synchronization between them.

– there is no way to know how blocks are scheduled during one single kernel execution.

– data must be kept in GPU memory, to reduce the overhead generated by copying between CPU and GPU.

– the total amount of threads running the same computation must be as large as possible.

– the number of execution branches inside one block should be as small as possible.

– global memory accesses should be coalescent, i.e. memory accesses done by physically parallel threads (2 x 16 at a time) must be consecutive and contained in a 128 Bytes range.

– shared memory is organized in 32x32 bit-wide banks. To avoid bank conflicts, each parallel thread (2 x 16 at a time) must access a different bank.

All the above characteristics always make designing efficient GPU code all the more constraining as nonsuited code would probably run even slower on GPU than on CPU.

5 Isolines

In the following, let I be the reference noiseless image (assuming we have one), I′ the noisy acquired image corrupted by Independent and Identically Distributed (IID) additive white gaussian noise of zero mean value and standard deviation σ. Let Î be the denoised image. Each pixel of I′ of coordinates (i, j) has its own gray level z(i, j).

As introduced above and since most common images are continuous and contain few edges, they can be decomposed into a set of constant gray level lines called isolines. Then our goal is to find, for each single pixel of a noisy image, the isoline it belongs to. The generalized likelihood criterion (GL) is used to select the best isoline among all the considered ones, all of which must have the same number of pixels in order to be compared.

5.1 Fixed-length isolines

For each pixel (i, j) of the corrupted image, we look for the gray level of the isoline it belongs to, inside a rectangular window ω centered on (i, j). Inside ω, let S^{n} be the isoline part which the center pixel belongs to. S^{n} is a set of n pixel positions (i_{q}, j_{q}) (q ∈ [0; n[). The gray levels z along S^{n} follow a gaussian probability density function whose parameters μ_{S^{n}} (mean value of isoline part) and σ (standard deviation brought by gaussian noise ) are unknown. Let \overline{S^{n}} be defined by ω = S^{n} ∪ \overline{S^{n}}. For each pixel, the mean values μ_{ij} of gray levels z over \overline{S^{n}} are unknown and supposed independant . Let Z be the gray levels of pixels in ω and \{μ_{ij}\}\overline{S^{n}} the mean values of pixels in \overline{S^{n}}. The likelihood is given by:

P [Z|S^{n}, μ_{S^{n}}, \{μ_{ij}\}\overline{S^{n}} , σ]

When separating contributions from regions S^{n} and \overline{S^{n}}, it becomes:

∏_{(i,j)∈S^{n}} P [z(i, j)|μ_{S^{n}}, σ]. ∏_{(i,j)∈\overline{S^{n}}} P [z(i, j)| \{μ_{ij}\}\overline{S^{n}} , σ](1)

The goal is then to estimate the value of the above expression, in order to find the boundaries of S^{n} that maximize expression (1). Let us consider that, on \overline{S^{n}}, the values z(i, j) are the likelihood estimations μ̂_{ij} for μ_{ij}. The second term of expression (1) becomes:

∏_{(i,j)∈\overline{S^{n}}} P [z(i, j)| \left\{{ \atop μ_{ij}}\right\}\overline{S^{n}} , σ] = 1  (2)

which leads to the generalized likelihood expression:

∏_{(i,j)∈S^{n}} P [z(i, j)|μ_{S^{n}}, σ]  (3)

As we know the probability density function on S^{n}, (3) can then be developped as

∏_{(i,j)∈S^{n}} \frac{1}{\sqrt{2πσ^{2}}} e^{− \frac{(z(i,j)−μ_{S^{n}})^{2}}{2σ^{2}}}  (4)

The log-likelihood is then given by:

−\frac{n}{2}log (2π) − \frac{n}{2}log (σ^{2}) − \frac{n}{2}  (5)

inside which the vector of parameters (μ_{S^{n}}, σ) is determined by maximum likelihood estimation

( \begin{matrix}μ_{S^{n}} = \frac{1}{n} ∑_{(i,j)∈S^{n}} z(i, j) \\ σ̂^{2} = \frac{1}{n} ∑_{(i,j)∈S^{n}} (z(i, j) − μ_{S^{n}})^{2} \end{matrix}

The selection of the best isoline is done by searching which one maximizes the expression of equation (5).

5.2 Lengthenable isolines

Searching for larger isolines should lead to better filtering as a larger number of pixels would be involved. However, processing all possible isolines starting from each pixel would be too costly in computing time, even in the case

Fig. 1 Determination and lengthening of an isoline: The gray level z of each pixel is seen as an elevation value. S^{n} is the n pixel length isoline for pixel of coordinates (i, j). The elongation of S^{n} by S^{p} (p pixel length) is submitted to the GLRT condition (see eq. (8)).

of a small GPU-processed 512x512 pixel image. Therefore, we chose to build large isolines inside an iterative process including a mandatory validation stage between each lengthening iteration, so as to reduce the number of pixel combinations to be examined and keep the estimation of deviation σ within a satisfactory range of values.

Let S^{n} be a previously selected isoline part and S^{p} connected to S^{n} in such a way that S^{p} could be seen as an addition to S^{n} so as to define a possible valid isoline S^{n+p}. Figure 1 illustrates this situation with a very simple example image. In this figure, the gray level of each pixel is used as its corresponding height (z) in order to visualize isolines easily. Some of the orthogonal isoline projections have been drawn in dotted line in the (i, j) plane. Both labeled parts S^{p} and S^{n} are represented in the (i, j) plane and in the 3D associated plot.

In order to decide whether S^{n+p} can be considered as an actual isoline, we compare the log-likelihood of both hypothesis below by using GLRT (Generalized Likelihood Ratio Test):

First, assuming that S^{n+p} is an isoline, the gray levels of its pixels share the same mean value μ_{n+p}. According to (5), its log-likelihood is

−\frac{(n + p)}{2} (log (2π) + 1) − \frac{(n + p)}{2}log (σ̂_{1}^{2})  (6)

where σ̂_{1} is the estimation of the standard deviation along S^{n}.

Second, considering S^{n} and S^{p} as two separate isoline parts connected together, the gray levels of their pixels have two different mean values μ_{n} and μ_{p}. The log-likelihood is the sum of both log-likelihoods, given by

−\frac{(n + p)}{2} (log (2π) + 1) − \frac{n}{2}log (σ̂_{2}^{2}) − \frac{p}{2}log (σ̂_{2}^{2}) (7)

where σ̂_{2} is the estimation of the standard deviation along S^{n} and S^{p}.

The difference between (6) and (7) leads to the expression of GLRT (S^{n+p}, S^{n}, S^{p}, T_{max}):

T_{max} − (n + p). [log (σ̂_{1}^{2}) − log (σ̂_{2}^{2})]  (8)

The decision to validate lengthening from S^{n} to S^{n+p} depends whether GLRT (S^{n+p}, S^{n}, S^{p}, T_{max}) is higher or lower than 0. Value T_{max} is the GLRT threshold.

6 Isoline models

The most obvious model considers isolines as polylines. Each isoline can then be curved by allowing a direction change at the end of each segment; we shall call such isolines poly-isolines.

In order to keep the number of candidate isolines within reasonable range, we chose to build them by combinating segments described by simple pre-computed patterns. Each pattern p_{l,d} describes a segment of length l and direction d. For one given l value, all p_{l,d} patterns are grouped into a matrix denoted P_{l}. Figure 8 shows an example of such a pattern matrix for l = 5.

To fit the GPU-specific architecture, we define regularly distributed D primary directions (D = 32 in our examples).

6.1 Poly-isolines with limited deviation angle (PI-LD)

At one stage we implemented an algorithm parsing the tree of all possible polyline configurations, but the process proved far too slow regarding our goal, even on GPU, because of the amount of memory involved (and consequent memory accesses) and because of the necessary reduction stage for which GPU efficiency is not maximum. So we focused on a variant inspired by Bertaux et al (2004) in which the selected direction of the next segment depends on the whole of the previously built and validated poly-isoline.

Let us consider a poly-isoline S^{n} under construction, starting from pixel (i, j) and made of K validated segments s_{k} (k ∈ [1; K]) of length l, each of them having its own direction d_{k}. The coordinates of the ending pixel of each segment s_{k} are denoted (i_{k}, j_{k}). Both of the following sums

C_{x} (Z(S^{n})) = ∑_{(i,j)∈S^{n}} z(i, j)  (9)

and C_{x^{2}} (Z(S^{n})) = ∑_{(i,j)∈S^{n}} z(i, j)^{2}  (10)

have been obtained during the previous lengthening steps.

Let us examine now how to decide wether to add a new segment to S^{n} or to stop the lenghtening process. The main idea is to apply each pattern p_{l,d} to the ending pixel (i_{k}, j_{k}), on the condition that its direction is contained within the limits of maximum deviation Δd_{max}.

Maximum deviation Δd_{max} prevents poly-isolines from beeing of circular shape (or backward-oriented) which would possibly generate supplementary artefacts in the output image. Another of its benefits is to reduce the number of combinations to be evaluated.

For each allowed pattern, GLRT is performed in order to decide if the corresponding segment could likely be added to the end of the poly-isoline S^{n}. If none is validated by GLRT, the poly-isoline S^{n} is stopped.

If at least one segment has been accepted by GLRT, the one that leads to the maximum likelihood (ML) value of the lengthened poly-isoline S^{n+l} is selected and integrated to S^{n+l} as s_{K+1}.

In order to avoid critical situations where the first selected segment would not share the primary direction of the actual poly-isoline, no selection is performed on the level of the first segment; D poly-isolines are kept and submitted to the lengthening process. To ensure isotropy, each of them shares the direction of one pattern p_{l,d} (d ∈ [0; D]).

Eventually, the poly-isoline with the maximum likelihood value is selected among the longest ones.

Figure 2 illustrates one stage of the lengthening process with the example of a two-segment poly-isoline at the beginning of stage (l = 5 and Δd_{max} = 2).

6.2 Poly-isolines with precomputed directions (PI-PD)

Though much faster, the PI-LD-based filter may be considered a bit weak compared to state-of-the-art filters like BM3D family algorithms Dabov et al (2009). Furthermore, we saw that this way of building poly-isolines requires the alternate use of two different types of validation at each lengthening stage: GLRT and maximum likelihood minimization. In order to be performed, each of them generates numerous branches during kernel execution, which does not fit GPU architecture well and leads to execution times that we hoped would be more impressive.

Within the PI-LD model, at each pixel (i, j), as no selection is done at the first stage, D poly-isolines are computed and kept as candidate though, obviously, only one follows the actual isoline at (i, j). So, if we assume we can achieve a robust determination of the direction at any given pixel of this isoline, it becomes unnecessary to perform the selection at each lenghtening step. Thus, at each pixel (i, j), only the first segment has to be determined in order to obtain the local direction of the isoline. This leads to an important reduction of the work complexity: the above PI-LD model needs to evaluate D. (2.Δ_{dmax} + 1)^{K−1} segments at each pixel position, while only D.K evaluations are needed in the second case. For example, with a maximum of K = 5 segments and a maximum deviation of Δ_{dmax} = 2, the PI-LD needs to evaluate up to 20000 segments per pixel where only 160 should be enough.

(a) Isoline with two validated segments s_{1} and s_{2}.

(b) First evaluated segment, corresponding to pattern p_{5,0}.

(c) Second evaluated segment, corresponding to pattern p_{5,1}.

(d) Third evaluated segment, corresponding to pattern p_{5,2}.

(e) Fourth evaluated segment, corresponding to pattern p_{5,3}.

(f) Fifth evaluated segment, corresponding to pattern p_{5,4}.

Fig. 2 Example of lengthening process starting with a twosegment poly-isoline (l = 5, Δd_{max} = 2). The initial situation is shown in 2a, while 2b to 2f represent the successive candidate segments. The direction index of the last validated segment is d_{2} = 2 (2a). It implies that direction indexes allowed for the third segment range from d_{2} − Δdmax = 0 to d_{2} + Δdmax = 4 (2b to 2f). The lengthening of the poly-isoline is accepted if at least one segment has a positive GLRT. If there are several, the one which minimizes the standard deviation of the whole poly-isoline is selected.

On the basis of these observations, we propose a new model that we shall call PI-PD, that completely separates the validation stages performed in the PI-LD model implementation mentioned above:

A first computation stage selects the best first segment s_{1} starting at each pixel (i, j) of the input image. Its direction index d_{1}(i, j) is then stored in a reference matrix denoted I_{Θ}; sums C_{x} and C_{x2} along s_{1}(i, j) are also computed and stored in a dedicated matrix I_{Σ}. It can be noticed that this selection method of s_{1} segments is a degraded version of PI-LD constrained by K = 1.

A second stage manages the now independant lengthening process. For one given state of a poly-isoline where the last added segment has been s_{K}, the pattern whose direction index is given by d = I_{Θ}(i_{K}, j_{K}) defines the only segment to be evaluated. Both corresponding sums C_{x} and C_{x2} are read from matrix I_{Σ} and used in GLRT evaluation. The last point is to prevent poly-isolines from turning back.

Figure 3 details this process, starting from the same initial state as in figure 2 with the noticeable difference that no deviation limit is needed.

Thus, as introduced above, work complexity is considerably reduced, as each pattern is only applied once at one given pixel (i, j), and associated values are computed only once; they are re-used every time one polyisoline’s segment ends at pixel (i, j). Also, this fits GPU constraints better, as it avoids multiple branches during kernel execution. It remains that, the building of polyisolines is done without global likelihood optimization.

Eventually, the model has been improved by adding to it the ability to thicken poly-isolines from one pixel up to three which allows to achieve higher PSNR values by increasing the number of pixels of poly-isolines in addition to the lengthening process. This may apply to large images which do not contain small relevant details, as it may blur small significant details or objects present in the noisy image. Still, this feature makes PI-PD more versatile than our reference BM3D, which has prohibitive computation times when processing large images (over 5 minutes for a 4096x4096 pixel image) and thus should require a slicing stage prior to processing them, causing some overhead.

6.3 Hybrid PI-PD

As the determination of each segment’s direction only involves a few pixels, the PI-PD model may not be robust enough in regions where the surface associated with Z has a low local slope value regarding power of noise σ^{2}. We shall call those regions Low Slope Regions (LSR). Figure 4 shows this lack of robustness with an example of two drawings of additive white gaussian noise applied on the same reference image (Figure 6). Within this image, we focused on a small 11x11 pixel window containing two LSR with one sharp edge between them.

Figures 4d and 4e show that the directions computed by PI-PD are identical from one drawing to the other near the edge (lines 5-7), while they vary in LSR (lines 1-4, 8-11).

Within such regions, our speed goals forbid us to compute isoline directions with the PI-LD model, more robust but far too slow. Instead we propose a fast solution which implies designing an edge detector whose principle is to re-use the segment patterns defined in section 6 and to combine them by pairs in order to detect any possible LSR around the center pixel. If a LSR is detected, the output gray-level value is the average value computed on the current square window, otherwise, the PI-PD output value is used.

In order to further simplify computation, only the patterns that do not share any pixel are used. These patterns have a direction which is a multiple of 45^{∘}.

Each base direction (Θ_{i}) and its opposite (Θ_{i}+π) [2π] define a line that separates the square window in two

(a) Poly-isoline with two validated segments.

(b) Next direction is read from element (i_{2}, j_{2}) of I_{Θ}.

(c) Pattern p_{l,d_{3}} is then applied at (i_{2}, j_{2}) and GLRT is performed. Both sums needed to perform GLRT are read from element (i_{2}, j_{2}) of I_{Σ}.

(d) If accepted by GLRT, segment s_{3} is added to polyisoline.

Fig. 3 Example of PI-PD lengthening process starting with a two-segment poly-isoline (l = 5). The initial situation is represented in 3a, while 3a to 3d represent the successive processing steps. The end pixel of the last validated segment is (i_{2}, j_{2}) (3a). Reference matrices I_{Θ} and I_{Σ} provide the values needed to select the pattern to be applied on (i_{2}, j_{2}) (3b and 3c). GLRT is performed to validated lengthening or not. This process goes on until one submitted segment does not comply with GLRT.

regions (top and and bottom regions, denoted T and B). We assume that segments on the limit belong to the T region which includes pixels of orientation from Θ_{i} to Θ_{i} + π. This region comprises three more segments of directions (Θ_{i} + \frac{π}{4} ), (Θ_{i} + \frac{2π}{4} ) and (Θ_{i} + \frac{3π}{4} ). The other region (B) only includes three segments of directions (Θ_{i} + \frac{5π}{4} ), (Θ_{i} + \frac{6π}{4} ) and (Θ_{i} + \frac{7π}{4} ).

Figure 5 illustrates this organization for Θ_{i} = Θ_{4} = 45^{∘}. Each bar represents a pixel in the detector’s window. Pixels with null height are not involved in the GLRT. Pixels represented by higher bars define the T region and those represented by shorter bars define the B region.

For each Θ_{i}, one GLRT is performed in order to decide whether the two regions T and B defined above are likely to be seen as a single region or as two different

(a) Reference image

(b) Image corrupted by random drawing n^{∘}1 (c) Image corrupted by random drawing n^{∘}2

(d) Isoline directions for random drawing n^{∘}1 (e) Isoline directions for random drawing n^{∘}2

Fig. 4 Zoom on a small square window of the airplane image. 4a reproduce the zoom on the window, taken from the reference image of Figure 6. 4b, 4c and 4a and are 3D views where each bar represents a pixel whose gray-level corresponds to the height of the bar. Figures 4d and 4e are 2D top views of the window. The chosen window shows an edge between two regions of low slope. The images 4b and 4c are corrupted with two different random drawings of the same additive white gaussian noise (AWGN) of power σ^{2} and mean value 0. 4d and 4e show, for each pixel of the window, the direction of the isoline found by PI-PD. In regions of low slope (the two regions at the top and the bottom), the determination of the direction is not robust. But near the edge, directions do not vary from one drawing to another.

ones, separated by an edge as shown in figure 5. The center pixel is located on the edge. Equations (6), (7) and (8) lead to a similar GLRT expression:

T 2_{max} − (8.l + 1). [log (σ̂_{3}^{2}) − log (σ̂_{4}^{2})]  (11)

where σ_{3} is the standard deviation considering that the two regions are likely to define a single one and σ_{4} the standard deviation if an edge is more likely to separate the two regions. T 2_{max} is the decision threshold. With equation (11), a negative result leads to an edge detection, oriented towards direction Θ_{i}. When GLRT is known for each Θ_{i}, we apply the following hybridation policy:

Fig. 5 Edge detector. 3D view representing an example square 11x11 pixel window (l = 5) used in the edge detector for Θ_{4} = 45^{∘} around a center pixel colored in black. Each pixel is represented by a bar. Bars of height value 0 are for pixels that are not involved in the detector. Top region is defined by five pattern segments and includes the center pixel. Bottomp region only includes three pattern segments. The different height values are meant to distinguish between each of the three different sets of pixels and their role.

a) more than one negative GLRT: the PI-PD output value is used.

b) only one negative GLRT: the center pixel is likely to be on a well-defined edge, and only the region it belongs to is considered. The average value of its pixel gray levels is then used.

c) no negative GLRT: the window around the center pixel is likely to be a LSR. The average value on the whole square window is used (11x11 pixels in the example of Figure 5).

(a) Reference noiseless airplane image (b) Location of the example window in the reference image.

Fig. 6 Location of the example window inside the reference image. Figure 6a shows the whole reference image and 6b zooms on the part where the example 11x11 pixel window is.

It must be noticed that point b) has been introduced in order to achieve smoother transitions between regions to which PI-PD is applied and those in which the plain average value is used. Figure 7 shows an example of such a classification achieved by the edge detector. The detector has been applied on the top noisy airplane image with a GLRT threshold value T 2_{max} = 2. Black pixels represent pixel classified as on an edge, while white ones are those which belong to LSR.

(a) Noisy airplane image (b) Pixel classification performed by the edge detector.

Fig. 7 Pixel classification inside the noisy image. Figure 7a shows the noisy input image and 7b reproduces the output classification of pixels, as a black and white image, obtained with threshold value T 2_{max} = 2. Black pixels are supposed to be near an edge, while white pixels belong to Low Slope Regions.

7 Hybrid PI-PD filter Implementation: details

All implementation details that will be given here are relative to the proposed PI-PD models and Nvidia^{\begin{matrix}c \\ ◯ \end{matrix}} GPU devices.

7.1 Segment patterns

The first kernel to be run is kernel genPaths() which generates matrix P_{l}. Its elements (Δi; Δj) are the relative coordinates of the pixels which define segment patterns p_{l,d}. The dimensions of matrix P_{l} are D rows × l columns. To fit GPU architecture as closely as possible, we chose D = 32 patterns. Each segment s_{k} of a polyisoline can then be seen as a pattern p_{l,d} applied on the starting pixel (i, j) of this segment, denoted p_{l,d}(i, j).

The example in figure 8 shows the first quarter of P_{5} and the corresponding eight discrete segment patterns in the first quadrant. The three remaining quarters of the matrix are easily deduced by applying successive rotations of angle \frac{π}{2} to the above elements.

7.2 Generation of reference matrices I_{Σ} and I_{Θ}

In order to generate both matrices, a GPU kernel kernel precomp() computes, in parallel for each pixel (i, j):

– the direction δ of the most likely segment s_{1} = p_{l,δ}(i, j) among the D possible ones. This value is stored in matrix I_{Θ} at position (i, j).

– values C_{x}(s_{1}) and C_{x^{2}}(s_{1}) defined in equations (9) and (10). This vector of values is stored in matrix I_{Σ} at position (i, j).

In order to reduce processing time, the input image is first copied into texture memory (see algorithm 1 for

P5 = [\begin{matrix}(0, 1) & (0, 2) & (0, 3) & (0, 4) & (0, 5) \\ (0, 1) & (0, 2) & (−1, 3) & (−1, 4) & (−1, 5) \\ (0, 1) & (−1, 2) & (−1, 3) & (−2, 4) & (−2, 5) \\ (−1, 1) & (−1, 2) & (−2, 3) & (−3, 4) & (−3, 5) \\ (−1, 1) & (−2, 2) & (−3, 3) & (−4, 4) & (−5, 5) \\ (−1, 1) & (−2, 1) & (−3, 2) & (−4, 3) & (−5, 3) \\ (−1, 0) & (−2, 1) & (−3, 1) & (−4, 2) & (−5, 2) \\ (−1, 0) & (−2, 0) & (−3, 1) & (−4, 1) & (−5, 1) \\ … & … & … & … & … \end{matrix}]

Fig. 8 Top: example segment patterns p_{5,d} for d ∈ [0; 7]; the black pixel represents the center pixel (i, j), which does not belong to the pattern. The gray ones define the actual pattern segments. Bottom: the first 8 lines of corresponding matrix P_{5} whose elements are the positions of segment pixels with respect to the center pixel.

initializations and memory transfer details), thus taking advantage of the 2D optimized caching mechanism.

This kernel follows the one thread per pixel rule. Consequently, each value of P_{l} has to be accessed by every thread of a block. That led us to load it from texture memory first, then copy it into all shared memory blocks. This has proved to be the fastest scheme.

Algorithm 2 summarizes the computations achieved by kernel precomp(). Vector (C_{x}, C_{x2}) stores the values of C_{x}(s_{1}) and C_{x^{2}}(s_{1}) associated with the current tested pattern. Vector (C_{x−best}, C_{x2−best}) stores the values of C_{x}(s_{1}) and C_{x^{2}}(s_{1}) associated with the best previously tested pattern.

In the same manner, σ and σ_{best} are deviation values for current and best tested patterns.

The selection of the best pattern is driven by the value of the standard deviation of candidate isolines. Lines 2 and 3 compute both sums for the first pattern to be evaluated. Line 4 computes its standard deviation. Then, lines 5 to 14 loop on each pattern and keep values associated with the best pattern found. These values are eventually stored in matrices I_{Θ} and I_{Σ} on lines 16 and 17.

7.3 PI-PD lengthening process: kernel PIPD()

This parallel kernel is run in order to obtain the image of the isolines. It is detailed in algorithm 3, (see section 6.2 for process description).

Lines from 2 to 11 perform allocations for the first lengthening to evaluate. More precisely, (i_{1}, j_{1}) represents the starting pixel of the current segment; (i_{2}, j_{2})

Algorithm 1: Initializations in GPU memory

1: l ← step size;

2: D ← number of primary directions;

3: I_{n} ← noisy image;

                                                       4: I_{ntex} ← I_{n};  /* copy to texture mem. */

5: P_{l} ← kernel genPaths ; /* pattern matrix */

                                                       6: P_{ltex} ← P_{l};  /* copy to texture mem. */

7: T_{max} ← GLRT threshold (lengthening); 8: T 2_{max} ← GLRT threshold (edge detection);

Algorithm 2: generation of reference matrices,

                                                  kernel kernel precomp()

1: foreach pixel (i, j) do /* in parallel */

2: C_{x−best} ← ∑_{(y,x)∈p_{l,0}(i,j)} I_{ntex}(i + y, j + x) ;

3: C_{x2−best} ← ∑_{(y,x)∈p_{l,0}(i,j)} I_{ntex}^{2}(i + y, j + x) ;

4: σ_{best} ← standard deviation along p_{l,0}(i, j) ;

                                                             /* loop on each pattern                   */

5: foreach d ∈ [1; D − 1] do

6: C_{x} ← ∑_{(y,x)∈p_{l,d}(i,j)} I_{ntex}(i + y, j + x);

7: C_{x2} ← ∑_{(y,x)∈p_{l,d}(i,j)} I_{ntex}^{2}(i + y, j + x);

8: σ ← standard deviation along p_{l,d}(i, j);

9:  if σ_{d} < σ_{best} then  /* keep the best */
10:  C_{x−best} ← C_{x} ; ; 11:  C_{x2−best} ← C_{x2} ; ; 12:  Θ_{best} ← d ;

13: end 14: end

15:  I_{Σ}(i, j) ← [C_{x−best}, C_{x2−best}] ;  /* stores */
16:  I_{Θ}(i, j) ← Θ_{best} ;  /* in matrices */

17: end

is both its ending pixel and the starting pixel of the next segment; d_{1} and d_{2} are their directions, read from precomputed matrix I_{Θ}. C_{x}^{1} and C_{x2}^{1} are the gray-level sums along the current poly-isoline; C_{x}^{2} and C_{x2}^{2} are the gray-level sums of the candidate segment. The current poly-isoline ends at (i_{1}, j_{1}) and is made of l_{1} pixels (already accepted segments); its standard deviation is σ_{1}. The loop extending from lines 12 to 21 performs the allocations needed to proceed one segment forward, as long as GLRT is true. If the lengthening has been accepted, the length of the poly-isoline is updated in line 13, and the same is done with C_{x} and C_{x2} which are read from precomputed matrix I_{Σ} (see equations (9) and (10) for definition). Finally, using direction value d_{2}, it translates the coordinates (i_{1}, j_{1}) to the end of the newly elongated poly-isoline, and (i_{2}, j_{2}) to the end of the next segment to be tested. As soon as the GLRT condition becomes false, line 23 eventually produces the output value of the denoised image at pixel (i, j), that is, the average graylevel value along the poly-isoline.

Algorithm 3: PI-PD lengthening process

 kernel PIPD()

1: foreach pixel (i, j) do /* in parallel */

2:  (C_{x}^{1}, C_{x2}^{1}) ← z(i, j) ;  /* starting pixel */

 3:  (i_{1}, j_{1}) ← (i, j) ;  /* first segment */

4:  (C_{x}^{1}, C_{x2}^{1}) ← I_{Σ}(i_{1}, j_{1}) ;  /* read matrix */

 5:  d_{1} ← I_{Θ}(i, j) ;  /* read matrix */
 6:  l_{1} ← l ;  /* isoline length */

7: σ_{1} ← (C_{x2}^{1}/l_{1} − C_{x}^{1})/l_{1}; 8: (i_{2}, j_{2}) ← end of first segment; 9: (C_{x}^{2}, C_{x2}^{2}) ← I_{Σ}(i_{2}, j_{2}) ; /* 2^{nd} segment */ 10: d_{2} ← I_{Θ}(i_{2}, j_{2}); 11: σ_{2} ← (C_{x2}^{2}/l − C_{x}^{2})/l ; 12: while GLRT (σ_{1}, σ_{2}, l_{1}, l) < T_{max} do

13:  l_{1} ← l_{1} + l ;  /* lengthening */

14: (C_{x}^{1}, C_{x2}^{1}) ← (C_{x}^{1}, C_{x2}^{1}) + (C_{x}^{2}, C_{x2}^{2});

15:  σ_{1} ← (C_{x2}^{1}/l_{1} − C_{x}^{1})/l_{1} ;  /* update */
16:  (i_{1}, j_{1}) ← (i_{2}, j_{2}) ;  /* step forward */

17: d_{1} ← d_{2}; 18: (i_{2}, j_{2}) ← end of next segment;

/* next segment */ (C_{x}^{2}, C_{x2}^{2}) ← I_{Σ}(i_{2}, j_{2});

19: d_{2} ← I_{Θ}(i_{2}, j_{2}); 20: σ_{2} ← (C_{s2}^{2}/l − C_{s}^{2})/l ; 21: end 22: end 23: Î(i, j) ← C_{x}^{1}/l_{1} ; /* isoline value */

7.4 Hybrid PI-PD : kernel edge detector()

As introduced in section 6.3, the aim of kernel kernel edge detector() is to divide pixels into two classes according to their belonging to a LSR or not. Algorithm 4 explains the detailled procedure. Lines 2 to 6 initialize values of the direction index (Θ), the number of edges detected (edgeCount), the gray-level sum along the pixels that defines the H half-plane (sumEdge) and the number of pixels that defines both half-planes H and L (nH, nL). Then the loop starting at line 7 performs the GLRT for every considered direction index Θ. Values sumH and sumL are vectors of two parameters x and y, parameter x being the sum of gray-level values and y the sum of square gray-level values. Value sumH is computed along the pixels of halfplane H and is obtained by loop at lines 10 to 14; Value sumL is computed along the pixels of half-plane L and is obtained by loop at lines 15 to 19. Value I_{ntex}(i, j) refers to the gray-level value at pixel (i,j) previously stored in texture memory. Eventually, the isoline level value is output at line 27, 30 or 33 depending on the situation (see 6.3 for details about the decision process).

8 Results

The proposed hybrid PI-PD model has been evaluated with the 512x512 pixel sample images used by den (2007) in order to make relevant comparisons with other filtering techniques. As we aim to address image processing Algorithm 4: edge detector and pixel classifier

                                                  kernel edge detector()

1: foreach pixel (i, j) do /* in parallel */

                                                       2:  Θ ← 0;  /* direction index */

3: edgeCount ← 0; 4: sumEdge ← 0; 5: nH ← 5l + 1; 6: nL ← 3l; 7: while (Θ < 32) do 8: sumH ← (I_{ntex}(i, j), I_{ntex}^{2}(i, j)); 9: sumL ← (0, 0); 10: for (α = Θ to α = Θ + 16 by step 4) do

11:  sP at ← ∑_{(y,x)∈P_{l,α}(i,j)} I_{ntex}(i + y, j + x); ; 12:  sP at2 ← ∑_{(y,x)∈P_{l,α}(i,j)} I_{ntex}^{2}(i + y, j + x); ; 13:  sumH ← sumH + (sP at, sP at2);

14: end

15: for (α = Θ + 20 to α = Θ + 28 by step 4) do

16:  sP at ← ∑_{(y,x)∈P_{l,α}(i,j)} I_{ntex}(i + y, j + x); ; 17:  sP at2 ← ∑_{(y,x)∈P_{l,α}(i,j)} I_{ntex}^{2}(i + y, j + x); ; 18:  sumL ← sumL + (sP at, sP at2);

19: end

20: if (GLRT (sumH, nH, sumL, nL) > T 2_{max})

then

{21:  edgeCount ← edgeCount + 1; \atop 22:  sumEdge ← sumH.x;}

23: end

24: Θ ← Θ + 4;

25: end

                                                             /* outputs isoline value                  */

26: if (edgeCount == 0) then

27:  Î(i, j) ← \frac{(sumH.x + sumL.x)}{nH + nL} ; /* LSR */

28: end

29: if (edgeCount == 1) then

30:  Î(i, j) ← \frac{(sumEdge)}{nH}

31: end

32: if (edgeCount > 1) then

33: Î(i, j) ← Î_{P IP D}(i, j); /* PI-PD */ 34: end 35: end

in very noisy conditions (as in Perrot et al (2011)), we focused on the noisiest versions, degraded by AWGN of standard deviation σ = 25.

Quality measurements of the denoised images in comparison with reference images have been obtained by the evaluation of:

a) Peak Signal to Noise Ratio (PSNR) that quantifies the mean square error between denoised and reference images: MSE(I, Î). We used the following expression:

P SNR = 10.log_{10} (\frac{max(Î)}{MSE(I, Î)})

PSNR values are given in dB and highest values mean best PSNR.

b) The Mean Structure Similarity Index (MSSIM, defined in Wang et al (2004)), which quantifies local similarities between denoised and reference images inside a sliding window. MSSIM values belong to an interval [0; 1]; the closer to 1 the better.

PSNR is widely used to measure image quality but can be misleading when used by itself: as demonstrated in Wang et al (2004), the processing of noisy images can bring a high PSNR value but very bad visual quality. This is avoided by the use of the MSSIM index along with the PSNR value: when both of them show high values, the overall quality can be considered high.

Result figure 9 provides the PSNR and MSSIM of every image, denoised with three different filters: average 5x5, hybrid PI-PD and BM3D. The noisy column shows the values for each image before denoising. BM3D (Dabov et al (2009)) is taken as a reference in terms of denoising quality, while the average filter is taken as a reference in terms of processing time. The window size of 5x5 pixels has been choosen to achieve PSNR values similar to those obtained by PI-PD.

BM3D code is run on a quad-core Xeon E31245 at 3.3GHz and 8GByte RAM under linux kernel 3.2 (64bits), while PI-PD as well as average filter codes is run on a Nvidia C2070 GPU hosted by a PC running linux kernel 2.6.18 (64bits). The average filter used is an efficient parallel GPU implementation that we developped. It is a generic and versatile separable convolution kernel that outputs more than 700MPixels per second in the 5x5 averaging configuration.

Hybrid PI-PD measurements were performed with n = 25, l = 5, T_{max} = 1 and T 2_{max} = 2. BM3D measurements have been performed with the freely available BM3D software proposed in Dabov et al (2009).

The hybrid PI-PD model proves much faster than BM3D and better than the average 5x5 filter. Processing the thirteen images of the database reveals that hybrid PI-PD brings an average improvement of 1.5dB (PSNR) and 7.2% (MSSIM) against the average filter at the cost of 35 times its computational duration. Against hybrid PI-PD, BM3D achieves an average improvement of 2.4dB and 4.6% at the cost of 350 times as much duration. Actually, the 5x5 average filter takes around 0.35 ms to process an image while hybrid PI-PD needs around 11 ms and BM3D around 4.3 s.

It must be noticed that experimental optimization show that the vector of parameter values T_{max} = 1 and T 2_{max} = 2 is optimal for 11 of the 13 images of the database. Better results are obtained with a slightly different value of T 2_{max} for peppers or zelda whose denoised

Image  Noisy average hybrid BM3D
5x5 PI-PD
airplane 19.49dB 26.39dB 28.46dB 30.88dB
0.58 0.84 0.88 0.93
barbara 20.04dB 22.76dB 24.26dB 30.60dB
0.70 0.76 0.83 0.94
boat  20.33dB 25.58dB 27.54dB 30.02dB
0.66 0.81 0.87 0.91
couple  20.28dB 25.25dB 27.33dB 29.77dB
0.69 0.79 0.87 0.91
elaine  19.85dB 28.71dB 28.94dB 30.60dB
0.59 0.86 0.87 0.91
fingerprint 20.34dB 23.33dB 26.07dB 27.93dB
0.93 0.87 0.95 0.96
goldhill  19.59dB 26.47dB 27.43dB 29.22dB
0.67 0.82 0.87 0.88
lena  19.92dB 27.99dB 29.14dB 31.80dB
0.60 0.84 0.88 0.93
man  20.38dB 24.74dB 26.74dB 28.14dB
0.71 0.80 0.86 0.87
mandrill 19.34dB 20.34dB 22.38dB 24.75dB
0.77 0.69 0.83 0.88
peppers 19.53dB 27.30dB 28.68dB 30.87dB
0.61 0.86 0.87 0.92
stream  20.35dB 23.23dB 25.35dB 26.34dB
0.80 0.78 0.87 0.88
zelda  17.71dB 23.13dB 27.71dB 30.49dB
0.58  0.87  0.88  0.93

Fig. 9 Comparison between hybrid PI-PD, average and BM3D filters. PI-PD parameter values: n = 25, l = 5, T_{max} = 1 and T 2_{max} = 2. The noisy column correspond to the noisy input images, before denoising. Timings: average filter in around 0.35 ms hybrid PI-PD in around 11.0 ms and BM3D in around 4.3 s

images can obtain a MSSIM index of 0.90. Most of the computational time of hybrid PI-PD is spent by the edge detector, which clearly does not fit GPU requirements to achieve good performance. For information, the simple PI-PD model runs in less than 4 ms in the same conditions.

Figure 10 shows denoised images produced by hybrid PI-PD model compared with the output of the BM3D and the average 5x5 filters. The figure illustrates the merits and drawbacks of each model: edges are well preserved by hybrid PI-PD, but a staircase effect is visible, a wellknown artefact inherent to this type of neighborhood filters. Our recent GPU-implementation of the regression method proposed in Buades et al (2006) brings a mean improvement of 1dB at the cost of 0.4 ms.

9 Conclusion, future work

From the start, our approach, unlike quite a few others, has been to base this study on the conception and characteristics of the targeted hardware (Nvidia Graphic cards).

So as to get high execution speeds, we chose, for example, to find a method that remains local (concentrating on the immediate neighborhood of the center pixel),

(a) Noisy image σ = 25 (b) Average 5x5 filter, in 0.35 ms

(c) PI-PD hybrid filter, n = 25, l = 5, T_{max} = 1, T 2_{max} = 2, in 11 ms

(d) BM3D filter, in 4.3s

Fig. 10 Comparison of 512x512 images denoised from noisy airplane image (10a) with a PI-PD filter (10b), PI-PD hybrid filter (10c) and BM3D filter (10d). Only zoomed parts of images are shown in order to ensure better viewing.

but still provides very significant benefits, using our technique of progressive lengthening.

Nevertheless, our method has proved slightly suboptimal and lacking robustness in flat regions (see above, Low Slope Regions), even if the actual visual effect may be considered quite satisfactory.

As a first step to address the above drawbacks, we have devised a hybrid method that detects and applies distinct processing to LSR regions (see above). Processing speeds remain fast, and much higher than the BM3D implementation taken as quality reference. This is very promising, and opens the perspective of real-time high definition image sequence processing at 25 fps, provided we improve the edge detector, which currently limits the HD frame rate at 16fps (High Definition: 1920x1080 pixels).

To further improve the quality of output images, we also implemented a efficient parallel implementation of the staircase effect reduction technique presented in Buades et al (2006). With this method, searching for best improvement factors leads to different parameters values for each image processed, which prompts to studying some way of overriding such parameters.

Our study so far has been based on additive noise; we are currently working on transposing criteria to various multiplicative noise types. We also extended the process to color images with very interesting visual results to be confirmed by the experimental measurement currently in progress.

References

(2007) DenoiseLab Philosophy: A Standard Test Set and

Evaluation Method to Compare Denoising Algorithms (2010) NVIDIA CUDA C Programming Guide v3.1.1.

NVIDIA Corporation Bertaux N, Frauel Y, Réfrégier P, Javidi B (2004) Speckle

removal using a maximum-likelihood technique with isoline gray-level regularization. J Opt Soc Am A 21(12):2283–2291, DOI 10.1364/JOSAA.21.002283, URL http://josaa.osa.org/abstract.cfm?URI=josaa-21-12-2283

Buades A, Coll B, Morel J (2005) A review of image denoising algorithms, with a new one. Multiscale Modeling and Simulation 4(2):490–530

Buades A, Coll B, Morel JM (2006) The staircasing effect in neighborhood filters and its solution. IEEE Transactions on Image Processing 15(6):1499–1505

Caselles V, Coll B, Morel JM (1997) Scale space versus to-

pographic map for natural images. Springer, pp 29–49 Chen W, Beister M, Kyriakou Y, Kachelries M (2009) High

performance median filtering using commodity graphics hardware. In: Nuclear Science Symposium Conference Record (NSS/MIC), 2009 IEEE, pp 4142–4147, DOI 10.1109/NSSMIC.2009.5402323

Coll B, Morel JM, Buades A (2011) Non-local means denois-

ing. Image Processing On Line Dabov K, Foi R, Katkovnik V, Egiazarian K (2009) Bm3d im-

age denoising with shape-adaptive principal component analysis. In: Proc. Workshop on Signal Processing with Adaptive Sparse Structured Representations (SPARS’09

Matheron G (1975) Random sets and integral geometry. Wi-

ley

McGuire  M  (2008)  A  fast,  small-radius
gpu  median  filter.  In:  ShaderX6,  URL

http://graphics.cs.williams.edu/papers/MedianShaderX6 Monasse P, Guichard F (1999) Scale-space from a level lines

tree. In: Nielsen M, Johansen P, Olsen O, Weickert J (eds) Scale-Space Theories in Computer Vision, Lecture Notes in Computer Science, vol 1682, Springer Berlin / Heidelberg, pp 175–186, URL http://dx.doi.org/10.1007/3-540-48236-9 16, 10.1007/3-540-48236-9 16

Palhano Xavier De Fontes F, Andrade Barroso G, Coupé P, Hellier P (2010) Real time ultrasound image denoising. Journal of Real-Time Image Processing DOI 10.1007/s11554-010-0158-5, URL http://hal.inria.fr/inria-00476122

Perrot G, Domas S, Couturier R, Bertaux N (2011) Gpu implementation of a region based algorithm for large images segmentation. In: Computer and Information Technology (CIT), 2011 IEEE 11th International Conference on, pp 291 –298, DOI 10.1109/CIT.2011.60

Sanchez RM, Rodriguez PA (2012) Bidimensional median filter for parallel computing architectures. In: Acoustics, Speech and Signal Processing (ICASSP), 2012 IEEE International Conference on, pp 1549–1552, DOI 10.1109/ICASSP.2012.6288187

Wang Z, Bovik AC, Sheikh HR, Member S, Simoncelli EP, Member S (2004) Image quality assessment: From error visibility to structural similarity. IEEE Transactions on Image Processing 13:600–612

Yang Q, Tan KH, Ahuja N (2009) Real-time o(1) bilateral filtering. In: CVPR, IEEE, pp 557–564, URL http://doi.ieeecomputersociety.org/10.1109/CVPRW.2009.5206542