Ba Chien Thai

Anissa Mokraoui
email: anissa.mokraoui@univ-paris13.frmatei@lipn.univ-paris13.fr

Basarab Matei

L2ti

HDR Image Tone Mapping Histogram Adjustment with Using An Optimized Contrast Parameter

Keywords: High Dynamic Range image, Low Dynamic Range image, Tone mapping, Lifting scheme, Optimized prediction

This paper proposes a Tone Mapping (TM) approach converting a High Dynamic Range (HDR) image into a Low Dynamic Range (LDR) image while preserving as much information of the HDR image as possible to ensure a good LDR image visual quality. Multiresolution approach with approximation and details are separated, and weighted by entropy, and used an optimized contrast parameter. The parameter relies on all accounts of luminance average themselves in order to enhance contrast as much as possible. Then, the pixel's distribution of the coarse reconstructed LDR image is then adjusted according to a perceptual quantizer with respect to the human visual system using a piecewise linear function. Simulation results provide good results, both in terms of visual quality and TMQI metric, compared to existing competitive TM approaches.

I. INTRODUCTION

The objective of a High Dynamic Range (HDR) image Tone Mapping (TM) approach is to find a trade-off between the relevant information (e.g. details, contrast, brightness...) to be preserved or discarded in the image ensuring a good visual quality of the displayed image on Low Dynamic Range (LDR) devices that would be appreciated by observers.

A state of the art on HDR image TM approaches is fairly complete in [START_REF] Banterle | Advanced High Dynamic Range Imaging: Theory and Practice[END_REF], [START_REF] Dufaux | High Dynamic Range Video 1st Edition : From Acquisition, to Display and Applications[END_REF] and [START_REF] Reinhard | High Dynamic Range Imaging 2nd Edition: Acquisition, Display, and Image-Based Lighting[END_REF] where a classification is provided. Among the developed TM strategies, this paper quickly reviews those that caught our attention because of their performance. In [START_REF] Durand | Fast Bilateral Filtering for The Display of High-dynamic-range Images[END_REF], the TM approach reduces the HDR contrast while preserving the image details. This work uses an edge-preserving bilateral filter to decompose the HDR image into two layers: a base layer encoding large-scale variations and a detail one. Contrast is then reduced only in the first layer while the details are kept unchanged. An adaptive logarithmic mapping method of luminance values is presented in [START_REF] Drago | Adaptive Logarithmic Mapping for Displaying High Contrast Scenes[END_REF]. It concerns the adjustment of the logarithmic basis depending on the radiance of the pixels. In [START_REF] Li | Compressing and Companding High Dynamic Range Images with Subband Architectures[END_REF], a subband architecture related on an oversampled Haar pyramid representation is proposed. Subband coefficients are re-scaled according to a gain control function reducing the high frequency magnitudes and boosting low ones. In [START_REF] Duan | Tone-mapping High Dynamic Range Images by Novel Histogram Adjustment[END_REF], a TM optimization approach using a histogram adjustment between linear mapping and the equalized histogram mapping is developed. A modification of this approach is made where revisited histogram equalization approaches are discussed in [START_REF] Husseis | Revisited Histogram Equalization as HDR Images Tone Mapping Operators[END_REF]. The latter considers both histogram equalization and human sensitivity to the light function.

In [START_REF] Fattal | Edge-Avoiding Wavelets and their Applications[END_REF], a second generation of wavelets based on the edge content of the image avoiding having pixels from both sides of an edge is proposed. In [START_REF] Thai | Performance Evaluation of High Dynamic Range Image Tone Mapping Operators Based on Separable Non-linear Multiresolution Families[END_REF], a separable non-linear multiresolution approach based on essentially non-oscillatory interpolation strategy has been investigated. These approaches take into account the singularities (e.g. edge points) in their mathematical models preserving then the structural information of the HDR images. In [START_REF] Thai | Image Tone Mapping Approach Using Essentially Non-Oscillatory Bi-quadratic Interpolations Combined with a Weighting Coefficients Strategy[END_REF], a non-separable non-linear multiresolution approach is proposed. The results provided in [START_REF] Fattal | Edge-Avoiding Wavelets and their Applications[END_REF], [START_REF] Thai | Performance Evaluation of High Dynamic Range Image Tone Mapping Operators Based on Separable Non-linear Multiresolution Families[END_REF] and [START_REF] Thai | Image Tone Mapping Approach Using Essentially Non-Oscillatory Bi-quadratic Interpolations Combined with a Weighting Coefficients Strategy[END_REF] show that the decomposition of the HDR image on different resolution levels would seem to be a good strategy. However, the choice of the decomposition filters is extremely important and decisive in the extractive power of the approximation and detail information. To do so, this paper proposes a separable near optimal lifting scheme. The "Predict" operation is performed by a new adaptive prediction operation to extract the finest detail coefficients to highlight the sharp transitions in the given HDR image.

The proposed HDR image TM approach has two main goals namely the preservation of the details and the adjustment of the contrast in accordance with the LDR display devices. It is composed of four stages. The first one judiciously decomposes the HDR image into different resolution levels (section II-A). The second one concerns the weighting strategy of the detail and approximation coefficients (section II-B). The third one reconstructs the coarse LDR image (section II-C). Finally the fourth one adjusts the contrast according to a perceptual linear quantizer (section II-D). Section III discusses the simulation results. Section IV concludes the paper.

II. PROPOSED HDR IMAGE TONE MAPPING APPROACH

This section concerns the first, second, third and fourth stages of the proposed HDR image TM mapping approach (i.e. decomposition, weighting, reconstruction and contrast adjustment stages).

