Keywords

1 Introduction

Nowadays, high dynamic range (HDR) image technology draws lots of attention from researchers and groups all around the world. The real scenes in the nature often have a range of light variance much broader than 255, thus causing the conventional cameras not able to capture the whole range in only a single shot. Scenes including both bright reflection and dark shadows are always difficult for photographers to capture. So, kinds of software develop the post-processing methods to synthesize HDR images for those photographers.

Promoted by the prosperity of HDR imaging, lots of commercial software and image signal processing (ISP) hardware develop the function to get HDR images. Among them, ‘software’ methods always recover the HDR scene radiances from a stack of low dynamic range images (LDRs) exposed at different periods. For example, [1, 2] solves camera response function (CRF) first, and then weights and sums the LDRs together, [3] fuses different LDRs together utilizing some characters such as contrast, saturation and exposure time, [4] averages float LDR pixel values by specially designed weighting function, [6] takes an original image with high resolution but low dynamic range as input, and generates an output image with extended dynamic range by applying the trade-off between the image resolution and the dynamic range, [14] develops a linear relation between with two LDRs and then fuse them together to a HDR one based on the linear function and [16] uses the gradient information of the visible regions and recovers the saturated regions by the energy minimization approach. However, always limited by time cost and delays, ‘ISP’ methods are less than ‘software’ methods and need specially design. One of them first gets the exposure ratio k by dividing the high exposure time by low exposure time as k = e1/e2, and then multiplies the low exposed pixels by the exposure ratio. The whole synthesis equation can be written as XHDR = Xlow*k + Xhigh. However, this method has some drawbacks, such as bad fidelity and low brightness.

In this paper, we propose an algorithm to synthesize HDR image based on the logarithm intensity mapping function (IMF), which denotes the relation between intensity values of same pixel location from LDRs. We first calibrate the logarithm IMF curve using least square method and then get HDR image by weighting and summing the mapped intensities from LDRs. To eliminate ghosting of HDR and avoid aligning LDR images, we capture different exposed images using spatially varying pixel exposures (SVEs) in one shot. Thus, we develop a low-computational-cost and ISP-suited method to synthesize HDR image. To illustrate the main method, the rest parts of our paper are organized as follows. In Sect. 2, the basic theories of HDR imaging are introduced such as camera response function, intensity mapping function, SVE and cross-histogram. Section 3 outlines theoretical framework of our main algorithm including the solving of IMF and the designing of weighting function. The experiment results are presented in Sect. 4 and the further discussions are included in Sect. 5.

2 Basic Knowledge of High Dynamic Imaging Work

In work [2, 10], the acquisition model of camera has been proposed. Given the scene radiance L, the irradiance E reach into the camera can be derived by E = epL, where e is the exposure and p is an optical factor. The exposure e can be further expressed as e = (πd2)t, where t is the exposure time and d is the aperture size. Then, the pixel value X of image and irradiance E are related by the non-linear camera response function (CRF) f expressed as:

$$ X = f\left( E \right) $$
(1)

And works [2, 4] hold the point view that obtaining the HDR image is a matter of recovering scene radiance E from the intensity values X by reverse response function g:

$$ g = f^{ - 1} $$
(2)

Finally, the obtained HDR pixel value Xi is expressed as follows in [2, 4].

$$ X_{i} = \frac{{\sum\limits_{j} {w(i,j)E(i,j)} }}{{\sum\limits_{j} {w(i,j)} }} $$
(3)

In the above equation, w (i, j) is a weighting function, i is the index of pixel and j is the index of differently exposed image.

To solve out the reverse camera response function g, kinds of methods are utilized. The primary method makes use of uniformly illuminated chart with known reflectance patches and computes the function g based on the relationship between illumination and intensity value. However this method appears inconvenient for most situations. So, chart-free methods are introduced such as sampling and fitting [2], solving intensity mapping function (IMF) τ instead [5] and so on. These CRF-based methods are practical when built in software, but hardly applicable to real time system due to its huge computational cost and time delay.

To establish our method, IMF τ is deployed instead of function g. In two differently exposed LDR images, the intensity value in the corresponding pixel location can be paired together as (X1, X2), where X1 is from low exposure and X2 denotes is from high exposure. The number of these pairs is expressed as J(X1, X2) and shown in the figure called cross-histogram (see Fig. 2). And function τ well describes the relation between values X1, and X2. In [5], it has the expression as

