A basic problem in optics concerns lenses that map a source into an image. If we imagine the source image as a collection of point sources of light, a "perfect" lens would map each unique point in the source plane onto a unique point in the image plane. Ignoring diffraction effects, this is precisely what a pinhole camera does.

Real lenses are not perfect. Light from a point source is spread out over a broad area in the image plane. This behavior is characterized by the so-called point-spread function, of PSF, of the lens. And when the PSF is known, it can, in principle, be inverted: the original image can be reconstructed from the "blurry" observed image.

Let us characterize the source plane by $I({\bf x})$ representing the intensity of light at the location ${\bf x}$. Similarly, let us characterize the image plane by $I'({\bf x}')$. Let the PSF be ${\tt PSF}({\bf x},{\bf x}')$. The relationship between source and image is, then, given by

$$I'({\bf x}')=\iint d^2{\rm x}~I({\bf x}){\tt PSF}({\bf x},{\bf x}').\tag{1}$$

This is obviously the definition of an integral operator, the convolution operator:

$$I'({\bf x}')={\tt PSF}({\bf x},{\bf x}')\otimes I({\bf x}).\tag{2}$$

Inverting this operator is a formidable task. While it is possible to represent the operator in matrix form, the resulting matrix is very large, and inverting it is not only computationally expensive but also fraught with numerical instabilities. Hence, a better method is needed. Fortunately, such a method is readily available, so long as the PSF satisfies a simple condition: namely, that it can be written in the form

$${\rm PSF}({\bf x},{\bf x}')={\tt PSF}(\beta{\bf x}+{\bf x}'),\tag{3}$$

with $\beta$ constant. In this case, we can Fourier-transform Eq. (1) as follows:

$$\begin{align}
\hat{I}'({\bf f}')=\iint d^2{\bf x}'e^{-2\pi i{\bf x}'\cdot{\bf f}'}I'({\bf x}')
&=\iint d^2{\bf x}' e^{-2\pi i{\bf x}'\cdot{\bf f}'}\iint d^2{\bf x}I({\bf x}){\tt PSF}(\beta{\bf x}+{\bf x}')\nonumber\\
&=\iint d^2{\bf x}e^{-2\pi i(-\beta{\bf x})\cdot{\bf f}'}I({\bf x})
\iint d^2{\bf x}'e^{-2\pi i(\beta{\bf x}+{\bf x}')\cdot{\bf f}'}{\tt PSF}({\beta\bf x}+{\bf x}')\nonumber\\
&=\widehat{\tt PSF}({\bf f}')\iint d^2{\bf x}e^{-2\pi i(-\beta{\bf x})\cdot{\bf f}'}I({\bf x})
=\widehat{\tt PSF}({\bf f}')\hat{I}(-\beta{\bf f}').\tag{4}
\end{align}$$

Consequently, if we know the Fourier-transform of the observed signal in the frequency domain, we can recover the original through simple division:

$$\hat{I}(-\beta{\bf f}')=\frac{\hat{I}'({\bf f}')}{\widehat{\tt PSF}({\bf f}')}.\tag{5}$$

This result is a specific application of the so-called convolution theorem.