Before describing the approach, introduce some notations. The original HDR image, at the finest resolution level J, is assumed to be of size N J × M J . The index j refers to the resolution level (with j = 0, ..., J -1). Denote l HDR the HDR image luminance. In the rest of this paper, the HDR image luminance is considered in the logarithm domain since it is well adapted to the human visual system. It is denoted I J and defined as follows I J = {I J (x n , y m ) = log 10 (l HDR (x n , y m )) for 1 ≤ n ≤ N J and 1 ≤ m ≤ M J where I J (x n , y m ) is the HDR logarithm luminance value of the pixel located at position (x n , y m ) on the image.

A. First stage: HDR image decomposition

The proposed HDR image decomposition is performed according to the forward process of a separable near optimal cell-average lifting scheme. This strategy is motivated by the fact that the relevant details are accurately predicted since the coefficients of the filters are adapted locally to the image to be processed.

The decomposition consists to go from the finest resolution level J to the coarse resolution level 0. At a given resolution level j -1 (with 0 ≤ j ≤ J -1), the algorithm deals with the approximation coefficients denoted I j (x n , y k ) (with 1 ≤ n ≤ N j and 1 ≤ k ≤ M j ) computed at resolution level j. For a given n belonging to [1, N j ], the algorithm starts with splitting the horizontal 1D-signal, i.e. I j (x n , y k ) for 1 ≤ k ≤ M j , into a set of odd and even indexes as follows:

{I j (x n , y k ) with 1 ≤ k ≤ M j } := {I j (x n , y 2k-1 ), I j (x n , y 2k ) with 1 ≤ k ≤ M j /2 }. ( 1
)
This process is then repeated for all n. Based on this split, the approximation coefficient located at position (x n , y k ), denoted V j-1 (x n , y k ), is computed on a Cell-Average (CA) scheme. For a given n, it is expressed as follows:

V j-1 (x n , y k ) = I j (x n , y 2k-1 ) + I j (x n , y 2k ) 2 for 1 ≤ k ≤ M j /2. ( 2
)
This process is then repeated for all n.

The detail coefficient, denoted d j-1 (x n , y 2k-1 ), is computed at odd indexes and is provided below for a given n:

d j-1 (x n , y 2k-1 ) = Îj (x n , y 2k-1 ) -I j (x n , y 2k-1 ) for 1 ≤ k ≤ M j /2, (3)
where Îj (x n , y 2k-1 ) is the predicted logarithm luminance value at odd position (x n , y 2k-1 ) and resolution level j. This predicted value is expressed as a linear weighted combination of the neighboring approximation coefficients at resolution level j -1:

Îj (x n , y 2k-1 ) = 2 i=0 u j-1 i (x n ) • V j-1 (x n , y k+i-1 ) for 1 ≤ k ≤ M j /2, (4)
where the weights u j-1 i (x n ) must preserve the initial 1Dsignal average which results in satisfying this condition:

2 i=0 u j-1 i (x n ) = 1. (5)
In what follows, the weights are rather denoted in a vector form

u j-1 = (u j-1 0 (x n ), u j-1 1 (x n ), u j-1 2 (x n ))
t to lighten the writing. These weights are deduced so that the Mean Squared Error (MSE) between Îj (x n , y 2k-1 ) and I j (x n , y 2k-1 ) is minimized:

E = arg min u j-1 i Îj (x n , y 2k-1 ) -I j (x n , y 2k-1 ) 2 2 for 1 ≤ k ≤ M j /2. (6)
This results in solving the following equation:

Γ j-1 • u j-1 = r j-1 , (7)
where u j-1 is the weight vector. r j-1 is the crosscorrelation vector r j-1 = (r j-1 (0), r j-1 (1), r j-1 (2)) t where r j-1 (i) represents the cross-correlation function between V j-1 (x n , y k+i-1 ) and

I j (x n , y 2k-1 ) for 1 ≤ k ≤ M j /2.
Γ j-1 is the autocorrelation matrix defined as:

Γ j-1 =   R j-1 (0) R j-1 (1) R j-1 (2) R j-1 (-1) R j-1 (0) R j-1 (1) R j-1 (-2) R j-1 (-1) R j-1 (0)   , (8)
where R j-1 (i) is the autocorrelation of V j-1 (x n , y k ). These weights, associated to the row (x n ), are then deduced so that the partial derivatives of the MSE (given by equation ( 6)) with respect to u j-1 i (i.e. i = 0, 1, 2) are equal to zero:

u j-1 (x n ) = (Γ j-1 ) -1 • r j-1 . (9)
The weight vectors are then stored to be used in the adaptive lifting scheme backward process to reconstruct the decomposed 1D-signal.

For sake of convenience, the approximation and details coefficients are organized as follows:

{W j (x k , y m ) for 1 ≤ k ≤ N j and ∀m ∈ [1, M j /2]} := {V j-1 (x n , y k ) for 1 ≤ k ≤ M j /2 and ∀n ∈ [1, N j ]}, {U j (x k , y m ) for 1 ≤ k ≤ N j , ∀m ∈ [1, M j /2]} := {d j-1 (x n , y 2k-1 ) for 1 ≤ k ≤ M j /2, ∀n ∈ [1, N j ]}.
(10) The split, approximation and detail operations are applied on W j (x k , y m ) and U j (x k , y m ) for a given m ∈ [1, M j /2] (i.e. on the vertical direction). Note that the approximation step requires the prediction of Ŵ j (x 2k-1 , y m ) (respectively

Û j (x 2k-1 , y m )) based on a set of weights v j-1 i (respectively w j-1 i
). These weights need to be stored for the backward process lifting scheme to reconstruct the decomposed image. Finally, the approximation resolution level

I j (x n , y m ) is divided into 4 blocks I j := (I j-1 , d j-1 HL , d j-1 LH , d j-1 HH ).
The decomposition is thus iterated on I j-1 until j = 0.

The finest HDR image I J is then represented by 3J + 1 resolution levels

I J := (I 0 , d 0 HL , d 0 LH , d 0 HH , ..., d J-1 HL , d J-1 LH , d J-1 HH ).

B. Second stage: Weighting strategy of the coefficients

To reduce the dynamic range of the HDR image, the TM approach proposes to weight the approximation and detail coefficients in an appropriate way before performing the adaptive lifting scheme backward process described in section II-C.

Denote N l the number of resolution levels equal to J; E a the entropy of the approximation coefficients at the coarsest resolution level (i.e. j = 0); and E j d the entropy of the detail coefficients at resolution level j (i.e. d j HL , d j LH and d j HH ) getting therefore

N l + 1 entropies (E a , E 0 d , E 1 d , ..., E j d ,..., E N l -1 d
). From these entropies, positive weights smaller than one are deduced as follows:

     α a = N l -1 i=0 E i d Ea+ N l -1 i=0 E i d for j = 0 α j d = Ea+ N l -1 i=0, =j E i d Ea+ N l -1 i=0 E i d for j = 0, ..., N l -1 (11)
α a (respectively α j d ) is the weight associated to the approximation (respectively detail) coefficients at resolution level j = 0 (respectively j).

The coefficients of the four coarsest resolution levels are first modified according to:

I 0 = α a × I 0 , d 0 HL = α 0 d × d 0 HL , d 0 LH = α 0 d × d 0 LH , d 0 HH = α 0 d × d 0 HH . (12)
The approximation subband, denoted I 1 , is then reconstructed (see section II-C) and the number of levels is reduced to N l -1. The N l entropies (associated to 3N l -2 subbands) are calculated again to update the weights α a , α 1 d (equation ( 11) with N l = N l -1). After that, these weights are applied on the coefficients I 1 , d 1 HL , d 1 HL , d 1 HL to build I 2 as explained above. This process is iterated until N l = 0 to reconstruct the coarse tone mapped HDR image denoted I J LDR , called coarse LDR image.

C. Third stage: Reconstruction of the coarse LDR image

The reconstruction stage is carried out inversely to the decomposition stage. Assume that the adaptive lifting scheme backward algorithm processed all resolution levels until j -1. The next step consists to recover the approximation coefficients I j of size N j × M j using the four weighted blocks I j-1 , d j-1 LH , d j-1 HL and d j-1 HH . The algorithm first deals with the coefficients in a vertical direction using I j-1 and d j-1 LH

(denoted d j-1 ). At a given m, the approximation coefficient at odd position, denoted W j (x 2k-1 , y m ), is deduced:

W j (x 2k-1 , y m ) = Ŵ j (x 2k-1 , y m ) -d j-1 (x k , y m ) for 1 ≤ k ≤ N j /2, ( 13
)
where Ŵ j (x 2k-1 , y m ) is the predicted coefficient, at odd position, deduced from the weighted combination v j-1 i (with i = 0, 1, 2) of the neighboring approximation coefficients I j-1 (x 2k-1 , y m ) as:

Ŵ j (x 2k-1 , y m ) = 2 i=0 v j-1 i (y m ) • I j-1 (x k+i-1 , y m ) for 1 ≤ k ≤ N j /2. ( 14
)
The weights v j-1 i (y m ) are those computed and stored in the decomposition process. The approximation coefficient W j (x 2k , y m ), at even position, is deduced in a CA scheme:

W j (x 2k , y m ) = 2I j-1 (x k , y m ) -W j (x 2k-1 , y m ) for 1 ≤ k ≤ N j /2. ( 15
)
The odd and even approximation coefficients are then merged to built 1D-signal:

{W j (x k , y m ) for 1 ≤ k ≤ N j } := {W j-1 (x 2k-1 , y m ), W j-1 (x 2k , y m ) for 1 ≤ k ≤ N j /2}. ( 16
) This process is repeated for all m to reconstruct W j of size N j × M j /2. These same steps are applied on d j-1 HL and d j-1 HH to generate the new block U j of size N j × M j /2. As in the decomposition strategy, W j and U j are renamed V j-1 and d j-1 . The same steps, as described above, are then performed but according to a horizontal direction. The reconstruction is iterated to finally build the image called coarse LDR image which is denoted I J LDR .

D. Fourth stage: Piecewise linear perceptual quantizer

This stage proposes to adjust locally the distribution of the coarse LDR image logarithm luminance I J LDR according to the HVS to enhance the contrast using a piecewise linear function. This strategy is inspired from [START_REF] Mai | Visually-favorable tone-mapping with high compression performance[END_REF] which has been developed for compression purpose. However, modifications are made mainly to avoid the problem of empty bins of equal size.

To do so, the I J LDR values are first sorted and classified into equal B bins defined by cutting points denoted c i uLDR (with 1 ≤ i ≤ B). A non-uniform histogram equalization is also performed. c i nuLDR (with 1 ≤ i ≤ B) cutting points, defining the bounds of the non-uniform consecutive B bins, are deduced. The lower bound (cutting point) of each bin is then adjusted as follows:

l i LDR (1) = c i uLDR + β i (c i nuLDR -c i uLDR ), (17)
where β i is a positive parameter smaller than 1 depending on the sub-interval [c i uLDR , c i nuLDR ]. It is computed as the cutting point l i LDR (1) as the average luminance on the subinterval [c i uLDR , c i nuLDR ] in [START_REF] Husseis | Revisited Histogram Equalization as HDR Images Tone Mapping Operators[END_REF]:

β i = mean( I J LDR [c i uLDR , c i nuLDR ]) -c i uLDR c i nuLDR -c i uLDR , (18)
Therefore the I J LDR values are classified into non-uniform B bins as follows:

I J LDR = { l LDR (k) for k = 1, ..., N J × M J } = {[ l 1 LDR (1), ..., l 1 LDR (K 1 )], ..., [ l i LDR (1), ..., l i LDR (K i )], ..., [ l B LDR (1), ..., l B LDR (K B )]}, (19)
depending on the quantization level set, where K i is the number of values in the i -th bin (i.e. 1 ≤ i ≤ B; K i > 0) and satisfying the following relation

B i=1 K i = N J × M J .
The "s-shape" TM perceptual curve, as discussed in [START_REF] Dowling | The Retina: An Approachable Part of the Brain[END_REF] and [START_REF] Geisler | Effects of Bleaching and Backgrounds on the Flash Response of the Cone System[END_REF], is modelled by a piecewise linear curve on each bin (see Fig. 1). Consider the i -th bin, defined by [ l i LDR (1), ..., l i LDR (K i )], the coarse LDR values are then modeled as follows:

li LDR (k) = a i l i LDR (k) + b i with k ∈ [1, K i ], (20)
where a i (with a i = 0) and b i are two unknown parameters depending on the i -th bin. This equation, after some mathematical manipulations, can be rewritten as follows:

l i LDR (k) = li LDR (k) -li LDR (1) a i + l i LDR (1), (21)
where the unknown parameter a i is deduced so that the the norm-M space error between the coarse LDR value and its quantized version, denoted l i LDR (k) (i.e. when li LDR (k) has been supported a ceiling process according to a scalar quantization), is minimized in the i -th bin: arg min

a i l i LDR (k) -l i LDR (k) M M , with M ≥ 1. ( 22
)
This equation is simplified as follows:

arg min

a i Ki k=1 l i LDR (k)-li LDR (k) a i M • p i , ( 23
)
where p i is the probability of the l i LDR (k) value in the i -th bin given by p i = Ki B i=1 Ki . Extending equation ( 23) to all bins involves the computation of a global norm-M space error Ψ deduced as follows:

Ψ = B i=1 Ki k=1 l i LDR (k)-li LDR (k) a i M • p i . ( 24
)
The expected value

Ki k=1 |l i LDR (k) -li LDR (k)| M on
each bin is assumed to be equal and is denoted ξ. Equation ( 24) is then simplified and becomes:

Ψ = B i=1 pi (a i ) M • ξ, with a i > 0. ( 25
)
Denote l LDRmax (respectively l LDRmin ) the maximum (respectively minimum) LDR luminance value. Introduce δ i as the difference between coarse LDR luminance in two consecutive bins:

δ i = li+1 LDR (1) -li LDR (1). (26)
A constraint related to the limit sum of the projected heights equal to the entire LDR range results in:

B i=1 a i • δ i = l LDRmax -l LDRmin . ( 27
)
Therefore the optimization problem is written as follows:

arg min

a i B i=1 p i (a i ) M •ξ, s.t B i=1 a i •δ i = l LDRmax -l LDRmin .
(28) This problem is solved analytically using the Lagrangian function. After some mathematical manipulations, the slope a i is deduced:

a i = (l LDRmax -l LDRmin ) • (p i ) 1/(M +1) B i=1 δ i • (p i ) 1/(M +1) . ( 29
)
Hence the unknown parameter b i is calculated (i.e. b i = li LDR (1)-a i × l i LDR (1)) and LDR mapped values are deduced according to equation (20). The global piecewise linear curve is continuous and strictly monotonic increasing according to the positive slopes (i.e a i > 0, or angles 0 • < atan(a i ) < 90 • ).

In both cases the parameter M → ∞ or p i = const, the limitation of slope

a i = l LDRmax -l LDRmin B i=1 δ i = l LDRmax -l LDRmin l B LDR (K B )-l 1 LDR (1)
= const. Therefore, the "s-shape" curve is linear.

III. SIMULATION RESULTS

This section provides the performance of the proposed tone mapped HDR image. The tone mapped image quality is measured with the TMQI (Tone-Mapped image Quality Index) metric [START_REF] Yeganeh | Objective Quality Assessment of Tonemapped Images[END_REF]. Simulations have been conducted under Matlab environnement using the HDR Toolbox ( [START_REF] Banterle | Advanced High Dynamic Range Imaging: Theory and Practice[END_REF]) with 274 test HDR images. For lack of space, we only present the results obtained with 8 HDR images ("Anturium", "Bottle Small", "Office", "Oxford Church", "Memorial", "Light", "WardFlowers" and "StreetLamp") with different dynamic range (or contrast ratio) from 8 f-stops to 19 f-stops.

The proposed approach is compared to: (i) NONSEP ENO-CA [13] ; (ii) SEP ENO-CA [12] with parameters α 1 = 0.3, α 2 = 0.7; (iii) Li TMO [8] with Haar multiscale; (iv) Fattal [11] using RBW method with parameters α = 0.8, β = 0.3, γ = 0.8; (v) Duan [9] using β = 0.5; (vi) TMOs in HDR Toolbox: Drago [7] , Reinhard [14] , Ward [15] , Durand [6] , Schlick [17] with the default parameters as given in the HDR Toolbox. The different parameters are chosen so as to give the best results in terms of TMQI metric in all methods.

Table I provides the TMQI metrics. The proposed TM approach namely "Proposed LJ NM" is deployed with B = 256, l LDRmax = 255, l LDRmin = 0, J = 1, ..., 5 and M = 1. Our approach is competitive to those developed in the literature. Combined to Table II, more the number of resolution levels increase as well as less norm space concerns, more the performance increase as the details around light can been clearly seen in Fig. 3. In addition, Fig. 7 and Fig. 11 show histograms of "Light" TM images with different levels and norm spaces. Fig. 13 give tone mapping curves with different norm spaces. The visual quality of our method in Fig. 8a is better than "Duan" method in Fig. 9b around the light.

Fig. 4 compares the visual quality of the "Church" tone mapped image using "Duan" method and our approach. The stained glass window at the church background presents a better contrast and details with our approach although the TQMI are identical.

Fig. 5 compares the visual quality of the "WardFlowers" tone mapped image using "Fattal" and our approach. Some details, on flowers and rocks, are lost on "Fattal" tone mapped image compared to our approach. Moreover, our tone mapped image if of better contrast. A similar result is provided by Fig. 6 where the HDR "StreetLamp" image has been mapped using "SEP ENO" method and our method. The brightness of our tone mapped is better.

The performance of our approach is confirmed on more than 274 test HDR images where the details and contrast are better represented than other competitive methods.

IV. CONCLUSION

This paper proposed a new HDR image TM approach able at the same time to extract the relevant details and enhance [7] 0.874 0.801 0.800 0.814 0.800 0.800 Reinhard [14] 0.778 0.807 0.826 0.789 0.791 0.794 Ward [15] 0.806 0.783 0.775 0.817 0.795 0.789 Durand [6] 0.811 0.892 0.825 0.929 0.814 0.800 Tumblin [16] 0.715 0.713 0.735 0.675 0.759 0.750 Schlick [17] 0.770 0.835 0.926 0.970 0.787 0.780 Duan [9] 0.964 0.916 0.955 0.986 0.935 0.916 Fattal [11] 0.889 0.928 0.943 0.889 0.927 0.971 Li [8] 0.964 0.954 0.854 0.877 0.834 0.888 SEP ENO [12] 0.896 0.934 0.943 0.895 0.932 0.970 NONSEP [13] 0

Fig. 1 :

 1 Fig.1: Piecewise linear curve modelization ("s-shape" curve).

Fig. 2 :Fig. 3 :

 23 Fig. 2: "Light" HDR test image (17.46 f-stops) -Left image: proposed optimized β (1 level Norm 1, TMQI=0.920); Right image: proposed β = 0.25 (TMQI=0.813).

Fig. 4 :

 4 Fig. 4: "Oxford Church" HDR test image (15.46 f-stops) -Left image: proposed (5 levels, TMQI=0.985); Right image: "Duan" (TMQI=0.986).

Fig. 5 :Fig. 7 :Fig. 8 :

 578 Fig. 5: "WardFlowers" HDR test image (14.01 f-stops) -Up image: proposed (5 levels, TMQI=0.930); Down image: "Fattal" (TMQI=0.875).

Fig. 9 :

 9 Fig. 9: Normalized histograms of "Light" tone mapped images (level 1) and "Duan".

Fig. 11 :Fig. 12 :

 1112 Fig. 11: Normalized histograms of "Light" tone mapped images (norm 1, 2) with optimized β.

Fig. 13 :

 13 Fig. 13: Tone-mapping curves with M norm space for "Light" HDR image with optimized β.

TABLE I :

 I Tone Mapped Image Quality Index (TMQI)

	TMOs	Anturium Bottle Office Church Memorial Light
	DR f-stops	8.73	16.03 16.29	15.46	18.38	17.46
	Drago

TABLE II :

 II Tone Mapped Image Quality Index (TMQI) for "Light" HDR test image (17.46 f-stops) for different norm spaces M = 1, 2, 3, 4, 6, 10 .920 0.891 0.876 0.867 0.856 0.847 Proposed L2 0.923 0.897 0.883 0.874 0.864 0.855 Proposed L3 0.931 0.907 0.894 0.885 0.876 0.887 Proposed L4 0.937 0.915 0.903 0.895 0.887 0.878 Proposed L5 0.942 0.922 0.911 0.903 0.895 0.888

	Norms	N1	N2	N3	N4	N6	N10
	Proposed L1 0


Ba Chien Thai

Anissa Mokraoui
email: anissa.mokraoui@univ-paris13.frmatei@lipn.univ-paris13.fr

Basarab Matei

L2ti

HDR Image Tone Mapping Approach based on Near Optimal Separable Adaptive Lifting Scheme

This paper proposes a Tone Mapping (TM) approach converting a High Dynamic Range (HDR) image into a Low Dynamic Range (LDR) image while preserving as much information of the HDR image as possible to ensure a good LDR image visual quality. This approach is based on a separable near optimal lifting scheme using an adaptive powerful prediction step. The latter relies on a linear weighted combination depending on the neighboring coefficients extracting then the relevant finest details in the HDR image at each resolution level. Moreover the approximation and detail coefficients are modified according to the entropy of each subband. The pixel's distribution of the coarse reconstructed LDR image is then adjusted according to a perceptual quantizer with respect to the human visual system using a piecewise linear function. Simulation results provide good results, both in terms of visual quality and TMQI metric, compared to existing competitive TM approaches.

I. INTRODUCTION

The objective of a High Dynamic Range (HDR) image Tone Mapping (TM) approach is to find a trade-off between the relevant information (e.g. details, contrast, brightness...) to be preserved or discarded in the image ensuring a good visual quality of the displayed image on Low Dynamic Range (LDR) devices that would be appreciated by observers.

A state of the art on HDR image TM approaches is fairly complete in [START_REF] Banterle | Advanced High Dynamic Range Imaging: Theory and Practice[END_REF], [START_REF] Dufaux | High Dynamic Range Video 1st Edition : From Acquisition, to Display and Applications[END_REF] and [START_REF] Reinhard | High Dynamic Range Imaging 2nd Edition: Acquisition, Display, and Image-Based Lighting[END_REF] where a classification is provided. Among the developed TM strategies, this paper quickly reviews those that caught our attention because of their performance. In [START_REF] Durand | Fast Bilateral Filtering for The Display of High-dynamic-range Images[END_REF], the TM approach reduces the HDR contrast while preserving the image details. This work uses an edge-preserving bilateral filter to decompose the HDR image into two layers: a base layer encoding large-scale variations and a detail one. Contrast is then reduced only in the first layer while the details are kept unchanged. An adaptive logarithmic mapping method of luminance values is presented in [START_REF] Drago | Adaptive Logarithmic Mapping for Displaying High Contrast Scenes[END_REF]. It concerns the adjustment of the logarithmic basis depending on the radiance of the pixels. In [START_REF] Li | Compressing and Companding High Dynamic Range Images with Subband Architectures[END_REF], a subband architecture related on an oversampled Haar pyramid representation is proposed. Subband coefficients are re-scaled according to a gain control function reducing the high frequency magnitudes and boosting low ones. In [START_REF] Duan | Tone-mapping High Dynamic Range Images by Novel Histogram Adjustment[END_REF], a TM optimization approach using a histogram adjustment between linear mapping and the equalized histogram mapping is developed. A modification of this approach is made where revisited histogram equalization approaches are discussed in [START_REF] Husseis | Revisited Histogram Equalization as HDR Images Tone Mapping Operators[END_REF]. The latter considers both histogram equalization and human sensitivity to the light function.

In [START_REF] Fattal | Edge-Avoiding Wavelets and their Applications[END_REF], a second generation of wavelets based on the edge content of the image avoiding having pixels from both sides of an edge is proposed. In [START_REF] Thai | Performance Evaluation of High Dynamic Range Image Tone Mapping Operators Based on Separable Non-linear Multiresolution Families[END_REF], a separable non-linear multiresolution approach based on essentially non-oscillatory interpolation strategy has been investigated. These approaches take into account the singularities (e.g. edge points) in their mathematical models preserving then the structural information of the HDR images. In [START_REF] Thai | Image Tone Mapping Approach Using Essentially Non-Oscillatory Bi-quadratic Interpolations Combined with a Weighting Coefficients Strategy[END_REF], a non-separable non-linear multiresolution approach is proposed. The results provided in [START_REF] Fattal | Edge-Avoiding Wavelets and their Applications[END_REF], [START_REF] Thai | Performance Evaluation of High Dynamic Range Image Tone Mapping Operators Based on Separable Non-linear Multiresolution Families[END_REF] and [START_REF] Thai | Image Tone Mapping Approach Using Essentially Non-Oscillatory Bi-quadratic Interpolations Combined with a Weighting Coefficients Strategy[END_REF] show that the decomposition of the HDR image on different resolution levels would seem to be a good strategy. However, the choice of the decomposition filters is extremely important and decisive in the extractive power of the approximation and detail information. To do so, this paper proposes a separable near optimal lifting scheme. The "Predict" operation is performed by a new adaptive prediction operation to extract the finest detail coefficients to highlight the sharp transitions in the given HDR image.

The proposed HDR image TM approach has two main goals namely the preservation of the details and the adjustment of the contrast in accordance with the LDR display devices. It is composed of four stages. The first one judiciously decomposes the HDR image into different resolution levels (section II-A). The second one concerns the weighting strategy of the detail and approximation coefficients (section II-B). The third one reconstructs the coarse LDR image (section II-C). Finally the fourth one adjusts the contrast according to a perceptual linear quantizer (section II-D). Section III discusses the simulation results. Section IV concludes the paper.

II. PROPOSED HDR IMAGE TONE MAPPING APPROACH

This section concerns the first, second, third and fourth stages of the proposed HDR image TM mapping approach (i.e. decomposition, weighting, reconstruction and contrast adjustment stages).

Before describing the approach, introduce some notations. The original HDR image, at the finest resolution level J, is assumed to be of size N J × M J . The index j refers to the resolution level (with j = 0, ..., J -1). Denote l HDR the HDR image luminance. In the rest of this paper, the HDR image luminance is considered in the logarithm domain since it is well adapted to the human visual system. It is denoted I J and defined as follows I J = {I J (x n , y m ) = log 10 (l HDR (x n , y m )) for 1 ≤ n ≤ N J and 1 ≤ m ≤ M J where I J (x n , y m ) is the HDR logarithm luminance value of the pixel located at position (x n , y m ) on the image.

A. First stage: HDR image decomposition

The proposed HDR image decomposition is performed according to the forward process of a separable near optimal cell-average lifting scheme. This strategy is motivated by the fact that the relevant details are accurately predicted since the coefficients of the filters are adapted locally to the image to be processed.

The decomposition consists to go from the finest resolution level J to the coarse resolution level 0. At a given resolution level j -1 (with 0 ≤ j ≤ J -1), the algorithm deals with the approximation coefficients denoted I j (x n , y k ) (with 1 ≤ n ≤ N j and 1 ≤ k ≤ M j ) computed at resolution level j. For a given n belonging to [1, N j ], the algorithm starts with splitting the horizontal 1D-signal, i.e. I j (x n , y k ) for 1 ≤ k ≤ M j , into a set of odd and even indexes as follows:

{I j (x n , y k ) with 1 ≤ k ≤ M j } := {I j (x n , y 2k-1 ), I j (x n , y 2k ) with 1 ≤ k ≤ M j /2 }. (1)
This process is then repeated for all n. Based on this split, the approximation coefficient located at position (x n , y k ), denoted V j-1 (x n , y k ), is computed on a Cell-Average (CA) scheme. For a given n, it is expressed as follows:

V j-1 (x n , y k ) = I j (x n , y 2k-1 ) + I j (x n , y 2k ) 2 for 1 ≤ k ≤ M j /2. ( 2
)
This process is then repeated for all n.

The detail coefficient, denoted d j-1 (x n , y 2k-1 ), is computed at odd indexes and is provided below for a given n:

d j-1 (x n , y 2k-1 ) = Îj (x n , y 2k-1 ) -I j (x n , y 2k-1 ) for 1 ≤ k ≤ M j /2, (3)
where Îj (x n , y 2k-1 ) is the predicted logarithm luminance value at odd position (x n , y 2k-1 ) and resolution level j. This predicted value is expressed as a linear weighted combination of the neighboring approximation coefficients at resolution level j -1:

Îj (x n , y 2k-1 ) = 2 i=0 u j-1 i (x n ) • V j-1 (x n , y k+i-1 ) for 1 ≤ k ≤ M j /2, (4)
where the weights u j-1 i (x n ) must preserve the initial 1Dsignal average which results in satisfying this condition:

2 i=0 u j-1 i (x n ) = 1. (5)
In what follows, the weights are rather denoted in a vector form

u j-1 = (u j-1 0 (x n ), u j-1 1 (x n ), u j-1 2 (x n ))
t to lighten the writing. These weights are deduced so that the Mean Squared Error (MSE) between Îj (x n , y 2k-1 ) and I j (x n , y 2k-1 ) is minimized:

E = arg min u j-1 i Îj (x n , y 2k-1 ) -I j (x n , y 2k-1 ) 2 2 for 1 ≤ k ≤ M j /2. (6)
This results in solving the following equation:

Γ j-1 • u j-1 = r j-1 , (7)
where u j-1 is the weight vector. r j-1 is the crosscorrelation vector r j-1 = (r j-1 (0), r j-1 (1), r j-1 (2)) t where r j-1 (i) represents the cross-correlation function between V j-1 (x n , y k+i-1 ) and I j (x n , y 2k-1 ) for 1 ≤ k ≤ M j /2. Γ j-1 is the autocorrelation matrix defined as:

Γ j-1 =   R j-1 (0) R j-1 (1) R j-1 (2) R j-1 (-1) R j-1 (0) R j-1 (1) R j-1 (-2) R j-1 (-1) R j-1 (0)   , (8)
where R j-1 (i) is the autocorrelation of V j-1 (x n , y k ). These weights, associated to the row (x n ), are then deduced so that the partial derivatives of the MSE (given by equation ( 6)) with respect to u j-1 i (i.e. i = 0, 1, 2) are equal to zero:

u j-1 (x n ) = (Γ j-1 ) -1 • r j-1 . (9)
The weight vectors are then stored to be used in the adaptive lifting scheme backward process to reconstruct the decomposed 1D-signal. For sake of convenience, the approximation and details coefficients are organized as follows:

{W j (x k , y m ) for 1 ≤ k ≤ N j and ∀m ∈ [1, M j /2]} := {V j-1 (x n , y k ) for 1 ≤ k ≤ M j /2 and ∀n ∈ [1, N j ]}, {U j (x k , y m ) for 1 ≤ k ≤ N j , ∀m ∈ [1, M j /2]} := {d j-1 (x n , y 2k-1 ) for 1 ≤ k ≤ M j /2, ∀n ∈ [1, N j ]}.
(10) The split, approximation and detail operations are applied on W j (x k , y m ) and U j (x k , y m ) for a given m ∈ [1, M j /2] (i.e. on the vertical direction). Note that the approximation step requires the prediction of Ŵ j (x 2k-1 , y m ) (respectively Û j (x 2k-1 , y m )) based on a set of weights v j-1 i (respectively w j-1 i ). These weights need to be stored for the backward process lifting scheme to reconstruct the decomposed image. Finally, the approximation resolution level

I j (x n , y m ) is divided into 4 blocks I j := (I j-1 , d j-1 HL , d j-1 LH , d j-1 HH ).
The decomposition is thus iterated on I j-1 until j = 0.

The finest HDR image I J is then represented by 3J + 1 resolution levels

I J := (I 0 , d 0 HL , d 0 LH , d 0 HH , ..., d J-1 HL , d J-1 LH , d J-1 HH ).

B. Second stage: Weighting strategy of the coefficients

To reduce the dynamic range of the HDR image, the TM approach proposes to weight the approximation and detail coefficients in an appropriate way before performing the adaptive lifting scheme backward process described in section II-C.

Denote N l the number of resolution levels equal to J; E a the entropy of the approximation coefficients at the coarsest resolution level (i.e. j = 0); and E j d the entropy of the detail coefficients at resolution level j (i.e. d j HL , d j LH and d j HH ) getting therefore

N l + 1 entropies (E a , E 0 d , E 1 d , ..., E j d ,..., E N l -1 d
). From these entropies, positive weights smaller than one are deduced as follows:

     α a = N l -1 i=0 E i d Ea+ N l -1 i=0 E i d for j = 0 α j d = Ea+ N l -1 i=0, =j E i d Ea+ N l -1 i=0 E i d for j = 0, ..., N l -1 (11)
α a (respectively α j d ) is the weight associated to the approximation (respectively detail) coefficients at resolution level j = 0 (respectively j).

The coefficients of the four coarsest resolution levels are first modified according to:

I 0 = α a × I 0 , d 0 HL = α 0 d × d 0 HL , d 0 LH = α 0 d × d 0 LH , d 0 HH = α 0 d × d 0 HH . (12)
The approximation subband, denoted I 1 , is then reconstructed (see section II-C) and the number of levels is reduced to N l -1. The N l entropies (associated to 3N l -2 subbands) are calculated again to update the weights α a , α 1 d (equation ( 11) with N l = N l -1). After that, these weights are applied on the coefficients I 1 , d 1 HL , d 1 HL , d 1 HL to build I 2 as explained above. This process is iterated until N l = 0 to reconstruct the coarse tone mapped HDR image denoted I J LDR , called coarse LDR image.

C. Third stage: Reconstruction of the coarse LDR image

The reconstruction stage is carried out inversely to the decomposition stage. Assume that the adaptive lifting scheme backward algorithm processed all resolution levels until j -1. The next step consists to recover the approximation coefficients I j of size N j × M j using the four weighted blocks I j-1 , d j-1 LH , d j-1 HL and d j-1 HH . The algorithm first deals with the coefficients in a vertical direction using I j-1 and d

j-1 LH

(denoted d j-1 ). At a given m, the approximation coefficient at odd position, denoted W j (x 2k-1 , y m ), is deduced:

W j (x 2k-1 , y m ) = Ŵ j (x 2k-1 , y m ) -d j-1 (x k , y m ) for 1 ≤ k ≤ N j /2, (13)
where Ŵ j (x 2k-1 , y m ) is the predicted coefficient, at odd position, deduced from the weighted combination v j-1 i (with i = 0, 1, 2) of the neighboring approximation coefficients I j-1 (x 2k-1 , y m ) as:

Ŵ j (x 2k-1 , y m ) = 2 i=0 v j-1 i (y m ) • I j-1 (x k+i-1 , y m ) for 1 ≤ k ≤ N j /2. (14)
The weights v j-1 i (y m ) are those computed and stored in the decomposition process. The approximation coefficient W j (x 2k , y m ), at even position, is deduced in a CA scheme:

W j (x 2k , y m ) = 2I j-1 (x k , y m ) -W j (x 2k-1 , y m ) for 1 ≤ k ≤ N j /2. (15)
The odd and even approximation coefficients are then merged to built 1D-signal:

{W j (x k , y m ) for 1 ≤ k ≤ N j } := {W j-1 (x 2k-1 , y m ), W j-1 (x 2k , y m ) for 1 ≤ k ≤ N j /2}. ( 16
) This process is repeated for all m to reconstruct W j of size N j × M j /2. These same steps are applied on d j-1 HL and d j-1 HH to generate the new block U j of size N j × M j /2. As in the decomposition strategy, W j and U j are renamed V j-1 and d j-1 . The same steps, as described above, are then performed but according to a horizontal direction. The reconstruction is iterated to finally build the image called coarse LDR image which is denoted I J LDR .

D. Fourth stage: Piecewise linear perceptual quantizer

This stage proposes to adjust locally the distribution of the coarse LDR image logarithm luminance I J LDR according to the HVS to enhance the contrast using a piecewise linear function. This strategy is inspired from [START_REF] Mai | Visually-favorable tone-mapping with high compression performance[END_REF] which has been developed for compression purpose. However, modifications are made mainly to avoid the problem of empty bins of equal size.

To do so, the I J LDR values are first sorted and classified into equal B bins defined by cutting points denoted c i uLDR (with 1 ≤ i ≤ B). A non-uniform histogram equalization is also performed. c i nuLDR (with 1 ≤ i ≤ B) cutting points, defining the bounds of the non-uniform consecutive B bins, are deduced. The lower bound (cutting point) of each bin is then adjusted as follows:

l i LDR (1) = c i uLDR + β(c i nuLDR -c i uLDR ), (17)
where β is a positive parameter smaller than 1. Therefore the I J LDR values are classified into non-uniform B bins as follows:

I J LDR = { l LDR (k) for k = 1, ..., N J × M J } = {[ l 1 LDR (1), ..., l 1 LDR (K 1 )], ..., [ l i LDR (1), ..., l i LDR (K i )], ..., [ l B LDR (1), ..., l B LDR (K B )]}, (18)
depending on the quantization level set, where K i is the number of values in the i -th bin (i.e. 1 ≤ i ≤ B; K i > 0) and satisfying the following relation

B i=1 K i = N J × M J .
The "s-shape" TM perceptual curve, as discussed in [START_REF] Dowling | The Retina: An Approachable Part of the Brain[END_REF] and [START_REF] Geisler | Effects of Bleaching and Backgrounds on the Flash Response of the Cone System[END_REF], is modelled by a piecewise linear curve on each bin (see Fig. 1). Consider the i -th bin, defined by

[ l i LDR (1), ..., l i LDR (K i )],
the coarse LDR values are then modeled as follows:

li LDR (k) = a i l i LDR (k) + b i with k ∈ [1, K i ], (19)
where a i (with a i = 0) and b i are two unknown parameters depending on the i -th bin. This equation, after some mathematical manipulations, can be rewritten as follows:

l i LDR (k) = li LDR (k) -li LDR (1) a i + l i LDR (1), (20)
where the unknown parameter a i is deduced so that the MSE between the coarse LDR value and its quantized version, denoted l i LDR (k) (i.e. when li LDR (k) has been supported a ceiling process according to a scalar quantization), is minimized in the i -th bin: arg min

a i l i LDR (k) -l i LDR (k) 2 2 . ( 21
)
This equation is simplified as follows:

arg min

a i Ki k=1 l i LDR (k)-li LDR (k) a i 2 • p i , (22)
where p i is the probability of the l i LDR (k) value in the i -th bin given by p i = Ki B i=1 Ki . Extending equation ( 22) to all bins involves the computation of a Global MSE (GMSE) deduced as follows:

GM SE = B i=1 Ki k=1 l i LDR (k)-li LDR (k) a i 2 • p i . (23)
The variance of (l i LDR (k) -li LDR (k)) on each bin is assumed to be equal and is denoted ξ. Equation ( 23) is then simplified and becomes:

GM SE = B i=1 pi (a i ) 2 • ξ. (24)
Denote l LDRmax (respectively l LDRmin ) the maximum (respectively minimum) LDR luminance value. Introduce δ i as the difference between coarse LDR luminance in two consecutive bins:

δ i = li+1 LDR (1) -li LDR (1). (25)
A constraint related to the limit sum of the projected heights equal to the entire LDR range results in:

B i=1 a i • δ i = l LDRmax -l LDRmin . (26)
Therefore the optimization problem is written as follows:

arg min

a i B i=1 p i (a i ) 2 •ξ, s.t B i=1 a i •δ i = l LDRmax -l LDRmin .
(27) This problem is solved analytically using the Lagrangian function. After some mathematical manipulations, the slope a i is deduced:

a i = (l LDRmax -l LDRmin ) • (p i ) 1/3 B i=1 δ i • (p i ) 1/3 . ( 28
)
Hence the unknown parameter b i is calculated (i.e. b i = li LDR (1)-a i × l i LDR (1)) and LDR mapped values are deduced according to equation [START_REF] Mai | Visually-favorable tone-mapping with high compression performance[END_REF]. The global piecewise linear curve is continuous and strictly monotonic increasing according to the positive slopes (i.e a i > 0, or angles 0 • < atan(a i ) < 90 • ).

III. SIMULATION RESULTS

This section provides the performance of the proposed tone mapped HDR image. The tone mapped image quality is measured with the TMQI (Tone-Mapped image Quality Index) metric [START_REF] Yeganeh | Objective Quality Assessment of Tonemapped Images[END_REF]. Simulations have been conducted under Matlab environnement using the HDR Toolbox ( [START_REF] Banterle | Advanced High Dynamic Range Imaging: Theory and Practice[END_REF]) with 274 test HDR images. For lack of space, we only present the results obtained with 8 HDR images ("Anturium", "Bottle Small", "Office", "Oxford Church", "Memorial", "Light", "WardFlowers" and "StreetLamp") with different dynamic range (or contrast ratio) from 8 f-stops to 19 f-stops.

The proposed approach is compared to: (i) NONSEP ENO-CA [13] ; (ii) SEP ENO-CA [12] with parameters α 1 = 0.3, α 2 = 0.7; (iii) Li TMO [8] with Haar multiscale; (iv) Fattal [11] using RBW method with parameters α = 0.8, β = 0.3, γ = 0.8; (v) Duan [9] using β = 0.5; (vi) TMOs in HDR Toolbox: Drago [7] , Reinhard [14] , Ward [15] , Durand [6] , Schlick [17] with the default parameters as given in the HDR Toolbox. The different parameters are chosen so as to give the best results in terms of TMQI metric in all methods.

Table I provides the TMQI metrics. The proposed TM approach namely "Proposed LJ" is deployed with B = 256, β = 0.25, l LDRmax = 255, l LDRmin = 0 and J = 1, ..., 5. Our approach is competitive to those developed in the literature. More the number of resolution levels increase, more the performance increase.

Fig. 2 compares the visual quality of the "Church" tone mapped image using "Duan" method and our approach. The stained glass window at the church background presents a better contrast and details with our approach although the TQMI are identical. Fig. 3 compares the "Memorial" tone mapped image using "Duan" and "Fattal" methods and our approach. The details on tills (see Fig. 3) and rosette (see Fig. 4) are better rendered by our approach.

Fig. 5 compares the visual quality of the "WardFlowers" tone mapped image using "Fattal" and our approach. Some details, on flowers and rocks, are lost on "Fattal" tone mapped image compared to our approach. Moreover, our tone mapped image if of better contrast. A similar result is provided by Fig. 6 where the HDR "StreetLamp" image has been mapped using "SEP ENO" method and our method. The brightness of our tone mapped is better.

The performance of our approach is confirmed on more than 274 test HDR images where the details and contrast are better represented than other competitive methods.

IV. CONCLUSION

This paper proposed a new HDR image TM approach able at the same time to extract the relevant details and enhance the contrast of LDR images. This is essentially related to : (i) the forward process of the near optimal local adaptive cell average lifting scheme where the filter coefficients are locally adapted to the content; (ii) the weighting operation depending on the information of each subband; (iii) the adjustment of the coarse LDR image luminance distribution according to the perceptual piecewise linear function. Simulation results confirm the relevance of the proposed approach both in terms of the TMQI metric and the visual quality of the displayed image. [7] 0.874 0.801 0.800 0.814 0.800 0.800 Reinhard [14] 0.778 0.807 0.826 0.789 0.791 0.794 Ward [15] 0.806 0.783 0.775 0.817 0.795 0.789 Durand [6] 0.811 0.892 0.825 0.929 0.814 0.800 Tumblin [16] 0.715 0.713 0.735 0.675 0.759 0.750 Schlick [17] 0.770 0.835 0.926 0.970 0.787 0.780 Duan [9] 0.964 0.916 0.955 0.986 0.935 0.969 Fattal [11] 0.889 0.928 0.943 0.889 0.927 0.971 Li [8] 0.964 0.954 0.854 0.877 0.834 0.888 SEP ENO [12] 0.896 0.934 0.943 0.895 0.932 0.970 NONSEP [13] 0

Fig. 1 :

 1 Fig. 1: Piecewise linear curve modelization ("s-shape" curve).

Fig. 2 :

 2 Fig. 2: "Oxford Church" HDR test image (15.46 f-stops) -Up image: proposed (5 levels, TMQI=0.985); Down image: "Duan" (TMQI=0.986).

Fig. 3 :

 3 Fig. 3: "Memorial" HDR test image (18.38 f-stops) -Left image: proposed (5 levels, TMQI=0.951); Middle image: "Duan" (TMQI=0.935); Right image: "Fattal" (TMQI=0.927).

Fig. 4 :

 4 Fig. 4: LDR luminance "Rosette" zoom -Left image: proposed (5 levels, TMQI=0.951); Middle image: "Duan" (TMQI=0.935); Right image: "Fattal" (TMQI=0.927).

Fig. 5 :

 5 Fig. 5: "WardFlowers" HDR test image (14.01 f-stops) -Up image: proposed (5 levels, TMQI=0.930); Down image: "Fattal" (TMQI=0.875).

Fig. 6 :

 6 Fig. 6: "StreetLamp" HDR test image (13.83 f-stops) -Up image: proposed (5 levels, TMQI=0.911); Down image: "SEP ENO" (TMQI=0.855).

TABLE I :

 I Tone Mapped Image Quality Index (TMQI)

	TMOs	Anturium Bottle Office Church Memorial Light
	DR f-stops	8.73	16.03 16.29	15.46	18.38	17.46
	Drago


Ba Chien
email: tbchien@dut.udn.vn

Anissa Mokraoui
email: anissa.mokraoui@univ-paris13.fr

Performance Analysis of the Revisited Tone Mapped Quality Index for Tone Mapped HDR Images Evaluation

Keywords: High Dynamic Range image, Image Tone Mapping Operators, Human Visual System, Mean Opinion Score, Tone Mapped Quality Index

This paper discusses the appropriate choice of the Tone Mapped Quality Index (TMQI) parameters. The TMQI metric has been developed to evaluate and compare the visual quality of the Tone Mapped (TM) High Dynamic Range (HDR) images. However, the parameters of this metric have been adjusted on a specific and reduced set of Tone Mapping Operators (TMO) and a restricted base of HDR images. The visual rendering evaluation of the TM images is then no longer consistent with the metric when new TMO (excluded from the initial training process) are used. To better adjust these parameters, discussions and experiments have been conducted on a large extension of the initial TMOs set and tone mapped HDR images dataset. The new parameters show a strong correlation between the revisited metric and the Mean Opinion Score (MOS).

I. INTRODUCTION

A High Dynamic Range (HDR) image contains both underexposed (i.e. very dark) and over-exposed (i.e. bright) areas in accordance with the real-world scene. Typical examples of HDR scenes are sunny outdoor scenes or views of an indoor room with objects visible outside the window. The dynamic ranges of several natural scenes have about 1 : 160 average contrast ratio. Outdoor scenes usually have a larger dynamic range, which can reach a contrast ratio of three orders of magnitude 1 : 1000 or sometimes even more. However, these HDR images cannot be visualized on standard Low Dynamic Range (LDR) display devices since their dynamic range is smaller than that of HDR images. The currently used display technologies are Cathode Ray Tube (CRT) or Liquid Crystal Displays (LCD), they are typically discretized at 8-bits and rarely at 10-bits per color channel. This means that colors shades are limited to 255. Moreover HDR display devices currently remain too expensive. Therefore many image Tone Mapping Operators (TMOs) have been developed. Their main concerns are to reduce the dynamic range (contrast, color gamut, details...) of those HDR images to the dynamic range of LDR display device while preserving, as much as possible, the appearance of the captured scene in terms of contrast and the overall impression of brightness and colors.

Many HDR image TMOs have been developed to convert HDR images into LDR images. However to measure the performance of these approaches, objective metrics are required not only to compare the visual rendering quality but also to provide a better understanding of the relationship between an operator and the image attributes. Indeed this can help in the development of TMOs more suited to a particular application such as an improved method for highlighting those parts of an image that are perceptually important to the Human Visual System (HVS).

The assessment of the Tone Mapped (TM) image is restricted either to no reference criteria, which were mostly trained and tested in different context, or to the specifically adjusted full reference metrics. All of the full reference image quality metrics assume that the dynamic range of the original and processed image are the same. However, in the case of TM HDR images, the dynamic range between the two versions is different (i.e. I HDR is HDR image while I LDR is LDR). That is why the development of objective metrics is not easy.

The main purpose of this paper is to choose one objective metric among many others that is best suited to our concerns, namely the visual quality evaluation of the TM HDR images with no HDR display device. To do so, this paper proceeds as follows. Section II introduces the well known objective metrics to assess the visual quality of TM HDR images. Section III analyzes the performance of these objective metrics when a large set of TMOs are applied to an extensive HDR image database. The metrics are compared to subjective assessments made by observers. The analysis of simulation results shows that none of the objective metrics is highly correlated to the Mean Opinion Score (MOS). However, the Tone Mapped Quality Index (TMQI) would seem to better meet our expectations although it is designed in a different context of experimentation of ours. It is then revised in section III-C with a large set of TMOs and an extensive HDR image database. Section IV confirms the performance of the revisited TMQI metric. Section V concludes this paper.

II. AN OVERVIEW OF OBJECTIVE VISUAL QUALITY

METRICS FOR TONE MAPPED HDR IMAGES This section introduces some objective metrics often used to assess the visual quality of the tone mapped HDR images.

A. Dynamic Range Independent Metric

The first metric, called Dynamic Range Independent Metric (DRIM) [START_REF] Aydin | Dynamic range independent image quality assessment[END_REF], exploits the detection model from HDR-VDP [START_REF] Mantiuk | Predicting Visible Differences in High Dynamic Range Images -Model and its Calibration[END_REF] which is calibrated on ModelFest dataset [START_REF] Watson | Visual Detection of Spatial Contrast Patterns: Evaluation of Five Simple Models In[END_REF] to indicate the regions containing visible contrast in the HDR image and its TM version. The DRIM metric proposes three distortion maps showing the regions where the contrast is either lost (i.e. contrast change is perceivable in HDR but imperceivable in LDR image), amplified (the opposite case), or reversed (the polarity is changed -mostly caused by halo artifacts) by tone mapping process. This DRIM metric inherits high computational requirements and the necessity to specify viewing conditions and display parameters. Note that this metric is designed for visualization of the perceived distortions regarding the image contrast. As a result, it does not provide a single value for quality assessment. For this reason, it will not be considered in the performance analysis section.

B. Tone Mapped image Quality Index

Another HDR image tone mapping metric has been introduced by Yeganeh and Wang known as Tone Mapped image Quality Index (called TMQI1) [START_REF] Yeganeh | Objective Quality Assessment of Tone Mapped Images[END_REF]. It is composed of two parts: the Structural Fidelity 1 (SF 1) part and Statistical Naturalness 1 (SN 1) part described below.

1) Structural Fidelity: SF 1 is a modified version of the Multi Scale Structural Similarity Index (MS-SSIM) [START_REF] Wang | Multiscale structural similarity for image quality assessment[END_REF]. For a patch u, it is given by:

SF 1 IHDR(u), ILDR(u) = 2σ IHDR(u) • σ ILDR(u) + C1 σ 2 IHDR(u) + σ 2 ILDR(u) + C2 × σ IHDR(u), ILDR(u) + C2 σ IHDR(u) • σ ILDR(u) + C2 , (1)
where σ (.) is the standard deviation of the patch after the nonlinear mapping; σ I HDR (u) , σ I LDR (u) and σ I HDR (u), I LDR (u) are the local standard deviations and cross correlation between the two corresponding patches in HDR and LDR images, respectively; and C 1 , C 1 are positive constants that make sure the equation not divided by zero.

Note that the luminance component is missing, compared the SSIM definition, but the structural element (i.e. the second fraction in equation ( 1)) remains the same.

The mapping is defined as follows:

σ (σ) = 1 √ 2πθ σ σ -∞ exp[- (x -τ σ ) 2 2θ 2 σ ], (2)
where

θ σ (f ) = τ σ (f ) k , (3)
with f being a spatial frequency and k representing a constant obtained from Croziers law [START_REF] Crozier | On the Variability of Critical Illumination for Flicker Fusion and Intensity Discrimination[END_REF], typically ranging from 2.3 to 4. The authors proposed to set k = 3. The threshold for the signals standard deviation is given by:

τ σ (f ) = μ √ 2 × λ × CSF (f ) , ( 4
)
where μ is the mean intensity value (set to 128 by the authors) and λ is a constant used to fit the physiological data. TMQI1 uses the CSF as introduced by Mannos and Sakrison [START_REF] Mannos | The Effects of a Visual Fidelity Criterion of the Encoding of Images[END_REF] and fits to the data measured by Kelly [START_REF] Kelly | Effects of Sharp Edges on the Visibility of Sinusoidal Gratings[END_REF]. The map of SF 1 is averaged on each scale and the final SF 1 index is obtained in the same way as in the case of MS-SSIM.

2) Statistical Naturalness: The second part of the metric, i.e. SN 1, does not need a reference. It is based on the assumption that naturalness of an image can be modeled by the probability distribution of the brightness and contrast (means and standard deviations) in natural gray-scale images. The means and standard deviations are modelled by a Gaussian (N ) and a Beta (B) distribution respectively. The distributions are then given by N (115.94, 27.99) and B(4.4, 10.1). Assuming that brightness and contrast are mutually independent, the probability that the image is natural is then expressed as:

