Deriving the Inverse Filter of Image Convolution K

Question

Does anyone has an idea how to calculate the Inverse of a 2-D filter?

Let's say I have a 3x3 filter:
0 1 0
1 1 1
0 1 0
I want to find it's inverse. It's easy to do using DFT. But let's say I want to do it by convolution. Now, that's the problem, Matlab symbolic isn't my specialty. Assuming there's a 3X3 Inverse Filter it means convolution of the two will result in:
0 0 0
0 1 0
0 0 0
The problem is to create the right set of equations for that and solving it. Doing it with symbols is easy to think yet I couldn't do it.

Any ideas? Thanks.

P.S. I'm not sure there an Inverse Filter for this one as it has zeros in its DTFT.

Moreover, someone should allow Latex in this forum the way MathOverflow has.

Answer 1

In practice, one useful solution is Wiener deconvolution http://en.wikipedia.org/wiki/Wiener_deconvolution which basically, in Fourier space, divides by the spectrum of the given filter. Zeros are handled by adding a fudge term: instead of 1/H(w), use H(w)/( abs(H(w))^2 + c) where H(w) is the discrete Fourier transform of the filter h(x) (add a 2nd dimension as you like) and "w" is supposed to be omega. The constant is chosen based on the noise level of the signal or image.

Answer 2

Let h[n] be the finite impulse response of a 1D filter. What does that imply about its inverse filter? The inverse cannot possibly be FIR.

Proof: Let H(omega) G(omega) = 1 for all omega, where H is the DTFT of h[n], and G is the DTFT of g[n]. If h[n] is FIR, then g[n] must be IIR.

Of course, there are ways to approximate the inverse IIR filter with an FIR filter. A basic method is adaptive filtering, e.g., Least Mean Squares (LMS) algorithm. Or just truncate the IIR filter. You still need to worry about stability, though.

For practical purposes, there probably is no desirable solution to your specific question. Especially if, for example, this is in image processing and you are trying to inverse an FIR blur filter with FIR sharpening filter. The final image will not look that good, period, unless your sharpening filter is really, really large.

EDIT: Let y[n] = b0 x[n-0] + b1 x[n-1] + ... + bN x[n-N]. Let this equation characterize the forward system, where y is the output and x is the input. By definition, the impulse response is the output when the input is an impulse: h[n] = b0 d[n-0] + b1 d[n-1] + ... + bN d[n-N]. This impulse response has finite length N+1.

Now, consider the inverse system where x is the output and y is the input. Then the impulse response is described by the recurrence equation d[n] = b0 h[n] + b1 h[n-1] + ... + bN h[n-N]. Equivalently, b0 h[n] = d[n] - b1 h[n-1] - ... - bN h[n-N].

Without loss of generality, assume that b0 and bN are both nonzero. For any m, if h[m] is nonzero, then h[m+N] is also nonzero. Because this system has feedback, its impulse response is infinitely long. QED.

Causality does not matter. The inverse of a delay is an advance, and vice versa. Neither a delay or an advance alter the finiteness of an impulse response. Shift an infinite impulse response left or right; it is still infinite.

EDIT 2: For clarification, this proof is not related to my original "proof". One was in frequency domain, the other in time domain.

Answer 3

One way is to write a matrix representation of convolution with your filter and then try to find some (regularized) inverse of this matrix.

For example the 1D [1,-1] filter is maybe the simplest discrete differential approximation that exists.

Its matrix will have two diagonals 1 on main diagonal and -1 on first off diagonal.

It's inverse will be integral filter. An IIR filter with infinite trail of 1s. Since our matrix has limited size this will in practice mean filling up the matrix with 1s everywhere on one side of the diagonal.

Answer 4

For a separable filter (i.e., one in which you can pull out a horizontal and vertical filter which can be applied in any order and give the same result as the composite 2-D filter), you could attempt to calculate the inverse of each on individually.

The inverse of a filter H(w) is G(w)=1/H(w), so one way to do it is to take the impulse response (the h[n] time-domain coeffs) and inverse-DFT them. It's not always easy to get an analytic expression for such a filter, so you could either compute it numerically (approximating to the desired precision) or do adaptive inverse filtering, as Steve suggested. See Widrow and Stearn's Adaptive Signal Processing for more info on this latter method.

Answer 5

Deriving the Inverse Kernel of a Given 2D Convolution Kernel

This is basically a generalization of the question - Deriving the Inverse Filter of Image Convolution Kernel.

Problem Formulation

Given a Convolution Kernel $ f \in \mathbb{R}^{m \times n} $ find its inverse kernel, $ g \in \mathbb{R}^{p \times q} $ such that $ f \ast g = h = \delta $.

Solution I

One could build the Matrix Form of the Convolution Operator.
The Matrix form can replicate various mode of Image Filtering (Exclusive):

• Boundary Conditions (Applied as padding the image and apply Convolution in Valid mode).
  • Zero Padding Padding the image with zeros.
  • Circular
    Using Circular / Periodic continuation of the image. Match Frequency Domain convolution.
  • Replicate
    Replicating the edge values (Nearest Neighbor).
    Usually creates the least artifacts in real world.
  • Symmetric
    Mirroring the image along the edges.
• Convolution Shape
  • Full The output size is to full extent - $ \left( m + p - 1 \right) \times \left( n + q - 1 \right) $.
  • Same
    Output size is the same as the input size ("Image").
  • Valid
    Output size is the size of full overlap of the image and kernel.

When using Image Filtering the output size matches the input size hence the Matrix Form is square and the inverse is defined.
Using Convolution Matrix the matrix might not be square (Unless "Same" is chosen) hence the Pseudo Inverse should be derived.

One should notice that while the input matrix should be sparse the inverse matrix isn't.
Moreover while the Convolution Matrix will have special form (Toeplitz neglecting the Boundary Conditions) the inverse won't be.
Hence this solution is more accurate than the next one yet it also use higher degree of freedom (The solution isn't necessarily in Toeplitz form).

Solution II

In this solution the inverse is derived using minimization of the following cost function:

$$ \arg \min_{g} \frac{1}{2} {\left\| f \ast g - h \right\|}_{2}^{2} $$

The derivative is given by:

$$ \frac{\partial \frac{1}{2} {\left\| f \ast g - h \right\|}_{2}^{2} }{\partial g} = f \star \left( f \ast g - h \right) $$

Where $ \star $ is the Correlation operation.

In practice the convolution in the Objective Function is done in full mode (MATLAB idiom).
In practice it is done:

hObjFun     = @(mG) 0.5 * sum((conv2(mF, mG, 'full') - mH) .^ 2, 'all');
mObjFunGrad = conv2(conv2(mF, mG, 'full') - mH, mF(end:-1:1, end:-1:1), 'valid');

Where the code flips mF for Correaltion and use valid (In matrix form it is the Adjoint / Transpose -> the output is smaller).

The optimization problem is Strictly Convex and can be solved easily using Gradient Descent:

for ii = 1:numIteraions
    mObjFunGrad = conv2(conv2(mF, mG, 'full') - mH, mF(end:-1:1, end:-1:1), 'valid');
    mG          = mG - (stepSize * mObjFunGrad);
end

Full code is available in my StackOverflow Q2080835 GitHub Repository.