# A Motion Detection Algorithm Using Local Phase Information

Authors: Aurel A. Lazar; Nikul H. Ukani; Yiyin Zhou
Journal: Computational Intelligence and Neuroscience
 (2016)

---

## Abstract

Previous research demonstrated thatglobal phase alone can be used to faithfully represent visual scenes. Here we provide a reconstruction algorithm by using onlylocal phase information. We also demonstrate that local phase alone can be effectively used to detect local motion. The local phase-based motion detector is akin to models employed to detect motion in biological vision, for example, the Reichardt detector. The local phase-based motion detection algorithm introduced here consists of two building blocks. The first building block measures/evaluates the temporal change of the local phase. The temporal derivative of the local phase is shown to exhibit the structure of a second order Volterra kernel with two normalized inputs. We provide an efficient, FFT-based algorithm for implementing the change of the local phase. The second processing building block implements the detector; it compares the maximum of the Radon transform of the local phase derivative with a chosen threshold. We demonstrate examples of applying the local phase-based motion detection algorithm on several video sequences. We also show how the locally detected motion can be used for segmenting moving objects in video scenes and compare our local phase-based algorithm to segmentation achieved with a widely used optic flow algorithm.


---

## Body

## 1. Introduction

Following Marr, the design of an information processing system can be approached on multiple levels [1]. Figure 1 illustrates two levels of abstraction, namely, the algorithmic level and the physical circuit level. On the algorithmic level, one studies procedurally how the information is processed independently of the physical realization. The circuit level concerns the actual realization of the algorithm in physical hardware, for example, a biological neural circuit or silicon circuits in a digital signal processor. In this paper, we put forth a simple motion detection algorithm that is inspired by motion detection models of biological visual systems (in vivo neural circuit) and provide an efficient realization that can easily be implemented on commodity (in silico) DSP chips.Figure 1
Algorithms can have different physical realizations. For example, an algorithm can be implemented by neural circuits in a biological system. Alternately, it can be implemented on a digital signal processor. Algorithms direct the implementation on the physical layer. Conversely, biological neural circuits inspire new algorithmic designs, which can, in turn, be expressed and improved upon by a realizationin silico.Visual motion detection is critical to the survival of animals. Many biological visual systems have evolved highly efficient/effective neural circuits to detect visual motion. Motion detection is performed in parallel with other visual coding circuits and starts already in the early stages of visual processing. In the retina of vertebrates, it is known that at least three types of Direction-Selective Ganglion Cells (DSGC) are responsible for signaling visual motion at this early stage [2]. In flies, direction-selective neurons are found in the optic lobe, 3 synapses away from the photoreceptors [3].The small number of synapses between photoreceptors and direction-selective neurons suggests that the processing involved in motion detection is not highly complex but still very effective. In addition, the biological motion detection circuits are organized in a highly parallel way to enable fast, concurrent computation of motion. It is also interesting to note that the early stages of motion detection are carried out largely in the absence of spiking neurons, indicating that initial stages of motion detection are preferably performed in the “analog” domain. Taking advantage of continuous time processing may be critical for quickly processing motion since motion intrinsically elicits fast and large changes in the intensity levels, that is, large amounts of data under stringent time constraints.Modern, computer-based motion detection algorithms often employ optic flow techniques to estimate spatial changes in consecutive image frames [4, 5]. Although, often time, optic flow estimation algorithms produce accurate results, the computational demand to perform many of these algorithms is too high for real-time implementation.Several models for biological motion detection are available and their architecture is quite simple [6]. The Reichardt motion detector [7] was thought to be the underlying model for motion detection in insects [8]. The model is based on a correlation method to extract motion induced by spatiotemporal information patterns of light intensity. Therefore, it relies on a correlation/multiplication operation. A second model is the motion energy detector [9]. It uses spatiotemporal separable filters and a squaring nonlinearity to compute motion energy and it was shown to be equivalent to the Reichardt motion detector. Earlier work in the rabbit retina was the foundation to the Barlow-Levick model [10] of motion detection. The model relies on inhibition to compensate motion in the null direction.In this paper, we provide an alternative motion detection algorithm based on local phase information of the visual scene. Similar to mechanisms in other biological models, it operates in continuous time and in parallel. Moreover, the motion detection algorithm we propose can be efficiently implemented on parallel hardware. This is, again, similar to the properties of biological motion detection systems. Rather than focusing on velocity of motion, we focus on localization, that is, where the motion occurs in the visual field as well as the direction of motion.It has been shown that images can be represented by their global phase alone [11]. Here we provide a reconstruction algorithm of visual scenes by only usinglocal phase information, thereby demonstrating the spectrum of the representational capability of phase information.The Fourier shift property clearly suggests the relationship between the global shift of an image and the global phase shift in the frequency domain. We elevate this relationship by computing the change of local phase to indicate motion that appears locally in the visual scene. The local phases are computed using window functions that tile the visual field with overlapping segments, making it amenable for a highly parallel implementation. In addition, we propose a Radon-transform-based motion detection index on the change of local phases for the readout of the relation between local phases and motion.Interestingly, phase information has been largely ignored in the field of linear signal processing and for good reason. Phase-based processing is intrinsically non-linear. Recent researches, however, showed that phase information can be smartly employed in speech processing [12] and visual processing [13]. For example, spatial phase in an image is indicative of local features such as edges when considering phase congruency [14]. The role of spatial phase in computational and biological vision, emergence of visual illusions, and pattern recognition is discussed in [15]. Together with our result for motion detection, these studies suggest that phase information has a great potential for achieving efficient visual signal processing.This paper is organized as follows. In Section2, we show that local phase information can be used to faithfully represent visual information in the absence of amplitude information. In Section 3, we develop a simple local phase-based algorithm to detect motion in visual scenes and provide an efficient way for its implementation. We then provide examples and applications of the proposed motion detection algorithm to motion segmentation. We also compare our results to those obtained using motion detection algorithms based on optic flow techniques in Section 4. Finally, we summarize our results in Section 5.


## 2. Representation of Visual Scenes Using Phase Information

The use of complex valued transforms is widespread, for both representing and processing images. When represented in polar coordinates, the output of a complex valued transform of a signal can be split into amplitude and phase. In this section, we define two types of phases of an image: global phase and local phase. We then argue that both types of phases can faithfully represent an image and we provide a reconstruction algorithm that recovers the image from local phase information. This indicates that phase information alone can largely represent an image or video signal.It will become clear in the sections that follow that the use of phase for representing of image and video information leads to efficient ways to implement certain types of image/video processing algorithms, for example, motion detection algorithms.


### 2.1. The Global Phase of Images

Recall that the Fourier transform of a real valued imageu
=
u
(
x
,
y
), (
x
,
y
)
∈
R
2, is given by(1)
U
^
ω
x
,
ω
y
=
∫
R
2
u
x
,
y
e
-
j
ω
x
x
+
ω
y
y
d
x
d
y with U
^
(
ω
x
,
ω
y
)
∈
C and (
ω
x
,
ω
y
)
∈
R
2.In polar coordinates, the Fourier transform ofu can be expressed as(2)
U
^
ω
x
,
ω
y
=
A
^
ω
x
,
ω
y
e
j
ϕ
^
ω
x
,
ω
y
,where A
^
(
ω
x
,
ω
y
)
∈
R is the amplitude and ϕ
^
(
ω
x
,
ω
y
) is the phase of the Fourier transform of u.Definition 1.
The amplitude of the Fourier transform of an imageu
=
u
(
x
,
y
), (
x
,
y
)
∈
R
2, is called the global amplitude of u. The phase of the Fourier transform of an image u
=
u
(
x
,
y
), (
x
,
y
)
∈
R
2, is called the global phase of u.It is known that the global phase of an image plays an important role in the representation of natural images [11]. A classic example is to take two images and to exchange their global phases before their reconstruction using the inverse Fourier transform. The resulting images are slightly smeared but largely reflect the information contained in the global phase.


### 2.2. The Local Phase of Images

The global phase indicates the offset of sinusoids of different frequencies contained in the entire image. However, it is not intuitive to relate the global phase to local image features such as edges and their position in an image. To study these local features, it is necessary to modify the Fourier transform such that it reflects properties of a restricted region of an image. It is natural to consider the Short-Time (-Space) Fourier Transform (STFT):(3)
U
ω
x
,
ω
y
,
x
0
,
y
0
=
∫
R
2
u
x
,
y
w
x
-
x
0
,
y
-
y
0
e
-
j
ω
x
x
-
x
0
+
ω
y
y
-
y
0
d
x
d
y
,where w
=
w
(
x
,
y
), (
x
,
y
)
∈
R
2, is a real valued window function centered at (
x
0
,
y
0
)
∈
R
2. Typical choices of window functions include the Hann window and the Gaussian window. The (effectively) finite support of the window restricts the Fourier transform to local image analysis.Similarly to the Fourier transform, the STFT can be expressed in polar coordinates as(4)
U
ω
x
,
ω
y
,
x
0
,
y
0
=
A
ω
x
,
ω
y
,
x
0
,
y
0
e
j
ϕ
ω
x
,
ω
y
,
x
0
,
y
0
,where A
(
ω
x
,
ω
y
,
x
0
,
y
0
)
∈
R is the amplitude and ϕ
(
ω
x
,
ω
y
,
x
0
,
y
0
) is the phase of the STFT.Definition 2.
The amplitude of the STFT of an imageu
=
u
(
x
,
y
), (
x
,
y
)
∈
R
2, is called the local amplitude of u. The phase of the STFT of an image u
=
u
(
x
,
y
), (
x
,
y
)
∈
R
2, is called the local phase of u.Note that whenw is a Gaussian window, the STFT of u evaluated at (
ω
x
,
ω
y
,
x
0
,
y
0
) can be equivalently viewed as the response of a complex-valued Gabor receptive field (5)
h
x
,
y
=
e
-
x
-
x
0
2
+
y
-
y
0
2
/
2
σ
2
e
-
j
ω
x
x
-
x
0
+
ω
y
y
-
y
0to u.In this case, the window of the STFT clearly is given by(6)
w
x
,
y
=
e
-
x
2
+
y
2
/
2
σ
2
.Therefore, the STFT can be realized by an ensemble of Gabor receptive fields that are common in modeling simple and complex cells (neurons) in the primary visual cortex [16].


### 2.3. Reconstruction of Images from Local Phase

Amplitude and phase can be interpreted as measurements/projections of images that are indicative of their information content. Classically, when both the global amplitude and phase are known, it is straightforward to reconstruct the image. The reconstruction calls for computing the inverse Fourier transform given the global amplitude and phase. Similarly, when using local amplitude and phase, if the sampling functions form a basis or frame in a space of images, the reconstruction is provided by the formalism of wavelet theory, for example, using Gabor wavelets [17].Amplitude or phase represents partial information extracted from visual scenes. They are obtained via nonlinear sampling, that is, a nonlinear operation for extracting the amplitude and phase information from images. The nonlinear operation makes reconstruction from either amplitude or phase alone difficult. Earlier studies and recent development in solving quadratic constraints, however, suggest that it is possible to reconstruct images from global or local amplitude information [18, 19].While computing the amplitude requires a second order (quadratic) nonlinearity, computing the phase calls for higher order nonlinear operators (Volterra kernels). It is possible, however, to reconstruct up to a constant scale an image from its global phase information alone without explicitly using the amplitude information. An algorithm was provided for solving this problem in the discrete signal processing domain in [11]. The algorithm smartly avoids using the inverse tangent function by reformulating the phase measurements as a set of linear equations that are easy to solve.Using a similar argument, we demonstrate in the following that, up to a constant scale, a bandlimited signalu
=
u
(
x
,
y
), (
x
,
y
)
∈
R
2, can be reconstructed from its local phase alone. We first formulate the encoding of an image u by local phase, using the Gabor receptive fields as a special case. It is straightforward then to formulate the problem with local phase computed from other types of STFTs.Formally, we consider an image,u
=
u
(
x
,
y
), on the domain R
2, to be an element of a space of trigonometric polynomials H of the form(7)
u
x
,
y
=
∑
l
x
=
-
L
x
L
x
∑
l
y
=
-
L
y
L
y
c
l
x
l
y
e
l
x
l
y
x
,
y
,where(8)
e
l
x
l
y
=
exp
⁡
j
l
x
Ω
x
L
x
x
exp
⁡
j
l
y
Ω
y
L
y
y
,
l
x
=
-
L
x
,
…
,
L
x
,
l
y
=
-
L
y
,
…
,
L
yare the set of basis functions of H, Ω
x, L
x are the bandwidth and the order, respectively, of the space in the x dimension, and Ω
y, L
y are the bandwidth and the order, respectively, of the space in the y dimension.H is a Reproducing Kernel Hilbert Space with inner product [20, 21](9)
u
1
,
u
2
=
∫
0
,
T
x
×
0
,
T
y
u
1
x
,
y
u
2
x
,
y
¯
d
x
d
y
,where T
x
=
2
π
L
x
/
Ω
x, T
y
=
2
π
L
y
/
Ω
y are the period of the space in the x and y dimensions, respectively.Consider a bank ofN Gabor receptive fields(10)
h
k
l
,
m
n
x
,
y
=
T
k
l
w
x
,
y
e
-
j
ω
x
m
x
+
ω
y
n
y
,where T
k
l is the translation operator with T
k
l
u
=
u
(
x
-
k
b
0
,
y
-
l
b
0
), k
,
l
∈
Z, b
0
>
0, 0
≤
k
b
0
≤
T
x, and 0
≤
l
b
0
≤
T
y, and ω
x
m
=
m
ω
0, ω
y
n
=
n
ω
0, m
,
n
∈
Z, ω
0
>
0, -
Ω
x
≤
m
ω
0
≤
Ω
x, and -
Ω
y
≤
n
ω
0
≤
Ω
y. The responses of the Gabor receptive fields to the input u
(
x
,
y
) are given by(11)
∫
R
2
u
x
,
y
h
k
l
,
m
n
x
,
y
d
x
d
y
=
∫
R
2
u
x
,
y
w
x
-
k
b
0
,
y
-
l
b
0
e
-
j
ω
x
m
x
-
k
b
0
+
ω
y
n
y
-
l
b
0
d
x
d
y
=
A
k
l
,
m
n
e
j
ϕ
k
l
,
m
n
,where A
k
l
,
m
n
≥
0, A
k
l
,
m
n
∈
R, is the local amplitude and ϕ
k
l
,
m
n
∈
[
0,2
π
) is the local phase.Dividing both sides of (11) by e
j
ϕ
k
l
,
m
n, we have(12)
∫
R
2
u
x
,
y
w
x
-
k
b
0
,
y
-
l
b
0
e
-
j
ω
x
m
x
-
k
b
0
+
ω
y
n
y
-
l
b
0
+
ϕ
k
l
,
m
n
d
x
d
y
=
A
k
l
,
m
n
.Since A
k
l
,
m
n
∈
R, we have(13)
∫
R
2
u
x
,
y
w
x
-
k
b
0
,
y
-
l
b
0
cos
⁡
ω
x
m
x
-
k
b
0
+
ω
y
n
y
-
l
b
0
+
ϕ
k
l
,
m
n
d
x
d
y
=
A
k
l
,
m
n
,
(14)
∫
R
2
u
x
,
y
w
x
-
k
b
0
,
y
-
l
b
0
sin
⁡
ω
x
m
x
-
k
b
0
+
ω
y
n
y
-
l
b
0
+
ϕ
k
l
,
m
n
d
x
d
y
=
0
.Note that w
x
-
k
b
0
,
y
-
l
b
0
sin
⁡
(
ω
x
m
(
x
-
k
b
0
)
+
ω
y
n
(
y
-
l
b
0
)
+
ϕ
k
l
,
m
n
) is a real-valued Gabor receptive field with a preferred phase at ϕ
k
l
,
m
n
+
π
/
2
-
ω
x
m
k
b
0
-
ω
y
n
l
b
0.Remark 3.
Assuming that the local phase informationϕ
k
l
,
m
n is obtained via measurements, that is, filtering the image u with pairs of Gabor receptive fields (10), the set of linear equations (14) has a simple interpretation: the image u is orthogonal to the space spanned by the functions(15)
w
x
-
k
b
0
,
y
-
l
b
0
sin
⁡
ω
x
m
x
-
k
b
0
+
ω
y
n
y
-
l
b
0
+
ϕ
k
l
,
m
n
,where (
k
,
l
,
m
,
n
)
∈
I with I
=
{
(
k
,
l
,
m
,
n
)
∈
Z
4
∣
0
≤
k
b
0
≤
T
x
,
0
≤
l
b
0
≤
T
y
,
-
Ω
x
≤
m
ω
0
≤
Ω
x, -
Ω
y
≤
n
ω
0
≤
Ω
y
}.We are now in the position to provide a reconstruction algorithm of the image from phaseϕ
k
l
,
m
n, (
k
,
l
,
m
,
n
)
∈
I.Lemma 4.
u can be reconstructed from ϕ
k
l
,
m
n, (
k
,
l
,
m
,
n
)
∈
I as(16)
u
x
,
y
=
∑
l
x
=
-
L
x
L
x
∑
l
y
=
-
L
y
L
y
c
l
x
l
y
e
l
x
l
y
x
,
y
,
with(17)
Φ
c
=
0
,where Φ is a matrix whose pth row and qth column entry are(18)
Φ
p
q
=
∫
R
2
w
x
-
k
b
0
,
y
-
l
b
0
sin
⁡
ω
x
m
x
-
k
b
0
+
ω
y
n
y
-
l
b
0
+
ϕ
k
l
,
m
n
e
l
x
l
y
x
,
y
d
x
d
y
.Here p traverses the set I, and q
=
(
2
L
y
+
1
)
(
l
x
+
L
x
)
+
(
l
y
+
L
y
+
1
). c is a vector of the form(19)
c
=
c
-
L
x
,
-
L
y
,
c
-
L
x
,
-
L
y
+
1
,
…
,
c
-
L
x
,
L
y
,
c
-
L
x
+
1
,
-
L
y
,
c
-
L
x
+
1
,
-
L
y
+
1
,
…
,
c
-
L
x
+
1
,
L
y
,
…
,
c
L
x
,
-
L
y
,
c
L
x
,
-
L
y
+
1
,
…
,
c
L
x
,
L
y
Tthat belongs to the null space of Φ. A necessary condition for perfect reconstruction of u, up to a constant scale, is that N
≥
(
2
L
x
+
1
)
(
2
L
y
+
1
)
-
1, where N is the number of phase measurements.Proof.
Substituting (7) into (14), we obtain(20)
∑
l
x
=
-
L
x
L
x
∑
l
y
=
-
L
y
L
y
c
l
x
l
y
∫
R
2
w
x
-
k
b
0
,
y
-
l
b
0
sin
⁡
ω
x
m
x
-
k
b
0
+
ω
y
n
y
-
l
b
0
+
ϕ
k
l
,
m
n
e
l
x
l
y
x
,
y
d
x
d
y
=
0
, for all k
,
l
,
m
,
n
∈
Z. Therefore, we have(21)
Φ
c
=
0
, inferring that c is in the null space of Φ.
IfN
<
(
2
L
x
+
1
)
(
2
L
y
+
1
)
-
1, it follows from the rank-nullity theorem that dim
⁡
(
n
u
l
l
(
Φ
)
)
>
1 leading to multiple linearly independent solutions to Φ
c
=
0.Example 5.
In Figure2 an example of reconstruction of an image is shown using only local phase information. The reconstructed signal was scaled to match the original signal. The SNR of the reconstruction is 44.48 [dB]. An alternative way to obtain a unique reconstruction is to include an additional measurement, for example, the mean value of the signal ∫
R
2
‍
u
x
,
y
d
x
d
y to the system of linear equations (14).Figure 2
An example of reconstruction of image from local phase information. (a) Original image, the dimension of the space is 16,129. (b) Reconstructed image from 20,102 phase measurements, scaled to match the original. (c) Error.
(a)
(b)
(c)


