Fourier-Domain Optimization for Image Processing

by   Majed El Helou, et al.

Image optimization problems encompass many applications such as spectral fusion, deblurring, deconvolution, dehazing, matting, reflection removal and image interpolation, among others. With current image sizes in the order of megabytes, it is extremely expensive to run conventional algorithms such as gradient descent, making them unfavorable especially when closed-form solutions can be derived and computed efficiently. This paper explains in detail the framework for solving convex image optimization and deconvolution in the Fourier domain. We begin by explaining the mathematical background and motivating why the presented setups can be transformed and solved very efficiently in the Fourier domain. We also show how to practically use these solutions, by providing the corresponding implementations. The explanations are aimed at a broad audience with minimal knowledge of convolution and image optimization. The eager reader can jump to Section 3 for a footprint of how to solve and implement a sample optimization function, and Section 5 for the more complex cases.



page 1

page 2

page 3

page 4


Algorithmic Optimisations for Iterative Deconvolution Methods

We investigate possibilities to speed up iterative algorithms for non-bl...

Phase and TV Based Convex Sets for Blind Deconvolution of Microscopic Images

In this article, two closed and convex sets for blind deconvolution prob...

Tuning Multigrid Methods with Robust Optimization

Local Fourier analysis is a useful tool for predicting and analyzing the...

Deep Learning for space-variant deconvolution in galaxy surveys

Deconvolution of large survey images with millions of galaxies requires ...

Deconvolving Images with Unknown Boundaries Using the Alternating Direction Method of Multipliers

The alternating direction method of multipliers (ADMM) has recently spar...

Fourier Image Transformer

Transformer architectures show spectacular performance on NLP tasks and ...

High Speed Tracking With A Fourier Domain Kernelized Correlation Filter

It is challenging to design a high speed tracking approach using l1-norm...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

A wide range of image processing applications are tackled in the literature by minimizing a carefully-defined loss function. These application-specific loss functions generally have one thing in common: the optimization variable is of the same order of magnitude in size as the image, and often it is the output image itself. This poses a big computational challenge when considering large images. Given that image sizes are consistently increasing, such optimization-based solutions are not suitable for deployment on low energy and mobile devices.

To remedy the computational problem, there exist different approaches that attempt to avoid explicitly formulating an optimization problem and settle instead for sub-optimal heuristics. This is because the current optimization schemes, which are essentially gradient descent based, are unable to converge fast enough with megapixel images. On the other hand, a high resolution is necessary for a multitude of computer vision applications as well as for better user experience. So an ideal solution would be capable of solving the formulated image optimization problem even for large image sizes. The bottom line is we desire high quality optimization without sacrificing computational efficiency or user experience.

What we propose in this paper is not a novel method, rather a thorough explanation, visualization, and proof of validity of a very powerful optimization solver based on the use of the Fourier transform. In fact, the technique we explain has been central to several recent works of research 

deconvHyperLap ; IRCNN ; l0smooth ; antiGlare . However, a thorough explanation of the method, including its mathematical background and implementation details, is lacking in the literature. Using a simple non-blind deconvolution problem as an example, we outline the details of the method, and elaborate upon the implementation details, in Section 2. After the theoretical and practical grounds are established, we present and explain all the steps involved in solving the image optimization problem in Section 3. We also provide publicly accessible code to complement the paper.

The presented method addresses optimization problems that have a closed-form solution. In cases where there are inter-pixel relations that cannot be translated to convolutions, for example resulting from pixel-position-dependent mappings or complex cost functions, the direct method does not apply. For such situations we additionally address "half-quadratic splitting" techniques, which allow to split (as the name suggests) the optimization into two steps, one of which can typically be solved in closed form, and the other with a traditional gradient descent approach. The final solution for such problems is obtained by iterating over these two steps nick ; deconvHyperLap ; kruse2017learning ; IRCNN . More details are provided in Section 5.

2 Simple Non-blind Deconvolution