SN 1(I LDR ) = 1 K(I LDR ) × P df N (I LDR ) × P df B (I LDR ), (5)
where K is a factor used for the normalization, thus:

K(I LDR ) = max{P df N (I LDR ), P df B (I LDR )}. (6)
The TMQI1 expression is given by the combination of the two measures defined as:

T M QI1(I HDR , I LDR ) = a × SF 1 α (I HDR , I LDR ) +(1 -a) × SN 1 β (I LDR ), (7)
where the values of the parameters have been deduce as follows: a = 0.8012, α = 0.3046 and β = 0.7088 in an empirical way using 8 TMOs ("Reinhard" [START_REF] Reinhard | Dynamic Range Reduction Inspired by Photoreceptor Physiology[END_REF], "Drago" [START_REF] Drago | Adaptive Logarithmic Mapping for Displaying High Contrast Scenes[END_REF], "Durand" [START_REF] Durand | Fast Bilateral Filtering for The Display of High-dynamic-range Images[END_REF], "Mantiuk" [START_REF] Mantiuk | A Perceptual Framework for Contrast Processing of High Dynamic Range Images[END_REF], "Pattanaik" [START_REF] Sumanta | Time-dependent Visual Adaptation for Fast Realistic Image Display[END_REF], and three remained TMOs in Adobe Photoshop namely "Exposure and Gamma", "Equalize Histogram" and "Local Adaptation" respectively) and a reduced database with 15 HDR test images whose dynamic range from 8 f-stops to 18 f-stops.