$$ X_{2} = \tau (X_{1} ) . $$
(4)

Some important properties of IMF τ are as follows from [5]:

  • τ(0) = 0

  • τ monotonically increases

  • X <=τ (X)

  • \( \mathop {\lim }\limits_{n \to \infty } \tau^{ - n} (X) = 0 \)

As explained in [5], due to quantization and saturation of the pixel value, the cross-histogram curve may seem not so smooth and monotonic increasing as the theory describes. So it introduces accumulative histogram function and obtains the IMF τ. Different with algorithm from [5], we assume IMF have a logarithm formation and propose a fitting algorithm to calibrate the intensity mapping function after a partition of cross-histogram. Further on, we develop a weighting function and utilize it with IMF to construct HDR image. In the end, we test our algorithm on computer and get satisfying HDR image. The detail algorithm will be included in the following chapters and we will show that our algorithm is a better alterative suiting real-time ISP hardware system.

3 Algorithm to Construct HDR Image

We divide our method into several sections. Section 3.1 introduces our technology to get LDRs. Section 3.2 includes the algorithm to solving the intensity mapping function from the cross-histogram. The methods to synthesis HDR image based on the IMF is shown in Sect. 3.3. Section 3.4 takes a tone-mapping strategy to show HDR image on LDR screen.

3.1 Acquisition of Spatially Varying Pixel Exposures

Spatially varying pixel exposure (SVE) has been introduced by many CMOS systems and used in differently-exposed frames acquisition. To obtain SVEs, optical masks with spatially varying transmittance are put adjacent to the CMOS sensor (see Fig. 1a). And different number of pixels in one pattern can be applied in different CMOS systems. The optical masks control the amount of light transmits through them and producing differently exposed frames in one shot [8, 9]. With SVEs, we easily acquire the different exposed images. And thanks to the short duration of one shot, the effect of object motion and scene motion can be eliminated and thus avoid the global alignment before synthesis. Here, we make use of SVE masks with two different-exposed pixels in one pattern and obtain two frames, one long exposure and one short exposure, in one shot (see Fig. 1b).

Fig. 1.
figure 1figure 1

SVE pattern with 4 exposures is shown on the left (marked as Fig. 1a). The right one (marked as Fig. 1b) shows the procedure of synthesis of HDR pixels with 2 exposures.

Full size image

3.2 Recovery of Intensity Mapping Function

The Grossberg’s Algorithm.

The cross-histogram of two different exposures contains all information to construct the intensity mapping function as proved in [5]. The function τ can be derived from the image’s cumulative histogram called H(x) which denotes the number of pixels with value between 0 and x. Given the continuous histogram function h(x) which indicates the number of pixels with intensity value x, H(x) can be solved by:

$$ H(X) = \int_{0}^{X} {h(u)du} $$
(5)

Actually the histogram function h(x) is discrete, so the above integration needs to be turned to accumulation. After the accumulation, the corresponding pixel values X1 and X2 from the same pixel location of two images have the following relation as (6) shows. It implies the numbers of pixels with smaller values than them is exactly the same [5].

$$ H_{1} (X_{1} ) = H_{2} (X_{2} ) $$
(6)

And from IMF, we get \( X_{2} = \tau (X_{1} ) \). Then, substituting it into (6) yields

$$ H_{1} (X_{1} ) = H_{2} (\tau (X_{1} )) $$
(7)

Supposing X is arbitrary value from the low exposure, IMF τ can be derived from reverse function of H as follows.

$$ \tau (X) = H_{2}^{ - 1} (H_{1} (X)) $$
(8)

The Partitioning and Fitting Algorithm.

We first partition the cross-histogram into different parts and propose to calibrate the logarithm IMF in some of the parts. In [14], it illustrates a linear J-function which has similar sub-parts and properties with parts of IMF. Here we extend the partition strategy of J-function and apply it to our IMF. Partition of the cross-histogram is shown in Fig. 2b and explanations are as follows.

Fig. 2.
figure 2figure 2

