pix/

# Scattering and Gathering for Spatially Varying Blurs

Nicholas Chimitt, , Xingguang Zhang, ,
Yiheng Chi, , and Stanley H. Chan The authors are with the School of
Electrical and Computer Engineering, Purdue University, West Lafayette, IN
47907, USA. Corresponding author: Stanley Chan, email<EMAIL_ADDRESS>work is supported in part by the Intelligence Advanced Research Projects
Activity (IARPA) under Contract No. 2022-21102100004, and in part by the
National Science Foundation under the grants CCSS-2030570 and IIS-2133032. The
views and conclusions contained herein are those of the authors and should not
be interpreted as necessarily representing the official policies, either
expressed or implied, of IARPA, or the U.S. Government. The U.S. Government is
authorized to reproduce and distribute reprints for governmental purposes
notwithstanding any copyright annotation therein.

###### Abstract

A spatially varying blur kernel $h(\mathbf{x},\mathbf{u})$ is specified by an
input coordinate $\mathbf{u}\in\mathbb{R}^{2}$ and an output coordinate
$\mathbf{x}\in\mathbb{R}^{2}$. For computational efficiency, we sometimes
write $h(\mathbf{x},\mathbf{u})$ as a linear combination of spatially
invariant basis functions. The associated pixelwise coefficients, however, can
be indexed by either the input coordinate or the output coordinate. While
appearing subtle, the two indexing schemes will lead to two different forms of
convolutions known as _scattering_ and _gathering_ , respectively. We discuss
the origin of the operations. We discuss conditions under which the two
operations are identical. We show that scattering is more suitable for
simulating how light propagates and gathering is more suitable for image
filtering such as denoising.

###### Index Terms:

Spatially varying blur, basis representation, scattering, gathering

##  I Introduction

In the two-dimensional space, the convolution between an input image
$J(\mathbf{x})$ and a shift-invariant kernel $h(\mathbf{x})$ produces an
output image $I(\mathbf{x})$ via the well-known integral

$I(\mathbf{x})=\int_{-\infty}^{\infty}h(\mathbf{x}-\mathbf{u})J(\mathbf{u})\;d\mathbf{u}.$
(1)

In this equation, $\mathbf{x}\in\mathbb{R}^{2}$ is a two-dimensional
coordinate in the output space and $\mathbf{u}\in\mathbb{R}^{2}$ is a
coordinate in the input space. This definition is ubiquitous in all shift-
invariant systems.

If the kernel $h$ is spatially _varying_ , then it is no longer a function of
the coordinate difference $\mathbf{x}-\mathbf{u}$ but a function of two
variables $\mathbf{x}$ and $\mathbf{u}$. The resulting kernel
$h(\mathbf{x},\mathbf{u})$ will give the input-output relationship via the
integral

$I(\mathbf{x})=\int_{-\infty}^{\infty}h(\mathbf{x},\mathbf{u})J(\mathbf{u})\;d\mathbf{u},$
(2)

also known as the superposition integral. That is, at every output coordinate
$\mathbf{x}$, there is a kernel $h(\mathbf{x},\mathbf{u})$ which is a function
of the input coordinate $\mathbf{u}$.

While spatially varying kernels are more difficult to analyze because they
cannot be directly handled by Fourier transforms [1], they are common in image
_formation_ and image _processing_ which are twins in many situations, such as
in the case of kernel estimation [2, 3, 4] or image restoration [5, 6, 7]. In
image formation, the spatially varying kernels are used to model how light
propagates from the object plane to the image plane. These kernels are known
as the point spread functions (PSF) which may be spatially varying due to
various degradations in the medium or aberrations in the imaging system (such
as spherical aberration). In image processing, the spatially varying kernels
are used to filter the input image for applications such as denoising or
interpolation. The spatially varying nature in these situations can come from
examples such as non-local edge-aware filters shapes and orientations of the
filters change depending on the image.

The theme of this paper is about the decomposition of the kernel in terms of
basis functions. If $h$ is _spatially invariant_ , we may express it via the
equation

$h(\mathbf{x}-\mathbf{u})=\sum_{m=1}^{M}a_{m}\varphi_{m}(\mathbf{x}-\mathbf{u}),$
(3)

where $\\{\varphi_{1},\varphi_{2},\ldots,\varphi_{M}\\}$ are orthogonal basis
functions. These functions could be as simples as the derivatives of the
Gaussians, or they can be learned from a dataset of kernels via the principal
component analysis. The scalars $\\{a_{1},a_{2},\ldots,a_{M}\\}$ are the basis
coefficients. They are often constructed according the local image statistics
or the underpinning physics.

The decomposition of $h$ into orthogonal basis functions can be useful from a
parameterization point of view. While $h$ is high-dimensional, the
decomposition allows us to represent $h$ in a low-dimensional space using the
set of $M$ coefficients. For applications such as blind deconvolution, this
low-dimensional representation can be effective for kernel estimation because
the search space is smaller.

In the case of _spatially varying_ kernels, the subject of this paper, the
basis representation in (3) needs to be modified so that it can take into
account of the two variables $\mathbf{x}$ and $\mathbf{u}$. However, there is
an ambiguity due to the existence of two options which we call _gathering_ and
_scattering_ :

$\displaystyle\text{(Gathering)}\qquad h(\mathbf{x},\mathbf{u})$
$\displaystyle=\sum_{m=1}^{M}a_{\mathbf{x},m}\varphi_{m}(\mathbf{x}-\mathbf{u}),$
(4) $\displaystyle\text{(Scattering)}\qquad h(\mathbf{x},\mathbf{u})$
$\displaystyle=\sum_{m=1}^{M}a_{\mathbf{u},m}\varphi_{m}(\mathbf{x}-\mathbf{u}).$
(5)

In both options, the spatially varying kernel $h(\mathbf{x},\mathbf{u})$ is
written as a combination of _invariant_ kernels
$\\{\varphi_{1},\varphi_{2},\ldots,\varphi_{M}\\}$. These $\varphi_{m}$’s are
spatially invariant, so they can be written as
$\varphi_{m}(\mathbf{x}-\mathbf{u})$. The difference between the two options
lies in the coefficient $a_{\mathbf{x},m}$ and $a_{\mathbf{u},m}$. In the
former case, the basis coefficient $a_{\mathbf{x},m}$ is indexed by the output
coordinate $\mathbf{x}$. For any pixel $\mathbf{x}$ in the output space, the
equation linearly combines the basis functions at that output pixel location
through $\\{a_{\mathbf{x},1},\ldots,a_{\mathbf{x},M}\\}$. In the latter case,
the basis coefficient $a_{\mathbf{u},m}$ is indexed by the input coordinate
$\mathbf{u}$. It is not immediately obvious why one would want to do so, but
now we do so via symmetry, opting to discuss this later in the paper.

Remark: Readers may wonder if we can define $h(\mathbf{x},\mathbf{u})$ using a
global $a_{m}$ instead of a pixelwise $a_{\mathbf{x},m}$ or
$a_{\mathbf{u},m}$. If we do so, for example by defining
$h(\mathbf{x},\mathbf{u})=\sum_{m=1}^{M}a_{m}\varphi_{m}(\mathbf{x}-\mathbf{u})$,
then $h(\mathbf{x},\mathbf{u})$ will be invariant because it is a linear
combination of invariant basis functions. This will defeat the purpose
studying a set of varying kernels. $\square$

At the first glance, the two choices above seem subtle to an extent that one
may expect a minor difference in terms image quality. However, the two
equations have two fundamentally different physical meanings. Even though the
resulting images may look similar, one of them is better suited for image
formation and the other is for image processing. To give readers a preview of
the main claims of the paper, we summarize them as follows:

  * •

Gathering: $a_{\mathbf{x},m}$ is for image _processing_ such as denoising
filter, interpolation filter, etc.

  * •

Scattering: $a_{\mathbf{u},m}$ is for image _formation_ such as modeling
atmospheric turbulence.

It should be noted that the study of time variant systems goes back to
classical time-variant filter banks, wavelets, or Kalman filtering/state-space
applications [8, 9, 10, 11, 12, 13]. The study of spatially varying kernels as
a sum of invariant ones has been performed in works such as [14, 15, 16]. In
these works, particular aspects of the problem are analyzed as an end-use
case. In this work, we focus on motivating the difference between the two
approximations from the side of modeling and describing where each one is more
appropriately applied.

##  II Computational Aspects of Gathering and Scattering

###  II-A Understanding Gathering

Gathering is a decomposition using the _output_ coordinates. Suppose that we
pass an image $J(\mathbf{u})$ through a spatially varying kernel
$h(\mathbf{x},\mathbf{u})$. Assuming that the spatially varying kernel
$h(\mathbf{x},\mathbf{u})$ has a basis representation shown in (4), by
substituting it into (2), we can show that

$\displaystyle I(\mathbf{x})$ $\displaystyle\overset{\text{by}\,\eqref{eq: I
varying}}{=}\int_{-\infty}^{\infty}h(\mathbf{x},\mathbf{u})J(\mathbf{u})\;d\mathbf{u}$
$\displaystyle\overset{\text{by}\,\eqref{eq:
x}}{=}\int_{-\infty}^{\infty}\left(\sum_{m=1}^{M}a_{\mathbf{x},m}\varphi_{m}(\mathbf{x}-\mathbf{u})\right)J(\mathbf{u})\;d\mathbf{u}$
$\displaystyle=\sum_{m=1}^{M}a_{\mathbf{x},m}\left(\int_{-\infty}^{\infty}\varphi_{m}(\mathbf{x}-\mathbf{u})J(\mathbf{u})\;d\mathbf{u}\right).$