## 2.1. The Global Phase of Images

Recall that the Fourier transform of a real valued imageu
=
u
(
x
,
y
), (
x
,
y
)
∈
R
2, is given by(1)
U
^
ω
x
,
ω
y
=
∫
R
2
u
x
,
y
e
-
j
ω
x
x
+
ω
y
y
d
x
d
y with U
^
(
ω
x
,
ω
y
)
∈
C and (
ω
x
,
ω
y
)
∈
R
2.In polar coordinates, the Fourier transform ofu can be expressed as(2)
U
^
ω
x
,
ω
y
=
A
^
ω
x
,
ω
y
e
j
ϕ
^
ω
x
,
ω
y
,where A
^
(
ω
x
,
ω
y
)
∈
R is the amplitude and ϕ
^
(
ω
x
,
ω
y
) is the phase of the Fourier transform of u.Definition 1.
The amplitude of the Fourier transform of an imageu
=
u
(
x
,
y
), (
x
,
y
)
∈
R
2, is called the global amplitude of u. The phase of the Fourier transform of an image u
=
u
(
x
,
y
), (
x
,
y
)
∈
R
2, is called the global phase of u.It is known that the global phase of an image plays an important role in the representation of natural images [11]. A classic example is to take two images and to exchange their global phases before their reconstruction using the inverse Fourier transform. The resulting images are slightly smeared but largely reflect the information contained in the global phase.


## 2.2. The Local Phase of Images

The global phase indicates the offset of sinusoids of different frequencies contained in the entire image. However, it is not intuitive to relate the global phase to local image features such as edges and their position in an image. To study these local features, it is necessary to modify the Fourier transform such that it reflects properties of a restricted region of an image. It is natural to consider the Short-Time (-Space) Fourier Transform (STFT):(3)
U
ω
x
,
ω
y
,
x
0
,
y
0
=
∫
R
2
u
x
,
y
w
x
-
x
0
,
y
-
y
0
e
-
j
ω
x
x
-
x
0
+
ω
y
y
-
y
0
d
x
d
y
,where w
=
w
(
x
,
y
), (
x
,
y
)
∈
R
2, is a real valued window function centered at (
x
0
,
y
0
)
∈
R
2. Typical choices of window functions include the Hann window and the Gaussian window. The (effectively) finite support of the window restricts the Fourier transform to local image analysis.Similarly to the Fourier transform, the STFT can be expressed in polar coordinates as(4)
U
ω
x
,
ω
y
,
x
0
,
y
0
=
A
ω
x
,
ω
y
,
x
0
,
y
0
e
j
ϕ
ω
x
,
ω
y
,
x
0
,
y
0
,where A
(
ω
x
,
ω
y
,
x
0
,
y
0
)
∈
R is the amplitude and ϕ
(
ω
x
,
ω
y
,
x
0
,
y
0
) is the phase of the STFT.Definition 2.
The amplitude of the STFT of an imageu
=
u
(
x
,
y
), (
x
,
y
)
∈
R
2, is called the local amplitude of u. The phase of the STFT of an image u
=
u
(
x
,
y
), (
x
,
y
)
∈
R
2, is called the local phase of u.Note that whenw is a Gaussian window, the STFT of u evaluated at (
ω
x
,
ω
y
,
x
0
,
y
0
) can be equivalently viewed as the response of a complex-valued Gabor receptive field (5)
h
x
,
y
=
e
-
x
-
x
0
2
+
y
-
y
0
2
/
2
σ
2
e
-
j
ω
x
x
-
x
0
+
ω
y
y
-
y
0to u.In this case, the window of the STFT clearly is given by(6)
w
x
,
y
=
e
-
x
2
+
y
2
/
2
σ
2
.Therefore, the STFT can be realized by an ensemble of Gabor receptive fields that are common in modeling simple and complex cells (neurons) in the primary visual cortex [16].


## 2.3. Reconstruction of Images from Local Phase

Amplitude and phase can be interpreted as measurements/projections of images that are indicative of their information content. Classically, when both the global amplitude and phase are known, it is straightforward to reconstruct the image. The reconstruction calls for computing the inverse Fourier transform given the global amplitude and phase. Similarly, when using local amplitude and phase, if the sampling functions form a basis or frame in a space of images, the reconstruction is provided by the formalism of wavelet theory, for example, using Gabor wavelets [17].Amplitude or phase represents partial information extracted from visual scenes. They are obtained via nonlinear sampling, that is, a nonlinear operation for extracting the amplitude and phase information from images. The nonlinear operation makes reconstruction from either amplitude or phase alone difficult. Earlier studies and recent development in solving quadratic constraints, however, suggest that it is possible to reconstruct images from global or local amplitude information [18, 19].While computing the amplitude requires a second order (quadratic) nonlinearity, computing the phase calls for higher order nonlinear operators (Volterra kernels). It is possible, however, to reconstruct up to a constant scale an image from its global phase information alone without explicitly using the amplitude information. An algorithm was provided for solving this problem in the discrete signal processing domain in [11]. The algorithm smartly avoids using the inverse tangent function by reformulating the phase measurements as a set of linear equations that are easy to solve.Using a similar argument, we demonstrate in the following that, up to a constant scale, a bandlimited signalu
=
u
(
x
,
y
), (
x
,
y
)
∈
R
2, can be reconstructed from its local phase alone. We first formulate the encoding of an image u by local phase, using the Gabor receptive fields as a special case. It is straightforward then to formulate the problem with local phase computed from other types of STFTs.Formally, we consider an image,u
=
u
(
x
,
y
), on the domain R
2, to be an element of a space of trigonometric polynomials H of the form(7)
u
x
,
y
=
∑
l
x
=
-
L
x
L
x
∑
l
y
=
-
L
y
L
y
c
l
x
l
y
e
l
x
l
y
x
,
y
,where(8)
e
l
x
l
y
=
exp
⁡
j
l
x
Ω
x
L
x
x
exp
⁡
j
l
y
Ω
y
L
y
y
,
l
x
=
-
L
x
,
…
,
L
x
,
l
y
=
-
L
y
,
…
,
L
yare the set of basis functions of H, Ω
x, L
x are the bandwidth and the order, respectively, of the space in the x dimension, and Ω
y, L
y are the bandwidth and the order, respectively, of the space in the y dimension.H is a Reproducing Kernel Hilbert Space with inner product [20, 21](9)
u
1
,
u
2
=
∫
0
,
T
x
×
0
,
T
y
u
1
x
,
y
u
2
x
,
y
¯
d
x
d
y
,where T
x
=
2
π
L
x
/
Ω
x, T
y
=
2
π
L
y
/
Ω
y are the period of the space in the x and y dimensions, respectively.Consider a bank ofN Gabor receptive fields(10)
h
k
l
,
m
n
x
,
y
=
T
k
l
w
x
,
y
e
-
j
ω
x
m
x
+
ω
y
n
y
,where T
k
l is the translation operator with T
k
l
u
=
u
(
x
-
k
b
0
,
y
-
l
b
0
), k
,
l
∈
Z, b
0
>
0, 0
≤
k
b
0
≤
T
x, and 0
≤
l
b
0
≤
T
y, and ω
x
m
=
m
ω
0, ω
y
n
=
n
ω
0, m
,
n
∈
Z, ω
0
>
0, -
Ω
x
≤
m
ω
0
≤
Ω
x, and -
Ω
y
≤
n
ω
0
≤
Ω
y. The responses of the Gabor receptive fields to the input u
(
x
,
y
) are given by(11)
∫
R
2
u
x
,
y
h
k
l
,
m
n
x
,
y
d
x
d
y
=
∫
R
2
u
x
,
y
w
x
-
k
b
0
,
y
-
l
b
0
e
-
j
ω
x
m
x
-
k
b
0
+
ω
y
n
y
-
l
b
0
d
x
d
y
=
A
k
l
,
m
n
e
j
ϕ
k
l
,
m
n
,where A
k
l
,
m
n
≥
0, A
k
l
,
m
n
∈
R, is the local amplitude and ϕ
k
l
,
m
n
∈
[
0,2
π
) is the local phase.Dividing both sides of (11) by e
j
ϕ
k
l
,
m
n, we have(12)
∫
R
2
u
x
,
y
w
x
-
k
b
0
,
y
-
l
b
0
e
-
j
ω
x
m
x
-
k
b
0
+
ω
y
n
y
-
l
b
0
+
ϕ
k
l
,
m
n
d
x
d
y
=
A
k
l
,
m
n
.Since A
k
l
,
m
n
∈
R, we have(13)
∫
R
2
u
x
,
y
w
x
-
k
b
0
,
y
-
l
b
0
cos
⁡
ω
x
m
x
-
k
b
0
+
ω
y
n
y
-
l
b
0
+
ϕ
k
l
,
m
n
d
x
d
y
=
A
k
l
,
m
n
,
(14)
∫
R
2
u
x
,
y
w
x
-
k
b
0
,
y
-
l
b
0
sin
⁡
ω
x
m
x
-
k
b
0
+
ω
y
n
y
-
l
b
0
+
ϕ
k
l
,
m
n
d
x
d
y
=
0
.Note that w
x
-
k
b
0
,
y
-
l
b
0
sin
⁡
(
ω
x
m
(
x
-
k
b
0
)
+
ω
y
n
(
y
-
l
b
0
)
+
ϕ
k
l
,
m
n
) is a real-valued Gabor receptive field with a preferred phase at ϕ
k
l
,
m
n
+
π
/
2
-
ω
x
m
k
b
0
-
ω
y
n
l
b
0.Remark 3.
Assuming that the local phase informationϕ
k
l
,
m
n is obtained via measurements, that is, filtering the image u with pairs of Gabor receptive fields (10), the set of linear equations (14) has a simple interpretation: the image u is orthogonal to the space spanned by the functions(15)
w
x
-
k
b
0
,
y
-
l
b
0
sin
⁡
ω
x
m
x
-
k
b
0
+
ω
y
n
y
-
l
b
0
+
ϕ
k
l
,
m
n
,where (
k
,
l
,
m
,
n
)
∈
I with I
=
{
(
k
,
l
,
m
,
n
)
∈
Z
4
∣
0
≤
k
b
0
≤
T
x
,
0
≤
l
b
0
≤
T
y
,
-
Ω
x
≤
m
ω
0
≤
Ω
x, -
Ω
y
≤
n
ω
0
≤
Ω
y
}.We are now in the position to provide a reconstruction algorithm of the image from phaseϕ
k
l
,
m
n, (
k
,
l
,
m
,
n
)
∈
I.Lemma 4.
u can be reconstructed from ϕ
k
l
,
m
n, (
k
,
l
,
m
,
n
)
∈
I as(16)
u
x
,
y
=
∑
l
x
=
-
L
x
L
x
∑
l
y
=
-
L
y
L
y
c
l
x
l
y
e
l
x
l
y
x
,
y
,
with(17)
Φ
c
=
0
,where Φ is a matrix whose pth row and qth column entry are(18)
Φ
p
q
=
∫
R
2
w
x
-
k
b
0
,
y
-
l
b
0
sin
⁡
ω
x
m
x
-
k
b
0
+
ω
y
n
y
-
l
b
0
+
ϕ
k
l
,
m
n
e
l
x
l
y
x
,
y
d
x
d
y
.Here p traverses the set I, and q
=
(
2
L
y
+
1
)
(
l
x
+
L
x
)
+
(
l
y
+
L
y
+
1
). c is a vector of the form(19)
c
=
c
-
L
x
,
-
L
y
,
c
-
L
x
,
-
L
y
+
1
,
…
,
c
-
L
x
,
L
y
,
c
-
L
x
+
1
,
-
L
y
,
c
-
L
x
+
1
,
-
L
y
+
1
,
…
,
c
-
L
x
+
1
,
L
y
,
…
,
c
L
x
,
-
L
y
,
c
L
x
,
-
L
y
+
1
,
…
,
c
L
x
,
L
y
Tthat belongs to the null space of Φ. A necessary condition for perfect reconstruction of u, up to a constant scale, is that N
≥
(
2
L
x
+
1
)
(
2
L
y
+
1
)
-
1, where N is the number of phase measurements.Proof.
Substituting (7) into (14), we obtain(20)
∑
l
x
=
-
L
x
L
x
∑
l
y
=
-
L
y
L
y
c
l
x
l
y
∫
R
2
w
x
-
k
b
0
,
y
-
l
b
0
sin
⁡
ω
x
m
x
-
k
b
0
+
ω
y
n
y
-
l
b
0
+
ϕ
k
l
,
m
n
e
l
x
l
y
x
,
y
d
x
d
y
=
0
, for all k
,
l
,
m
,
n
∈
Z. Therefore, we have(21)
Φ
c
=
0
, inferring that c is in the null space of Φ.
IfN
<
(
2
L
x
+
1
)
(
2
L
y
+
1
)
-
1, it follows from the rank-nullity theorem that dim
⁡
(
n
u
l
l
(
Φ
)
)
>
1 leading to multiple linearly independent solutions to Φ
c
=
0.Example 5.
In Figure2 an example of reconstruction of an image is shown using only local phase information. The reconstructed signal was scaled to match the original signal. The SNR of the reconstruction is 44.48 [dB]. An alternative way to obtain a unique reconstruction is to include an additional measurement, for example, the mean value of the signal ∫
R
2
‍
u
x
,
y
d
x
d
y to the system of linear equations (14).Figure 2
An example of reconstruction of image from local phase information. (a) Original image, the dimension of the space is 16,129. (b) Reconstructed image from 20,102 phase measurements, scaled to match the original. (c) Error.
(a)
(b)
(c)


## 3. Visual Motion Detection from Phase Information

In this section we consider visual fields that change as a function of time. For notational simplicityu will denote here the space-time intensity of the visual field.


### 3.1. The Global Phase Equation for Translational Motion

