\section{Introduction}\label{sec:introduction}}

\IEEEPARstart{F}{iltering} techniques, especially these methods with so-called ``feature-preserving'' property, have been used as a basic tool in many signal/image processing fields. This favor may be largely due to the fact that many natural signals, as well as visual clues (color, lighting, depth and surfaces, etc), tend to be locally piece-wise smoothing, among which are discontinuous boundaries that convey an important proportion of visual information, despite they are formed by a minority of sparse features --- singularities and edges.
On the basis of this prior knowledge, these sparse features usually need to be preserved in practical applications. It has been also witnessed the prevailing use of ``feature-preserving'' filtering methods in a wide range of tasks, including denoise~\cite{donoho1995noising, perona1990scale, tschumperle2006fast}, image tone mapping and high dynamic range (HDR) image compression~\cite{farbman2008edge, he2012guided, tomasi1998bilateral}, details enhancement~\cite{farbman2008edge, he2012guided, paris2011local, xu2011image},  stylization~\cite{kyprianidis2011image, xu2011image, xu2012structure}, dehazing~\cite{he2012guided}, depth filtering~\cite{barron2016fast, li2016deep, shen2015mutual, Young_2019_ICCV, zhang2014rolling} and so on, as well as the extensions in 3D cases: mesh and point clouds refinements~\cite{fleishman2003bilateral,he2013mesh,meyer2003discrete,oztireli2009feature,wang2016mesh,zhang2015guided}.

In the literature, a large number of filtering algorithms have been proposed under different purposes or research backgrounds. Despite the varying forms and properties, they always share a very similar goal --- exploring, implicitly or explicitly, a special signal/image decomposition to easily extract or analyze the useful information. In general, these filtering methods can be divided into two broad categories: local and global filters. In local filters, each output can be simply computed as a weighted sum of local neighbors, in which the weights can be pre-defined or computed using local information. This strategy can be found in many typical instances including median filter~\cite{weiss2006fast}, bilateral filter~\cite{tomasi1998bilateral}, guided filter~\cite{he2012guided} and their extensions~\cite{barron2016fast,chen2007real,kou2015gradient,li2014weighted,weiss2006fast}. These local methods have been well studied for their simplification and ease of implementation. One bottleneck of these local filtering methods is that they may cause some unexpected effects in some applications, for example, halo artifacts in image detail enhancement~\cite{tomasi1998bilateral,chen2007real} and HDR image compression~\cite{fattal2002gradient,tomasi1998bilateral} because of the ignorance of the global attributes. In contrast, the global methods are usually based on an optimization framework to find a globally optimal solution to an objective function involving a data constraint term and a prior smoothness term. Due to the global property, the optimization-based approaches, including total variation (TV)~\cite{rudin1992nonlinear}, weighted least squares (WLS) filter~\cite{farbman2008edge,min2014fast} and $L_0$-norm gradient minimization~\cite{ono2017l,xu2011image}, to some extent, overcome the artifacts with better performance in a range of image processing and computer graphics tasks. Despite the success of these optimization-based methods, the high performance is generally achieved at the much-increased cost of solving a large-scale linear or no-linear system~\cite{ono2017l, xu2011image}, which is always not negligible even with the recent endeavor in accelerations.

In this paper, we propose a novel semi-sparse minimization framework for smoothing filters. This new model is inspired by the sparsity-inducing models proposed in the cutting-edge filtering methods~\cite{ono2017l,xu2011image,ye2013sparse}. As illustrated therein, sparsity prior knowledge is usually formulated into an optimization-based problem associated with some sparsity-induced regularizers. It has shown that these methods can produce high-quality filtering results and are particularly suitable for preserving sparse features (singularities and edges). However, as shown in Fig. 1, they may fail to produce high-quality fitting/filtering results in the regions where the sparse features and polynomial smoothing surfaces coexist. It fails because the sparse gradient prior is not fully admitted in polynomial smoothing surfaces and imposing sparse regularization in these regions will produce severe (step-like) artifacts. In contrast, we here show that a so-called semi-sparsity has more favorable properties for fitting the sharpening edges and polynomial-smoothing surfaces, which, as a result, can be utilized to alleviate the drawbacks existing in the sparsity-induced models. We interpret how to reformulate this kind of prior knowledge and show its benefits to smoothing filters under a semi-sparse optimization model. The major contributions of this paper are summarized as follows:

\begin{itemize}
	\setlength{\itemsep}{0pt}
	\setlength{\parsep}{0pt}
	\setlength{\parskip}{0pt}

	\item A semi-sparsity prior knowledge is introduced with the verification on natural images, and its benefits to smoothing filters are discussed under an optimization-based framework.

	\item A new minimization is designed for the semi-sparsity prior knowledge with $L_0$-norm regularization in the higher-order gradient domains. This gives rise to a new filtering method with the powerful simultaneous-fitting ability in both sparse features  and non-sparse regions.

	\item A half-quadratic splitting algorithm is introduced to solve the $L_0$-norm minimization problem, in which the fast Fourier transform (FFT) is also employed for acceleration.

	\item Experimental results on 2D images and 3D meshes are provided to show the favorable properties of the semi-sparsity model. A series of visual results are presented in comparison of the sparse-induced $L_0$ gradient smoothing filters, as well as the cutting-edge ``feature-aware'' filtering methods.

\end{itemize}

We remark here that the simultaneous-fitting ability makes it possible to remove the low-amplitude structures in the signals and preserve the minority but significant edges and singularities without introducing notorious stair-case artifacts in polynomial-smoothing surfaces. That is to say,  our new model can be used as a ``feature-aware'' filter in a wide range of fields as many existing filtering methods. The demonstration of its versatility is discussed in more detail substantially in cooperation with its benefits to many image processing and computer graphical applications.

The paper is organized as follows. In Section 2, we discuss the related approaches under a generalized optimization framework. In Section 3, we derive our approach in detail.  The $L_0$ minimization and its limitations are briefly interpreted in Sections 3.1. The advantages of semi-sparsity prior knowledge to $L_0$ gradient minimization and its universal applicability in natural images are respectively presented in Section 3.2 and 3.3. In Section 4, we formulate the semi-sparsity prior as an optimization-based problem with $L_0$-norm regularization in the higher-order gradient domains. The solver and its acceleration are discussed as well. In Section 5, a more detailed analysis is discussed. In Section 6, we provide the experimental results for smoothing filters. In Section 7, we first show the results on 2D images on several applications, and then extend our semi-sparse model to 3D cases and show the advantages to mesh and point cloud denoising. In Section 8, we draw our concluding remarks.

\section{Related Work}
Within the signal/image processing fields, a number of filtering algorithms have been proposed to reduce noise, or to perform a special signal decomposition. Mathematically, many of them can be, explicitly or implicitly, formulated into an optimization-based framework with the following form:
\begin{equation}
	\begin{aligned}
		\mathop{\min}_{u} {{\left\Vert {u}-Af\right\Vert }_2^2}+\lambda\mathcal{R}(u),
	\end{aligned}
	\label{f2}
\end{equation}
where, ${f}$ and ${u}$ are the observed and output signals\footnote{We, throughout, use the symbols ${f}$ and $u$ as functions or operators and the bold ones ${\boldsymbol{f}}$ and $\boldsymbol{u}$ as their discrete counterparts.}, the weight $A$ specifies a spatially-varying confidence map of smoothness of $f$, $\mathcal{R}(u)$ is known as the regularization term; and the global parameter ${\lambda}>0$ specifies a trade-off between fidelity term and regularization term. We review the internal relationships of the existing works based on Eq. (1) despite their varying forms and generalizations.

Firstly, we interpret a type of filtering methods in the context of Eq. (1), where only weight ${A}$ is considered and the regularization term $\mathcal{R}(u)$ is ignored, or set ${\lambda}=0$ equivalently. As interpreted in~\cite{milanfar2012tour}, the general construction of $A$ begins by specifying a kernel function $K$ in each sample. The kernel can in fact depend on the given data samples, the locations of these samples and even pre-computed features.
For example, it is easy to define a general Gaussian filter kernel $K$ with the elements,
\begin{equation}
	\begin{aligned}
		{\boldsymbol{K}_{i,j}={\text{exp}( - {(\boldsymbol{x}_i-\boldsymbol{x}_j)}^T \boldsymbol{Q}_{i,j} {(\boldsymbol{x}_i-\boldsymbol{x}_j)}) }},
	\end{aligned}
	\label{f2}