Recognizing that the integral is a spatially invariant convolution, we can
show that

$\displaystyle I(\mathbf{x})$ $\displaystyle=\underset{\text{linear
combination of invariant
blurs}}{\underbrace{\sum_{m=1}^{M}a_{\mathbf{x},m}\;\;\underset{\text{invariant
blur}}{\underbrace{(\varphi_{m}\circledast J)(\mathbf{x})}}}}.$ (6)

Therefore, the spatially varying convolution is now replaced by a sum of $M$
spatially _invariant_ convolutions. Each term in the sum
$(\varphi_{m}\circledast J)(\mathbf{x})$ is the convolution between the image
$J(\mathbf{u})$ and the basis functions $\varphi_{m}(\mathbf{u})$.

In terms of computation, if we start with the image $J(\mathbf{u})$, we can
pre-compute the filtered images $(\varphi_{m}\circledast J)(\mathbf{x})$, and
form a weighted sum of these images to construct the output $I(\mathbf{x})$.
The weights are pixelwise, and so the resulting blur is spatially varying.
Figure 1 illustrates how this idea can be implemented.

[width=]./pix/Figure01_Gathering_Block_Diagram

Figure 1: Gathering is for Image Processing. In gathering, we construct a set
of convolved images $(\varphi_{m}\circledast J)(\mathbf{x})$ and _gather_ them
via pixelwise multiplication. The process offers a significant saving in terms
of computation.

The illustration in Figure 1 should remind readers of many classical signal
processing techniques. For example, the steerable filters by Freeman and
Adelson [17] employed the idea to detect image edges using a bank of Laplacian
filters. The idea was also used in applied mathematics literature. For
example, Nagy and O’Leary [18] used the idea to decompose the spatially
varying kernels using non-overlapped blocks. Nagy and O’Leary’s formulation is
a special case of our model by making the coefficients $a_{\mathbf{x},m}$ as a
binary mask. In [19], Popkin et al. used the idea to perform fast computation
of spatially varying blurs with the basis functions computed via principal
component analysis.

To readers who prefer the perspective of deep learning, the gathering given in
(4) is what we called _convolution_ in deep neural networks. (This should be
contrasted with scattering which will be discussed shortly. Scattering is
equivalent to _transposed convolution_.) To see why gathering is identical to
the convolution (in deep-learning), we can refer to Figure 2. Starting with a
local neighborhood in the input space, we assign each pixel a weight which is
specified by the kernel and map the sum to the output coordinate.

[width=0.65]./pix/Figure02_Gathering_Illustration

Figure 2: The operation derived in (6) follows the intuitive argument that a
neighboring input pixels are mapped to an output. Thus, each
$h(\mathbf{x},\mathbf{u})$ is defined based on the output coordinate
$\mathbf{x}$: At every $\mathbf{x}$, there is a set of basis functions that
define $h(\mathbf{x},\mathbf{u})$.

We summarize our findings here:  What is gathering? We apply spatially
invarying kernels first, and then combine the results with weights. •
Consistent with the convolution we learned in Oppenheim and Wilsky [20],
“flip, shift, and integrate”. • Equivalent to the “convolution” in deep neural
networks.

###  II-B Understanding Scattering

We now discuss the scattering equation. For now, we follow the same approach
as we did for the gathering case by analyzing the computation. Substituting
(5) into (2), we see that

$\displaystyle I(\mathbf{x})$ $\displaystyle\overset{\text{by}\,\eqref{eq: I
varying}}{=}\int_{-\infty}^{\infty}h(\mathbf{x},\mathbf{u})J(\mathbf{u})\;d\mathbf{u}$
$\displaystyle\overset{\text{by}\,\eqref{eq:
u}}{=}\int_{-\infty}^{\infty}\left(\sum_{m=1}^{M}a_{\mathbf{u},m}\varphi_{m}(\mathbf{x}-\mathbf{u})\right)J(\mathbf{u})\;d\mathbf{u}$
$\displaystyle=\sum_{m=1}^{M}\left(\int_{-\infty}^{\infty}\varphi_{m}(\mathbf{x}-\mathbf{u})\cdot\Big{(}a_{\mathbf{u},m}J(\mathbf{u})\Big{)}\;d\mathbf{u}\right).$

The last equation is trickier to understand. For a fixed index $m$, the
product term $a_{\mathbf{u},m}J(\mathbf{u})$ is a pixelwise multiplication of
the coefficient map $a_{\mathbf{u},m}$ and the input image $J(\mathbf{u})$.
Therefore, if we express it using the elementwise multiplication $\odot$, we
can write the equation as

$\displaystyle I(\mathbf{x})$ $\displaystyle=\underset{\text{linear
combination of invariant
blurs}}{\underbrace{\sum_{m=1}^{M}\bigg{(}\varphi_{m}\quad\circledast\underset{\text{pixelwise
weighted input}}{\underbrace{\Big{(}a_{m}\odot
J\Big{)}}}\bigg{)}(\mathbf{x})}}.$ (7)

What (7) implies is interesting. We need to first construct the $M$
coefficient maps $\\{a_{\mathbf{u},1},\ldots,a_{\mathbf{u},M}\\}$ where each
$a_{\mathbf{u},m}$ is a value in an array of the same size as the image $J$.
We multiply $J$ with each coefficient map to obtain $M$ “weighted images”
$\widetilde{J}_{m}(\mathbf{u})=a_{\mathbf{u},m}J(\mathbf{u})$. Then we
convolve these $M$ intermediate images through
$(\varphi_{m}\circledast\widetilde{J}_{m})(\mathbf{u})$. Finally, summing over
the $M$ images will give us the final result.

The process is summarized in Figure 3. Comparing it with Figure 1, we notice
that the difference is at where the pixelwise multiplication is located. In
Figure 1, the pixelwise multiplication is done _after_ the convolutions
whereas in Figure 3, it is done _before_ the convolutions. Since the two
processes involve the same number of convolutions and multiplications, the
complexity of executing the computation is identical.

[width=]./pix/Figure03_Scattering_Block_Diagram

Figure 3: Scattering is for Image Formation. Here we show a modified version
of Mao et al’s phase-to-space (P2S) transform [21] by switching the order of
multiplication and convolution. The new formulation follows the derivation in
(7). We argue that this formulation is the better (and the physically valid)
approach when modeling how light propagates from the object plane to the image
plane.

The situation described in (7) is the mirror of Figure 2. Instead of having a
local neighborhood whose pixels are mapped to a common output coordinate, the
input pixel distributes its influence to a local neighborhood in the output
space, as shown in Figure 4. In deep learning, we call it the _transposed
convolution_ , which is typically used to upsample the features in tasks such
as image super-resolution.

At this point, we can give readers a brief intuition of where the scattering
equation comes from. It is based on how light propagates from the source to
the image plane. A more detailed analysis will be performed in the next
section.

Consider a point source in the object plane. We call the point source as a
delta function. When the point source emits light, scalar diffraction theory
[22] tells us that the wave intensity formed on the image plane is the result
of Fresnel diffraction. Assuming a circular aperture and a vacuum medium, the
point spread function (aka blur kernel $h$ or the impulse response) is the
Airy disc. If there are perturbations along the propagation path such as
atmospheric turbulence, then the resulting point spread function will have a
deformed shape caused by the phase distortions. In any case, the impulse
response $h$ tells us how the energy of delta function is redistributed in the
object plane. Thus, for every fixed $\mathbf{u}$, $h(\mathbf{x},\mathbf{u})$
is the value of the point spread function at a running coordinate
$\mathbf{x}$.

[width=0.65]./pix/Figure04_Scattering_Illustration

Figure 4: The operation behind the Scattering equation: A point source
$\delta(\mathbf{u})$ emits light which propagates to the object plane. The
point spread function (aka the blur kernel $h$) seen in the object plane is
centered at $\mathbf{u}$. At each fixed $\mathbf{u}$, the pixel values
$h(\mathbf{x},\mathbf{u})$ are determined by how the propagation medium
affects the amplitude and phase at coordinate $\mathbf{x}$.

The emission of the light from the source to its neighbors is exactly
described by Figure 4. If there is a delta function $\delta(\mathbf{u})$, it
will be mapped to an impulse response $h(\mathbf{x},\mathbf{u})$ centered at
$\mathbf{u}$. (The variable $\mathbf{x}$ is the dummy variable for any
integration to calculate the intensity on the image plane.) Therefore, when
defining the basis representation for $h(\mathbf{x},\mathbf{u})$, the basis
functions and the basis coefficients must be centered at $\mathbf{u}$. This
will give us (5) where again we emphasize that the coefficient
$a_{\mathbf{u},m}$ is indexed by $\mathbf{u}$. The interpretation here is that
for every input coordinate $\mathbf{u}$, there is a set of basis functions
that defines the point spread function $h$. As we can tell from the name, the
point spread function _spreads_ the influence of a point outward. Light does
not gather the pixels and map the sum to an output coordinate.