Letu
=
u
(
x
,
y
,
t
), (
x
,
y
)
∈
R
2, t
∈
R, be a visual stimulus. If the visual stimulus is a pure translation of the signal at u
(
x
,
y
,
0
), that is,(22)
u
x
,
y
,
t
=
u
x
-
s
x
t
,
y
-
s
y
t
,
0
,where(23)
s
x
t
=
∫
0
t
v
x
s
d
s
,
s
y
t
=
∫
0
t
v
y
s
d
sare the total length of translation at time t in each dimension and v
x
(
t
) and v
y
(
t
) are the corresponding instantaneous velocity components, then the only difference between u
(
x
,
y
,
t
) and u
(
x
,
y
,
0
) in the Fourier domain is captured by their global phase. More formally, consider the following.Lemma 6.
The change (derivative) of the global phase is given by(24)
d
ϕ
^
ω
x
,
ω
y
,
t
d
t
=
-
ω
x
v
x
t
-
ω
y
v
y
t
=
-
ω
x
,
ω
y
v
x
t
,
v
y
t
T
,where, by abuse of notation, ϕ
^
ω
x
,
ω
y
,
t denotes the global phase of u
x
,
y
,
t and ϕ
^
ω
x
,
ω
y
,
0 is the initial condition.Proof.
IfF
u
(
·
,
·
,
0
)  (
ω
x
,
ω
y
) is the 2D (spatial) Fourier transform of u
=
u
(
x
,
y
,
0
), (
x
,
y
)
∈
R
2, by the Fourier shift theorem, we have(25)
F
u
·
,
·
,
t
ω
x
,
ω
y
=
F
u
·
,
·
,
0
ω
x
,
ω
y
e
-
j
ω
x
s
x
t
+
ω
y
s
y
t
.For a certain frequency component (
ω
x
,
ω
y
), the change of its global phase over time amounts to(26)
d
ϕ
^
ω
x
,
ω
y
,
t
d
t
=
-
ω
x
d
s
x
t
d
t
-
ω
y
d
s
y
t
d
t
=
-
ω
x
v
x
t
-
ω
y
v
y
t
=
-
ω
x
,
ω
y
v
x
t
,
v
y
t
T
.Therefore, in the simple case where the entire visual field is shifting, the derivative of the phase of Fourier components indicates motion, and it can be obtained by the inner product between the component frequency and the velocity vector.


### 3.2. The Change of Local Phase

#### 3.2.1. The Local Phase Equation for Translational Motion

The analysis in Section3.1 applies toglobal motion. This type of motion occurs most frequently when the imaging device, either an eye or a camera, moves. Visual motion in the natural environment, however, is more diverse across the screen since it is, often time, produced by multiple moving objects. The objects can be small and the motion more localized in the visual field.Taking the global phase ofu will not simply reveal where motion of independent objects takes place or their direction/velocity of motion. The ease of interpretation of motion by using the Fourier transform, however, motivates us to reuse the same concept in detecting local motion. This can be achieved by restricting the domain of the visual field where the Fourier transform is applied.To be able to detect local motion, we consider the local phase ofu
(
x
,
y
,
t
) by taking the STFT with window function w
(
x
,
y
). Note that, the STFT and its ubiquitous implementation in DSP chips can be extensively used in any dimension. For simplicity and without loss of generality, we consider the window to be centered at (
0,0
). The STFT is given by(27)
∫
R
2
u
x
,
y
,
t
w
x
,
y
e
-
j
ω
x
x
+
ω
y
y
d
x
d
y
=
A
00
ω
x
,
ω
y
,
t
e
j
ϕ
00
ω
x
,
ω
y
,
t
,where, by abuse of notation, A
00
(
ω
x
,
ω
y
,
t
) is the amplitude and ϕ
00
(
ω
x
,
ω
y
,
t
) the local phase.Before we move on to the mathematical analysis, we can intuitively explain the relation between the change in local phase and visual motion taking place across the window support. First, if the stimulus undergoes a uniform change of intensity or it changes proportionally over time due to lighting conditions, for example, the local phase does not change since the phase is invariant with respect to intensity scaling. Therefore, the local phase does not change for such nonmotion stimuli. Second, a rigid edge moving across the window support will induce a phase change.For a strictly translational signal within the window support (footprint), for example,(28)
u
x
,
y
,
t
=
u
x
-
s
x
t
,
y
-
s
y
t
,
0
,
f
o
r
x
,
y
∈
supp
⁡
w
x
,
y
, where s
x
(
t
) and s
y
(
t
) are as defined in (23), we have the following result.Lemma 7.
Consider the following:(29)
d
ϕ
00
d
t
ω
x
,
ω
y
,
t
=
-
d
s
x
t
d
t
ω
x
-
d
s
y
t
d
t
ω
y
+
v
00
ω
x
,
ω
y
,
t
,where, by abuse of notation, ϕ
00
(
ω
x
,
ω
y
,
t
) is the local phase of u
(
x
,
y
,
t
) and ϕ
00
(
ω
x
,
ω
y
,
0
) is the initial condition. The functional form of the term v
00
=
v
00
(
ω
x
,
ω
y
,
t
) is provided in the Appendix.Proof.
The derivation of (29) and the functional form of v
00 are given in the Appendix.Remark 8.
We notice that the derivative of the local phase has similar structure to that of the global phase, but for the added termv
00. Through simulations, we observed that the first two terms in (29) dominate over the last term for an ON or OFF moving edge [3]. For example, Figure 3 shows the derivative of the local phase given in (29) for an ON edge moving with velocity (
40,0
) pixels/sec.Figure 3
Derivative of the local phase for an ON edge moving in the positivex direction.Remark 9.
Note thatϕ
00
(
ω
x
,
ω
y
,
t
) may not be differentiable even if u
(
x
,
y
,
t
) is differentiable, particularly when A
00
(
ω
x
,
ω
y
,
t
)
=
0. For example, the spatial phase can jump from a positive value to zero when A
00
(
ω
x
,
ω
y
,
t
) diminishes. This also suggests that the instantaneous local spatial phase is less informative about a region of a visual scene whenever A
00
(
ω
x
,
ω
y
,
t
) is close to zero. Nevertheless, the time derivative of the local phase can be approximated by applying a high-pass filter to ϕ
00
(
ω
x
,
ω
y
,
t
).


#### 3.2.2. The Block Structure for Computing the Local Phase

We construct Gaussian windows along thex, y dimensions. The Gaussian windows are defined as(30)
T
k
l
w
x
,
y
=
e
-
x
-
x
k
2
+
x
-
y
l
2
/
2
σ
2
,where x
k
=
k
b
0, y
l
=
l
b
0, in which b
0
∈
Z
+ is the distance between two neighboring windows and 1
≤
k
b
0
≤
P
x, 1
≤
l
b
0
≤
P
y, where P
x
,
P
y
∈
Z
+ are the number of pixels of the screen in x and y directions, respectively.We then take the 2D Fourier transform of the windowed video signalu
(
x
,
y
,
t
)
(
T
k
l
w
)
(
x
,
y
) and write in polar form(31)
∫
R
2
u
x
,
y
,
t
T
k
l
w
x
,
y
e
-
j
ω
x
x
-
x
k
+
ω
y
y
-
y
l
d
x
d
y
=
A
k
l
ω
x
,
ω
y
,
t
e
j
ϕ
k
l
ω
x
,
ω
y
,
t
.The above integral can be very efficiently evaluated using the 2D FFT in discrete domain defined on M
×
M blocks approximating the footprint of the Gaussian windows. For example, the standard deviation of the Gaussian windows we use in the examples in Section 4 is 4 pixels. A block of 32
×
32 pixels (M
=
32) is sufficient to cover the effective support (or footprint) of the Gaussian window. At the same time, the size of the block is a power of 2, which is most suitable for FFT-based computation. The processing of each block (
k
,
l
) is independent of all other blocks; thereby, parallelism is readily achieved.Note that the size of the window is informed by the size of the objects one is interested in locating. Measurements of the local phase using smaller window functions are less robust to noise. Larger windows would enhance object motion detection if the object size is comparable to the window size. However, there would be an increased likelihood of independent movement of multiple objects within the same window, which is not modeled here and thereby may not be robustly detected.Therefore, for each block(
k
,
l
), we obtain M
2 measurements of the phase ϕ
k
l
(
ω
x
m
,
ω
y
m
,
t
) at every time instant t, with (
ω
x
m
,
ω
y
m
)
∈
D
2, where(32)
D
2
=
ω
x
m
=
m
ω
0
,
ω
y
n
=
n
ω
0
,
m
,
n
=
-
M
/
2
,
-
M
/
2
+
1
,
…
,
M
/
2
-
1
, with ω
0
=
2
π
/
M. We then compute the temporal derivative of the phase, that is, d
ϕ
k
l
/
d
t
(
ω
x
,
ω
y
,
t
) for (
ω
x
,
ω
y
)
∈
D
2.We further illustrate an example of the block structure in Figure4. Figure 4(a) shows an example of an image of 64
×
64 pixels. Four Gaussian windows are shown each with a standard deviation of 4 pixels. The distance between the centers of two neighboring Gaussian windows is 6 pixels. The red solid square shows a 32
×
32-pixel block with k
=
3, l
=
3, which encloses effective support of the Gaussian window on top-left (k
=
0, l
=
0 is the block with Gaussian window centered at pixel (
1,1
)). The green dashed square shows another 32
×
32-pixel block with k
=
3, l
=
7. The two Gaussian windows on the right are associated with the blocks k
=
7, l
=
3 and k
=
7, l
=
7, respectively. Cross section of all Gaussian windows with k
=
7, l
∈
[
0,10
], that is, those centered on the magenta line, are shown in Figure 4(b). The red and blue curve in Figure 4(b) correspond to the two Gaussian windows shown in Figure 4(a). Figure 4(b) also suggests that some of the Gaussian windows are cut off on the boundaries. This is, however, equivalent to assuming that the pixel values outside the boundary are always zero, and it will not significantly affect motion detection based on the change of local phase.Figure 4
Example of block structure. (a) Four yellow disks show 4 Gaussian window functions with translation parameterk
=
3, l
=
3 (top-left), k
=
3, l
=
7 (bottom-left), k
=
7, l
=
3 (top-right), and k
=
7, l
=
7 (bottom-right). The red solid square shows a 32
×
32-pixel block approximating the Gaussian window with k
=
3, l
=
3, and the green dashed square shows a 32
×
32-pixel block approximating the Gaussian window with k
=
3, l
=
7. (b) Cross section of all 11 Gaussian windows with translation parameters k
=
7, l
∈
[
0,10
]. The cross section is taken as indicated by the magenta line in (a), and the red and blue curve correspond to cross sections of the two Gaussian windows shown in (a) centered on the magenta line.
(a)
(b)Since the phase, and thereby the phase change, is noisier when the local amplitude is low, an additional denoising step can be employed to discount the measurements ofd
ϕ
k
l
/
d
t
(
ω
x
,
ω
y
,
t
) for low amplitude values A
k
l
(
ω
x
,
ω
y
,
t
). The denoising is given by(33)
d
ϕ
k
l
d
t
ω
x
,
ω
y
,
t
A
k
l
ω
x
,
ω
y
,
t
1
/
M
2
∑
ω
x
,
ω
y
∈
D
2
A
k
l
ω
x
,
ω
y
,
t
+
ϵ
,where ϵ
>
0 is a constant, and (
ω
x
,
ω
y
)
∈
D
2.


### 3.3. The Phase-Based Detector

We propose here a block FFT based algorithm to detect motion using phase information. Such an algorithm is, due to its simplicity and parallelism, highly suitable for anin silico implementation.


#### 3.3.1. Radon Transform on the Change of Phases

We exploit the approximately linear structure of the phase derivative for blocks exhibiting motion by computing the Radon transform ofd
ϕ
k
l
/
d
t
(
ω
x
,
ω
y
,
t
) over a circular bounded domain C
=
{
(
ω
x
,
ω
y
)
∣
(
ω
x
,
ω
y
)
∈
D
2
,
ω
x
2
+
ω
y
2
<
π
2
}.The Radon transform of the change of phase in the domainC is given by(34)
R
d
ϕ
k
l
d
t
ρ
,
θ
,
t
=
∫
R
d
ϕ
k
l
d
t
ρ
·
cos
⁡
θ
-
s
·
sin
⁡
θ
,
ρ
·
sin
⁡
θ
+
s
·
cos
⁡
θ
,
t
1
C
ρ
·
cos
⁡
θ
-
s
·
sin
⁡
θ
,
ρ
·
sin
⁡
θ
+
s
·
cos
⁡
θ
d
s
, where(35)
1
C
ω
x
,
ω
y
=
1
if
ω
x
,
ω
y
∈
C
0
o
t
h
e
r
w
i
s
e
.The Radon transform (
R
d
ϕ
k
l
/
d
t
)
(
ρ
,
θ
,
t
) evaluated at a particular point (
ρ
0
,
θ
0
,
t
0
) is essentially an integral of d
ϕ
k
l
/
d
t
(
ω
x
,
ω
y
,
t
0
) along a line oriented at angle π
/
2
+
θ
0 with the ω
x axis and at distance |
ρ
0
| along the cos
⁡
θ
0
,
sin
⁡
(
θ
0
) direction from (
0,0
).If, for a particulark and l, we have d
ϕ
k
l
/
d
t
(
ω
x
,
ω
y
,
t
)
=
-
v
x
(
t
)
ω
x
-
v
y
(
t
)
ω
y, we have(36)
R
d
ϕ
k
l
/
d
t
ρ
,
θ
,
t
c
ρ
,
θ
=
ρ
-
v
x
t
cos
⁡
θ
-
v
y
t
sin
⁡
θ
,where(37)
c
ρ
,
θ
=
∫
R
1
C
ρ
·
cos
⁡
θ
-
s
·
sin
⁡
θ
,
ρ
·
sin
⁡
θ
+
s
·
cos
⁡
θ
d
s
.
c
(
ρ
,
θ
) is a correction term due to different length of line integrals for different values of (
ρ
,
θ
) in the bounded domain C.After computing the Radon transform ofd
ϕ
k
l
/
d
t
(
ω
x
,
ω
y
,
t
) for every block (
k
,
l
) at time t
0, we compute the Phase Motion Indicator (PMI), defined as(38)
P
M
I
k
l
=
max
θ
∈
0
,
π
⁡
∑
ρ
R
d
ϕ
k
l
/
d
t
ρ
,
θ
,
t
0
c
ρ
,
θ
.If the P
M
I
k
l is larger than a chosen threshold, motion is deemed to occur in block (
k
,
l
) at time t
0.Using the Radon transform makes it easier to separate rigid motion from noise. Since the phase is quite sensitive to noise, particularly when the amplitude is very small, the change of phase under noise may have comparable magnitude to that due to motion as mentioned earlier. The change of phase under noise, however, does not possess the structure suggested by (29) in the (
ω
x
,
ω
y
) domain. Instead, it appears to be more randomly distributed. Consequently, the PMI value is comparatively small for these blocks (see also Section 3.3.3).Moreover, the direction of motion, for block(
k
,
l
) where motion is detected, can be easily computed as(39)
θ
^
k
l
=
π
sign
⁡
∑
ρ
>
0
R
d
ϕ
k
l
/
d
t
ρ
,
α
k
l
,
t
0
/
c
ρ
,
α
k
l
+
1
2
+
α
k
l
,where(40)
α
k
l
=
argmax
θ
∈
0
,
π
⁡
∑
ρ
R
d
ϕ
k
l
/
d
t
ρ
,
θ
,
t
0
c
ρ
,
θ
.This follows from (36).


#### 3.3.2. The Phase-Based Motion Detection Algorithm

We formally summarize the above analysis as Algorithm1. Figure 5 shows a schematic diagram of the proposed phase-based motion detection algorithm.Algorithm 1:Phase-based motion detection algorithm using the FFT.
Input: Visual stream or Video data u
(
x
,
y
,
t
)
Output: Detected Motion
Construct Gaussian window of sizeM
×
M, denote as w;
for  t
=
Δ
T
,
2
Δ
T
,
…
,
N
Δ
T  do
Divide the screen of frame at time t into multiple M
×
M blocks u
k
l
(
t
), with overlaps;
foreach  block  u
k
l
(
t
)  do
Multiply pixel-wise the blocku
k
l
(
t
) with the Gaussian window u
k
l
(
t
)
w;
TakeM
×
M 2D FFT of u
k
l
(
t
)
w, denote phase as ϕ
k
l
(
ω
x
,
ω
y
,
t
) and amplitude as A
k
l
(
ω
x
,
ω
y
,
t
);
foreach  frequency  (
ω
x
,
ω
y
)  do
Highpass filter, temporally,ϕ
k
l
(
ω
x
,
ω
y
,
t
), denote ϕ
k
l
′
(
ω
x
,
ω
y
,
t
);
Optionally, denoise byϕ
k
l
′
(
ω
x
,
ω
y
,
t
)
←
ϕ
k
l
′
(
ω
x
,
ω
y
,
t
)
A
k
l
(
ω
x
,
ω
y
,
t
)
/
1
/
M
2
∑
(
ω
x
,
ω
y
)
∈
D
2
‍
A
k
l
(
ω
x
,
ω
y
,
t
)
+
ϵ;
end
Compute the radon transform(
R
ϕ
k
l
′
)
(
ρ
,
θ
,
t
);
ComputeP
M
I
k
l
(
t
) according to (38);
if  P
M
I
k
l
(
t
) > threshold  then
Motion Detected at block(
k
,
l
) at time t;
Compute direction of motion according to (39);
end
end
endFigure 5
Schematic diagram of proposed phase-based motion detection algorithm.The algorithm is subdivided into two parts. The first part computes local phase changes and the second part is the phase-based motion detector.In the first part, the screen is divided into overlapping blocks. For example, the red, green, and blue blocks in the plane “divide into overlapping blocks” correspond to the squares of the same color covering the video stream. A Gaussian window is then applied on each block, followed by a 2D FFT operation that is used to extract the local phase. A temporal high-pass filter is then employed to extract phase changes.In the second part, the PMI is evaluated for each block based on the Radon transform of the local phase changes in each block. Motion is detected for blocks with PMI larger than a preset threshold, and the direction of motion is computed as in (39).It is easy to notice that the algorithm can be highly parallelized.


