An Image Sharpness Enhancement Algorithm Based on Green Function

2.1 Image enhancement

Many image enhancement algorithms can improve the subjective visual effect of the original image to a certain extent. But each algorithm has its own defects. For instance, histogram equalization [3] can greatly enhance the global contrast of the image, but at the expense of the contrast in local areas. Through histogram equalization, the gray values are easily combined excessively, resulting in the loss of some details of the original image. As a result, the enhanced image has a reduced quality. With the aid of histogram segmentation, Cao et al. [4] performed histogram equalization of each sub histogram to improve the brightness of low-illumination images; the histogram segmentation alleviates the problem of excessive enhancement to a certain extent. Qadar et al. [5] proposed a recursive separation weighted histogram equalization algorithm, which decomposes the histogram into several sub histograms, weighs each sub histogram, and then implement histogram equalization; however, the quality of the enhanced image often degrades, due to the loss of statistical information.

Some researchers held that the spatial information of the image is more suitable for image enhancement than the first-order statistics. For example, Zhou et al. [6] put forward the context and variation contrast enhancement algorithm, which adjusts the image with two-dimensional (2D) histogram. Under the best contrast tone mapping framework, Song et al. [7] proposed a new image enhancement algorithm that overcomes the limitations of contrast enhancement histogram, using the joint probability of spatial information.

2.2 Retinex model

The Retinex model [8], a popular model in the field of image enhancement, mimics human’s visual perception of brightness and colors. The model is built on the constancy of colors: Under the irradiation of natural light, different objects reflect lights of different wavelengths to human eyes, creating different colors; in other words, the colors are not intrinsic properties of the objects.

Currently, there are four types of Retinex models: path model [9], partial differential equation model [10], center surround model [11] and variational model [12]. The random spray Retinex model proposed by Provenzi et al. [13] has too many parameters and a huge computing load. The primary defect of the path model is that the conditions of local information processing no longer apply, when the path increases and the scope expands. Wu and Wu [14] analyzed Retinex model from a new perspective, and constructed a partial differential equation model. The model can achieve accurate results with only a few parameters. Xie et al. [15] proposed single scale Retinex (SSR) algorithm and multi-scale Retinex (MSR) algorithm, but the two algorithms have problems like unnatural color. Ma et al. [16] designed a multi-scale Retinex algorithm with color restoration ability, and corrected the enhanced image with the color recovery factor. Fine and MacLeod [17] presented a natural reservation algorithm based on prior brightness statistics, and obtained satisfactory results of quality enhancement. The variational model is based on prior construction constraints. Hou et al. [18] developed a Retinex algorithm under the variational framework, under the assumption that the illumination components vary smoothly in space, and calculated the illumination by gradient descent. Xie et al. [19] proposed a method to estimate the probability of illumination and reflection at the same time, and demonstrated its good results on detail enhancement. Li and Yin [20] evaluated the enhanced image by adding a coefficient matrix. This evaluation method does well in image enhancement, but the enhanced image is often noisy owing to the excessive brightness enhancement.

The above analysis shows that the Retinex model has a serious halo phenomenon in the enhanced image, despite adopting a hierarchical structure. The halo covers the edge details of the high-contrast image. The contrast can be obviously enhanced through nonlinear grayscale stretching. Nevertheless, this approach fails to consider the subjective visual effect of human eyes, and underestimates the image information entropy through the combination of gray values, leading to a poor visual effect of the enhanced image.

2.3 Thesis statement

To solve the defects of existing algorithms, this paper presents an image sharpness enhancement algorithm based on the Green function is proposed. To retain as much image details as possible, the illumination layer separated by Retinex model was enhanced separately without changing the detail layer. The former layer contains the details of the original image, while the latter covers the contrast information of the image.

The existence of halo is a sign of the serious gray mutation near the edge of the illumination image, for the algorithm emphasizes the spatial relationship over the gray relation between pixels. The halo can be effectively suppressed through bilateral filtering, which synthesizes the gray value change of pixels based on their relative positions.

