1.

Introduction

The dynamic range of natural scenes varies accordingly with different lighting conditions, from lower than 40 dB for ordinary uniformly lit scenes to higher than 120 dB for mixed indoor-outdoor scenes. However, human observers can cope well with such tremendous difference through local adaptation embedded in retina and cortex. Powered by sophisticated mechanisms, the human visual systems (HVS) can perceive a very large intensity range simultaneously, and even larger after long-time adaptation.^1,2 On the other hand, the dynamic range of current image displays and other reproduction media is rather limited. Consequently, there is a great demand for high-fidelity dynamic range reduction in many media-related applications, such as remote imaging, medical imaging, virtual reality, and digital photography.^3–5

Algorithms that map high-dynamic-range (HDR) input intensities to low-dynamic-range (LDR) displayable signals are called tone mapping operators (TMOs). Tone mapping has become a hot topic during the past two decades, and abundant publications about this topic can be found. Comprehensive reviews of the state of the art in this field are given in Refs. 1, 4, and 6. Among the tone-mapping literature, some TMOs apply a certain method to decompose the original image into several layers with different scales, and produce results by manipulating the dynamic range of these layers independently.^3,7–10

In this paper, similar to decomposition-based TMOs, we propose an illumination-reflectance decomposition method based on the iterative Retinex algorithm for HDR tone mapping. The iterative Retinex algorithm¹¹ is combined with the anisotropic diffusion¹² in computing the illumination. The discontinuities in the illumination can be well described using the proposed estimator. Moreover, a jumping-spiral iteration is proposed to improve the symmetry of the edge response over aligned directions. The proposed illumination estimation is later adopted in HDR tone mapping. Results show that the proposed algorithm is able to produce good results. Meanwhile, undesirable artifacts are effectively suppressed, which are usually not avoided using the standard Retinex algorithms. The proposed algorithm also outperforms some similar TMOs, evaluated by the state of the art objective metric.¹³

The rest of this article is arranged as follows. In Sec. 2, the related work is briefly reviewed. Then the computational background of the proposed method is presented in Sec. 3. In Sec. 4, the proposed illumination-reflectance decomposition method is derived in detail. It is then applied in HDR tone mapping later in this section. Experiments are done in Sec. 5, where the effectiveness of the proposed algorithm in tone mapping is approved, and encouraging results are produced compared with some other similar algorithms. Finally, conclusions are drawn in Sec. 6.

2.

Related Works

3.

Computational Background

4.

Adapting the Iterative Retinex Computation for Tone Mapping

The basic assumption applied by Retinex theory is that the illumination is spatially smooth. In most cases, however, the abrupt edges in HDR scenes are caused by the illumination discontinuities. As a consequence, in reproducing these scenes, the results will inevitably disrupted by halos using the standard Retinex algorithms.

In order to obtain better results in tone mapping, we make several modifications based on the McCann–Sobel Retinex algorithm. First, in order to better compensate the asymmetric effect of the spiral iteration without increasing the computational complexity, we slightly modify the order of the computation process. Second, an edge-stopping function is introduced into the iterative smoothing computation. Such modification enables the representation of the abrupt edges of the illumination. Third, the estimated reflectance is not directly used as the tone-mapped result. Instead, it is relit by the range-compressed version of the estimated illumination to produce visually more pleasing results. Moreover, the proposed algorithm is only performed in the luminance component of the original image rather than individually applied in three color channels as done in the standard Retinex algorithms. The color information is reproduced according to the original image using an existing technique.

4.1.

Jumping-Spiral Iteration

The original iterative Retinex algorithm proceeds in a spiral manner and produces asymmetric spatial response. This can easily cause distortions. Although this kind of asymmetric response can be reduced by increasing the number of iterations, the computational complexity will increase, and the strength of the detail enhancement will weaken as well.

The dual-spiral iteration manner proposed in Ref. 37 is axial-symmetric. However, the edge responses in horizontal directions are different from that in vertical directions. Furthermore, the dual-spiral iteration manner is less effective than the single-spiral one with comparative computational burden using the McCann–Sobel algorithm.