We summarize our findings here:  What is scattering? We apply weighted
averaging first, and then filter the weighted averages, and finally add. •
Consistent with how light propagates. See Goodman’s _Fourier Optics_ [22] and
our next section. • Equivalent to the “transposed convolution” in deep neural
networks.

###  II-C Conditions for Equivalence

After elaborating on the computations of gathering and scattering, we now
explain the conditions under which the two are equivalent. As expected, the
two are equivalent if and only if the underlying blur is spatially invariant.

The first observation is about the very first equation of convolution defined
in (1). When comparing this equation to the operation illustrated in Figure 2,
we actually made an implicit assumption about the interchangeability of $J$
and $h$:

$\displaystyle\underset{\text{how scattering operates}}{\underbrace{\int
h(\mathbf{x}-\mathbf{u})J(\mathbf{u})d\mathbf{u}}}=\underset{\text{how
gathering operates}}{\underbrace{\int
h(\mathbf{u})J(\mathbf{x}-\mathbf{u})d\mathbf{u}}}.$ (8)

While this may appear subtle (and perhaps we always take it for granted), it
is this interchangeability that allows us to shift-and-crop a region
$J(\mathbf{x}-\mathbf{u})$ from the image $J(\mathbf{u})$ to carry out the
gathering integration efficiently. In fact, the operation defined on the left
hand side of (8) matches with how scattering operates, whereas the operation
on the right hand side matches with how gathering operates. The
interchangeability implies that the two operations are valid as long as the
convolution is shift-invariant so that $h(\mathbf{x},\mathbf{u})$ can be
written as $h(\mathbf{x}-\mathbf{u})$. When $h$ is spatially varying, then we
cannot interchange the order.

Let us look at the equations more carefully through the lens of matrices and
vectors. Let $\mathbf{H}\in\mathbb{R}^{N\times N}$ be the matrix
representation of the spatially varying blur kernel
$h(\mathbf{x},\mathbf{u})$. We assume that there is a set of _circulant_
matrices $\mathbf{H}_{1},\mathbf{H}_{2},\ldots,\mathbf{H}_{M}$ representing
the set of $M$ spatially invariant basis functions
$\\{\varphi_{1},\varphi_{2},\ldots,\varphi_{M}\\}$.

We consider two sets of diagonal matrices. For every index $m$, we define

$\displaystyle\mathbf{D}^{\mathbf{x}}_{m}=\text{diag}\left\\{\begin{bmatrix}a_{\mathbf{x}_{1},m}\\\
\vdots\\\
a_{\mathbf{x}_{N},m}\end{bmatrix}\right\\},\;\;\mathbf{D}^{\mathbf{u}}_{m}=\text{diag}\left\\{\begin{bmatrix}a_{\mathbf{u}_{1},m}\\\
\vdots\\\ a_{\mathbf{u}_{N},m}\end{bmatrix}\right\\}$

Then, the gathering and scattering equations are

$\displaystyle(\text{Gathering}):\qquad\mathbf{H}^{\mathbf{x}}$
$\displaystyle=\sum_{m=1}^{M}\mathbf{D}_{m}^{\mathbf{x}}\mathbf{H}_{m}$ (9)
$\displaystyle(\text{Scattering}):\qquad\mathbf{H}^{\mathbf{u}}$
$\displaystyle=\sum_{m=1}^{M}\mathbf{H}_{m}\mathbf{D}_{m}^{\mathbf{u}}$ (10)

In other words, the difference lies in how we order the diagonal matrices and
the spatially invariant convolution matrices.

It is not difficult to show that the two constructions can never be the same
unless $\mathbf{D}^{\mathbf{x}}_{m}$ and $\mathbf{D}^{\mathbf{u}}_{m}$ are
identity matrices up to a scalar multiple. To see this, we consider the case
where $M=1$.

###### Theorem 1

Let $\mathbf{H}\in\mathbb{R}^{N\times N}$ be a square matrix. Let
$\mathbf{A}=\text{diag}[a_{1},\ldots,a_{N}]$ and
$\mathbf{B}=\text{diag}[b_{1},\ldots,b_{N}]$ be two diagonal matrices. Then,
$\mathbf{A}\mathbf{H}=\mathbf{H}\mathbf{B}$ if and only if
$\mathbf{A}=\mathbf{B}=\lambda\mathbf{I}$ for some constant $\lambda$ where
$\mathbf{I}$ is the identity matrix.

###### Proof:

We just need to write out the terms. For $\mathbf{A}\mathbf{H}$, we can show
that

$\displaystyle\mathbf{A}\mathbf{H}=\begin{bmatrix}a_{1}h_{11}&a_{1}h_{12}&\ldots&a_{1}h_{1N}\\\
a_{2}h_{21}&a_{2}h_{22}&\ldots&a_{2}h_{2N}\\\ \vdots&\vdots&\ddots&\vdots\\\
a_{N}h_{N1}&a_{N}h_{N2}&\ldots&a_{N}h_{NN}\\\ \end{bmatrix},$

and for $\mathbf{H}\mathbf{B}$, we can show that

$\displaystyle\mathbf{H}\mathbf{B}=\begin{bmatrix}b_{1}h_{11}&b_{2}h_{12}&\ldots&b_{N}h_{1N}\\\
b_{1}h_{21}&b_{2}h_{22}&\ldots&b_{N}h_{2N}\\\ \vdots&\vdots&\ddots&\vdots\\\
b_{1}h_{N1}&b_{2}h_{N2}&\ldots&b_{N}h_{NN}\\\ \end{bmatrix},$

By comparing terms, we can see that the only possibility for
$\mathbf{A}\mathbf{H}=\mathbf{H}\mathbf{B}$ is to require
$\mathbf{A}=\mathbf{B}=\lambda\mathbf{I}$. ∎

The result of the previous theorem implies that if we have a convolution
matrix $\mathbf{H}_{m}$ (which is circulant), for the scattering and gathering
operations to be equivalent, we need

$\mathbf{D}^{\mathbf{x}}_{m}\mathbf{H}_{m}=\mathbf{H}_{m}\mathbf{D}^{\mathbf{u}}_{m},\quad\text{for
all }m.$

Theorem 1 asserts that we need
$\mathbf{D}^{\mathbf{x}}_{m}=\mathbf{D}^{\mathbf{u}}_{m}=\lambda\mathbf{I}$.
But if
$\mathbf{D}^{\mathbf{x}}_{m}=\mathbf{D}^{\mathbf{u}}_{m}=\lambda\mathbf{I}$,
then the underlying blur must be spatially invariant.

Another observation of Theorem 1 is that in general,

$\sum_{m=1}^{M}\mathbf{D}_{m}^{\mathbf{x}}\mathbf{H}_{m}\not=\sum_{m=1}^{M}\mathbf{H}_{m}\mathbf{D}_{m}^{\mathbf{u}}.$
(11)

Therefore, the gathering and scattering equations (4) and (5) are _mutually
exclusive_. If we say that $h(\mathbf{x},\mathbf{u})$ can be exactly
represented by the gathering equation, then there will be an approximation
error when representing $h(\mathbf{x},\mathbf{u})$ using the scattering
equation, and vice versa.

Under what conditions would scattering = gathering? When the underlying blur
is spatially invariant. • Scattering and gathering are mutually exclusive. We
cannot simultaneously have (4) and (5) for a spatially varying blur. If one is
the correct representation, the other will have approximation error.

###  II-D Normalization

When performing a convolution, sometimes it is necessary to ensure that the
image intensity is not amplified or attenuated due to an improper
normalization. For example, in a spatially invariant blur, we almost always
require that

$\int h(\mathbf{u})d\mathbf{u}=1,$

assuming that $h(\mathbf{u})\geq 0$ for all $\mathbf{u}$. Otherwise, if the
integral is less than unity, the resulting (convolved) image will appear to be
dimmer. Translated to matrices and vectors, this is equivalent to
$\mathbf{H}\boldsymbol{1}=\boldsymbol{1}$ for an all-one vector
$\boldsymbol{1}$, assuming that $\mathbf{H}$ is circulant.

Suppose that we have a sequence of spatially invariant blurs
$\mathbf{H}_{1},\mathbf{H}_{2},\ldots,\mathbf{H}_{M}$ satisfying the property
that $\mathbf{H}_{m}\boldsymbol{1}=\boldsymbol{1}$ for all $m$. We want the
diagonal matrices
$\mathbf{D}_{1}^{\mathbf{x}},\mathbf{D}_{2}^{\mathbf{x}},\ldots,\mathbf{D}_{M}^{\mathbf{x}}$
to be defined in such a way that the gathering equation (9) will give us

$\displaystyle\boldsymbol{1}\quad\overset{\text{(we
want)}}{=}\quad\mathbf{H}^{\mathbf{x}}\boldsymbol{1}$
$\displaystyle=\left(\sum_{m=1}^{M}\mathbf{D}_{m}^{\mathbf{x}}\mathbf{H}_{m}\right)\boldsymbol{1}=\sum_{m=1}^{M}\mathbf{D}_{m}^{\mathbf{x}}\boldsymbol{1}.$

