1 April 2013 Adapting iterative retinex computation for high-dynamic-range tone mapping
Author Affiliations +
J. of Electronic Imaging, 22(2), 023006 (2013). doi:10.1117/1.JEI.22.2.023006
Retinex algorithms have been widely applied in many aspects of image processing. Based on the iterative Retinex algorithm, we propose an edge-preserving illumination estimation method. Inspired by the anisotropic diffusion, an edge-stopping function is introduced in the iterative computation. This modification enables the preservation of abrupt edges when computing the upper envelope of a given image. Based on the illumination-reflectance decomposition, a high-dynamic-range (HDR) radiance map can be easily tone-mapped to be a low-dynamic-range image by compressing the range of the estimated illumination. Artifacts are effectively suppressed using the proposed method. Meanwhile, we also propose a jumping-spiral iteration manner to improve the symmetry of the edge response. Experimental results show that the proposed tone mapping algorithm is very effective in reproducing HDR scenes, and has a better performance compared with some similar operators.
Pan, An, and He: Adapting iterative retinex computation for high-dynamic-range tone mapping



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.34.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,78.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 algorithm11 is combined with the anisotropic diffusion12 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.


Related Works


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.1415.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.2324.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,2829.30.31 and shadow removal,3233.34 among others.


Retinex-Based TMOs

Sobol31 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 Susstrunk30 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.


Decomposition-Based TMOs

There are a number of previous works for HDR tone mapping based on multiscale image decomposition.3,78.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 Dorsey3 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 L0 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.


Computational Background



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:


where i=log(I), r=log(R), and l=log(L), respectively. Since objects never reflect more energy than the incident light, many Retinex algorithms discount the illumination estimated using the upper envelope of the original image,5,11,23,36 where the estimated illumination satisfies li for all pixels in the image domain.


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:


where p is the neighborhood center pixel index and q is the index of one of the neighboring points. R(·) is the reset operation, which truncates the data exceeding a preset upper bound. The computation proceeds in a inward spiral manner, where only the neighboring pixels in four aligned directions are processed. Equivalently, Eq. (2) can be rewritten as


where Cw is the preset upper bound, which is referred as the perceived lightness of the brightest points in the scene.

Cooper and Baqai37 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 Cw=0 in Eq. (3). Consequently, the McCann–Sobel algorithm is a special case of the Frankle–McCann algorithm.


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 max{·} 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.


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.


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×128 monochrome image. The size of the bright box in the center of the image is 32×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:


where I and If 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



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



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 Malik12 can be formulated as follows:


where λ is the step size of the discrete time steps. ηp denotes the spatial neighborhood of the pixel at p, and |η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:



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

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



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:



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:



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



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


where C{R,G,B} denotes the three independent color components in RGB space. Lin and Lout denote the luminance before and after tone mapping, respectively. And S is a parameter controlling the saturation of the reproduced color.


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, γ=5, σ=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, ×10, ×100, and ×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 γ 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 γ in Eq. (17). The value of σ 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. 910. 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) γ=2, 4, 8, and 12, respectively, where σ=1, n=1 and S=0.8. Larger γ leads to stronger compression of the dynamic range, while smaller γ 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 σ=0.1, 0.5, 2.5, and 12.5, respectively, where γ=5, n=1 and S=0.8. Halos are eliminated with smaller σ, but details may diminish in the tone-mapped result and the ringing artifacts may be elevated at the same time. On the contrary, larger σ produces better overall contrast, but the halos may become visible as σ 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 γ=5, σ=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.3839.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 γ and σ, 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 El, Ea, Er, and Eavg, 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.

S0.8Color saturation
γ4Illumination discounting strength
σ1.5Sigma of the edge-stopping function
n1Number of iterations

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 TMOs3,78.9.10 are very similar to the proposed one. Two noticeable methods among them are selected for comparsion: that of Durand and Dorsey3 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.

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.31 Third column: results produced by Durand and Dorsey.3 Fourth column: results produced by Farbman et al.7 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.




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.


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.


1. S. Ferradanset al., “An analysis of visual adaptation and contrast perception for tone mapping,” IEEE Trans. Pattern Anal. Mach. Intell. 33(10), 2002–2012 (2011).ITPIDJ0162-8828 http://dx.doi.org/10.1109/TPAMI.2011.46 Google Scholar

2. P. LeddaL. P. SantosA. Chalmers, “A local model of eye adaptation for high dynamic range images,” in Proc. 3rd Int. Conf. Computer Graphics, Virtual Reality, Visualisation and Interaction in Africa, AFRIGRAPH ‘04, pp. 151–160, ACM, New York (2004). Google Scholar

3. F. DurandJ. Dorsey, “Fast bilateral filtering for the display of high-dynamic-range images,” ACM Trans. Graph. 21(3), 257–266 (2002).ATGRDF0730-0301 http://dx.doi.org/10.1145/566654.566574 Google Scholar