#### 3.3.3. Example

We provide an illustrative example in Figure6 showing how motion is detected using Algorithm 1. The full video of this example can be found in Supplementary Video S1; see Supplementary Material available online at http://dx.doi.org/10.1155/2016/7915245. Figure 6(a) depicts a still from the “highway video” in the Change Detection 2014 dataset [22] evaluated at a particular time t
0. As suggested by the algorithm, the screen in Figure 6(a) is divided into 26
×
19 overlapping blocks and the window functions are applied to each block. Local phases can then be extracted from the 2D FFT of each windowed block, and the local phase changes are obtained by temporal high-pass filtering. The phase change is shown in Figure 6(b) for all blocks, with block (
12,11
) enlarged in Figure 6(c) and block (
23,6
) enlarged in Figure 6(d) (see also the plane “2D FFT and extract phase change” in Figure 5). Note that at the time of the video frame, block (
12,11
) covers a part of the vehicle in motion in the front, and block (
23,6
) corresponds to an area of the highway pavement where no motion occurs.Figure 6
Evaluating the motion detection algorithm at timet
0 (see also Supplementary Video S1). (a) A still from a 156
×
112-pixel video at time t
0. (b) d
ϕ
k
l
/
d
t
(
ω
x
,
ω
y
,
t
0
) for all 26
×
19 blocks, each of size 32
×
32. The blocks are concatenated to create a large “phase image” with their relative neighbors kept. (c) d
ϕ
12,11
/
d
t
(
ω
x
,
ω
y
,
t
0
) shown in the (
ω
x
,
ω
y
)
∈
D
2 space. (d) d
ϕ
23,6
/
d
t
(
ω
x
,
ω
y
,
t
0
) shown in the (
ω
x
,
ω
y
)
∈
D
2 space. (e) Phase Motion Indicator (PMI) for all blocks at t
0. Each pixel (
k
,
l
) is color coded for P
M
I
k
l. (f) max
ω
x
,
ω
y
∈
D
2
⁡
d
ϕ
k
l
/
d
t
ω
x
,
ω
y
,
t
0 for all blocks at t
0. Each pixel (
k
,
l
) is color coded for this maximum value of block (
k
,
l
).Figure6(f) depicts, for each block (
k
,
l
), the maximum phase change over all (
ω
x
,
ω
y
)
∈
D
2; that is,(41)
max
ω
x
,
ω
y
∈
D
2
⁡
d
ϕ
k
l
d
t
ω
x
,
ω
y
,
t
0
. We observe from the figure that, for regions with low amplitude, such as the region depicting the road, when the normalization constant is absent, the derivative of the phase can be noisy. For these blocks the maximum of d
ϕ
k
l
/
d
t
ω
x
,
ω
y
,
t
0 over all (
ω
x
,
ω
y
)
∈
D
2 is comparable to the maximum obtained for blocks that cover the vehicles in motion.However, (29) suggests that the local phase change from multiple filter pairs centered at the same spatial position (
k
,
l
) can provide a constraint to robustly estimate motion and its direction. Given the block structure employed in the computation of the local phase, it is natural to utilize phase change information from multiple sources.Indeed, if, for a particular block(
k
,
l
), d
ϕ
k
l
/
d
t
(
ω
x
,
ω
y
,
t
)
=
-
v
x
(
t
)
ω
x
-
v
y
(
t
)
ω
y, then it is easy to see that d
ϕ
k
l
/
d
t
(
ω
x
,
ω
y
,
t
) will be zero on the line v
x
(
t
)
ω
x
+
v
y
(
t
)
ω
y
=
0 and have opposite sign on either side of this line. For example, in Figures 6(b) and 6(c), d
ϕ
k
l
/
d
t
(
ω
x
,
ω
y
,
t
0
) clearly exhibits this property for blocks that cover a vehicle in motion. The PMI is a tool to evaluate this property.Finally, the PMIs for all blocks are shown compactly in a heat map in Figure6(e). The figure shows clearly that the blocks corresponding to the two moving vehicles have a high PMI value while the stationary background areas have a low PMI value, allowing one to easily detect motion by employing simple thresholding (see also the plane “Radon transform and extract strength orientation of plane” in Figure 5). In addition, the orientation of motion in each block is readily observable even by inspection in Figure 6(b) by a line separating the yellow part and blue part in each block. Further results about the direction of motion are presented in Section 4.


### 3.4. Relationship to Biological Motion Detectors

A straightforward way to implement local motion detectors is to apply a complex-valued Gabor receptive field (5) to the video signal u, and then take the derivative of the phase with respect to time or apply a high-pass filter on the phase to approximate the derivative.We present here an alternate implementation without explicitly computing the phase. This will elucidate the relation between the phase-based motion detector presented in the previous section and some elementary motion detection models used in biology, such as the Reichardt motion detector [7] and motion energy detector [9].Assuming that the local phaseϕ
00
(
ω
x
,
ω
y
,
t
) is differentiable, we have (see also the Appendix)(42)
d
ϕ
00
ω
x
,
ω
y
,
t
d
t
=
d
b
ω
x
,
ω
y
,
t
/
d
t
a
ω
x
,
ω
y
,
t
-
d
a
ω
x
,
ω
y
,
t
/
d
t
b
ω
x
,
ω
y
,
t
a
ω
x
,
ω
y
,
t
2
+
b
ω
x
,
ω
y
,
t
2
,where a
(
ω
x
,
ω
y
,
t
) and b
(
ω
x
,
ω
y
,
t
) are, respectively, the real and imaginary parts of A
00
(
ω
x
,
ω
y
,
t
)
e
j
ϕ
00
(
ω
x
,
ω
y
,
t
).We notice that the denominator of (42) is the square of the local amplitude of u, and the numerator is of the form of a second order Volterra kernel. This suggests that the time derivative of the local phase can be viewed as a second order Volterra kernel that processes twonormalized spatially filtered inputs v
1 and v
2.We consider an elaborated Reichardt motion detector as shown in Figure7. It is equipped with a quadrature pair of Gabor filters whose outputs are r
1
(
t
)
=
a
(
ω
x
,
ω
y
,
t
) and r
2
(
t
)
=
b
(
ω
x
,
ω
y
,
t
), respectively, for a particular value of (
ω
x
,
ω
y
). The pair of Gabor filters that provide these outputs are the real and imaginary parts of w
(
x
,
y
)
e
-
j
(
ω
x
x
+
ω
y
y
). It also consists of a temporal high-pass g
1
(
t
) and temporal low-pass filter g
2
(
t
) [9]. The output of the elaborated Reichardt detector follows the diagram in Figure 7 and can be expressed as(43)
r
2
∗
g
1
t
r
1
∗
g
2
t
-
r
1
∗
g
1
t
r
2
∗
g
2
t
.The response can also be characterized by a second order Volterra kernel. We notice the striking similarity between (43) and the numerator of (42). In fact, the phase-based motion detector shares some properties with the Reichardt motion detector. For example, it is straightforward to see that a single phase-based motion detector is tuned to the temporal frequency of a moving sinusoidal grating.Figure 7
Diagram of an elaborated Reichardt motion detector. It consists of a quadrature pair of Gabor filters whose outputsr
1 and r
2 provide inputs to the high-pass filters g
1
(
t
) and the low-pass filters g
2
(
t
).Since the motion energy detector is formally equivalent to an elaborated Reichardt motion detector [9], the structure of the motion energy detector with divisive normalization is also similar to the phase-based motion detector.


## 3.1. The Global Phase Equation for Translational Motion

Letu
=
u
(
x
,
y
,
t
), (
x
,
y
)
∈
R
2, t
∈
R, be a visual stimulus. If the visual stimulus is a pure translation of the signal at u
(
x
,
y
,
0
), that is,(22)
u
x
,
y
,
t
=
u
x
-
s
x
t
,
y
-
s
y
t
,
0
,where(23)
s
x
t
=
∫
0
t
v
x
s
d
s
,
s
y
t
=
∫
0
t
v
y
s
d
sare the total length of translation at time t in each dimension and v
x
(
t
) and v
y
(
t
) are the corresponding instantaneous velocity components, then the only difference between u
(
x
,
y
,
t
) and u
(
x
,
y
,
0
) in the Fourier domain is captured by their global phase. More formally, consider the following.Lemma 6.
The change (derivative) of the global phase is given by(24)
d
ϕ
^
ω
x
,
ω
y
,
t
d
t
=
-
ω
x
v
x
t
-
ω
y
v
y
t
=
-
ω
x
,
ω
y
v
x
t
,
v
y
t
T
,where, by abuse of notation, ϕ
^
ω
x
,
ω
y
,
t denotes the global phase of u
x
,
y
,
t and ϕ
^
ω
x
,
ω
y
,
0 is the initial condition.Proof.
IfF
u
(
·
,
·
,
0
)  (
ω
x
,
ω
y
) is the 2D (spatial) Fourier transform of u
=
u
(
x
,
y
,
0
), (
x
,
y
)
∈
R
2, by the Fourier shift theorem, we have(25)
F
u
·
,
·
,
t
ω
x
,
ω
y
=
F
u
·
,
·
,
0
ω
x
,
ω
y
e
-
j
ω
x
s
x
t
+
ω
y
s
y
t
.For a certain frequency component (
ω
x
,
ω
y
), the change of its global phase over time amounts to(26)
d
ϕ
^
ω
x
,
ω
y
,
t
d
t
=
-
ω
x
d
s
x
t
d
t
-
ω
y
d
s
y
t
d
t
=
-
ω
x
v
x
t
-
ω
y
v
y
t
=
-
ω
x
,
ω
y
v
x
t
,
v
y
t
T
.Therefore, in the simple case where the entire visual field is shifting, the derivative of the phase of Fourier components indicates motion, and it can be obtained by the inner product between the component frequency and the velocity vector.


## 3.2. The Change of Local Phase

### 3.2.1. The Local Phase Equation for Translational Motion

The analysis in Section3.1 applies toglobal motion. This type of motion occurs most frequently when the imaging device, either an eye or a camera, moves. Visual motion in the natural environment, however, is more diverse across the screen since it is, often time, produced by multiple moving objects. The objects can be small and the motion more localized in the visual field.Taking the global phase ofu will not simply reveal where motion of independent objects takes place or their direction/velocity of motion. The ease of interpretation of motion by using the Fourier transform, however, motivates us to reuse the same concept in detecting local motion. This can be achieved by restricting the domain of the visual field where the Fourier transform is applied.To be able to detect local motion, we consider the local phase ofu
(
x
,
y
,
t
) by taking the STFT with window function w
(
x
,
y
). Note that, the STFT and its ubiquitous implementation in DSP chips can be extensively used in any dimension. For simplicity and without loss of generality, we consider the window to be centered at (
0,0
). The STFT is given by(27)
∫
R
2
u
x
,
y
,
t
w
x
,
y
e
-
j
ω
x
x
+
ω
y
y
d
x
d
y
=
A
00
ω
x
,
ω
y
,
t
e
j
ϕ
00
ω
x
,
ω
y
,
t
,where, by abuse of notation, A
00
(
ω
x
,
ω
y
,
t
) is the amplitude and ϕ
00
(
ω
x
,
ω
y
,
t
) the local phase.Before we move on to the mathematical analysis, we can intuitively explain the relation between the change in local phase and visual motion taking place across the window support. First, if the stimulus undergoes a uniform change of intensity or it changes proportionally over time due to lighting conditions, for example, the local phase does not change since the phase is invariant with respect to intensity scaling. Therefore, the local phase does not change for such nonmotion stimuli. Second, a rigid edge moving across the window support will induce a phase change.For a strictly translational signal within the window support (footprint), for example,(28)
u
x
,
y
,
t
=
u
x
-
s
x
t
,
y
-
s
y
t
,
0
,
f
o
r
x
,
y
∈
supp
⁡
w
x
,
y
, where s
x
(
t
) and s
y
(
t
) are as defined in (23), we have the following result.Lemma 7.
Consider the following:(29)
d
ϕ
00
d
t
ω
x
,
ω
y
,
t
=
-
d
s
x
t
d
t
ω
x
-
d
s
y
t
d
t
ω
y
+
v
00
ω
x
,
ω
y
,
t
,where, by abuse of notation, ϕ
00
(
ω
x
,
ω
y
,
t
) is the local phase of u
(
x
,
y
,
t
) and ϕ
00
(
ω
x
,
ω
y
,
0
) is the initial condition. The functional form of the term v
00
=
v
00
(
ω
x
,
ω
y
,
t
) is provided in the Appendix.Proof.
The derivation of (29) and the functional form of v
00 are given in the Appendix.Remark 8.
We notice that the derivative of the local phase has similar structure to that of the global phase, but for the added termv
00. Through simulations, we observed that the first two terms in (29) dominate over the last term for an ON or OFF moving edge [3]. For example, Figure 3 shows the derivative of the local phase given in (29) for an ON edge moving with velocity (
40,0
) pixels/sec.Figure 3
Derivative of the local phase for an ON edge moving in the positivex direction.Remark 9.
Note thatϕ
00
(
ω
x
,
ω
y
,
t
) may not be differentiable even if u
(
x
,
y
,
t
) is differentiable, particularly when A
00
(
ω
x
,
ω
y
,
t
)
=
0. For example, the spatial phase can jump from a positive value to zero when A
00
(
ω
x
,
ω
y
,
t
) diminishes. This also suggests that the instantaneous local spatial phase is less informative about a region of a visual scene whenever A
00
(
ω
x
,
ω
y
,
t
) is close to zero. Nevertheless, the time derivative of the local phase can be approximated by applying a high-pass filter to ϕ
00
(
ω
x
,
ω
y
,
t
).


### 3.2.2. The Block Structure for Computing the Local Phase

We construct Gaussian windows along thex, y dimensions. The Gaussian windows are defined as(30)
T
k
l
w
x
,
y
=
e
-
x
-
x
k
2
+
x
-
y
l
2
/
2
σ
2
,where x
k
=
k
b
0, y
l
=
l
b
0, in which b
0
∈
Z
+ is the distance between two neighboring windows and 1
≤
k
b
0
≤
P
x, 1
≤
l
b
0
≤
P
y, where P
x
,
P
y
∈
Z
+ are the number of pixels of the screen in x and y directions, respectively.We then take the 2D Fourier transform of the windowed video signalu
(
x
,
y
,
t
)
(
T
k
l
w
)
(
x
,
y
) and write in polar form(31)
∫
R
2
u
x
,
y
,
t
T
k
l
w
x
,
y
e
-
j
ω
x
x
-
x
k
+
ω
y
y
-
y
l
d
x
d
y
=
A
k
l
ω
x
,
ω
y
,
t
e
j
ϕ
k
l
ω
x
,
ω
y
,
t
.The above integral can be very efficiently evaluated using the 2D FFT in discrete domain defined on M
×
M blocks approximating the footprint of the Gaussian windows. For example, the standard deviation of the Gaussian windows we use in the examples in Section 4 is 4 pixels. A block of 32
×
32 pixels (M
=
32) is sufficient to cover the effective support (or footprint) of the Gaussian window. At the same time, the size of the block is a power of 2, which is most suitable for FFT-based computation. The processing of each block (
k
,
l
) is independent of all other blocks; thereby, parallelism is readily achieved.Note that the size of the window is informed by the size of the objects one is interested in locating. Measurements of the local phase using smaller window functions are less robust to noise. Larger windows would enhance object motion detection if the object size is comparable to the window size. However, there would be an increased likelihood of independent movement of multiple objects within the same window, which is not modeled here and thereby may not be robustly detected.Therefore, for each block(
k
,
l
), we obtain M
2 measurements of the phase ϕ
k
l
(
ω
x
m
,
ω
y
m
,
t
) at every time instant t, with (
ω
x
m
,
ω
y
m
)
∈
D
2, where(32)
D
2
=
ω
x
m
=
m
ω
0
,
ω
y
n
=
n
ω
0
,
m
,
n
=
-
M
/
2
,
-
M
/
2
+
1
,
…
,
M
/
2
-
1
, with ω
0
=
2
π
/
M. We then compute the temporal derivative of the phase, that is, d
ϕ
k
l
/
d
t
(
ω
x
,
ω
y
,
t
) for (
ω
x
,
ω
y
)
∈
D
2.We further illustrate an example of the block structure in Figure4. Figure 4(a) shows an example of an image of 64
×
64 pixels. Four Gaussian windows are shown each with a standard deviation of 4 pixels. The distance between the centers of two neighboring Gaussian windows is 6 pixels. The red solid square shows a 32
×
32-pixel block with k
=
3, l
=
3, which encloses effective support of the Gaussian window on top-left (k
=
0, l
=
0 is the block with Gaussian window centered at pixel (
1,1
)). The green dashed square shows another 32
×
32-pixel block with k
=
3, l
=
7. The two Gaussian windows on the right are associated with the blocks k
=
7, l
=
3 and k
=
7, l
=
7, respectively. Cross section of all Gaussian windows with k
=
7, l
∈
[
0,10
], that is, those centered on the magenta line, are shown in Figure 4(b). The red and blue curve in Figure 4(b) correspond to the two Gaussian windows shown in Figure 4(a). Figure 4(b) also suggests that some of the Gaussian windows are cut off on the boundaries. This is, however, equivalent to assuming that the pixel values outside the boundary are always zero, and it will not significantly affect motion detection based on the change of local phase.Figure 4
Example of block structure. (a) Four yellow disks show 4 Gaussian window functions with translation parameterk
=
3, l
=
3 (top-left), k
=
3, l
=
7 (bottom-left), k
=
7, l
=
3 (top-right), and k
=
7, l
=
7 (bottom-right). The red solid square shows a 32
×
32-pixel block approximating the Gaussian window with k
=
3, l
=
3, and the green dashed square shows a 32
×
32-pixel block approximating the Gaussian window with k
=
3, l
=
7. (b) Cross section of all 11 Gaussian windows with translation parameters k
=
7, l
∈
[
0,10
]. The cross section is taken as indicated by the magenta line in (a), and the red and blue curve correspond to cross sections of the two Gaussian windows shown in (a) centered on the magenta line.
(a)
(b)Since the phase, and thereby the phase change, is noisier when the local amplitude is low, an additional denoising step can be employed to discount the measurements ofd
ϕ
k
l
/
d
t
(
ω
x
,
ω
y
,
t
) for low amplitude values A
k
l
(
ω
x
,
ω
y
,
t
). The denoising is given by(33)
d
ϕ
k
l
d
t
ω
x
,
ω
y
,
t
A
k
l
ω
x
,
ω
y
,
t
1
/
M
2
∑
ω
x
,
ω
y
∈
D
2
A
k
l
ω
x
,
ω
y
,
t
+
ϵ
,where ϵ
>
0 is a constant, and (
ω
x
,
ω
y
)
∈
D
2.


