A Convex Model for Edge-Histogram Specification with Applications to Edge-Preserving Smoothing

Abstract

The goal of edge-histogram specification is to find an image whose edge image has a histogram that matches a given edge-histogram as much as possible. Mignotte has proposed a non-convex model for the problem in 2012. In his work, edge magnitudes of an input image are first modified by histogram specification to match the given edge-histogram. Then, a non-convex model is minimized to find an output image whose edge-histogram matches the modified edge-histogram. The non-convexity of the model hinders the computations and the inclusion of useful constraints such as the dynamic range constraint. In this paper, instead of considering edge magnitudes, we directly consider the image gradients and propose a convex model based on them. Furthermore, we include additional constraints in our model based on different applications. The convexity of our model allows us to compute the output image efficiently using either Alternating Direction Method of Multipliers or Fast Iterative Shrinkage-Thresholding Algorithm. We consider several applications in edge-preserving smoothing including image abstraction, edge extraction, details exaggeration, and documents scan-through removal. Numerical results are given to illustrate that our method successfully produces decent results efficiently.

1. Introduction

Histogram specification is a process where the image histogram is altered such that the histogram of the output image follows a prescribed distribution. It is one of the many important tools in image processing with numerous applications such as image enhancement [1,2,3], segmentation [4,5,6] among many others. To match the histogram of an image I to a target histogram $\mathbf{h}$, one can use the following procedure:

• Compute the probability density function of I, $p\left(r_j\right)=\frac{n_j}{n}$, $j=1,\cdots ,L$, where $r_j$ denotes the intensity values, $n_j$ denotes the number of pixels whose value equals to $r_j$, and n is the total number of pixels in I.
• With the probability density function $p\left(r_j\right)$, compute the cumulative density function $S\left(r_j\right)$ of I.
• With $\mathbf{h}$ serving as a probability density function, compute the cumulative density function $G\left(s_j\right)$.
• For $p=1,2,\cdots ,L$, replace $r_p$ by $s_q$, where $S\left(r_p\right)=G\left(s_q\right)$.

In the typical 8-bit image setting, $L=256$ and $r_p=p-1$, for $p=1,2,\cdots ,256$.

However, the above procedure gives an inexact histogram matching in discrete images, since the number of pixels is far greater than the number of intensity values in discrete images. Therefore, to ensure exact matching of histograms, different algorithms are introduced, see for instance [7,8,9]. The main idea of these algorithms is to obtain a total order of the pixel values before setting pixel values according to the target histogram.

The goal of edge-histogram specification is to find an image whose edge image has a histogram that matches a given edge-histogram as much as possible. An edge-histogram of an image counts the frequencies of different gradient values in an image and plots the distribution of gradient values. Edge-histograms of natural images have been studied in previous literature [10,11]. It is known that the edge-histograms of natural images contain a sharp peak at 0. This phenomenon is due to the abundance of flat regions in natural images. In natural images, flat regions are seen in the sky, water, walls, and many others. The flat regions in the image produce zero gradients and hence the resulting edge-histogram contains a sharp peak at 0. In addition, as stated in [10], the edge-histograms of natural images have a higher kurtosis and heavier tails compared to a Gaussian distribution with the same mean and variance.

Histogram of gradients has also been used in other fields of image processing. Histogram of oriented gradients (HOG) is used as a feature descriptor for object detection [12,13,14]. This technique first counts the frequencies of gradient orientation in terms of a histogram in a local region, and combines the histograms of gradients in different regions to a single vector as output.

Unlike HOG, our focus is on matching edge-histograms, and we do not consider local edge-histograms. To match the edge-histogram, the author in [11] proposed the following non-convex model for the problem. Given a discrete input image I of size m-by-n, let $\mathbf{r}\in {\mathbb{R}}^k$ be a vector storing the pairwise differences

$$r_{s,t}:=|I_s-I_t|,s=1,2,\cdots mn,t\in {\mathcal{N}}_s,$$

(1)