In what follows, we assume that the reader has some background in basic linear algebra and image convolution-deconvolution. We assume that the number of pixels in an image is . We treat one-channel images, and assume the pixel values to lie in .

2.1 Mathematical Framework

In the simple case of deconvolution, we observe an image , which is the result of the convolution between an image and a convolution kernel . The size of the observed image is reduced with respect to the input image such that and

. This size reduction is commonly avoided by zero-padding. The goal is to recover an accurate approximation of the input image

. The kernel can be the point spread function (PSF) of an acquisition system, omnipresent in coded aperture systems for instance. In non-blind deconvolution, we assume to be known.

We start with an example convolution between a kernel and a matrix (our image),


where stands for 2D convolution. This convolution can be written in terms of a matrix multiplication. For this, let us define the flipped and zero-filled kernel which in our example is given by:


where the mirroring is needed to mimic the convolution operation (mirroring and sliding), and the zero-padding is necessary to fit the image dimensions, as we will see shortly. The length of corresponds to the number of pixels in our image (

). We also need to reshape the matrix into a column vector in row-wise order, yielding:


We define the mapping which takes an input vector of length and maps it to a matrix of size such that the matrix is circulant, and made up of integer shifts of the input vector. The operator applied to a 3 dimensional vector yields, for instance:


This allows us to rewrite the convolution equation in terms of a matrix multiplication:


where is the vectorized result of the nonvalid convolution and is, in this particular example:


The appellation nonvalid is due to boundary effects: unlike which is , is , with the extra column appended as an artifact of the convolution (dependent on the kernel size). For instance, in our original example, the result is:


Notice, however, that is exactly equal to the first columns of . This is what is referred to as the "circular boundary conditions" when using this framework in the literature deconvHyperLap .

We could solve Eq. (5) for by inverting the circulant matrix. Note however that the circulant matrix is of size , and such an inversion is very costly. Therefore, we move to the Fourier domain, exploiting the following well-known theorem.

Theorem 1 The discrete Fourier transform (DFT) matrix diagonalizes any circulant matrix. In other words, any matrix in the frequency domain, which is the transform of a circulant matrix, is a diagonal matrix.

Based on Theorem 1, in Eq. (5) becomes a diagonal matrix in the Fourier domain (denoted by ). The multiplication in the time/space domain naturally becomes a convolution in the Fourier domain:


Two important things to note here are that, first, we are consistently utilizing the DFT across Eq. (8). And second, the convolution is between a diagonal matrix and a column vector, making it effectively nothing more than a point-wise multiplication of their entries. This point-wise multiplication is seen in Eq. (12) but where the non-trivial entries are compactly collected into two small matrices instead of one large diagonal matrix and a corresponding-length vector. This means that every entry in is multiplied by one diagonal element in to obtain the corresponding element in . Therefore, a simple point-wise division is enough to recover , and applying the inverse Fourier transform yields , that is, the deconvolved image in vectorized form.

2.2 Kernels with Two Dimensions

The generalization to 2D convolution kernels is straight-forward. The only difference relative to the previous section is in the construction of . Taking as an example a 2D kernel :


the corresponding , for the same image size of , is given by:


In general, always contains the unrolled version of a doubly-mirrored or flipped version of called :


The way is unrolled is row by row: the first row is filled into the first entries of , then the second row is filled in but beginning at entry , where is the number of columns in the image (i.e. the size of an image row). In this example, is filled in at position . This procedure is repeated until all the rows of have been traversed and filled in . The remainder of is left as the original zero filling.

For completeness, it is theoretically possible, although practically not relevant, to have a convolution kernel wider than the image and thus not fitting in this construct. In that case, however, there is not a single valid convolution being carried out, since the size of the valid output normally is where and . Thus the valid convolution is the empty set whenever or .

2.3 Practical Implementation