## 3.2.1. The Local Phase Equation for Translational Motion

The analysis in Section3.1 applies toglobal motion. This type of motion occurs most frequently when the imaging device, either an eye or a camera, moves. Visual motion in the natural environment, however, is more diverse across the screen since it is, often time, produced by multiple moving objects. The objects can be small and the motion more localized in the visual field.Taking the global phase ofu will not simply reveal where motion of independent objects takes place or their direction/velocity of motion. The ease of interpretation of motion by using the Fourier transform, however, motivates us to reuse the same concept in detecting local motion. This can be achieved by restricting the domain of the visual field where the Fourier transform is applied.To be able to detect local motion, we consider the local phase ofu
(
x
,
y
,
t
) by taking the STFT with window function w
(
x
,
y
). Note that, the STFT and its ubiquitous implementation in DSP chips can be extensively used in any dimension. For simplicity and without loss of generality, we consider the window to be centered at (
0,0
). The STFT is given by(27)
∫
R
2
u
x
,
y
,
t
w
x
,
y
e
-
j
ω
x
x
+
ω
y
y
d
x
d
y
=
A
00
ω
x
,
ω
y
,
t
e
j
ϕ
00
ω
x
,
ω
y
,
t
,where, by abuse of notation, A
00
(
ω
x
,
ω
y
,
t
) is the amplitude and ϕ
00
(
ω
x
,
ω
y
,
t
) the local phase.Before we move on to the mathematical analysis, we can intuitively explain the relation between the change in local phase and visual motion taking place across the window support. First, if the stimulus undergoes a uniform change of intensity or it changes proportionally over time due to lighting conditions, for example, the local phase does not change since the phase is invariant with respect to intensity scaling. Therefore, the local phase does not change for such nonmotion stimuli. Second, a rigid edge moving across the window support will induce a phase change.For a strictly translational signal within the window support (footprint), for example,(28)
u
x
,
y
,
t
=
u
x
-
s
x
t
,
y
-
s
y
t
,
0
,
f
o
r
x
,
y
∈
supp
⁡
w
x
,
y
, where s
x
(
t
) and s
y
(
t
) are as defined in (23), we have the following result.Lemma 7.
Consider the following:(29)
d
ϕ
00
d
t
ω
x
,
ω
y
,
t
=
-
d
s
x
t
d
t
ω
x
-
d
s
y
t
d
t
ω
y
+
v
00
ω
x
,
ω
y
,
t
,where, by abuse of notation, ϕ
00
(
ω
x
,
ω
y
,
t
) is the local phase of u
(
x
,
y
,
t
) and ϕ
00
(
ω
x
,
ω
y
,
0
) is the initial condition. The functional form of the term v
00
=
v
00
(
ω
x
,
ω
y
,
t
) is provided in the Appendix.Proof.
The derivation of (29) and the functional form of v
00 are given in the Appendix.Remark 8.
We notice that the derivative of the local phase has similar structure to that of the global phase, but for the added termv
00. Through simulations, we observed that the first two terms in (29) dominate over the last term for an ON or OFF moving edge [3]. For example, Figure 3 shows the derivative of the local phase given in (29) for an ON edge moving with velocity (
40,0
) pixels/sec.Figure 3
Derivative of the local phase for an ON edge moving in the positivex direction.Remark 9.
Note thatϕ
00
(
ω
x
,
ω
y
,
t
) may not be differentiable even if u
(
x
,
y
,
t
) is differentiable, particularly when A
00
(
ω
x
,
ω
y
,
t
)
=
0. For example, the spatial phase can jump from a positive value to zero when A
00
(
ω
x
,
ω
y
,
t
) diminishes. This also suggests that the instantaneous local spatial phase is less informative about a region of a visual scene whenever A
00
(
ω
x
,
ω
y
,
t
) is close to zero. Nevertheless, the time derivative of the local phase can be approximated by applying a high-pass filter to ϕ
00
(
ω
x
,
ω
y
,
t
).


## 3.2.2. The Block Structure for Computing the Local Phase

We construct Gaussian windows along thex, y dimensions. The Gaussian windows are defined as(30)
T
k
l
w
x
,
y
=
e
-
x
-
x
k
2
+
x
-
y
l
2
/
2
σ
2
,where x
k
=
k
b
0, y
l
=
l
b
0, in which b
0
∈
Z
+ is the distance between two neighboring windows and 1
≤
k
b
0
≤
P
x, 1
≤
l
b
0
≤
P
y, where P
x
,
P
y
∈
Z
+ are the number of pixels of the screen in x and y directions, respectively.We then take the 2D Fourier transform of the windowed video signalu
(
x
,
y
,
t
)
(
T
k
l
w
)
(
x
,
y
) and write in polar form(31)
∫
R
2
u
x
,
y
,
t
T
k
l
w
x
,
y
e
-
j
ω
x
x
-
x
k
+
ω
y
y
-
y
l
d
x
d
y
=
A
k
l
ω
x
,
ω
y
,
t
e
j
ϕ
k
l
ω
x
,
ω
y
,
t
.The above integral can be very efficiently evaluated using the 2D FFT in discrete domain defined on M
×
M blocks approximating the footprint of the Gaussian windows. For example, the standard deviation of the Gaussian windows we use in the examples in Section 4 is 4 pixels. A block of 32
×
32 pixels (M
=
32) is sufficient to cover the effective support (or footprint) of the Gaussian window. At the same time, the size of the block is a power of 2, which is most suitable for FFT-based computation. The processing of each block (
k
,
l
) is independent of all other blocks; thereby, parallelism is readily achieved.Note that the size of the window is informed by the size of the objects one is interested in locating. Measurements of the local phase using smaller window functions are less robust to noise. Larger windows would enhance object motion detection if the object size is comparable to the window size. However, there would be an increased likelihood of independent movement of multiple objects within the same window, which is not modeled here and thereby may not be robustly detected.Therefore, for each block(
k
,
l
), we obtain M
2 measurements of the phase ϕ
k
l
(
ω
x
m
,
ω
y
m
,
t
) at every time instant t, with (
ω
x
m
,
ω
y
m
)
∈
D
2, where(32)
D
2
=
ω
x
m
=
m
ω
0
,
ω
y
n
=
n
ω
0
,
m
,
n
=
-
M
/
2
,
-
M
/
2
+
1
,
…
,
M
/
2
-
1
, with ω
0
=
2
π
/
M. We then compute the temporal derivative of the phase, that is, d
ϕ
k
l
/
d
t
(
ω
x
,
ω
y
,
t
) for (
ω
x
,
ω
y
)
∈
D
2.We further illustrate an example of the block structure in Figure4. Figure 4(a) shows an example of an image of 64
×
64 pixels. Four Gaussian windows are shown each with a standard deviation of 4 pixels. The distance between the centers of two neighboring Gaussian windows is 6 pixels. The red solid square shows a 32
×
32-pixel block with k
=
3, l
=
3, which encloses effective support of the Gaussian window on top-left (k
=
0, l
=
0 is the block with Gaussian window centered at pixel (
1,1
)). The green dashed square shows another 32
×
32-pixel block with k
=
3, l
=
7. The two Gaussian windows on the right are associated with the blocks k
=
7, l
=
3 and k
=
7, l
=
7, respectively. Cross section of all Gaussian windows with k
=
7, l
∈
[
0,10
], that is, those centered on the magenta line, are shown in Figure 4(b). The red and blue curve in Figure 4(b) correspond to the two Gaussian windows shown in Figure 4(a). Figure 4(b) also suggests that some of the Gaussian windows are cut off on the boundaries. This is, however, equivalent to assuming that the pixel values outside the boundary are always zero, and it will not significantly affect motion detection based on the change of local phase.Figure 4
Example of block structure. (a) Four yellow disks show 4 Gaussian window functions with translation parameterk
=
3, l
=
3 (top-left), k
=
3, l
=
7 (bottom-left), k
=
7, l
=
3 (top-right), and k
=
7, l
=
7 (bottom-right). The red solid square shows a 32
×
32-pixel block approximating the Gaussian window with k
=
3, l
=
3, and the green dashed square shows a 32
×
32-pixel block approximating the Gaussian window with k
=
3, l
=
7. (b) Cross section of all 11 Gaussian windows with translation parameters k
=
7, l
∈
[
0,10
]. The cross section is taken as indicated by the magenta line in (a), and the red and blue curve correspond to cross sections of the two Gaussian windows shown in (a) centered on the magenta line.
(a)
(b)Since the phase, and thereby the phase change, is noisier when the local amplitude is low, an additional denoising step can be employed to discount the measurements ofd
ϕ
k
l
/
d
t
(
ω
x
,
ω
y
,
t
) for low amplitude values A
k
l
(
ω
x
,
ω
y
,
t
). The denoising is given by(33)
d
ϕ
k
l
d
t
ω
x
,
ω
y
,
t
A
k
l
ω
x
,
ω
y
,
t
1
/
M
2
∑
ω
x
,
ω
y
∈
D
2
A
k
l
ω
x
,
ω
y
,
t
+
ϵ
,where ϵ
>
0 is a constant, and (
ω
x
,
ω
y
)
∈
D
2.


## 3.3. The Phase-Based Detector

We propose here a block FFT based algorithm to detect motion using phase information. Such an algorithm is, due to its simplicity and parallelism, highly suitable for anin silico implementation.


### 3.3.1. Radon Transform on the Change of Phases

We exploit the approximately linear structure of the phase derivative for blocks exhibiting motion by computing the Radon transform ofd
ϕ
k
l
/
d
t
(
ω
x
,
ω
y
,
t
) over a circular bounded domain C
=
{
(
ω
x
,
ω
y
)
∣
(
ω
x
,
ω
y
)
∈
D
2
,
ω
x
2
+
ω
y
2
<
π
2
}.The Radon transform of the change of phase in the domainC is given by(34)
R
d
ϕ
k
l
d
t
ρ
,
θ
,
t
=
∫
R
d
ϕ
k
l
d
t
ρ
·
cos
⁡
θ
-
s
·
sin
⁡
θ
,
ρ
·
sin
⁡
θ
+
s
·
cos
⁡
θ
,
t
1
C
ρ
·
cos
⁡
θ
-
s
·
sin
⁡
θ
,
ρ
·
sin
⁡
θ
+
s
·
cos
⁡
θ
d
s
, where(35)
1
C
ω
x
,
ω
y
=
1
if
ω
x
,
ω
y
∈
C
0
o
t
h
e
r
w
i
s
e
.The Radon transform (
R
d
ϕ
k
l
/
d
t
)
(
ρ
,
θ
,
t
) evaluated at a particular point (
ρ
0
,
θ
0
,
t
0
) is essentially an integral of d
ϕ
k
l
/
d
t
(
ω
x
,
ω
y
,
t
0
) along a line oriented at angle π
/
2
+
θ
0 with the ω
x axis and at distance |
ρ
0
| along the cos
⁡
θ
0
,
sin
⁡
(
θ
0
) direction from (
0,0
).If, for a particulark and l, we have d
ϕ
k
l
/
d
t
(
ω
x
,
ω
y
,
t
)
=
-
v
x
(
t
)
ω
x
-
v
y
(
t
)
ω
y, we have(36)
R
d
ϕ
k
l
/
d
t
ρ
,
θ
,
t
c
ρ
,
θ
=
ρ
-
v
x
t
cos
⁡
θ
-
v
y
t
sin
⁡
θ
,where(37)
c
ρ
,
θ
=
∫
R
1
C
ρ
·
cos
⁡
θ
-
s
·
sin
⁡
θ
,
ρ
·
sin
⁡
θ
+
s
·
cos
⁡
θ
d
s
.
c
(
ρ
,
θ
) is a correction term due to different length of line integrals for different values of (
ρ
,
θ
) in the bounded domain C.After computing the Radon transform ofd
ϕ
k
l
/
d
t
(
ω
x
,
ω
y
,
t
) for every block (
k
,
l
) at time t
0, we compute the Phase Motion Indicator (PMI), defined as(38)
P
M
I
k
l
=
max
θ
∈
0
,
π
⁡
∑
ρ
R
d
ϕ
k
l
/
d
t
ρ
,
θ
,
t
0
c
ρ
,
θ
.If the P
M
I
k
l is larger than a chosen threshold, motion is deemed to occur in block (
k
,
l
) at time t
0.Using the Radon transform makes it easier to separate rigid motion from noise. Since the phase is quite sensitive to noise, particularly when the amplitude is very small, the change of phase under noise may have comparable magnitude to that due to motion as mentioned earlier. The change of phase under noise, however, does not possess the structure suggested by (29) in the (
ω
x
,
ω
y
) domain. Instead, it appears to be more randomly distributed. Consequently, the PMI value is comparatively small for these blocks (see also Section 3.3.3).Moreover, the direction of motion, for block(
k
,
l
) where motion is detected, can be easily computed as(39)
θ
^
k
l
=
π
sign
⁡
∑
ρ
>
0
R
d
ϕ
k
l
/
d
t
ρ
,
α
k
l
,
t
0
/
c
ρ
,
α
k
l
+
1
2
+
α
k
l
,where(40)
α
k
l
=
argmax
θ
∈
0
,
π
⁡
∑
ρ
R
d
ϕ
k
l
/
d
t
ρ
,
θ
,
t
0
c
ρ
,
θ
.This follows from (36).


### 3.3.2. The Phase-Based Motion Detection Algorithm

