Semi-Global Weighted Least Squares in Image Filtering

TL;DR

This paper introduces Semi-Global Weighted Least Squares (SG-WLS), an efficient approximation method for image filtering that significantly reduces computational costs while maintaining performance close to the original WLS model.

Contribution

We propose SG-WLS, a novel iterative approach using one-dimensional WLS subsystems with a special neighborhood construction, enabling faster and more memory-efficient image filtering.

Findings

SG-WLS is approximately 20 times faster than the original WLS model.

Memory usage of SG-WLS is significantly lower, at most a small fraction of WLS.

SG-WLS achieves comparable filtering performance to the full 2D WLS model.

Abstract

Solving the global method of Weighted Least Squares (WLS) model in image filtering is both time- and memory-consuming. In this paper, we present an alternative approximation in a time- and memory- efficient manner which is denoted as Semi-Global Weighed Least Squares (SG-WLS). Instead of solving a large linear system, we propose to iteratively solve a sequence of subsystems which are one-dimensional WLS models. Although each subsystem is one-dimensional, it can take two-dimensional neighborhood information into account due to the proposed special neighborhood construction. We show such a desirable property makes our SG-WLS achieve close performance to the original two-dimensional WLS model but with much less time and memory cost. While previous related methods mainly focus on the 4-connected/8-connected neighborhood system, our SG-WLS can handle a more general and larger neighborhood system thanks to the proposed fast solution. We show such a generalization can achieve better performance than the 4-connected/8-connected neighborhood system in some applications. Our SG-WLS is $\sim20$ times faster than the WLS model. For an image of $M\times N$, the memory cost of SG-WLS is at most at the magnitude of $max\{\frac{1}{M}, \frac{1}{N}\}$ of that of the WLS model. We show the effectiveness and efficiency of our SG-WLS in a range of applications. The code is publicly available at: https://github.com/wliusjtu/Semi-Global-Weighted-Least-Squares-in-Image-Filtering.

Tables

Table 1. Table 1: Computation time of different methods for smoothing a 1-megapixel RGB image.

Method	WLS [6]	FGS [20]	Ours ($r=1,\tau =1$)
Time (seconds)	8.02	0.1	0.35

Table 2. Table 2: Mean Absolute Difference (MAD) of guided depth upsampling errors for different methods. Best results are in bold.

		Art			Book			Dolls			Laundry			Moebius			Reindeer
	$2\times$	$4\times$	$8\times$	$2\times$	$4\times$	$8\times$	$2\times$	$4\times$	$8\times$	$2\times$	$4\times$	$8\times$	$2\times$	$4\times$	$8\times$	$2\times$	$4\times$	$8\times$
JBU [11]	1.59	2.06	3.18	0.87	1.24	2.04	0.91	1.2	1.98	0.94	1.38	2.15	0.89	1.28	2.05	0.95	1.36	2.24
GF [10]	1.91	2.23	3.11	0.84	1.19	1.86	0.87	1.17	1.89	1.01	1.31	2.25	0.92	1.19	1.88	1.06	1.34	1.98
WLS [6]	1.58	2.52	3.96	0.9	1.25	1.85	0.93	1.3	1.84	1.03	1.5	2.33	0.94	1.34	1.97	1.09	1.6	2.42
FGS [20]	1.36	2.01	3.71	1.25	1.73	2.58	1.33	1.89	2.62	1.11	1.62	2.61	1.35	1.95	2.81	1.46	2.08	3.07
Ours ($r=1,\tau =1$)	1.6	2.4	3.75	0.95	1.31	1.97	0.96	1.37	2.05	1.05	1.6	2.45	0.94	1.37	2.02	1.12	1.68	2.49
Ours ($r=4,\tau =4$)	1.26	1.9	3.07	0.82	1.12	1.73	0.87	1.11	1.81	0.86	1.17	2	0.82	1.08	1.79	0.9	1.32	2.01

Equations

$$E(U)=\sum _{i\in \Omega }(U_i-F_i{)}^2+\lambda \sum _{i\in \Omega }\sum _{j\in N(i)}{\omega }_{i,j}(U_i-U_j{)}^2$$

$$\begin{array}{rl}& {\omega }_{i,j}={\omega }_{i,j}^{frac}=\frac{1}{|i-j{|}^{{\alpha }_s}+\epsilon }\cdot \frac{1}{|G_i-G_j{|}^{{\alpha }_r}+\epsilon }\\ & {\omega }_{i,j}={\omega }_{i,j}^{exp}=\exp \left(-\frac{|i-j{|}^2}{2{\sigma }_s^2}\right)\cdot \exp \left(-\frac{|G_i-G_j{|}^2}{2{\sigma }_r^2}\right)\end{array}$$