Therefore, as long as we can ensure that the sum of the $M$ diagonal matrices
$\\{\mathbf{D}_{m}^{\mathbf{x}}\,|\,m=1,\ldots,M\\}$ is a vector of all one’s,
we are guaranteed to have $\mathbf{H}^{\mathbf{x}}$ to have unit rows.
Converting this into the basis representation, it is equivalent to asking

$\sum_{m=1}^{M}a_{\mathbf{x},m}=1,\qquad\text{for all }m,$ (12)

which is reasonably easy to satisfy. For implementation, if
$\mathbf{H}^{\mathbf{x}}$ does not have rows sum to the unity such that
$\mathbf{H}^{\mathbf{x}}\boldsymbol{1}\not=\boldsymbol{1}$, the simplest
approach is to define a diagonal matrix $\mathbf{D}$ such that
$\mathbf{D}^{-1}\mathbf{H}^{\mathbf{x}}\boldsymbol{1}=\boldsymbol{1}$. From
the derivations above, it is clear that the diagonal matrix should have the
elements

$\mathbf{D}=\text{diag}\left\\{\sum_{m=1}^{M}\mathbf{D}_{m}^{\mathbf{x}}\mathbf{H}_{m}\boldsymbol{1}\right\\}.$
(13)

Therefore, the overall operation applied to an image is

$\widehat{\mathbf{I}}=\text{diag}\left\\{\sum_{m=1}^{M}\mathbf{D}_{m}^{\mathbf{x}}\mathbf{H}_{m}\boldsymbol{1}\right\\}^{-1}\left(\sum_{m=1}^{M}\mathbf{D}_{m}^{\mathbf{x}}\mathbf{H}_{m}\mathbf{J}\right),$
(14)

where $\mathbf{J}\in\mathbb{R}^{N}$ is the vector representation of the input,
and $\widehat{\mathbf{I}}\in\mathbb{R}^{N}$ is the output.

The normalization of the scattering equation is more complicated because the
diagonal matrices do not commute with the spatially invariant blurs, and so we
are not able to simplify the equation:

$\displaystyle\boldsymbol{1}\quad\overset{\text{(we
want)}}{=}\quad\mathbf{H}^{\mathbf{u}}\boldsymbol{1}$
$\displaystyle=\sum_{m=1}^{M}\mathbf{H}_{m}(\mathbf{D}_{m}^{\mathbf{u}}\boldsymbol{1}).$

The only work around solution is to define a diagonal matrix $\mathbf{D}$ such
that $\mathbf{D}^{-1}\mathbf{H}^{\mathbf{u}}\boldsymbol{1}=\boldsymbol{1}$.
This would require that

$\mathbf{D}=\text{diag}\left\\{\sum_{m=1}^{M}\mathbf{H}_{m}\mathbf{D}_{m}^{\mathbf{u}}\boldsymbol{1}\right\\}.$
(15)

In other words, when computing the scattering equation and if normalization is
required, then we need to compute the set of convolutions _twice_ : once for
the input, and once for the normalization constant:

$\widehat{\mathbf{I}}=\text{diag}\left\\{\sum_{m=1}^{M}\mathbf{H}_{m}\mathbf{D}_{m}^{\mathbf{x}}\boldsymbol{1}\right\\}^{-1}\left(\sum_{m=1}^{M}\mathbf{H}_{m}\mathbf{D}_{m}^{\mathbf{x}}\mathbf{J}\right).$
(16)

To summarize the normalization:  If normalization is needed, then For
gathering: Either ensure $\sum_{m=1}^{M}a_{\mathbf{x},m}=1$ for all
$\mathbf{x}$, or perform the calculation in (14). • For scattering: Perform
the calculation in (16). The cost is twice the number of convolutions than
gathering.

##  III Origin of the Scattering Equation

In this section, we explain the origin of the scattering equation from a
physics point of view. To make our discussions concrete, the running example
we use is a point source propagating a random medium as shown in Figure 5.

[width=0.95]./pix/Figure05_Fourier_Optics

Figure 5: The coordinate system of a typical optical system, in the presence
of a random medium.

When we discuss optics, the coordinate system can sometimes be confusing. For
the purpose of our discussion, we follow the coordinate system defined Figure
5. The object plane uses the coordinate $\mathbf{u}\in\mathbb{R}^{2}$. We can
think of it as the _input_ coordinate. As light propagates through the random
medium, it reaches at the aperture of a lens. The coordinate on the lens is
denoted by $\boldsymbol{\rho}\in\mathbb{R}^{2}$. When the light reaches the
image plane, the coordinate becomes $\mathbf{x}\in\mathbb{R}^{2}$, which is
also the _output_ coordinate.

Deriving the PSF equation from the Rayleigh-Sommerfeld integral will be too
lengthy for our paper. Therefore, we skip the derivation and refer the readers
to [22, Ch. 4 & 5]. What we shall do here is to highlight the four components
of Figure 5.

Source to Aperture: The propagation of a point source from the object plane to
the aperture, in the absence of the random medium, is characterized by the
free-space propagation. The electromagnetic field defined upon the aperture is
given by [22, Eq 5-25]

$U(\mathbf{u},\boldsymbol{\rho})=\frac{1}{j\lambda
z_{1}}\exp\left\\{j\frac{k}{2z_{1}}|\mathbf{u}-\boldsymbol{\rho}|^{2}\right\\},$
(17)

where $k=2\pi/\lambda$ is the wave number, and $z_{1}$ is the distance from
the source to the aperture. The notation $|\mathbf{u}-\boldsymbol{\rho}|$
denotes the Euclidean distance between the two coordinates $\mathbf{u}$ and
$\boldsymbol{\rho}$. This equation describes a parabolic wave propagating
outward from $\mathbf{u}$. The farther apart $\mathbf{u}$ and
$\boldsymbol{\rho}$ is, the weaker the field $U(\mathbf{u},\boldsymbol{\rho})$
will become.

Aperture and lens: Right at the lens, the incident field will be imparted by
the pupil function of the lens and its phase response. For a lens with a focal
length of $f$, the field at the exit of the aperture is [22, Eq 5-26]

$U^{\prime}(\mathbf{u},\boldsymbol{\rho})=U(\mathbf{u},\boldsymbol{\rho})P(\boldsymbol{\rho})\exp\left\\{j\frac{k}{2f}|\boldsymbol{\rho}|^{2}\right\\},$
(18)

where $P(\boldsymbol{\rho})$ is the pupil function, typically chosen to be a
circular indicator function.

Aperture to image: When the incident field exits the lens, it propagates via
Fresnel diffraction to the image plane. Referring to [22, Eq 5-27], we can
show that

$\underset{=h(\mathbf{x},\mathbf{u})}{\underbrace{U^{\prime\prime}(\mathbf{x},\mathbf{u})}}=\frac{1}{j\lambda
z_{2}}\int_{-\infty}^{\infty}U^{\prime}(\mathbf{u},\boldsymbol{\rho})\exp\left\\{j\frac{k}{2z_{2}}|\mathbf{x}-\boldsymbol{\rho}|^{2}\right\\}\,d\boldsymbol{\rho}.$
(19)

Notice that the final electromagnetic field
$U^{\prime\prime}(\mathbf{x},\mathbf{u})$ arriving at the image plane is
originated from a point source. As such,
$U^{\prime\prime}(\mathbf{x},\mathbf{u})$ is the point spread function
$h(\mathbf{x},\mathbf{u})$.

The PSF $h(\mathbf{x},\mathbf{u})$ can be expressed (with some approximation)
as [22, Eq 5-36]:

$h(\mathbf{x},\mathbf{u})=\underset{=h(\mathbf{x}-\mathbf{u})\,\,\text{if
$S=1$}}{\underbrace{\frac{1}{\lambda^{2}z_{1}z_{2}}\int_{-\infty}^{\infty}P(\boldsymbol{\rho})\exp\left\\{-j\frac{k}{z_{2}}(\mathbf{x}-S\mathbf{u})^{T}\boldsymbol{\rho}\right\\}d\boldsymbol{\rho}}},$
(20)

where $S=-z_{2}/z_{1}$ is the magnification factor. If, for simplicity, we
assume $z_{1}=-z_{2}$ so that $S=1$, then $h(\mathbf{x},\mathbf{u})$ is
completely characterized by the coordinate difference $\mathbf{x}-\mathbf{u}$.
This will give us $h(\mathbf{x},\mathbf{u})=h(\mathbf{x}-\mathbf{u})$, and so
$h(\mathbf{x},\mathbf{u})$ represents a spatially _invariant_ kernel.

Random medium. The fourth element we need to discuss, which is also the source
of the problem, is the random medium. The random medium introduces a random
amplitude and phase distortion as

$R_{\mathbf{u}}(\boldsymbol{\rho})=\underset{\text{amplitude}}{\underbrace{A_{\mathbf{u}}(\boldsymbol{\rho})}}\;\times\;\underset{\text{phase}}{\underbrace{\exp\left\\{-j\phi_{\mathbf{u}}(\boldsymbol{\rho})\right\\}}}.$
(21)