We propose a jumping-spiral iteration manner as illustrated in Fig. 1. After every loop where neighboring pixels with the same shift distance are processed once, the process jumps to the opposite side of the starting direction of the previous loop. The dashed arrows denote the jumps.

Fig. 1

The proposed jumping-spiral iteration manner: This figure depicts the case that only one iteration is taken.

In order to test the symmetry of different iteration manners over directions, a computer-generated image is used in this experiment, as shown in Fig. 2(a). The test image is a $128\times 128$ monochrome image. The size of the bright box in the center of the image is $32\times 32$. Then, following the McCann–Sobel algorithm, the three different iteration manners are adopted in computing the illumination estimation of the test image. The results are respectively shown as Fig. 2(b)–2(d), which are all produced with two iterations.

Fig. 2

The spatial behavior of the iterative smoothing with different iteration manners; (a) is the input image, (b), (c) and, (d) are, respectively the results of the single-spiral iteration manner, and the proposed jumping-spiral iteration manner.

To give an objective evaluation of the nonsymmetry of the iteration manners, we define an asymmetric error measure of a square matrix as follows:

Eq. (8)

$$m=\frac{E(I)+E(I_f)}{2},$$

where $I$ and $I_f$ are the square matrix to be tested and its flipped version, respectively. $E(·)$ is an error norm measuring the diagonal symmetry of a square matrix, which is defined as

Eq. (9)

$$E(I)={(\int {|I(x,y)I(y,x)|}^2)}^{\frac{1}{2}}.$$

Since the test image is a diagonal-symmetric square matrix, the asymmetric error of the results over directions can be predicted by the corresponding value of $m$. The smaller $m$ indicates the better symmetry of the result.

The asymmetric errors of the three iteration manners decrease as the number of iteration increases, as plotted in Fig. 3. It is shown that the proposed jumping-spiral iteration manner produces less asymmetric error than the other two when the number of iterations $n\ge 2$. Consequently, the symmetry of the spatial response of the spiral iteration can be further improved using the jumping-spiral manner.

Fig. 3

Asymmetric errors versus number of iterations.

4.2.

Iterative Data-Dependent Smoothing

Without nonlinear operations, the iterative Retinex computation can be interpreted as an asymmetric isotropic diffusion, as depicted in Eq. (7). The diffusion procedure is data independent, and thus smooths over edges. Inspired by the anisotropic diffusion, an iterative data-dependent smoothing method is proposed.

The discrete form of the anisotropic diffusion proposed by Perona and Malik¹² can be formulated as follows:

Eq. (10)

$$I_p^{k+1}-I_p^k=\frac{\lambda }{|{\eta }_p|}\sum _{{q\in \eta }_p}g(|I_q^k-I_p^k|)(I_q^k-I_p^k),$$

where $\lambda$ is the step size of the discrete time steps. ${\eta }_p$ denotes the spatial neighborhood of the pixel at $p$, and $|{\eta }_p|$ is the number of neighbors. $g(·)$ is the edge-stopping function, which is generally a Gaussian function.

The anisotropic diffusion is very effective in edge-preserving local smoothing. Comparing Eq. (10) with Eq. (7), the corresponding asymmetric anisotropic diffusion is formulated as follows:

Eq. (11)

$$l_p^{k+1}-l_p^k=\frac{1}{2}g(|l_q^k-l_p^k|)(l_q^k-l_p^k).$$

The above formulation is a nonlinear data-dependent smoothing. In order to reduce the blocking and ringing effects, it is further modified to be a linear data-dependent smoothing method:

Eq. (12)

$$l_p^{k+1}-l_p^k=\frac{1}{2}g(|i_q-i_p|)(l_q^k-l_p^k).$$

The data-dependent smoothing formulation is then adopted in the iterative procedure. Given a natural monochrome image as shown in Fig. 4(a), the data-independent smoothing, as given by Eq. (7), produces an overall smoothed image, as shown in Fig. 4(b). On the contrary, the data-dependent smoothing, as described by Eq. (12), is able to preserve sharp edges, as shown in Fig. 4(c).

Fig. 4

Data-dependent smoothing versus data-independent smoothing: (a) is the original natural image, which is downloaded from http://image.baidu.com, (b) and (c) are respectively the spatial behavior of the data-dependent smoothing method depicted in Eq. (7) and the data-independent smoothing method depicted in Eq. (12).

4.3.

Edge-Preserving Illumination Estimation

Similar to Eq. (6), the asymmetric linear anisotropic diffusion described in Eq. (12) can be formulated as a weighted sum of the intermediate values at the neighborhood center and at one of its neighboring point as follows:

Eq. (13)

$$l_p^{k+1}=[1-\frac{g(|i_q-i_p|)}{2}]l_p^k+\frac{g(|i_q-i_p|)}{2}{l}_q^k.$$

In order to characterize the discontinuities in the illumination, the data-dependent smoothing is adopted in computing the upper envelope of the input image. Following the McCann–Sobel algorithm, the proposed edge-preserving estimator for the upper envelope is formulated as follows:

Eq. (14)

$$l_p^{k+1}=[1-\frac{g(|i_q-i_p|)}{2}]l_p^k+\frac{g(|i_q-i_p|)}{2}\max \{i_p,l_q^k\}.$$

It is experienced that local structures can be over-smoothed in the illumination estimation with a small number of iterations, and thus the corresponding reflectance is over-exaggerated. When increasing the number of iteration, however, the depth of the reflectance map is enlarged, and thus the effectiveness and the efficiency of the algorithm are reduced as well. In order to better handle these problems, the final illumination estimation is obtained by averaging the output with a small number of iterations and the original intensity:

Eq. (15)

$$l=\frac{i+\widehat{l}}{2},$$

where $\widehat{l}$ is the edge-preserving upper envelope estimation by Eq. (14).

Note that the result of Eq. (15) is still an upper envelope of the original image. Given the same natural image as shown in Fig. 4(a), the illumination estimation and the corresponding reflectance are produced by the proposed algorithm, which are respectively shown in Fig. 5(a) and 5(b). Notice that the illumination estimation inherits the edge-preserving property of the asymmetric anisotropic diffusion depicted in Eq. (12). Furthermore, the illumination discontinuities are well described in the estimated illumination. As a result, the corresponding reflectance image is a halo-free LDR image.

Fig. 5

The illumination-reflectance decomposition using the proposed edge-preserving illumination estimation; (a) is the estimated illumination and (b) is the corresponding reflectance.

5.

Experimental Results

In order to illustrate the effectiveness of the proposed algorithm in HDR tone mapping, we experiment on the method in three aspects. First, the validity of the proposed method is tested using a common HDR radiance map. Second, the effects of varying the parameters are illustrated and discussed. Third, the proposed algorithm is objectively evaluated and compared with some closely related TMOs.

First of all, given an HDR radiance map as shown in Fig. 6, It is first converted into a monochrome intensity map according to Eq. (16). Then it is decomposed into two portions using the proposed algorithm, an HDR illumination map and an LDR reflectance map, as shown in Fig. 7(a) and 7(b), respectively. Then the dynamic range of the illumination estimation is reduced using a Gamma mapping. Afterwards, the reflectance map is relit by the range-compressed version of the estimated illumination, thus leads to a new intensity map with relatively much lower dynamic range. Finally, the color information is reproduced according to the original radiance map as discussed in the previous section. The corresponding tone-mapped result using the proposed algorithm is shown as Fig. 7(c), which is produced with $S=0.8$, $\gamma =5$, $\sigma =1$, and $n=1$. Compared with the global Gamma-mapped result as shown in Fig. 7(d), all details are preserved while the overall dynamic range is significantly reduced without visible artifacts.

Fig. 6

Window scene: the original HDR radiance map is displayed in linear scale, from (a) to (d), normalized, $\times 10$, $\times 100$, and $\times 1000$, respectively. The details all around the scene cannot be simultaneously represented in any one of the left four linearly scaled images, (e) is the color-coded intensity map. The numbers on the right side of the color bar denote the log10 units of the intensity. It is shown that the dynamic range of the original radiance map is about 100 dB.

Fig. 7

Validation of the proposed tone mapping algorithm: (a) Illumination estimation, in which the abrupt variations of the original intensity are mainly encoded. (b) Corresponding reflectance map in linear space, which has a much lower dynamic range than the original map and the illumination estimation. (c) Tone-mapped result produced by the proposed algorithm, where fine details can be perceived without visible artifacts. (d) Globally mapped image using the same $\gamma$ as in producing (c). However, details are diminished.

Second, several parameters are used in the proposed algorithm, which can be tuned to obtain better tone-mapping results. The strength of the dynamic range reduction can be tuned by the parameter $\gamma$ in Eq. (17). The value of $\sigma$ indicates the threshold determining the magnitude of edges to be preserved in the illumination. And the number of iteration can also slightly affect the appearance of the output image. The results obtained using different parameter settings can be illustrated in Figs. 8Fig. 9–10. Although the saturation of the reproduced color can be tuned by varying $S$ in Eq. (18), it is not discussed in detail in this paper. Users can easily experience the performance of different color saturation and tune the value of $S$ as needed.

Fig. 8

Belgium house: results illustrating the effect of varying the strength of the illumination range reduction. From (a) to (d) $\gamma =2$, 4, 8, and 12, respectively, where $\sigma =1$, $n=1$ and $S=0.8$. Larger $\gamma$ leads to stronger compression of the dynamic range, while smaller $\gamma$ produces a result closer to the original image. HDR radiance map courtesy of Dani Lischinski.

Fig. 9

Atrium night: results illustrating the effect of varying the sigma of the edge-stopping function. From (a) to (d), where $\sigma =0.1$, 0.5, 2.5, and 12.5, respectively, where $\gamma =5$, $n=1$ and $S=0.8$. Halos are eliminated with smaller $\sigma$, but details may diminish in the tone-mapped result and the ringing artifacts may be elevated at the same time. On the contrary, larger $\sigma$ produces better overall contrast, but the halos may become visible as $\sigma$ increases. HDR radiance map courtesy of Karol Myszkowski.

Fig. 10

Desk scene: results illustrating the effect of varying the number of iteration. From (a) to (d), $n=1$, 4, 16, and 128, respectively, where $\gamma =5$, $\sigma =1$ and $S=0.8$. Smaller number of iteration leads to better photographic look of the scene. Larger number of iteration will increase the depth of highlight, and results in better fidelity regarding to the original scene. HDR radiance map courtesy of ILM.

Third, as stated in Ref. 1, a good TMO should reproduce realistically the color sensation and detail visibility of the original scene. Thus, in order to evaluate the performance of a TMO, one should quantitatively measure its ability in reproducing these two features. Color reproduction of the proposed algorithm is done simply following the existing works. The quality of the reproduced color sensation is not further discussed in this paper.

With respect to contrast reproduction, several subjective metrics can be found in the literature.^38–40 The authors ask subjects to tell the differences between the tone-mapped images and the corresponding HDR real-world scenes or the corresponding linear reproductions on an HDR display. However, sometimes it is not applicable to ask subjects to evaluate the computer-generated HDR images, or real-time HDR imaging pipelines. Meanwhile, the prototypes of the HDR display are now still too expensive for low-cost applications, and they are not yet commonly available. Fortunately, objective evaluation metrics are proposed. Cadik et al.⁴¹ propose a method based on a basic set of attributes, which is built on the statistics of subjective experiments. Aydin et al.¹³ proposed a quantitative measure of the quality of the tone-mapped images based on psychophysical experiments. In this remarkable work, the original HDR radiance map can be directly used as the reference, and the corresponding tone-mapped image obtained by some TMO is compared using a carefully calibrated perceptual model. The visible differences between the two images are categorized into three types. As noted by the authors, they are loss of visible contrast, amplification of invisible contrast, and contrast reversal, respectively. A distortion map of the same size as the original image is given as the output. The distortion type and the corresponding magnitude are respectively coded using certain color and its saturation.