C. Modified Tone Mapped image Quality Index

Based on the TMQI1, Ma et al. [START_REF] Ma | High Dynamic Range Image Compression by Optimizing Tone Mapped Image Quality Index[END_REF] revised both of the terms and proposed a new version of the TMQI1 called TMQI2. Namely, the contrast visibility model for HDR images has been adapted to the local luminance. The estimate of the contrast in the HDR reference is therefore computed as the standard deviation in a patch divided by the local mean. However, the second term is modified much more severely.

D. Feature Similarity Index for TM images

The last full-reference metric is the Feature Similarity Index for TM images (called FSITM) [START_REF] Nafchi | FSITM: A Feature Similarity Index For Tone-Mapped Images[END_REF]. It is based on the phase congruency features to calculate the difference between original and TM version of the the HDR image. More specifically, it computes the Locally Weighted Mean Phase Angle (LWMPA) to deduce the phase congruency. The main advantage of this feature is its robustness against noise.

III. PERFORMANCE ANALYSIS OF THE REVISITED VISUAL QUALITY OBJECTIVE METRIC

This section first compares the visual quality objective metric of the TM HDR images, described in the previous section, to the MOS. Discussions are provided in section III-B. Then, to improve the evaluation performance in accordance with the discussions provided in section III-B, section III-C revisits the selected objective metric so that to be consistent with the MOS.

