1 \documentclass[a4paper,10pt]{article}
2 \usepackage[utf8]{inputenc}
5 \title{Fast GPU-based denoising filter using isoline levels}
13 In this study, we describe a GPU-based filter for image denoising whose principle rests on Matheron's level sets theory, first introduced in 1975 but rarely implemented in digital image processing because of high computation costs. We use the fact that, within a natural image, the significant contours of objects coincide with parts of the image level-lines.
14 The presented algorithm assumes an \emph{a priori} knowledge of the corrupting noise type and uses the polygonal level-line modelling constraint to estimate each pixel's gray-level by local maximum likelihood optimization.
15 Over the 11 tested images, results show a high quality/runtime ratio with an average improvement of 7.14 dB of the SNR ratio and 30\%\ of the structural similarity index, performed in no more than 0.011, which is 350 faster than the BM3D filter, and potentially allows processing high definition images at over 16 fps.
19 \section{introduction}
20 The fast growth in digital devices able to take pictures
21 or make movies has made digital processing become more
22 important. However, the wide range of applications requiring noise
23 removal makes it difficult to propose a universal filter; concurrently, the increase in pixel density of the CCD
24 or CMOS sensors used to measure light intensity leads
25 to higher noise levels and requires higher data rates in
26 processing algorithms.
27 To quantify the quality of an image processing algorithm one should also take into account the fact that human perception varies significantly from person to person.
29 To date, many researchers have successfully speeded up GPU image processing
30 algorithms. For example Mc Guire [11], Chen
31 et al [7] and Sanchez [15] reported quite fast median filters, able to output up to 300 million pixels per second, i.e. 150 fps of high definition video
32 images. Bilateral filltering has also been successfully proposed by Yang [17].
34 Still, most high quality algorithms,
35 like NL-means [8] or BM3D [9] make use of non-local similarities and/or frequency domain transforms and the speedups achieved by their current GPU implementations, though quite significant (as shown for example with
36 NL-means in [13]), do not come near those achieved by local methods such as gaussian, median or neighborhood filters as they have not \emph{originally} been designed for the specific GPU architecture.
39 It is established that a natural image, i.e any photograph of an outdoor or indoor scene taken by a standard camera, is of bounded variation and thus can be
40 decomposed into a set of level-lines [?]. Accordingly, researchers have succeeded in implementing algorithms (mostly dedicated to
41 to segmentation and classification) [6, 12? ]. In that respect, \emph{Bertaux
42 et al} described a method which significantly reduces speckle noise inside coherent images by using level lines to constrain the minimization process [3]. Those level lines, more precisely \emph{iso-gray-level lines}, are called
43 \emph{isolines} and, as explained in [3], consist in neighborhoods of polyline shapes determined by use of the maximum likelihood optimization.
44 This method has proved to bring good enhancement and to preserve edges between regions, but its costs in computation time (a little less than one minute on a PIII-1Ghz cpu processing a 2-Mpixel image) does not allow real-time image processing.
46 We started by designing a set of GPU implementations with various optimization heuristics, in order to find out how to both minimize loss in quality and preserve admissible execution times. Algorithms have been tested with reference images taken from the site of \emph{DenoiseLab} [1] and for which results have been published in [? ?]. As will be detailed further down, our statistical observation of the output images produced by the \emph{Bertaux et al} method enabled us to propose a very fast and simple parallel GPU-based denoising method with good results in terms of average gray-level error and edge preservation. In terms of processing speed, on the basis of the BM3D timings listed in [9] and with our own measurements, our GPU-based filter runs around 350 times faster and thus is able to process high definition images at over 16fps.
48 In the following, section 2 briefly focuses on the recent
49 characteristics of recent Nvidia GPU's. Section 3 introduces the
50 theory and notations used to define isolines. Section 4 describes the two isoline-based models that lead to the final hybrid model, while section 5 details the parallel implementation of the proposed algorithm. Finally, we present our results in section 6 before drawing our conclusions and outlining our future work in section 7.
52 \section{NVidia’s GPU architecture}
54 GPUs are multi-core, multi-threaded processors, optimized for highly parallel computation. Their design focuses on the SIMT (Single Instruction Multiple Threads) model that devotes more transistors to data processing, rather than resorting to data-caching and flow control. We only have had access to two Nvidia cards, respectively models C1060 and the more recent C2070, designed for general purpose computing. An Nvidia C2070, for example, features 6 GBytes global \textit{off-chip} memory and a total of \textit{448 cores} (1.15 GHz) bundled in several Streaming Multi-processors (SM), along with an amount of much faster \textit{on-chip} memory : \textit{thread registers} (63 per thread)
55 and \textit{shared memory} (up to 48 KB per thread block). The memory bandwidth claimed by the manufacturer is 144 GB per second.
57 Moreover, Nvidia provides CUDA (for Compute Unified Device Architecture) which makes it quite easy to write parallel code
58 for their target GPUs. However, writing efficient code is not trivial as serialization must be avoided as much
59 as possible. Thus, code design requires that one pays attention to a number of points, among which:
62 CUDA organizes threads by a) thread blocks in which
63 synchronization is possible, b) a grid of blocks with
64 no possible synchronization between them.
66 the order in which threads are scheduled is not defined. Threads are run in parallel by groups called\textit{ warps}. For each, distinct executions engines run odd-indexed threads and even-indexed threads.
69 data must be kept in GPU memory, to reduce the overhead generated by copying between CPU and
72 the total amount of threads running the same computation must be as large as possible.
75 the number of execution branches inside one block should be as small as possible.
77 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 an 128 Byte aligned block.
79 shared memory is organized in 32x32 bit-wide banks.
80 Bank-conflict occurs when two threads in a half-warp
81 map to the same bank but at different locations.
86 In the following, let \textit{I} be the reference noiseless image
87 (assuming we have one), \textit{I'} the noisy acquired image corrupted by one additive white gaussian noise of zero mean
88 value and standard deviation σ. Let \textit{Î} be the denoised image. Each pixel of \textit{I} of coordinates (i, j) has its own gray level z(i, j).
89 Within the noisy image, our goal is to select, for each pixel, the set of neighbor pixels that locally best fits the underlying level line of the reference image (in Matheron's meaning). In the rest of this article, we call \textit{isoline} one such set of pixels. An \textit{isoline} can be one single \textit{isoline-segment} or one \textit{polyline} composed by at least two connected \textit{isoline-segments}. The generalized likelihood criterion (GL) is used to select the best \textit{isoline(-segment)} among all the considered ones.
91 \subsection{Isoline-segments}
93 \subsection{Extendable isolines}
96 \section{isoline models}
98 \subsection{Poly-isolines with limited deviation angle (PI-LD)}
100 \subsection{Poly-isolines with computed directions (PI-PD)}
102 Though much faster, the PI-LD-based filter may be considered inferior to state-of-the-art filters such as BM3D family algorithms [9]. In addition, it has to be noted that this way of building poly-isolines requires the alternate use of two different types of validation at each extension stage, namely GLRT and maximum likelihood minimization, each of them generating numerous divergent branches during kernel execution that have to be separately executed by each thread warp, leading to unwanted overhead.
104 For each pixel (i, j), as no selection is done at the first stage within the PI-LD model, \textit{D} poly-isolines are computed and kept as candidates but, since only one follows the actual isoline at (i, j), it becomes unnecessary to perform the selection at each lenghtening step to achieve a robust determination of the direction at any given pixel of this isoline and only the first segment has to be determined to obtain the local direction of the isoline, which leads to an important reduction of the work complexity :
105 the above PI-LD model needs to evaluate (2.∆dmax + 1) segments at each extension stage (∆dmax segments on both side of the previous orientation angle). As we allow a maximum of K extension stages and as the first D orientation are mandatorily evaluated, it brings K−1
106 to the count of D.(2.∆dmax + 1) segments to evaluate at each pixel position, while only D.K evaluations are needed in the second case.
107 For example, with a maximum of K = 5 segments and a maximum orientation change of ∆dmax = 2, the PI-LD needs to evaluate up
108 to 20000 segments per pixel where only 160 would be enough.
113 The PI-PD model brings significant simplification of the computing process, as each pixel pattern is only needs to be applied once along with the computation of its associated values, that are re-used every time one poly-isoline’s segment ends (at pixel (i, j). In addition, multiple branches are avoided during kernel execution, which better fits the GPU's architecture.constraints.
114 It remains that the building of poly-
115 isolines is done without global likelihood optimization.
116 In addition to the extension process, the model has been further improved by adding the ability to thicken poly-isolines from one to three pixels, allowing to achieve higher PSNR values by increasing the number of pixels of poly-isolines: This may apply preferably to
117 large images that do not contain small relevant details liable to be slightly blurred. Still, this feature makes PI-PD
118 more versatile than our reference BM3D which has prohibitive computation times when processing large images (over 5 minutes for a 4096x4096 pixel image).
121 \subsection{Hybrid PI-PD}
123 As the determination of each segment’s direction only involves 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\textit{ Low Slope Regions} (LSR). Figure 7 shows that lack of robustness with two successive drawings of additive white gaussian noise applied to our reference image of figure 6. We focused on a small 11x11 pixel window containing two LSR with one sharp edge between them. Figures 7d and 7e 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 LSRs (lines 1-4, 8-11).
127 It must be noticed that point b) has been introduced
128 in order to achieve smoother transitions between regions
129 to which PI-PD is applied and those in which the plain
130 average value is used. Figure 9 shows an example of such
131 classification achieved by the edge detector. The detector has been applied to the noisy airplane image
132 with a GLRT threshold value T 2max = 2. Darker dots represent pixels classified as belonging to an edge, whiter dots pixels classified as belonging to a LSR.
136 \section{conclusion - work in progress}
138 Instead of only searching to obtain speedups from existing CPU algorithms, we have tried to design specific GPU-based high-speed and robust elementary processing cells, aimed at achieving both satisfactory noise reduction and high execution speeds. For that purpose, we chose to design a method that remains local (concentrating on each center pixel) and provides very significant benefits thanks to our technique of progressive \textit{isoline} extension. Nevertheless, despite the fact the actual visual effect may be considered quite satisfactory, it proved slightly sub-optimal and lacking robustness in flat regions (see above, Low Slope Regions), which led us to devise 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
139 HD frame rate to 16fps-high Definition (1920x1080 pixels).
142 To further improve the quality of output images, we also implemented a parallel version of the staircase effect reduction technique presented in [5]. With this method, searching for best improvement factors leads to different parameters values for each image processed, which
143 prompts to studying some way of overriding such parameters.
144 So far, our research has been focussed on additive noise; we are currently working on transposing criteria to various multiplicative noise types. We have also started extending the process to color images with very interesting visual results to be confirmed by the experimental measurement currently in progress.