\end{equation}
where, $\boldsymbol{x}$ is a feature vector and ${\boldsymbol{Q}}$ is symmetric positive definite (SPD) matrix. Let ${\boldsymbol{x=u}}$, ${\boldsymbol{Q}_{i,j}} = \frac{1}{h^2}\boldsymbol{I}$ (${h}$ is a control parameter related to the variance of noise). In matrix language, one has $\boldsymbol{A}=\boldsymbol{\Lambda}^{-1}\boldsymbol{K}$, where $\boldsymbol{\Lambda}$ is nontrivial diagonal matrix with diagonal elements $\boldsymbol{\Lambda}_{j,j} = \sum_{i}{\boldsymbol{K}_{i,j}}$\footnote{The reader is referred to the surveys~\cite{milanfar2013symmetrizing} for more detailed discussion of their relationship.}. As we can see that $\boldsymbol{A}$ is a (row-) stochastic matrix as its rows sum to one. This is always necessary because it keeps the mean value of the signal unchanged during the filtering process. The above configuration can be found in many representative instances. For example, the median filters~\cite{weiss2006fast,  zhang2014100} use the (weighted) median value of samples in a local area; the famous bilateral filter~\cite{barron2016fast,  elad2005retinex, li2016fast,   porikli2008constant, tomasi1998bilateral} takes both spatial and data-wise distances into account. By analogy, non-local means (NLM) filter~\cite{buades2011non} is obtained when extending the model into patch-wise data. It is also possible to compute the corresponding filter kernel for locally adaptive regression kernel (LARK)~\cite{takeda2007kernel}. Moreover, we claim that domain transform filter~\cite{gastal2011domain} and guided image filter~\cite{he2012guided,  kou2015gradient, li2014weighted, ochotorena2019anisotropic} belong to this case as well, as they explicitly define a Laplacian kernel matrix depending on the local information.

It is also easy to see from Eq. (2) that the kernel ${K}$ decays exponentially, therefore the contribution of the samples far from the center can be neglected in practice. This
helps to reduce the computational cost, because only a few neighbors are sufficient to provide the high-quality approximate results. It is worthy to notice that many accelerations and extensions of the above local filters also benefit from this local property. The low computational cost is one of main reasons that the local filters always attract widespread attention in practical applications.

Simultaneously, there is another way for filtering methods based on the optimization-based framework of Eq. (1). These filtering methods have also emerged as powerful tools for their remarkable feature-aware filtering results. In the context of smoothing filters, the weight (matrix) function $\boldsymbol {A}$ is always assumed to be identity matrix ($\boldsymbol{A = I}$). The filter result is highly dependent on the regularization term $\mathcal{R}(u)$. The regularization technique can be traced back to image restoration, where a quadratic data-fidelity term is used to measure the distance between input and output, and an extra regularization~\cite{rudin1992nonlinear} is added to penalize the strength of smoothness of output.  Differing from the above local filters, the global solution of Eq. (1) is emphasized. The methods using this strategy are usually classified as global filters. The related works can be found from the range of total variation (TV)~\cite{rudin1992nonlinear} to weighted least squares (WLS) frameworks~\cite{farbman2008edge, kim2017fast, min2014fast},  $L_0$ gradient minimization ~\cite{nguyen2015fast, ono2017l, xu2011image} and so on. The difference mainly lies in the usage of prior knowledge for regularization. For example, the TV-based regularization~\cite{rudin1992nonlinear} has a tendency to preserve the singularities and strong edges and produce piece-wise smoothing results. In the WLS filter,  a novel weighted $L_{2}$ gradient norm is exploited to provide visual-friendly smoothing outputs, which can be used in image enhancement with halo-free results. When the penalty is changed into a sparse gradient regularization, it gives $L_0$ gradient minimization smoothing methods; and it has turned out that such a sparse regularization favors piece-wise constant solutions with excellent approximation ability in the flattening regions.

Despite the different forms and generalizations, these methods all strive to explore an adaptive filtering strategy --- smoothing data but allowing to preserve singularities and discontinuities. They have been used as powerful tools for the remarkable feature-aware filtering results. However, these optimization-based methods need to solve a linear or non-linear system, which is time-consuming due to the large-scale data samples in practice. Despite the recent endeavor in developing acceleration techniques, they are typically still an order of magnitude slower than the local ones, limiting the use in many real-time applications.

\section{Preliminaries}

We start our work with a brief review of $L_0$ gradient minimization and then extend it into higher-order gradient domains to give our semi-sparse smoothing model. Moreover, we also verify its universally applicable in natural images.

\subsection{$L_0$ Gradient Minimization}
The $L_{0}$ norm\footnote{It is not a true norm in the rigorous sense, but intuitively gives a useful tool for comparing two vectors. } counts the number of its non-zero entries, which directly measures the degree of sparsity of the vector. As a metric, it has been proved to be capable of capturing the singularities and discontinues of signals. In signal and image processing, a corresponding $L_0$ regularization is always a favorable and powerful mathematical tool to formulate the sparsity-induced priors. It has been used in a huge amount of existing works, including sparse signal/image recovering~\cite{donoho2003optimally, elad2010sparse} sparse representation and compressed sensing~\cite{donoho2006compressed} and so on.

In the context of filtering techniques, it has proved that $L_0$-norm minimization can be directly used for optimizing the gradients of an image ~\cite{nguyen2015fast, ono2017l, xu2011image}. More specifically, let $\nabla u$ be the gradient field of image $u$, it is suggested in ~\cite{xu2011image} to minimize the objective function\footnote{In regard to Eq. (1), ${W}$ is reduced into Kronecker delta function, and the smoothness of filter is determined by the regularization ${\left\Vert \nabla u \right\Vert }_0$.},
\begin{equation}
	\begin{aligned}
		\mathop{\min}_{u} {{\left\Vert {u}-{f}\right\Vert }_2^2}+\lambda{\left\Vert \nabla u \right\Vert }_0,
	\end{aligned}
	\label{f2}
\end{equation}
where, ${{\left\Vert \cdot \right\Vert}_0}$ is the so-called $L_0$-norm.
It turns out that the $L_0$ norm regularization provides a powerful tool for creating piece-wise constant results. This idea underpins several sparse gradient models~\cite{nguyen2015fast, ono2017l, xu2011image}, in which the gradients of signals are usually assumed to be sparse---only a tiny part of the gradient magnitudes is relatively large and the remainder is negligible for approximating to zero.

\subsection{Semi-sparsity Priors}

Despite of the advantages of $L_0$ gradient minimization, it may fail in some cases when sparse gradient prior is not fully admitted in such as polynomial smoothing surfaces. As shown in Fig. 1, a 1D noisy signal $u$ is composed of singularities, sharp edges, and slope line (or, polynomial surfaces). Clearly, ${{\nabla} u}$ is not fully sparse within the area of slope line, because the gradient values are relatively larger than zero. Thus, the sparse gradient prior fails here and produces stair-case results (blue). This failure is largely due to the fact that sparsity prior in the gradient domain is no longer fully satisfied within the dense counterparts. This motivates us to resort to a new regularization for better fitting results.
\begin{figure}[t]
	\begin{center}
		\begin{minipage}{0.5\textwidth}
			{\includegraphics[width=\textwidth]{./imgs/smooth_1d.jpg}}
		\end{minipage}
	\end{center}
\caption{An example of 1D signal smoothing. This signal is composed of the spikes (singularities), sharp edges, and slope line (or, polynomial surface). It is clear that $L_0$ gradient minimization~\cite{xu2011image} preserves spikes and sharp edges successfully but produces stair-case artifacts within the slope region (\textcolor{blue}{blue line}). In contrast, our model attains much better fitting (filtering) results in polynomial surfaces with comparable results in spikes and sharp edges (\textcolor{red}{red line}).}
\end{figure} %

As interpreted in Fig. 1, the non-zero elements of gradient ${{\nabla} u}$ is not fully sparse but are densely distributed in some region. We say $u$ is a so-called semi-sparse signal in gradient domains in this case. On the other hand, it is easy to verify that the second order gradient  ${{\Delta} u}$ is sparse. Without loss of generality, we define that $u$ is a semi-sparse signal, if the higher-order gradients ${\nabla}^{n-1} u$  and ${\nabla}^{n} u$ satisfy,
\begin{equation}
	\begin{aligned}
		\begin{cases}
			& {\left\Vert {\nabla}^{n} u\right\Vert }_0 < M,\\
			&{\left\Vert {\nabla}^{n-1} u\right\Vert }_0 > N,
		\end{cases}
	\end{aligned}
	\label{f2}
