3

I want to deconvolve an image $h$ by a kernel $f$. More precisely, let

$$G = \operatorname*{argmin}_g \|f \ast g - h\|_2$$

be the set of least-squares solutions. I want to find the least-norm solution in this set:

$$ \hat{g} \in \operatorname*{argmin}_{g \in G} \|g\|_2 $$

which is

$$ \hat{g} = (f \mathrel{\ast})^+ h $$

where $^+$ is the pseudoinverse and $f \mathrel{\ast}$ is the Toeplitz matrix constructed from $f$.

What algorithm can efficiently compute $\hat{g}$, particularly in the setting of 2D images? Richardson–Lucy deconvolution doesn't work because I get division-by-zero errors (e.g. when $f$ is the discrete Laplace operator). Are there other, possibly gradient-based, iterative methods that don't suffer from this problem?

Royi
  • 10,050
user76284
  • 6,408
  • The eigenvalues of the convolution matrix form a sequence that converges to zero, at least in the limit as your image gets more and more pixels. Therefore, you need some way to deal with eigenvalues that are small, but nonzero. The pseudoinverse only deals with the eigenvalues that are exactly zero. A conventional approach which performs well is to regularize, and choose the regularization parameter via the L-curve method. – Nick Alger Nov 04 '24 at 21:22

1 Answers1

2

The problem can be written in terms of Linear Operator, a Matrix, and vectors. Where the vectors are the column stack form of the image.

The column stack of an image is the result of applying the Vectorization Operator on the image matrix. The problem becomes:

$$ \arg \min_{\boldsymbol{x}} \frac{1}{2} {\left\| \boldsymbol{H} \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} $$

Where the blurred image in a column stack form is given by $\boldsymbol{y}$, the matrix form of the blurring operator is given by $\boldsymbol{H}$.

There are 2 main challenges to the problem, even if the blurring operator is known:

  1. The Matrix Condition Number
    The blurring operator is basically an LPF filter.
    As such, there are frequencies which it might nearly zero out.
    Those will affect the condition number of the matrix operator (See the effect in 1D Deconvolution with Gaussian Kernel (MATLAB) which shows the connection).
    The solution is regularization. Be it in the form of Lucy Richardson, Tikhonov Regularization or a more effective model like the Total Variation Model (See How to Solve ADMM for TV Minimization Problem For Different Sizes $A$ and $x$ in $Ax=b$).
  2. The Scale of the Problem
    The dimensions of the matrix operator, $\boldsymbol{H} \in \mathbb{R}^{{n}^{2} \times {n}^{2}}$ where $n$ is the number of pixels in the image.
    Although the matrix is sparse, it makes it hard to use direct solvers.
    The correct way to approach the problem is define the operator using a wrapper within an iterative solver (Conjugate Gradient, Gradient Descent, etc...).
Royi
  • 10,050