We formally summarize the above analysis as Algorithm1. Figure 5 shows a schematic diagram of the proposed phase-based motion detection algorithm.Algorithm 1:Phase-based motion detection algorithm using the FFT.
Input: Visual stream or Video data u
(
x
,
y
,
t
)
Output: Detected Motion
Construct Gaussian window of sizeM
×
M, denote as w;
for  t
=
Δ
T
,
2
Δ
T
,
…
,
N
Δ
T  do
Divide the screen of frame at time t into multiple M
×
M blocks u
k
l
(
t
), with overlaps;
foreach  block  u
k
l
(
t
)  do
Multiply pixel-wise the blocku
k
l
(
t
) with the Gaussian window u
k
l
(
t
)
w;
TakeM
×
M 2D FFT of u
k
l
(
t
)
w, denote phase as ϕ
k
l
(
ω
x
,
ω
y
,
t
) and amplitude as A
k
l
(
ω
x
,
ω
y
,
t
);
foreach  frequency  (
ω
x
,
ω
y
)  do
Highpass filter, temporally,ϕ
k
l
(
ω
x
,
ω
y
,
t
), denote ϕ
k
l
′
(
ω
x
,
ω
y
,
t
);
Optionally, denoise byϕ
k
l
′
(
ω
x
,
ω
y
,
t
)
←
ϕ
k
l
′
(
ω
x
,
ω
y
,
t
)
A
k
l
(
ω
x
,
ω
y
,
t
)
/
1
/
M
2
∑
(
ω
x
,
ω
y
)
∈
D
2
‍
A
k
l
(
ω
x
,
ω
y
,
t
)
+
ϵ;
end
Compute the radon transform(
R
ϕ
k
l
′
)
(
ρ
,
θ
,
t
);
ComputeP
M
I
k
l
(
t
) according to (38);
if  P
M
I
k
l
(
t
) > threshold  then
Motion Detected at block(
k
,
l
) at time t;
Compute direction of motion according to (39);
end
end
endFigure 5
Schematic diagram of proposed phase-based motion detection algorithm.The algorithm is subdivided into two parts. The first part computes local phase changes and the second part is the phase-based motion detector.In the first part, the screen is divided into overlapping blocks. For example, the red, green, and blue blocks in the plane “divide into overlapping blocks” correspond to the squares of the same color covering the video stream. A Gaussian window is then applied on each block, followed by a 2D FFT operation that is used to extract the local phase. A temporal high-pass filter is then employed to extract phase changes.In the second part, the PMI is evaluated for each block based on the Radon transform of the local phase changes in each block. Motion is detected for blocks with PMI larger than a preset threshold, and the direction of motion is computed as in (39).It is easy to notice that the algorithm can be highly parallelized.


### 3.3.3. Example

We provide an illustrative example in Figure6 showing how motion is detected using Algorithm 1. The full video of this example can be found in Supplementary Video S1; see Supplementary Material available online at http://dx.doi.org/10.1155/2016/7915245. Figure 6(a) depicts a still from the “highway video” in the Change Detection 2014 dataset [22] evaluated at a particular time t
0. As suggested by the algorithm, the screen in Figure 6(a) is divided into 26
×
19 overlapping blocks and the window functions are applied to each block. Local phases can then be extracted from the 2D FFT of each windowed block, and the local phase changes are obtained by temporal high-pass filtering. The phase change is shown in Figure 6(b) for all blocks, with block (
12,11
) enlarged in Figure 6(c) and block (
23,6
) enlarged in Figure 6(d) (see also the plane “2D FFT and extract phase change” in Figure 5). Note that at the time of the video frame, block (
12,11
) covers a part of the vehicle in motion in the front, and block (
23,6
) corresponds to an area of the highway pavement where no motion occurs.Figure 6
Evaluating the motion detection algorithm at timet
0 (see also Supplementary Video S1). (a) A still from a 156
×
112-pixel video at time t
0. (b) d
ϕ
k
l
/
d
t
(
ω
x
,
ω
y
,
t
0
) for all 26
×
19 blocks, each of size 32
×
32. The blocks are concatenated to create a large “phase image” with their relative neighbors kept. (c) d
ϕ
12,11
/
d
t
(
ω
x
,
ω
y
,
t
0
) shown in the (
ω
x
,
ω
y
)
∈
D
2 space. (d) d
ϕ
23,6
/
d
t
(
ω
x
,
ω
y
,
t
0
) shown in the (
ω
x
,
ω
y
)
∈
D
2 space. (e) Phase Motion Indicator (PMI) for all blocks at t
0. Each pixel (
k
,
l
) is color coded for P
M
I
k
l. (f) max
ω
x
,
ω
y
∈
D
2
⁡
d
ϕ
k
l
/
d
t
ω
x
,
ω
y
,
t
0 for all blocks at t
0. Each pixel (
k
,
l
) is color coded for this maximum value of block (
k
,
l
).Figure6(f) depicts, for each block (
k
,
l
), the maximum phase change over all (
ω
x
,
ω
y
)
∈
D
2; that is,(41)
max
ω
x
,
ω
y
∈
D
2
⁡
d
ϕ
k
l
d
t
ω
x
,
ω
y
,
t
0
. We observe from the figure that, for regions with low amplitude, such as the region depicting the road, when the normalization constant is absent, the derivative of the phase can be noisy. For these blocks the maximum of d
ϕ
k
l
/
d
t
ω
x
,
ω
y
,
t
0 over all (
ω
x
,
ω
y
)
∈
D
2 is comparable to the maximum obtained for blocks that cover the vehicles in motion.However, (29) suggests that the local phase change from multiple filter pairs centered at the same spatial position (
k
,
l
) can provide a constraint to robustly estimate motion and its direction. Given the block structure employed in the computation of the local phase, it is natural to utilize phase change information from multiple sources.Indeed, if, for a particular block(
k
,
l
), d
ϕ
k
l
/
d
t
(
ω
x
,
ω
y
,
t
)
=
-
v
x
(
t
)
ω
x
-
v
y
(
t
)
ω
y, then it is easy to see that d
ϕ
k
l
/
d
t
(
ω
x
,
ω
y
,
t
) will be zero on the line v
x
(
t
)
ω
x
+
v
y
(
t
)
ω
y
=
0 and have opposite sign on either side of this line. For example, in Figures 6(b) and 6(c), d
ϕ
k
l
/
d
t
(
ω
x
,
ω
y
,
t
0
) clearly exhibits this property for blocks that cover a vehicle in motion. The PMI is a tool to evaluate this property.Finally, the PMIs for all blocks are shown compactly in a heat map in Figure6(e). The figure shows clearly that the blocks corresponding to the two moving vehicles have a high PMI value while the stationary background areas have a low PMI value, allowing one to easily detect motion by employing simple thresholding (see also the plane “Radon transform and extract strength orientation of plane” in Figure 5). In addition, the orientation of motion in each block is readily observable even by inspection in Figure 6(b) by a line separating the yellow part and blue part in each block. Further results about the direction of motion are presented in Section 4.


## 3.3.1. Radon Transform on the Change of Phases

We exploit the approximately linear structure of the phase derivative for blocks exhibiting motion by computing the Radon transform ofd
ϕ
k
l
/
d
t
(
ω
x
,
ω
y
,
t
) over a circular bounded domain C
=
{
(
ω
x
,
ω
y
)
∣
(
ω
x
,
ω
y
)
∈
D
2
,
ω
x
2
+
ω
y
2
<
π
2
}.The Radon transform of the change of phase in the domainC is given by(34)
R
d
ϕ
k
l
d
t
ρ
,
θ
,
t
=
∫
R
d
ϕ
k
l
d
t
ρ
·
cos
⁡
θ
-
s
·
sin
⁡
θ
,
ρ
·
sin
⁡
θ
+
s
·
cos
⁡
θ
,
t
1
C
ρ
·
cos
⁡
θ
-
s
·
sin
⁡
θ
,
ρ
·
sin
⁡
θ
+
s
·
cos
⁡
θ
d
s
, where(35)
1
C
ω
x
,
ω
y
=
1
if
ω
x
,
ω
y
∈
C
0
o
t
h
e
r
w
i
s
e
.The Radon transform (
R
d
ϕ
k
l
/
d
t
)
(
ρ
,
θ
,
t
) evaluated at a particular point (
ρ
0
,
θ
0
,
t
0
) is essentially an integral of d
ϕ
k
l
/
d
t
(
ω
x
,
ω
y
,
t
0
) along a line oriented at angle π
/
2
+
θ
0 with the ω
x axis and at distance |
ρ
0
| along the cos
⁡
θ
0
,
sin
⁡
(
θ
0
) direction from (
0,0
).If, for a particulark and l, we have d
ϕ
k
l
/
d
t
(
ω
x
,
ω
y
,
t
)
=
-
v
x
(
t
)
ω
x
-
v
y
(
t
)
ω
y, we have(36)
R
d
ϕ
k
l
/
d
t
ρ
,
θ
,
t
c
ρ
,
θ
=
ρ
-
v
x
t
cos
⁡
θ
-
v
y
t
sin
⁡
θ
,where(37)
c
ρ
,
θ
=
∫
R
1
C
ρ
·
cos
⁡
θ
-
s
·
sin
⁡
θ
,
ρ
·
sin
⁡
θ
+
s
·
cos
⁡
θ
d
s
.
c
(
ρ
,
θ
) is a correction term due to different length of line integrals for different values of (
ρ
,
θ
) in the bounded domain C.After computing the Radon transform ofd
ϕ
k
l
/
d
t
(
ω
x
,
ω
y
,
t
) for every block (
k
,
l
) at time t
0, we compute the Phase Motion Indicator (PMI), defined as(38)
P
M
I
k
l
=
max
θ
∈
0
,
π
⁡
∑
ρ
R
d
ϕ
k
l
/
d
t
ρ
,
θ
,
t
0
c
ρ
,
θ
.If the P
M
I
k
l is larger than a chosen threshold, motion is deemed to occur in block (
k
,
l
) at time t
0.Using the Radon transform makes it easier to separate rigid motion from noise. Since the phase is quite sensitive to noise, particularly when the amplitude is very small, the change of phase under noise may have comparable magnitude to that due to motion as mentioned earlier. The change of phase under noise, however, does not possess the structure suggested by (29) in the (
ω
x
,
ω
y
) domain. Instead, it appears to be more randomly distributed. Consequently, the PMI value is comparatively small for these blocks (see also Section 3.3.3).Moreover, the direction of motion, for block(
k
,
l
) where motion is detected, can be easily computed as(39)
θ
^
k
l
=
π
sign
⁡
∑
ρ
>
0
R
d
ϕ
k
l
/
d
t
ρ
,
α
k
l
,
t
0
/
c
ρ
,
α
k
l
+
1
2
+
α
k
l
,where(40)
α
k
l
=
argmax
θ
∈
0
,
π
⁡
∑
ρ
R
d
ϕ
k
l
/
d
t
ρ
,
θ
,
t
0
c
ρ
,
θ
.This follows from (36).


## 3.3.2. The Phase-Based Motion Detection Algorithm

We formally summarize the above analysis as Algorithm1. Figure 5 shows a schematic diagram of the proposed phase-based motion detection algorithm.Algorithm 1:Phase-based motion detection algorithm using the FFT.
Input: Visual stream or Video data u
(
x
,
y
,
t
)
Output: Detected Motion
Construct Gaussian window of sizeM
×
M, denote as w;
for  t
=
Δ
T
,
2
Δ
T
,
…
,
N
Δ
T  do
Divide the screen of frame at time t into multiple M
×
M blocks u
k
l
(
t
), with overlaps;
foreach  block  u
k
l
(
t
)  do
Multiply pixel-wise the blocku
k
l
(
t
) with the Gaussian window u
k
l
(
t
)
w;
TakeM
×
M 2D FFT of u
k
l
(
t
)
w, denote phase as ϕ
k
l
(
ω
x
,
ω
y
,
t
) and amplitude as A
k
l
(
ω
x
,
ω
y
,
t
);
foreach  frequency  (
ω
x
,
ω
y
)  do
Highpass filter, temporally,ϕ
k
l
(
ω
x
,
ω
y
,
t
), denote ϕ
k
l
′
(
ω
x
,
ω
y
,
t
);
Optionally, denoise byϕ
k
l
′
(
ω
x
,
ω
y
,
t
)
←
ϕ
k
l
′
(
ω
x
,
ω
y
,
t
)
A
k
l
(
ω
x
,
ω
y
,
t
)
/
1
/
M
2
∑
(
ω
x
,
ω
y
)
∈
D
2
‍
A
k
l
(
ω
x
,
ω
y
,
t
)
+
ϵ;
end
Compute the radon transform(
R
ϕ
k
l
′
)
(
ρ
,
θ
,
t
);
ComputeP
M
I
k
l
(
t
) according to (38);
if  P
M
I
k
l
(
t
) > threshold  then
Motion Detected at block(
k
,
l
) at time t;
Compute direction of motion according to (39);
end
end
endFigure 5
Schematic diagram of proposed phase-based motion detection algorithm.The algorithm is subdivided into two parts. The first part computes local phase changes and the second part is the phase-based motion detector.In the first part, the screen is divided into overlapping blocks. For example, the red, green, and blue blocks in the plane “divide into overlapping blocks” correspond to the squares of the same color covering the video stream. A Gaussian window is then applied on each block, followed by a 2D FFT operation that is used to extract the local phase. A temporal high-pass filter is then employed to extract phase changes.In the second part, the PMI is evaluated for each block based on the Radon transform of the local phase changes in each block. Motion is detected for blocks with PMI larger than a preset threshold, and the direction of motion is computed as in (39).It is easy to notice that the algorithm can be highly parallelized.


## 3.3.3. Example

We provide an illustrative example in Figure6 showing how motion is detected using Algorithm 1. The full video of this example can be found in Supplementary Video S1; see Supplementary Material available online at http://dx.doi.org/10.1155/2016/7915245. Figure 6(a) depicts a still from the “highway video” in the Change Detection 2014 dataset [22] evaluated at a particular time t
0. As suggested by the algorithm, the screen in Figure 6(a) is divided into 26
×
19 overlapping blocks and the window functions are applied to each block. Local phases can then be extracted from the 2D FFT of each windowed block, and the local phase changes are obtained by temporal high-pass filtering. The phase change is shown in Figure 6(b) for all blocks, with block (
12,11
) enlarged in Figure 6(c) and block (
23,6
) enlarged in Figure 6(d) (see also the plane “2D FFT and extract phase change” in Figure 5). Note that at the time of the video frame, block (
12,11
) covers a part of the vehicle in motion in the front, and block (
23,6
) corresponds to an area of the highway pavement where no motion occurs.Figure 6
Evaluating the motion detection algorithm at timet
0 (see also Supplementary Video S1). (a) A still from a 156
×
112-pixel video at time t
0. (b) d
ϕ
k
l
/
d
t
(
ω
x
,
ω
y
,
t
0
) for all 26
×
19 blocks, each of size 32
×
32. The blocks are concatenated to create a large “phase image” with their relative neighbors kept. (c) d
ϕ
12,11
/
d
t
(
ω
x
,
ω
y
,
t
0
) shown in the (
ω
x
,
ω
y
)
∈
D
2 space. (d) d
ϕ
23,6
/
d
t
(
ω
x
,
ω
y
,
t
0
) shown in the (
ω
x
,
ω
y
)
∈
D
2 space. (e) Phase Motion Indicator (PMI) for all blocks at t
0. Each pixel (
k
,
l
) is color coded for P
M
I
k
l. (f) max
ω
x
,
ω
y
∈
D
2
⁡
d
ϕ
k
l
/
d
t
ω
x
,
ω
y
,
t
0 for all blocks at t
0. Each pixel (
k
,
l
) is color coded for this maximum value of block (
k
,
l
).Figure6(f) depicts, for each block (
k
,
l
), the maximum phase change over all (
ω
x
,
ω
y
)
∈
D
2; that is,(41)
max
ω
x
,
ω
y
∈
D
2
⁡
d
ϕ
k
l
d
t
ω
x
,
ω
y
,
t
0
. We observe from the figure that, for regions with low amplitude, such as the region depicting the road, when the normalization constant is absent, the derivative of the phase can be noisy. For these blocks the maximum of d
ϕ
k
l
/
d
t
ω
x
,
ω
y
,
t
0 over all (
ω
x
,
ω
y
)
∈
D
2 is comparable to the maximum obtained for blocks that cover the vehicles in motion.However, (29) suggests that the local phase change from multiple filter pairs centered at the same spatial position (
k
,
l
) can provide a constraint to robustly estimate motion and its direction. Given the block structure employed in the computation of the local phase, it is natural to utilize phase change information from multiple sources.Indeed, if, for a particular block(
k
,
l
), d
ϕ
k
l
/
d
t
(
ω
x
,
ω
y
,
t
)
=
-
v
x
(
t
)
ω
x
-
v
y
(
t
)
ω
y, then it is easy to see that d
ϕ
k
l
/
d
t
(
ω
x
,
ω
y
,
t
) will be zero on the line v
x
(
t
)
ω
x
+
v
y
(
t
)
ω
y
=
0 and have opposite sign on either side of this line. For example, in Figures 6(b) and 6(c), d
ϕ
k
l
/
d
t
(
ω
x
,
ω
y
,
t
0
) clearly exhibits this property for blocks that cover a vehicle in motion. The PMI is a tool to evaluate this property.Finally, the PMIs for all blocks are shown compactly in a heat map in Figure6(e). The figure shows clearly that the blocks corresponding to the two moving vehicles have a high PMI value while the stationary background areas have a low PMI value, allowing one to easily detect motion by employing simple thresholding (see also the plane “Radon transform and extract strength orientation of plane” in Figure 5). In addition, the orientation of motion in each block is readily observable even by inspection in Figure 6(b) by a line separating the yellow part and blue part in each block. Further results about the direction of motion are presented in Section 4.


## 3.4. Relationship to Biological Motion Detectors