A. Evaluation context

Note that during the visual evaluation the observers have no HDR display device (no-reference). The TM HDR images are displayed on the ColorEdge CG242W Color LCD monitor available in the lab with a good calibration in a random order during the testing period.

An observer assigns to each TM HDR image a mark ranging from 0 until 5: 5 for excellent, 4 very good, 3 good, 2 accepted, 1 unsatisfactory and 0 failed. The MOS table is obtained upon 26 observers. The chosen criteria is based on ability to recover details with natural ways.

Pearsons Linear Correlation Coefficient (PLCC) is selected to measure the relation between the objective visual quality metric of the TM image (i.e. Metric) and the MOS as follows:

P LCC = N n=1 M m=1 (MOS(n, m) -MOS) • (Metric(n, m) -Metric) N n=1 M m=1 (MOS(n, m) -MOS) 2 × 1 N n=1 M m=1 (Metric(n, m) -Metric) 2 , ( 8
)
where MOS and Metric are the mean values of two matrices MOS and Metric, respectively. P LCC ranges from -1 and +1.

Fifteen TMOs are used for the training such as "Drago" [START_REF] Drago | Adaptive Logarithmic Mapping for Displaying High Contrast Scenes[END_REF], "Reinhard" [START_REF] Reinhard | Dynamic Range Reduction Inspired by Photoreceptor Physiology[END_REF], "Ward" [START_REF] Ward | A Visibility Matching Tone Reproduction Operator for High Dynamic Range Scenes[END_REF], "Durand" [START_REF] Durand | Fast Bilateral Filtering for The Display of High-dynamic-range Images[END_REF], "Tumblin" [START_REF] Tumblin | Tone Reproduction for Realistic Images[END_REF], "Schlick" [START_REF] Schlick | Quantization Techniques for Visualization of High Dynamic Range Pictures[END_REF], "Duan" [START_REF] Duan | Tone-mapping High Dynamic Range Images by Novel Histogram Adjustment[END_REF], "Fattal" WRB [START_REF] Fattal | Edge-Avoiding Wavelets and their Applications[END_REF], "Li" [START_REF] Li | Compressing and Companding High Dynamic Range Images with Subband Architectures[END_REF], "Husseis" [START_REF] Husseis | Revisited Histogram Equalization as HDR Images Tone Mapping Operators[END_REF], "NUHA" [START_REF] Thai | Piecewise linear perceptual quantizer as a non-uniform histogram equalization adjustment for contrast enhancement of tone mapped HDR images[END_REF], "SEPENO" [START_REF] Thai | Performance Evaluation of High Dynamic Range Image Tone Mapping Operators Based on Separable Non-linear Multiresolution Families[END_REF], "NONSEPENO" [START_REF] Thai | Image Tone Mapping Approach Using Essentially Non-Oscillatory Bi-quadratic Interpolations Combined with a Weighting Coefficients Strategy[END_REF], "CEDP Lin" [START_REF] Thai | Contrast Enhancement and Details Preservation of Tone Mapped High Dynamic Range Images[END_REF] with β = 0.25 and "CEDP Opt" [START_REF] Thai | Contrast Enhancement and Details Preservation of Tone Mapped High Dynamic Range Images[END_REF] with adaptive β. The training TM HDR image database is based on 24 HDR images with different dynamic range (or contrast ratio) from 7 f-stops to 29 f-stops, namely "Lausanne1", "CraterLake1", "Shasta2", "Synagogue", "Anturium", "BowRiver", "Bridges", "Stairway1", "ArchRock", "DollDoll", "ClockBuilding", "OxfordChurch", "Bot-tlesSmall", "Montreal", "SmallOffice", "Light", "BridgeStudios2", "Memorial", "ClaridgeHotel", "Mistaya1", "BrookHouse", "Peace-Rocks", "GGpark2" and "AtriumNight". Therefore the training database contains 360 TM HDR.