In the previous section we have leveraged two properties of the Fourier transform: a time/space-domain convolution corresponds to a multiplication in the Fourier domain, and the DFT diagonalizes circulant matrices. However, an obvious issue in Eq. (5) is the data complexity. For an image containing (possibly millions of) pixels, the matrix is of size , making it very expensive to take its Fourier transform or even to simply store it. However, once in the Fourier domain, only diagonal entries are non-zero, so the matrix is sparse. In both MATLAB and python there exist libraries that can handle relatively large sparse matrices and can compute their inverses, however, these are (as of this date) limited to too small image sizes. In practical implementations, it is possible to bypass these steps and to go directly to the Fourier domain, which is what we explain in this section.

So instead of converting our images into vectorized forms, we can keep them in their original matrix form, and compute and arrange the diagonal entries of into a single matrix of same size as the images. It is also possible to work with the two-dimensional Fourier transform of the images, if the corresponding entries in the diagonal matrix ( Eq. (8)) are adjusted. And by adjusted, we mean that a one-dimensional Fourier transform needs to be applied to the matrix-rearranged diagonal entries. This is exactly what does, by first distributing the entries of the input kernel over a matrix of same size as the image, then applying the two-dimensional Fourier transform to that resulting matrix. In the provided implementations, the convolution operation corresponds to:


where and are built-in functions, is point-wise multiplication, holds the dimensions of the image, and is the mirrored version of . All what does is compute the diagonal entries of the diagonal matrix by reshaping into a matrix and taking the DFT of that matrix. The deconvolution simply becomes:


where the division is again point-wise. The final deconvolved image can be obtained by going back to the spatial domain using the inverse Fourier transform ().

(a) Visualization of the operation of . In colors are from top to bottom: the convolution kernel, the same kernel with entries rearranged, the two-dimensional Fourier transform of the latter matrix. In grey are similarly the input image, and its two-dimensional Fourier transform. (b) Visualization of the large convolution in the Fourier domain between the Fourier transform of the circulant matrix and the Fourier transform of the vectorized image. This effectively amounts to a point-wise multiplication.
Figure 1: Visualizations of python implementations of both (a) the approach and (b) the large matrix multiplication on vectorized images.

2.4 Remarks

1) Note that the one-dimensional kernel we used for illustration implements an image gradient in the horizontal direction (vertical edge detector). This will be used in the next section.

2) Circulant matrices have some well-defined properties, which are exploited in the proposed solution in Section 3. These properties are outlined in the theorems below.

Theorem 2 The transpose of a circulant matrix is also a circulant matrix.
Proof The proof is trivial and this can easily be visualized. Let be the bottom row of , and its mirrored version, we then have , and thus is circulant.

Theorem 3 Multiplying a circulant matrix by a circulant matrix yields a circulant matrix .
Proof We know from Theorem 2 that is also circulant; . (a) This implies that not only the rows but also the columns of are circulant variants of a single vector. Let denote the circular shifting of the vector by steps. (b) The inner product only depends on the difference . (a) & (b) imply by construction that where and is the number of columns of .

Theorem 4 The sum of a circulant matrix and a circulant matrix is also a circulant matrix. And is a linear function.
Proof The sum is trivially the circulant matrix . And also trivially for any constant .

3 Convex Optimization Solution

3.1 Mathematical Framework

In this section, we explain how closed-form solutions to convex image optimization problems can be computed efficiently (), instead of running a descent algorithm. Such approaches have been used in deconvHyperLap , l0smooth , and antiGlare , however, a clear explanation and implementation details are missing in the literature.

A common denominator of many image processing problems is that a solution can be found by minimizing a cost function, tailored to the problem, over a high-dimensional argument, typically the image itself. One example of such cost functions is:


where the first and second cost terms enforce data consistency between the given guiding matrices and , and the image on one side and the image gradient on the other side, respectively, weighted by the regularization factor . This cost function can be augmented with any number of cost terms including additional kernels and guiding matrices. Remember that the gradient is also a convolution with the kernel discussed in Section 2. Note that we cannot allow point-wise operations on the optimization argument if we want a closed-form solution, which is discussed in the following section, to exist. Otherwise, half-quadratic splitting is required for the optimization solution to benefit from the proposed approach.

