# Estimate the Image Using Multi Realizations of Its Convolution with a Known Filters Using Wiener Filter

Suppose we have a corrupted image $$Y = H*X + \epsilon$$ formed by taking an image $$X$$, convolving it with a point-spread function $$H$$, and adding gaussian noise $$\epsilon$$. Then we know that the Wiener Filter can compute the MMSE estimator of $$X$$ given $$H$$ and the signal-to-noise ratio (SNR).

Given a set of $$n$$ images $$Y_i = H_i * X + \epsilon_i$$, is there a generalized Wiener Filter estimate of $$X$$ given the $$H_i$$'s and the SNRs?

(I'm cross-posting this from stats.SE, after someone suggested this community might be a better fit. Original link: https://stats.stackexchange.com/questions/540244/multi-image-wiener-filter)

• Yes, but I'm still working out the exact math. It's pretty obvious, but tedious -- given that, plus the fact that Wiener came out with his version in 1949, I suspect there's a paper on it published before 1966. Try searching on "multivariate Wiener filter", or "Wiener filter with multiple measurements". Aug 16, 2021 at 15:21
• @TimWescott, The math is pretty easy and elegant. You may see my answer.
– Royi
Aug 17, 2021 at 6:06
• Sorry for the delay, I just marked it! Aug 19, 2021 at 19:03

This is a nice question.

The math is actually pretty simple once you embrace the method I derived the Wiener Filter in - How Is the Formula for the Wiener Deconvolution Derived?

So, here is the model:

$$\boldsymbol{y}_{i} = \boldsymbol{h}_{i} \ast \boldsymbol{x} + \boldsymbol{w}, \; i = 1, 2, \ldots, n$$

Where $$\boldsymbol{w}$$ is an additive white gaussian noise which is independent of the signal and $$\boldsymbol{x} \sim N \left( 0, {\sigma}_{x}^{2} \right)$$.

Then the optimization model becomes (The MMSE Estimator):

$$\arg \min_{\boldsymbol{x}} \frac{1}{ 2 {\sigma}_{n}^{2} } \sum_{i = 1}^{n} {\left\| {H}_{i} \boldsymbol{x} - \boldsymbol{y}_{i} \right\|}_{2}^{2} + \frac{1}{2 {\sigma}_{x}^{2}} {\left\| \boldsymbol{x} \right\|}_{2}^{2} = \arg \min_{\boldsymbol{x}} \frac{1}{ { 2 \sigma}_{n}^{2} } {\left\| H \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} + \frac{1}{2 {\sigma}_{x}^{2}} {\left\| \boldsymbol{x} \right\|}_{2}^{2}$$

Where $$H = \begin{bmatrix} {H}_{1} \\ {H}_{2} \\ \vdots \\ {H}_{n} \end{bmatrix}$$ is the model matrix where $${H}_{i}$$ is the matrix form of the $$i$$ -th filter and $$\boldsymbol{y} = \begin{bmatrix} \boldsymbol{y}_{1} \\ \boldsymbol{y}_{2} \\ \vdots \\ \boldsymbol{y}_{n} \end{bmatrix}$$ is a concatenation of the output images in a column stack form.

Then the solution is given by:

$$\hat{\boldsymbol{x}} = {\left( {H}^{T} H + \frac{ {\sigma}_{n}^{2} }{ {\sigma}_{x}^{2} } I \right)}^{-1} {H}^{T} \boldsymbol{y}$$

• Oh that's nice. I need to brush up on my multivariate statistics. Aug 17, 2021 at 14:58
• Note you could also generalize it for measurement noise that differs by observation i. That just introduces weights to the least squares solution. Aug 17, 2021 at 17:21
• Also note that circulant matrices (circular convolutions) are diagonalized by the Fourier basis. Aug 17, 2021 at 17:22
• @MarkBorgerding, Mark, Weighing will be different not if the realizations of the noise are different (They in the model above) but if four each image we have different noise statistics (Different variance in this model). Also, I'm aware of the diagonalization of Circulant Matrices. This is the matrix form of the convolution theorem. Pay attention that usually in Image Processing $H$ isn't circulant as periodic boundary conditions introduce artifacts we don't want.
– Royi
Aug 17, 2021 at 17:37
• For Gaussian Distribution the MAP and the MMSE are the same (Since the maximum point is also the mean).
– Royi
Aug 19, 2021 at 21:26