## 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}4.^{–}^{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}8.9.^{–}^{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^{11} is combined with the anisotropic diffusion^{12} 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.^{13}

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

## 2.1.

### Retinex Theory and Its Variants

Land’s Retinex theory is one of the most famous attempts to explain and simulate the perception of lightness and color of the HVS.^{14}15.^{–}^{16} After about 50 years of evolution, many different variants are proposed. According to their implementation forms, the standard Retinex algorithms can be roughly categorized as path-based algorithms,^{15}^{,}^{17}^{,}^{18} iterative algorithms,^{11}^{,}^{19}^{,}^{20} center/surround algorithms,^{21}^{,}^{22} and PDE-type algorithms.^{23}24.25.^{–}^{26}

Together with many extensions for other purposes, the huge family of Retinex-type algorithms are now widely applied in many aspects of image processing, such as color constancy,^{17}^{,}^{18}^{,}^{20}^{,}^{25} general image enhancement,^{5}^{,}^{21}^{,}^{27} HDR image rendering,^{28}29.30.^{–}^{31} and shadow removal,^{32}33.^{–}^{34} among others.

## 2.2.

### Retinex-Based TMOs

Sobol^{31} introduces several effective improvements based on the Frankle–McCann Retinex algorithm, including the ratio modification operator (RMO), the spatially varying contrast gain, and the partial contrast strength mask. The RMO truncates the intermediate values of the contrast ratios to a limited range. Consequently, the dynamic range of the original image is reduced, while fine details are preserved. Similar to Sobol’s RMO, Drago et al.^{28} propose a soft contrast clipping function to automatically manipulate local contrast and suppress artifacts for tone mapping. However, the contrast clipping operation of these two TMOs changes the nature of the contrast of the original intensity. As a result, the appearance of the scene may be changed significantly.

Based on the center/surround Retinex, Meylan and Susstrunk^{30} introduce the adaptive filtering in computing a surround-like contrast mask. More recently, Kim et al.^{29} apply a single-scale Retinex on intermediate data to enhance the naturalness of the result images. In both methods, a large support of the surround function is used. Though methods are applied to reduce the complexity, the computational burden of these TMOs is relatively heavy.

## 2.3.

### Decomposition-Based TMOs

There are a number of previous works for HDR tone mapping based on multiscale image decomposition.^{3}^{,}^{7}8.9.^{–}^{10}^{,}^{35} The original image is often decomposed into several layers corresponding to different scales. The dynamic range of the output is manipulated by reweighting these layers.

In order to avoid artifacts in tone mapping using this type of operators, edge-preserving local smoothing is often employed in pursuing the base layer. Durand and Dorsey^{3} propose a very efficient approximate of the bilateral filter to compute the piece-wise smooth base layer. Farbman et al.^{7} propose an effective edge-preserving operator based on the weighted-least-square optimization. Xu et al.^{10} employ ${\mathcal{L}}_{0}$ gradient minimization in computing the base layer. However, the number of layers is arbitrary, and no indicative selection of the weights is presented for perceptual matching in anyone of these works.

## 3.

## Computational Background

## 3.1.

### Preliminaries

The Retinex theory is employed to compensate the effects of uneven illumination. Following some previous works,^{11}^{,}^{23}^{,}^{24} a given image $I$ can be decomposed into the reflectance $R$ and the illumination $L$, such that $I=RL$. The decomposition is often done in logarithmic space, such that:

^{5}

^{,}

^{11}

^{,}

^{23}

^{,}

^{36}where the estimated illumination satisfies $l\ge i$ for all pixels in the image domain.

## 3.2.

### Iterative Retinex Computation

The first iterative Retinex algorithm is an efficient variant of the Retinex algorithm, which was developed by Frankle and McCann.^{19} It is generalized by Funt et al.^{20} to handle images with arbitrary resolution. The formulation of the generalized Frankle–McCann algorithm is given as follows:

Cooper and Baqai^{37} propose some extensions to the Frankle–McCann Retinex algorithm. They introduce the distance weighting to enhance local contrast, the soft reset to reduce halos, and the dual-spiral iteration manner to improve the symmetry of the spatial response.

McCann and Sobel adapt the iterative scheme to first estimate the illumination using the upper envelope of the original image.^{11} The iterative update operation of the McCann–Sobel Retinex algorithm is given as below:

Combined with Eq. (1), the above equation can be rewritten as

Note that Eq. (5) is exactly the formulation of the Frankle–McCann Retinex algorithm when ${C}_{w}=0$ in Eq. (3). Consequently, the McCann–Sobel algorithm is a special case of the Frankle–McCann algorithm.

## 3.3.

### Another Look into the Iterative Retinex Computation

Consider the formulation of the McCann–Sobel Retinex algorithm given as Eq. (4). Without the nonlinear operation for the upper envelope, it can be simplified as follows:

Alternatively, Eq. (6) can be rewritten in a differential form:

This can be interpreted as directional isotropic diffusion. Its asymmetric effect is compensated by the spiral iteration manner. Combining the progressively shortened shift distance, Eq. (7) can be considered as an efficient approximate of the discrete isotropic diffusion with isometric step size. The initial shift distance, together with the number of iteration, controls the number of diffusion steps, and further controls the support size and the scale of the smoothing kernel.

Furthermore, the nonlinear operations are introduced to truncate the intermediate data to fit a preset interval. The reset operation restricts that the value of the estimated reflectance does not exceed the maximum brightness. Similarly, the operator $\mathrm{max}\{\xb7\}$ of the McCann–Sobel algorithm corresponds to the physical constraint that the objects never reflect more energy than the incident light. Thus the illumination estimation is a Gaussian-like upper envelope of the original image.

## 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.

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.

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

where $I$ and ${I}_{f}$ are the square matrix to be tested and its flipped version, respectively. $E(\xb7)$ is an error norm measuring the diagonal symmetry of a square matrix, which is defined asSince 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.

## 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^{12} can be formulated as follows:

## (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}),$$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:

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:

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).

## 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:

## (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:

## (14)

$${l}_{p}^{k+1}=[1-\frac{g(|{i}_{q}-{i}_{p}|)}{2}]{l}_{p}^{k}+\frac{g(|{i}_{q}-{i}_{p}|)}{2}\mathrm{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:

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.

## 4.4.

### Tone Mapping Based on the Proposed Algorithm

In order to avoid gray-world violation problems when handling the color images, the original image in RGB color space is first converted into a grayscale image. The color-to-gray conversion model employed in this paper is given as follows:

Then the illumination estimation of $I$ is computed using the proposed upper envelope estimator. As illustrated in Fig. 5, the proposed variant of the iterative Retinex algorithm is able to decouple the HDR illumination and the LDR reflectance. Unlike the standard Retinex algorithms, the estimated reflectance is not directly taken as the output, but relit using the range-compressed version of the estimated illumination. For simplicity, the range of the illumination is reduced using a Gamma mapping. The modified intensity can be computed as follows:

where $\gamma >1$, which can be tuned to obtain various strength of contrast reduction.Finally, the color information is reproduced according to the original HDR radiance map. In color reproduction, we follow the method employed in Refs. 3, 9, and 29, which can be formulated as follows:

## (18)

$${C}_{\mathrm{out}}={\left(\frac{{C}_{\mathrm{in}}}{{L}_{\mathrm{in}}}\right)}^{S}{L}_{\mathrm{out}},$$## 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.

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.

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}39.^{–}^{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.^{41} propose a method based on a basic set of attributes, which is built on the statistics of subjective experiments. Aydin et al.^{13} 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.

Parameter | Value | Description |
---|---|---|

S | 0.8 | Color saturation |

γ | 4 | Illumination discounting strength |

σ | 1.5 | Sigma of the edge-stopping function |

n | 1 | Number of iterations |

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}8.9.^{–}^{10} are very similar to the proposed one. Two noticeable methods among them are selected for comparsion: that of Durand and Dorsey^{3} and that of Farbman et al.^{7} 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.

## Table 2

Error percentages, averaged over the test set.

Algorithms | El | Ea | Er | Eavg |
---|---|---|---|---|

Ours | 20.43 | 3.21 | 7.18 | 10.27 |

Sobol04 | 39.69 | 20.61 | 19.83 | 26.71 |

Durand02 | 53.85 | 6.53 | 21.52 | 27.30 |

Farbman08 | 13.48 | 12.67 | 17.62 | 14.59 |

## 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.

## Acknowledgments

The authors express their warm thanks to T. Aydin and R. Mantiuk for their efforts on providing the online access to the dynamic range independent image quality assessment metric. This work is partially supported by the National Natural Science Foundation of China, under Grant Nos. 61075043 and 90820302.

## References

## Biography

**Shengdong Pan** received his BE and ME degrees in control science and engineering from National University of Defense Technology (NUDT), Changsha, Hunan, China, in 2006 and 2008, respectively. He is currently working toward his PhD degree in NUDT. His research interests include image acquisition and processing, high dynamic range compression, computer vision, and pattern recognition.

**Xiangjing An** received his BS degree in automatic control from the Department of Automatic Control, National University of Defense Technology (NUDT), Changsha, China, in 1995 and his PhD degree in control science and engineering from the College of Mechatronics and Automation (CMA), NUDT in 2001. He has been a visiting scholar for cooperation research in the Boston University during 2009 to 2010. Currently, he is an associate professor at the Institute of Automation, CMA, NUDT. His research interests include image processing, computer vision, biologically inspired feature extraction.

**Hangen He** received his BSc degree in nuclear physics from Harbin Engineering Institute, Harbin, China, in 1968. He was a visiting professor at the University of the German Federal Armed Forces in 1996 and 1999, respectively. He is currently a professor in the College of Mechatronics and Automation (CMA), National University of Defense Technology (NUDT), Changsha, Hunan, China. His research interests include artificial intelligence, reinforcement learning, learning control, and robotics. He has served as a member of editorial boards of several journals and has cochaired many professional conferences. He is a joint recipient of more than a dozen academic awards in China.