4. E. Reinhardet al., High Dynamic Range Imaging: Acquisition, Display, and Image-Based Lighting, Morgan Kaumann, The Netherlands (2006). Google Scholar

5. D. ShakedR. Keshet, “Robust recursive envelope operators for fast Retinex,” Technical Report No. HPL-2002-74 (2002). Google Scholar

6. K. Devlin, “A review of tone reproduction techniques,” Technical Report No. CSTR-02-005 (2002). Google Scholar

7. Z. Farbmanet al., “Edge-preserving decompositions for multi-scale tone and detail manipulation,” ACM Trans. Graph. 27(3), 67:1–67:10 (2008).ATGRDF0730-0301 http://dx.doi.org/10.1145/1360612 Google Scholar

8. K. SubrC. SolerF. Durand, “Edge-preserving multiscale image decomposition based on local extrema,” ACM Trans. Graph. 28(5), 147:1–147:9 (2009).ATGRDF0730-0301 http://dx.doi.org/10.1145/1618452 Google Scholar

9. J. TumblinG. Turk, “LCIS: a boundary hierarchy for detail-preserving contrast reduction,” in Proc. 26th Annual Conf. Computer Graphics and Interactive Techniques, SIGGRAPH ‘99, pp. 83–90, ACM Press/Addison-Wesley Publishing Co., New York (1999). Google Scholar

10. L. Xuet al., “Image smoothing via L0 gradient minimization,” ACM Trans. Graph. 30(6), 174:1–174:12 (2011).ATGRDF0730-0301 http://dx.doi.org/10.1145/2070781.2024208 Google Scholar

11. J. J. McCannI. Sobel, “Experiments with Retinex,” HPL Color Summit, Hewlett Packard Laboratories, Israel (1998). Google Scholar

12. P. PeronaJ. Malik, “Scale-space and edge detection using anisotropic diffusion,” IEEE Trans. Pattern Anal. Mach. Intell. 12(7), 629–639 (1990).ITPIDJ0162-8828 http://dx.doi.org/10.1109/34.56205 Google Scholar

13. T. O. Aydinet al., “Dynamic range independent image quality assessment,” ACM Trans. Graph. 27(3), 69:1–69:10 (2008). http://dx.doi.org/10.1145/1360612.1360668 Google Scholar

14. E. H. Land, “The Retinex theory of color vision,” Sci. Am. 237(6), 108–128 (1977).SCAMAC0036-8733 http://dx.doi.org/10.1038/scientificamerican1277-108 Google Scholar

15. E. H. LandJ. J. McCann, “Lightness and the Retinex theory,” J. Opt. Soc. Am. 61(1), 1–11 (1971).JOSAAH0030-3941 http://dx.doi.org/10.1364/JOSA.61.000001 Google Scholar

16. E. Provenziet al., “Mathematical definition and analysis of the retinex algorithm,” J. Opt. Soc. Am. A 22(12), 2613–2621 (2005).JOAOD60740-3232 http://dx.doi.org/10.1364/JOSAA.22.002613 Google Scholar

17. E. H. Land, “Recent advances in Retinex theory and some implications for cortical computations: color vision and the natural image,” Proc. Natl. Acad. Sci. 80(16), 5163–5169 (1983).PIPSBD0370-0046 http://dx.doi.org/10.1073/pnas.80.16.5163 Google Scholar

18. E. Provenziet al., “Random spray Retinex: a new Retinex implementation to investigate the local properties of the model,” IEEE Trans. Image Process. 16(1), 162–171 (2007).IIPRE41057-7149 http://dx.doi.org/10.1109/TIP.2006.884946 Google Scholar

19. J. A. FrankleJ. J. McCann, “Method and apparatus for lightness imaging,” U.S. Patent No. 4384336 (1983). Google Scholar

20. B. FuntF. CiureaJ. McCann, “Retinex in MATLAB™,” J. Electron. Imag. 13(1), 48–57 (2004).JEIME51017-9909 http://dx.doi.org/10.1117/1.1636761 Google Scholar

21. D. JobsonZ. RahmanG. Woodell, “A multiscale Retinex for bridging the gap between color images and the human observation of scenes,” IEEE Trans. Image Process. 6(7), 965–976 (1997).IIPRE41057-7149 http://dx.doi.org/10.1109/83.597272 Google Scholar

22. E. H. Land, “An alternative technique for the computation of the designator in the Retinex theory of color vision,” Proc. Natl. Acad. Sci. 83(10), 3078–3080 (1986).PIPSBD0370-0046 http://dx.doi.org/10.1073/pnas.83.10.3078 Google Scholar

23. R. Kimmelet al., “A variational framework for Retinex,” Int. J. Comput. Vis. 52(1), 7–23 (2003).IJCVEQ0920-5691 http://dx.doi.org/10.1023/A:1022314423998 Google Scholar

24. W. Maet al., “An l1-based variational model for Retinex theory and its application to medical images,” in Proc. IEEE Conf. Comput. Vis. and Pattern Recognit., pp. 153–160, IEEE Computer Society, Providence, Rhode Island (2011). Google Scholar

