Mobile-end Tone Mapping based on Integral Image and Integral Histogram

The integral image is prominently used within the Viola-Jones object detection framework from 2001 [19]. It can significantly reduce the computation burden when an accumulative sum of image area is required. The following equation defines an integral image.

where $i(x,y)$ is the value of the pixel at $(x,y)$. A great feature of integral image $I$ is that summation of any rectangular region in the original image $i$ can be computed efficiently in a single pass. For example, if there are four points $A(x_{0},y_{0})$, $B(x_{1},y_{0})$, $C(x_{1},y_{1})$ and $D(x_{0},y_{1})$ $(x_{0}<x_{1}$, $y_{0}<y_{1})$ in image $i$, the accumulative summation of rectangular $w$ that is enclosed by the four points is equal to:

It can be seen from Eq. 6 that with the help of integral image, only one addition and two subtractions are required for calculating an accumulative sum.

A natural extension of integral image is the integral histogram. It is a fast way to extract histogram of a specific area [20]. If the pixel values of an image $i$ are equally segmented into $n$ different bins of a histogram, then we can build an $n$-channel map $h$ where the $k$-th channel $h_{k}$ is defined as follow:

$b_{k}$ indicates the $k$-th bin. For each channel $h_{k}$, we can have an integral image $H_{k}$

Then, we can build a $n$-channel integral histogram $H$ where the $k$-th bin of the histogram is calculated with

It is evident that with the help of integral histogram. An $n$-bin histogram of any local area can be calculated easily with $n$ addition and 2$n$ subtractions.

With the help of integral histogram, the histogram calculation of Eq. 1 becomes straightforward. We first equally divide image $l$ into $n$ bins based on pixel values. $n$ is a user defined parameter. The integral histogram of the logarithmic image $l$ can be built according to Eq. 7 and Eq. 8. Then for any receptive field $w_{i}$, the pixel count of $f_{w_{i}}(b_{i})$ can be solved by Eq. 9. Hence, Eq. 1 and Eq. 2 can be derived easily. The local variance value $\sigma_{{w_{i}},(x,y)}^{2}$ has the following form:

where $\sum_{w_{i},(x,y)}i^{2}$ and $\sum_{w_{i},(x,y)}i$ are the accumulative sums of receptive field $w_{i},(x,y)$ centered at pixel location $(x,y)$, and $|w_{i,(x,y)}|^{2}$ is the number of pixels within the receptive field. The two accumulative sums in Eq. 10 can be computed with integral image of $i$ and $i^{2}$.

The majority of the computational burden of Eq. 1 to Eq. 4 are optimized with the help of integral image and integral histogram. The remaining computation is just some basic operations. After the formation of the integral image and integral histogram, the entire tone mapping process becomes pixel-parallel. Which means all pixels can be computed in parallel without data dependency. This feature is significantly helpful when it comes to GPU implementation where large amounts of single-instruction-multiple-data (SIMD) processing units are available. We will show the GPU pixel-parallel implementation in Section III.

In this part, we make a detailed analysis of how the proposed algorithm compresses the dynamic range while preserving local details.

The main reason that the algorithm adopts a simple tone mapping method (Eq. 2) rather than a more complicated one to tone map every receptive field $w_{i}$ is because the local dynamic range of a single receptive field $w_{i}$ is mostly much lower than the dynamic range of the entire WDR image. Fig. 1 shows an example that a simple tone mapping method is already sufficient to produce good results. The image on top of Fig. 1 is the original WDR image. We tone map this image with different algorithms. Two patches of the same size and at the same location are selected from the corresponding tone mapped images. From left to right on the two bottom rows of Fig. 1, the image patches are from Mantiuk et al. [8], Durand et al. [10], Reinhard et al. [21] and Drago et al. [3], respectively. The last column shows image patches from the original WDR image and then tone mapped with a logarithmic response function. The image patches using local processing show comparable or even better results when compared with other algorithms. In the lamp area, the image is not as saturated as the other four images, and in the car area, the image is brighter than the other four images. This example demonstrates that with local adaptation, even a simple tone mapping method can reveal the local details of wide dynamic range scenes.