\end{equation}
where, ${\nabla}^{n}$ is a $n$-th (partial) differential operator, and $M, N$ are some appropriate natural number with ${M\ll N}$. The idea of Eq. (4) is easy to understand, that is, the energy of ${\nabla}^{n} u$ is much smaller than that of  ${\nabla}^{n-1} u$. Intuitively, this happens if and only if the support\footnote{The support of vector $u$ is its number non-zeros elements.} of ${\nabla}^{n} u$ is much less than that ${\nabla}^{n-1} u$; in other words, ${\nabla}^{n} u$ is much more sparse that ${\nabla}^{n-1} u$. Taking a piece-wise polynomial function with $n$ degree into account, it is clear that the $n$-th order gradient ${\nabla}^{n} u$ is sparse, but that is not the case for $k$-th $(k<n)$ order gradient ${\nabla}^{k} u$. This observation thus motivates us to extend Eq. (3) into higher-order gradient domains to regularize the semi-sparse signals with singularities and polynomial surfaces. In practice, the choice of $ n $ is usually determined by the property of signals.

\subsection{Verification on Natural Images}
In this section, we verify the semi-sparsity prior knowledge and show its universal applicability in natural images. It is well-known that the gradient of natural images has a heavy-tailed distribution is usually characterized as mixture of Gaussians, Laplacian or hyper Laplacian distributions in image denoising~\cite{rudin1992nonlinear, zuo2014gradient}, debluring~\cite{xu2013unnatural}, or enhancement~\cite{gong2014image}. We here interpret that our semi-sparsity prior knowledge can be derived from the distributions of higher-order gradients of images.

For the sake of this purpose, we use Kodak image dataset~\cite{kodak} for interpretation and show the normalized distributions of higher-order gradients ${\nabla}^{n} u$ in Fig. 2.  Clearly, the distributions reveal heavy tails in the log-scale domain, which is consistent with the conclusions in many sparsity-induced models. Moreover, it is clear that more and more values of high-order gradients approach to  $0$, as the order $n$ increases. It means that the sparsity is enhanced when the order $n$ goes higher. On the other hand, the gap of sparsity between the $1^{st}$ and $2^{nd}$ gradients is much bigger than that between the $2^{nd}$ and $3^{rd}$ gradient. This trend also holds in the higher-order cases and implicit that the sparsity tends to be bounded when increases the order $n$. This trend also implies that the areas belonging to polynomial surfaces of degree $ n> 2 $ are relatively small and can be ignored in natural images.

\begin{figure}[h]
	\begin{center}
		{\includegraphics[width=0.4\textwidth]{./imgs/gradient_distribution.jpg}}
	\end{center}
	\caption{The statistic distributions of the higher-order ($n=1\scriptsize{\sim}5$) gradients learned from the Kodak dataset~\cite{kodak}.}
\end{figure} %

 The above claims intuitively motivate us to chose $n=2$ as the highest order in most of our subsequent experiments. This choice has extra two-fold benefits: (1) $n=2$ is the simplest semi-sparsity case, and (2) the computational cost of the proposed model is also moderate when $n$ is small. The more profound theoretical justification of the choice of the order $n$ is out of the scope here.

\section{Formulation}
According to Eq. (4), we here propose an optimization-based framework for semi-sparse filtering as follows:
\begin{equation}
	\begin{aligned}
		\mathop{\min}_{u} {{\left\Vert {u}-{f}\right\Vert }_2^2}+\alpha \sum_{k=1}^{n-1}  {\left\Vert {\nabla}^{k}u-{\nabla}^{k} f \right\Vert }_2^2+\beta{\left\Vert {\nabla}^{n} u \right\Vert }_0
	\end{aligned}
	\label{f2}
\end{equation}
where $\alpha$ and $\beta$ are the weights. The first term in Eq. (5) is data fidelity. The second term measures the similarity of higher-order gradients ${\nabla}^{k}u$ and ${\nabla}^{k}f$ in consideration of the piece-wise polynomial surfaces. The third term  ${\left\Vert {\nabla}^{n} u \right\Vert }_0$ favors the highest-order gradient ${ {\nabla}^{n} u}$ to be fully sparse.  The idea is straightforward, that is, a sparse-induced $L_0$-norm constraint is only imposed on the highest order $n-th$ gradient domain, as the ones less than $n$ is not fully sparse. Mathematically, Eq. (5) can be viewed as a generalized extension of Eq. (3) in higher-order gradient cases.

Due to the non-convexity of $L_0$-norm, a direct solver to Eq. (5) is not available. In general, it can be solved with greedy algorithms~\cite{tropp2004greed} to iteratively select the sub-optimal solutions, or be reduced with some appropriate relaxations, for instance, replacing $L_0$-norm by its convex relaxation $L_1$-norm~\cite{chen2001atomic, donoho2003optimally, wang2008new}. We here use an alternating minimization algorithm to split it into two sub-problems based on variable-splitting and penalty techniques~\cite{wang2008new, xu2011image}.

For the sake of simplicity, we take the second-order gradient of Eq. (5) into account. By introducing an auxiliary variable $w$, we rewrite Eq. (5) as,
\begin{equation}
	\begin{aligned}
		\mathop{\min}_{u,w} {{\left\Vert {u}\!-\!{f}\right\Vert }_2^2}+\alpha {\left\Vert {\nabla}u \!- \!{\nabla}f\right\Vert }_2^2 +\beta {c(w)}  + \lambda {\left\Vert \Delta u \!-\! w \right\Vert }_2^2
	\end{aligned}
	\label{f2}
\end{equation}
where $\Delta$ is the Laplace operator and $c(w) = {\left\Vert \Delta w \right\Vert}_0$, and $\lambda$ controls the similarity between the variables $w$ and ${\Delta}u$. It is also well-known that the solution of Eq. (6) globally converges to that of Eq. (5) when $\lambda \to \infty$. As explained in~\cite{wang2008new, xu2011image},  the above relaxation has two-side benefits: (1) by introducing the variable $w$, the objective function can be iteratively minimized with respect to one and the other, and (2) each problem has a closed-form formula with low computational complexity by fixing either of two variables $u$ and $w$ during the minimization scheme.

\textbf{Sub-problem 1:}  Clearly, the objective function of Eq. (6) reduces to a quadratic function with respect to $u$ by fixing $w$, giving an equivalent minimization problem,
\begin{equation}
	\begin{aligned}
		\mathop{\min}_{u} {{\left\Vert {u}-{f}\right\Vert }_2^2}+\alpha {\left\Vert {\nabla}u -{\nabla} f \right\Vert }_2^2 +\lambda{\left\Vert \Delta u -w \right\Vert }_2^2.
	\end{aligned}
	\label{f2}
\end{equation}

It is easy to solve $u$ for the quadratic property of Eq. (7). For simplicity, we consider the discrete case and denote  $\boldsymbol{D}$ and $\boldsymbol{L}$ as discrete gradient operator and Laplace operator in matrix-form, respectively. The optimal solution $u$ is then given by the following linear system:
\begin{equation}
	\small
	\begin{aligned}
		(\boldsymbol{I} + \alpha \boldsymbol{D}^T \boldsymbol{D} + \lambda  \boldsymbol{L}^T \boldsymbol{L}) u =  	(\boldsymbol{I} + \alpha \boldsymbol{D}^T \boldsymbol{D}) f + \lambda  \boldsymbol{L}^T w,
	\end{aligned}
	\label{f2}
\end{equation}
where $\boldsymbol{I}$ is identity matrix. Clearly, the left-hand
side of Eq. (8) is symmetric and semi-positive, several solvers such as Gauss-Seidel method and preconditioned conjugate gradients (PCG) method~\cite{saad2003iterative} can be used directly. However, these solvers are still time-consuming when the size of problem is huge, millions variables for example, in practical tasks.