B. Comparison of the Pearson Linear Correlation Coefficient using default parameters of the quality metrics

Based on the evaluation context described above, one can notice that the PLCC associated with the TMQI1 [START_REF] Yeganeh | Objective Quality Assessment of Tone Mapped Images[END_REF] is highest (0.588) than those with TMQI2 [START_REF] Ma | High Dynamic Range Image Compression by Optimizing Tone Mapped Image Quality Index[END_REF] (0.512) and FSITM [START_REF] Nafchi | FSITM: A Feature Similarity Index For Tone-Mapped Images[END_REF] (0.018). According to the provided results, the TMQI1 [START_REF] Yeganeh | Objective Quality Assessment of Tone Mapped Images[END_REF] (called TMQI D) is retained but requires some changes. The parameters as tuned in the paper [START_REF] Yeganeh | Objective Quality Assessment of Tone Mapped Images[END_REF] are not well adapted since they have been trained only on a small set of TMOs (only 8) and a reduced test HDR images (only 15). Indeed this metric does not take into account the new distortions introduced by the new TMOs (e.g. "Duan" [START_REF] Duan | Tone-mapping High Dynamic Range Images by Novel Histogram Adjustment[END_REF], "Fattal" WRB [START_REF] Fattal | Edge-Avoiding Wavelets and their Applications[END_REF], "Li" [START_REF] Li | Compressing and Companding High Dynamic Range Images with Subband Architectures[END_REF], "Husseis" [START_REF] Husseis | Revisited Histogram Equalization as HDR Images Tone Mapping Operators[END_REF], "NUHA" [START_REF] Thai | Piecewise linear perceptual quantizer as a non-uniform histogram equalization adjustment for contrast enhancement of tone mapped HDR images[END_REF], "SEPENO" [START_REF] Thai | Performance Evaluation of High Dynamic Range Image Tone Mapping Operators Based on Separable Non-linear Multiresolution Families[END_REF], "NON-SEPENO" [START_REF] Thai | Image Tone Mapping Approach Using Essentially Non-Oscillatory Bi-quadratic Interpolations Combined with a Weighting Coefficients Strategy[END_REF], "CEDP" [START_REF] Thai | Contrast Enhancement and Details Preservation of Tone Mapped High Dynamic Range Images[END_REF]). That is why this paper addresses, in section III-C, the question of re-evaluating the TMQI parameters to be in accordance with the MOS.

C. Revised TMQI quality metric

This section computes the best values of the parameters of the TMQI metric (i.e. a, α, β) in an iterative way as described in the The Absolute Mean Error (AME) between the metric and normalized MOS (using the same training database) has been reduced from 0.3758 (TMQI D) to 0.2626 (TMQI P).

Fig. 2a and Fig. 2b show that the parameters are better designed on our training database ((a) "DollDoll" HDR test image (each group has 15 corresponding TMOs) as a column matrix, (b) "Ward" method (each group has 24 corresponding images) as a row matrix). The absolute mean errors are reduced for the image "DollDoll" and method "Ward".

Fig. 3 presents the TMQI metric versus the normalized MOS. One can observe that the values related to the points (TMQI P, normalized MOS) are more gathered around the bisector line than the points (TMQI D, normalized MOS). The more population of points near the bisector line the better the quality evaluation metric.

IV. SIMULATION RESULTS

This section evaluates the visual quality of the TM HDR images using the default TMQI and the revisited TMQI (i.e TMQI P) metrics. The resulting metrics are then compared to the MOS. Simulation results have been conducted on a large set of TM HDR images not belonging to the TM image training database to show the evaluation performance of our metric. TM "Lausanne", "WreathBuilding", "AtriumNight", "LondonChapel", "AdobeLobby" and "DomeBuilding" HDR images are selected to discuss the results. The respective TM images are provided by Fig. 8, Fig. 9, Fig. 10. One can observe that The visual quality evaluated by the revisited metric is closer to the MOS than the default TMQI metric.

V. CONCLUSION

This paper revised the TMQI objective metric to evaluate the visual quality of the tone mapped HDR images. The new values of the

Fig. 1 .

 1 After the training process, the best values of the parameters are set to: a = 0.1, α = 0.1, β = 0.2 and P LCC = 0.7120. One can notice that the Pearsons Linear Correlation Coefficient with the corresponding MOS increases from 0.588 (denoted TMQI D for Default TMQI) to 0.7120 (denoted TMQI P for revisited TMQI).

Fig. 1 .

 1 Fig. 1. Iterative algorithm to tune the TMQI parameters using PLCC criterion.

Fig. 4 ,

 4 Fig. 5, Fig. 6, Fig. 7 confirm the good choice of the tuned parameters when the visual evaluation is carried out on the TM HDR images taken from the training database.

Fig. 2 .

 2 Fig. 2. Absolute Mean Errors (AME) as one of our well-designed examples.

Fig. 3 .

 3 Fig. 3. TMQI D and TMQI P metrics versus the normalized MOS scores.

