SURE-Type Functionals as Criteria for Parametric PSF Estimation

Abstract

Point spread function (PSF) estimation plays an important role in blind image deconvolution. This paper proposes two novel criteria for parametric PSF estimation, based on Stein’s unbiased risk estimate (SURE), namely, prediction-SURE and its variant. We theoretically prove the SURE-type functionals incorporating exact (complementary) smoother filtering as the valid criteria for PSF estimation. We also provide the theoretical error analysis for the regularizer approximations, by which we show that the proposed frequency-adaptive regularization term yields more accurate PSF estimate than others. In particular, the proposed SURE-variant enables us to avoid estimation of noise variance, which is a key advantage over the traditional SURE-like functional. Finally, we propose an efficient algorithm for the minimizations of the criteria. Not limited to the examples we show in this paper, the proposed SURE-based framework has a great potential for other imaging applications, provided the parametric PSF form is available.

Appendices

Appendix 1: Proof of Theorem 3

Proof

Based on Lemma 1, the prediction-SURE can be written as

$$\begin{aligned} \epsilon =\frac{1}{N} \big \Vert \mathbf {U}\mathbf {y}-\mathbf {y}\big \Vert ^2 +\frac{2\sigma ^2}{N} {\mathrm {div}}_{\mathbf {y}}\big (\mathbf {U}\mathbf {y}\big ) -\sigma ^2. \end{aligned}$$

We omit the subscript $\mathbf {U}_{\mathbf {H},\lambda }$ as $\mathbf {U}$ for brevity in this proof.

Now, we are going to compute the divergence term—${\mathrm {div}}_{\mathbf {y}}\big (\mathbf {U}\mathbf {y}\big )$. Notice that the matrix $\mathbf {U}$ is the function of $\mathbf {y}$, hence, ${\mathrm {div}}_{\mathbf {y}}\big (\mathbf {U}\mathbf {y}\big )\ne {\mathrm {Tr}}(\mathbf {U})$.

First, by definition of divergence term, we have

$$\begin{aligned} {\mathrm {div}}_{\mathbf {y}}\big (\mathbf {U}\mathbf {y}\big ) ={\mathrm {Tr}}(\mathbf {U})+\bigg ( \underbrace{\frac{\partial \mathrm {diag}(\mathbf {U})}{\partial \mathbf {y}}}_{\alpha }\bigg )^\mathrm {T}\mathbf {y}, \end{aligned}$$

(28)

where the vector $\mathrm {diag}(\mathbf {U})\in \mathbf {R}^N$ consists of the diagonal element $U_{n,n}$ of matrix $\mathbf {U}$:

$$\begin{aligned} \frac{\partial \mathrm {diag}(\mathbf {U})}{\partial \mathbf {y}} =\left[ \frac{\mathrm {d}U_{1,1}(\mathbf {y})}{\mathrm {d}y(0)}, \frac{\mathrm {d}U_{2,2}(\mathbf {y})}{\mathrm {d}y(1)}, \ldots , \frac{\mathrm {d}U_{N,N}(\mathbf {y})}{\mathrm {d}y(N-1)}\right] ^\mathrm {T}\end{aligned}$$

Notice that under periodic boundary condition, $\mathbf {U}$ is circulant matrix, whose diagonal element is a constant given by

$$\begin{aligned}&U_{n,n}(\mathbf {y})=\underbrace{u_N(0)}_\text {filter}\\&\quad =\frac{1}{N}\sum _{k=0}^{N-1} \underbrace{\frac{|H(k)|^2}{|H(k)|^2 +\frac{\lambda }{|Y(k)|^2}}}_{U(k)},\qquad \mathrm {for}\ n=0,1,...,N-1. \end{aligned}$$

The second equality is from inverse Fourier transform. Now, we consider the n-th element of $\alpha$:

$$\begin{aligned}&\alpha _n=\frac{\partial U_{n,n}(\mathbf {y})}{\partial y(n)} =\frac{\partial U_{n,n}(\mathbf {y})}{\partial |Y(k)|^2} \cdot \frac{\partial |Y(k)|^2}{\partial y(n)} \\&\quad =\frac{1}{N}\sum _{k=0}^{N-1} \frac{\partial \frac{|H(k)|^2}{|H(k)|^2 +\frac{\lambda }{|Y(k)|^2}}}{\partial |Y(k)|^2} \cdot \frac{\partial |Y(k)|^2}{\partial y(n)}, \end{aligned}$$

where

$$\begin{aligned} \frac{\partial \displaystyle \bigg (\frac{|H(k)|^2}{|H(k)|^2 +\frac{\lambda }{|Y(k)|^2}}\bigg )}{\partial |Y(k)|^2} =\frac{|H(k)|^2\lambda }{\Big (|H(k)|^2 +\frac{\lambda }{|Y(k)|^2}\Big )^2\cdot |Y(k)|^4} \end{aligned}$$

and for any fixed n:

$$\begin{aligned}&\frac{\partial |Y(k)|^2}{\partial y(n)}= \frac{\partial \displaystyle \bigg |\sum _{l=0}^{N-1}y(l) e^{-j\frac{2\pi kl}{N}}\bigg |^2}{\partial y(n)}\\&\quad = \frac{\partial \displaystyle \Bigg ( \bigg |\sum _{l=0}^{N-1}y(l)\cos \frac{2\pi kl}{N}\bigg |^2 + \bigg |\sum _{l=0}^{N-1}y(l)\sin \frac{2\pi kl}{N}\bigg |^2 \Bigg )}{\partial y(n)}, \end{aligned}$$

where the two terms are

$$\begin{aligned}&\frac{\partial \displaystyle \bigg (\sum _{l=0}^{N-1}y(l)\cos \frac{2\pi kl}{N}\bigg )^2}{\partial y(n)}\\&\quad =2\cos \frac{2\pi nk}{N}\bigg ( \sum _{l=0}^{N-1}y(l)\cos \frac{2\pi lk}{N}\bigg ) \\&\quad =2\cos \frac{2\pi nk}{N}{\mathcal {R}}\{Y(k)\} \end{aligned}$$

and

$$\begin{aligned}&\frac{\partial \displaystyle \bigg (\sum _{l=0}^{N-1}y(l)\sin \frac{2\pi kl}{N}\bigg )^2}{\partial y(n)}\\&=2\sin \frac{2\pi nk}{N}\bigg ( \sum _{l=0}^{N-1}y(l)\sin \frac{2\pi lk}{N}\bigg ) \\&=-2\sin \frac{2\pi nk}{N}{\mathcal {I}}\{Y(k)\} \end{aligned}$$

Hence, $\alpha _n$ becomes

$$\begin{aligned} \alpha _n= & {} \frac{2}{N}\sum _{k=0}^{N-1} \underbrace{\frac{|H(k)|^2\lambda }{\Big (|H(k)|^2 +\frac{\lambda }{|Y(k)|^2}\Big )^2\cdot |Y(k)|^4}} _{V(k): \text {real-valued}}\\&\quad \cdot \bigg ({\mathcal {R}}\big \{Y(k)\big \}\cos \frac{2\pi nk}{N} \\&-{\mathcal {I}}\big \{Y(k)\big \}\sin \frac{2\pi nk}{N}\bigg )\\= & {} \frac{2}{N}\sum _{k=0}^{N-1} {\mathcal {R}}\big \{\underbrace{V(k)Y(k)}_{G(k)}\big \} \cos \frac{2\pi nk}{N} \\&-\frac{2}{N}\sum _{k=0}^{N-1} {\mathcal {I}}\big \{\underbrace{V(k)Y(k)}_{G(k)}\big \} \sin \frac{2\pi nk}{N}\\= & {} 2{\mathcal {R}}\{g(n)\}, \end{aligned}$$

where g(n) is the inverse Fourier transform of $G(k)=V(k)Y(k)$: $g(n)=\frac{1}{N}\sum _{k=0}^{N-1} G(k)e^{j\frac{2\pi nk}{N}}$. Since V(k) is real-valued, and Y(k) is Hermitian symmetric, g(n) is real-valued. Thus, we have $\alpha _n=2{\mathcal {R}}\{g(n)\}=2g(n)$.

