{'fraction_non_alphanumeric': 0.028154967395473726, 'fraction_numerical': 0.0076332949750671265, 'mean_word_length': 4.940077466393256, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Adapting iterative retinex computation for high-dynamic-range tone mapping

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


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.

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.

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

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: 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 l ≥ i 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 C w is the preset upper bound, which is referred as the perceived lightness of the brightest points in the scene. 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)

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 maxf·g 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.
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.
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ð·Þ 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 n ≥ 2. Consequently, the symmetry of the spatial response of the spiral iteration can be further improved using the jumping-spiral manner.

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: where λ is the step size of the discrete time steps. η p denotes the spatial neighborhood of the pixel at p, and jη p j is the number of neighbors. gð·Þ is the edge-stopping function, which is generally a Gaussian function. The anisotropic diffusion is very effective in edgepreserving 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).

Edge-Preserving Illumination Estimation
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: gðji q − i p jÞ 2 maxfi p ; l k q g: (14) It is experienced that local structures can be oversmoothed in the illumination estimation with a small number of iterations, and thus the corresponding reflectance is overexaggerated. 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: wherel 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.

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 colorto-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 ∈ fR; G; Bg denotes the three independent color components in RGB space. L in and L out 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 tonemapped 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.
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. [8][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 realtime HDR imaging pipelines. Meanwhile, the prototypes of the HDR display are now still too expensive for lowcost 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, 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.
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 E l , E a , E r , and E 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.
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 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.

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