where $I_i$ denotes the value of the i-th pixel of I, and ${\mathcal{N}}_s$ denotes the set of indices of a fixed shape neighborhood of the s-th pixel, for example, the north, east, south, and west neighboring points of the s-th pixel. Here $k=mn|{\mathcal{N}}_s|$ with $|{\mathcal{N}}_s|$ the number of elements in ${\mathcal{N}}_s$. That is, each entry of the vector $\mathbf{r}$ contains the absolute difference between the values of two neighboring pixels in I. Note that the order of storing the pairwise differences does not matter, as it does not affect the output of histogram specification.

Given a target edge-histogram $\mathbf{h}$, one can apply histogram specification on $\mathbf{r}$ to match $\mathbf{h}$. Denote the output by $\mathbf{d}$:

$$\mathbf{r}=(r_{s_1,t_1},\cdots ,r_{s_k,t_k})\underset{\mathrm{specification}}{\overset{\mathrm{histogram}}{\to }}(d_{s_1,t_1},\cdots ,d_{s_k,t_k})=\mathbf{d}.$$

(2)

To ensure exact histogram matching, the author of [11] used the algorithm in [7]. After applying histogram specification as in (2), one gets the output image X by solving the minimization problem

$${\displaystyle \underset{X}{min}\sum _{s=1}^{mn}\sum _{t\in {\mathcal{N}}_s}{\left({(X_s-X_t)}^2-d_{s,t}^2\right)}^2.}$$

(3)

Model (3) is solved by a conjugate gradient procedure followed by a stochastic local search, please refer to [7,11] for the full details of the algorithm. The non-convex nature of the model hinders the computations, and it is difficult to include additional constraints. In this paper, we propose a convex model that can include additional constraints based on different applications in edge-preserving smoothing.

Edge-preserving smoothing is a popular topic in image processing and computer graphics. Its aim is to suppress insignificant details and keep important edges intact. As an example, the input image in Figure 3a contains textures on the slate and the goal of edge-preserving smoothing is to remove such textures and keep only the object boundaries as in Figure 3c. Numerous methods have been introduced to perform the task. Anisotropic diffusion [15,16] performs smoothing by means of solving a non-linear partial differential equation. Bilateral filtering [17,18,19,20] is a method combining domain filters and range filters. They are widely used because of their simplicity. Optimization frameworks such as the weighted least squares (WLS) [21] and TV regularization [22,23] are also introduced. In WLS, a regularization term is added to minimize the horizontal and vertical gradients with corresponding smoothness weights. Recently, models based on $l_0$-gradient minimization [24,25,26,27,28,29] have become popular. These models focus on the $l_0$-norm of the image gradients.

One application of edge-preserving smoothing is scan-through removal. Written or printed documents are usually subjected to various degradations. In particular, two-sided documents can be suffered from the effect of back-to-front interference, known as “see-through”, see Figure 10. The problem is especially severe in old documents, which is caused by the bad quality of the paper or ink-bleeding. These effects greatly reduce the readability and hinder optical character recognition. Therefore it is of great importance to remove such interference. However, physical restoration is difficult as it may damage the original contents of the documents, which is clearly undesirable as the contents may be important. Consequently, different approaches in the field of image processing are considered to restore the images digitally.

These approaches can be mainly classified into two classes: Blind and Non-blinds methods. Non-blind methods [30,31,32,33,34,35,36,37,38] require the information of both sides. These methods usually consist of two steps. First, the two sides of the images are registered. Then, the output image is computed based on the registered images. It is obvious that these methods strongly depend on the quality of registration; therefore, highly accurate registration is needed. However, perfect registration is hard to achieve in practice due to numerous sources of errors including local distortions and scanning errors. Furthermore, information from the back page is not available on some occasions. Therefore, blind methods which do not assume the availability of the back page are also developed in solving the problem, see [39,40,41,42,43,44].

In this paper, we propose a convex model for applications in edge-preserving smoothing. In our work, we modify the objective function in the non-convex model in [11] so that we only need to solve a convex minimization problem to obtain the output. The simplicity of our model allows us to incorporate different useful constraints such as the dynamic range constraint; the convexity of our model allows us to compute the output efficiently by Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) [45] or Alternating Direction Method of Multipliers (ADMM) [46,47]. We introduce different edge-histograms and suitable constraints in our model, and apply them to different imaging tasks in edge-preserving smoothing, including image abstraction, edge extraction, details exaggeration, and scan-through removal.