However, revealing local detail is not enough for good tone mapping, because artifacts usually emerge if there is no consideration of global statistics. Figure. 2 shows a tone mapping example of our algorithm with only one receptive field ($s=1$ in Eq. 3). It is apparent that the smaller the receptive field, the better the details, but there are many disturbing artifacts especially at uniform areas. On the other hand, the larger the receptive field, the fewer the details as well as the artifacts. Actually, the size of the receptive field can balance between the global and local statistics during tone mapping. Consider the most extreme situation when the receptive field is as large as the WDR image itself, then the tone mapping becomes global. In this circumstance, there will be no artifacts because the tone mapping is global, it can preserve the brightness consistency of the image based on the monotonic tone mapping curve. But, the details are also concealed due to the global compression feature. If the receptive field is so small that there are only several pixels, the details will be greatly exaggerated because every receptive field will be given the maximum display range between $min(L_{d})$ and $max(L_{d})$ regardless of its size. If the pixels have similar values, for instance, they are all background pixels of clear sky, then, some pixel will still be tone mapped to $min(L_{d})$ and become artifacts because of Eq. 2. The local compression could sabotage the image brightness consistency, especially at uniform areas because the tone mapping function of the entire WDR image is no longer monotonic. Although the results are significantly different when tone mapped with large and small receptive fields, it is also evident that the displayed visual appearances are complimentary which contain local detail and global brightness consistency. Since every pixel is included in $s$ different receptive fields, one solution to obtain these local and global appearances is to fuse all receptive fields using a weighted summation as shown in Eq. 3.

To reveal local detail and maintain brightness consistency, we should consider the two possible cases of local receptive fields. In a “flat” receptive field where pixel values are mostly the same, we should tone map it more globally to maintain its image consistency and reduce possible artifacts. In a “high variance” receptive field where pixel values fluctuate significantly, we should tone map it more locally to reveal details. Eq. 4 is taken from guided image filter [12] and it is a very efficient way to measure the “flatten” or “variance” degree of a local image area. If the variance value $\sigma_{{w_{i}},(x,y)}^{2}$ is much larger than $\epsilon$, the weight value will be close to 1, which indicates that the centre pixel of receptive field ${w_{i}},(x,y)$ is “high variance”. If the variance is much less than $\epsilon$, the weight value will be close to 0, and the center pixel is regarded as belonging to a “flat” receptive field. Eq. 4 gives pixel-wise and receptive-field-wise weight for each pixel and it has edge preserving feature that is helpful for artifacts reduction. Fig. 3 shows the weight map calculated with different receptive field. Bright pixels mean the computed weight value is close to one and dark pixels mean that it is close to zero. It can be seen that in smaller scales, the weight can better preserve details such as edges and textures.

There are three free parameters involved in our algorithm, namely, the number of different receptive fields $s$, the number of the histogram bins $n$, and the regularization term $\epsilon$. The number of receptive fields $s$ is the most important parameter in the proposed algorithm. As we have discussed previously, to maintain the image brightness consistency, the largest receptive field should be the same as the image size. Yet to obtain possible details in different receptive fields, we need smaller receptive fields as well. We adopt the popular image pyramid method and define the relationship for any two adjacent receptive fields $w_{i+1}$ and $w_{i}$ as $w_{i}/w_{i+1}=2$. Consider a $2048\times 2048$ WDR image, $s=5$ will yield the smallest receptive field of $128\times 128$ which is enough for the tone mapping function of Eq. 2 to preserve local details. The effect of the number of bins $n$ and the regularization term $\epsilon$ parameters for a tone mapped image are shown in Fig. 4. Nine tone mapped images are presented in a matrix with $n$ varying vertically and $\epsilon$ varying horizontally. It is apparent that the image contrast increases as $n$ increases. The decrease of $\epsilon$ in Eq. 4 will result image with better local details. In Fig. 4, the texture of the window glass is best preserved with smallest $\epsilon$ value. We find values for $n>5$ and $\epsilon=0.1$ to usually produce satisfactory results, i.e. good brightness while preserving local contrast and details.