Cross-histogram maps are shown. The top one (marked as Fig. 1a) is a 3D histogram, in which the Z-coordinate denotes the number of pixels have value (X1, X2) respectively. The bottom one (marked as Fig. 2b) is a 2D cross-histogram which has been partition into 4 parts.

Full size image

  • Part I. This part is named the low exposed zone. In this zone, the pixels’ abscissa values are almost the same and near a constant. These low abscissa values mean they are in the dark scene and lack of exposure in low exposed image. However, their ordinate values are not constants because they are exposed in a long period in the high exposed image.

  • Part II. We name this part the appropriate exposed zone. In this part, the abscissa value and ordinate value are both variables. Because this part of both images expose appropriately and the pixels here contains most information of the scene.

  • Part III. We call this part the noise zone. In this part, the pixel values are disperse from the ‘main curve’ due to the photoelectron noise and intensity quantization.

  • Part IV. We name this part the over exposed zone. Pixels in this part have the same ordinate values nearly 255. However, because exposed appropriately in low exposed image, the abscissa values are less than saturation.

We see that the low-exposed zone and over-exposed zone both loss some information from the scene. In reality, when exposed at a short period, the dark scene can’t be recognized, and when exposed at long period, the bright scene may show saturation and noising. The aim of HDR imaging is to include these dark and bright parts in one image by recover the intensity values from other zones. In one sense, the pixel values from zone I and IV are not reliable due to noise, so we chose zone II and III to recover the intensity mapping function. Different from [14], we extend the noise zone and then use a denoising method to make the data useful in this area.

Because all the values are discrete in the cross-histogram image, so we apply a fitting strategy to calibrate the curve that best satisfy the IMF in a least square error sense. And to recover the function means recovering a finite number of pixel values that satisfy X2 = τ(X1) in the range of (0, 255). To achieve it, we assume the IMF has a logarithm formation and solve the least squares problem in zone II and III. Here, we give the form of the function that

$$ \tau (x) = a_{0} *\log (x) + a_{1} $$
(9)

And we have the objective function written as:

$$ O = \sum\limits_{{(x_{i} ,y_{i} ) \in S}} {(\tau (x_{i} ) - y_{i} )^{2} } $$
(10)

In the objective function, xi denotes the abscissa value and yi denotes the ordinate value.

$$ S = \{ (x,y)|x \ge \varepsilon_{1} ,y \le \varepsilon_{2} \} $$
(11)

In it, ε1 and ε2 are the threshold values of zone II and III. To solve the function, we need to denoise the values in the zone III. Here we take a ‘denoising procedure’ by getting together the pixels have the same abscissa value X and then arranging their ordinate values from small to large. To get the “denoised” value Y, we get the middle 10 ordinate values and average them. This way, we eliminate the noised pixels in zone III and at the same time reduce the number of pixels involved in the later procedure from thousands to less than 255 (see Fig. 3a).

Fig. 3.
figure 3figure 3

Cross-histogram maps after denosing are shown. The top one (marked as Fig. 3a) shows the relation between pixels in two LDR images. The bottom one (marked as Fig. 3b) shows a logarithm IMF line after fitting (red continues line) and cross-histogram (blue circles) (Color figure online).

Full size image
$$ Y = \frac{1}{10}\sum\limits_{{i \in \{ i|x_{i} = X\} \, }} {y_{i} } $$
(12)

Then, to minimize the objective function in (10), we apply derivation on both sides and then solve the linear equations to get the two parameters. And the curve and the pixel values after fitting are shown in Fig. 3b. We see that in zone II and III, the logarithm curve fit the pixel values well.