The contribution of this paper is twofold. First, we propose a convex variant of the model in [11] which can be solved efficiently. In addition, our model can easily include additional constraints based on different applications. Second, we demonstrate applications of edge-histogram specification in edge-preserving smoothing.

The outline of the paper is as follows: Section 2 describes the proposed convex model, Section 3 presents the applications of our model with numerical results, and conclusions are then presented in Section 4.

2. Our Model

In our model, we do not consider the edge magnitudes as in (1). Instead, we directly consider the image gradients and define $\mathbf{r}$ with entries

$$r_{s,t}:=I_s-I_t.$$

Similar to [11], our model consists of two parts. First, given $\mathbf{r}$, we apply histogram specification on $\mathbf{r}$ to obtain $\mathbf{d}$ as in (2). In the second part, we solve a minimization problem to obtain the output X. Instead of solving the non-convex model (3), we propose a convex model. In Section 2.1, we present our convex model and its solvers. In Section 2.2, Section 2.3 and Section 2.4, we apply our model to specific applications in edge-preserving smoothing.

In the following discussions, we consider only grayscale images. For colored images, we apply our method to R, G, B channels separately.

2.1. Proposed Convex Model

Instead of (3), we propose the convex model

$${\displaystyle \underset{X}{min}\sum _{s=1}^{mn}\sum _{t\in {\mathcal{N}}_s}{\left|(X_s-X_t)-d_{s,t}\right|}^p+{\iota }_{\mathbf{C}}\left(X\right),}$$

(4)

where $p=1$ or 2, $\mathbf{C}$ is a convex set to be discussed in Section 2.4, and ${\iota }_{\mathbf{C}}$ denotes the indicator function of $\mathbf{C}$. The choices of p and $\mathbf{C}$ depend on the applications. Let $\mathbf{x}$ be a vector such that its s-th entry is $X_s$. Then we can rewrite (4) as

$$\underset{\mathbf{x}}{min}\left|\right|G\mathbf{x}-{\mathbf{d}\left|\right|}_p^p+{\iota }_{\mathbf{C}}\left(\mathbf{x}\right).$$

(5)

In our tests, we use ${\mathcal{N}}_s=\{s_v,s_h\}$, where $s_v$ and $s_h$ denote the indices of the north of north and west neighboring points of the s-th pixel respectively. Hence we can write $G={(G_h,G_v)}^T$, where $G_h$, $G_v$ are the horizontal and vertical backward difference operators. We use periodic boundary condition for pixels outside the boundaries, see ([48], p. 258).

For $p=2$, model (5) can be solved by FISTA [45]. For $p=1$, we rewrite (5) as

$$\begin{array}{c}\underset{\mathbf{x},\mathbf{y}}{min}{\parallel \mathbf{y}-\mathbf{d}\parallel }_1+{\iota }_{\mathbf{C}}\left(\mathbf{x}\right)\hfill \\ \mathrm{s}.\mathrm{t}.G\mathbf{x}=\mathbf{y},\hfill \end{array}$$

which can be solved by ADMM [46,47].

2.2. Construction of Target Edge-Histogram

For the applications we considered in this paper, one objective is to remove textures in the images where their edges have a small magnitude. As an example, the textures on the slate in Figure 1a produce smaller edge magnitudes compared to the boundaries of the slate and the letters on the slate, see Figure 1b,c. To eliminate those textures, we could set the values of edges with a small magnitude to zero. Hence, in this paper, we propose to use edge-histograms similar to that shown in Figure 2b as target edge-histogram which is obtained by thresholding the input edge-histogram in Figure 2a. In particular, the target edge-histogram is dependent on the input image. We remark here that it is not uncommon to construct the target histogram based on the input. For example, such construction is used in image segmentation [6].

Since we are just thresholding the edges with small values to zero, the edge-histogram specification (2) can be done easily as follows. Given any input image Y, we first compute its gradients $y_{s,t}=Y_s-Y_t$. Then we set