$$A\cdot U=F$$

$$\displaystyle\small A_{i,j}=\left\{\begin{array}[]{l}1+\lambda\sum_{j\in N(i)}\omega_{i,j}\ \ \ \text{for}\ i=j\\ -\lambda\omega_{i,j}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \text{for}\ j\in N(i)\\ 0\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \text{otherwise}\end{array}\right.$$

$$A_s\cdot u=f$$

$$A_s=\left[\begin{array}{c}\begin{array}{cccccccccc}{a}_1& b_{1,1}& b_{1,2}& \cdots & b_{1,r}& & & & & \\ c_{1,1}& a_2& b_{2,1}& \cdots & b_{2,r-1}& b_{2,r}& & & & \\ c_{1,2}& c_{2,1}& a_3& \cdots & b_{3,r-2}& b_{3,r-1}& b_{3,r}& & & \\ ⋮& ⋮& ⋮& \ddots & \ddots & \ddots & \ddots & \ddots & & \\ c_{1,r}& c_{2,r-1}& c_{3,r-2}& \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \\ & c_{2,r}& c_{3,r-1}& \ddots & \ddots & a_{s-r}& b_{s-r,1}& b_{s-r,2}& \cdots & b_{s-r,r}\\ & & c_{3,r}& \ddots & \ddots & b_{s-r,1}& \ddots & \ddots & \ddots & ⋮\\ & & & \ddots & \ddots & b_{s-r,2}& \ddots & a_{s-2}& b_{s-2,1}& b_{s-2,2}\\ & & & & \ddots & ⋮& \ddots & c_{s-2,1}& a_{s-1}& b_{s-1,1}\\ & & & & & c_{s-r,r}& \cdots & c_{s-2,2}& c_{s-1,1}& a_s\end{array}\end{array}\right]$$

$$A_s=P\cdot Q$$

$$P\cdot y=f$$

$$Q\cdot u=y$$

$$P=\left[\begin{array}{c}\text{setcounter}MaxMatrixCols11\begin{array}{ccccccccccc}{\alpha }_1& & & & & & & & & & \\ {\gamma }_{1,1}& {\alpha }_2& & & & & & & & & \\ {\gamma }_{1,2}& {\gamma }_{2,1}& {\alpha }_3& & & & & & & & \\ ⋮& ⋮& ⋮& \ddots & & & & & & & \\ {\gamma }_{1,r}& {\gamma }_{2,r-1}& {\gamma }_{3,r-2}& \ddots & \ddots & & & & & & \\ & {\gamma }_{2,r}& {\gamma }_{3,r-1}& \ddots & \ddots & {\alpha }_{s-r}& & & & & \\ & & {\gamma }_{3,r}& \ddots & \ddots & {\gamma }_{s-r,1}& \ddots & & & & \\ & & & \ddots & \ddots & ⋮& \ddots & \ddots & & & \\ & & & & \ddots & ⋮& \ddots & \ddots & {\alpha }_{s-1}& & \\ & & & & & {\gamma }_{s-r,r}& \cdots & {\gamma }_{s-2,2}& {\gamma }_{s-1,1}& {\alpha }_s& \end{array}\end{array}\right]$$

$$Q=\left[\begin{array}{c}\begin{array}{ccccccccc}1& {\beta }_{1,1}& {\beta }_{1,2}& \cdots & {\beta }_{1,r}& & & & \\ & 1& {\beta }_{2,1}& \cdots & {\beta }_{2,r-1}& {\beta }_{1,r}& & & \\ & & 1& \cdots & {\beta }_{3,r-2}& {\beta }_{2,r-1}& {\beta }_{1,r}& & \\ & & & \ddots & \ddots & \ddots & \ddots & \ddots & \\ & & & & 1& {\beta }_{s-r,1}& \cdots & \cdots & {\beta }_{s-r,r}\\ & & & & & \ddots & \ddots & \ddots & ⋮\\ & & & & & & 1& {\beta }_{s-2,1}& {\beta }_{s-2,2}\\ & & & & & & & 1& {\beta }_{s-1,1}\\ & & & & & & & & 1\end{array}\end{array}\right]$$

$$\small\begin{array}[]{l}\alpha_{1}=a_{1},\ \gamma_{1,i}=c_{1,i},\ \beta_{1,i}=b_{1,i}/\alpha_{1}\end{array}$$

$$\displaystyle\small\begin{array}[]{l}\alpha_{k}\ \ =a_{k}-\sum_{t=1}^{k-1}\gamma_{k-t,t}\beta_{k-t,t},\\ \gamma_{k,i}=c_{k,i}-\sum_{t=1}^{k-1}\gamma_{k-t,i+t}\beta_{k-t,t},\\ \beta_{k,i}=\frac{1}{\alpha_{k}}(b_{k,i}-\sum_{t=1}^{k-1}\beta_{k-t,i+t}\gamma_{k-t,t})\end{array}$$

$$\displaystyle\small\begin{array}[]{l}\gamma_{k,i}=c_{k,i}-\sum_{t=1}^{r-i}\gamma_{k-t,i+t}\beta_{k-t,t},\\ \beta_{k,i}=\frac{1}{\alpha_{k}}(b_{k,i}-\sum_{t=1}^{r-i}\beta_{k-t,i+t}\gamma_{k-t,t})\end{array}$$

$$\displaystyle\small\begin{array}[]{l}\alpha_{k}\ \ =a_{k}-\sum_{t=1}^{r}\gamma_{k-t,t}\beta_{k-t,t},\\ \end{array}$$

$$\displaystyle\small\begin{array}[]{l}\gamma_{k,i}=c_{k,i}-\sum_{t=1}^{r-i}\gamma_{k-t,i+t}\beta_{k-t,t},\\ \beta_{k,i}=\frac{1}{\alpha_{k}}(b_{k,i}-\sum_{t=1}^{r-i}\beta_{k-t,i+t}\gamma_{k-t,t})\end{array}$$

$$\displaystyle\small\begin{array}[]{l}\gamma_{k,r}=c_{k,r},\beta_{k,r}=\frac{b_{k,r}}{\alpha_{k}}\end{array}$$

$$\small\begin{array}[]{l}y_{1}=\frac{f_{1}}{\alpha_{1}}\end{array}$$

$$\displaystyle\small\begin{array}[]{l}y_{k}=\frac{1}{\alpha_{k}}(f_{k}-\sum_{t=1}^{k-1}\gamma_{t,k-t}y_{t})\end{array}$$

$$\displaystyle\small\begin{array}[]{l}y_{k}=\frac{1}{\alpha_{k}}(f_{k}-\sum_{t=1}^{r}\gamma_{k-t,t}y_{k-t})\end{array}$$

$$\small\begin{array}[]{l}u_{s}=y_{s}\end{array}$$

$$\small\begin{array}[]{l}u_{k}=y_{k}-\sum_{t=1}^{s-k}\beta_{k,t}u_{k+t}\end{array}$$

$$\small\begin{array}[]{l}u_{k}=y_{k}-\sum_{t=1}^{r}\beta_{k,t}u_{k+t}\end{array}$$

$$U(m)=\frac{A^{-1}F(m)}{A^{-1}H(m)},m\in \Omega$$

Full text

1 Introduction

Image smoothing is an important operation in both image processing and computer graphics. Many applications require decomposing an image into a piecewise smooth base layer which contains the main structure information and a detail layer which captures the residual details. To achieve the decomposition, Edge-Preserving Smoothing (EPS) is required. EPS can be achieved with local filters which compute the output as a weighted average of the input. Bilateral filter [24] is one of the well-known filters which has been widely used in various applications such as image upsampling [11], flash/no flash image filtering [21] and HDR tone mapping [5]. There are also other local filters based on different theories and computational models [8, 10, 19, 29]. Most local filters can be efficiently computed. However, they can cause gradient reversals and halo artifacts [6, 10] which are their main drawbacks.

There are also EPS methods based on global filters [6, 17, 26, 27]. These methods seek a globally optimal solution to the objective function. The objective function usually contains a data constraint term and a prior smoothness term and needs to be solved globally in a principled manner. Global filters can overcome the limitations of local filters such as gradient reversals and halo artifacts [6]. However, most of global filters are time-consuming and some of them are also memory-consuming. The Weight Least Squares (WLS) model has been widely used in image processing and computer graphics [1, 4, 6, 17]. The solution to the model typically needs to solve a large linear system which is both time- and memory-consuming.

In this paper, we present a new approximation to the WLS model which is denoted as Semi-Global Weighted Least Squares (SG-WLS). The main contributions of this paper are as follows.

$-$ Our SG-WLS can achieve close performance to the original WLS model in several challenging tasks. Yet, our SG-WLS is $\sim 20$ times faster than the original WLS model. For an image of $M\times N$, the memory cost of SG-WLS is at most at the magnitude of $\max\{\frac{1}{M},\frac{1}{N}\}$ of that of the original WLS model.

$-$ The 1D filters in our SG-WLS can take two-dimensional neighborhood information into account each time due to the newly proposed neighborhood construction. This is different from previous methods [8, 22, 29, 20] that can only consider neighbors in one dimension each time. Our neighborhood construction thus overcomes the limitation of previous methods and shows improved results in several applications.

$-$ We propose a fast solution to each 1D filter in our SG-WLS, which is denoted as r-band LU decomposition, to handle a more general and larger neighborhood system other than the 4-connected/8-connected neighborhood system. Such a generalization can achieve better performance than the 4-connected/8-connected neighborhood system in some applications such as guided depth upsampling.

Background The WLS model is a fundamental optimization framework that has been widely used [1, 4, 6, 17]. Given a target image $F$ to be filtered and a guidance image $G$, the formulation of WLS is defined as:

[TABLE]

where $\Omega$ represents the set of all the coordinates. $\lambda$ is a parameter that balances the data term and smoothness term. A larger $\lambda$ results in a larger smoothing effect on $F$. $N(i)$ is the neighborhood of the pixel with coordinate $i$ which is a square patch of $(2r+1)\times(2r+1)$ centered at $i$. $\omega_{i,j}$ is the guidance weight based on $G$. Based on different applications which we will detail in Sec. 3, we adopt two types of guidance weight for different applications:

[TABLE]

where $\alpha_{s}$, $\alpha_{r}$, $\sigma_{s}$ and $\sigma_{r}$ are constants defined by the user. $\varepsilon$ is a small constant that prevents division by zero in areas where $G$ is constant. In this paper, we set $\varepsilon=0.0001$.

The unique minimum of Eq. (1) can be obtained by solving the following large linear system:

[TABLE]

here $U$ and $F$ are vector representations of two-dimensional images. If an image is of size $M\times N$, then its vector representation is a $S\times 1$ vector where $S=M\times N$. $A$ is a $S\times S$ matrix and is defined as:

[TABLE]

There are two challenging issues in solving Eq. (3): (I) Large memory cost. Despite the memory cost of the solver, the memory cost of storing $A$ in Eq. (7) is $\mathcal{O}(MNr^{2})$. As $r$ becomes larger, the memory cost increases very fast. (II) Large time cost. Solving Eq. (3) is an inverse problem which is time-consuming since Eq. (3) is very large. Note that the time cost also increases as $r$ becomes larger.

Eq. (3) can be directly solved by modern linear solvers such as Preconditioned Conjugate Gradient (PCG) [12, 13]. However, the convergence of PCG strongly depends on a good choice of the preconditioner [17]. Besides, the memory cost is not reduced. To solve the above two issues, various approximate methods have been proposed in the literature. Barron et al. [3] proposed to first project the original image into a bilateral space. Then a much smaller linear system is solved with PCG [23]. The output is projected back into the original image space which is the final filtered result. The method in [2] shares the similar idea. Both of their methods are post-processed by the domain transform filtering [8] to smooth out the blocky artifacts introduced by the bilateral grid. Xu et al. [25] proposed to first cluster the image in an affinity space with kd tree and then solve another linear system with the clustered image. In summary, the key idea of these methods is to reduce the dimension of either the large matrix ($``A"$ in Eq. (3)) or the input image. Thus the inverse operation of the very large system can be reduced to the matrix inverse of much smaller matrixes which can reduce both time and memory cost.

The Fast Global Smoother (FGS) proposed by Min et al. [20] is closely related to our work. FGS divides the WLS model into a sequence of subsystems in each row and column which is different from the previous methods [1, 2, 3, 25]. As FGS is locally global in each row and column, we denote it as a semi-global method. Our work is also a semi-global method. It is different from FGS [20] in the following two aspects: (I) Our construction of subsystems is different from that in FGS. FGS only considers neighbors in one dimension each time (row direction or column direction) which has largely destroyed the two-dimensional neighborhood system in the original model. This can cause noticeable artifacts as illustrated in Fig. 1(c). On the contrary, our method can take two-dimensional neighborhood information into account. Thus, our method can overcome the limitation of FGS and achieves close performance to the original WLS model as illustrated in Fig. 1(d). (II) The FGS only adopts the 4-connected/8-connected neighborhood system while our method can handle a more general $(2r+1)\times(2r+1)$ neighborhood system thanks to the proposed fast solution. We show that such a generalization makes our method achieve better performance than the 4-connected/8-connected one in some applications such as guided depth upsampling [7, 16].

2 Semi-Global Weighted Least Squares

2.1 Neighborhood Construction of Subsystems

Like many previous separate filters [8, 22, 29], the FGS proposed by Min et al. [20] separates a two-dimensional filtering process into an iterative one-dimensional filtering process. FGS is fast and can process a 1-megapixel RGB image in $0.1$ seconds on a standard desktop. However, the main limitation of this separation is that only neighborhood in one dimension is considered each time. In such cases, the original two-dimensional neighbors has been largely destroyed as illustrated in Fig. 3. In some applications, such limitation can cause noticeable artifacts. We show some examples of image detail enhancement and HDR tone mapping obtained by FGS in Fig. 1(c) and Fig. 2(b) respectively.

In this section, we show a new type of neighborhood construction that can handle neighbors in two dimensions each time within one-dimensional separate filters. For a pixel at row $s$ and column $t$ in image $I$ of size $M\times N$, we denote it as $I_{(s,t)}$. For a WLS model with a $(2r+1)\times(2r+1)$ neighborhood system, we first extract $2r+1$ columns around the $k$th ($k=r+1,\cdots,N-r$) column in the image, which are denoted as $[I_{(*,k-r)},I_{(*,k-r+1)},\cdots,I_{(*,k)},\cdots,I_{(*,k+r-1)},I_{(*,k+r)}]$. Here $I_{(*,k)}$ denotes all the pixels in the $k$th column. If arranged in their original order, these $2r+1$ columns totally have $M$ row vectors of size $2r+1$. Then for each $i$th ($i=1,\cdots,M$) row vector, if $i$ is even, we reverse the order of the $i$th row vector. Finally, these row vectors are connected head to end to form a $(2r+1)M\times 1$ column vector. We denote this process as neighborhood construction along column direction. It can also be performed along row direction for each $k$th ($k=r+1,\cdots,M-r$) row. In this way, the formed 1D vector is a $(2r+1)N\times 1$ one. An example of $r=1$ is illustrated in Fig. 3.

Now we show how the above neighborhood construction can handle neighbors in two dimensions. First, note that for any $2r+1$ neighboring pixels in the formed 1D vector, they are also neighbors in the original $(2r+1)\times(2r+1)$ neighborhood system. Then take the neighborhood construction along column direction for example, the final 1D vector also contains pixels from neighboring $r$ columns on each side of the current column, which also contains neighbors in the row direction. In this way, the final vector can handle neighbors from both row and column directions. As illustrated in Fig. 3, pixel $``4"$ and pixel $``6"$ are neighbors of pixel $``5"$ in row direction in the original neighborhood system. In the neighborhood construction along column direction, pixel $``4"$ and pixel $``6"$ are still neighbors of pixel $``5"$. This newly designed neighborhood construction can well overcome the limitation of previous separate filters that only handle neighbors in one dimension [8, 22, 29]. Fig. 1(d) and Fig. 2(c) show results obtained with our newly designed neighborhood construction. Our results are indistinguishable to the ones obtained with WLS [6] and well overcome the limitation of FGS [20]. Details of our method will be described in Sec. 2.3.

2.2 Fast and Exact Solution to Subsystems

Since any neighboring $2r+1$ pixels in the 1D vector formed in the neighborhood construction step are neighboring pixels in the original two-dimensional neighborhood system, we can solve another WLS model using this 1D vector with a $2r+1$ neighborhood system. Similarly, for each formed 1D vector, we can obtain a much smaller linear system as follows:

[TABLE]

where $f$ is the formed 1D vector in the neighborhood construction step. $A_{s}$ is formed in a way similar to Eq. (7). Eq. (8) is denoted as a subsystem of our method. Note that the neighborhood system $N(i)$ here is one-dimensional $2r+1$ pixels. To be explicit, $A_{s}$ has the following structure:

[TABLE]

here $s=(2r+1)M$ for neighborhood construction along column direction and $s=(2r+1)N$ for neighborhood construction along row direction.

Matrix $A_{s}$ in Eq. (8) is positive definite and it is also a diagonal matrix with bandwidth of $r$. For $r=1$, solving Eq. (8) has been well studied and has classical solution such as LU decomposition [9]. However, here we need to solve a more general case with $r\geq 1$ which is seldom studied. In this section, we propose a solution to Eq. (8) to handle the case for any $r\geq 1$, which is denoted as $r$-band LU decomposition. Moreover, such a decomposition can be completed efficiently. To be more explicit, this process can be formulated as follows:

[TABLE]

[TABLE]

[TABLE]

Assuming $P$ and $Q$ have the following structures:

[TABLE]

[TABLE]

According to Eq. (9), for $i=1,\cdots,r$, we have:

[TABLE]

For $k=2,\cdots,r$ and $i=1,\cdots,r-k+1$, we have:

[TABLE]

For $k=3,\cdots,r$ and $i=r-k+2,\cdots,r-1$, we have:

[TABLE]

For $k=r+1,\cdots,s$, we have:

[TABLE]

For $k=r+1,\cdots,s-1$ and $i=1,\cdots,\min\{r-1,s-k\}$, we have:

[TABLE]

For $k=2,\cdots,s-r$, we have:

[TABLE]

Eqs. (12)$\sim$(26) are the r-band LU decomposition in Eq. (9). When this is completed, we can solve Eq. (8) through Eq. (10) and Eq. (11). When solving Eq. (10), we first have:

[TABLE]

For $k=2,\cdots,r$, we have:

[TABLE]

For $k=r+1,\cdots,s$, we have:

[TABLE]

When solving Eq. (11), we first have:

[TABLE]

For $k=s-1,\cdots,s-r+1$, we have:

[TABLE]

For $k=s-r,\cdots,1$, we have:

[TABLE]

The subsystem in Eq. (8) can be solved exactly through Eqs. (12)$\sim$(34). Note that for $r=1$, only Eqs. (12), (21), (26), (27), (31), (32) and (34) are needed. For $r=2$, only Eq. (19) is not needed.

2.3 Semi-global Edge-preserving Smoothing

Based on the subsystems described in Sec. 2.1 and Sec. 2.2, we can perform edge-preserving smoothing which can benefit numerous applications [6, 20]. For a $M\times N$ image with a $(2r+1)\times(2r+1)$ two-dimensional neighborhood system, our goal is to divide the original two-dimensional WLS model described in Eq. (1) into a sequence of one-dimensional WLS models with a $2r+1$ one-dimensional neighborhood system in Eq. (8).

Our smoothing process contains four steps. The first step is the one-dimensional neighborhood construction described in Sec. 2.1. $2r+1$ columns (along column direction) or rows (along row direction) in the original image are needed each time. This step results in 1D vectors of size $(2r+1)M$ along column direction or $(2r+1)N$ along row direction. Each vector has a $2r+1$ neighborhood system. The second step is solving the linear system with the formed 1D vector which is described in Sec. 2.2 . The third step is transforming the solution in the second step into an image patch of $2r+1$ rows/columns. This is a simple inverse operation of the neighborhood construction in the first step. The fourth step is averaging pixel values. This is because one pixel can be involved in several subsystems. These values in different subsystems of the same pixel are averaged as the final output. The above four steps are performed $T$ times along column direction and row direction alternatively to get the final smoothed image. In this paper, we find $T=2\sim 4$ is appropriate for most applications. As each subsystem is a globally optimized one while it is only performed in a local region of an image, we call our method as Semi-Global Weighted Least Squares (SG-WLS).

When SG-WLS is used for sparse interpolation such as guided depth upsampling [7, 16] and colorization [15], it cannot be directly applied to the sparse input data due to the unstable result. Instead, we perform guided sparse interpolation in a way similar to the one in [14, 16, 20]. Denote the index map of input as $H$, then the output is computed as:

[TABLE]

The WLS smoothing is applied to both the input $F$ and the index map $H$. This procedure is approximated with our SG-WLS.

Different from the FGS [20] that only adopts the 4-connected/8-connected neighborhood system, our SG-WLS focuses on a more general case where the neighborhood system is $(2r+1)\times(2r+1)$. In fact, the 8-connected neighborhood system in the FGS [20] is a special case of our neighborhood system with $r=1$. In their work, the smoothing strength is increased by enlarging $\lambda$. For our SG-WLS, we show that the smoothing strength can also be increased by enlarging the radius $r$ of the neighborhood system. Particularly, in some cases, enlarging $r$ can achieve the smoothing property that cannot be achieved by simply enlarging $\lambda$. Fig. 4 shows one example of guided depth upsampling. When enlarging $\lambda$, depth edges have been blurred while some parts are still noisy as highlighted. However, when we use a larger $r$ but a smaller $\lambda$, depth edges are well preserved while the noise is also smoothed. Note that Yang et al. [28] also adopted an $11\times 11$ neighborhood system other than a 4-connected/8-connected one for guided depth upsampling. This also validates the effectiveness of larger neighborhood systems. Quantitative measurement of Mean Absolute Difference (MAD) between the result and the groundtruth also shows the effectiveness of using larger neighborhood systems.

The neighborhood construction in the first step can be performed in a sliding manner. Take the neighborhood construction along column direction for example, when a 1D vector is formed with $2r+1$ columns centered at the $k$th column, i.e., $[I_{(*,k-r)},I_{(*,k-r+1)},\cdots,I_{(*,k)},\cdots,I_{(*,k+r-1)},I_{(*,k+r)}]$. Then the next 1D vector is formed with $2r+1$ columns centered at the $(k+1)$th column. However, we find this kind of neighborhood construction is redundant. We can slide it with a step $\tau$, i.e., for the above case, the next 1D vector can formed with $2r+1$ columns centered at the $(k+\tau)$th column. According to our experimental results, for $\tau\leq r$, our SG-WLS has similar performance to that of $\tau=1$ but with smaller computational cost. Fig. 5 shows an example of detail enhancement and guided depth upsampling. The MAD of guided depth upsampling results in Fig. 5(c) also shows increasing $\tau$ seldom decreases the performance. Fig. 6 shows computation time comparison among different values of $\tau$. As shown in Fig. 6, increasing $\tau$ can greatly reduce the computational cost.

Computation Complexity

Assume that the image is of size $M\times N$. The neighborhood system is $(2r+1)\times(2r+1)$ in the original two-dimensional WLS model. Thus, the neighborhood system of our SG-WLS is $2r+1$. Despite the few elements near the boundary, when solving a subsystem in Eq. (8), there are $r$ times multiplication for each $\alpha_{k}$ and $\gamma_{k,i}$, $r+1$ times multiplication for each $\beta_{k,i}$, $r+1$ times addition for $\alpha_{k}$, $\gamma_{k,i}$ and $\beta_{k,i}$. Thus, Eqs. (12)$\sim$(26) require $\mathcal{O}(sr+sr^{2}+sr(r+1))$ times multiplication operations and $\mathcal{O}(s(r+1)+2sr(r+1))$ times addition operations. Since $s=(2r+1)M$ along column direction and $s=(2r+1)N$ along row direction, the computational complexity is $\mathcal{O}(Mr^{3})$ or $\mathcal{O}(Nr^{3})$ for both multiplication and addition. Similarly, solving Eq. (10) and Eq. (11) through Eqs. (27)$\sim$(34) requires $\mathcal{O}(Mr^{2})$ or $\mathcal{O}(Nr^{2})$ multiplication and addition operations. Thus, the computational complexity of solving a subsystem in Eq. (8) is $\mathcal{O}(Mr^{3})$ or $\mathcal{O}(Nr^{3})$. There are total $\frac{N}{\tau}$ subsystems when applying SG-WLS along column direction and $\frac{M}{\tau}$ subsystems along row direction. Thus, when applied $T$ times, the final computational complexity of SG-WLS is $\mathcal{O}(\frac{T}{\tau}MNr^{3})$. Fig. 6 shows the computation time of different $r$ and $\tau$ for filtering a 1-megapixel RGB image on a computer with an Intel i7 3.40GHz CPU and 8GB memory. The iteration number is fixed with $T=4$. In particular, we compare the computation time of our SG-WLS of $r=1,\tau=1$ with that of the WLS model [6] and FGS [20] in Table 1. Our method is over $20\times$ faster than the WLS model [6]. Although our SG-WLS is a little slower than FGS [20], it can overcome several limitations of FGS [20] and achieves close performance to the WLS model [6].

Storing matrixes $A_{s},P,Q$ in Eq. (9) is the main memory cost of our SG-WLS which is $\mathcal{O}(sr)$. Since the memory cost of Eqs. (12)$\sim$(34) is $\mathcal{O}(s)$, the final memory cost of our SG-WLS is $\mathcal{O}(sr+s)$ which is $\mathcal{O}(Mr^{2})$ or $\mathcal{O}(Nr^{2})$. Note that despite the memory cost of the solver which may be different among different solvers, the memory cost of storing the matrix $A$ in the original WLS model in Eq. (3) is $\mathcal{O}(MNr^{2})$. Thus, the memory cost of our SG-WLS is at most at the magnitude of $max\{\frac{1}{M},\frac{1}{N}\}$ of that of the WLS model.

3 Applications and Experimental Results

We test our SG-WLS on four applications including image detail enhancement, HDR tone mapping, guided depth upsampling and image colorization. For the first two applications which represent applications with dense input data, we adopt $\omega_{i,j}=\omega_{i,j}^{frac}$ in Eq. (2). For the rest two applications, we adopt $\omega_{i,j}=\omega_{i,j}^{exp}$ in Eq. (2). These two applications represent applications with sparse input data which can be regular (guided depth upsampling) or irregular (image colorization). For more experimental results, please refer to our supplementary materials.

Image detail enhancement aims at enhancing the details of an image while avoiding artifacts such as gradient reversals and halos [6, 10]. In our experiments, an image is decomposed into a base layer and a detail layer through edge-preserving smoothing such as FGS [20] and WLS [6]. Parameters of our SG-WLS are set as follows: $r=1,\tau=1,\alpha_{s}=\alpha_{r}=1.2,\lambda=900$. Similar to the result in Fig. 1, we illustrate another example of our detail enhancement results and comparison with results of other methods in Fig. 7. All the results are obtained by enlarging $3$ times of their corresponding detail layer. Visual comparison shows that our SG-WLS can achieve close performance to the WLS [6]. Note that the result of FGS appears pale while results of our SG-WLS and WLS [6] are more colorful. Please zoom in for better visual comparison.

HDR tone mapping is another application that needs edge-preserving smoothing. Based on the multi-scale tone mapping framework proposed by Farbman et al. [6] 111The source code can be downloaded here http://www.cs.huji.ac.il/~danix/epd/, the input image is decomposed into a base layer and three detail layers. The base layer is nonlinearly mapped to a low dynamic range and is re-combined with detail layers. Filters are applied to the logarithmic HDR images. Parameters of our SG-WLS are set as follows: $r=1,\tau=1,\alpha_{s}=\alpha_{r}=1.2$ and $\lambda=5/40/320$ for the first/second/third detail layer. Results are illustrate in Fig. 2 as well as Fig. 8. There are noticeable blocky artifacts in the results of FGS [20] while our SG-WLS can well overcome this limitation and shows close performance to WLS [6].

Guided depth upsampling aims at enlarging the resolution and smoothing the noise of a small noisy input depth map with the guidance of a color image. The input depth map is firstly projected onto the high resolution coordinate. We adopt Eq. (35) for this task where our SG-WLS is applied to both the sparse input and the index map. The guidance weight is based on the guidance image. Results of FGS [20] are obtained in a similar manner. The parameter setting of our SG-WLS is as follows: $r=4,\tau=4,\sigma_{s}=4,\sigma_{r}=3$ and $\lambda=100/200/400$ for $2\times/4\times/8\times$ upsampling. We adopt the parameters used in [16] for FGS [20] to produce the results. Upsampling results of different methods are shown in Fig. 9. The main challenges of guided depth upsampling are texture copy artifacts and blurring edges [18]. As illustrated in highlighted regions, our SG-WLS shows better performance in handing the challenges than compared methods. We further show MAD of different methods in Table 2. Note that the results of our SG-WLS with $r=4,\tau=4$ clearly outperform the results of WLS [6], FGS [20] and our SG-WLS with $r=1,\tau=1$. This also validates the effectiveness of using large neighborhood systems.

Image colorization is colorizing a gray image given user specified scribbles. Two chrominance channels $U$ and $V$ extracted from the input color scribbles are used as sparse input which is propagated with the guidance of the gray image. Similar to guided depth upsampling, we use Eq. (35) for the propagation of the $U$ and $V$ channel. The parameters of our SG-WLS are set as follows: $r=4,\tau=2,\sigma_{s}=4,\sigma_{r}=2,\lambda=900$. We show experimental results of different methods in Fig. 10. Note that the hair in the red circle in the result of FGS [20] is seldom colorized while the ones in our result and the result of WLS [6] are properly colorized.

Conclusion In this paper, we have presented a fast alternative approximation to the Weighted Least Squares (WLS) model, termed Semi-Global Weighted Least Squares (SG-WLS). Both the time cost and the memory cost of our SG-WLS are much more efficient than that of the WLS model while it can achieve close performance to the WLS model in several applications. Our SG-WLS can overcome several limitations of previous related work due to the newly designed neighborhood system construction. Thanks to the proposed fast solution to 1D filters, our SG-WLS is also capable of a more general and larger neighborhood system other than the 4-connected/8-connected neighborhood system adopted by previous work. We show such a generalization can achieve better performance in some applications such as guided depth upsampling. Through experiments of several applications, we show the effectiveness and efficiency of our SG-WLS.

Acknowledgements Correspondence should be addressed to C. Shen or J. Yang. This work was in part supported by NSFC, China ($\#$61572315, 61503250, 61471230) and 973 Plan, China ($\#$2015CB856004).