Notice that in this definition, the distortion has a coordinate pair
$(\mathbf{u},\boldsymbol{\rho})$. In the case of a random medium, the position
of the source and aperture will parameterize the slice of the atmosphere the
wave can be said to propagate through. Therefore, the distortion must be
defined by these two coordinates. Taking the position of our imaging system as
our reference, the distortions can be completely parameterized by
$\mathbf{u}$.

The impact of the random medium takes place at source-to-aperture, i.e.,
appending (21) to (17). Thus, the electromagnetic field incident upon the
aperture is

$\displaystyle U(\mathbf{u},\boldsymbol{\rho})$
$\displaystyle=\frac{1}{j\lambda
z_{1}}A_{\mathbf{u}}(\boldsymbol{\rho})\exp\left\\{-j\phi_{\mathbf{u}}(\boldsymbol{\rho})\right\\}$
$\displaystyle\qquad\times\exp\left\\{j\frac{k}{2z_{1}}|\mathbf{u}-\boldsymbol{\rho}|^{2}\right\\}.$

Consequently, the point spread function is

$\displaystyle h(\mathbf{x},\mathbf{u})$
$\displaystyle=\frac{1}{\lambda^{2}z_{1}z_{2}}\int_{-\infty}^{\infty}\underset{=g_{\mathbf{u}}(\boldsymbol{\rho})}{\underbrace{A_{\mathbf{u}}(\boldsymbol{\rho})\exp\left\\{-j\phi_{\mathbf{u}}(\boldsymbol{\rho})\right\\}P(\boldsymbol{\rho})}}$
$\displaystyle\qquad\times\exp\left\\{-j\frac{k}{z_{2}}(\mathbf{x}-\mathbf{u})^{T}\boldsymbol{\rho}\right\\}d\boldsymbol{\rho}$
(22)

where again we assumed $S=1$. Defining the constant
$\kappa=1/(\lambda^{2}z_{1}z_{2})$, we see that (22) takes the form of

$h(\mathbf{x},\mathbf{u})=\kappa\int_{-\infty}^{\infty}g_{\mathbf{u}}(\boldsymbol{\rho})e^{-j2\pi\boldsymbol{\xi}^{T}\boldsymbol{\rho}}d\boldsymbol{\rho},$

where we defined $\boldsymbol{\xi}=(\mathbf{x}-\mathbf{u})/(\lambda z_{2})$
and

$\displaystyle
g_{\mathbf{u}}(\boldsymbol{\rho})\overset{\text{def}}{=}A_{\mathbf{u}}(\boldsymbol{\rho})\exp\left\\{-j\phi_{\mathbf{u}}(\boldsymbol{\rho})\right\\}P(\boldsymbol{\rho}).$

Therefore, (22) can be seen as the Fourier transform of the random medium and
the pupil function via

$\displaystyle h(\mathbf{x},\mathbf{u})$
$\displaystyle=\text{Fourier}\Big{\\{}g_{\mathbf{u}}(\boldsymbol{\rho})\Big{\\}}\Big{|}_{\frac{\mathbf{x}-\mathbf{u}}{\lambda
z_{2}}},$ (23)

where we specify that the transform is evaluated at the coordinate
$(\mathbf{x}-\mathbf{u})/(\lambda z_{2})$.

At this point, it should become clear that the PSF $h(\mathbf{x},\mathbf{u})$
must be a function of $\mathbf{x}-\mathbf{u}$. However, since
$g_{\mathbf{u}}(\boldsymbol{\rho})$ is also indexed by $\mathbf{u}$, the
resulting PSF $h(\mathbf{x},\mathbf{u})$ should inherent the index
$\mathbf{u}$. This will give us

$\displaystyle h(\mathbf{x},\mathbf{u})$ $\displaystyle=\text{some function of
$(\mathbf{x}-\mathbf{u})$, index by $\mathbf{u}$},$
$\displaystyle=\sum_{m=1}^{M}a_{\mathbf{u},m}\varphi_{m}(\mathbf{x}-\mathbf{u}).$
(24)

where in the last step we use the linear combination of basis functions as the
model for such $h(\mathbf{x},\mathbf{u})$. The basis function $\varphi_{m}$
here captures the spatial invariance, whereas the coefficients
$a_{\mathbf{u},m}$ capture the spatially varying indices.

The implication of our derivation is important as it explains why in optics
simulation, such as imaging through random medium, must follow the
_scattering_ equation if we choose to represent the PSF using a set of
spatially invariant basis functions. This is due to the fact that nature
relies on the _source_ location to determine the response of the system.
Though our derivation relied upon a random medium, this concept extends to
cases such as spherical aberrations or defocus blur in a scene which has
objects which are beyond the depth of field, both of which may be formulated
as a phase error. In both of these examples, the source location is what
parameterizes the phase error.

##  IV Do These Actually Matter?

The question to ask now is: Given the gathering equation and the scattering
equation, does it really matter if we choose the “wrong” one? The goal of this
section is to answer the question through a few examples.

###  IV-A Scattering Works for Optical Simulation

In the first example, we consider the problem of simulating the resulting
image given a specific light source. The light source $J(\mathbf{u})$ we
consider here consists of two delta functions:

$J(\mathbf{u})=\delta(\mathbf{u}+\Delta)+\delta(\mathbf{u}-\Delta),$

where $\Delta$ is a small displacement from $\mathbf{u}$. For convenience, we
consider a plane with two halves. The separation is located at the origin; any
pixel that is on the left is denoted as “$\mathbf{u}<0$” (this notation means
that the _horizontal_ component of $\mathbf{u}$ is less than zero). Similarly,
any pixel that is on the right is denoted as “$\mathbf{u}\geq 0$”. Therefore,
$\delta(\mathbf{u}+\Delta)$ is on the left, and $\delta(\mathbf{u}-\Delta)$ is
on the right.

Imagine that in front of the light source, we put two transparent sheets with
different phase profiles (which can be engineered using a meta-material.) This
will give us a spatially varying blur kernel $h(\mathbf{x},\mathbf{u})$, and
for simplicity, if light is emitted on the left hand side, then the blur uses
a smaller radius; if the light is emitted on the right hand side, then the
blur uses a larger radius. Thus, we write

$h(\mathbf{x},\mathbf{u})=\begin{cases}\frac{1}{2\pi\sigma_{1}^{2}}\exp\left\\{-\frac{\|\mathbf{x}-\mathbf{u}\|^{2}}{2\sigma_{1}^{2}}\right\\}\overset{\text{def}}{=}\varphi_{1}(\mathbf{x}-\mathbf{u}),&\,\mathbf{u}<0,\\\
\frac{1}{2\pi\sigma_{2}^{2}}\exp\left\\{-\frac{\|\mathbf{x}-\mathbf{u}\|^{2}}{2\sigma_{2}^{2}}\right\\}\overset{\text{def}}{=}\varphi_{2}(\mathbf{x}-\mathbf{u}),&\,\mathbf{u}\geq
0,\end{cases}$ (25)

where $\sigma_{1}<\sigma_{2}$. Figure 6(b) illustrates these spatially varying
blur kernels. For visualization purposes, we show only the PSFs at a grid of
points. In reality, the PSFs is defined continuously over $\mathbf{u}$.

[width=0.48]./pix/Figure06_Gathering_Example_J |  [width=0.48]./pix/Figure06_Gathering_Example_h
---|---
(a) Input $J(\mathbf{u})$ | (b) Blurs $h(\mathbf{x},\mathbf{u})$
| (illustrated over a grid of points)
Figure 6: Visualization of an example with (a) $J(\mathbf{u})$ and (b) a grid
of spatially varying blur kernels.

Before we do any calculus, we can perform a thought experiment. Figure 7
illustrates a hypothetical experimental setup. On the object plane there are
two points emitting light through a meta surface with two different phase
profiles. As the light propagates outward from the source through diffraction,
the electromagnetic fields superimpose over each other. When the light reaches
the aperture, the two diffraction patterns overlap. Therefore, the resulting
image $I(\mathbf{x})$, without any calculation, should be one big diffraction
pattern. It is impossible to obtain a sharp cutoff and two diffraction
patterns.

[width=]./pix/Figure07_Gathering_Example_Illustration

Figure 7: Thought experiment with two points on the object plane, diffracting
through two different metasurfaces. The resulting image should, in principle,
be one superimposed diffraction pattern.

With this thought experiment in mind, we can now talk about the gathering and
the scattering equations. For the gathering equation, since the PSF has a
simple binary structure, we can define it as

$\displaystyle h^{\text{gather}}(\mathbf{x},\mathbf{u})$
$\displaystyle=\underset{=a_{\mathbf{x},1}}{\underbrace{\mathbb{I}\\{\mathbf{x}\in\text{left}\\}}}\times\varphi_{1}(\mathbf{x}-\mathbf{u})$
$\displaystyle\qquad\underset{=a_{\mathbf{x},2}}{\underbrace{\mathbb{I}\\{\mathbf{x}\in\text{right}\\}}}\times\varphi_{2}(\mathbf{x}-\mathbf{u}),$
(26)

where $\mathbb{I}\\{\cdot\\}$ is an indicator function For scattering, the
equation takes a similar form

