Efficient algorithm to impose boundary conditions while applying a 2d filter

Question

It is possible to filter an image by using a 2d DFT. This will result in periodic boundary conditions. Now, is it possible to use the DFT to filter an image, while maintaining predefined boundary conditions (either Dirichlet or Neumann type).

I could implement a filter directly in the spatial domain. For example, solving the diffusion equation and pinnig pixel values, is equivalent to apply gaussian blur (strength depends on the number of iterations). However, this technique does not generalize very well, because one needs to find a PDE (and discretize it) whose Green's function has the desired behaviour in the frequency domain. Another approach would be to implement a convolution, and not touch the boundary itself but let it affect the surrounding pixels. However, this will becomes computationally expensive for larger images or filter sizes.

Answer 1

Another approach would be to implement a convolution, and not touch the boundary itself but let it affect the surrounding pixels. However, this will becomes computationally expensive for larger images or filter sizes.

You can do this by padding the image on all sides with pixels equal to your desired boundary. Make the padding one larger than the length of your desired filter (I'm assuming you're using a filter of finite extent). Then filter using the FFT method. When you're done, trim off the padding.

Edit:

Alternatively, based on Thomas's answer: there's a fast almost-Gaussian algorithm* out there that involves implementing a 2x2 or 3x3 "Gaussian" filter repeatedly. By the central limit theorem, the result tends nicely toward a real Gaussian, and because the kernel is teeny the math is fast.

It should be easy to either virtually or actually pad the edges of your image with one or two rows/columns of pixels with your boundary conditions. In fact, it should be much easier to satisfy Dirichlet, Neumann, or your choice of mixed boundary conditions with this approach.

* I should be honoring the creator by naming it, but I can't remember -- names fall out of my head as soon as I hear them.

Answer 2

First, make sure your problem is indeed faster in frequency domain. While asymptotically DFT based convolution should be faster in practice for small kernels it is not.

One way to do it is to pad (See MATLAB's padarray()) the array and then use valid equivalent convolution as in Replicate MATLAB's conv2() in Frequency Domain.

Answer 3

One such technique is the p+s or periodic plus smooth proposed by Lionel Moisan in Periodic plus smooth image decomposition, Journal of Mathematical Imaging and Vision, 2011, with code (perdecomp) and details. It basically uses an FFT.

This decomposition permits to avoid edge effects arising when the Discrete Fourier Transform of an image is computed, that are due to the implicit periodization. A discrete image u is decomposed under the form u = p + s, where p is called the periodic component of u and s is the smooth component of u. The periodic component p looks like the original image u, except that its Discrete Fourier Transform avoids the classical "cross structure" artifact. The smooth component s is an image that presents very slow variations inside the image domain.