The loss function in Eq. (14) is a summation of squared Frobenius norms and is thus convex with a unique solution that can be found in closed form. The optimal solution can be found at the zero-gradient location, so we differentiate the loss function with respect to our optimization variable image . To do this, we write the loss in terms of matrix multiplications (with no convolution operations) and vectorized images. Based on our analysis in the previous section, this leads to the following, equivalent loss:


where and are defined as in Section 2. Note that they do not have to be computed in practice, as we can directly use their Fourier representation, as explained in Section 2.3. Taking the partial derivative of our loss function, we get:


Expanding and setting the derivative to zero yields:


which can be rearranged in the standard form that is analyzed in Section 2:


where is a circulant matrix based on Theorems 2, 3 and 4, and is a known vector.

3.2 Practical Implementation

We can solve Eq. (18) in the Fourier domain as derived in Section 2.3. The Fourier domain equivalent of is given by:


where denotes the conjugate or Hermitian transpose of a matrix. As before, we obtain a non-sparse representation of the Fourier coefficients by leveraging :


where is the size of the (non-vectorized) input image. To sum up with, the optimal solution becomes:


where the division and all multiplications are point-wise, and is the MATLAB symbol for the conjugate or Hermitian transpose.

4 GitHub demos

On the GitHub page are made available a MATLAB and a python demo that can be downloaded and run without any additional setup or toolboxes needed. The demos are split into two parts A and B as follows.

Part A illustrates the simple convolution operation described in Section 2. The convolution is carried with the built-in functions, but also with the large matrix multiplication as well as in the Fourier domain. One can visualize the results, which are automatically printed, and notice that the nonvalid entries (last column in the chosen example) are different, but all valid entries are consistent.

Part B illustrates how to solve a guided deblurring image optimization in the Fourier domain. The script reads a guide RGB image, and a blurry near-infrared (NIR) image. The deblurred output is the image minimizing the loss function that is given by:


where is a regularization factor (set to 1), is the blurry NIR input, and

the estimated blur. For simplicity, it is assumed to be Gaussian and constant across the entire image.

is the luminance of the RGB guide and is the optimization variable. The interested reader can refer to the following works for the more advanced and better performing versions of this multispectral deblurring loss function ICIP ; zahra .

5 Half-Quadratic Splitting

The previous sections explain how to write-out the closed-form solution of the optimization, and how to solve for it efficiently in the Fourier domain. However, as we discuss in our introduction, not any optimization can be solved directly with this approach, without half-quadratic splitting. Let us consider a general loss function we aim to minimize:


where is non-quadratic in and a solution cannot be written-out in closed form. This is where half-quadratic splitting comes into play. The loss can be minimized iteratively by optimizing for and one after the other. For that, the optimization loss is rewritten in terms of and , where is a latent variable used to split the two functions in for the iterative optimization. We obtain the new loss, which we minimize with respect to both and :


With this new definition, can be minimized by minimizing over , which can be done efficiently using the approach we described in previous sections, since is independent of , then minimizing over , for instance with descent approaches. However, this method fails to converge to the same solution as , since and are not constrained to being equal. To enforce this and converge to the minimizer of , we add another term to the loss:


where is a small positive constant that is increased with every minimization iteration. As , we obtain (as a result of minimizing , since otherwise their difference makes the loss infinitely large) thus making and equivalent. The solution is obtained by alternately minimizing over and with the following steps (shown for time step ):


where is an increasing function, and can be initialized randomly or with an educated guess when possible. Eq. (26) can be solved with the method presented in this paper, and Eq. (27) needs to be solved with descent algorithms.

The half-quadratic splitting approach we present belongs to the family of variable splitting optimization techniques, with the particularity that the splitting approach is designed specifically to leverage the closed-form Fourier solution. More precisely, the variable splitting optimization is carried with the penalty approach, where the penalty term added in is responsible for constraining to be equal to . This constraint is made less and less loose as we increase the value of with every iteration. This behavior can be controlled by tuning the parameter .