$\displaystyle h^{\text{scatter}}(\mathbf{x},\mathbf{u})$
$\displaystyle=\underset{=a_{\mathbf{u},1}}{\underbrace{\mathbb{I}\\{\mathbf{u}\in\text{left}\\}}}\times\varphi_{1}(\mathbf{x}-\mathbf{u})$
$\displaystyle\qquad\underset{=a_{\mathbf{u},2}}{\underbrace{\mathbb{I}\\{\mathbf{u}\in\text{right}\\}}}\times\varphi_{2}(\mathbf{x}-\mathbf{u}),$
(27)

where we replaced $a_{\mathbf{x},m}$ by $a_{\mathbf{u},m}$.

By comparing the gathering equation (26) and the scattering equation (27) with
the original spatially varying $h(\mathbf{x},\mathbf{u})$ in (25), it is clear
that only the scattering equation will match with the original
$h(\mathbf{x},\mathbf{u})$ because they are both indexed by $\mathbf{u}$.
However, to confirm that this is indeed the case, it would be helpful to look
at the resulting image, as illustrated in Figure 8. Let us explain how these
figures are obtained.

[width=0.48]./pix/Figure08_Gathering_Example_Igather |  [width=0.48]./pix/Figure08_Gathering_Example_Iscatter
---|---
(a) $I^{\text{gather}}(\mathbf{x})$ | (b) $I^{\text{scatter}}(\mathbf{x})$
Figure 8: Comparison between gathering and scattering for the setup in Figure
7. Notice that for this optical experiment, we should expect the resulting
image to contain one big diffraction pattern. However, only the scattering
equation demonstrates this.

The resulting image of the gathering equation can be shown as follows:

$\displaystyle
I^{\text{gather}}(\mathbf{x})=\int\left(a_{\mathbf{x},1}\varphi_{1}(\mathbf{x}-\mathbf{u})+a_{\mathbf{x},2}\varphi_{2}(\mathbf{x}-\mathbf{u})\right)J(\mathbf{u})\,d\mathbf{u}.$

Since $a_{\mathbf{x},1}=1$ if $\mathbf{x}\geq 0$ and $a_{\mathbf{x},1}=0$ if
$\mathbf{x}<0$ (similarly for $a_{\mathbf{x},2}$), the coefficients
$a_{\mathbf{x},1}$ and $a_{\mathbf{x},2}$ will create two cases as

$\displaystyle I^{\text{gather}}(\mathbf{x})$
$\displaystyle=\begin{cases}\int_{-\infty}^{\infty}\varphi_{1}(\mathbf{x}-\mathbf{u})J(\mathbf{u})\,d\mathbf{u},&\quad\mathbf{x}<0,\\\
\int_{-\infty}^{\infty}\varphi_{2}(\mathbf{x}-\mathbf{u})J(\mathbf{u})\,d\mathbf{u},&\quad\mathbf{x}\geq
0.\end{cases}$
$\displaystyle=\begin{cases}\int_{-\infty}^{\infty}\varphi_{1}(\mathbf{x}-\mathbf{u})(\delta(\mathbf{u}+\Delta)+\delta(\mathbf{u}-\Delta))\,d\mathbf{u},&\quad\mathbf{x}<0,\\\
\int_{-\infty}^{\infty}\varphi_{2}(\mathbf{x}-\mathbf{u})(\delta(\mathbf{u}+\Delta)+\delta(\mathbf{u}-\Delta))\,d\mathbf{u},&\quad\mathbf{x}\geq
0.\end{cases}$
$\displaystyle=\begin{cases}\varphi_{1}(\mathbf{x}+\Delta)+\varphi_{1}(\mathbf{x}-\Delta),&\quad\mathbf{x}<0,\\\
\varphi_{2}(\mathbf{x}+\Delta)+\varphi_{2}(\mathbf{x}-\Delta),&\quad\mathbf{x}\geq
0.\end{cases}$ (28)

If we draw $I^{\text{gather}}(\mathbf{x})$, we will obtain the figure shown in
Figure 8(a).

For scattering, we can carry out the same derivation and show that

$\displaystyle
I^{\text{scatter}}(\mathbf{x})=\int_{-\infty}^{\infty}\left(a_{\mathbf{u},1}\varphi_{1}(\mathbf{x}-\mathbf{u})+a_{\mathbf{u},2}\varphi_{2}(\mathbf{x}-\mathbf{u})\right)J(\mathbf{u})\,d\mathbf{u}$
$\displaystyle=\int_{-\infty}^{\infty}\varphi_{1}(\mathbf{x}-\mathbf{u})[a_{\mathbf{u},1}J(\mathbf{u})]\,d\mathbf{u}$
$\displaystyle\qquad\qquad+\int_{-\infty}^{\infty}\varphi_{2}(\mathbf{x}-\mathbf{u})[a_{\mathbf{u},2}J(\mathbf{u})]\,d\mathbf{u}$
$\displaystyle=\int_{-\infty}^{\infty}\varphi_{1}(\mathbf{x}-\mathbf{u})\delta(\mathbf{u}+\Delta)\,d\mathbf{u}$
$\displaystyle\qquad\qquad+\int_{-\infty}^{\infty}\varphi_{2}(\mathbf{x}-\mathbf{u})\delta(\mathbf{u}-\Delta)\,d\mathbf{u}$
$\displaystyle=\varphi_{1}(\mathbf{x}+\Delta)+\varphi_{2}(\mathbf{x}-\Delta),$
(29)

where we used the fact that
$a_{\mathbf{u},1}=\mathbb{I}\\{\mathbf{u}\in\text{left}\\}$ and so
$a_{\mathbf{u},1}\delta(\mathbf{u}+\Delta)=\delta(\mathbf{u}+\Delta)$ whereas
$a_{\mathbf{u},2}\delta(\mathbf{u})=0$. Similarly we have
$a_{\mathbf{u},1}\delta(\mathbf{u}-\Delta)=0$ and
$a_{\mathbf{u},2}\delta(\mathbf{u}-\Delta)=\delta(\mathbf{u}-\Delta)$.

If we draw the resulting image, we will obtain the figure shown in Figure
8(b). This is consistent with what we expect in Figure 7 and the theoretical
derivation in Section III.

###  IV-B Gathering Works for Image Filtering

In the second example, we consider the problem of _image filtering_. The
scenario is that we are given a noisy image $J(\mathbf{u})$ that contains two
regions:

$\displaystyle J(\mathbf{u})$
$\displaystyle=\begin{cases}\theta_{1}+\text{Gauss}(0,\sigma_{1}^{2}),&\qquad\mathbf{u}<0,\\\
\theta_{2}+\text{Gauss}(0,\sigma_{2}^{2}),&\qquad\mathbf{u}\geq 0.\end{cases}$
$\displaystyle=(\theta_{1}+W_{1}(\mathbf{u}))\times(1-\text{Step}(\mathbf{u}))$
$\displaystyle\qquad\qquad\qquad+(\theta_{2}+W_{2}(\mathbf{u}))\times\text{Step}(\mathbf{u})$
(30)

with two signal levels $\theta_{1}$ and $\theta_{2}$, and two noise standard
deviations $\sigma_{1}$ and $\sigma_{2}$ such that $\sigma_{1}>\sigma_{2}$. In
this equation, $W_{1}(\mathbf{u})\sim\text{Gauss}(0,\sigma_{1}^{2})$ and
$W_{2}(\mathbf{u})\sim\text{Gauss}(0,\sigma_{2}^{2})$ denote the Gaussian
noise. The function $\text{Step}(\mathbf{u})$ denotes the horizontal step
function where $\text{Step}(\mathbf{u})=1$ for any $\mathbf{u}\geq 0$, i.e.,
residing on the right, and $\text{Step}(\mathbf{u})=0$ for any $\mathbf{u}<0$.
For illustration, we show in Figure 9(a) the case where $\theta_{1}=0.8$,
$\theta_{2}=0.2$, $\sigma_{1}=0.1$ and $\sigma_{2}=0.02$.

[width=0.48]./pix/Figure09_Scattering_Example_J |  [width=0.48]./pix/Figure09_Scattering_Example_h
---|---
(a) Input $J(\mathbf{u})$ | (b) Denoising filter $h(\mathbf{x},\mathbf{u})$
| (illustrated over a grid of points)
Figure 9: Thought experiment of two noisy regions in an image. To denoise this
image, ideally we would want to apply to different filters with a sharp
boundary at the transition.

To denoise this image, we consider the simplest approach assuming that we
_knew_ the partition of the two regions. Suppose that we want to denoise the
left side. Since we know that the noise is stronger, we shall apply a stronger
filter. As illustrated in Figure 10, for this filter to be effective along the
boundary, we should apply a mask _after_ the filtering is done.

[width=]./pix/Figure10_Scattering_Example_Illustration

Figure 10: Thought experiment with two noisy half-planes. As we perform the
denoising step, we would hope that the sharp boundary is preserved.

The spatially varying filter we propose here takes the form

$\displaystyle h(\mathbf{x},\mathbf{u})$
$\displaystyle=\begin{cases}\frac{1}{\sqrt{2\pi
s_{1}^{2}}}\exp\left\\{-\frac{\|\mathbf{x}-\mathbf{u}\|^{2}}{2s_{1}^{2}}\right\\}\overset{\text{def}}{=}\varphi_{1}(\mathbf{x}-\mathbf{u}),&\mathbf{x}<0,\\\
\frac{1}{\sqrt{2\pi
s_{2}^{2}}}\exp\left\\{-\frac{\|\mathbf{x}-\mathbf{u}\|^{2}}{2s_{2}^{2}}\right\\}\overset{\text{def}}{=}\varphi_{2}(\mathbf{x}-\mathbf{u}),&\mathbf{x}\geq
0,\end{cases}$