A straightforward way to implement local motion detectors is to apply a complex-valued Gabor receptive field (5) to the video signal u, and then take the derivative of the phase with respect to time or apply a high-pass filter on the phase to approximate the derivative.We present here an alternate implementation without explicitly computing the phase. This will elucidate the relation between the phase-based motion detector presented in the previous section and some elementary motion detection models used in biology, such as the Reichardt motion detector [7] and motion energy detector [9].Assuming that the local phaseϕ
00
(
ω
x
,
ω
y
,
t
) is differentiable, we have (see also the Appendix)(42)
d
ϕ
00
ω
x
,
ω
y
,
t
d
t
=
d
b
ω
x
,
ω
y
,
t
/
d
t
a
ω
x
,
ω
y
,
t
-
d
a
ω
x
,
ω
y
,
t
/
d
t
b
ω
x
,
ω
y
,
t
a
ω
x
,
ω
y
,
t
2
+
b
ω
x
,
ω
y
,
t
2
,where a
(
ω
x
,
ω
y
,
t
) and b
(
ω
x
,
ω
y
,
t
) are, respectively, the real and imaginary parts of A
00
(
ω
x
,
ω
y
,
t
)
e
j
ϕ
00
(
ω
x
,
ω
y
,
t
).We notice that the denominator of (42) is the square of the local amplitude of u, and the numerator is of the form of a second order Volterra kernel. This suggests that the time derivative of the local phase can be viewed as a second order Volterra kernel that processes twonormalized spatially filtered inputs v
1 and v
2.We consider an elaborated Reichardt motion detector as shown in Figure7. It is equipped with a quadrature pair of Gabor filters whose outputs are r
1
(
t
)
=
a
(
ω
x
,
ω
y
,
t
) and r
2
(
t
)
=
b
(
ω
x
,
ω
y
,
t
), respectively, for a particular value of (
ω
x
,
ω
y
). The pair of Gabor filters that provide these outputs are the real and imaginary parts of w
(
x
,
y
)
e
-
j
(
ω
x
x
+
ω
y
y
). It also consists of a temporal high-pass g
1
(
t
) and temporal low-pass filter g
2
(
t
) [9]. The output of the elaborated Reichardt detector follows the diagram in Figure 7 and can be expressed as(43)
r
2
∗
g
1
t
r
1
∗
g
2
t
-
r
1
∗
g
1
t
r
2
∗
g
2
t
.The response can also be characterized by a second order Volterra kernel. We notice the striking similarity between (43) and the numerator of (42). In fact, the phase-based motion detector shares some properties with the Reichardt motion detector. For example, it is straightforward to see that a single phase-based motion detector is tuned to the temporal frequency of a moving sinusoidal grating.Figure 7
Diagram of an elaborated Reichardt motion detector. It consists of a quadrature pair of Gabor filters whose outputsr
1 and r
2 provide inputs to the high-pass filters g
1
(
t
) and the low-pass filters g
2
(
t
).Since the motion energy detector is formally equivalent to an elaborated Reichardt motion detector [9], the structure of the motion energy detector with divisive normalization is also similar to the phase-based motion detector.


## 4. Exploratory Results

In this section, we apply the phase-based motion detection algorithm on several video sequences and demonstrate its efficiency and effectiveness in detecting local motion. The motion detection algorithm is compared to two well-known biological motion detectors, namely, the Reichardt motion detector and the Barlow-Levick motion detector. We then show that the detected local motion can be used in motion segmentation tasks. The effectiveness of the segmentation is compared to segmentation using motion information obtained from a widely used optic flow based algorithm available in the literature [4].


### 4.1. Efficient Parallel Implementation

The algorithm was implemented in PyCUDA [23] and tested on an NVIDIA GeForce GTX TITAN GPU. All computations use single precision floating points. The processing speeds of the algorithm for several screen sizes are listed in Table 1. Clearly, the proposed phase-based motion detection algorithm has real-time capability to process video even with full High Definition screen size.Table 1
Processing Speeds of the proposed algorithm, the Reichardt motion detector, and the Barlow-Levick Motion detector.
Screen size
Processing capability (frames per second)
Proposed
Reichardt
Barlow-Levick
320
×
240
420
790
800
720
×
576
135
685
690
1280
×
720
60
274
293
1920
×
1080
27
191
194For comparison, we implemented the Reichardt motion detector and the Barlow-Levick motion detector. Their respective diagrams are shown in Figure8. Note that we moved the high-pass filter h
1
(
t
) to the front of the low-pass filter h
2
(
t
). This configuration provides a superior performance to the one in Figure 7. For both the Reichardt and the Barlow-Levick detectors, the videos are first blurred by a Gaussian filter (with the same variance as the Gaussian window in (30) used for the phase-based motion detector) and subsampled at the center of each overlapping block in the phase-based motion detector. The subsampled video then provides inputs to two 2D arrays of the circuits shown in Figure 8, one for the horizontal direction and one for the vertical direction. The outputs of the horizontal and vertical motion circuits form an array of motion vectors that indicate the strength and direction of motion. The three tested motion detectors have the same number of outputs as a result. The processing speeds of the Reichardt motion detector and the Barlow-Levick motion detector, both implemented in PyCUDA and tested on the same GPU, are shown in Table 1.Figure 8
(a) A Reichardt motion detector detecting motion between the two pointsu
k
l and u
k
+
1
,
l in the blurred and downsampled video. (b) A Barlow-Levick motion detector detecting motion between the two points u
k
l and u
k
+
1
,
l in the blurred and downsampled video. g
1
(
t
) denotes a high-pass filter and g
2
(
t
) denotes a low-pass filter.
(a)
(b)Note that the Reichardt motion detector and the Barlow-Levick motion detector are highly efficient due to the simplicity of their algorithms. The phase-based motion detection algorithm, however, is a much more sophisticated algorithm, and yet it can be implemented in real-time using parallel computing devices.The fast GPU implementation is based on the FFT and Matrix-Matrix multiplication. It is expected that those operations can be efficiently implemented in hardware, for example, FPGA.


### 4.2.  Examples of Phased-Based Motion Detection

We applied our motion detection algorithm on video sequences of the Change Detection 2014 dataset [22] that did not exhibit camera egomotion. For these video sequences, the standard deviation of the Gaussian window functions was set to 4 pixels and the block size was chosen to be 32
×
32 pixels. Threshold and normalization parameters were kept the same with the exception of the “thermal video” (in order to deal with larger background noise levels, see below). We also tested the same video sequences using the Reichardt motion detector and the Barlow-Levick motion detector. For the Reichardt detector, the high-pass filters were chosen as first order filters with a time constant of 200 milliseconds, and the low-pass filters were chosen as first order filters with a time constant of 300 milliseconds (assuming that the frame rate is 50 frames per second). Threshold was set to 2. For the Barlow-Levick motion detector, the time constant of the first high-pass filters was set to 250 milliseconds. The low-pass filters were the same as in the Reichardt detector. Threshold was set to 2.The first video was taken from a highway surveillance camera (“highway video”) under good illumination conditions and high contrast. The video had moderate noise, particularly on the road surface. The detected motion is shown in the top left panel of Figure9 (see Supplementary Video S2 for full video). The phase-based motion detection algorithm captured both the moving cars and the tree leaves moving (due to the wind). In the time interval between the 9th and 10th second, the camera was slightly moved left- and right-wards within 5 frames, again possibly due to the wind. Movement due to this shift was captured by the motion detection algorithm and the algorithm was fast enough to correctly determine the direction of this movement. In this video, we already noted that this algorithm suffers from the aperture problem. For example, in front of the van where a long, horizontal edge is present, the detected motion is mostly pointing downwards. In addition to moving downwards, the edge is also moving to the left, however. This is expected since the algorithm only detects motion locally and does not take into account the overall shape of any object.Figure 9
Top: motion detection of the “highway video” using the phase-based algorithm (a), the Reichardt motion detector (b), and the Barlow-Levick motion detector (c). Red arrows indicate detected motion and its direction. Bottom: the contrast of the video was artificially reduced by 5-fold and the mean was reduced to3
/
5 of the original. The red arrows are duplicated from the motion detection result on the original video as in top. Blue arrows are the result of motion detection on the video with reduced contrast. If motion is detected both in the original video and in the video with reduced contrast, then the arrow is shown in magenta (as a mix of blue and red) (see Supplementary Video S2 for full video).
(a)
(b)
(c)For comparison, the motion detection results for the Reichardt motion detector and the Barlow-Levick motion detector are shown in the top middle and top right panel of Figure9 (see Supplementary Video S2 for full video). The Reichardt detector performed relatively well when vehicles are moving faster, but the direction it predicts for vehicles moving slower, for example, on the back of the image is not accurate. In addition, motion is still detected in some parts of the screen where the vehicles have just passed by. The detection result for the Barlow-Levick motion detector was poorer. In particular, the response to OFF edge movement is always the opposite to the actual movement direction.We then squeezed the range of the screen intensity from[
0,1
] to [
0.2,0.4
], resulting in a video with a lower mean luminance and lower contrast. The motion detection results on the low-contrast video are shown in the bottom 3 panels of Figure 9 (see Supplementary Video S2 for full video). For reference, the motion detection results for the original video are shown in red arrows. Motion detected in the low-contrast video is shown in blue arrows if no motion is detected for the block in the original video. If motion is detected in both the original and the low-contrast video, the arrows are shown in magenta. Figure 9 clearly shows that while the motion detection performance is degraded for all three detectors, the phase-based motion detector performed still quite well in the low-contrast video, detecting most of the moving vehicles. The other two detectors missed many of the blocks where motion was detected in the original movie.To quantify how well each detector works under different contrast conditions, we computed the ratio of unthresholded output values of each motion detector between lower contrast and full contrast video. For example, in the case of phase-based motion detector, we computed, for each block, the ratio between the PMI index for the lower contrast video and that for the full contrast video. The ratios are then averaged across all blocks where motion is detected in the full contrast video. This average for different contrasts is shown in Figure10 by the blue curve. In the ideal case when the normalization constant ϵ is 0, the phase detector should produce invariant responses to videos with different contrast. The curve shown here is mainly due to a nonzero ϵ. The ratios for the Reichardt motion detector and the Barlow-Levick motion detector are computed similarly and are shown in red and yellow, respectively. It is clear that the phase-based motion detector has a superior performance across a range of contrast values. At 20
% contrast, the phase-based detector still has 50
% of the PMI index value for full contrast video. As expected, the response of Barlow-Levick motion detector is linear with respect to contrast, and the Reichardt motion detector has a quadratic relation to contrast. For a fixed threshold value a larger ratio equates to a more consistent performance at lower contrast.Figure 10
The ratio for the phase-based motion detector is defined as follows: we set the motion detected in the100
% contrast case as the baseline; we compute the ratio between the PMI index for the lower contrast video and that for the full contrast video for each block; the ratios are averaged across all the blocks where motion is detected in the full contrast video; the average is given as the ratio for the phase-based motion detector. The ratio for the Reichardt motion detector and the Barlow-Levick motion detector is similarly defined.The second video was captured by a surveillance camera in a train station (“train station video”). The video was under moderate room light with a low noise level. The front side of the video had high contrast; illumination on the back side was quite low. The detected motion is shown in the video of Figure11 (see Supplementary Video S3 for full video). Movements of people were successfully captured by the motion detection algorithm.Figure 11
Motion detection of the “train station video” using the phase-based algorithm (a), the Reichardt motion detector (b), and the Barlow-Levick motion detector (c) (see Figure9 for description of each panel) (see Supplementary Video S3 for full video).
(a)
(b)
(c)The third video was a “thermal video” with large amount of background noise (thermal video). The threshold for detecting motion was raised by 60% in order to mitigate the increased level of noise. The detected motion is shown in the video of Figure12 (see Supplementary Video S4 for full video).Figure 12
Motion detection of the “thermal video” using the phase-based algorithm (a), the Reichardt motion detector (b), and the Barlow-Levick motion detector (c) (see Figure9 for description of each panel) (see Supplementary Video S4 for full video).
(a)
(b)
(c)The last example we show here was taken from a highway surveillance camera at night (“winterstreet video”). The overall illumination on the lower-left side was low whereas illumination was moderate on the upper-right side where the road was covered by snow. The detected motions are shown in the video in Figure13 (see Supplementary Video S5 for full video). We note that, overall, car movements were successfully detected. Car movements on the lower-left side, however, suffered from low illumination and some parts of the car were not detected well due to the trade-off employed for noise suppression.Figure 13
Motion detection of the “winterstreet video” using the phase-based algorithm (a), the Reichardt motion detector (b), and the Barlow-Levick motion detector (c) (see Figure9 for description of each panel) (see Supplementary Video S5 for full video).
(a)
(b)
(c)With a higher threshold, the phase-base motion detection algorithm is able to detect motion under noisy conditions. We added to the original “highway video” and “train station video” Gaussian white noise with standard deviation5
% of the maximum luminance range. The results are shown, respectively, in Figures 14 and 15 (see Supplementary Videos S6 and S7 for full videos).Figure 14
Phase-based motion detection applied on the “highway video” with added Gaussian white noise (see Supplementary Video S6 for full video).Figure 15
Phase-based motion detection applied on the “train station video” with added Gaussian white noise (see Supplementary Video S7 for full video).


### 4.3. Examples of Motion Segmentation

We asked whether the detected motion signals in the video sequences can be useful for segmenting moving objects from the background. We applied a larger threshold to only signal motion for salient objects. The32
×
32 blocks, however, introduce large boundaries around the moving objects. To reduce the boundary and to segment the moving object more closely to the actual object boundary, we applied the motion detection algorithm around the detected boundary with 16
×
16 blocks. If 16
×
16 blocks did not indicate motion, then the corresponding area was removed from the segmented object area.For comparison, we performed motion segmentation based on local motion detection based on an optic flow algorithm [4]. The segmentation in this case was implemented by comparing the length of the optic flow vectors with an appropriate threshold.We employed a simple thresholding for both phase-based motion detection algorithms and optic flow based motion detection. More sophisticated algorithms may produce better segmentation results. Thus, the segmentation was purely based on motion cues and no postprocessing at pixel level was performed. The state-of-the-art results for the Change Detection 2014 dataset utilize multiple cues such as motion, color, and background extraction to segment objects and thereby achieve better results. We are exploring here, however, only the case where motion is the only cue for segmentation. Therefore, the ground truth information from the dataset was not applicable to our test. We will, therefore, only show the effectiveness of motion segmentation visually.We first applied motion based segmentation on 2-second segment of the “highway video.” The result using the local phase-based motion detection algorithm is shown in the video of Figure16(a) (see Supplementary Video S8a for full video) and that using optic flow based motion detection algorithm is shown in the video of Figure 16(b) (see Supplementary Video S8b for full video). Both videos are played back at 1
/
4 of speed. With a higher threshold, the movement of the leaves was no longer picked up by the phase-based motion detector. Therefore, only the cars were identified as moving objects and they are indicated in red. Although the moving objects were not perfectly segmented on their boundary, they were mostly captured. For the optic flow based segmentation, since the regions of interest are set by thresholding the length of the velocity vector, objects moving at lower speed, for example, the cars on the top, were not always picked up by the segmentation algorithm.Figure 16
Segmentation of moving cars based on detected motion in the “highway video.” (a) Motion detected by the phase-based motion detection algorithm. (b) Motion detected using optic flow algorithm [4] (see Supplementary Videos S8a and S8b for full videos).
(a)
(b)We then applied the motion segmentation to a 2-second segment of the “train station video,” the “thermal video,” and the “winterstreet video.” The results are shown, respectively, in the videos of Figures17, 18, and 19 (see Supplementary Videos S9a and S9b, S10a and S10b, and S11a and S11b, resp., for full videos).Figure 17
Segmentation of moving people in the “train station video” based on detected motion. (a) Motion detected by the phase-based motion detection algorithm. (b) Motion detected using optic flow algorithm [4] (see Supplementary Videos S9a and S9b for full videos).
(a)
(b)Figure 18
Segmentation of moving people in the “thermal video” based on detected motion. (a) Motion detected by the phase-based motion detection algorithm. (b) Motion detected using optic flow algorithm [4] (see Supplementary Videos S10a and S10b for full videos).
(a)
(b)Figure 19
Segmentation of moving cars in the “winter street video” based on detected motion. (a) Motion detected by the phase-based motion detection algorithm. (b) Motion detected using optic flow algorithm [4] (see Supplementary Videos S11a and S11b for full videos).
(a)
(b)These results show that for the purpose of detecting local motion and its use as a motion segmentation cue, the local phase-based motion detector works as good, if not better than a simple thresholding segmentation using an optic flow based algorithm.


## 4.1. Efficient Parallel Implementation