Fig. 4 .

 4 Fig. 4. "DollDoll" HDR test image (13.89 f-stops) -CEDP Opt with level 1, Norm 1 with MOS=0.813, TMQI D=0.843 and TMQI P=0.808 (close to the MOS).

Fig. 5 .

 5 Fig. 5. "BottlesSmall" HDR test image (16.03 f-stops) -NONSEP with level 2 with MOS=0.563, TMQI D=0.903 and TMQI P=0.889 (close to the MOS).

Fig. 6 .

 6 Fig. 6. "OxfordChurch" HDR test image (15.43 f-stops) -CEDP Opt with level 1, norm 1 with MOS=0.738, TMQI D=0.878 and TMQI P=0.788 (close to the MOS).

Fig. 7 .

 7 Fig. 7. Left: "Synagogue" HDR test image (8.57 f-stops) -CEDP Lin with level 4, norm 1 with MOS=0.750, TMQI D=0.908 and TMQI P=0.869 (close to the MOS), Right: "PeaceRocks" HDR test image (24.13 fstops) -SEPENOCA with level 1 with MOS=0.613, TMQI D=0.841 and TMQI P=0.788 (close to the MOS).

Fig. 8 .

 8 Fig. 8. Left: "Lausanne1" HDR test image (7.71 f-stops) -SEPENOCA with level 2 with MOS=0.625, TMQI D=0.835 and TMQI P=0.754 (close to the MOS), Right: "WreathBuilding" HDR test image (23.88 f-stops) -Fattal RBW with MOS=0.560, TMQI D=0.908 and TMQI P=0.589 (close to the MOS).

Fig. 9 .

 9 Fig. 9. Left: "AtriumNight" HDR test image (28.68 f-stops) -CEDP Opt with level 5, norm 1 with MOS=0.412, TMQI D=0.768 and TMQI P=0.562 (close to the MOS); Right: "LondonChapel" HDR test image (14.19 f-stops) -Husseis with MOS=0.850, TMQI D=0.722 and TMQI P=0.853 (close to the MOS).


Ba Chien Thai
email: bachien.thai@univ-paris13.fr

Anissa Mokraoui
email: anissa.mokraoui@univ-paris13.fr

Basarab Matei
email: matei@lipn.univ-paris13.fr

Mahdia Tabal

IMAGE TONE MAPPING APPROACH USING ESSENTIALLY NON-OSCILLATORY BI-QUADRATIC INTERPOLATIONS COMBINED WITH A WEIGHTING COEFFICIENTS STRATEGY

This paper proposes a Tone Mapping Operator (TMO) to convert High Dynamic Range (HDR) images into Low Dynamic Range (LDR) images. It is related on two consecutive stages. The first one decomposes the HDR image on several resolution levels according to a nonseparable multiresolution approach using Essentially Non-Oscillatory (ENO) strategy where adapted bi-quadratic interpolations are exploited. The underlying idea of this decomposition is to better represent the details of the complex HDR image preserving as much as possible the HDR image quality. In the second stage, a weighting coefficients process is deployed to reduce judiciously the dynamic range of each subband's coefficient. The weights are deduced from Gaussian filtering operations combined with a monotonic decreasing function where its parameters depend on the statistical properties of each HDR subband, neighborhood and resolution level. Simulation results show that the proposed image TMO has the ability to better represent the details compared to other operators.

INTRODUCTION

High Dynamic Range (HDR) images are considered as high-quality images on which different levels of exposure areas of the real-world scene are rendered at the same time faithfully to the Human Visual System (HVS). However these images cannot be viewed on standard Low Dynamic Range (LDR) display devices since their dynamic range is smaller than those of HDR images. Unfortunately HDR display devices are still too expensive to replace all LDR devices. To overcome this problem, Tone Mapping Operators (TMOs) have been developed [START_REF] Banterle | Advanced High Dynamic Range Imaging : Theory and Practice[END_REF]. Their main purpose is to reduce the dynamic range (contrast, color gamut, details...) of HDR images to LDR dynamic range device while preserving as much as possible the captured scene appearance. For more details, a complete state of the art is available in reference [START_REF] Banterle | Advanced High Dynamic Range Imaging : Theory and Practice[END_REF] where a classification into local, global, segmentation, frequency and perceptual operators is proposed. Although each TMO approach has its own specificity, it is not possible to review the wide range of the developed work in this paper. Only those selected (according to their good performance) to be compared to our approach are summarized below.

In [START_REF] Durand | Fast Bilateral Filtering for The Display of High-dynamic-range Images[END_REF], Durand and Dorsey proposed an edgepreserving bilateral filter to decompose the HDR image into two layers : a base layer encoding large-scale variations and a detail one. Contrast is then reduced only in the first layer while the details are kept unchanged. The combination of these tone mapping layers produces an LDR image. In [START_REF] Drago | Adaptive Logarithmic Mapping for Displaying High Contrast Scenes[END_REF], Drago et al. presented an adaptive logarithmic mapping method of luminance values. It concerns the adaptive adjustment of the logarithmic basis depending on the radiance of the pixels. A set of logarithmic functions ranging from log 2 to log 10 were used to preserve scene details and improve the rendering contrast. A bias power function was used to ensure a smooth interpolation between the different logarithm bases. In [START_REF] Li | Compressing and Companding High Dynamic Range Images with Subband Architectures[END_REF], Li et al. proposed a subband architecture related on an oversampled Haar pyramid representation. Subband coefficients are rescaled using a gain control function reducing the high frequency magnitudes and boosting low ones. The modified subbands are then convolved with the Haar synthesis filters and summed to reconstruct the final LDR image. In [START_REF] Duan | Tonemapping High Dynamic Range Images by Novel Histogram Adjustment[END_REF], Duan et al. proposed an optimization approach based on a histogram adjustment between linear mapping and equalized histogram mapping. In [START_REF] Fattal | Edge-Avoiding Wavelets and their Applications[END_REF], Fattal et al. proposed a second generation of wavelets based on the edge content of the image avoiding having pixels from both sides of an edge. This approach was then exploited to map an HDR image into an LDR image. In [START_REF] Thai | Performance Evaluation of High Dynamic Range Image Tone Mapping Operators Based on Separable Nonlinear Multiresolution Families[END_REF], the authors proposed non-linear separable multiresolution families using data dependent interpolation where a set of operators have been derived. These latter have the ability to consider in their mathematical models the image singularity points. However the weakness of this approach is that the detail and approximation coefficients are weighed in the same way without taking into account the statistical properties of the coefficients. To improve the performance of this TMO, this paper investigates on non-separable and nonlinear Essentially Non-Oscillatory (ENO) multiresolution where each subband's coefficient is weighted according to the statistical properties of the subband.

This paper is organized as follows. Section 2 intro-duces the basic concepts of the non-linear multiresolution approach using data dependent interpolation. Section 3 presents the non-separable and non-linear multiresolution image TMO. Section 4 evaluates the TMO approach performance. Section 5 concludes this paper.

BASIC CONCEPTS OF HARTEN'S NON-LINEAR MULTIRESOLUTION

This section introduces the non-linear multiresolution approach based on Essentially Non-Oscillatory (ENO) interpolation strategy developed by Harten [START_REF] Harten | Multiresolution Representation of Data : A General Framework[END_REF], [START_REF] Aràndiga | Nonlinear Multiscale Decompositions : The Approach of A.Harten[END_REF]. Cell-Average (CA) multiresolution is considered since it has the ability to introduce in its mathematical model the isolated singularities such as edge points in the image avoiding the Gibbs phenomenon particularly harmful in HDR tone mapped images.

Denote j the resolution level. Harten's multiresolution strategy is based on two discrete interlevel operators, called "projection" D j-1 j and "prediction" P j j-1 operators. The pixel value v j k , located at position k = (k 1 , k 2 ) and resolution level j, is defined by the cell average of the underlying bidimensional function v modeling the image, as follows :

v j k = 2 2j C j k v(x, y)dxdy, (1)
where

C j k is the cell defined as [2 -j k 1 , 2 -j (k 1 + 1)] × [2 -j k 2 , 2 -j (k 2 + 1)] with 0 ≤ k 1 , k 2 ≤ 2 j -1.
The projection operator D j-1 j , acting from the fine to coarse level, computes the dataset v j-1 = D j-1 j v j corresponding to the coarser version of v j . In the cell average context, the operator D j-1 j is the averaging operator and is always linear :

v j-1 k = 1 4 (v j 2k + v j 2k+e1 + v j 2k+e2 + v j 2k+e1+e2 ), (2)
where e 1 and e 2 are unit vectors oriented to the right and upward, respectively. The prediction operator P j j-1 acts from the coarse to fine levels by computing an approximation vj of v j from v j-1 : vj = P j j-1 v j-1 . This operator may be nonlinear. Both prediction and projection operators satisfy the following consistency condition : D j-1 j P j j-1 = I, where I is the identity operator. The null space of D j-1 j represents the detail space of dimension N j -N j-1 . The prediction error j = v j -vj belongs to the detail space and then expanding it on a basis of that space getting the detail vector (d j-1 k ) k=1,..,Nj -Nj-1 . Therefore, one can equivalently write v j in the form (v j-1 , d j-1 ). Iterating this process, one obtains a multiresolution representation of the representation v J at the scale J into (v 0 , d 0 , d 1 , ..., d J-1 ).

NON-SEPARABLE MULTIRESOLUTION IMAGE TMO USING BI-QUADRATIC INTERPOLATIONS

This section deals with the proposed HDR image TMO given by the scheme of Figure 1. The first stage is based on the non-separable decomposition developed in [START_REF] Matei | Nonlinear Cell-Average Multiscale Signal Representations : Application to Signal Denoising[END_REF], [START_REF] Matei | Nonlinear and Nonseparable Bidimensional Multiscale Representation Based on Cell-Average Representation[END_REF]. The prediction operator is related on the smoothness properties of the function v, considered as a bi-dimensional smooth function and represented by a piecewise polynomial function of degree 2, modeling the complex HDR image (i.e. smooth regions separated by edges). To adapt the prediction near the singularities, ENO interpolation techniques are then deployed.

Prediction step

Consider a shift vector defined as r ∈ {-e 1e 2 , -e 1 , -e 1 + e 2 , -e 2 , 0, e 2 , e 1 -e 2 , e 1 , e 1 + e 2 }. Define for a given cell C j k and a shift r the set of cells :

S j (k, r) = {C j k+l+r ; l = (l 1 , l 2 ), -1 ≤ l 1 , l 2 ≤ 1 and k = (k 1 , k 2 ), -1 ≤ k 1 , k 2 ≤ 1}, (3)
called stencil of the cell C j k associated to the shift r. The stencil S j (k, r) has 3 × 3 cells and contains the cell C j k . An example of such stencil is the center stencil obtained for k = 0 and shift r = 0, i.e. :

S j (0, 0) =    C j (-1,-1) C j (0,-1) C j (+1,-1) C j (-1,0) C j (0,0) C j (+1,0) C j (-1,+1) C j (0,+1) C j (+1,+1)    . (4)
Therefore 9 stencils of size 3 × 3, containing the cell C j k associated to the different shifts r, exist. To emulate the ENO strategy, the cost function ζ for the stencil S j (k, r) is defined as :

ζ(S j (k, r)) = l∈S j (k,r) |∆v j-1 l |, (5)
where

∆v j-1 l = v j-1 l+1 -v j-1 l
and l = (l 1 , l 2 ). This cost function measures the sum of oscillations in S j (k, r). The least oscillatory stencil, for the cell C j k , is retained as follows :

r * = argmin -M +1≤r1,r2≤M -1 ζ(S j (k, r)). (6)
The least oscillatory stencil S j (k, r * ) is then retained. The predicted values vj are deduced from the values at the coarse level v j-1 k on S j (k, r * ). For this, a bi-quadratic polynomial p k,r * , defined by its 9 unknown coefficients p mn , is selected :

