## 3 Proposed Methods
### \(\ell_{0}\)-Smoothing Based on Box-Type Gradient Constraint
#### 3.2.2 Optimization
\begin{cases}\mathbf{A}=\lambda\mathbf{I}+\gamma_{1}^{-1}\mathbf{ D}^{\top}\mathbf{D}+\gamma_{2}^{-1}\mathbf{D}^{\top}\mathbf{P}^{\top}\mathbf{PD}+ \gamma_{3}^{-1}\mathbf{D}^{\top}\mathbf{M}^{\top}\mathbf{MD},\\ \mathbf{q}=\lambda\mathbf{s}_{\text{in}}+\gamma_{1}^{-1}\mathbf{D}^{\top}( \mathbf{z}_{1}-\mathbf{b}_{1})+\gamma_{2}^{-1}\mathbf{D}^{\top}\mathbf{P}^{ \top}(\mathbf{z}_{2}-\mathbf{b}_{2})+\gamma_{3}^{-1}\mathbf{D}^{\top} \mathbf{M}^{\top}(\mathbf{z}_{3}-\mathbf{b}_{3}),\end{cases} \]  
where \(\mathbf{A}\) is a matrix of size \(\{3N\times 3N\}\), \(\mathbf{q}\) is a \(3N\)-dimensional vector, and \(\mathbf{I}\in\mathbb{R}^{3N\times 3N}\) is an identity matrix. The optimal solution is obtained by the inverse problem \(\mathbf{s}^{*}=\mathbf{A}^{-1}\mathbf{q}\). If we set \(\gamma_{3}=\gamma_{2}\), the third term of the right side in the formula \(\mathbf{A}\) can be rewritten as \(\gamma_{2}^{-1}\mathbf{D}^{\top}\mathbf{D}=\gamma_{2}^{-1}\mathbf{D}^{\top} \mathbf{P}^{\top}\mathbf{PD}+\gamma_{3}^{-1}\mathbf{D}^{\top}\mathbf{M}^{ \top}\mathbf{MD}\). Since the boundary condition of \(\mathbf{D}\) is periodic and the matrix \(\mathbf{A}\) is a block circulant matrix with circulant blocks (BCCB), it is diagonalized by the 2D fast discrete Fourier transform (2 DFFT). Thus, the inverse problem of (20) can be calculated as  
\[\mathbf{s}^{*}=\mathbf{F}^{*}\Big{\{}\lambda\mathbf{I}+(\gamma_{1}^{-1}+ \gamma_{2}^{-1})\mathbf{\Sigma}\Big{\}}^{-1}\mathbf{F}\mathbf{q}, \tag{21}\]  
where \(\mathbf{F}\) and \(\mathbf{F}^{*}\) are the 2D-FFT and its inverse matrices, respectively, and \(\mathbf{\Sigma}\) is the diagonal matrix with its entries being the Fourier-transformed Laplacian filter kernel. Therefore, its inversion is reduced to entry-wise division. For the update of \(\mathbf{z}_{1}\), we need to compute the pseudo-proximity operator of the \(\ell_{0,1}\)-norm. By dividing \(\mathbf{z}_{1}\) to sub-vectors as \(\mathbf{z}_{1}:=[\mathbf{z}_{1}^{\mathcal{G}_{1}\top}\dots\mathbf{z}_{1}^{ \mathcal{G}_{N}\top}]^{\top}\), in which \(\mathcal{G}_{n}\)\((n=1,\dots,N)\) is introduced in Section 3.1.2 and \(\mathbf{z}_{1}^{\mathcal{G}_{n}}\) can be regarded as the auxiliary variable for \((\mathbf{D}\mathbf{s})^{\mathcal{G}_{n}}\) (see (15)), the pseudo-proximity operator equals to a group hard-thresholding operation for each sub-vector [48; 61]. Thus, the optimal solution of the \(n\)-th sub-vector \(\mathbf{z}_{1}^{\mathcal{G}_{n}}\) is obtained by  
\[\mathbf{z}_{1}^{\mathcal{G}_{n}*}\leftarrow\text{prox}_{\gamma_{1}\|\cdot\|_{0 _{1}}}(\mathbf{d}_{n})=\left\{ \begin{array}{ll}\mathbf{d}_{n},&\text{if }\sqrt{\sum_{m=1}^{6}d_{n,m}^{2}}\geq\sqrt{2\gamma_{1}},\\ \mathbf{0},&\text{otherwise},\end{array} \right. \tag{22}\]  
where \(d_{n,m}\) is the \(m\)-th element of \(\mathbf{d}_{n}:=(\mathbf{D}\mathbf{s}+\mathbf{b}_{1})^{\mathcal{G}_{n}}\)\((\in\mathbb{R}^{6})\). To update \(\mathbf{z}_{2}\) and \(\mathbf{z}_{3}\), we need to compute the proximity operator of the indicator functions \(\iota_{\mathcal{C}_{[\mathbf{g}^{+}]}}\) and \(\iota_{\mathcal{C}_{[\mathbf{g}^{-},\mathbf{g}^{+}]}}\). By (9), these are simply projections onto each set, respectively, i.e., for \(i=1,\ldots,m_{+}\) by  
\[\mathbf{z}_{2}\leftarrow\left[\mathcal{P}_{\mathcal{C}_{[\mathbf{g}^{+}]}}( \mathbf{PD}\mathbf{s}+\mathbf{b}_{2})\right]_{i}=\left\{ \begin{array}{ll} 0,&\text{if }(\mathbf{PD}\mathbf{s}+\mathbf{b}_{2})_{i}<0,\\ g_{i}^{+},&\text{if }(\mathbf{PD}\mathbf{s}+\mathbf{b}_{2})_{i}>g_{i}^{+},\\ (\mathbf{PD}\mathbf{s}+\mathbf{b}_{2})_{i},&\text{if }0\leq(\mathbf{PD} \mathbf{s}+\mathbf{b}_{2})_{i}\leq g_{i}^{+},\end{array} \right. \tag{23}\]  
and for \(j=1,\ldots,m_{-}\) by  
\[\mathbf{z}_{3}\leftarrow\left[\mathcal{P}_{\mathcal{C}_{[\mathbf{g}^{-}, \mathbf{g}]}}(\mathbf{MD}\mathbf{s}+\mathbf{b}_{3})\right]_{j}=\left\{