To save the computational cost, we simply set $n=1$ and $S=0.8$. After experiencing the visual appearance of the results with different $\gamma$ and $\sigma$, we come to a parameter set as listed in Table 1. In the rest of this paper, the results of the proposed algorithm are produced using these parameters, unless otherwise noted. Afterward, the performance of the proposed algorithm is evaluated using the objective metric proposed in Ref. 13. Following the original work of this metric, we use green for loss of visible contrast, blue for amplification of invisible contrast, and red for contrast reversal, when illustrating the distortion maps of the tested images. These types of distortion errors, together with the average error, are denoted by $E_l$, $E_a$, $E_r$, and $E_{\mathrm{avg}}$, respectively, in the following texts. Some of the tone-mapping results by the proposed algorithm on commonly available HDR radiance maps with the corresponding distortion maps are illustrated in Fig. 11.

Table 1

Default parameter setting.

Fig. 11

Some results obtained by the proposed algorithm: Images in the first column are the tone-mapped results produced by the proposed algorithm, while the images in the second column are the distortion maps corresponding to the first column, respectively. As shown in the tone-mapped images, details can be observed without visible artifacts. Thus only a few visible differences are reported by the evaluation metric in the corresponding distortion maps. First row: Foggy night, HDR radiance map courtesy of Jack Tumblin. Second row: Cathedral, HDR radiance map courtesy of Dani Lischinski. Third row: Stanford memorial, HDR radiance map courtesy of Paul Debevec.

The iterative Retinex algorithms in Refs. 11, 19, 20, and 37 are not designed for HDR tone mapping, thus their performance is not discussed here. Nevertheless, the improvements in Ref. 31 are very effective in HDR rendering. This remarkable variant of the iterative Retinex algorithm is then involved in the comparison. Moreover, decomposition-based TMOs^3,7–10 are very similar to the proposed one. Two noticeable methods among them are selected for comparsion: that of Durand and Dorsey³ and that of Farbman et al.⁷ The evaluation metric by Aydin et al. is once again employed to evaluate the performance of the four algorithms.

Although parameters can greatly affect the final results, we use the recommended parameters given by the authors when producing the results of the other three algorithms. Some results of the proposed TMO and the other three, together with the corresponding distortion maps, are illustrated in Fig. 12. In the comparison experiment, a test set is formed comprising 10 HDR radiance maps mentioned in the previous texts. Finally, a perspective summation of the numerical results on the test set is given in Table 2, where the numbers are obtained by averaging the error percentages over the test set. The numbers in boldface denote the best ones in the corresponding columns. Although the averaged loss error percentage of the proposed algorithm is greater than that of Ref. 7, the other three terms are all smaller than that of the other methods. The normalized average error rates of the four algorithms are plotted together in Fig. 13 as well. Results show that the proposed algorithm outperforms the other three methods.

Fig. 12

Comparison results using the dynamic-range-independent metric in Ref. 13. First column: tone-mapped results produced by the proposed algorithm. Second column: results produced by Sobol.³¹ Third column: results produced by Durand and Dorsey.³ Fourth column: results produced by Farbman et al.⁷ Images in the first, the third and the fifth rows are tone-mapped results, while their corresponding distortion maps are illustrated in the second, the fourth and the sixth rows, respectively. First row: Nave, HDR radiance map courtesy of Paul Debevec. Third row: Rosette, HDR radiance map courtesy of Paul Debevec. Fifth row: Tree, HDR radiance map courtesy of ILM.

Table 2

Error percentages, averaged over the test set.

Fig. 13

Average error rate: the horizontal axis denotes the average error rate. The vertical axis is the index of the HDR radiance maps, where the names of the corresponding HDR radiance maps are given on the left side. It is shown that the proposed operator produces less error than the other three methods.

6.

Conclusions

In this paper, inspired by anisotropic diffusion, an edge-preserving upper envelope estimator is proposed based on the McCann–Sobel iterative Retinex algorithm. With several modifications, the iterative Retinex computation is made more effective for HDR tone mapping. The proposed modifications include the jumping-spiral iteration manner for better symmetry of edge response over directions, the data-dependent smoothing for edge preservation when estimating the illumination, and illumination-reflectance decomposition for HDR tone mapping. The effectiveness of the proposed method was evaluated using established protocols, and encouraging results are produced compared with some related methods.