where we assume that $s_{1}>s_{2}$. We are careful about the index in this
equation, remarking that the conditions are applied to $\mathbf{x}$ instead of
$\mathbf{u}$. We will illustrate what will happen if the conditions are
applied to $\mathbf{u}$.

The gathering and the scattering equation for this example are identical to
those in (26) and (27). Most importantly, the coefficients $a_{\mathbf{x},m}$
and $a_{\mathbf{u},m}$ are binary masks indicating whether the pixel
$\mathbf{x}$ (or $\mathbf{u}$) is on the left / right hand side.

The resulting images $I^{\text{gather}}(\mathbf{x})$ and
$I^{\text{scatter}}(\mathbf{x})$ follow a similar derivation as in (28) and
(29). For the gathering equation, we recognize that

$\displaystyle
I^{\text{gather}}(\mathbf{x})=\int_{-\infty}^{\infty}\left(a_{\mathbf{x},1}\varphi_{1}(\mathbf{x}-\mathbf{u})+a_{\mathbf{x},2}\varphi_{2}(\mathbf{x}-\mathbf{u})\right)J(\mathbf{u})\,d\mathbf{u}$
$\displaystyle=a_{\mathbf{x},1}\int_{-\infty}^{\infty}\varphi_{1}(\mathbf{x}-\mathbf{u})\,\Bigl{\\{}\left[\theta_{1}+W_{1}(\mathbf{u})\right](1-\text{Step}(\mathbf{u}))$
$\displaystyle\qquad\qquad\qquad\qquad+\left[\theta_{2}+W_{2}(\mathbf{u})\right]\text{Step}(\mathbf{u})\Bigr{\\}}\,d\mathbf{u}$
$\displaystyle\quad+a_{\mathbf{x},2}\int_{-\infty}^{\infty}\varphi_{2}(\mathbf{x}-\mathbf{u})\,\Bigl{\\{}\left[\theta_{1}+W_{1}(\mathbf{u})\right](1-\text{Step}(\mathbf{u}))$
$\displaystyle\qquad\qquad\qquad\qquad+\left[\theta_{2}+W_{2}(\mathbf{u})\right]\text{Step}(\mathbf{u})\Bigr{\\}}\,d\mathbf{u}$
$\displaystyle=\begin{cases}[((\theta_{1}+W_{1}(\mathbf{u}))(1-\text{Step}(\mathbf{u}))\\\
\qquad+(\theta_{2}+W_{2}(\mathbf{u}))\text{Step}(\mathbf{u}))\circledast\varphi_{1}(\mathbf{u})](\mathbf{x}),&\quad\mathbf{x}<0,\\\
[((\theta_{1}+W_{1}(\mathbf{u}))(1-\text{Step}(\mathbf{u}))\\\
\qquad+(\theta_{2}+W_{2}(\mathbf{u}))\text{Step}(\mathbf{u}))\circledast\varphi_{2}(\mathbf{u})](\mathbf{x}),&\quad\mathbf{x}\geq
0,\end{cases}$ (31)

where the last equality holds because $a_{\mathbf{x},1}$ and
$a_{\mathbf{x},2}$ are indicator functions.

For the scattering equation, we recognize that

$\displaystyle
I^{\text{scatter}}(\mathbf{x})=\int_{-\infty}^{\infty}\left(a_{\mathbf{u},1}\varphi_{1}(\mathbf{x}-\mathbf{u})+a_{\mathbf{u},2}\varphi_{2}(\mathbf{x}-\mathbf{u})\right)J(\mathbf{u})\,d\mathbf{u}$
$\displaystyle\quad=\int_{-\infty}^{\infty}a_{\mathbf{u},1}\varphi_{1}(\mathbf{x}-\mathbf{u})[\theta_{1}+W_{1}(\mathbf{u})](1-\text{Step}(\mathbf{u}))\,d\mathbf{u}$
$\displaystyle\quad\quad+\int_{-\infty}^{\infty}a_{\mathbf{u},2}\varphi_{2}(\mathbf{x}-\mathbf{u})[\theta_{2}+W_{2}(\mathbf{u})]\text{Step}(\mathbf{u})\,d\mathbf{u},$

where we used the facts $a_{\mathbf{u},1}=(1-\text{Step}(\mathbf{u}))$ and
$a_{\mathbf{u},2}=\text{Step}(\mathbf{u})$, and hence
$a_{\mathbf{u},1}\text{Step}(\mathbf{u})=0$ and
$a_{\mathbf{u},2}(1-\text{Step}(\mathbf{u}))=0$. As a result, we can simplify
the above equations as

$\displaystyle I^{\text{scatter}}(\mathbf{x})$
$\displaystyle=[((\theta_{1}+W_{1}(\mathbf{u}))(1-\text{Step}(\mathbf{u})))\circledast\varphi_{1}(\mathbf{u})](\mathbf{x})$
$\displaystyle\quad+[((\theta_{2}+W_{2}(\mathbf{u}))\text{Step}(\mathbf{u}))\circledast\varphi_{2}(\mathbf{u})](\mathbf{x}).$
(32)

So, the result is the _sum_ of two convolutions with the input image
$J(\mathbf{u})$.

We visualize the results in Figure 11. As we can see here, while both of them
can offer denoising to some extent, the gathering approach handles the
boundary much better because the masking is performed _after_ the filtering.
If the masking is performed before the filtering, then (32) tells us that we
are summing two convolutions of the same image. Therefore, the edge is
blurred.

[width=0.48]./pix/Figure11_Scattering_Example_Igather |  [width=0.48]./pix/Figure11_Scattering_Example_Iscatter
---|---
(a) $I^{\text{gather}}(\mathbf{x})$ | (b) $I^{\text{scatter}}(\mathbf{x})$
Figure 11: Comparison between gathering and scattering for the setup in Figure
9. Notice that for this denoising experiment, the better method should produce
a sharp transition along the boundary.

##  V Discussions

As we write this paper, two underpinning questions are constantly asked: (1)
what utility does it bring? (2) How does it affect how we solve an inverse
problem? In this section, we briefly share our findings.

###  V-A Utility: Simulating atmospheric turbulence

One of the biggest utilities (for the scattering approach) is the simulation
of atmospheric turbulence. The latest turbulence simulators, based on phase-
over-aperture [23], phase-to-space (P2S) transform [21], and dense-field P2S
[24], are all using the _gathering_ equation. This present paper explains why
these prior simulators have a potential hindsight when decomposing the
spatially varying blur into spatially invariant blurs, which should have been
using the scattering equation instead.

Table I shows a comparison between the dense-field P2S simulator [24] (which
is based on the gathering equation), and a new simulator implemented using the
scattering equation.111The new simulator contains a few other modifications
including expanding the parameter space, and expanding the kernel supports.
However, the biggest change is the adoption of the scattering equation. Our
testing dataset is based on the text recognition dataset released to the UG2+
challenge [25]. We tested three image reconstruction models: TSRWGAN [26],
ESTRNN [27], and TMT [28]. For each model, we report the recognition accuracy
in terms of CRNN / DAN / ASTER from the restored images.

TABLE I: Comparison on the turbulence-text dataset. The metrics are
recognition accuracy (%) by text recognition models CRNN / DAN / ASTER from
the restored images.

0.46!  Model Simulator in [24] New simulator TSRWGAN 57.75 / 71.45 / 73.10
60.30 / 73.90 / 74.40 ESTRNN 84.50 / 96.25 / 95.45 87.10 / 97.80 / 96.95 TMT
79.25 / 85.20 / 88.00 80.90 / 87.25 / 88.55

As we can see in this table, the new simulator indeed offers a sizable amount
of improvement over the previous simulator. Since the change from the
gathering to the scattering equation attributes to a substantial portion of
the simulator update, the utility of our study is evident.

###  V-B We have an inverse problem, which model to use?

Solving an inverse often requires an optimization. For a spatially varying
blur problem, the typical formulation is

$\widehat{\mathbf{J}}=\mathop{\underset{\mathbf{J}}{\mbox{argmin}}}\;\;\|\mathbf{I}-\mathbf{H}\mathbf{J}\|^{2}+\lambda
R(\mathbf{J}),$ (33)

for some regularization functions $R(\mathbf{J})$. However, if we want to be
exact about the forward model, we need to go back to the light propagation
physics and use the scattering equation to represent $\mathbf{H}$. This will
give us

$\widehat{\mathbf{J}}=\mathop{\underset{\mathbf{J}}{\mbox{argmin}}}\;\;\left\|\mathbf{I}-\left(\sum_{m=1}^{M}\mathbf{H}_{m}\mathbf{D}_{m}^{\mathbf{u}}\right)\mathbf{J}\right\|^{2}+\lambda
R(\mathbf{J}).$ (34)