The algorithm is implemented using Apple’s GPU application programming interface, Metal. Metal allows for the parallel processing of data much like OpenGL or other shader programming languages. In Metal, kernel function is the basic function that runs in SIMD fashion in GPU. Each instance of kernel function is called a thread. Threads are further organized into threadgroups that are executed together and can share a common block of memory. Since the proposed algorithm is pixel-parallel, then we can assign every pixel a thread to tone map itself and all threads share the common memory which is the calculated integral image and integral histogram. The overall implementation flow of the algorithm is shown in Fig. 5. The input and output of the algorithm are the calculated lumniance image and tone mapped lumniance image, respectively. The blue blocks indicate functions that execute only once during the computation. They are used to calculate the logarithmic image, integral histogram $H$ and integral image $I(i)$ and $I(i^{2})$, respectively. MPSImageIntegral function that is provided by Metal Performance Shaders is used to calculate the integral image and integral histogram. The core functions are labeled with green color in Fig. 5.

Our implementation calculates different scales in sequence, the functions that are required in this procedure are labeled as green, and they will be called for $s$ times in a main loop (shaded gray). The pseudocode code of Algorithm 1 shows how the Eq. 4 is implemented in GPU as kernel function. $I_{\sum}$ and $I^{2}_{\sum}$ are the computed integral image and of $i$ and $i^{2}$, respectively. They can be accessed by all threads. $grid.x$ and $grid.y$ are the vertical and horizontal coordinates of current thread. Line 2 to line 5 are used to calculate coordinates of the four corner pixels of Eq. 6 and Eq. 9. $|w|$ calculates the number of pixels in the receptive fields. Eq. 6 and Eq. 9 are executed in line 13 and line 18. Line 19 implement Eq. 10 and Eq. 4 is carried out in line 20. The implementation of algorithm 1 only employs basic addition, subtraction and indexing operations of GPU.

Eq. 1, Eq. 2 and Eq. 3 are also implemented as kernel functions because they only depend on the pixel value and the integral histogram. The parallel processing feature of the algorithm greatly reduce the computation complexity and cost of the implementation. Evaluation results of the GPU implementation will be shown in the following section.

It can be understood from the earlier explanation that our algorithm utilizes all available display levels in each local processing. This significantly increases image brightness, especially in dark regions. Our experiments tested different WDR images to evaluate the image quality of the tone mapped images. One example is shown in Fig. 6. It shows one image tone mapped with different algorithms. In reading order, the images are tone mapped with Durand et al., Gu et al., Paris et al., and the proposed algorithm, respectively. A visual impression can tell that our algorithm gives an image with higher brightness value over the other three images. In fact the average brightness value for the four images is 93.78, 120.92, 98.01 and 122.14 respectively. Despite the fact that the image generated with Gu et al. has similar overall brightness, it lacks global contrast and many unwanted details are enhanced which makes the image look unnatural.

For objective assessment, we use the tone-mapped image quality index (TMQI) [24] to calculate an overall quality score that combines a multi-scale structural fidelity measure and a measure of image naturalness. The structural fidelity measure is a full-reference assessment based on the structural similarity (SSIM) index. The naturalness measure is a no-reference assessment based on statistics of good-quality natural images. The results of the naturalness score, structural fidelity score, and the overall score are listed in Table I, Table II and Table III. The winner algorithm’s score is shown in bold font. In the naturalness score, our algorithm scores highest in 5 images and achieves an average value of 0.8221 for which is highest among the tested algorithms. In terms of structural similarity, our algorithm wins 9 out of 10 images and achieves an average score of 0.9213. In terms of overall quality, our algorithm produces the best scores for 9 images.

Field assessment of the image quality has also been conducted using the developed iOS application. Fig. 6 shows two sets of examples. The first scene is taken in a commonly seen extreme light condition. The four images with different exposure are shown on the two left columns and the tone mapped image is shown on the right-most column. Our algorithm produces a bright image with clear details — For example, the outdoor cloud and the frame of skylight window are clearly visible. The inner structure of the building such as stairs and railing are also clearly shown in the tone mapped image. Another scene is a portrait image because taking portraits is the most common use of phone camera. The bottom image of Fig. 6 shows an example image of a portrait taken by our application. This image is taken under backlight condition, it is a challenging situation because the strong background light could make the portrait dark or even invisible. The tone mapped image in Fig. 6 shows the portrait clearly and the detail in the background is also well preserved. This experiment indicates that our application is very suitable for imaging under extreme light conditions.

Ideally, we hope to compare the application performance with some other works which are implemented also in mobile-end. However, most tone mapping algorithms reported are implemented on desktop platforms. It would be unfair to compare implementations on different platforms. Hence, we focus our analysis on the mobile-end implementation.

The proposed algorithm has three parameters, namely the number of scales $s$, the number of bins $n$ and the regularization term $\epsilon$. A detailed performance analysis should fully consider the variation of the three parameters. The regularization term $\epsilon$ is a scalar value which does not affect neither the memory requirement nor the total amount of computation, so it will not affect the performance of our application. However, the number of scales affects the computation burden and the number of bins $n$ affects the required memory usage. To evaluate the two parameters $s$ and $n$ ’s effect on the processing speed, we tested $9$ parameter sets, varying $s$ from 3 to 5 and $n$ from 3 to 5. The evaluation results of three commonly used resolutions $480\times 640$, $720\times 1280$ and $1920\times 1080$ are shown in Table IV. For identical parameter settings, the processing time increases mostly linearly as the image resolution increases. For example, the processing times are about $200$ ms, $400$ ms and $1000$ ms for the three different resolutions when $n=5,s=5$. The influence of the $n$ parameter on the processing time is limited. Under same image resolution with a fixed number of scales, the processing time only fluctuates in a very limited range. In $720\times 1280$ resolution, the processing time are $402$ ms, $422$ ms and $432$ ms when $s$ equals to $5$ and $n$ equals to 3, 4 and 5, respectively. In $1080\times 1920$ resolution, the processing time is $788$ ms, $793$ ms and $886$ ms when $s$ equals to $4$ and $n$ equals to $3$, $4$ and $5$, respectively. From Table IV, we can see that the scale parameter $s$ has greater influence on the application processing time than the $n$ parameter. In $480\times 640$ resolution, the application processing time is $122$, $145$ and $166$ ms when the $n$ is equal to 3 an $s$ is equal to 3, 4, 5, respectively. As for the $720\times 1280$ and $1080\times 1920$ resolution, the processing time can increase about $50\%$ when scale $s$ changes from 3 to 5. It is easy to conclude that the greatest factor that affects the application processing time is the image resolution. Parameter $s$ affects the application time secondarily while parameter $n$ has only very little influence.

A processing time breakdown using two different parameter settings is shown in Fig. 7. Fig. 7 (a) is the processing time breakdown when $n=5$, $s=5$ and image resolution is $1920\times 1080$. Fig. 7 (b) is the processing time breakdown when $n=3$, $s=3$ and image resolution equals to $480\times 640$. We choose the two extreme cases to show the rough proportion of each function. Since Eq. 3, Eq. 4 and Eq. 1 & Eq. 2 will be looped for several times, hence, they consumes the majority of processing time. They account for 23%, 28% and 24% of the total processing time in Fig. 7 (a) and 32%, 28% and 19% in Fig. 7(b). The integral histogram and integral images only need to be computed once in our algorithm and they consume 19% and 6% of processing time in Fig. 7 (a), 18% 3% in Fig. 7(b) respectively.