## 2 Materials and Methods
### Image Restoration  
With the estimated blur kernel \(\@vec{k}\) from the above method, the latent image can be revealed by nonblind deconvolution algorithms, such as Wiener filtering and the Lucy-Richardson method, which were introduced decades ago. These deconvolution algorithms are classic and have been widely implemented in many fields. However, they are sensitive to noise, especially when the blur kernels are estimated incorrectly, and the ringing effect is unavoidable. Therefore, inspired by the kernel estimation discussed above, the sparse representation techniques can also be applied in the restoration of the latent image. Based on the statistical character of natural images in recent studies, the distributions of image gradient domains tend to have heavy tails and sharp tops, which means that most of their distribution is on small values, but the probability of large values is much greater compared to the Gaussian distribution [14; 15; 16; 17; 18; 19]. Therefore, the hyper-Laplacian distribution-based prior can be used as a regularization term in the cost function as follows:  
\[\@vec{f}=\arg\min_{f}\parallel\@vec{k}\otimes\@vec{f}-\@vec{g}\parallel_{2}^{ 2}+\alpha_{2}\@vec{\parallel}\left(\
abla\@vec{f},\@vec{V}\@vec{g}\right) \parallel_{p} \tag{17}\]  
where \(p\) is derived from the hyper-Laplacian model, which signifies the slope of the exponential term of the density function. In Equation (17), it becomes the constraint term as a quasi-norm form, \(\|\cdot\|_{p}\), as shown above. \(\
abla_{x}\@vec{f}\) and \(\
abla_{y}\@vec{f}\) are the image gradient domains in the horizontal (\(x\)-axis), and vertical (\(y\)-axis), respectively. The optimization can be solved by half-quadratic penalty method similar to that of Equation (7). With auxiliary variable \(\@vec{v}=\left(v_{x},\ v_{y}\right)\), the alternating scheme is shown as follows:  
\[\left\{ \begin{array}{l}\hat{\
u}=\arg\min_{\@vec{v}}\parallel v\parallel_{p }+\frac{\beta_{2}}{\alpha_{2}}\@vec{\parallel}\left(\
abla_{x}\@vec{f},\ \
abla_{y}\@vec{f}\right)-\left(v_{x},\ v_{y}\right)\parallel_{2}^{2}\\ \hat{\gamma}=\arg\min_{f}\parallel\@vec{k}\otimes\@vec{f}-\@vec{g}\parallel_{ 2}^{2}+\beta_{2}\@vec{\parallel}\left(\
abla_{x}\@vec{f},\ \
abla_{y}\@vec{f}\right)-\left(v_{x},\ v_{y}\right)\parallel_{2}^{2}\end{array} \right. \tag{18}\]  
The solutions for the \(\hat{\
u}\)-subproblem vary with different \(p\) values. In particular, when \(p=1\), the analytical solution can be obtained in Equation (19) by the soft thresholding algorithm [17] as follows:  
\[\hat{\
u}=sgn(\@vec{v})\max\{\@vec{0},\ |\@vec{v}|-T\} \tag{19}\]where \(T=\alpha_{2}/(2\beta_{2})\) is the threshold and \(sgn(\cdot)\) represents the signum function. It should be noted that the solutions in the horizontal and vertical orientation have the same form in the following Equation (20)-(22), therefore the subscripts related to orientations are omitted for briefness and convenience. In general, for an arbitrary \(p\), the minimization problem can be solved by setting the derivative of the \(\hat{v}\)-subproblem to zero. Therefore, for the \(i\)th element of \(v\) and \(\
abla f\) in horizontal orientation, we get:  
\[p|v_{i}|^{p-1}sgn(v_{i})+\frac{2\beta_{2}}{\alpha_{2}}\Big{(}v_{i}-\
abla f_{ i}\Big{)}=0 \tag{20}\]  
In particular, when \(p=0.5\), the \(\hat{v}\)-subproblem in Equation (18) can be simplified to a cubic function:  
\[v_{i}^{3}-2\
abla f_{i}v_{i}^{2}+\
abla f_{i}^{2}\mathbf{u}_{i}-\Bigg{(}\frac{ \alpha_{2}}{4\beta_{2}}\Bigg{)}^{2}sgn(v_{i})=0 \tag{21}\]  
which can be solved using [PERSON]'s formula [31]. When \(p=2/3\), the \(\hat{v}\)-subproblem in Equation (18) can be expanded as a quartic function:  
\[v_{i}^{4}-3\
abla f_{i}v_{i}^{3}+3\
abla f_{i}^{2}v_{i}^{2}-\
abla f_{i}^{2} v_{i}-\Bigg{(}\frac{\alpha_{2}}{3\beta_{2}}\Bigg{)}^{3}=0 \tag{22}\]  
which can be solved using [PERSON]'s and [PERSON] solutions [31]. For some special cases where \(p>1\), analytical solutions can be referred to in [32].