Alternatively, the fast Fourier transforms (FFTs) can be used for acceleration. As suggested in~\cite{wang2008new, xu2011image},  $\boldsymbol{D}^T \boldsymbol{D}$ and $\boldsymbol{L}^T \boldsymbol{L}$ are all block circulant under the periodic boundary condition\footnote{In our experiments, we use periodic boundary conditions and compute the FFTs of constants and higher-order operators in the initialization step.}; therefore, we solve the corresponding linear equation with Fast Fourier Transforms (FFTs), yielding the solution
\begin{equation}
	\small
	\begin{aligned}
		u =
		\mathcal{F}^{-1}\!\Big(\frac{\mathcal{F}(f)\!+\!\alpha{\mathcal{F}(\nabla)}^\ast\!\circ\!{\mathcal{F}(f)}\!+\!\lambda{\mathcal{F}(\Delta)}^\ast\!\circ{\mathcal{F}(w)}}
		{\mathcal{F}(1)\!+\!\alpha {\mathcal{F}(\nabla)}^\ast\!\circ\!\mathcal{F}(\nabla)\!+\!\lambda{\mathcal{F}(\Delta)}^\ast\!\circ\!\mathcal{F}(\Delta)}\Big),
	\end{aligned}
	\label{f2}
\end{equation}
where  $\mathcal{F}$ is fast Fourier transform(FFT) operator and $\mathcal{F}(1)$ is the FFT of the delta function; while ``$\ast$" denotes complex conjugacy and ``$\circ$" is component-wise multiplication. The division is a component-wise operator as well. The benefit of Eq. (9) is that the formulation only involves one FFT and inverse FFT and component-wise operators. It is much faster than solving the huge matrix equations.
\begin{figure*}[t]
	\begin{center}
		\subfigure[Noisy Input]
		{\includegraphics[width=0.19\textwidth]{./imgs/Stripe/Noisy.png}}

		\subfigure[Guided Filter~\cite{he2012guided}]
		{\includegraphics[width=0.19\textwidth]{./imgs/Stripe/gf.png}}

		\subfigure[WLS Filter~\cite{farbman2008edge}]
		{\includegraphics[width=0.19\textwidth]{./imgs/Stripe/WLS.png}}
		\subfigure[$L_0$ Smoothing~\cite{xu2011image}]
		{\includegraphics[width=0.19\textwidth]{./imgs/Stripe/L0.png}}
		\subfigure[Our Result ($n=2$)]
		{\includegraphics[width=0.19\textwidth]{./imgs/Stripe/SAP.png}}

	\end{center}
	\caption{Noisy image smoothing. (a) Noisy input, (b) Bilateral filter ($\delta_s = 30, \delta_r = 0.2$)~\cite{tomasi1998bilateral}, (c) Guided filter ($r=15, \epsilon =0.1^2$)~\cite{he2012guided}, (d) Total variation (TV) ($\lambda = 1.35$)~\cite{rudin1992nonlinear},  (e) WLS filtering ($\alpha = 2.0, \lambda = 1.2$)~\cite{farbman2008edge}, (f) $L_0$ smoothing ($\lambda = 0.01, \kappa =1.05$)~\cite{xu2011image}, (g) Our result ($\alpha = 0.5, \beta =0.01, \kappa =1.05, \tau=0.8$), (h) Ground truth.  PSNR, (a)$\scriptsize{\sim}$(g): 30.12, 32.54, 37.24, 37.80, 42.64, 40.61, 41.64, 44.65.}
\end{figure*} %

\textbf{Sub-problem 2:}  By analogy, $w$ can be solved with a fixed $u$. In this case, the first two terms in Eq. (6) are constant with respect to $w$, and it is equivalent to solving the following problem,
\begin{equation}
	\begin{aligned}
		\mathop{\min}_{w} {\left\Vert \Delta u -w \right\Vert }_2^2 + \frac{\beta}{\lambda} {c(w)}.
	\end{aligned}
	\label{f2}
\end{equation}
where $c(w)$ counts the number of non-zero elements in $w$. As demonstrated in~\cite{xu2011image}, Eq. (10) can be spatially decomposed into elements, in which each variable $w_i$ is estimated individually. The optimal solution of Eq. (11) is attained when the $i$-th variable $w_i$ satisfies,
\begin{equation}
	\begin{aligned}
		{w_i} =
		\begin{cases}
			0, &\quad {\left\Vert \Delta u_i \right\Vert }_2^2<\frac{\beta}{\lambda},\\
			\Delta u_i, &\quad {\left\Vert \Delta u_i \right\Vert }_2^2 \geq \frac{\beta}{\lambda}.
		\end{cases}
	\end{aligned}
	\label{f2}
\end{equation}
It is clear that Eq. (11) is a hard-threshold operator and can be solved efficiently. This simplification is the main benefit of the splitting scheme, as it makes the altered problem empirically solvable even with a huge scale of variables.

\begin{algorithm}[!t]
	\caption{Semi-Sparse Smoothing Filtering}
	\label{euclid}
	\begin{algorithmic}[1]
		\REQUIRE Signal ${f}$, weights ${\alpha, \beta, \lambda}$, parameters $\lambda_0$, $\lambda_{max}$ \\and rates $\kappa$ and $\tau$.
		\STATE $u \leftarrow f$, $\lambda \leftarrow \lambda_{0}$, ${i \leftarrow 0}$.
		\WHILE{$\lambda<\lambda_{max}$}
		\STATE With $u^i$, solve $w_i$  for Eq. (6).
		\STATE With $w^i$, solve $u_{i+1}$ for  Eq. (8).
		\STATE $\lambda \leftarrow \kappa \lambda$, $\alpha  \leftarrow \tau \alpha$, $i++$.
		\ENDWHILE
		\ENSURE Smoothing result ${u}$.
	\end{algorithmic}
\end{algorithm}

\section{More Analysis}
\textbf{Iterations$\And$Convergence:} Notice that the solution of Eq. (6) converges to that of Eq. (5) on the premise that $\lambda$ is large enough. However, it is generally difficult to determine the best value under a given accuracy. We take an adaptive strategy as explained in~\cite{wang2008new,xu2011image}, where $\lambda$ is increased with $\kappa$ $(\kappa>1)$ times in each iteration. A similar configuration is also introduced for $\alpha$, yet, decreasing with $\tau$ $(0 \le \tau \le 1)$ times in each iteration, as the second term in Eq. (5) plays the same data-fidelity role and has little influence on smoothing results.  The whole alternative minimization scheme is summarized in Algorithm 1.
\begin{figure*}[!t]
	\begin{center}
		\subfigure[Input Image]
		{\includegraphics[width=0.245\textwidth]{./imgs/woman/Noisy_local.png}}
		\subfigure[BF ($\sigma_s = 5, \sigma_r = 0.2$)~\cite{weiss2006fast}]
		{\includegraphics[width=0.245\textwidth]{./imgs/woman/BF_local.png}}
		\subfigure[$L_0$ Smoothing($\lambda = 0.02$)~\cite{xu2011image}]
		{\includegraphics[width=0.245\textwidth]{./imgs/woman/L0_local.png}}
		\subfigure[Ours ($\alpha = 0.1, \beta = 0.01$)]
		{\includegraphics[width=0.245\textwidth]{./imgs/woman/SP_local.png}}
	\end{center}
	\caption{Visual comparison on a natural face image. As suggested in~\cite{weiss2006fast}, the filters are employed in the logarithmic intensity for better visual results.}