Finally, the second term of (28) becomes

$$\begin{aligned} \alpha ^\mathrm {T}\mathbf {y}=2(\mathbf {F}G)^\mathrm {T}\mathbf {F}Y =2G^\mathrm {T}\mathbf {F}^\mathrm {T}\mathbf {F}Y=2G^\mathrm {T}Y\\ =\sum _{k=0}^{N-1}V(k)Y^*(k)Y(k) =\sum _{k=0}^{N-1}V(k)|Y(k)|^2\\ =\sum _{k=0}^{N-1}\underbrace{ \frac{|H(k)|^2\lambda }{\Big (|H(k)|^2 +\frac{\lambda }{|Y(k)|^2}\Big )^2\cdot |Y(k)|^2}}_{Q(k)}, \end{aligned}$$

where $\mathbf {F}$ is 1-D DFT matrix. The proof is completed. $\square$

Appendix 2: Proof of Corollary 1

Proof

The proof is very similar to that of Theorem 3 (see Appendix 1). The only difference from one-dimensional case is how to compute $\frac{\partial |Y(k_1,k_2)|^2}{\partial y(n_1,n_2)}$:

$$\begin{aligned} \frac{\partial |Y(k_1,k_2)|^2}{\partial y(n_1,n_2)}= \frac{\partial \displaystyle \Bigg |\sum _{l_1=0}^{M-1}\sum _{l_2=0}^{N-1} y(l_1,l_2)e^{-j\frac{2\pi k_1l_1}{M}} e^{-j\frac{2\pi k_2l_2}{N}}\Bigg |^2}{\partial y(n_1,n_2)}. \end{aligned}$$

The numerator is

$$\begin{aligned}&\Bigg |\sum _{l_1=0}^{M-1}\sum _{l_2=0}^{N-1} y(l_1,l_2)e^{-j\frac{2\pi k_1l_1}{M}} e^{-j\frac{2\pi k_2l_2}{N}}\Bigg |^2 \nonumber \\&=\Bigg [\sum _{l_1=0}^{M-1}\sum _{l_2=0}^{N-1} y(l_1,l_2)\Big (\underbrace{\cos \frac{2\pi k_1l_1}{M} \cos \frac{2\pi k_2l_2}{N} -\sin \frac{2\pi k_1l_1}{M} \sin \frac{2\pi k_2l_2}{N}}_{A(l_1,l_2)}\Big )\Bigg ]^2 \nonumber \\&\quad +\Bigg [\sum _{l_1=0}^{M-1}\sum _{l_2=0}^{N-1} y(l_1,l_2)\Big (\underbrace{\cos \frac{2\pi k_1l_1}{M} \sin \frac{2\pi k_2l_2}{N} +\sin \frac{2\pi k_1l_1}{M} \cos \frac{2\pi k_2l_2}{N}}_{B(l_1,l_2)}\Big )\Bigg ]^2.\nonumber \\ \end{aligned}$$

(29)

The derivatives of both the terms of (29) are

$$\begin{aligned}&\frac{\partial \displaystyle \bigg [\sum _{l_1=0}^{M-1}\sum _{l_2=0}^{N-1} y(l_1,l_2)A(l_1,l_2)\bigg ]^2}{\partial y(n_1,n_2)} \\&\quad =2\bigg ( \underbrace{\sum _{l_1=0}^{M-1}\sum _{l_2=0}^{N-1} y(l_1,l_2)A(l_1,l_2)}_{{\mathcal {R}}\{Y(k_1,k_2)\}}\bigg )\\&\qquad \times \Big ( \underbrace{\cos \frac{2\pi k_1n_1}{M} \cos \frac{2\pi k_2n_2}{N} -\sin \frac{2\pi k_1n_1}{M} \sin \frac{2\pi k_2n_2}{N}}_{A(n_1,n_2)}\Big ) \end{aligned}$$

and