$$ \left\{ {\begin{array}{*{20}l} {\frac{\partial O}{{\partial a_{0} }} = \sum\limits_{{(x_{i} ,y_{i} )}} {2*(\tau (x_{i} ) - y_{i} )\log (x_{i} ) = 0} } \hfill \\ {\frac{\partial O}{{\partial a_{1} }} = \sum\limits_{{(x_{i} ,y_{i} )}} {2*(\tau (x_{i} ) - y_{i} ) = 0} } \hfill \\ \end{array} } \right. $$
(13)

3.3 HDR Imaging Based on the Intensity Mapping Function

After the recovery of IMF, we apply a weighting algorithm to the dark and bright scene in LDR images and synthesize the different exposures together to HDR image. As we know, the low exposure includes the information of saturated scene and the high exposure includes the information of dark scene, so we proposed to fuse these parts together based on the intensity mapping function. In [14], it proposed the zone I, II and III of high exposure contain the needed information of scene and try to fuse the zone IV of low exposure to high exposure linearly.

Here, we first map the zone IV from low exposure to the high exposure by substituting the pixel values to IMF. Given the intensity value of low exposure Xlow, we have the value after mapping written as:

$$ X^{\prime}_{low} = \tau (x_{low} ) $$
(14)

Second, we develop a weighting function to fuse the different zone together. In the low exposure, the information in zone IV is weighted 1, while information in other zones is weighted 0. In the high exposure, information in zone I, II and III are weighted 1, while the other is weighted 0. In order to make the transition from zone III to IV smoother, we give the weighting function an S-shape increasing (and decreasing) between the intervals. The weighting functions W1 and W2 are shown below, in which \( \xi \) is the dividing value of zone III and IV. ε and α are the parameters which define the width of transition zone and can be manually changed to improve the visual sense.

$$ W_{1} = \left\{ {\begin{array}{*{20}l} 0 \hfill & {x_{1} < \xi - \varepsilon } \hfill \\ 1 \hfill & {x_{1} > \xi + \varepsilon } \hfill \\ {\frac{{x_{1} - \xi + \varepsilon }}{2*\varepsilon } *\frac{1}{1 + \alpha *(\exp ( - (x - \xi ))/2)}} \hfill & {others} \hfill \\ \end{array} } \right. $$
(15)
$$ W_{2} = \left\{ {\begin{array}{*{20}l} 1 \hfill & {x_{2} < \xi - \varepsilon } \hfill \\ 0 \hfill & { \, x_{2} > \xi + \varepsilon } \hfill \\ {1 - \frac{{x_{2} - \xi + \varepsilon }}{2*\varepsilon } *\frac{1}{1 + \alpha *(\exp ( - (x - \xi ))/2)}} \hfill & {others} \hfill \\ \end{array} } \right. $$
(16)

Then, we take the weighted average of the intensity value from the high exposure and “mapped” low exposure to get the result, HDR intensity value. This procedure can be expressed as follows:

$$ X_{HDR} (i,j) = \frac{{W_{1} (x^{\prime}_{1} (i,j))x^{\prime}(i,j) + W_{2} (x_{2} (i,j))x_{2} (i.j)}}{{W_{1} (x^{\prime}_{1} (i,j)) + W_{2} (x_{2} (i,j))}} $$
(17)

By the above procedures, we build up a pixel level synthesis algorithm to construct a HDR image from different SVEs. The LDR images and corresponding HDR images (tone-mapping operators are applied on HDRs) are shown in Figs. 5 and 6.

To extend our algorithm to synthesize more LDR images, we would like to choose a mid-exposed image as a reference map. Then, we obtain the IMF between the LDR stacks with the reference respectively and synthesize high dynamic stacks. Last, we weight and accumulate the HDR stack to a final HDR image (see Fig. 4). Although more LDRs are needed to get satisfied HDR image, we choose two LDRs for use because the time delay of hardware increases greatly due to the large LDR stack.

Fig. 4.
figure 4figure 4

The overall procedure of dealing with more pictures is shown in this figure.

Full size image

3.4 A Contrast-Preserving Tone-Mapping Operator

The aim of tone-mapping is to preserve contrast, brightness and details of a HDR image which is shown on a typical 255-intensity-level screen [11]. Here we use a global contrast-preserving mapping operator to reduce the intensity levels to 255.

First we normalize the pixel intensity from (0, Xmax) to (0, 1) by applying the equation:

$$ X_{f} \frac{{X_{i} }}{{X_{\hbox{max} } }} $$
(18)

After normalization, the pixel values all become float numbers. Then we take a few steps to apply the multipeak histogram equalization operator from [12] to enhance the contrast while compress the intensity value.

The algorithm in [12] detects N breakpoints {d1, d2 … dN} and N + 1 peaks of the histogram. Then the histogram B is divided into N + 1 subsections based on the breakpoints which can be expressed as:

$$ B = B_{1} \cup B_{2} \cup \ldots \cup B_{N + 1} $$
(19)

For each Bj, the probability density function can be written as follows [15], where i denotes the index of peak and j denotes the index of intensity level:

$$ p_{i} (j) = \frac{{X_{f} (j)}}{{\sum\limits_{{d_{i} \le j < d_{i + 1} }} {X_{f} (j)} }} $$
(20)

Then the histogram equalization process is applied separately to each of the peaks of histogram. The cumulative intensity function is written as follows:

$$ C_{i} (k) = \sum\limits_{{j = d_{i} }}^{k} {p_{i} (j)} $$
(21)

After the multipeak histogram equalization, the intensity levels of each peak are spread out between the breakpoints. Lastly, we extend the range of intensity value from (0, 1) to (0, 255). This way, we can preserve the contrast and compress the intensity levels.

4 Experiment Results

We test our algorithm on variety of scenes and then compare the HDR images with results from method of [2] and [6]. The respective low dynamic stacks are taken in different exposure time while other settings of camera are kept static during a shot.

In Fig. 5 the ‘office’ pictures which are taken under different exposures are shown in the first line. The low exposed image has dark areas in indoor scene like the book shelf and desk, while the high exposed image has saturated areas in the outdoor scene such as the tree and road. Shown in the second line, our HDR image has appropriately recovered missing information in LDRs. To contrast with the HDR image from our algorithm, the third line and fourth line show the ‘software’ results from [2] and [6]. Figure 6 shows another stack of images which include castles in bright light. In the first line are two LDRs, the second line is our HDR image while the third and last line are the ‘software’ results. In all our and [2] ’s HDR images, the missing areas clearly reappear and remain high contrast. However, some images based on [6] ’s method reappear with high quality, while some other images with excessive bright parts may have an even brighter appearance due to the summation of the intensity without any distinguish with pixels, and thus cause losing of some information.

Fig. 5.
figure 5figure 5

Four images show scenes of offices. In the first line is LDR images. In the second line is HDR image and details got by our algorithm. The contrast HDR images with details resulted from ‘software’ method [2] and [6] are shown in the third and fourth line.

Full size image
Fig. 6.
figure 6figure 6

These images show scenes of castles. In the first line is LDR images. In the second line is HDR image and details got by our algorithm. The contrast HDR images with details resulted from ‘software’ method [2] and [6] are shown in the third and last line.

Full size image

The experiment results suggest that our HDR images can keep good or even better visual quality than those ‘software’ results. In fact, our method needs not to know the exact exposure time of the LDR stacks and thus becomes more convenient to photographers. Another advantage is the low time cost, mostly because our technology needs not to resolve the camera response function, which is always the most time-costly part. (In fact, the method from [6] has similar running time as ours, but it trades the quality of HDR picture (see Figs. 5 and 6’s details) for fast speed. So we didn’t satisfy with the method even the running time is short.) We simulate both algorithm on the same computer and show the running time in Table 1. (We run the MATLAB code of our own and from [2] respectively.) Table 1 demonstrates different size of images varies greatly in the computational time. To construct a better HDR image with more texture and higher contrast, higher-resolution LDR images are needed. Totally, our IMF-based method shows a time-cost superiority towards ‘software’ method. Thus, our algorithm is expected to be more suitable to build in the graphic hardware in the future. In addition, it makes capturing of HDR videos possible with an ordinary camera.

Table 1. Computational time of our algorithm and ‘software’ methods.
Full size table

5 Conclusion and Discussion

HDR imaging has been a popular subject in the machine vision field in the past years. Lots of software techniques have been proposed to synthesize HDR image using a bracket of low dynamic exposures. However, most of them need to resolve the overdetermined camera response function and thus cost a lot of time. So we propose to synthesize HDR images based on the intensity mapping function. And then we apply a contrast-preserving tone-mapping operator to turn the intensity level range to (0, 255). Totally, our algorithm shows advantages of less computational cost to the ‘software’ methods and thus more suitable to be built in real-time ISP hardware.

In the future work, we would like to develop better denoising algorithm which can eliminate unreliable pixels and increase HDR’s SNR in the dark and saturated areas. Also, we would like to investigate different tone-mapping techniques to maintain better texture and more details.