(Quasi)Periodicity Quantification in Video Data, Using Topology

TL;DR

This paper presents a new topology-based framework for quantifying periodic and quasiperiodic dynamics in videos without segmentation or training, demonstrating robustness and novel applications like detecting biphonation.

Contribution

It introduces a topology-driven method to measure recurrence in videos directly, without segmentation or object tracking, advancing the analysis of complex recurrent behaviors.

Findings

Scores are robust against noise levels.

Periodicity score outperforms existing methods.

Successfully detects biphonation in vocal fold videos.

Abstract

This work introduces a novel framework for quantifying the presence and strength of recurrent dynamics in video data. Specifically, we provide continuous measures of periodicity (perfect repetition) and quasiperiodicity (superposition of periodic modes with non-commensurate periods), in a way which does not require segmentation, training, object tracking or 1-dimensional surrogate signals. Our methodology operates directly on video data. The approach combines ideas from nonlinear time series analysis (delay embeddings) and computational topology (persistent homology), by translating the problem of finding recurrent dynamics in video data, into the problem of determining the circularity or toroidality of an associated geometric space. Through extensive testing, we show the robustness of our scores with respect to several noise models/levels, we show that our periodicity score is superior to other methods when compared to human-generated periodicity rankings, and furthermore, we show that our quasiperiodicity score clearly indicates the presence of biphonation in videos of vibrating vocal folds, which has never before been accomplished end to end quantitatively.

Figures

Tables

Table 1. Table 1: AUROC values for different levels of noise, from the binary classification task: periodic (bird flapping, heart beating, pendulum) vs. non-recurrent (driving - left subcell, explosions - right subcell) based on the periodicity score (Equation  19 ). Also two synthetic quasiperiodic videos (sideways disks, modulated pulses) are compared to the same two non-recurrent videos based on the quasiperiodicity score (Equation  20 ).

	Blur 20		Blur 40		Blur 80		Awgn $\sigma =1$		Awgn $\sigma =2$		Awgn $\sigma =3$		Bit Err 5%		Bit Err 10%		Bit Err 15%
Bird Flapping	1	1	1	1	1	1	1	1	1	1	1	1	1	0.97	1	0.92	0.75	0.79
Heart Beat	1	1	1	1	0.91	0.94	1	1	1	1	0.9	0.89	1	0.91	0.98	0.87	0.56	0.62
Pendulum	1	0.94	0.84	0.85	0.43	0.68	1	1	1	1	1	1	1	1	1	0.99	0.98	0.97
QuasiPeriodic Disks	1	0.85	1	0.83	0.75	0.39	1	1	1	1	1	1	1	0.92	0.99	0.93	0.82	0.82
QuasiPeriodic Pulses	1	0.9	1	0.87	1	0.83	1	1	1	1	1	1	1	0.96	0.99	0.95	0.95	0.95

Table 2. Table 2: The Kendall τ 𝜏 \tau scores between all of the machine rankings and the Hodge aggregated human rankings.

	Human	SW+TDA	Freq[6]	CDscore[6]	Clarity[7]
Human	1	0.663	-0.295	0.347	0.284
SW+TDA	0.663	1	-0.316	0.221	0.516
Freq[6]	-0.295	-0.316	1	-0.0842	-0.189
CDscore[6]	0.347	0.221	-0.0842	1	0.411
Clarity[7]	0.284	0.516	-0.189	0.411	1

Table 3. Table 3: Average runtimes, in milliseconds, per video for all of the algorithms

SW+TDA	Freq[6]	CDscore[6]	Clarity[7]
3101ms	73ms	176ms	154ms

Table 4. Table 4: Results of our sliding window pipeline on videos of periodic vocal folds, biphonation, and irregularities. We give the max persistence periodicity score (PS), the modified periodicity score (MPS), the harmonic score (HS), and quasiperiodic score (QPS) presented in Section  3.4 . We also show the window size (Win) that the autocorrelation technique in Section  3.6 gives. We have bolded the top three MPS and QPS scores across all videos. The max modified periodic scores include the two periodic videos and one of the biphonation videos. The max quasiperiodic scores are all of the biphonation videos, which means the one with a high periodicity score could be ruled out of the periodicity category.

Video Name	Win	PS	MPS	QPS
Periodic 1 ([13])	16	0.816	0.789	0.011
Periodic 2 ([26], Figure 21)	32	0.601	0.533	0.009
Biphonation 1 ([28])	53	0.638	0.294	0.292
Biphonation 2 ([28])	42	0.703	0.583	0.116
Biphonation 3 ([28], Figure 22)	67	0.515	0.076	0.426
Mucus Perturbed Periodic ([49])	94	0.028	0.019	0.004
Irregular ([13], Figure 23)	232	0.18	0.097	0.04

Equations

$$\begin{array}[]{rccl}X:\mathbb{R}^{+}\longrightarrow\mathbb{R}^{W\times H}\end{array}$$

$$SW_{d,\tau}X(t)=\left[\begin{array}[]{c}X(t)\\ X(t+\tau)\\ \vdots\\ X(t+d\tau)\end{array}\right]\in\mathbb{R}^{W\times H\times(d+1)}$$

$$f_h(t)=\cos \left(\frac{\pi }{5}t\right)+\cos \left(\frac{\pi }{15}t\right)$$

$$f_q(t)=\cos \left(\frac{\pi }{5}t\right)+\cos \left(\frac{1}{5}t\right)$$

$${\mathbb{T}}^N=\underset{N-times}{\underbrace{S^1\times \cdots \times S^1}},S^1=\{z\in \mathbb{C}:|z|=1\}$$

$$X_i(t)=\sum _{n=1}^N{a}_{in}\cos ({\omega }_{n}t)$$

$$X(t)=\sum _{n=1}^N{A}^{n}cos({\omega }_{n}t)$$