$$\begin{aligned}&\frac{\partial \displaystyle \bigg [\sum _{l_1=0}^{M-1}\sum _{l_2=0}^{N-1} y(l_1,l_2)B(l_1,l_2)\bigg ]^2}{\partial y(n_1,n_2)} \\&\quad =2\bigg ( \underbrace{ \sum _{l_1=0}^{M-1}\sum _{l_2=0}^{N-1} y(l_1,l_2)B(l_1,l_2)}_{-{\mathcal {I}}\{Y(k_1,k_2)\}}\bigg )\\&\qquad \times \Big ( \underbrace{\cos \frac{2\pi k_1n_1}{M} \sin \frac{2\pi k_2n_2}{N} +\sin \frac{2\pi k_1n_1}{M} \cos \frac{2\pi k_2n_2}{N}}_{B(n_1,n_2)}\Big ). \end{aligned}$$

Hence, $\alpha _{n_1,n_2}$ is

$$\begin{aligned} \alpha _{n_1,n_2}= & {} \frac{2}{MN}\sum _{k_1=0}^{M-1}\sum _{k_2=0}^{N-1} \underbrace{ \frac{|H(k_1,k_2)|^2\lambda }{\Big (|H(k_1,k_2)|^2 +\frac{\lambda }{|Y(k_1,k_2)|^2}\Big )^2 \cdot |Y(k_1,k_2)|^4} }_{V(k_1,k_2): \text {real-valued}}\\&\cdot \Big ({\mathcal {R}}\{Y(k_1,k_2)\}A(n_1,n_2) -{\mathcal {I}}\{Y(k_1,k_2)\}B(n_1,n_2)\Big )\\= & {} \frac{2}{MN}\sum _{k_1=0}^{M-1}\sum _{k_2=0}^{N-1} {\mathcal {R}}\big \{\underbrace{V(k_1,k_2)Y(k_1,k_2)} _{G(k_1,k_2)}\big \}A(n_1,n_2)\\&-\frac{2}{MN}\sum _{k_1=0}^{M-1}\sum _{k_2=0}^{N-1} {\mathcal {I}}\big \{\underbrace{V(k_1,k_2)Y(k_1,k_2)} _{G(k_1,k_2)}\big \}B(n_1,n_2)\\= & {} 2{\mathcal {R}}\big \{g(n_1,n_2)\big \}, \end{aligned}$$

where $g(n_1,n_2)$ is the inverse Fourier transform of $G(k_1,k_2)=V(k_1,k_2)Y(k_1,k_2)$: $g(n_1,n_2)=\frac{1}{MN}\sum _{k_1=0}^{M-1}\sum _{k_2=0}^{N-1} G(k_1,k_2)e^{j\frac{2\pi n_1k_1}{M}} e^{j\frac{2\pi n_2k_2}{N}}$. Since $V(k_1,k_2)$ is real-valued, and $Y(k_1,k_2)$ is Hermitian symmetric, $g(n_1,n_2)$ is real-valued. Thus, we have $\alpha _{n_1,n_2}=2{\mathcal {R}}\{g(n_1,n_2)\} =2g(n_1,n_2)$.

Finally, the second term of (28) becomes

$$\begin{aligned} \alpha ^\mathrm {T}\mathbf {y}&=2(\mathbf {F}G)^\mathrm {T}\mathbf {F}Y =2G^\mathrm {T}\mathbf {F}^\mathrm {T}\mathbf {F}Y=2G^\mathrm {T}Y\\&=\sum _{k_1=0}^{M-1}\sum _{k_2=0}^{N-1} V(k_1,k_2)Y^*(k_1,k_2)Y(k_1,k_2)\\&=\sum _{k_1=0}^{M-1}\sum _{k_2=0}^{N-1} V(k_1,k_2)|Y(k_1,k_2)|^2\\&=\sum _{k_1=0}^{M-1}\sum _{k_2=0}^{N-1} \underbrace{ \frac{|H(k_1,k_2)|^2\lambda }{\Big (|H(k_1,k_2)|^2 +\frac{\lambda }{|Y(k_1,k_2)|^2}\Big )^2 \cdot |Y(k_1,k_2)|^2}}_{Q(k_1,k_2)}, \end{aligned}$$

where $\mathbf {F}$ is 2-D DFT matrix. The proof is completed. $\square$