p(x, y) = 0≤m,n≤2 p mn x m y n , (7)
to predict the values as follows : where w is a vector basis in two dimensions w ∈ {0, e 1 , e 2 , e 1 + e 2 }. The 9 unknown coefficients are deduced from the 9 interpolation conditions on the stencil S j (k, r * ) defined by :

vj 2k+w,r = 2 2j I j w p k,r * dxdy, ( 8
) Split v j (x,y) D j j-1 ENO Prediction P j j-1 v j 2k v j 2k+e1 v j 2k+e2 v j 2k+e1+e2 v j 2k v j 2k+e1 v j 2k+e2 v j 2k+e1+e2 h d hv h h d d j-1 (x,y) dv j-1 (x,y) d h j-1 (x,y) v j-1 (x,y) d d j-1 (x,y) dv j-1 (x,y) d h j-1 (x,y) g(s) F v j-1 (x,y)
v kl = k k-1 l l-1
p k,r * (x, y)dxdy.

(9)

Detail step

The prediction errors, with respect to v j , are by construction :

j w = v j w -vj w with w ∈ {0, e 1 , e 2 , e 1 + e 2 }. (10)
Due to the consistency property D j-1 j P j j-1 = I, we get that w∈{0,e1,e2,e1+e2} j w = 0. These errors can be written in a non-redundant way by representing them in a basis of of KerD j-1 j . The matrix of this change of basis H has the columns

h lo = [ 1 4 , 1 4 , 1 4 , 1 4 ], h h = [ 1 4 , -1 4 , 1 4 , -1 4 ], h v = [ 1 4 , 1 4 , -1 4 , -1 4 ],
and

h d = [ 1 4 , -1 4 , -1 4 , 1 4
]. Therefore the details vectors are deduced from j 2k , j 2k+e1 , j 2k+e2 , j 2k+e1+e2 as follows :

d j-1 h = {h h × [ j 2k , j 2k+e1 , j 2k+e2 , j 2k+e1+e2 ] }, (11
) for horizontal detail block ;

d j-1 v = {h v × [ j 2k , j 2k+e1 , j 2k+e2 , j 2k+e1+e2 ] }, (12)
for vertical detail block ; and

d j-1 d = {h d × [ j 2k , j 2k+e1 , j 2k+e2 , j 2k+e1+e2 ] }, (13)
for diagonal detail block.

Note that the size of the block is half of each direction of v j . The multiresolution representation is then obtained

v j := (v j-1 , (d j-1 d , d j-1 h , d j-1 v )).
Iterating the process until resolution level J involves a set of approximation and detail coefficients :

(v 0 , (d 0 d , d 0 h , d 0 v ), (d 1 d , d 1 h , d 1 v ), ..., (d J-1 d , d J-1 h , d J-1 v
)).

Adaptive weighting coefficients strategy

This section proposes to independently weight each subband's coefficient to reduce in an appropriate way the HDR image dynamic range.

The proposed weight is a positive value between 0 and 1 depending on the statistical properties of each subband and the neighborhood. The main operations are presented below.

First, the absolute value of all subband coefficients is taken :

A 0 0 = |v 0 |, A j d = |d j d |, A j h = |d j h |, A j v = |d j v | for j = 0, ..., J -1.
Denote A j sb (x, y), with sb = 0, d, h, v and j = 0, ..., J -1, the positive coefficient at position (x, y) of the subband sb at resolution level j. Define g(x, y) the Gaussian filter given by :

g(x, y) = 1 2πσ 2 sb e -x 2 +y 2 2σ 2 sb , ( 14
)
where σ sb is the standard deviation of the Gaussian distribution of the subband sb deduced from the median component as in [START_REF] Matei | Nonlinear Cell-Average Multiscale Signal Representations : Application to Signal Denoising[END_REF] :

σ sb = median(A j sb ) 0.6745 . ( 15
)
This filter is used to compute the weighting coefficient a j sb (x, y) associated to each coefficient A j sb (x, y) in the corresponding domain as follows :

a j sb (x, y) = A j sb (x, y) * g(x, y), ( 16
)
where * is the convolutional product. After this, a monotonic nonlinear decreasing function F (as in [START_REF] Li | Compressing and Companding High Dynamic Range Images with Subband Architectures[END_REF]) is applied on the coefficients a j sb , resulting in the positive weighting coefficient q j sb (x, y) between 0 and 1, as follows :

q j sb (x, y) = F (a j sb (x, y)) = 1 ( a j sb (x,y)+ε δj ) 1-γ , ( 17
)
where γ ∈ [0, 1] is the compressive factor ; j the resolution level (i.e. j = 0, ..., J -1) associated to the subband sb. ε is used to avoid singularities when the power 1 -γ is a positive value. δ j is considered as a gain control stability level depending on the resolution level j and the subband coefficients :

δ j = [1 - j(1 -ξ) J ] × Āj sb , (18)
where Āj sb is the average of A j sb . Finally, the subband coefficients are weighted as follows :

v 0 = q 0 0 ⊗ v 0 and d j sb = q j sb ⊗ d j sb (19)
with sb = d, h, v ; j = 0, ..., J -1 ; and ⊗ is the product term by term.

SIMULATION RESULTS

This section evaluates the performance of the proposed TMO using the TMQI (Tone-Mapped image Quality Index) metric developed in [START_REF] Yeganeh | Objective Quality Assessment of Tone-mapped Images[END_REF]. It evaluates the quality of the HDR image tone mapped compared to the original HDR image. Note that this metric is upper-bounded by 1.

Simulations have been conducted under Matlab environnement using the HDR Toolbox with its HDR test images : "Anturium", "Bottle Small" renamed "Bottle", "Office", "Oxford Church", "Memorial" and "Light" [START_REF] Banterle | Advanced High Dynamic Range Imaging : Theory and Practice[END_REF]. The different parameters are chosen so as to give the best results in terms of TMQI metric in all methods. These latter are performed on the logarithmic transformation of the Luma values L w (i.e log 10 (L w + ε), ε avoiding the logarithmic singularities).

Simulation results are compared to competitive methods : (i) SEP [10] with β = 0.3, γ = 0.7, J=2, upd=0, b = e ; (ii) Li TMO [11] with Haar multiscale ; (iii) Fattal [9] using RBW method with α = 0.8, β = 0.3, γ = 0.8, J=2, upd=0 ; (iv) Drago [2] , Reinhard [3] , Ward [4] , Durand [5] , Tumblin [6] , Schlick [7] with the default parameters as given in the HDR Toolbox ; and (v) Duan [8] using β = 0.5. For the proposed method, denoted "Proposed", simulations are conducted with the following parameters : ε = 10 -9 , J = 2, γ = 0.6, ξ = 0.1, σ i deduced from equation [START_REF] Matei | Nonlinear and Nonseparable Bidimensional Multiscale Representation Based on Cell-Average Representation[END_REF]. Figure 2 clearly shows the impact of the second stage on the visual quality of the LDR image. To further improve the visual quality of some HDR tone mapped images, it is possible to add the histogram adjustment, developed in [START_REF] Duan | Tonemapping High Dynamic Range Images by Novel Histogram Adjustment[END_REF], as a post-processing. This postprocessing can be applied on the HDR tone mapped image in a local way on equal non-overlapping blocks of size 32×32 or in a global way. Indeed this approach adjusts the contrast between linear mapping and histogram equalized mapping thanks to a control parameter β ∈ [0, 1] (we set β = 0.3). The two approaches are named "Proposed G" for global and "Proposed L" for local post-processing.

Table 1 provides the TMQI of the HDR tone mapped images according to the above listed methods. One can observe that our approach is competitive to those developed in the literature.

Figure 3 shows the "Bottle" HDR tone mapped images according to the "Proposed", "Proposed G" and "Proposed L" approaches. Figure 4 presents the "Bottle" HDR tone mapped images using [START_REF] Duan | Tonemapping High Dynamic Range Images by Novel Histogram Adjustment[END_REF] with local ("Duan L") and global ("Duan G") methods. The visual rendering of our tone mapped images is better than those developed in [START_REF] Duan | Tonemapping High Dynamic Range Images by Novel Histogram Adjustment[END_REF], even when the post-processing is not used. This is observed even if the TMQI is higher for "Duan G" (0.91) and "Duan L" (0.92) than our "Proposed" (0.87) (see Table 1) where "Proposed L" TMQI remains high (0.93). Note that the TMQI metric does not always reflect the visual quality.

Figures 5,6, 7 and 8 provide the visual quality of the "Memorial" HDR tone mapped images. One can observe that the details are well rendered in favor of our method (e.g. tiles, rosette). In Figure 6, the tiles on the top right on the tone mapped image are better rendered in our method.

Based on these simulation results, the proposed operator has the ability to better represent the details compared to other operators. Its strength points are inherent not only to the non-separable multiresolution decomposition but also to the ENO strategy. It suffices to compare its performance with that of separable multiresolution ENO strategy developed in [START_REF] Thai | Performance Evaluation of High Dynamic Range Image Tone Mapping Operators Based on Separable Nonlinear Multiresolution Families[END_REF].

TMOs

Anturium Bottle Office Church Memorial Light Drago [2] 0.874 0.801 0.801 0.814 0.800 0.800 Reinhard [3] 0.778 0.807 0.826 0.789 0.792 0.794 Ward [4] 0.806 0.783 0.775 0.817 0.793 0.789 Durand [5] 0.811 0.892 0.825 0.929 0.814 0.800 Tumblin [6] 0.715 0.713 0.735 0.675 0.761 0.750 Schlick [7] 0.770 0.835 0.926 0.970 0.788 0.780 Duan G [8] 0.964 0.916 0.955 0.986 0.935 0.969 Duan L [8] 0.910 0.924 0.927 0.906 0.839 0.924 Fattal [9] 0.889 0.928 0.928 0.889 0.927 0.971 SEP [10] 0.888 0.931 0.928 0.885 0.903 0.960 Li [11] 0

CONCLUSION

This paper focused on HDR image tone mapping operator. As a first stage, the operator decomposes the HDR image into several resolution levels according to nonseparable orthogonal transform using ENO bi-quadratic interpolation. This latter has the advantage to take into its mathematical model the possible singularities in all directions. Therefore the details are well represented in the subbands. As a second stage, a weighting strategy adapted to each subband's coefficient is developed in order to reduce adequately the dynamic range of the HDR image while preserving the more important details in the HDR image. Simulation results show good performance achieved not only in terms of TMQI metrics but also in terms of visual quality.

Fig. 1 .

 1 Fig. 1. Forward non-separable multiresolution and adaptive weighting strategy

Fig. 2 .

 2 Fig. 2. (Up) Without adaptive weighting ; (Down) With adaptive weighting (ε = 10 -9 ; J = 2 ; γ = 0.6 ; ξ = 0.1).

Fig. 3 .

 3 Fig. 3. (Up) Proposed (ε = 10 -9 ; J = 2 ; γ = 0.6 ; ξ = 0.1) ; (Middle) Proposed with global post-processing ; (Down) Proposed with local post-processing.

Fig. 4 .

 4 Fig. 4. (Up) Duan with global processing ; (Down) Duan with local processing.

Fig. 5 .

 5 Fig. 5. (Left) Duan (β = 0.5, local) ; (Right) Proposed with local post-processing.

Fig. 6 .

 6 Fig. 6. (Left) Fattal (RBW, α = 0.8, β = 0.3, γ = 0.8, J=2, upd=0) ; (Right) Proposed with local postprocessing.

Fig. 7 .

 7 Fig. 7. (Left) Li (Haar multiscale) with post-processing ; (Right) Proposed with local post-processing.

Table 1 .

 1 Tone Mapped Image Quality Index (TMQI)

		.964	0.954 0.854	0.877	0.834	0.888
	Proposed	0.938	0.873 0.948	0.820	0.846	0.932
	Proposed G	0.975	0.927 0.828	0.968	0.896	0.946
	Proposed L	0.920	0.930 0.861	0.939	0.881	0.909