Another variable splitting strategy is similar to optimizing under the constraint . The constrained optimization can be carried out with the augmented Lagrangian method, or ADMM (first appearing in ADMM ) when the constraint has more than one component vector and it is possible to alternate the optimization. For the purpose of this paper, we only present the half-quadratic splitting technique because its formulation allows the leveraging of the closed-form Fourier-domain solution that we present.

6 Final Remarks

  1. It is preferable, when a closed-form solution is readily available, to make use of it and avoid gradient descent approaches. The latter are far less efficient, and also require hyper-parameter tweaking. Obtaining the closed-form solution in the Fourier domain is itself more efficient both computationally ( improvement) and memory-wise ( instead of ) compared to matrix inversion.

  2. When a closed-form solution cannot be formulated, half-quadratic splitting can be used to divide the optimization into two steps, one having a closed-form solution and the other requiring traditional descent-based optimization.

  3. In some (simplistic) scenarios where the optimization has a "null-space", the optimal solution is non-unique. This is simpler to reason about with an example. One such case is when the optimization variable is always convolved in the loss function with a gradient kernel (ex: [-1 1]) that is symmetric in absolute value but of opposite signs on each side of the center. Such an operator is not affected by image-wide constant shifts, since they get added then removed by the convolution with [-1 1]. This means that any constant shift added to the optimal solution does not affect the loss value, and will itself also be an optimal solution too. This issue is a well-known problem in image deconvolution.

7 Conclusion

In this paper, we present the framework for efficiently solving image deconvolution and, more generally, a variety of image optimization problems, in the Fourier domain. We explain the mathematical background in simple terms, as well as the steps involved in the derivation of the optimization solution. We provide an open-source MATLAB and python demo showing the presented theory in practice for two selected examples. Researchers can leverage this framework for solving other image-related problems, by formulating the matching optimization and solving it efficiently for very large image sizes. We also explain briefly the approach of half-quadratic splitting, which makes it possible to leverage the proposed method in the case of more general loss functions.

We invite the readers who implement this approach in other programming languages to contribute a sample implementation or demo similar to ours to the GitHub repository.


The authors would like to thank Dr. Zahra Sadeghipoor and Dr. Nikolaos Arvanitopoulos for valuable discussions and advice.


  • (1) N. Arvanitopoulos, R. Achanta, and S. Susstrunk. Single image reflection suppression. In

    Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)

    , pages 4498–4506, 2017.
  • (2) M. El Helou, Z. Sadeghipoor, and S. Süsstrunk. Correlation-based deblurring leveraging multispectral chromatic aberration in color and near-infrared joint acquisition. In IEEE International Conference on Image Processing (ICIP), 2017.
  • (3) R. Glowinski and A. Marroco. Sur l’approximation, par éléments finis d’ordre un, et la résolution, par pénalisation-dualité d’une classe de problèmes de dirichlet non linéaires. Revue française d’automatique, informatique, recherche opérationnelle. Analyse numérique, 9(R2):41–76, 1975.
  • (4) D. Krishnan and R. Fergus. Fast image deconvolution using hyper-laplacian priors. In Advances in Neural Information Processing Systems (NIPS), pages 1033–1041, 2009.
  • (5) J. Kruse, C. Rother, and U. Schmidt. Learning to push the limits of efficient fft-based image deconvolution. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 4586–4594, 2017.
  • (6) Y. Li and M. S. Brown. Single image layer separation using relative smoothness. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 2752–2759, 2014.
  • (7) Z. Sadeghipoor Kermani, Y. Lu, E. Mendez, and S. Süsstrunk. Multiscale guided deblurring: Chromatic aberration correction in color and near-infrared imaging. In European Conference on Signal Processing (Eusipco), 2015.
  • (8) T. Sandhan and J. Y. Choi. Anti-glare: Tightly constrained optimization for eyeglass reflection removal. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 1241–1250, 2017.
  • (9) K. Zhang, W. Zuo, S. Gu, and L. Zhang. Learning deep cnn denoiser prior for image restoration. In IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 2808–2817, 2017.