25. J. MorelA. PetroC. Sbert, “A PDE formalization of Retinex theory,” IEEE Trans. Image Process. 19(11), 2825–2837 (2010).IIPRE41057-7149 http://dx.doi.org/10.1109/TIP.2010.2049239 Google Scholar

26. H. TakahashiT. SaitoT. Komatsu, “Variational Retinex algorithm with its application to a high-quality chroma key,” in Proc. IEEE Int. Conf. Image Process., pp. 977–980, IEEE Signal Processing Society, Atlanta, Georgia (2006). Google Scholar

27. C.-T. ShenW.-L. Hwang, “Color image enhancement using Retinex with robust envelope,” in Proc. 16th IEEE Int. Conf. Image Process., pp. 3141–3144, IEEE Signal Processing Society, Cairo (2009). Google Scholar

28. F. Dragoet al., “Design of a tone mapping operator for high-dynamic-range images based upon psychophysical evaluation and preference mapping,” Proc. SPIE 5007, 321–331 (2003). http://dx.doi.org/10.1117/12.473919 Google Scholar

29. K. KimJ. BaeJ. Kim, “Natural HDR image tone mapping based on Retinex,” IEEE Trans. Consum. Electron. 57(4), 1807–1814 (2011).ITCEDA0098-3063 http://dx.doi.org/10.1109/TCE.2011.6131157 Google Scholar

30. L. MeylanS. Susstrunk, “High dynamic range image rendering with a Retinex-based adaptive filter,” IEEE Trans. Image Process. 15(9), 2820–2830 (2006).IIPRE41057-7149 http://dx.doi.org/10.1109/TIP.2006.877312 Google Scholar

31. R. Sobol, “Improving the Retinex algorithm for rendering wide dynamic range photographs,” J. Electron. Imag. 13(1), 65–74 (2004).JEIME51017-9909 http://dx.doi.org/10.1117/1.1636762 Google Scholar

32. G. D. FinlaysonS. D. HordleyM. Drew, “Removing shadows from images using Retinex,” in Proc. of IS&T/SID 10th Color Imaging Conf., pp. 73–79, IS&T, Scottsdale, Arizona (2002). Google Scholar

33. G. MaJ. Yang, “Shadow removal using Retinex theory,” in Proc. 3rd Chinese Conf. Intelligent Visual Surveillance pp. 25–28, CASIA, Beijing, China (2011). Google Scholar

34. J. SunY. DuY. Tang, “Shadow detection and removal from solo natural image based on retinex theory,” in Proc. 1st Int. Conf. on Intelligent Robotics and Applications: Part I, ICIRA ‘08, pp. 660–668, Springer-Verlag, Berlin, Heidelberg (2008). Google Scholar

35. P. ChoudhuryJ. Tumblin, “The trilateral filter for high contrast images and meshes,” in Proc. 14th Eurographics Workshop on Rendering, EGRW ‘03, pp. 186–196, Aire-la-Ville, Switzerland (2003). Google Scholar

36. M. Elad, “Retinex by two bilateral filters,” in Scale Space and PDE Methods in Computer Vision, R. KimmelN. SochenJ. Weickert, Eds., vol. 3459 of Lecture Notes in Computer Science, pp. 217–229, Springer, Berlin/Heidelberg (2005). Google Scholar

37. T. J. CooperF. A. Baqai, “Analysis and extensions of the Frankle-McCann Retinex algorithm,” J. Electron. Imag. 13(1), 85–92 (2004).JEIME51017-9909 http://dx.doi.org/10.1117/1.1636182 Google Scholar

38. J. Kuanget al., “Evaluating HDR rendering algorithms,” ACM Trans. Appl. Percept. 4(2), 1–27 (2007).1544-3558 http://dx.doi.org/10.1145/1265957 Google Scholar

39. P. Leddaet al., “Evaluation of tone mapping operators using a high dynamic range display,” ACM Trans. Graph. 24(3), 640–648 (2005).ATGRDF0730-0301 http://dx.doi.org/10.1145/1073204 Google Scholar

40. A. Yoshidaet al., “Perceptual evaluation of tone mapping operators with real-world scenes,” Proc. SPIE 5666, 192–203 (2005). Google Scholar

41. M. Cadiket al., “Evaluation of HDR tone mapping methods using essential perceptual attributes,” Comput. Graph. 32(3), 330–349 (2008).0097-8493 http://dx.doi.org/10.1016/j.cag.2008.04.003 Google Scholar



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.

© The Authors. Published by SPIE under a Creative Commons Attribution 3.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.
Shengdong Pan, Xiangjing An, Hangen He, "Adapting iterative retinex computation for high-dynamic-range tone mapping," Journal of Electronic Imaging 22(2), 023006 (1 April 2013). https://doi.org/10.1117/1.JEI.22.2.023006

Back to Top