The algorithm was implemented in PyCUDA [23] and tested on an NVIDIA GeForce GTX TITAN GPU. All computations use single precision floating points. The processing speeds of the algorithm for several screen sizes are listed in Table 1. Clearly, the proposed phase-based motion detection algorithm has real-time capability to process video even with full High Definition screen size.Table 1
Processing Speeds of the proposed algorithm, the Reichardt motion detector, and the Barlow-Levick Motion detector.
Screen size
Processing capability (frames per second)
Proposed
Reichardt
Barlow-Levick
320
×
240
420
790
800
720
×
576
135
685
690
1280
×
720
60
274
293
1920
×
1080
27
191
194For comparison, we implemented the Reichardt motion detector and the Barlow-Levick motion detector. Their respective diagrams are shown in Figure8. Note that we moved the high-pass filter h
1
(
t
) to the front of the low-pass filter h
2
(
t
). This configuration provides a superior performance to the one in Figure 7. For both the Reichardt and the Barlow-Levick detectors, the videos are first blurred by a Gaussian filter (with the same variance as the Gaussian window in (30) used for the phase-based motion detector) and subsampled at the center of each overlapping block in the phase-based motion detector. The subsampled video then provides inputs to two 2D arrays of the circuits shown in Figure 8, one for the horizontal direction and one for the vertical direction. The outputs of the horizontal and vertical motion circuits form an array of motion vectors that indicate the strength and direction of motion. The three tested motion detectors have the same number of outputs as a result. The processing speeds of the Reichardt motion detector and the Barlow-Levick motion detector, both implemented in PyCUDA and tested on the same GPU, are shown in Table 1.Figure 8
(a) A Reichardt motion detector detecting motion between the two pointsu
k
l and u
k
+
1
,
l in the blurred and downsampled video. (b) A Barlow-Levick motion detector detecting motion between the two points u
k
l and u
k
+
1
,
l in the blurred and downsampled video. g
1
(
t
) denotes a high-pass filter and g
2
(
t
) denotes a low-pass filter.
(a)
(b)Note that the Reichardt motion detector and the Barlow-Levick motion detector are highly efficient due to the simplicity of their algorithms. The phase-based motion detection algorithm, however, is a much more sophisticated algorithm, and yet it can be implemented in real-time using parallel computing devices.The fast GPU implementation is based on the FFT and Matrix-Matrix multiplication. It is expected that those operations can be efficiently implemented in hardware, for example, FPGA.


## 4.2.  Examples of Phased-Based Motion Detection

We applied our motion detection algorithm on video sequences of the Change Detection 2014 dataset [22] that did not exhibit camera egomotion. For these video sequences, the standard deviation of the Gaussian window functions was set to 4 pixels and the block size was chosen to be 32
×
32 pixels. Threshold and normalization parameters were kept the same with the exception of the “thermal video” (in order to deal with larger background noise levels, see below). We also tested the same video sequences using the Reichardt motion detector and the Barlow-Levick motion detector. For the Reichardt detector, the high-pass filters were chosen as first order filters with a time constant of 200 milliseconds, and the low-pass filters were chosen as first order filters with a time constant of 300 milliseconds (assuming that the frame rate is 50 frames per second). Threshold was set to 2. For the Barlow-Levick motion detector, the time constant of the first high-pass filters was set to 250 milliseconds. The low-pass filters were the same as in the Reichardt detector. Threshold was set to 2.The first video was taken from a highway surveillance camera (“highway video”) under good illumination conditions and high contrast. The video had moderate noise, particularly on the road surface. The detected motion is shown in the top left panel of Figure9 (see Supplementary Video S2 for full video). The phase-based motion detection algorithm captured both the moving cars and the tree leaves moving (due to the wind). In the time interval between the 9th and 10th second, the camera was slightly moved left- and right-wards within 5 frames, again possibly due to the wind. Movement due to this shift was captured by the motion detection algorithm and the algorithm was fast enough to correctly determine the direction of this movement. In this video, we already noted that this algorithm suffers from the aperture problem. For example, in front of the van where a long, horizontal edge is present, the detected motion is mostly pointing downwards. In addition to moving downwards, the edge is also moving to the left, however. This is expected since the algorithm only detects motion locally and does not take into account the overall shape of any object.Figure 9
Top: motion detection of the “highway video” using the phase-based algorithm (a), the Reichardt motion detector (b), and the Barlow-Levick motion detector (c). Red arrows indicate detected motion and its direction. Bottom: the contrast of the video was artificially reduced by 5-fold and the mean was reduced to3
/
5 of the original. The red arrows are duplicated from the motion detection result on the original video as in top. Blue arrows are the result of motion detection on the video with reduced contrast. If motion is detected both in the original video and in the video with reduced contrast, then the arrow is shown in magenta (as a mix of blue and red) (see Supplementary Video S2 for full video).
(a)
(b)
(c)For comparison, the motion detection results for the Reichardt motion detector and the Barlow-Levick motion detector are shown in the top middle and top right panel of Figure9 (see Supplementary Video S2 for full video). The Reichardt detector performed relatively well when vehicles are moving faster, but the direction it predicts for vehicles moving slower, for example, on the back of the image is not accurate. In addition, motion is still detected in some parts of the screen where the vehicles have just passed by. The detection result for the Barlow-Levick motion detector was poorer. In particular, the response to OFF edge movement is always the opposite to the actual movement direction.We then squeezed the range of the screen intensity from[
0,1
] to [
0.2,0.4
], resulting in a video with a lower mean luminance and lower contrast. The motion detection results on the low-contrast video are shown in the bottom 3 panels of Figure 9 (see Supplementary Video S2 for full video). For reference, the motion detection results for the original video are shown in red arrows. Motion detected in the low-contrast video is shown in blue arrows if no motion is detected for the block in the original video. If motion is detected in both the original and the low-contrast video, the arrows are shown in magenta. Figure 9 clearly shows that while the motion detection performance is degraded for all three detectors, the phase-based motion detector performed still quite well in the low-contrast video, detecting most of the moving vehicles. The other two detectors missed many of the blocks where motion was detected in the original movie.To quantify how well each detector works under different contrast conditions, we computed the ratio of unthresholded output values of each motion detector between lower contrast and full contrast video. For example, in the case of phase-based motion detector, we computed, for each block, the ratio between the PMI index for the lower contrast video and that for the full contrast video. The ratios are then averaged across all blocks where motion is detected in the full contrast video. This average for different contrasts is shown in Figure10 by the blue curve. In the ideal case when the normalization constant ϵ is 0, the phase detector should produce invariant responses to videos with different contrast. The curve shown here is mainly due to a nonzero ϵ. The ratios for the Reichardt motion detector and the Barlow-Levick motion detector are computed similarly and are shown in red and yellow, respectively. It is clear that the phase-based motion detector has a superior performance across a range of contrast values. At 20
% contrast, the phase-based detector still has 50
% of the PMI index value for full contrast video. As expected, the response of Barlow-Levick motion detector is linear with respect to contrast, and the Reichardt motion detector has a quadratic relation to contrast. For a fixed threshold value a larger ratio equates to a more consistent performance at lower contrast.Figure 10
The ratio for the phase-based motion detector is defined as follows: we set the motion detected in the100
% contrast case as the baseline; we compute the ratio between the PMI index for the lower contrast video and that for the full contrast video for each block; the ratios are averaged across all the blocks where motion is detected in the full contrast video; the average is given as the ratio for the phase-based motion detector. The ratio for the Reichardt motion detector and the Barlow-Levick motion detector is similarly defined.The second video was captured by a surveillance camera in a train station (“train station video”). The video was under moderate room light with a low noise level. The front side of the video had high contrast; illumination on the back side was quite low. The detected motion is shown in the video of Figure11 (see Supplementary Video S3 for full video). Movements of people were successfully captured by the motion detection algorithm.Figure 11
Motion detection of the “train station video” using the phase-based algorithm (a), the Reichardt motion detector (b), and the Barlow-Levick motion detector (c) (see Figure9 for description of each panel) (see Supplementary Video S3 for full video).
(a)
(b)
(c)The third video was a “thermal video” with large amount of background noise (thermal video). The threshold for detecting motion was raised by 60% in order to mitigate the increased level of noise. The detected motion is shown in the video of Figure12 (see Supplementary Video S4 for full video).Figure 12
Motion detection of the “thermal video” using the phase-based algorithm (a), the Reichardt motion detector (b), and the Barlow-Levick motion detector (c) (see Figure9 for description of each panel) (see Supplementary Video S4 for full video).
(a)
(b)
(c)The last example we show here was taken from a highway surveillance camera at night (“winterstreet video”). The overall illumination on the lower-left side was low whereas illumination was moderate on the upper-right side where the road was covered by snow. The detected motions are shown in the video in Figure13 (see Supplementary Video S5 for full video). We note that, overall, car movements were successfully detected. Car movements on the lower-left side, however, suffered from low illumination and some parts of the car were not detected well due to the trade-off employed for noise suppression.Figure 13
Motion detection of the “winterstreet video” using the phase-based algorithm (a), the Reichardt motion detector (b), and the Barlow-Levick motion detector (c) (see Figure9 for description of each panel) (see Supplementary Video S5 for full video).
(a)
(b)
(c)With a higher threshold, the phase-base motion detection algorithm is able to detect motion under noisy conditions. We added to the original “highway video” and “train station video” Gaussian white noise with standard deviation5
% of the maximum luminance range. The results are shown, respectively, in Figures 14 and 15 (see Supplementary Videos S6 and S7 for full videos).Figure 14
Phase-based motion detection applied on the “highway video” with added Gaussian white noise (see Supplementary Video S6 for full video).Figure 15
Phase-based motion detection applied on the “train station video” with added Gaussian white noise (see Supplementary Video S7 for full video).


## 4.3. Examples of Motion Segmentation

We asked whether the detected motion signals in the video sequences can be useful for segmenting moving objects from the background. We applied a larger threshold to only signal motion for salient objects. The32
×
32 blocks, however, introduce large boundaries around the moving objects. To reduce the boundary and to segment the moving object more closely to the actual object boundary, we applied the motion detection algorithm around the detected boundary with 16
×
16 blocks. If 16
×
16 blocks did not indicate motion, then the corresponding area was removed from the segmented object area.For comparison, we performed motion segmentation based on local motion detection based on an optic flow algorithm [4]. The segmentation in this case was implemented by comparing the length of the optic flow vectors with an appropriate threshold.We employed a simple thresholding for both phase-based motion detection algorithms and optic flow based motion detection. More sophisticated algorithms may produce better segmentation results. Thus, the segmentation was purely based on motion cues and no postprocessing at pixel level was performed. The state-of-the-art results for the Change Detection 2014 dataset utilize multiple cues such as motion, color, and background extraction to segment objects and thereby achieve better results. We are exploring here, however, only the case where motion is the only cue for segmentation. Therefore, the ground truth information from the dataset was not applicable to our test. We will, therefore, only show the effectiveness of motion segmentation visually.We first applied motion based segmentation on 2-second segment of the “highway video.” The result using the local phase-based motion detection algorithm is shown in the video of Figure16(a) (see Supplementary Video S8a for full video) and that using optic flow based motion detection algorithm is shown in the video of Figure 16(b) (see Supplementary Video S8b for full video). Both videos are played back at 1
/
4 of speed. With a higher threshold, the movement of the leaves was no longer picked up by the phase-based motion detector. Therefore, only the cars were identified as moving objects and they are indicated in red. Although the moving objects were not perfectly segmented on their boundary, they were mostly captured. For the optic flow based segmentation, since the regions of interest are set by thresholding the length of the velocity vector, objects moving at lower speed, for example, the cars on the top, were not always picked up by the segmentation algorithm.Figure 16
Segmentation of moving cars based on detected motion in the “highway video.” (a) Motion detected by the phase-based motion detection algorithm. (b) Motion detected using optic flow algorithm [4] (see Supplementary Videos S8a and S8b for full videos).
(a)
(b)We then applied the motion segmentation to a 2-second segment of the “train station video,” the “thermal video,” and the “winterstreet video.” The results are shown, respectively, in the videos of Figures17, 18, and 19 (see Supplementary Videos S9a and S9b, S10a and S10b, and S11a and S11b, resp., for full videos).Figure 17
Segmentation of moving people in the “train station video” based on detected motion. (a) Motion detected by the phase-based motion detection algorithm. (b) Motion detected using optic flow algorithm [4] (see Supplementary Videos S9a and S9b for full videos).
(a)
(b)Figure 18
Segmentation of moving people in the “thermal video” based on detected motion. (a) Motion detected by the phase-based motion detection algorithm. (b) Motion detected using optic flow algorithm [4] (see Supplementary Videos S10a and S10b for full videos).
(a)
(b)Figure 19
Segmentation of moving cars in the “winter street video” based on detected motion. (a) Motion detected by the phase-based motion detection algorithm. (b) Motion detected using optic flow algorithm [4] (see Supplementary Videos S11a and S11b for full videos).
(a)
(b)These results show that for the purpose of detecting local motion and its use as a motion segmentation cue, the local phase-based motion detector works as good, if not better than a simple thresholding segmentation using an optic flow based algorithm.


## 5. Discussion

Previous research demonstrated thatglobal phase information alone can be used to faithfully represent visual scenes. Here we provided a reconstruction algorithm of visual scenes by only usinglocal phase information. More importantly, local phase information can be effectively used to detect local motion. Through a simple temporal derivative of the phase, we obtained a second order Volterra kernel that is applied on two normalized inputs. The structure of the second order Volterra kernel in the phase-based motion detector is akin to models employed to detect motion in biological vision, for example, the Reichardt detector [7] and the motion energy detector [9].We then proposed an efficient, FFT-based algorithm employing the change in local phase for detecting motion. In order to exploit the special structure of the change in phase in the frequency domain that is due to rigid motion, the phase-based motion detection algorithm also incorporates the Radon transform, a transform closely related to the Fourier transform. Based on the Radon transform, a motion indicator was proposed to robustly detect whether the phase change is caused by motion. Therefore, the algorithm can be efficiently implemented whenever the FFT is available/supported. We showed examples of applying the phase-based motion detection algorithm on several video sequences. We also showed that the locally detected motion can be used for segmenting moving objects in video scenes. We compared the segmentation of moving objects using our local phase-based algorithm to segmentation achieved using a widely used optic flow based algorithm. Our results suggest that spatial phase information may provide an efficient alternative to perform many visual tasksin silico as well as in modelingin vivo biological vision systems. This is consistent with other recent findings [15].Note that phase information has been used for solving various visual tasks in the past. In fact, phase has been successfully employed in optic flow algorithms [24] and image registration for translation [25], both applied to motion related tasks. The phase-based optic flow algorithm applies the optic flow equation on the local phase of images rather than the intensity itself to achieve better resolution and robustness in estimating motion velocity. The phase correlation method computes the normalized cross-power spectrum of two images to extract the phase difference. It provides better accuracy and robustness as compared to the classical cross correlation approach applied to two consecutive images in a sequence. However, it has limitations when dealing with images with a repetitive structure.Our method differs from the above two cases in the following ways. First, it employs a simple temporal derivative/high-pass filtering on the phase to extract local phase changes. The structure of thephase change is exploited for better detection of motion. On the contrary, the phase correlation method considers the structure of thephase itself. This also allowed us to make the motion detection a more local, continuous process rather than purely operating globally on discrete frames. Second, it explores the structure of the change of phase due to motion in the frequency domain rather than in the spatial domain, in which the key constraints of optic flow equations are based. Third, instead of focusing on estimating the exact velocity or shift, our method is centered on the detection of motion with a coarse local estimate of direction. This is the case in the first steps of biological motion detection in the retina or optic lobe of insects; we discussed the resemblance of the method presented here to those of biological models of motion detection.The proposed motion detection algorithm, however, shares several advantages with the other phase-based methods. For example, it is sensitive to motion that only induces a subpixel shift between frames and for very small differences in intensity. In addition, when compared to the amplitude, the local phase is robust under different contrast and illumination conditions. Consequently, the algorithm presented here can operate in a wide range of contrast/illumination conditions.Furthermore, once the local phase is extracted from each block, motion detection becomes a localized, temporal operation on the local phase of each block. This forms the basis for the highly parallel structure in the phase-based motion detection algorithm. By contrast, traditional optic flow techniques often rely on explicit comparisons across spatial locations, which increase, for example, memory complexity since all states must be made available to their neighbors.We also notice that, for each block, a large number of measurements of phase changes are obtained. From Figure6, we see that this number is much higher than that of the original pixel space. In other words, in order to detect motion in a robust way, the local phase-based motion detection algorithm undergoes an expansion of measurements before settling down onto a single motion indicator value. This number of measurements, however, does not incur additional computational demand thanks to the highly efficient FFT algorithm. Biological visual systems often have a similar structure. For example, in the vertebrate retina, computations are carried out by an extraordinarily large number of neurons until measurements of the visual scene are projected by a fraction of the neurons onto the cortex [26]. Similar expansion takes place in the primary visual cortex as well.We also highlight the ease of implementation of the intrinsically parallel algorithm proposed here. The algorithm introduced in Section3.3 is based on the FFT algorithm and does not require solving an optimization problem. It can be efficiently implemented in hardware, for example, FPGAs, or in software, for example, using GPU accelerators. We note that extending the FFT to higher dimensions is straightforward and the implementations of higher dimensional FFTs are also highly efficient. Clearly, the methodology can be applied to motion detection of data in 3D or higher dimensional space, where the Radon transform operates over planes or hyperplanes.Finally, we argue that the change of phase, although highly nonlinear, can be obtained through a normalization (gain control) followed by a second order Volterra kernel. This separation of a higher order nonlinearity into gain control block and a lower order nonlinear filter can be used for modeling motion detection circuits in biological systems.



---

*Source: 7915245-2016-01-10.xml*