The only difference between bilateral filtering and the SSR algorithm lies in the selection of convolution kernel. The convolution kernel of bilateral filtering can be defined as:

$\begin{aligned} L_{x, y}(i, j) &=\gamma \exp \left[-\frac{(i-x)^{2}+(j-x)^{2}}{2 \sigma_{s}^{2}}\right.\\ & \left.-\frac{[S(i, j)-S(x, y)]^{2}}{2 \sigma_{r}^{2}}\right] \end{aligned}$      (15)

where, $\sigma_{S}^{2}$ and $\sigma_{r}^{2}$ are empty value and gray value, respectively.

The shape of the bilateral filter depends on the gray value and position of pixels. Therefore, the time complexity of bilateral filtering is O(r²), far greater than the O(r) of the SSR.

To keep the enhanced image consistent with the original image in the gradient domain, image enhancement was transformed into solving the Poisson equation as shown in formula (13).

When the gradient image of the illumination layer of the original image is known, the enhanced image $\tilde{I}$ can be obtained, after the Poisson equation of formula (13) has been solved. Next, $\tilde{I}$ should be replaced by the original illumination image I of formula (1), in order to produce the enhanced image with preserved details.

4.1 Solving Poisson equation with Green function

The Poisson equation can be solved by many methods, such as the Green function, fast Fourier transform, central difference method, and multigrid method [21]. In this paper, the second Green function method, which is based on Dirichlet boundary condition, is adopted to solve the Poisson equation of formula (13). The Dirichlet boundary condition controls the brightness of the entire planar region. That is, the enhanced image can be obtained by the Green function, as long as the boundary pixel brightness in the region with low brightness or high exposure are mapped to the brightness suitable for human eyes.

Let ϵ denote a planar closed region with σϵ as the boundary, ϵ^* denote the boundary condition, and r denote any point in ϵ. Then, formula (14) satisfies the Dirichlet boundary condition:

$\left\{\begin{array}{c}\Delta \tilde{I}(r)=\operatorname{div} \nabla F(r) \\ \left.\tilde{I}(r)\right|_{\sigma \epsilon}=\tilde{I}^{*}(r)\end{array}\right.$     (16)

According to the second Green function, the solution of formula (16) can be expressed as:

$\tilde{I}(r)=-\int_{\sigma \epsilon} \tilde{I}^{*}\left(r_{0}\right) \frac{\partial G\left(r, r_{0}\right)}{\partial n_{0}} d s$

$-\iint_{\epsilon} \operatorname{div} \nabla F(r) \times G\left(r, r_{0}\right) d \sigma$      (17)

where, n₀ is the outer normal direction of σϵ; r₀ is the position of point source charge in the planar region ϵ; G(r,r₀) is the guarantee that the Green function satisfies the point source function and the boundary homogeneity:

$\left\{\begin{array}{c}\left.G\left(r, r_{0}\right)\right|_{\sigma \epsilon}=0 \\ \Delta G\left(r, r_{0}\right)=-\vartheta\left(r-r_{0}\right)\end{array}\right.$     (18)

where, $\tilde{I}^{*}\left(r_{0}\right)$ and div∇F(r) are known; the Green function G is the only unknown. Formula (13) can be solved if the Green function satisfies the point source function and the boundary homogeneity.

The Green function is usually constructed by the electric image method, which can provide the Green function of a regular and symmetrical region. For a planar circular region ϵ, it is assumed that the dielectric constant ε equals 1, a standard punctual charge with q is placed at any point in ϵ, and the boundary is composed of m non-repeated points that satisfies σϵ∈{b₁,b₂,…,b_m}, ∀b_i∈R². By the principle of the electric image method, the equivalent charge q' of the induced charge formed by the point charge at any boundary point falls on the extension line between the center of the circle and point q, and lies outside the circle.

3.png

Figure 3. Diagram of point charge, boundary point, and equivalent charge

The point charge, boundary point, and equivalent charge is shown in Figure 3, where c is the center of the circle, r₀ is the point charge, r₁ is the equivalent charge, and b_i is any boundary point. These parameters satisfy the following relationship:

$\frac{\left\|c-b_{i}\right\|}{\left\|c-r_{1}\right\|}=\frac{\left\|c-r_{0}\right\|}{\left\|c-b_{i}\right\|}$     (19)

Since ∠r₀cb_i=∠b_icr₁, we have ∆r₀cb_i=∆b_icr₁. Then, the following relationship can be derived:

$\frac{\left\|b_{i}-p\right\|}{\left\|x_{n}-b_{i}\right\|}=\frac{R}{\left\|c-x_{0}\right\|}$      (20)

According to the theory of electromagnetic field, the Green function describes the potential of point charge r₀ in r. For any point in unbounded space, the Green function can be expressed as:

$G\left(r, r_{0}\right)=\frac{q}{2 \pi} \cdot \frac{1}{\ln \left\|r-r_{0}\right\|}$      (21)

For the circular region in Figure 3, the potential of point r is the superposition of the potential generated by r₀ in the unbounded space and the potential generated by the equivalent charge of the boundary. Thus, G(r,r₀) can be defined as:

$G\left(r, r_{0}\right)=\frac{q}{2 \pi} \cdot \frac{1}{\ln \left\|r-r_{0}\right\|}+\frac{q^{\prime}}{2 \pi} \cdot \frac{1}{\ln \left\|r-r_{1}\right\|}$     (22)

Because the Green function satisfies boundary homogeneity, that is, the potential at any boundary point b_i equals zero, the Green function of any boundary b_i must be equal to zero. Thus, we have:

$G\left(b_{i}, r_{0}\right)=\frac{q}{2 \pi} \cdot \frac{1}{\ln \left\|b_{i}-r_{0}\right\|}+\frac{q^{\prime}}{2 \pi} \cdot \frac{1}{\ln \left\|b_{i}-r_{1}\right\|}$     (23)

Sorting out formulas (19) and (20), and substituting the result into formula (23):

$q^{\prime}=-q \frac{\ln R}{\ln \left\|c-r_{\Omega}\right\|}$      (24)

Substituting formula (24 into formula (22), the Green function in circular region can be expressed as:

$G\left(r, r_{0}\right)=\frac{q}{2 \pi} \cdot \frac{1}{\ln \left\|r-r_{0}\right\|}$

$-\frac{q \ln R}{2 \pi \ln \left\|c-r_{0}\right\| l n\left\|r-r_{1}\right\|}$      (25)

Substituting formula (25) into formula (17), the solution of Poisson equation, i.e., the enhanced image, can be obtained.

According to formula (17), solving the gray value of pixels in the region is equivalent to solving the second type curve integral of the Green function for the boundary of the whole region. The complexity of the algorithm is O(mn), with m being the number of boundary points, and n being the total number of pixels in the region. The high complexity is not conducive to real-time processing. Thus, it is necessary to reduce the complexity of the algorithm by the boundary sampling with equal interval.

The boundary sampling with equal interval can be implemented by the following principle: First, the set of boundary points can be expressed as $\sigma \tilde{\epsilon}=\left\{\tilde{b}_{1}, \tilde{b}_{2}, \ldots, \tilde{b}_{s}\right\}$ based on the boundary points sampled with equal interval $\sigma \in$. After sampling, the Euclidean distance from any point x_i in the region to the boundary point $\tilde{b}_{j}$ after sampling is $\left\|x_{i}-\tilde{b}_{j}\right\|$. As can be seen from formula (17), when the point charge is close to the boundary point, the boundary point contributes greatly to the curve integral. The inverse is also true. Hence, a threshold T_τ can be set to measure the distance. The threshold must satisfy:

$\left\|x_{i}-\tilde{b}_{j}\right\|<T_{\tau}$     (26)

Figure 4 shows a set of points sampled from a point charge at any position in the region and the corresponding boundary. It can be seen that the effectiveness of the boundary points being sampled increases with the proximity to the point. This means, when formula (26) is satisfied, the number of boundary pixels is reduced to lower the complexity of the algorithm.

4.png

Figure 4. Boundary sampling diagram

4.2 Adaptive brightness mapping

Formula (17) suggests that it is possible to obtain the gray values of the pixels in the planar region, after the determination of the Dirichlet boundary condition of the region acted by Green function. Thus, the global brightness of the pixels in the planar region is controlled by the gray values of the boundary pixels.

The human eyes have different sensitivities to different brightness. For example, with the decrease of brightness, the local contrast is reduced, and the human eyes are less able to distinguish between details. Then, the image quality will be rather poor. The strong light, which greatly stimulates the human eyes, would submerge the image details.

If the planar region is dark, the boundary pixels should be made brighter; if the planar region is bright, the boundary pixels should be made darker. Then, the boundary conditions need to be brought into the Green function to solve the Poisson equation, and thus improve the visual effect of the human eyes.

Table 1. Piecewise adaptive brightness mapping

Note: $\widehat{K}$ represents brightness (unit: cd/m²).

This paper chooses the piecewise adaptive brightness mapping to map the brightness of boundary pixels to the level suitable for the human eyes. Table 1 shows the procedure of this mapping strategy.

In the gray image, the gray value measures the brightness of the corresponding pixel. Therefore, the piecewise adaptive brightness mapping cannot be directly applied. Instead, the gray value should be firstly converted into the corresponding brightness. Let K_min and K_max be the minimum and maximum brightness values of the display, respectively. Suppose the gray values of pixels change between 0 and 255. Then, the gray value of a pixel can be converted into brightness by:

$K(x, y)=I(x, y) / 255 \times\left(K_{\max }-K_{\min }\right)$     (27)

where, I is the illumination image obtained by fast bilateral filtering; K is the corresponding brightness image. Then, the brightness image needs to be brought into Table 1 to get the adaptive brightness of each pixel. Finally, the brightness should be converted into gray value by formula (26), thereby producing the corresponding gray image.

To sum up, the proposed image sharpness enhancement algorithm is implemented in the following procedure: Firstly, the original image is divided into a group of symmetrical closed planar regions. Next, the Poisson equation is solved according to the Dirichlet boundary condition of the region. Since the electric image method only applies to regular and symmetrical regions, the image is divided into several regular circular sub blocks through secondary partition. The solution of Poisson equation corresponding to each sub block is obtained by the electric image method. Finally, the solutions of all sub blocks are stitched into the enhanced image. The specific steps of the algorithm are as follows:

Step 1. Separate the illumination layer from the detail layer by fast bilateral filtering.

Step 2. Divide the illumination image into several tangent circular regions, according to the size of the original image.

Step 3. Segment the brightness of all boundary points to obtain the first boundary conditions of all regions.

Step 4. Extract a region, and place a point charge q at any point in that region. Apply the boundary sampling algorithm with interval d to the region, and solve the Green function of the region.

Step 5. Substitute the Green function into formula (17) to calculate the gray value of each pixel in the region, and mark the pixel that has been solved.

Step 6. Repeat Steps 4-5 until all partially cut circular regions have been traversed.

Step 7. Divide the remaining pixels into sub blocks by secondary partition.

Step 8. Place a point charge at any point in the unlabeled pixels of a region being divided into sub blocks, perform the boundary sampling algorithm with interval d  on the region, and solve the Green function of the region.

Step 9. Substituting the Green function into formula (17) to calculate the gray value of each unlabeled pixel.

Step 10. Repeat Steps 8-9 until all the sub blocks have been traversed.

Step 11. Bring the enhanced image into Retinex model to get the enhanced image.

To solve the defects of traditional image enhancement algorithm, this paper proposes an image sharpness enhancement algorithm. The proposed algorithm transforms image enhancement into solving Poisson equation in the gradient domain, solves the Green function by the electric image method, and reduces the complexity of the algorithm by the boundary sampling with equal interval. On this basis, the Dirichlet boundary condition was obtained through gray mapping, and the Poisson equation was solved by Green function to get the final enhanced image. Experimental results show that our algorithm can effectively improve the contrast, details, and clarity of the original image.