This problem is extremely difficult to solve because even if
$\mathbf{D}_{m}^{\mathbf{u}}$ is binary so that it forms a partition of the
image, we cannot solve for individual regions and combine them. Moreover, we
cannot solve for individual blurs because the $\mathbf{H}_{m}$’s are summed.
As a result, the decomposition offers little computational benefit although it
accurately reflects the physics.

On the other hand, the gathering equation will give us

$\widehat{\mathbf{J}}=\mathop{\underset{\mathbf{J}}{\mbox{argmin}}}\;\;\left\|\mathbf{I}-\left(\sum_{m=1}^{M}\mathbf{D}_{m}^{\mathbf{x}}\mathbf{H}_{m}\right)\mathbf{J}\right\|^{2}+\lambda
R(\mathbf{J}).$ (35)

If $\mathbf{D}_{m}^{\mathbf{x}}$ is binary (as in Nagy and O’Leary [18]), we
can partition the image into $M$ smaller regions and solve them individually.
The parallelism offered by the model is computationally appealing, and it has
been proven useful to some extent. The caveat is that, of course, the
gathering equation does not match the physics, and so we only solve a proxy to
the original optimization. However, in the case of a binary
$\mathbf{D}_{m}^{\mathbf{x}}$, there will be regions in which the two forms
will be identical provided that the support of the impulse response is
contained within the region spanned by $\mathbf{D}_{m}^{\mathbf{x}}$. In other
words, if a region can be approximated thought of as spatially invariant, the
two solutions in this region will be identical. In a general sense for
arbitrary $\mathbf{D}_{m}^{\mathbf{x}}$, any solution we obtain, regardless if
we can prove convergence of the optimization algorithm, will not resemble the
true solution to (33).

Our advice to practitioners, when solving an inverse problem related to the
spatially varying blur, is to have the sense of awareness for this kind of
mismatch. For deep neural networks, the mismatch is often less of an issue
when the capacity of the network is large enough. However, the training data
needs to capture enough of the physics in order to generalize well.

##  VI Conclusion

The gathering and the scattering representation of a spatially varying blur
are twins of the same blur kernel. They are mutually exclusive, in the sense
that if one is the exact representation of the original blur, the other one
can only be an approximation. They become identical if the underlying blur
kernel is spatially invariant. As summary, we recognize the following key
points:

Gathering has an origin of image filtering. It is an effective approach to
speed up the spatially varying blur via a small set of invarying blurs. The
approach is to filter the image first, and then combine the filtered images
through a pixelwise mask. Gathering offers better edge awareness in tasks such
as image denoising.

Scattering is originated from light propagation physics. It is an accurate
description according to the scalar diffraction theory in Fourier optics. The
approach is to weigh the images and then perform filtering afterward.
Scattering is _the_ model for describing how light propagates through a random
medium.

Our hope in this paper is to raise awareness about the different meanings and
implications of scattering and gathering. At the very least, gathering and
scattering are not just the convolution and transposed convolution in deep
neural networks.

## References

  * [1] S. H. Chan,  “Constructing a sparse convolution matrix for shift varying image restoration problems,”  in 2010 IEEE International Conference on Image Processing, 2010, pp. 3601–3604.
  * [2] A. Chakrabarti, T. Zickler, and W. T. Freeman,  “Analyzing spatially-varying blur,”  in 2010 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2010, pp. 2512–2519.
  * [3] X. Zhu, S. Cohen, S. Schiller, and P. Milanfar,  “Estimating spatially varying defocus blur from a single image,”  IEEE Transactions on Image Processing, vol. 22, no. 12, pp. 4879–4891, 2013.
  * [4] Y. Huang, É. Chouzenoux, and V. Elvira,  “Probabilistic modeling and inference for sequential space-varying blur identification,”  IEEE Transactions on Computational Imaging, vol. 7, pp. 531–546, 2021.
  * [5] M. M. Sondhi,  “Image restoration: The removal of spatially invariant degradations,”  Proceedings of the IEEE, vol. 60, no. 7, pp. 842–853, 1972.
  * [6] M. Suin, K. Purohit, and A. N. Rajagopalan,  “Degradation aware approach to image restoration using knowledge distillation,”  IEEE Journal of Selected Topics in Signal Processing, vol. 15, no. 2, pp. 162–173, 2021.
  * [7] D. O’Connor and L. Vandenberghe,  “Total variation image deblurring with space-varying kernel,”  Computational Optimization and Applications, vol. 67, no. 3, pp. 521–541, Jul 2017.
  * [8] B. Saleh and N. Subotic,  “Time-variant filtering of signals in the mixed time frequency domain,”  IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 33, no. 6, pp. 1479–1485, 1985.
  * [9] I. Sodagar, K. Nayebi, and T. P. Barnwell,  “Time-varying filter banks and wavelets,”  IEEE Transactions on Signal Processing, vol. 42, no. 11, pp. 2983–2996, 1994.
  * [10] I. Sodagar, K. Nayebi, T. P. Barnwell, and M. J. T. Smith,  “Time-varying analysis-synthesis systems based on filter banks and post filtering,”  IEEE Transactions on Signal Processing, vol. 43, no. 11, pp. 2512–2524, 1995.
  * [11] F. Wang and V. Balakrishnan,  “Robust kalman filters for linear time-varying systems with stochastic parametric uncertainties,”  IEEE Transactions on Signal Processing, vol. 50, no. 4, pp. 803–813, 2002.
  * [12] M. S. Arulampalam, S. Maskell, N. Gordon, and T. Clapp,  “A tutorial on particle filters for online nonlinear/non-gaussian bayesian tracking,”  IEEE Transactions on Signal Processing, vol. 50, no. 2, pp. 174–188, 2002.
  * [13] G. Wang,  “The most general time-varying filter bank and time-varying lapped transforms,”  IEEE Transactions on Signal Processing, vol. 54, no. 10, pp. 3775–3789, 2006.
  * [14] T. W. Mitchel, B. Brown, D. Koller, T. Weyrich, S. Rusinkiewicz, and M. Kazhdan,  “Efficient spatially adaptive convolution and correlation,”  Tech. Rep. 2006.13188, arXiv preprint, Jun 2020.
  * [15] K. Yanny, K. Monakhova, R. W. Shuai, and L. Waller,  “Deep learning for fast spatially varying deconvolution,”  Optica, vol. 9, no. 1, pp. 96–99, Jan 2022.
  * [16] M. Arigovindan, J. Shaevitz, J. McGowan, J. W. Sedat, and D. A. Agard,  “A parallel product-convolution approach for representing depth varying point spread functions in 3d widefield microscopy based on principal component analysis,”  Opt. Express, vol. 18, no. 7, pp. 6461–6476, Mar 2010.
  * [17] W. T. Freeman and E. H. Adelson,  “The design and use of steerable filters,”  IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 13, no. 9, pp. 891–906, 1991.
  * [18] J. G. Nagy and D. P. O’Leary,  “Restoring images degraded by spatially variant blur,”  SIAM Journal on Scientific Computing, vol. 19, no. 4, pp. 1063–1082, 1998.
  * [19] T. Popkin, A. Cavallaro, and D. Hands,  “Accurate and efficient method for smoothly space-variant gaussian blurring,”  IEEE Transactions on Image Processing, vol. 19, no. 5, pp. 1362–1370, 2010.
  * [20] A. V. Oppenheim, A. S. Willsky, and S. H. Nawab,  Signals & Systems (2nd Ed.),  Prentice-Hall, Inc., USA, 1996.
  * [21] Z. Mao, N. Chimitt, and S. H. Chan,  “Accelerating atmospheric turbulence simulation via learned phase-to-space transform,”  in Proc. IEEE/CVF International Conference on Computer Vision (ICCV), October 2021, pp. 14759–14768.
  * [22] J. W. Goodman,  Introduction to Fourier Optics,  Roberts and Company, Englewood, Colorado, 3 edition, 2005.
  * [23] N. Chimitt and S. H. Chan,  “Simulating anisoplanatic turbulence by sampling intermodal and spatially correlated Zernike coefficients,”  Optical Engineering, vol. 59, no. 8, pp. 1 – 26, 2020.
  * [24] N. Chimitt, X. Zhang, Z. Mao, and S. H. Chan,  “Real-time dense field phase-to-space simulation of imaging through atmospheric turbulence,”  IEEE Transactions on Computational Imaging, vol. 8, pp. 1159–1169, 2022.
  * [25] “Bridging the gap between computational photography and visual recognition: $5th$ $\text{UG2}+$ prize challenge,” http://cvpr2022.ug2challenge.org/dataset22_t3.html,  Track 3.
  * [26] D. Jin, Y. Chen, Y. Lu, J. Chen, P. Wang, Z. Liu, S. Guo, and X. Bai,  “Neutralizing the impact of atmospheric turbulence on complex scene imaging via deep learning,”  Nature Machine Intelligence, vol. 3, pp. 876–884, 2021.
  * [27] Z. Zhong, Y. Gao, Y. Zheng, B. Zheng, and I. Sato,  “Real-world video deblurring: A benchmark dataset and an efficient recurrent neural network,”  International Journal of Computer Vision, pp. 1–18, 2022.
  * [28] X. Zhang, Z. Mao, N. Chimitt, and S. H. Chan,  “Imaging through the atmosphere using turbulence mitigation transformer,”  Tech. Rep. 2207.06465, arXiv preprint, Jul 2022.