$$SW_{d,\tau}X(t)=\sum_{n=1}^{N}\left[\begin{array}[]{c}A^{n}\cos(\omega_{n}t)\\ \vdots\\ A^{n}\cos(\omega_{n}(t+d\tau)\end{array}\right]$$

$$S{W}_{d,\tau }X(t)=\sum _{n=1}^N{\vec{u}}_n\cos ({\omega }_{n}t)-{\vec{v}}_n\sin ({\omega }_{n}t)$$

$$f_{i}(t)=\left\{\begin{array}[]{cc}X_{i}(t)&0\leq t\leq\ell\\ 0&\text{otherwise}\end{array}\right\}$$

$$X_i(t)=\sum _{m=-\infty }^{\infty }{f}_i(t-m\ell )$$

$$\mathcal{F}\left\{X_i(t)\right\}(k)\propto \sum _{m=-\infty }^{\infty }\mathcal{F}(f_i(t))\left(\frac{m}{\ell }\right)\delta \left(\frac{m}{\ell }-k\right)$$

$$X_i(t)\propto \sum _{m=-\infty }^{\infty }\mathcal{F}\left(f_i(t)\right)\left(\frac{m}{\ell }\right)e^{i\frac{2\pi m}{\ell }t}$$

$$\sin \left(\frac{2\pi }{\ell }t\right)+\frac{1}{3}\sin \left(\frac{6\pi }{\ell }t\right)+\frac{1}{5}\sin \left(\frac{10\pi }{\ell }t\right)+\dots$$

$$K=\big{\{}\{a\},\{b\},\{c\},\{a,b\},\{b,c\},\{a,c\}\big{\}}$$

$$R_{\epsilon}(\mathbb{X}):=\big{\{}\sigma\subset\mathbb{X}:d_{\mathbb{X}}(x_{1},x_{2})\leq\epsilon,\;\forall x_{i},x_{j}\in\sigma\big{\}}$$

$$\hat{A}=U^{T}A=U^{T}USV=SV$$

$$SW_{d,\tau}(t)=\left[\begin{array}[]{c}U^{T}X(t)\\ \vdots\\ U^{T}X(t+d\tau)\end{array}\right]$$

$$\Vert S{W}_{d,\tau }(t)-S{W}_{d,\tau }(t^{\prime })\Vert =\Vert S{W}_{d,\tau }X(t)-S{W}_{d,\tau }X(t^{\prime })\Vert$$

$$||S{W}_{d,1}X(i)-S{W}_{d,1}X(j)|{|}_2^2=\sum _{m=0}^d||X(i+m)-X(j+m)|{|}_2^2$$

$${\widetilde{SW}}_{d,\tau }(t)=\frac{S{W}_{d,\tau }(t)-(S{W}_{d,\tau }(t{)}^{T}1)1}{||S{W}_{d,\tau }(t)-(S{W}_{d,\tau }(t{)}^{T}1)1|{|}_2}$$

$${\hat{X}}_i(t)=X_i(t)\ast -at{\exp }^{-t^2/(2{\sigma }^2)}$$

$$m{p}_1({\mathsf{dgm}}_n)=\max \{d-b:(b,d)\in {\mathsf{dgm}}_n\}$$

$$PS=\frac{1}{\sqrt{3}}m{p}_1({\mathsf{dgm}}_1)$$

$$QPS=\sqrt{\frac{m{p}_2({\mathsf{dgm}}_1)m{p}_1({\mathsf{dgm}}_2)}{3}}$$

$$MPS=\frac{1}{\sqrt{3}}\left(m{p}_1({\mathsf{dgm}}_1)-m{p}_2({\mathsf{dgm}}_1)\right)$$

$$d\tau =\frac{\pi k}{L}\left(\frac{d}{d+1}\right)$$

$$r_t(\tau )=\sum _{j=t}^{t+N-1-\tau }{x}_j{x}_{j+\tau }$$

$$d_t(\tau )=\sum _{j=t}^{t+N-1-\tau }(x_j-x_{j+\tau }{)}^2$$

$$m_t(\tau )=\sum _{j=t}^{t+N-1-\tau }(x_j^2+x_{j+\tau }^2)$$

Full text

1 Introduction

Periodicity characterizes many natural motions including animal locomotion (walking/wing flapping/slithering), spinning wheels, oscillating pendulums, etc. Quasiperiodicity, thought of as the superposition of non-commensurate frequencies, occurs naturally during transitions from ordinary to chaotic dynamics [11]. The goal of this work is to automate the analysis of videos capturing periodic and quasiperiodic motion. In order to identify both classes of motion in a unified framework, we generalize 1-dimensional (1D) sliding window embeddings [38] to reconstruct periodic and quasiperiodic attractors from videos111Some of the analysis and results appeared as part of the Ph.D. thesis of the first author [tralie2017geometric].222Code to replicate results: https://github.com/ctralie/SlidingWindowVideoTDA333Supplementary material and videos: https://www.ctralie.com/Research/SlidingWindowVideoQuasi. We analyze the resulting attractors using persistent homology, a technique which combines geometry and topology (Section 2.2), and we return scores in the range $[0,1]$ that indicate the degree of periodicity or quasiperiodicity in the corresponding video. We show that our periodicity measure compares favorable to others in the literature when ranking videos (Section 4.2). Furthermore, to our knowledge, there is no other method able to quantify the existence of quasiperiodicity directly from video data.

Our approach is fundamentally different from most others which quantify periodicity in video. For instance, it is common to derive 1D signals from the video and apply Fourier or autocorrelation to measure periodicity. By contrast, our technique operates on raw pixels, avoiding common video preprocessing and tracking entirely. Using geometry over Fourier/autocorrelation also has advantages for our applications. In fact, as a simple synthetic example shows (Figure 3), the Fourier Transform of quasiperiodic signals is often very close to the Fourier transform of periodic signals. By contrast, the sliding window embeddings we design yield starkly different geometric structures in the periodic and quasiperiodic cases. We exploit this to devise a quasiperiodicity measurement, which we use to indicate the degree of “biphonation” in videos of vibrating vocal folds (Section 4.3), which is useful in automatically diagnosing speech pathologies.

In the context of applied topology, our quasiperiodicity score is one of the first applications of persistent $H_{2}$ to high dimensional data, which is largely possible due to recent advancements in the computational feasibility of persistent homology [3].

1.1 Prior Work on Recurrence in Videos

1.1.1 1D Surrogate Signals

One common strategy for detecting periodicity in video is to derive a 1D function to act as a surrogate for its dynamics, and then to use either frequency domain (Fourier transform) or time domain (autocorrelation, peak finding) techniques. One of the earliest works in this genre finds level set surfaces in a spatiotemporal “XYT” volume of video (all frames stacked on top of each other), and then uses curvature scale space on curves that live on these “spatiotemporal surfaces” as the 1D function [1]. [34] use Fourier Transforms on pixels which exhibit motion, and define a measure of periodicity based on the energy around the Fourier peak and its harmonics. [10] extract contours and find eigenshapes from the contours to classify and parameterize motion within a period. Frequency estimation is done by using Fourier analysis and peak detection on top of other 1D statistics derived from the contours, such as area and center of mass. Finally, [47] derive a 1D surrogate function based on mutual information between the first and subsequent frames, and then look for peaks in the similarity function with the help of a watershed method.

1.1.2 Self-Similarity Matrices

Another class of techniques relies on self-similarity matrices (SSMs) between frames, where similarity can be defined in a variety of ways. [37] track a set of points on a foreground object and compare them with an affine invariant similarity. Another widely recognized technique for periodicity quantification [6], derives periodicity measures based on self-similarity matrices of L1 pixel differences. This technique has inspired a diverse array of applications, including analyzing the cycles of expanding/contracting jellyfish [33], analyzing bat wings [2], and analyzing videos of autistic spectrum children performing characteristic repetitive motions such as “hand flapping” [22]. We compare to this technique in Section 4.2.

1.1.3 Miscellaneous Techniques for Periodic Video Quantification

There are also a number of works that don’t fall into the two categories above. Some works focus solely on walking humans, since that is one of the most common types of periodic motion in videos of interest to people. [29] look at the “braiding patterns” that occur in XYT slices of videos of walking people. [17] perform blob tracking on the foreground of a walking person, and use the ratio of the second and first eigenvalues of PCA on that blob.

For more general periodic videos, [43] make a codebook of visual words and look for repetitions within the resulting string. [23] take a deep learning approach to counting the number of periods that occur in a video segment. They use a 3D convolutional neural network on spatially downsampled, non-sequential regions of interest, which are uniformly spaced in time, to estimate the length of the cycle. Finally, perhaps the most philosophically similar work to ours is the work of [41], who use cohomology to find maps of MOCAP data to the circle for parameterizing periodic motions, though this work does not provide a way to quantify periodicity.

1.1.4 Our Work

We show that geometry provides a natural way to quantify recurrence (i.e. periodicity and quasiperiodicity) in video, by measuring the shape of delay embeddings. In particular, we propose several optimizations (section 3) which make this approach feasible. The resulting measure of quasiperiodicity, for which quantitative approaches are lacking, is used in section 4.3 to detect anomalies in high-speed videos of vibrating vocal folds. Finally, in contrast to both frequency and time domain techniques, our method does not rely on the period length being an integer multiple of the sampling rate.

2 Background

2.1 Delay Embeddings And Their Geometry

Recurrence in video data can be captured via the geometry of delay embeddings; we describe this next.

2.1.1 Video Delay Embeddings

We will regard a video as a sequence of grayscale444For color videos we can treat each channel independently, yielding a vector in $\mathbb{R}^{W\times H\times 3}$. In practice, there isn’t much of a difference between color and grayscale embeddings in our framework for the videos we consider. image frames indexed by the positive real numbers. That is, given positive integers $W$ (width) and $H$ (height), a video with $W\times H$ pixels is a function

[TABLE]

In particular, a sequence of images $X_{1},X_{2},\ldots\in\mathbb{R}^{W\times H}$ sampled at discrete times $t_{1}<t_{2}<\cdots$ yields one such function via interpolation. For an integer $d\geq 0$, known as the dimension, a real number $\tau>0$, known as the delay, and a video $X:\mathbb{R}^{+}\longrightarrow\mathbb{R}^{W\times H}$, we define the sliding window (also referred to as time delay) embedding of $X$ – with parameters $d$ and $\tau$ – at time $t\in\mathbb{R}^{+}$ as the vector

[TABLE]

The subset of $\mathbb{R}^{W\times H\times(d+1)}$ resulting from varying $t$ will be referred to as the sliding window embedding of $X$. We remark that since the pixel measurement locations are fixed, the sliding window embedding is an “Eulerian” view into the dynamics of the video. Note that delay embeddings are generally applied to 1D time series, which can be viewed as 1-pixel videos ($W=H=1$) in our framework. Hence equation (1) is essentially the concatenation of the delay embeddings of each individual pixel in the video into one large vector. One of the main points we leverage in this paper is the fact that the geometry of the sliding window embedding carries fundamental information about the original video. We explore this next.

2.1.2 Geometry of 1-Pixel Video Delay Embeddings

As a motivating example, consider the harmonic (i.e. periodic) signal

[TABLE]

and the quasiperiodic signal

[TABLE]

We refer to $f_{h}$ as harmonic because its constitutive frequencies, $\frac{1}{10}$ and $\frac{1}{30}$, are commensurate; that is, they are linearly dependent over the rational numbers $\mathbb{Q}\subset\mathbb{R}$. By way of contrast, the underlying frequencies of the signal $f_{q}$, $\frac{1}{10}$ and $\frac{1}{10\pi}$, are linearly independent over $\mathbb{Q}$ and hence non-commensurate. We use the term quasiperiodicity, as in the non-linear dynamics literature [19], to denote the superposition of periodic processes whose frequencies are non-commensurate. This differs from other definitions in the literature (e.g. [43, robinson2016topological]) which regard quasiperiodic as any deviation from perfect repetition.

A geometric argument from [31] (see equation 7 below and the discussion that follows) shows that given a periodic function $f:[0,2\pi]\longrightarrow\mathbb{R}$ with exactly $N$ harmonics, if $d\geq 2N$ and $0<\tau<\frac{2\pi}{d}$ then the sliding window embedding $SW_{d,\tau}f$ is a topological circle (i.e. a closed curve without self-intersections) which wraps around an $N$-dimensional torus

[TABLE]

As an illustration, we show in Figure 1 a plot of $f_{h}$ and of its sliding window embedding $SW_{d,\tau}f_{h}$, via a PCA (Principal Component Analysis) 3-dimensional projection.

However, if $g:\mathbb{R}\longrightarrow\mathbb{R}$ is quasiperiodic with $N$ distinct non-commensurate frequencies then, for appropriate $d$ and $\tau$, $SW_{d,\tau}g$ is dense in (i.e. fills out) $\mathbb{T}^{N}$ [30]. Figure 2 shows a plot of the quasiperiodic signal $f_{q}(t)$ and a 3-dimensional projection, via PCA, of its sliding window embedding $SW_{d,\tau}f_{q}$.

The difference in geometry of the delay embeddings is stark compared to the difference between their power spectral densities, as shown in Figure 3.

Moreover, as we will see next, the interpretation of periodicity and quasiperiodicity as circularity and toroidality of sliding window embeddings remains true for videos with higher resolution (i.e. $\max\{W,H\}>1$). The rest of the paper will show how one can use persistent homology, a tool from the field of computational topology, to quantify the presence of (quasi)periodicity in a video by measuring the geometry of its associated sliding window embedding. In short, we propose a periodicity score for a video $X$ which measures the degree to which the sliding window embedding $SW_{d,\tau}X$ spans a topological circle, and a quasiperiodicity score which quantifies the degree to which $SW_{d,\tau}X$ covers a torus. This approach will be validated extensively: we show that our (quasi)periodicity detection method is robust under several noise models (motion blur, additive Gaussian white noise, and MPEG bit corruption); we compare several periodicity quantification algorithms and show that our approach is the most closely aligned with human subjects; finally, we provide an application to the automatic classification of dynamic regimes in high-speed laryngeal video-endoscopy.

2.1.3 Geometry of Video Delay Embeddings

Though it may seem daunting compared to the 1D case, the geometry of the delay embedding shares many similarities for periodic videos, as shown in [39]. Let us argue why sliding window embeddings from (quasi)periodic videos have the geometry we have described so far. To this end, consider an example video $X$ that contains a set of $N$ frequencies $\omega_{1},\omega_{2},...,\omega_{N}$. Let the amplitude of the $n^{\text{th}}$ frequency and $i^{\text{th}}$ pixel be $a_{in}$. For simplicity, but without loss of generality, assume that each is a cosine with zero phase offset. Then the time series at pixel $i$ can be written as

[TABLE]

Grouping all of the coefficients together into a $(W\times H)\times N$ matrix $A$, we can write

[TABLE]

where $A^{n}$ stands for the $n^{\text{th}}$ column of $A$. Constructing a delay embedding as in Equation 1:

[TABLE]

and applying the cosine sum identity, we get

[TABLE]

where $\vec{u}_{n},\vec{v}_{n}\in\mathbb{R}^{W\times H\times(d+1)}$ are constant vectors. In other words, the sliding window embedding of this video is the sum of linearly independent ellipses, which lie in the space of $d+1$ frame videos at resolution $W\times H$. As shown in [31] for the case of commensurate frequencies, when the window length is just under the length of the period, all of the $\vec{u}_{n}$ and $\vec{v}_{n}$ vectors become orthogonal, and so they can be recovered by doing PCA on $SW_{d,\tau}X(t)$. Figure 4 shows the components of the first 8 PCA vectors for a horizontal line of pixels in a video of an oscillating pendulum. Note how the oscillations are present both temporally and spatially.

2.1.4 The High Dimensional Geometry of Repeated Pulses

Using Eulerian coordinates has an important impact on the geometry of delay embeddings of natural videos. As Figure 5 shows, pixels often jump from foreground to background in a pattern similar to square waves.

These types of abrupt transitions require higher dimensional embeddings to reconstruct the geometry. To see why, first extract one period of a signal with period $\ell$ at a pixel $X_{i}(t)$:

[TABLE]

Then $X_{i}(t)$ can be rewritten in terms of the pulse as

[TABLE]

Since $X_{i}(t)$ repeats itself, regardless of what $f_{i}(t)$ looks like, periodic summation discretizes the frequency domain [32]

[TABLE]

Switching back to the time domain, we can write $X_{i}(t)$ as

[TABLE]

In other words, each pixel is the sum of some constant offset plus a (possibly infinite) set of harmonics at integer multiples of $\frac{1}{\ell}$. For instance, applying Equation 11 to a square wave of period $\ell$ centered at the origin is a roundabout way of deriving the Fourier Series

[TABLE]

by sampling the sinc function $\sin(\pi\ell f)/(\pi f)$ at intervals of $m/2\ell$ (every odd $m$ coincides with $\pi/2+k\pi$, proportional to $1/k$, and every even harmonic is zero conciding with $\pi k$). In general, the sharper the transitions are in $X_{i}(t)$, the longer the tail of $\mathcal{F}\{f_{i}(t)\}$ will be, and the more high frequency harmonics will exist in the embedding, calling for a higher delay dimension to fully capture the geometry, since every harmonic lives on a linearly independent ellipse. Similar observations about harmonics have been made in images for collections of patches around sharp edges ([48], Figure 2).

2.2 Persistent Homology

Informally, topology is the study of properties of spaces which do not change after stretching without gluing or tearing. For instance, the number of connected components and the number of (essentially different) 1-dimensional loops which do not bound a 2-dimensional disk, are both topological properties of a space. It follows that a circle and a square are topologically equivalent since one can deform one onto the other, but a circle and a line segment are not because that would require either gluing the endpoints of the line segment or tearing the circle. Homology [12] is a tool from algebraic topology designed to measure these types of properties, and persistent homology [50] is an adaptation of these ideas to discrete collections of points (e.g., sliding window embeddings). We briefly introduce these concepts next.

2.2.1 Simplicial Complexes

A simplicial complex is a combinatorial object used to represent and discretize a continuous space. With a discretization available, one can then compute topological properties by algorithmic means. Formally, a simplicial complex with vertices in a nonempty set $V$ is a collection $K$ of nonempty finite subsets $\sigma\subset V$ so that $\emptyset\neq\tau\subset\sigma\in K$ always implies $\tau\in K$. An element $\sigma\in K$ is called a simplex, and if $\sigma$ has $(n+1)$ elements then it is called an $n$-simplex. The cases $n=0,1,2$ are special, 0-simplices are called vertices, 1-simplices are called edges and 2-simplices are called faces. Here is an example to keep in mind: the circle $S^{1}=\{z\in\mathbb{C}:|z|=1\}$ is a continuous space but its topology can be captured by a simplicial complex $K$ with three vertices $a,b,c$, and three edges $\{a,b\}$, $\{b,c\}$, $\{a,c\}$. That is, in terms of topological properties, the simplicial complex

[TABLE]

can be regarded as a combinatorial surrogate for $S^{1}$: they both have 1 connected component, one 1-dimensional loop which does not bound a 2-dimensional region, and no other features in higher dimensions.

2.2.2 Persistent Homology of Point Clouds

The sliding window embedding of a video $X$ is, in practice, a finite set $\mathbb{SW}_{d,\tau}X=\{SW_{d,\tau}X(t):t\in T\}$ determined by a choice of $T\subset\mathbb{R}$ finite. Moreover, since $\mathbb{SW}_{d,\tau}X\subset\mathbb{R}^{W\times H\times(d+1)}$ then the restriction of the ambient Euclidean distance endows $\mathbb{SW}_{d,\tau}X$ with the structure of a finite metric space. Discrete metric spaces, also referred to as point clouds, are trivial from a topological point of view: a point cloud with $N$ points simply has $N$ connected components and no other features (e.g., holes) in higher dimensions. However, when a point cloud has been sampled from/around a continuous space with non-trivial topology (e.g., a circle or a torus), one would expect that appropriate simplicial complexes with vertices on the point cloud should reflect the topology of the underlying continuous space. This is what we will exploit next.

Given a point cloud $(\mathbb{X},d_{\mathbb{X}})$ – where $\mathbb{X}$ is a finite set and $d_{\mathbb{X}}:\mathbb{X}\times\mathbb{X}\longrightarrow[0,\infty)$ is a distance function – the Vietoris-Rips complex (or Rips complex for short) at scale $\epsilon\geq 0$ is the collection of non-empty subsets of $\mathbb{X}$ with diameter less than or equal to $\epsilon$:

[TABLE]

That is, $R_{\epsilon}(\mathbb{X})$ is the simplicial complex with vertex set equal to $\mathbb{X}$, constructed by adding an edge between any two vertices which are at most $\epsilon$ apart, adding all 2-dimensional triangular faces (i.e. 2-simplices) whose bounding edges are present, and more generally, adding all the $k$-simplices whose $(k-1)$-dimensional bounding facets have been included. We show in Figure 6 the evolution of the Rips complex on a set of points sampled around the unit circle.

The idea behind persistent homology is to track the evolution of topological features of complexes such as $R_{\epsilon}(\mathbb{X})$, as the scale parameter $\epsilon$ ranges from [math] to some maximum value $\epsilon_{\mathsf{max}}\leq\infty$. For instance, in Figure 6 one can see that $R_{0}(\mathbb{X})=\mathbb{X}$ has 40 distinct connect components (one for each point), $R_{0.30}(\mathbb{X})$ has three connected components and $R_{0.35}(\mathbb{X})$ has only one connected component; this will continue to be the case for every $\epsilon\geq 0.35$. Similarly, there are no closed loops in $R_{0}(\mathbb{X})$ or $R_{0.30}(\mathbb{X})$ bounding empty regions, but this changes when $\epsilon$ increases to $0.35$. Indeed, $R_{0.35}(\mathbb{X})$ has three 1-dimensional holes: the central prominent hole, and the two small ones to the left side. Notice, however, that as $\epsilon$ increases beyond $0.35$ these holes will be filled by the addition of new simplices; in particular, for $\epsilon>2$ one has that $R_{\epsilon}(\mathbb{X})$ will have only one connected component and no other topological features in higher dimensions.

The family $\mathcal{R}(\mathbb{X})=\{R_{\epsilon}(\mathbb{X})\}_{\epsilon\geq 0}$ is known as the Rips filtration of $\mathbb{X}$, and the emergence/dissapearence of topological features in each dimension (i.e., connected components, holes, voids, etc) as $\epsilon$ changes, can be codified in what are referred to as the persistence diagrams of $\mathcal{R}(\mathbb{X})$. Specifically, for each dimension $n=0,1,\ldots$ (0 = connected components, 1 = holes, 2 = voids, etc) one can record the value of $\epsilon$ for which a particular $n$-dimensional topological feature of the Rips filtration appears (i.e. its birth time), and when it disappears (i.e. its death time). The birth-death times $(b,d)\in\mathbb{R}^{2}$ of $n$-dimensional features for $\mathcal{R}(\mathbb{X})$ form a multiset $\mathsf{dgm}_{n}(\mathcal{R}(\mathbb{X}))$ — i.e. a set whose elements can come with repetition — known as the $n$-dimensional persistence diagram of the Rips filtration on $\mathbb{X}$. Since $\mathsf{dgm}_{n}(\mathcal{R}(\mathbb{X}))$ is just a collection of points in the region $\{(x,y)\in\mathbb{R}^{2}:0\leq x<y\}$, we will visualize it as a scatter plot. The persistence of a topological feature with birth-death times $(b,d)$ is the quantity $d-b$, i.e. its lifetime. We will also include the diagonal $y=x$ in the scatter plot in order to visually convey the persistence of each birth-death pair. In this setting, points far from the diagonal (i.e. with large persistence) represent topological features which are stable across scales and hence deemed significant, while points near the diagonal (i.e. with small persistence) are often associated with unstable features. We illustrate in Figure 7 the process of going from a point cloud to the 1-dimensional persistence diagram of its Rips filtration.

We remark that the computational task of determining all non-equivalent persistent homology classes of a filtered simplicial complex can, surprisingly, be reduced to computing the homology of a single simplicial complex [50, 3]. This is in fact a problem in linear algebra that can be solved via elementary row and column operations on appropriate boundary matrices.

The persistent homology of \mathcal{R}\big{(}\mathbb{SW}_{d,\tau}X\big{)}, and in particular its $n$-dimensional persistence diagrams for $n=1,2$, are the objects we will use to quantify periodicity and quasiperiodicity in a video $X$. Figures 8 and 9 show the persistence diagrams of the Rips filtrations, on the sliding window embeddings, for the commensurate and non-commensurate signals from Figures 1 and 2, respectively. We use fast new code from the “Ripser” software package to make persistent $H_{2}$ computation feasible [3].

3 Implementation Details

3.1 Reducing Memory Requirements with SVD

Suppose we have a video which has been discretely sampled at $N$ different frames at a resolution of $W\times H$, and we do a delay embedding with dimension $d$, for some arbitrary $\tau$. Assuming 32 bit floats per grayscale value, storing the sliding window embedding requires $4WHN(d+1)$ bytes. For a low resolution $200\times 200$ video only 10 seconds long at 30fps, using $d=30$ already exceeds 1GB of memory. In what follows we will address the memory requirements and ensuing computational burden to construct and access the sliding window embedding. Indeed, constructing the Rips filtration only requires pairwise distances between different delay vectors, this enables a few optimizations.

First of all, for $N$ points in $\mathbb{R}^{WH}$, where $N\ll WH$, there exists an $N$-dimensional linear subspace which contains them. In particular let $A$ be the $(W\times H)\times N$ matrix with each video frame along a column. Performing a singular value decomposition $A=USV^{T}$, yields a matrix $U$ whose columns form an orthonormal basis for the aforementioned $N$-dimensional linear subspace. Hence, by finding the coordinates of the original frame vectors with respect to this orthogonal basis

[TABLE]

and using the coordinates of the columns of $SV$ instead of the original pixels, we get a sliding window embedding of lower dimension

[TABLE]

for which

[TABLE]

Note that $SV$ can be computed by finding the eigenvectors of $A^{T}A$; this has a cost of $O(W^{2}H^{2}+N^{3})$ which is dominated by $W^{2}H^{2}$ if $WH\gg N$. In our example above, this alone reduces the memory requirements from 1GB to 10MB. Of course, this procedure is the most effective for short videos where there are actually many fewer frames than pixels, but this encompasses most of the examples in this work. In fact, the break-even point for a 200x200 30fps video is 22 minutes. A similar approach was used in the classical work on Eigenfaces [40] when computing the principal components over a set of face images.

3.2 Distance Computation via Diagonal Convolutions

A different optimization is possible if $\tau=1$; that is, if delays are taken exactly on frames and no interpolation is needed. In this case, the squared Euclidean distance between $SW_{d,1}X(i)$ and $SW_{d,1}X(j)$ is

[TABLE]

Let $D^{2}_{X}$ be the $N\times N$ matrix of all pairwise squared Euclidean distances between frames (possibly computed with the memory optimization in Section 3.1), and let $D^{2}_{Y}$ be the $(N-d)\times(N-d)$ matrix of all pairwise distances between delay frames. Then Equation 16 implies that $D^{2}_{Y}$ can be obtained from $D^{2}_{X}$ via convolution with a “rect function”, or a vector of $1$s of length $d+1$, over all diagonals in $D^{2}_{X}$ (i.e. a moving average). This can be implemented in time $O(N^{2})$ with cumulative sums. Hence, regardless of how $d$ is chosen, the computation and memory requirements for computing $D^{2}_{Y}$ depend only on the number of frames in the video. Also, $D_{Y}$ can simply be computed by taking the entry wise square root of $D^{2}_{Y}$, another $O(N^{2})$ computation. A similar scheme was used in [16] when comparing distances of 3D shape descriptors in videos of 3D meshes.

Figure 10 shows self-similarity matrices on embeddings of the pendulum video with no delay and with a delay approximately matching the period. The effect of a moving average along diagonals with delay eliminates the anti-diagonals caused by the video’s mirror symmetry.

Even for videos without mirror symmetries, such as a video of a running dog (Figure 11), introducing a delay brings the geometry into focus, as shown in Figure 12.

3.3 Normalization

A few normalization steps are needed in order to enable fair comparisons between videos with different resolutions, or which have a different range in periodic motion either spatially or in intensity. First, we perform a “point-center and sphere normalize” vector normalization which was shown in [31] to have nice theoretical properties.

That is,

[TABLE]

where $\boldsymbol{1}$ is a $WH(d+1)\times 1$ vector of all ones. In other words, one subtracts the mean of each component of each vector, and each vector is scaled so that it has unit norm (i.e. lives on the unit sphere in $\mathbb{R}^{WH(d+1)}$). Subtracting the mean from each component will eliminate additive linear drift on top of the periodic motion, while scaling addresses resolution / magnitude differences. Note that we can still use the memory optimization in Section 3.1, but we can no longer use the optimizations in Section 3.2 since each window is normalized independently.

Moreover, in order to mitigate nonlinear drift, we implement a simple pixel-wise convolution by the derivative of a Gaussian for each pixel in the original video before applying the delay embedding:

[TABLE]

This is a pixel-wise bandpass filter which could be replaced with any other bandpass filter leveraging application specific knowledge of expected frequency bounds. This has the added advantage of reducing the number of harmonics, enabling a smaller embedding dimesion $d$.

3.4 Periodicity/Quasiperiodicity Scoring

Once the videos are normalized to the same scale, we can score periodicity and quasiperiodicity based on the geometry of sliding window embeddings. Let $\mathsf{dgm}_{n}$ be the $n$-dimensional persistence diagram for the Rips filtration on the sliding window embedding of a video, and define $mp_{i}(\mathsf{dgm}_{n})$ as the $i$-th largest difference $d-b$ for $(b,d)\in\mathsf{dgm}_{n}$. In particular

[TABLE]

and $mp_{i}(\mathsf{dgm}_{n})\geq mp_{i+1}(\mathsf{dgm}_{n})$. We propose the following scores:

Periodicity Score (PS)

[TABLE]

Like [31], we exploit the fact that for the Rips filtration on $S^{1}$, the 1-dimensional persistence diagram has only one prominent birth-death pair with coordinates $\left(0,\sqrt{3}\right)$. Since this is the limit shape of a normalized perfectly periodic sliding window video, the periodicity score is between 0 (not periodic) and 1 (perfectly periodic). 2. 2.

Quasiperiodicity Score (QPS)

[TABLE]

This score is designed with the torus in mind. We score based on the second largest 1D persistence times the largest 2D persistence, since we want a shape that has two core circles and encloses a void to get a large score. Based on the K$\ddot{u}$nneth theorem of homology, the 2-cycle (void) should die the moment the smallest 1-cycle dies. 3. 3.

Modified Periodicity Score (MPS)

[TABLE]

We design a modified periodicity score which should be lower for quasiperiodic videos, than what the original periodicity score would yield.

Note that we use $\mathbb{Z}_{3}$ field coefficients for all persistent homology computations since, as shown by [31], this works better for periodic signals with strong harmonics. Before we embark on experiments, let us explore the choice of two crucial parameters for the sliding window embedding: the delay $\tau>0$ and the dimension $d\in\mathbb{N}$. In practice we determine an equivalent pair of parameters: the dimension $d$ and the window size $d\tau$.

3.5 Dimension and Window Size

Takens’ embedding theorem is one of the most fundamental results in the theory of dynamical systems [38]. In short, it contends that (under appropriate hypotheses) there exists an integer $D$, so that for all $d\geq D$ and generic $\tau>0$ the sliding window embedding $SW_{d,\tau}X$ reconstructs the state space of the underlying dynamics witnessed by the signal $X$. One common strategy for determining a minimal such $D$ is the false nearest-neighbors scheme [21]. The idea is to keep track of the $k$-th nearest neighbors of each point in the delay embedding, and if they change as $d$ is increased, then the prior estimates for $d$ were too low. This algorithm was used in recent work on video dynamics [42], for instance.

Even if we can estimate $d$, however, how does one choose the delay $\tau$? As shown in [31], the sliding window embedding of a periodic signals is roundest (i.e. so that the periodicity score $PS$ is maximized) when the window size, $d\tau$, satisfies the following relation:

[TABLE]

Here $L$ is number of periods that the signal has in $[0,2\pi]$ and $k\in\mathbb{N}$. To verify this experimentally, we show in Figure 13 how the periodicity score $PS$ changes as a function of window size for the pendulum video, and how the choice of window size from Equation 22 maximizes $PS$. To generate this figure we fixed a sufficiently large $d$ and varied $\tau$. Let us now describe the general approach: Given a video we perform a period-length estimation step (see section 3.6 next), which results in a positive real number $\ell$. For a given $d\in\mathbb{N}$ large enough we let $\tau>0$ be so that $d\tau=\ell\cdot\frac{d}{d+1}$.

3.6 Fundamental Frequency Estimation

Though Figure 13 suggests robustness to window size as long as the window is more than half of the period, we may not know what that is in practice. To automate window size choices, we do a coarse estimate using fundamental frequency estimation techniques on a 1D surrogate signal. To get a 1D signal, we extract the first coordinate of diffusion maps [4] using $10\%$ nearest neighbors on the raw video frames (no delay) after taking a smoothed time derivative. Note that a similar diffusion-based method was also used in recent work by [46] to analyze the frequency spectrum of a video of an oscillating 2 pendulum + spring system in a quasiperiodic state. Once we have the diffusion time series, we then apply the normalized autocorrelation method of [25] to estimate the fundamental frequency. In particular, given a discrete signal $x$ of length $N$, define the autocorrelation as

[TABLE]

However, as observed by [7], a more robust function for detecting periodicities is the squared difference function

[TABLE]

which can be rewritten as $d_{t}(\tau)=m_{t}(\tau)-2r_{t}(\tau)$ where

[TABLE]

Finally, [25] suggest normalizing this function to the range $[-1,1]$ to control for window size and to have an interpretation akin to a Pearson correlation coefficient:

[TABLE]

The fundamental frequency is then the inverse period of the largest peak in $n_{t}$ which is to the right of a zero crossing. The zero crossing condition helps prevent an offset of [math] from being the largest peak. Defining the normalized autocorrelation as in Equation 26 has the added advantage that the value of $n_{t}(\tau)$ at the peak can be used to score periodicity, which the authors call clarity. Values closer to 1 indicate more perfect periodicities. This technique will sometimes pick integer multiples of the period, so we multiply $n_{t}(\tau)$ by a slowly decaying envelope which is 1 for 0 lag and 0.9 for the maximum lag to emphasize smaller periods. Figure 14 shows the result of this algorithm on a periodic video, and Figure 15 shows the algorithm on an irregular video.

4 Experimental Evaluation

Next we evaluate the effectiveness of the proposed (Modified) Periodicity and Quasiperiodicity scores on three different tasks. First, we provide estimates of accuracy for the binary classifications periodic/not-periodic or quasiperiodic/not-quasiperiodic in the presence of several noise models and noise levels. The results illustrate the robustness of our method. Second, we quantify the quality of periodicity rankings from machine scores, as compared to those generated by human subjects. In a nutshell, and after comparing with several periodicity quantification algorithms, our approach is shown to be the most closely aligned with the perception of human subjects. Third, we demonstrate that our methodology can be used to automatically detect the physiological manifestations of certain speech pathologies (e.g., normal vs. biphonation), directly from high-speed videos of vibrating vocal folds.

4.1 Classification Under Varying Noise Levels/Models

As shown empirically in [8], a common source of noise in videos comes from camera shake (blur); this is captured by point spread functions resembling directed random walks [8, Figure 1] and the amount of blur (i.e. noise level) is controlled by the extent in pixels of the walk. Other sources are additive white Gaussian noise (awgn), controlled by the standard deviation of the Gaussian kernel, and MPEG bit errors quantified by the percentage of corrupted information. Figure 16 shows examples of these noise types.

For classification purposes we use three main recurrence classes. Three types of periodic videos (True periodic, TP): an oscillating pendulum, a bird flapping its wings, and an animation of a beating heart. Two types of quasiperiodic videos (True quasiperiodic, TQ): one showing two solid disks which oscillate sideways at non-commensurate rates, and the second showing two stationary Gaussian pulses with amplitudes non-commensurately modulated by cosine functions. Two videos without significant recurrence (True non-recurrent, TN): a video of a car driving past a landscape, and a video of an explosion. Each one of these seven videos is then corrupted by the three noise models at three different noise levels (blur = 20, 40, 80, awgn $\sigma=$ 1, 2, 3, bit error = 5, 10, 20%) as follows: given a particular video, a noise model and noise level, 600 instances are generated by sampling noise independently at random.

Results: We report in Table 1 the area under the Receiver Operating Characteristic (ROC) curve, or AUROC for short, for the classification task TP vs. TN (resp. TQ vs. TN) and binary classifier furnished by Periodicity (resp. Quasiperiodicity) Score.

For instance, for the Blur noise model with noise level of $80\times 80$ pixels, the AUROC from using the Periodicity Score to classify the 600 instances of the Heartbeat video as periodic, and the 600 instances of the Driving video as not periodic is 0.91. Similarly, for the MPEG bit corruption model with 5% of bit error, the AUROC from using the Quasiperiodicity score to classify the 600 instances of the Quasiperiodic Sideways Disks as quasiperiodic and the 600 instances of the Explosions video as not quasiperiodic is 0.92. To put these numbers in perspective, AUROC = 1 is associated with a perfect classifier and AUROC = 0.5 corresponds to classification by a random coin flip.

Overall, the type of noise that degrades performance the most across videos is the bit error, which makes sense, since this has the effect of randomly freezing, corrupting, or even deleting frames, which all interrupt periodicity. The blur noise also affects videos where the range of motion is small. The pendulum video, for instance, only moves over a range of 60 pixels at the most extreme end, so an 80x80 pixel blur almost completely obscures the motion.

4.2 Comparing Human and Machine Periodicity Rankings

Next we quantify the extent to which rankings obtained from our periodicity score (Equation 19), as well as three other methods, agree with how humans rank videos by periodicity. The starting point is a dataset of 20 different creative commons videos, each 5 seconds long at 30 frames per second. Some videos appear periodic, such as a person waving hands, a beating heart, and spinning carnival rides. Some of them appear nonperiodic, such as explosions, a traffic cam, and drone view of a boat sailing. And some of them are in between, such as the pendulum video with simulated camera shake.

It is known that humans are notoriously bad at generating globally consistent rankings of sets with more than 5 or 7 elements [27]. However, when it comes to binary comparisons of the type “should A be ranked higher than B?” few systems are as effective as human perception, specially for the identification of recurrent patterns in visual stimuli. We will leverage this to generate a globally consistent ranking of the 20 videos in our initial data set.

We use Amazon’s Mechanical Turk (AMT) [5] to present each pair of videos in the set of 20, $\binom{20}{2}=190$, each to three different users, for a total of 570 pairwise rankings. 15 unique AMT workers contributed to our experiment, using an interface as the one shown in Figure 17.

In order to aggregate this information into a global ranking which is as consistent as possible with the pairwise comparisons, we implement a technique known as Hodge rank aggregation [18]. Hodge rank aggregation finds the closest consistent ranking to a set of preferences, in a least squares sense. More precisely, given a set of objects $X$, and given a set of comparisons $P\subset X\times X$, we seek a scalar function $s$ on all of the objects that minimizes the following sum

[TABLE]

where $v_{ab}$ is a real number which is positive if $b$ is ranked higher than $a$ and negative otherwise. Thus, $s$ is a function whose discrete gradient best matches the set of preferences with respect to an $L^{2}$ norm. Note that the preferences that we feed the algorithm are based on the pairwise rankings returned from AMT. If video $b$ is greater than video $a$, then we assign $v_{ab}=1$, or -1 otherwise. Since we have 3 rankings for each video, we actually assign weights of +3, +1, -1, or -3. The +/- 3 are if all rankings agree in one direction, and the +/- 1 are if one of the rankings disagrees with the other two. Figure 18 shows a histogram of all of the weighted scores from users on AMT. They are mostly in agreement, though there are a few +/- 1 scores.

As comparison to the human scores, we use three different classes of techniques for machine ranking of periodicity.

Sliding Windows (SW): We sort the videos in decreasing order of Periodicity Score (Equation 19). We fix the window size at 20 frames and the embedding dimension at 20 frames (which is enough to capture 10 strong harmonics). We also apply a time derivative of width 10 to every frame.

Cutler-Davis [6]: The authors of this work present two different techniques to quantify periodicity from a self-similarity matrix (SSM) of video frames. The first is a frequency domain technique based on the peak of the average power spectral density over all columns (rows) of the SSM after linearly de-trending and applying a Hann window. To turn this into a continuous score, we report the ratio of the peak minus the mean over the standard deviation. This method will be referred to as Frequency Score.

As the authors warn, the frequency peak method has a high susceptibility to false positives. This motivated the design of a more robust technique in [6], which works by finding peaks in the 2D normalized autocorrelation of the Gaussian smoothed SSMs. For videos with mirror symmetry, the peaks will lie on a diamond lattice, while for videos without mirror symmetry, they will lie on a square lattice. After peak finding within neighborhoods, one simply searches over all possible lattices at all possible widths to find the best match with the peaks. Since each lattice is centered at the autocorrelation point $(0,0)$, no translational checks are necessary.

To turn this into a continuous score, let $E$ be the sum of Euclidean distances of the matched peaks in the autocorrelation image to the best fit lattice, let $r_{1}$ be the proportion of lattice points that have been matched, and let $r_{2}$ be the proportion of peaks which have been matched to a lattice point. Then we give the final periodicity score as

[TABLE]

A lattice which fits the peaks perfectly ($r_{1}=1$) with no error ($E=0$) and no false positive peaks ($r_{2}=1$) will have a score of 1, and any video which fails to have a perfectly matched lattice will have a score greater than 1. Hence, we sort in increasing order of the score to get a ranking.

As we will show, this technique agrees the second best with humans after our Periodicity Score ranking. One of the main drawbacks is numerical stability of finding maxes in non-isolated critical points around nearly diagonal regions in square lattices, which will erroneously inflate the score. Also, the lattice searching only occurs over an integer grid, but there may be periods that aren’t integer number of frames, so there will always be a nonzero $E$ for such videos. By contrast, our sliding window scheme can work for any real valued period length.

Diffusion Maps + Normalized Autocorrelation “Clarity”: Finally, we apply the technique from Section 3.6 to get an autocorrelation function, and we report the value of the maximum peak of the normalized autocorrelation to the right of a zero crossing, referred to as “clarity” by [25]. Values closer to 1 indicate more perfect repetitions, so we sort in descending order of clarity to get a ranking.

Figure 19 shows an example of these three different techniques on a periodic video. There is a dot which rises above the diagonal in the persistence diagram, a lattice is found which nearly matches the critical points in the autocorrelation image, and autocorrelation function on diffusion maps has a nice peak.

By contrast, for a nonperiodic video (Figure 20), there is hardly any persistent homology, there is no well matching lattice, and the first diffusion coordinate has no apparent periodicities.

Results: Once we have the global human rankings and the global machine rankings, we can compare them using the Kendall $\tau$ score [20]. Given a set of objects $N$ objects $X$ and two total orders $>_{1}$ and $>_{2}$, where $>(x_{a},x_{b})=1$ if $x_{a}>x_{b}$ and $>(x_{a},x_{b})=-1$ if $x_{a}<x_{b}$, the Kendall $\tau$ score is defined as

[TABLE]

For two rankings which agree exactly, the Kendall $\tau$ score will be 1. For two rankings which are exactly the reverse of each other, the Kendall $\tau$ score will be -1. In this way, it analogous to a Pearson correlation between rankings.

Table 2 shows the Kendall $\tau$ scores between all of the different machine rankings and the human rankings. Our sliding window video methodology (SW+TDA) agrees with the human ranking more than any other pair of ranking types. The second most similar are the SW and the diffusion clarity, which is noteworthy as they are both geometric techniques. Table 3 also shows the average run times, in milliseconds, of the different algorithms on each video on our machine. This does highlight one potential drawback of our technique, since TDA algorithms tend to be computationally intensive. However, at this scale (videos with at most several hundred frames), performance is reasonable.

4.3 Periodicity And Biphonation in High Speed Videos of Vocal Folds

In this final task we apply our methodology to a real world problem of interest in medicine. We show that our method can automatically detect certain types of voice pathologies from high-speed glottography, or high speed videos (4000 fps) of the left and right vocal folds in the human vocal tract [45, 44, 9]. In particular, we detect and differentiate quasiperiodicity from periodicity by using our geometric sliding window pipeline. Quasiperiodicity is a special case of what is referred to as “biphonation” in the biological context, where nonlinear phenomena cause a physical process to bifurcate into two different periodic modes, often during a transition to chaotic behavior [15]. The torus structure we sketched in Figure 9 has long been recognized in this context [28], but we provide a novel way of quantifying it.

Similar phenomena exist in audio [14, 15], but the main reason for studying laryngeal high speed video is understanding the biomechanical underpinnings of what is perceived in the voice. In particular, this understanding can potentially lead to practical corrective therapies and surgical interventions. On the other hand, the presence of biphonation in sound is not necessarily the result of a physiological phenomenon; it has been argued that it may come about as the result of changes in states of arousal [briefer2015segregation].

In contrast with our work, the existing literature on video-based techniques usually employs an inherently Lagrangian approach, where different points on the left and right vocal folds are tracked, and coordinates of these points are analyzed as 1D time series (e.g. [28, 35, 24, 13], [26]). This is a natural approach, since those are the pixels where all of the important signal resides, and well-understood 1D signal processing technique can be used. However, edge detectors often require tuning, and they can suddenly fail when the vocal folds close [24]. In our technique, we give up the ability to localize the anomalies (left/right, anterior/posterior) since we are not tracking them, but in return we do virtually no preprocessing, and our technique is domain independent.

Results: We use a collection of 7 high-speed videos for this analysis, drawn from a variety of different sources [49], [26], [28], [13]. There are two videos which correspond to “normal” periodic vocal folds, three which correspond to biphonation [28], and two which correspond to irregular motion555Please refer to supplementary material for an example video from each of these three classes. We manually extracted 400 frames per video (100milliseconds) and autotuned the window size based on autocorrelation of 1D diffusion maps (Section 3.6). We then chose an appropriate $\tau$ and chose a time spacing so that each point cloud would have 600 points. As shown in Table 4, our technique is able to differentiate between the four classes. We also show PCA and persistence diagrams for one example for each class. In Figure 21, we see what appears to be a loop in PCA, and one strong 1D persistent dot confirms this. In Figure 22, we see a prominent torus in the persistence diagram. In Figure 23, we don’t see any prominent structures in the persistence diagram, even though PCA looks like it could be a loop or a torus. Note, however, that PCA only preserves 13.7% of the variance in the signal, which is why high dimensional techniques are important to draw quantitative conclusions.

5 Discussion

We have shown in this work how applying sliding window embeddings to videos can be used to translate properties of the underlying dynamics into geometric features of the resulting point cloud representation. Moreover, we also showed how topological/geometric tools such as persistence homology can be leveraged to quantify the geometry of these embeddings. The pipeline was evaluated extensively showing robustness to several noise models, high quality in the produced periodicity rankings and applicability to the study of speech conditions form high-speed video data.

Moving forward, an interesting avenue related to medical applications is the difference between biphonation which occurs from quasiperiodic modes and biphonation which occurs from harmonic modes. [31] shows that $\mathbb{Z}_{3}$ field coefficients can be used to indicate the presence of a strong harmonic, so we believe a geometric approach is possible. This could be used, for example, to differentiate between subharmonic anomalies and quasiperiodic transitions [44].

Acknowledgments

The authors would like to thank Juergen Neubauer, Dimitar Deliyski, Robert Hillman, Alessandro de Alarcon, Dariush Mehta, and Stephanie Zacharias for providing videos of vocal folds. We also thank Matt Berger at ARFL for discussions about sliding window video efficiency, and we thank the 15 anonymous workers on the Amazon Mechanical Turk who ranked periodic videos.