$$z_{s,t}=\left\{\begin{array}{cc}{y}_{s,t}\hfill & \mathrm{if}|y_{s,t}|\ge \lambda ,\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(6)

The thresholded $z_{s,t}$, where its histogram is shown in Figure 2b, will be used as the vector $\mathbf{d}$ in (5) to obtain the output $\mathbf{x}$. It is obvious that different $\lambda$ gives different outputs, see Figure 3. We see that the smoothness of the output increases with $\lambda$.

2.3. Gaussian Smoothing and Iterations

Strong textures produce edges with large magnitude, which cannot be eliminated using a thresholded edge-histogram as in Figure 2b. To suppress them, the input image I will first pass through a Gaussian filter with standard deviation $\sigma$ to get the initial guess $X^{\left(0\right)}$. Larger $\sigma$ will have a greater effect in suppressing strong textures, but at the same time blur the image. Hence, $\sigma$ should be chosen small enough so that the Gaussian-filtered image is visually equal to I. In our tests, $\sigma$ is chosen to be less than 1. However, there is not an automatic way to compute the optimal $\sigma$. Hence in this model, we leave $\sigma$ as a parameter to be chosen manually by users. Let $X^{\left(0\right)}$ be the Gaussian-filtered image. Whenever such suppression is unnecessary, we set $\sigma =0$ and hence $X^{\left(0\right)}=I$.

As mentioned in Section 2.2, one of our objectives is to map weak edges to zero. This can be done by changing the $\lambda$ in the thresholded edge-histogram or by solving (5) repeatedly. More specifically, given $X^{\left(0\right)}$, we construct $\mathbf{d}$ using (6) and solve (5) to obtain $X^{\left(1\right)}$. Then we repeat the process to obtain $X^{\left(2\right)}$ and so on. Figure 4 shows a comparison; while we see in Figure 4b that the result after one iteration still contains textures in the grasses, almost all of them are removed after three iterations, see Figure 4c.

2.4. Convex Set $\mathbf{C}$

The model (3) does not consider the dynamic range constraint

$$I_s\in [0,255],\forall s=1,2,\cdots ,mn.$$

(7)

For example, consider the case when one defines $\mathbf{h}$ such that every pixel of $\mathbf{r}$ is doubled. In the absence of (7), it is easy to get an exact solution X of (3) if any one of the pixel values is given. However, there is no guarantee that the pixel values of X lie within $[0,255]$. Therefore, when X is converted back to the desired dynamic range, either by stretching or clipping, the edge-histogram is no longer preserved. To avoid this, it is better to include the dynamic range constraint in the objective function. Therefore, we use the following constraint in all our applications:

$$\mathbf{C}=\{\mathbf{x}:x_i\in [0,255],\forall i\}.$$

(8)

In scan-through removal, we assume the background in books and articles have a lighter intensity than the ink in all color channels. Therefore, in addition to the dynamic range constraint, we also keep the value of the background pixels unchanged. Hence, we set

$$\mathbf{C}=\{\mathbf{x}:x_i\in [0,255],\forall i\mathrm{and}{x}_i=X_i^{\left(0\right)}\mathrm{if}{x}_i\ge \alpha \},$$

(9)

where $\alpha$ is the approximate intensity of the background to be defined in Section 3.4.

3. Applications and Comparisons

Edge-preserving smoothing includes many different applications. In this section, we consider four applications, namely image abstraction, edge extraction, details exaggeration, and scan-through removal. For the first three applications, we use $p=2$ in (5) and solve it by FISTA with the input image as an initial guess. We compare with four existing methods: bilateral filtering [18], weighted-least square [21], $l_0$-smoothing [24], and $l_0$-projection [29]. For the scan-through removal, we use $p=1$ in (5) and solve it by ADMM with the input image and $\mathbf{d}$ as the initial guesses. We compare with one blind method [43] and three non-blind methods [32,44,49]. In all applications, the number of iterations is fixed at 3. The values $\lambda$ and $\sigma$ vary for different images and will be stated separately.

For the tests below, we select the parameters which give the output image with the best visual quality. Some of the comparison results are obtained directly from the authors’ work and some are done by ourselves. For the results done by ourselves, we list out the parameters we have used.

3.1. Image Abstraction

The goal of image abstraction is to remove textures and fine details so that the output looks un-photorealistic. This can be done by solving (5) with constraint (8). As shown in Figure 5, the textures of the objects in the photorealistic input image in Figure 5a is removed and our output in Figure 5f becomes un-photorealistic. We see that our model successfully eliminates almost all object textures and keeps the object boundaries intact. As we see in Figure 5f, the details in the basketball net in our output are kept intact, while it disappears, or almost disappears, in the outputs of other models.

3.2. Edge Extraction

Object textures are sometimes misclassified as edges during the edge detection process. In order to reduce misclassifications, image abstraction as discussed in the last section can be used to suppress object textures. Given an input image as shown in Figure 6a, objects of less importance such as clouds and grasses can be eliminated by image abstraction. Using our method, a smooth image as in Figure 6f is obtained. Edge detection or segmentation can then be applied to the output image to obtain a result with much fewer distortions. Figure 7 shows the results of applying the Canny edge detector to the grayscale version of Figure 6. We see that while the outputs of other models keep unnecessary details, our model produces a result containing only salient edges, and removes unimportant details.

3.3. Details Exaggeration

Details exaggeration is to enhance the fine details in an image as much as possible. Given an input image I, we obtain a smooth image X by our method where the textures in I are removed, see in Figure 8b. As seen in Figure 8c, the image $|I-X|$ has small values in regions with insignificant textures and large values in the parts containing strong textures. By enhancing ($I-X$) and adding it back to X, a details-exaggerated image J can be obtained, see Figure 8d. Mathematically, we have $J=X+s(I-X)$, where $s>1$ is a parameter controlling the extent of exaggeration. Figure 9 shows a comparison with the results by other methods. In Figure 9f we see that our model successfully produces a better result with more exaggerated details.

3.4. Scan-Through Removal

Two-sided documents can be suffered from the effect of back-to-front interference, known as “see-through”. Usually, “see-through” produces relatively small gradient fluctuations than the main content we want to preserve, see Figure 10. By considering the edges, one can identify interferences and eliminate them.

Recall in (9), we also impose a constraint that background pixels will not be modified. Here background pixels refer to the pixels with values not less than $\alpha$. To find a suitable $\alpha$, we first need to locate background regions—regions which contain only insignificant intensity change, i.e., the standard deviation of the intensity within the region should be small. Motivated by this, we design a multi-scale sliding window method to compute a suitable $\alpha$. A sliding window with size w is used to scan through an input image Y with stride $\lceil w/5\rceil$. At each location p, the mean intensity $m_p$ of the sliding window is computed and if its standard deviation ${\sigma }_p$ is smaller than a parameter $\widehat{\sigma }$, $m_p$ will be stored for future selection. After scanning through the whole image, we set $\alpha$ to the largest stored value to avoid choosing regions with purely foreground or interference. If ${\sigma }_p\ge \widehat{\sigma }$ for all p, we replace w by $w/2$ and scan through the image again. At the worst case when $w=1$, it is equivalent to setting $\alpha$ to the maximum intensity of the image. In our tests, we use $\widehat{\sigma }=3$.

The reason for using a varying window size is that a small window will have a chance of capturing extreme values and a large window will have a chance of failing in capturing pure background. Therefore, we start from a large window and stop once we find at least one region with a small standard deviation. The initial window we use is the largest square window of length $w=2^{\ell }$ that can fit inside the given image. Figure 10 shows the windows (red-colored squares) obtained by the procedure above. We see that it successfully locates a pure background. The background level $\alpha$ is the mean intensity of the corresponding square.

With $\alpha$ found, we solve our model (5) with constraint (9) to obtain the output. We test our method using the first image in Figure 10. Our output is shown in Figure 11b, where we see that the contents are kept and the back-page interferences are removed. Figure 11a shows a comparison with the blind method from [43]. We also compare our result with three non-blind methods [32,44,49]. For copyright reasons, we can only refer readers to the papers [32,44] to see the resulting images from the three methods. Our method outperforms the blind method and is comparable to the non-blind methods, while these non-blind methods require information from both sides.

4. Conclusions

We have proposed a convex model with suitable constraints for edge-preserving smoothing tasks including image abstraction, edge extraction, details exaggeration, and documents scan-through removal. Our convex model allows us to solve it efficiently by existing algorithms.

In this paper, because of the special applications we considered, we use only the thresholded histograms as target edge-histograms. In the future, we would investigate more general shapes of edge-histograms and apply them to a wider class of problems.