\end{figure*}
\begin{figure*}[!b
	\begin{center}
		\subfigure[Input]
		{\includegraphics[width = 0.245\textwidth]{./imgs/girl/SAP_src.png}}
		\subfigure[$\alpha = 10, \beta =0.002$]
		{\includegraphics[width = 0.245\textwidth]{./imgs/girl/SAP_alpha_10_lambda_005_res_local.png}}
		\subfigure[$\alpha = 10, \beta =0.01$]
		{\includegraphics[width = 0.245\textwidth]{./imgs/girl/SAP_alpha_10_lambda_01_res_local.png}}
		\subfigure[$\alpha = 10, \beta =0.1$]
		{\includegraphics[width = 0.245\textwidth]{./imgs/girl/SAP_alpha_10_lambda_05_res_local.png}}
		\subfigure[Close-ups]
		{\includegraphics[width = 0.245\textwidth]{./imgs/girl/SAP_contrast_local.png}}
		\subfigure[$\alpha = 0.1, \beta =0.002$]
		{\includegraphics[width = 0.245\textwidth]{./imgs/girl/SAP_alpha_01_lambda_005_res_local.png}}
		\subfigure[$\alpha = 0.1, \beta =0.01$]
		{\includegraphics[width = 0.245\textwidth]{./imgs/girl/SAP_alpha_01_lambda_01_res_local.png}}
		\subfigure[$\alpha = 0.1, \beta =0.1$]
		{\includegraphics[width = 0.245\textwidth]{./imgs/girl/SAP_alpha_01_lambda_05_res_local.png}}
	\end{center}
	\caption{Our semi-sparsity smoothing results of an 2D image with different  parameters $\alpha$ and $\beta$. In each row, $\alpha$ is fixed and $\beta$ is gradually increased to give more smoothing results; while, in each column, $\beta$ is fixed and $\alpha$ is increased to remove more strong local details. Note that all filtering results are generated in the logarithmic domain for better visual effects.}
\end{figure*}

The numerical results demonstrate the efficiency of iterative hard-thresholding algorithm, however the theoretical convergence guarantees are quite high as pointed in~\cite{sun2017convergence}.
The half-quadratic splitting algorithm is analogical to that of $L_1$-norm regularization~\cite{wang2008new}, while the convergence result there can not be extended to $L_0$-norm case for its non-convexity. Convergence analysis of iterative hard-thresholding for compressed sensing is also proposed under the restricted isometry property (RIP)~\cite{blumensath2009iterative}, in which the $L_0$-norm regularization is for signal, not its higher-order gradients. The convergence analysis is possible to construct based on iterative thresholding algorithms~\cite{bredies2015minimization,sun2017convergence,zeng2016sparse}. The reader familiar with them or more interested in the convergence result may refer therein to the details.

\textbf{Parameter Settings:}  The global parameters $\alpha, \beta$  and $\lambda$ are important to control the smoothness of results. We, in general, have $\alpha \in [0.01, 10]$, $\beta \in [0.001, 0.1]$ and $\lambda_0 =2\beta$. Notice that increasing $\lambda$ gradually removes the local tiny details until only large-scale or dense features are reserved. We here empirically set $\lambda_{max} = 1.0e5$ in most of our cases. Additionally, the scale factors $\kappa$ and $\tau$ are introduced to speed up the convergence of our alternating minimization framework~\cite{he2013mesh, wang2008new, xu2011image}. We set $\kappa = 2.0$ and a smaller $\kappa$ produces more near-ideal results.
In contrast, we also set $\tau = 0.5$ to ensures the smoothing results to be more and more sparse in the highest-order gradient domain. Such configurations provide a good choice to keep balance between efficiency and performance.

\textbf{Computational Complexity:} Benefiting from the alternating minimization scheme, the semi-sparse model, as a result, can be solved efficiently. In Algorithm 1, it has two main phases: linear system solver and hard-threshold shrinking in each iteration. The computation is dominated by the FFTs and inverse transforms. Each is at the cost of $O(N log(N))$, where N is the total number of pixels in an image. The hard-threshold shrinking can be computed in-place for its separability. The number of iterations controls the total time of the algorithm and it roughly needs 15 iterations to produce high level visual results  in most cases. Our Matlab implementation runs on the PC with Intel Core2 Duo CPU 2.13G. It takes roughly 3 seconds to process a $600 \times 400$ resolution color image.  The code is available online with the link \url{https://github.com/QingXin96/semi-sparsity-smoothing-filters}.
\begin{figure*}[!t
	\begin{center}
				\subfigure[Input (PSNR = 34.19)]
				{\includegraphics[width=0.3\textwidth]{./imgs/order/input_mesh.jpg}}
				\subfigure[Second-order (PSNR = 43.27)]
				{\includegraphics[width=0.3\textwidth]{./imgs/order/smoothing_order2_mesh.png}}
				\subfigure[Third-order (PSNR = 44.05)]
				{\includegraphics[width=0.3\textwidth]{./imgs/order/smoothing_order3_mesh.png}}
	\end{center}
		\caption{Comparison of the smoothing results on a synthesized 2D surface with the slope and $3^{rd}$ degree polynomial regions. The sparse regularization is imposed on the second-order ($n=2$) and higher-order ($n=3$) gradient domains, respectively.}
\end{figure*}

\section{Experimental Results}
In this section, we show the property of simultaneous-fitting ability of the semi-sparse minimization and the benefits to produce the ``edge-preserving'' smoothing results. Notice that  ``edge-preserving''  is very important property for image smoothing tools to remove the small scale details but preserve the sharp variations in image intensities. Specifically, we show the smoothing results on both synthesis and natural images. We also compare the proposed approach with several well-known methods, including bilateral filter~\cite{tomasi1998bilateral}, total variation method~\cite{rudin1992nonlinear}, guided image filter~\cite{he2012guided}, weighted least square (WLS) filter~\cite{farbman2008edge} and $L_0$ gradient minimization~\cite{xu2011image}. For the sake of comparison, the parameters of these filtering methods are either carefully hand-tuned to reach a similar level of smoothness or are selected based on a greedy search strategy for better visual results.

\textbf{Smoothing Filters:}  We first simulate an 1D noisy signal as shown in Fig. 1 and compare the performance of preserving the sharpening edges. As we can see, the cutting-edge $L_0$ gradient minimization~\cite{xu2011image} attains high-quality results in spikes and step-wise edges, but also produces stair-case results in the interval of slope line. As explained earlier, this failure is caused by the invalidity of sparse gradient regularization in the corresponding polynominal regions. In contrast, our semi-sparse model gives more accurate fitting results in the slope and also attains comparable performance in singularities and step-wise edges. In Fig. 3, we also show the results on a synthesized color image consisting of a polynomial surface ({\color{green} green rectangle}) and constant regions. As we can see, our method attains the same performance as that in 1D case, that is, a high-quality fitting result in the polynomial surface; while $L_0$ minimization~\cite{xu2011image} preserves the edges but produce stair-case artifacts in green region and (WLS) filter~\cite{farbman2008edge} also gives a high-quality result, yet, with slight deviations in green region. Other methods, to some extend, blur the strong edges or produce  piecewise constant results in green region.

To demonstrate the virtues of our semi-sparse smoothing filter, we further show the results on natural images. As shown in Fig. 4, the given image contains multiple polynomial surfaces (face region) and strong edges (hair and face contour). As suggested in~\cite{weiss2006fast}, we apply the filters in the logarithmic intensity for more visual-pleasant results. It is clear that our method gives rise to a more pleasurable smoothing result, as it removes the details in face region but also retains high-quality edges in hair region. While the L0 gradient minimization tends to produce the stripe-like results in the polynomial (face) surfaces.

For complementary, we briefly explain the roles of parameters ${\alpha, \beta}$ in our semi-sparse filtering algorithm. A visual comparison with varying parameters $\alpha$ and $\beta$ is presented in Fig. 5. In theory, these parameters need to be carefully configured to weigh the balance of three terms in Eq. (6). We conclude that $\beta$ plays a more important role in controlling the smoothness of results and the outputs are more vulnerable to the change of its value. In general, the outputs tend to be more and more smoothing when $\beta$ is increased. For the parameter $\alpha$, a similar trend is observed but with decreasing $\alpha$. It is worth noting that the second term in Eq. (5) and (6) plays a similar role as the data-fidelity term, thus $\alpha$ has a relatively small impact on global smoothing results but mainly adjusts the local details.

\textbf{Higher-order Regularization:} In general, the choice of $n$ of Eq. (5) is determined by the property of signals. Firstly, we show the smoothing results on a simulated surface consisting of a sloping plane and a $3^{rd}$ degree polynomial surface in Fig. 6. We use the Peak Signal-Noise-Ratio (PSNR) index to measure the smoothing results. In this case, a better result is produced when the sparsity-inducing regularization is employed on the $3^{rd}$ order gradient domain. In Fig. 7, we compare the difference on natural images. As we can see, the smoothing filter gives very similar results under the equivalent level of smoothness except in some tiny regions. This is maybe due to the fact that these areas of these rapid-increasing higher degree ($n \leq 3$) polynomial surfaces in natural images is too small and can be approximated as spikes.
 This can be also explained from the higher-order gradient distributions of natural images in Fig. 2, where the sparsity is enhanced when the order $n$ increases but the gap of sparsity between $1^{st}$ and $2^{nd}$ orders is much bigger than that between $2^{nd}$ and $3^{rd}$ orders. A similar trend can be observed in much higher cases. As a consequence, we claim that the polynomial smoothing surfaces with degree ($ n \leq 3$) not frequently occur in natural images, in other words, the supports (areas) of which is relatively smaller than that of the $n=2$ case. That is to say, it is usually enough for the natural images to impose the sparse regularization on second-order gradients.

\begin{figure}[!t
	\begin{center}
		\subfigure[Input]
		{\includegraphics[width=0.155\textwidth]{./imgs/order/input.jpg}}
		\subfigure[Second-order]
		{\includegraphics[width=0.155\textwidth]{./imgs/order/smoothing_order2_local.png}}
		\subfigure[Third-order]
		{\includegraphics[width=0.155\textwidth]{./imgs/order/smoothing_order3_local.png}}
	\end{center}
	\caption{Comparison of the smoothing results on a natural image with the sparse regularization on the second-order ($n=2$, PSNR = 33.96) and third-order ($n=3$, PSNR = 33.16) gradient domains, respectively.}
\end{figure}
\begin{figure*}[t]
	\begin{center}
		{\includegraphics[width=0.16\textwidth]{./imgs/flower0/img_src.png}}
		{\includegraphics[width=0.16\textwidth]{./imgs/flower0/img_bf.png}}
		{\includegraphics[width=0.16\textwidth]{./imgs/flower0/img_gf.png}}
		{\includegraphics[width=0.16\textwidth]{./imgs/flower0/img_wls.png}}
		{\includegraphics[width=0.16\textwidth]{./imgs/flower0/img_l0.png}}
		{\includegraphics[width=0.16\textwidth]{./imgs/flower0/img_sap.png}}
		{\includegraphics[width=0.16\textwidth]{./imgs/flower0/img_src_local.png}}
		{\includegraphics[width=0.16\textwidth]{./imgs/flower0/img_en_bf}}
		{\includegraphics[width=0.16\textwidth]{./imgs/flower0/img_en_gf.png}}
		{\includegraphics[width=0.16\textwidth]{./imgs/flower0/img_en_wls.png}}
		{\includegraphics[width=0.16\textwidth]{./imgs/flower0/img_en_l0.png}}
		{\includegraphics[width=0.16\textwidth]{./imgs/flower0/img_en_sap.png}}
		\subfigure[Input]
		{\includegraphics[width=0.16\textwidth]{./imgs/flower0/local_1.png}}
		\subfigure[BF~\cite{tomasi1998bilateral}]
		{\includegraphics[width=0.16\textwidth]{./imgs/flower0/local_2.png}}
		\subfigure[GF~\cite{he2012guided}]
		{\includegraphics[width=0.16\textwidth]{./imgs/flower0/local_3.png}}
		\subfigure[WLS~\cite{farbman2008edge}]
		{\includegraphics[width=0.16\textwidth]{./imgs/flower0/local_5.png}}
		\subfigure[$L_0$~\cite{xu2011image}]
		{\includegraphics[width=0.16\textwidth]{./imgs/flower0/local_6.png}}
		\subfigure[Ours]
		{\includegraphics[width=0.16\textwidth]{./imgs/flower0/local_7.png}}
	\end{center}
	\caption{Image enhancement with ($\times$2.5) detail layers. (a). Input image, (b). Bilateral filter ($r = 15, \delta =[7, 0.2^2]$)~\cite{tomasi1998bilateral}, (c). Guided filter ($r = 15, \epsilon = 0.15^2$)~\cite{farbman2008edge}, (d). WLS filter ($\alpha = 1.2, \lambda = 0.2$)~\cite{farbman2008edge}, (e). $L_0$ Smoothing ($\lambda = 0.01$)~\cite{xu2011image}, (f). Our result ($\alpha = 1.0, \beta =0.01$). All filters are configured to produce the similar level smoothing results (PSNR : (b)$\scriptsize{\sim}$(f): 28.83, 28.82, 28.67, 28.28, 28.36). }
\end{figure*} %

\section{Applications}

Similar to many existing filtering methods, it is always possible to use our semi-sparsity model in various signal processing fields. We describe some examples involving 2D images and 3D surfaces.

\subsection{2D Images}

\textbf{Details Manipulation:} \quad In many imaging systems, the captured images may suffer from degradation in details for the inappropriate focus or parameter-settings in image tone mapping, denoising and so on. In such cases, ``edge-aware'' filtering methods can be used to help recover or exaggerate the degraded details. In general, it is popular to use the smoothing filters to decompose an image into two components: a piece-wise smoothing base layer and a detail layer; then either of them can be remapped and combined for reconstructing better visual results. As shown in Fig. 8, we compare the results with the widely-used methods~\cite{farbman2008edge, fattal2002gradient, tomasi1998bilateral,  xu2011image}. To be fairness, we introduce a greedy search strategy for the parameter-settings to force each smoothing filter to give the similar level of smoothing results. Then the detail layer is obtained by subtracting the smoothing base-layer from the original image, which, in turn, are exaggerated by a scale ($\times 2.5$) factor and are recombined with the base-layers to form the enhanced images. As we can see, our semi-sparse smoothing filter helps to recover rich local details with less level artifacts compared with the cutting-edge methods~\cite{farbman2008edge,xu2011image}.  Besides, we also take a multi-scale  strategy for the detail enhancement, the configuration of which is also adopted from the WLS filter~\cite{farbman2008edge}. The different smoothing levels are applied to generate the base layers, and the detail layers are then classified into different bands by subtract one from the former smoothing result. As shown in Fig. 9,  the small, median and coarse levels are applied for producing smoothing base layers and the final merged result.

\textbf{HDR Image Compression:} \quad This technique arises from the limited capabilities of displays in processing HDR images.  This process is similar to the above image enhancement but on high dynamic range images. In this framework, the base layer is usually assumed to be piece-wise smoothing property with high dynamic range, which can be generated by a smoothing filter and need to be compressed into a low dynamic range (LDR) in non-linear way. The compressed output is then re-combined with the detail layer to reconstruct a new image for display or visual applications. The difficulties here is to keep a balance of the smoothness of piece-wise areas and sharp discontinuities to avoid halos for reasonable contrast maintenance in dynamic range compression. The weighted least square (WLS) filter~\cite{farbman2008edge} is one of powerful tools for its  capability of suppressing halos artifacts in HDR image compression. We here only replace the WLS filter with the other ones and keep the other configurations unchanging for fair comparison. As shown Fig. 10, WLS filter gives an LDR image with fine balance global contrast but lose some details in under-exposure region, which can be more or less recovered by  guided filter~\cite{ he2012guided} but it has a limitation in bringing out local tiny details; while $L0$ smoothing~\cite{xu2011image} tends to produce over-exaggerated local details, which is mainly caused by the stair-case fitting results in polynomial smoothing surfaces. In contrast, our semi-sparse filtering method provides a result with appropriate contrast in local and fine-balanced global image illumination, which can be also found in Fig. 11. We here employ a quantitative evaluation for HDR image compression, which is based on the Tone Mapped image Quality Index (TMQI)~\cite{yeganeh2012objective}. In this TMQI index, Structural Fidelity (SF) and Statistical Naturalness (SN) are considered to provide an objective quality assessment between the HDR images and the corresponding LDR output results.

\textbf{Image Stylization$\And$Abstraction} \quad This interesting task aims to produce different stylized results or non-photorealistic rendering (NPR) of an image. Here, we first take a smoothing strategy to remove the local details by using our semi-sparse filtering method. Then, as interpreted in~\cite{winnemoller2012xdog}, the main structure of images is extracted by applying the Difference-of-Gaussian (DoG) operator on smoothing outputs. Due to the simultaneously-fitting property of the semi-sparse model, it avoids to introducing extra edge lines in polynomial surfaces. Despite the isotropic property of the DoG operator,  it is also enough to produce visual-friendly results. As shown in~\cite{winnemoller2012xdog},
one can also produce different artistic effects such as hatching, woodcut, watercolor and so on. We show two simple examples: pen and color-pencil drawing effects in Fig. 12. One can also combine our semi-sparse smoothing filter model with more complex configurations as illustrated in~\cite{kang2008flow, lu2012combining} to produce more reasonable and aesthetic stylized results.

\begin{figure*}[!t]
	\begin{center}
		\subfigure[Input]
		{\includegraphics[width=0.19\textwidth]{./imgs/enhancer/orig.png}}
		\subfigure[Small]
		{\includegraphics[width=0.19\textwidth]{./imgs/enhancer/fine.png}}
		\subfigure[Median]
		{\includegraphics[width=0.19\textwidth]{./imgs/enhancer/medium.png}}
		\subfigure[Coarse]
		{\includegraphics[width=0.19\textwidth]{./imgs/enhancer/coarse.png}}
		\subfigure[Merged Result]
		{\includegraphics[width=0.19\textwidth]{./imgs/enhancer/combine.png}}
	\end{center}
	\caption{Multi-scale details manipulation. Three scales: (b) small ($\alpha = 0.01, \beta =0.001$) and (c) median ($\alpha = 0.1, \beta =0.005$), (d) coarse ($\alpha = 1.0, \beta =0.02$), and (e) is obtained by scaling ($\times 2$) detail-layers.}
\end{figure*} %
\begin{figure*}[!h]
	\begin{center}
		\subfigure[WLS~\cite{farbman2008edge}]
		{\includegraphics[width=0.24\textwidth]{./imgs/desk/desk_wls.png}}
		\subfigure[Guided Filter~\cite{he2012guided}]
		{\includegraphics[width=0.24\textwidth]{./imgs/desk/desk_gf.png}}
		\subfigure[$L_0$ Smoothing~\cite{xu2011image}]
		{\includegraphics[width=0.24\textwidth]{./imgs/desk/desk_l0_01_2_log5.png}}
		\subfigure[Our Result ($n=2$)]
		{\includegraphics[width=0.24\textwidth]{./imgs/desk/desk_sp.png}}
	\end{center}
	\caption{HDR image tone mapping. From left to right: (a) WLS filter ($\alpha = 2.0, \lambda = 1.2$)~\cite{farbman2008edge}, (b) Guided filter ($r =15, \epsilon = 0.1^2$)~\cite{he2012guided}, (c) $L_0$ smoothing ($\lambda = 0.1, \kappa =2.0$)~\cite{xu2011image}, and (d) Our result ($\alpha = 0.5, \beta =0.01, \kappa =1.05$). TMQI: 0.9295, 0.9305, 0.8965, 0.9562.}
\end{figure*} %

\begin{figure*}[!h
	\begin{center}
		\subfigure[Bilateral Filter (BF)~\cite{tomasi1998bilateral}]
		{\includegraphics[width = 0.32\textwidth]{./imgs/hdr/vinesunset_bf_local.png}}
		\subfigure[Guided Filter (GF)~\cite{he2012guided}]
		{\includegraphics[width = 0.32\textwidth]{./imgs/hdr/vinesunset_gf_local.png}}
		\subfigure[WLS~\cite{farbman2008edge}]
		{\includegraphics[width = 0.32\textwidth]{./imgs/hdr/vinesunset_wls_local.png}}
		\subfigure[L0 Minimization~\cite{xu2011image}]
		{\includegraphics[width = 0.32\textwidth]{./imgs/hdr/vinesunset_l0_local.png}}
		\subfigure[Relative Total Variation (RTV)~\cite{xu2012structure}]
		{\includegraphics[width = 0.32\textwidth]{./imgs/hdr/vinesunset_rtv_local.png}}
		\subfigure[Our Result]
		{\includegraphics[width = 0.32\textwidth]{./imgs/hdr/vinesunset_sap_local.png}}
	\end{center}
	\caption{An example of HDR image compression. We use different  filtering methods to provide the base layers, and the results are then obtained by combining these compressed base layers with detailed layers. We keep the same parameters except for the filters: (a) Bilateral filter ($\sigma_s = 15, \sigma_r =0.2$)~\cite{tomasi1998bilateral},  (b) Guided filter ($r = 15, \epsilon = 0.1^2$)~\cite{he2012guided}, (c) WLS filter ($\alpha = 1.2, \lambda = 0.5$)~\cite{farbman2008edge}, (d) L0 smoothing ($\lambda = 0.02$)~\cite{xu2011image}, (e) Relative total variation (RTV) ($\lambda = 0.02, \sigma =0.02$)~\cite{xu2012structure}, (f) Our result ($\alpha = 1.0, \beta =0.01$). TMQI: 0.7823, 0.7837, 0.7849, 0.7723, 0.7864, 0.7918.}
\end{figure*}

\begin{figure}[!t]
	\begin{center}
		\subfigure[Input]
		{\includegraphics[width=0.155\textwidth]{./imgs/Stylization/input.png}}
		\subfigure[Pen Drawing]
		{\includegraphics[width=0.155\textwidth]{./imgs/Stylization/pen.png}}
		\subfigure[Color-pencil]
		{\includegraphics[width=0.155\textwidth]{./imgs/Stylization/color_pencil.png}}
	\end{center}
	\caption{Image stylization.  From left to right: (a) Input image, (b) Pen drawing, and (c) Color pencil effect.}
\end{figure} %

\subsection {Mesh Denoising}
Similar to many filtering methods~\cite{fleishman2003bilateral, he2013mesh, sun2007fast, zhang2018static, zhang2015guided,  zheng2010bilateral}, we illustrate that our semi-sparse model can be extended to 3D geometry.  As depicted in~\cite{he2013mesh}, an $L_0$ minimization scheme is deduced on meshes with the sparse regularization of the first-order difference of edges. We extend this model into higher-order cases and show it enables us to smooth 3D meshes without causing stair-case artifacts as well. For the 3D meshes, the $2^{nd}$ order gradient operator corresponds to the curvature of local geometry. The one more than $2^{nd}$ is seldom used in practice. We here also choose $n=2$ in Eq. (5) for the simplicity of the model and the purpose of saving the computational cost.

\begin{figure}[!h]
	\begin{center}
		\subfigure[]
		{\includegraphics[width=0.11\textwidth]{./imgs/Mesh/op1.png}}
		\subfigure[]
		{\includegraphics[width=0.12\textwidth]{./imgs/Mesh/op2.png}}
		\subfigure[]
		{\includegraphics[width=0.125\textwidth]{./imgs/Mesh/op3.png}}
		\subfigure[]
		{\includegraphics[width=0.115\textwidth]{./imgs/Mesh/op4.png}}
	\end{center}
	\caption{The illustration of discrete operators on meshes.}
\end{figure}
\begin{figure*}[!b]
	\begin{center}
		\subfigure[Noisy Input]
		{\includegraphics[width=0.165\textwidth]{./imgs/Fandisk/Fandisk_Input_local.png}}
		\subfigure[BF~\cite{fleishman2003bilateral}]
		{\includegraphics[width=0.165\textwidth]{./imgs/Fandisk/Fandisk_BF_local.png}}
		\subfigure[SRML~\cite{sun2007fast}]
		{\includegraphics[width=0.16\textwidth]{./imgs/Fandisk/Fandisk_SRML_local.png}}
		\subfigure[$L_0$~\cite{he2013mesh}]
		{\includegraphics[width=0.16\textwidth]{./imgs/Fandisk/Fandisk_L0_local.png}}
		\subfigure[GNF~\cite{zhang2015guided}]
		{\includegraphics[width=0.16\textwidth]{./imgs/Fandisk/Fandisk_GF_local.png}}
		\subfigure[Ours]
		{\includegraphics[width=0.16\textwidth]{./imgs/Fandisk/Fandisk_Ours_local.png}}
	\end{center}
	\caption{3D mesh denoising, surface corrupted by Gaussian noise in random directions with standard deviation $\sigma = 0.2\bar{L}_e$ ($\bar{L}_e$ is the average length of edges). (a) Noisy input, (b) Bilateral filter (BF)~\cite{fleishman2003bilateral},  (c) SRML~\cite{sun2007fast}, (d) $L_0$ minimization~\cite{he2013mesh},(e) Guided normal filter (GNF)~\cite{zhang2015guided}, (f) Our result.}
\end{figure*} %

For completeness, we introduce a measure of mean curvature as interpreted in~\cite{oztireli2009feature}. Let $M$ be a non-degenerate triangulated surface with the vertices, edges, and triangles denoted as $v_{i(i=0,1,2\cdots,I)}$, $e_{j(j=0,1,2\cdots,J)}$ and $\tau_{k(k=0,1,2\cdots,K)}$, respectively. We introduce the relative orientation of an edge $e$ to a triangle $\tau$, which is denoted by $s(e, \tau)$. We assume that all triangles have the counter-clockwise orientations and all edges are with randomly chosen fixed orientations. If an edge $e$ and a triangle $\tau$ have the same
orientation, then $s(e, \tau) = 1$; otherwise, $s(e, \tau) = -1$.

Considering the barycentric coordinates of triangles in Fig. 10 (a), the outward-facing normal of triangle $\tau$ is denoted $e_{\tau}^\bot$. By definition, it is perpendicular to edge $e=(v_{\scaleto{i}{4pt}}\!-\!v_{\scaleto{i+1}{4pt}})$ and satisfies~\cite{pinkall1993computing},
\begin{equation}
	\small
	\begin{aligned}
		e_{\tau}^\bot \!=\! (v_{\scaleto{i}{4pt}}\!-\!v_{\scaleto{i+1}{4pt}})^\bot \!=\! \cot(\theta_1)(v_{\scaleto{i+1}{4pt}}\!-\!v_{\scaleto{0}{4pt}})\!+\!\cot(\theta_2)(v_{\scaleto{i}{4pt}}\!-\!v_{\scaleto{0}{4pt}})
	\end{aligned}
	\label{f2}
\end{equation}
where $\theta_1$ and $\theta_2$ are the angles as shown in Fig. 13 (a). We here define the first-order difference operators of edges based on the outward-facing normals. As depicted in Fig. 13 (b), two triangles $\tau^+$ and $\tau^-$ share the common edge $e$ and the orientation of $e$ is assumed to be consistent with $\tau^+$.

Let $u^e_{\tau}$ (or, $u_{\tau}$ for simplicity) be a vector restricted on edge $e$ associated to $u_{\tau}$, for example, $u^e_{\tau}=-e_{\tau}^\bot s(e,\tau)$. That's to say, it is opposite to the outward-facing normal $e_{\tau}^\bot$ with consistent orientation of edge $e$. Then,  it is natural to define the first-order difference over the edge $e$ as,
\begin{equation}
	\small
	\begin{aligned}
		\left[u\right]_e
		& =  u_{\tau^{+}}s(e, {\tau}^+)  \!+ \! u_{\tau^{-}}s(e, {\tau}^-)  \!= \!u_{\tau^{+}} \!-\! u_{\tau^{-}}.
	\end{aligned}
	\label{f2}
\end{equation}
It is clear that $\left[u\right]_e$ is equal to the sum of outward-facing normals of the  edges ({\color{red} red}) in Fig. 13 (b). This result is consistent with the conclusion in~\cite{he2013mesh}. Clearly, it degenerates to $0$, if $\tau^{+}$ and $\tau^{-}$ are in the same plane.

By analogy, we can define the second-order difference operator to construct our high order regularization models. Let $l$ be the line connecting the barycenter and one vertex of the triangle $\tau$ as depicted in Fig. 13 (c),  we then define the difference over the line $l$ as,
\begin{equation}
	\small
	\begin{aligned}
		\left[\left[u\right]\right]_l
		&= \left[u\right]_{e^+}s(e^+, {\tau}^{+})  \!+ \! \left[u\right]_{e^-}s(e^-,{\tau}^{-})\\
		&= (u_{\tau^{+}} \!-\! u_{\tau})-(u_{\tau} \!-\! u_{\tau^{-}}) = u_{\tau^{+}} \!-\! 2u_{\tau}\!+\! u_{\tau^{-}}.
	\end{aligned}
	\label{f2}
\end{equation}
We can see that the definition of $\left[\left[u\right]\right]_l$ is invariant under the choice of  orientation of edges. With above definitions, it is easy to compute the gradient operator $\nabla u \vert_e$ and Laplace operator $\Delta u\vert_e$ over the edge $e$. Using the divergence theorem~\cite{botsch2010polygon}, we have $\nabla u\vert_e = \left[u\right]_e$ and $\Delta u\vert_e$,
\begin{equation}
	\small
	\begin{aligned}
		\Delta u\vert_e = \sum_{l \in H(e)} \frac{L_l}{L_e} \left[\left[u\right]\right]_l s(e, {\tau}_l),
	\end{aligned}
	\label{f2}
\end{equation}
where $\tau_l$ is the triangle containing line $l$, ${L_l}$ and ${L_e}$ are length of line $l$ and edge $e$, $H(e)$ is the set of lines ({\color{green} green}) associated with the edge $e$ as shown in Fig. 13 (d). The reader is referred to~\cite{botsch2010polygon, crane2018discrete, liu2019triangulated, liu2020mesh, meyer2003discrete} for more details.
\begin{figure*}[!t
	\begin{center}
		\begin{tabular}{c}
			\begin{minipage}[t]{\textwidth}
				{\includegraphics[width=0.24\textwidth]{./imgs/man/noisy_local.png}}
				{\includegraphics[width=0.24\textwidth]{./imgs/man/L0_local.png}}
				{\includegraphics[width=0.24\textwidth]{./imgs/man/our_local.png}}
				{\includegraphics[width=0.24\textwidth]{./imgs/man/src_local.png}}
			\end{minipage}
		\end{tabular}
		\begin{tabular}{c}
			\begin{minipage}[t]{\textwidth}
				\subfigure[Noisy Input]
				{\includegraphics[width=0.24\textwidth]{./imgs/elephant/noisy_local.png}}
				\subfigure[L0 Minimization~\cite{he2013mesh}]
				{\includegraphics[width=0.24\textwidth]{./imgs/elephant/L0_local.png}}
				\subfigure[Our Result]
				{\includegraphics[width=0.24\textwidth]{./imgs/elephant/our_local.png}}
				\subfigure[Ground Truth]
				{\includegraphics[width=0.24\textwidth]{./imgs/elephant/src_local.png}}
			\end{minipage}
		\end{tabular}
	\end{center}
	\caption{Mesh denoising results. The first row shows a 3D mesh surface corrupted by Gaussian noise in random directions with standard deviation $\delta$= 0.2 average length of edges.}
\end{figure*}
\begin{figure}[!t]
	\begin{center}
		{\includegraphics[width=0.115\textwidth]{./imgs/Mesh/pyramid_input.png}}
		{\includegraphics[width=0.115\textwidth]{./imgs/Mesh/pyramid_cnr.png}}
		{\includegraphics[width=0.115\textwidth]{./imgs/Mesh/pyramid_l0.png}}
		{\includegraphics[width=0.115\textwidth]{./imgs/Mesh/pyramid_our.png}}
		\subfigure[]
		{\includegraphics[width=0.115\textwidth]{./imgs/Mesh/boy_input.png}}
		\subfigure[]
		{\includegraphics[width=0.115\textwidth]{./imgs/Mesh/boy_cnr.png}}
		\subfigure[]
		{\includegraphics[width=0.115\textwidth]{./imgs/Mesh/boy_l0.png}}
		\subfigure[]
		{\includegraphics[width=0.115\textwidth]{./imgs/Mesh/boy_our.png}}
	\end{center}
	\caption{Mesh denoising on scanned surfaces. (a) Noisy input, (b) Cascaded normal regression~\cite{wang2016mesh}, (c) $L_0$ minimization~\cite{he2013mesh}, (d) Our result.}
\end{figure}

Once $\nabla u\vert_e $ and $\Delta u\vert_e$ are obtained, it is easy to substitute them into Eq. (6) and solve the problem based on the Alg.1\footnote{We here solve Eq. (7) using the Gauss-Seidel method~\cite{saad2003iterative}, as the FFT solver is not unavailable on triangulated 3D meshes.}. We here also compare our semi-sparsity model with the feature-preserving denoising~\cite{sun2007fast}, bilateral filter~\cite{fleishman2003bilateral}, guided mesh filter~\cite{zhang2015guided}  and $L_0$ minimization~\cite{he2013mesh}. In Fig. 14, the original surface contains sparse sharpening features (corners and edges) and polynomial smoothing surfaces. Obviously, bilateral filter removes noise in smoothing areas but blurs sharp features; the SRML and GNF methods produce much better results with distortions in some edges, and $L_0$ minimization retains sharp features but cause stair-case artifacts in smoothing regions. While, our semi-sparse smoothing method produces noticeable better results in both sharpening and polynomial-smoothing regions. This is further demonstrated by the smoothing results in Fig. 15 and the real scanned surfaces in Fig. 16.

\section{Conclusion}
This paper has described a new semi-sparse smoothing method, which is formed as a general extension of $L_0$-norm minimization in higher-order gradient domains. We  demonstrate the virtue of preserving the singularities and fitting ability in polynomial surfaces. We show its benefits and advantages with a series of experiments in computer vision and graphics. The drawback of our method is the high computational cost in solving the iterative optimization problem, which could be partially solved by resorting to some approximate solutions like separable solvers. We leave
this for future studies.

\ifCLASSOPTIONcompsoc

  \section*{Acknowledgments}
  	The authors would like to thank to the on-line resources. This work was supported in part by the FWO Odysseus Project by the Leverhulme Grant RPG-2017-151 and EPSRC grant EP/R003025/1.
\else

  \section*{Acknowledgment}
\fi

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi