Learning V1 Simple Cells with Vector Representation of Local Content and Matrix Representation of Local Motion

TL;DR

This paper introduces a model for representing local image content and motion in V1, using vector and matrix representations to capture local displacements, and demonstrates its ability to learn Gabor-like filters and infer local motions.

Contribution

The paper presents a novel vector-matrix representation model for local content and motion in V1, capable of learning simple cell-like filters and estimating optical flow.

Findings

Learned Gabor-like filters matching V1 simple cells

Model can infer local motions in supervised and unsupervised settings

Achieved competitive optical flow estimation results

Abstract

This paper proposes a representational model for image pairs such as consecutive video frames that are related by local pixel displacements, in the hope that the model may shed light on motion perception in primary visual cortex (V1). The model couples the following two components: (1) the vector representations of local contents of images and (2) the matrix representations of local pixel displacements caused by the relative motions between the agent and the objects in the 3D scene. When the image frame undergoes changes due to local pixel displacements, the vectors are multiplied by the matrices that represent the local displacements. Thus the vector representation is equivariant as it varies according to the local displacements. Our experiments show that our model can learn Gabor-like filter pairs of quadrature phases. The profiles of the learned filters match those of simple cells in Macaque V1. Moreover, we demonstrate that the model can learn to infer local motions in either a supervised or unsupervised manner. With such a simple model, we achieve competitive results on optical flow estimation.

Figures

Tables

Table 1. Table 1: Average endpoint error of the inferred displacement and number of parameters. Abbreviation FN2 refers to FlowNet 2.0. ‘Ours-ref’ indicates that the results are post-processed by the refinement CNN.

	V1Deform		V1FlyingObjects		MPI-Sintel		# params (M)
	pre-trained	trained	pre-trained	trained	pre-trained	trained
FN2-C	1.324	1.130	0.852	1.034	0.363	0.524	39.18
FN2-S	1.316	0.213	0.865	0.261	0.410	0.422	38.68
FN2-CS	0.713	0.264	0.362	0.243	0.266	0.346	77.87
FN2-CSS	0.629	0.301	0.299	0.303	0.234	0.450	116.57
FN2	0.686	0.205	0.285	0.265	0.146	0.278	162.52
Ours	-	0.258	-	0.442	-	0.337	1.82
Ours-ref	-	0.156	-	0.202	-	0.140	1.84

Table 2. Table 2: Average endpoint error of the inferred displacement in unsupervised learning.

	Unsup	UnFlow-C	UnFlow-CS	UnFlow-CSS	Ours
V1FlyingObjects (train)	0.379	0.336	0.374	0.347	0.245
V1FlyingObjects (test)	0.811	0.399	0.394	0.453	0.316
MPI-Sintel	0.440	0.198	0.248	0.202	0.101

Table 3. Table 3: Ablation study measured by average endpoint error (AEE) of motion inference and mean squared error (MSE) of multi-step animation, learned from V1Deform dataset.

	Sub-vector dimension
	2	4	6	8	12
AEE	0.258	0.246	0.238	0.241	0.0.246
MSE	8.122	7.986	7.125	7. 586	8.017

Equations

$$v(x)=W\mathbf{I}[x],x\in D_{-},$$

$$\mathbf{I}={\mathbf{W}}^{\top }V=\sum _{x\in D_{-}}{W}^{\top }(x)v(x),$$

$$\displaystyle\begin{array}[c]{ccc}{v}_{t}(x)&\stackrel{{\scriptstyle M(\delta(x))\times}}{{\xrightarrow{\hskip 28.45274pt}}}&{v}_{t+1}(x)\\ &&\\ W\uparrow&\uparrow&\uparrow W\\ {\bf I}_{t}[x]&\stackrel{{\scriptstyle\delta(x)}}{{\xrightarrow{\hskip 28.45274pt}}}&{\bf I}_{t+1}[x]\end{array}$$

$$v_{t+1}(x)=M(\delta (x))v_t(x),\forall x\in D_{-}.$$

$$M(\delta )=I+B_1{\delta }_1+B_2{\delta }_2+B_{11}{\delta }_1^2+B_{22}{\delta }_2^2+B_{12}{\delta }_1{\delta }_2.$$

$$v_{t+1}(x)=\sum _{\text{d}x\in \mathcal{S}}M(\delta (x),\text{d}x)v_t(x+\text{d}x)$$

$$v^{(k)}(x)=W^{(k)}\mathbf{I}[x],k=1,...,K,$$

$$M(\delta )=\text{diag}(M^{(k)}(\delta ),k=1,...,K)$$

$$v_{t+1}^{(k)}(x)=M^{(k)}(\delta )v_t^{(k)}(x),k=1,...,K.$$

$$L_1=\sum _{k=1}^K\sum _{x\in D_{-}}{\Vert W^{(k)}{\mathbf{I}}_{t+1}[x]-M^{(k)}(\delta (x))W^{(k)}{\mathbf{I}}_t[x]\Vert }^2.$$

$$L_2=\sum _{s\in t,t+1}{\Vert {\mathbf{I}}_s-{\mathbf{W}}^{\top }\mathbf{W}{\mathbf{I}}_s\Vert }^2.$$

$$L_{1,x}(\delta )=\sum _{k=1}^K{\Vert W^{(k)}{\mathbf{I}}_{t+1}[x]-M^{(k)}(\delta )W^{(k)}{\mathbf{I}}_t[x]\Vert }^2.$$

$${\sum }_x{L}_{1,x}(\delta )+\Vert ▽\delta {\Vert }^2+{\Vert {\mathbf{I}}_{t+1}-\mathrm{warp}({\mathbf{I}}_t,\delta )\Vert }^2.$$

$$u=\sum _{i=0}^{m}M(\delta {)}^{m-i}{v}_{t+i,x}=\sum _{i=0}^{m}M(\delta {)}^{m-i}W{\mathbf{I}}_{t+i}[x]$$

$${\hat{v}}_{t+1}(x,\delta )$$

$${\delta }_{t+1}(x)$$

Full text

1 Introduction

Our understanding of the primary visual cortex or V1 (Hubel and Wiesel 1959) is still very limited (Olshausen and Field 2005). In particular, mathematical and representational models for V1 are still in short supply. Two prominent examples of such models are sparse coding (Olshausen and Field 1997) and independent component analysis (ICA) (Bell and Sejnowski 1997). Although such models may not provide detailed explanations at the level of neuronal dynamics, they help us understand the computational problems being solved by V1.

In this paper, we propose a model of this sort. It is a representational model of natural image pairs that are related by local pixel displacements. The image pairs can be consecutive frames of a video sequence, where the local pixel displacements are caused by the relative motions between the agent and the objects in the 3D environment. Perceiving such local motions can be crucial for inferring ego-motion, object motions, and 3D depth information.

As is the case with existing models, we expect our model to explain only partial aspects of V1, including: (1) The receptive fields of V1 simple cells resemble Gabor filters (Daugman 1985). (2) Adjacent simple cells have quadrature phase relationship (Pollen and Ronner 1981; Emerson and Huang 1997). (3) The V1 cells are capable of perceiving local motions. While existing models can all explain (1), our model can also account for (2) and (3) naturally. Compared to models such as sparse coding and ICA, our model has a component that serves a direct purpose of perceiving local motions.

Our model consists of the following two components.

(1) Vector representation of local image content. The local content around each pixel is represented by a high dimensional vector. Each unit in the vector is obtained by a linear filter. These local filters or wavelets are assumed to form a normalized wavelet tight frame, i.e., the image can be reconstructed from the vectors using the linear filters as the basis functions.

(2) Matrix representation of local displacement. The change of the image from the current time frame to the next time frame is caused by the displacements of the pixels. Each possible displacement is represented by a matrix that operates on the vector. When the image changes according to the displacements, the vector at each pixel is multiplied by the matrix that represents the local displacement, in other words, the vector at each pixel is rotated by the matrix representation of the displacement of this pixel. Thus the vector representation is equivariant as it varies according to the local displacements. See Fig. 1 for an illustration.

We train the model on image pairs where in each pair, the second image is a deformed version of the first image, and the deformation is either known or inferred in training. We learn the encoding matrices for vector representation and the matrix representation of local displacements from the training data.

Experiments show that our model learns V1-like units that can be well approximated by Gabor filters with quadrature phase relationship. The profiles of learned units match those of simples cells in Macaque V1. After learning the encoding matrices for vector representation and the matrix representations of the displacements, we can infer the displacement field using the learned model. Compared to popular optical flow estimation methods (Dosovitskiy et al. 2015; Ilg et al. 2017), which use complex deep neural networks to predict the optical flows, our model is much simpler and is based on explicit vector and matrix representations. We demonstrate comparable performance to these methods in terms of the inference of displacement fields.

In terms of biological interpretation, the vectors can be interpreted as activities of groups of neurons, and the matrices can be interpreted as synaptic connections. See Sections 5.1 and 5.2 for details.

2 Contributions and related work

This paper proposes a simple representational model that couples the vector representations of local image contents and matrix representations of local pixel displacements. The model explains certain aspects of V1 simple cells such as Gabor-like receptive fields and quadrature phase relationship, and adds to our understanding of V1 motion perception in terms of representation learning.

The following are two themes of related work.

(1) V1 models. Most well known models for V1 are concerned with statistical properties of natural images or video sequences. Examples include sparse coding model (Olshausen and Field 1997; Lewicki and Olshausen 1999; Olshausen 2003), independent component analysis (ICA) (Hyvärinen, Karhunen, and Oja 2004; Bell and Sejnowski 1997; van Hateren and Ruderman 1998), slowness criterion (Hyvärinen, Hurri, and Väyrynen 2003; Wiskott and Sejnowski 2002), and prediction (Singer et al. 2018). While these models are very compelling, they do not account for perceptual inference explicitly, which is the most important functionality of V1. On the other hand, our model is learned for the direct purpose of perceiving local motions caused by relative motion between the agent and the surrounding 3D environment. In fact, our model is complementary to the existing models for V1: similar to existing models, our work assumes a linear generative model for image frames, but our model adds a relational component with matrix representation that relates the consecutive image frames. Our model is also complementary to slowness criterion in that when the vectors are rotated by matrices, the norms of the vectors remain constant.

(2) Matrix representation. In representation learning, it is a common practice to encode the signals or states as vectors. However, it is a much less explored theme to represent the motions, actions or relations by matrices that act on the vectors. An early work in this theme is (Paccanaro and Hinton 2001), which learns matrices to represent relations. More recently, (Jayaraman and Grauman 2015) learns equivariant representation with matrix representation of ego-motion. (Zhu et al. 2021) learns generative models of posed images based on invariant representation of 3D scene with matrix representation of ego-motion (Gao et al. 2019, 2021) learn vector representation of self-position and matrix representation of self-motion in a representational model of grid cells. Our work constitutes a new development along this theme.

3 Representational model

3.1 Vector representation

Let $\{{\bf I}(x),x\in D\}$ be an image observed at a certain instant, where $x=(x_{1},x_{2})\in D$ is the 2D coordinates of pixel. $D$ is the image domain (e.g., $128\times 128$). We represent the image ${\bf I}$ by vectors $\{v(x),x\in D_{-}\}$, where each $v(x)$ is a vector defined at pixel $x$, and $D_{-}$ may consist of a sub-sampled set of pixels in $D$. ${V}=\{v(x),x\in D_{-}\}$ forms a vector representation of the whole image.

Conventionally, we assume the vector encoding is linear and convolutional. Specifically, let ${\bf I}[x]$ be a squared patch (e.g., $16\times 16$) of ${\bf I}$ centered at $x$. We can flatten ${\bf I}[x]$ into a vector (e.g., 256 dimensional) and let

[TABLE]

be the linear encoder, where $W$ is the encoding matrix that encodes ${\bf I}[x]$ into a vector $v(x)$. The rows of $W$ are the linear filters and can be displayed as local image patches of the same size as the image patch ${\bf I}[x]$. We can further write ${V}={\bf W}{\bf I}$ if we treat ${\bf I}$ as a vector, and the rows of ${\bf W}$ are the translated versions of $W$.

3.2 Tight frame auto-encoder

We assume that ${\bf W}$ is an auto-encoding tight frame. Specifically, let $W(x)$ denote the translation of filter $W$ to pixel $x$ and zero-padding the pixels outside the filters, so that each row of $W(x)$ is of the same dimension as ${\bf I}$. We then assume

[TABLE]

i.e., the linear filters for bottom-up encoding also serve as basis functions for top-down decoding. Both the encoder and decoder can be implemented by convolutional linear neural networks.

The tight frame assumption can be justified by the fact that under that assumption, for two images ${\bf I}$ and ${\bf J}$, we have $\langle{\bf W}{\bf I},{\bf W}{\bf J}\rangle={\bf I}^{\top}{\bf W}^{\top}{\bf W}{\bf J}=\langle{\bf I},{\bf J}\rangle,$ i.e., the vector representations preserve the inner product. As a result, we also have $\|{\bf W}{\bf I}\|=\|{\bf I}\|$ and $\|{\bf W}{\bf J}\|=\|{\bf J}\|$, so that the vector representations ${\bf W}{\bf I}$ and ${\bf W}{\bf J}$ preserve the angle between the images ${\bf I}$ and ${\bf J}$ and has the isometry property. That is, when the image ${\bf I}$ changes from ${\bf I}_{t}$ to ${\bf I}_{t+1}$, its vector representation $V$ changes from $V_{t}$ to $V_{t+1}$, and the angle between ${\bf I}_{t}$ and ${\bf I}_{t+1}$ is the same as the angle between $V_{t}$ and $V_{t+1}$.

In this paper, we assume a tight frame auto-encoder for computational convenience. A more principled treatment is to treat the decoder as a top-down generative model, and treat the encoder as approximate inference. Sparsity constraint can be imposed on the top-down decoder model.

3.3 Matrix representation

Let ${\bf I}_{t}$ be the image at time frame $t$. Suppose the pixels of ${\bf I}_{t}$ undergo local displacements $(\delta(x),\forall x)$, where $\delta(x)$ is the displacement at pixel $x$. The image transforms from ${\bf I}_{t}$ to ${\bf I}_{t+1}$. We assume that $\delta(x)$ is within a squared range $\Delta$ (e.g., $[-6,6]\times[-6,6]$ pixels) that is inside the range of the local image patch ${\bf I}_{t}[x]$ (e.g., $16\times 16$ pixels). We assume that the displacement field $(\delta(x),\forall x)$ is locally smooth, i.e., pixels within each local image patch undergoes similar displacements. Let $v_{t}(x)$ and $v_{t+1}(x)$ be the vector representations of ${\bf I}_{t}[x]$ and ${\bf I}_{t+1}[x]$ respectively. The transformation from ${\bf I}_{t}[x]$ to ${\bf I}_{t+1}[x]$ is illustrated by the following diagram:

[TABLE]

Specifically, we assume that

[TABLE]

That is, when ${\bf I}$ changes from ${\bf I}_{t}$ to ${\bf I}_{t+1}$, ${v}(x)$ undergoes a linear transformation, driven by a matrix $M(\delta(x))$, which depends on the local displacement $\delta(x)$. Thus $v(x)$ is an equivariant representation as it varies according to $\delta(x)$.

One motivation for modeling the transformation as a matrix representation operating on the vector representation of image patch comes from Fourier analysis. Specifically, an image patch ${\bf I}[x]$ can be expressed by the Fourier decomposition ${\bf I}[x]=\sum_{k}c_{k}e^{i\langle\omega_{k},x\rangle}$. Assuming all the pixels in the patch are shifted by a constant displacement $\delta(x)$, then the shifted image patch becomes ${\bf I}(x-\delta(x))=\sum_{k}c_{k}e^{-i\langle\omega_{k},\delta(x)\rangle}e^{i\langle\omega_{k},x\rangle}$. The change from the complex number $c_{k}$ to $c_{k}e^{-i\langle\omega_{k},dx\rangle}$ corresponds to rotating a 2D vector by a $2\times 2$ matrix. However, we emphasize that our model does not assume Fourier basis or its localized version such as Gabor filters. The model figures it out with generic vector and matrix representations.

Parametrization.

We consider two ways to parametrize the matrix representation of local displacement $M(\delta(x))$. First, we can discretize the displacement $\delta(x)$ into a finite set of possible values $\{\delta\}$, and we learn a separate $M(\delta)$ for each $\delta$. Second, We can learn a parametric version of $M(\delta)$ as the second order Taylor expansion of a matrix-valued function of $\delta=(\delta_{1},\delta_{2})$,

[TABLE]

In the above expansion, $I$ is the identity matrix, and $B=(B_{1},B_{2},B_{11},B_{22},B_{12})$ are matrices of coefficients of the same dimensionality as $M(\delta)$.

A more careful treatment is to treat $M(\delta)$ as forming a matrix Lie group, and write $M(\delta)$ as an exponential map of generator matrices that span the matrix Lie algebra. By assuming skew-symmetric generator matrices, $M(\delta)$ will be automatically a rotation matrix. The exponential map can be approximated by Taylor expansion similar to 9.

Local mixing.

If $\delta(x)$ is large, $v_{t+1}(x)$ may contain information from adjacent image patches of ${\bf I}_{t}$ in addition to ${\bf I}_{t}[x]$. To address this problem, we can generalize the motion model (Eq. (8)) to allow local mixing of encoded vectors. Specifically, let $\mathcal{S}$ be a local support centered at [math]. We assume that

[TABLE]

In the learning algorithm, we discretize $\text{d}x$ and learn a separate $M(\delta,\text{d}x)$ for each $\text{d}x$.

3.4 Sub-vectors and disentangled rotations

The vector $v(x)$ can be high-dimensional. For computational efficiency, we further divide $v(x)$ into $K$ sub-vectors, $v(x)=(v^{(k)}(x),k=1,...,K)$. Each sub-vector is obtained by an encoding sub-matrix $W^{(k)}$, i.e.,

[TABLE]

where $W^{(k)}$ consists of the rows of $W$ that correspond to $v^{(k)}$. In practice, we find that this assumption is necessary for the emergence of V1-like receptive field.

Correspondingly, the matrix representation

[TABLE]

is block diagonal. Each sub-vector $v^{(k)}(x)$ is transformed by its own sub-matrix $M^{(k)}(\delta)$:

[TABLE]

The linear transformations of the sub-vectors $v^{(k)}(x)$ can be considered as rotations. $v(x)$ is like a multi-arm clock, with each arm $v^{(k)}(x)$ rotated by $M^{(k)}(\delta(x))$. The rotations of $v^{(k)}(x)$ for different $k$ and $x$ are disentangled, meaning that the rotation of a sub-vector does not depend on other sub-vectors.

The assumption of sub-vectors and block-diagonal matrices is necessary for learning Gabor-like filters. For ${v}_{t+1}(x)=M(\delta(x)){v}_{t}(x)$, it is equivalent to $\tilde{v}_{t+1}(x)=\tilde{M}(\delta(x))\tilde{v}_{t}(x)$ if we let $\tilde{v}_{t}(x)=Pv_{t}(x)$, $\tilde{v}_{t+1}(x)=Pv_{t+1}(x)$, and $\tilde{M}(\delta(x))=PM(\delta(x))P^{-1}$, for an invertible $P$. Assuming block-diagonal matrices helps eliminates such ambiguity. More formally, a matrix representation is irreducible if it cannot be further diagonalized into smaller block matrices. Assuming block-diagonal matrices helps to make the matrix representation close to irreducible.

4 Learning and inference

4.1 Supervised learning

The input data consist of the triplets $({\bf I}_{t},(\delta(x),x\in D_{-}),{\bf I}_{t+1})$, where $(\delta(x),x\in D_{-})$ is the given displacement field. The unknown parameters to learn consist of matrices $(W,M(\delta),\delta\in\Delta)$, where $\Delta$ is the range of $\delta$. In the case of parametric $M$, we learn the $B$ matrices in the second order Taylor expansion (Eq. (9)). We assume that there are $K$ sub-vectors so that $M$ or $B$ are block-diagonal matrices. We learn the model by optimizing a loss function defined as a weighted sum of two loss terms, based on the linear transformation model (Eq. (8)) and tight frame assumption (Eq. (2)) respectively:

(1) Linear transformation loss

[TABLE]

For local mixing generalization, we substitute $M^{(k)}(\delta(x))$ by $\sum_{\text{d}x\in\mathcal{S}}M^{(k)}(\delta(x),\text{d}x)$.

(2) Tight frame auto-encoder loss

[TABLE]

4.2 Inference of motion

After learning $(W,M(\delta),\delta\in\Delta)$, given a testing pair $({\bf I}_{t},{\bf I}_{t+1})$, we can infer the pixel displacement field $(\delta(x),x\in D_{-})$ by minimizing the linear transformation loss: $\delta(x)=\arg\max_{\delta\in\Delta}L_{1,x}(\delta)$, where

[TABLE]

This algorithm is efficient in nature as it can be parallelized for all $x\in D_{-}$ and for all $\delta\in\Delta$.

If we use a parametric version of $(M(\delta),\delta\in\Delta)$ (Eq. (9)), we can minimize $\sum_{x}L_{1,x}(\delta)$ using gradient descent with $\delta$ initialized from random small values. To encourage the smoothness of the inferred displacement field, we add a penalty term $\|\triangledown\delta(x)\|^{2}$ for this setting.

4.3 Unsupervised learning

We can easily adapt the learning of the model to an unsupervised manner, without knowing the pixel displacement field $(\delta(x),x\in D_{-})$. Specifically, we can iterate the following two steps: (1) update model parameters by loss functions defined in Section 4.1; (2) infer the displacement field as described in 4.2. To eliminate the ambiguity of $M(\delta)$ with respect to $\delta$, we add a regularization term $\left\|{\bf I}_{t+1}-{\rm warp}({\bf I}_{t},\delta)\right\|^{2}$ in the inference step, where ${\rm warp}(\cdot,\cdot)$ is a differentiable warping function, and we use the parametric version of $(M(\delta),\delta\in\Delta)$. To summarize, we infer the displacement field $(\delta(x),x\in D_{-})$ by minimizing:

[TABLE]

In practice, for each image pair at each iteration, we start the inference by running gradient descent on the inferred displacement field from the previous iteration.

5 Discussions about model

5.1 Biological interpretations of cells and synaptic connections

The learned $(W,M(\delta)),\delta)$ can be interpreted as synaptic connections. Specifically, for each block $k$, $W^{(k)}$ corresponds to one set of connection weights. Suppose $\delta\in\Delta$ is discretized, then for each $\delta$, $M^{(k)}(\delta)$ corresponds to one set of connection weights, and $(M^{(k)}(\delta),\delta\in\Delta)$ corresponds to multiple sets of connection weights. For motion inference in a biological system, after computing $v_{t,x}^{(k)}=W^{(k)}{\bf I}_{t}[x]$, $M^{(k)}(\delta)v_{t,x}^{(k)}$ can be computed simultaneously for every $\delta\in\Delta$. Then $\delta(x)$ is inferred by max pooling according to Eq. (16).

$v_{t,x}^{(k)}$ can be interpreted as activities of simple cells, and $\|v_{t,x}^{(k)}\|^{2}$ can be interpreted as activity of a complex cell. If $M^{(k)}(\delta)$ is close to a rotation matrix, then we have norm stability so that $\|v_{t,x}^{(k)}\|\approx\|v_{t+1,x}^{(k)}\|$, which is closely related to the slowness property (Hyvärinen, Hurri, and Väyrynen 2003; Wiskott and Sejnowski 2002).

5.2 Spatiotemporal filters and recurrent implementation

If we enforce norm stability or the orthogonality of $M^{(k)}(\delta)$, then minimizing $\|v_{t+1,x}-M(\delta)v_{t,x}\|^{2}$ over $\delta\in\Delta$ is equivalent to maximizing $\langle v_{t+1,x},M(\delta)v_{t,x}\rangle$, which in turn is equivalent to maximizing $\|v_{t+1,x}+M(\delta)v_{t,x}\|^{2}$ so that $v_{t+1,x}$ and $M(\delta)v_{t,x}$ are aligned. This alignment criterion can be conveniently generalized to multiple consecutive frames, so that we can estimate the velocity at $x$ by maximizing the $m$-step alignment score $\|u\|^{2}$, where

[TABLE]

consists of responses of spatiotemporal filters or “animated” filters $(M(\delta)^{m-i}W,i=0,...,m)$, and $\|u\|^{2}$ corresponds to the energy of motion $\delta$ in the motion energy model (Adelson and Bergen 1985) for direction selective cells. Thus our model is connected with the motion energy model. Moreover, our model enables a recurrent network for computing $u$ by $u_{i}=v_{t+i,x}+M(\delta)u_{i-1}$ for $i=0,...,m$, with $u_{-1}=0$, and $u=u_{m}$. This recurrent implementation is much more efficient and biologically plausible than the plain implementation of spatiotemporal filtering which requires memorizing all the ${\bf I}_{t+i}$ for $i=0,...,m$. See (Pachitariu and Sahani 2017) for a discussion of biological plausibility of recurrent implementation of spatiotemporal filtering in general.

The spatiotemporal filters can also serve as spatiotemporal basis functions for the top-down decoder model.

6 Experiments

The code, data and more results can be found at http://www.stat.ucla.edu/~ruiqigao/v1/main.html

We learn our model $(W,M(\delta),\delta\in\Delta)$ from image pairs $({\bf I}_{t},{\bf I}_{t+1})$ with its displacement field $(\delta(x))$ known or unknown. The number of sub-vectors $K=40$, and the number of units in each sub-vector $v^{(k)}(x)$ is 2. We use Adam (Kingma and Ba 2014) optimizer for updating the model. To demonstrate the efficacy of the proposed model, we conduct experiments on two new synthetic datasets (V1Deform and V1FlyingObjects) and two public datasets (MPI-Sintel and MUG Facial Expression). The motivation to generate the synthetic datasets is that we find existing datasets such as Flying Chairs (Dosovitskiy et al. 2015), FlyingThings3D (Mayer et al. 2016), and KITTI flow (Geiger, Lenz, and Urtasun 2012) contain image pairs with fairly large motions, which are unlikely consecutive frames perceived by V1. Thus we generate two synthetic datasets: V1Deform and V1FlyingObjects, which contains image pairs with only small local displacements and therefore better serve our purpose of studying motion perception in V1. See Fig. 4 for some examples from the synthetic datasets.

6.1 Datasets

In this subsection, we elaborate the generation process of the two new synthetic datasets, and introduce the public datasets we use in this work.

V1Deform. For this dataset, we consider random smooth deformations for natural images. Specifically, We obtain the training data by collecting static images for $({\bf I}_{t})$ and simulate the displacement field $(\delta(x))$. The simulated displacement field is then used to transform ${\bf I}_{t}$ to obtain ${\bf I}_{t+1}$. We retrieve natural images as ${\bf I}_{t}$ from MIT places365 dataset (Zhou et al. 2016). The images are scaled to 128 $\times$ 128. We sub-sample the pixels of images into a $m\times m$ grid ($m=4$ in the experiments), and randomly generate displacements on the grid points, which serve as the control points for deformation. Then $\delta(x)$ for $x\in D$ can be obtained by spline interpolation of the displacements on the control points. We get ${\bf I}_{t+1}$ by warping ${\bf I}_{t}$ using $\delta(x)$ (Jaderberg et al. 2015). When generating a displacement $\delta=(\delta_{1},\delta_{2})$, both $\delta_{1}$ and $\delta_{2}$ are randomly sampled from a range of $[-6,+6]$. We synthesize $20,000$ pairs for training and $3,000$ pairs for testing.

V1FlyingObjects. For this dataset, we consider separating the displacement field into motions of the background and foreground, to jointly simulate the self-motion of the agent and the motion of the objects in the natural 3D scenes. To this end, we apply affine transformations to background images collected from MIT places365 (Zhou et al. 2016) and foreground objects from a public 2D object dataset COIL-100 (Nene et al. 1996). The background images are scaled to $128\times 128$, and the foreground images are randomly rescaled. To generate motion, we randomly sample affine parameters of translation, rotation, and scaling for both the foreground and background images. The motions of the foreground objects are relative to the background images, which can be explained as the relative motion between the moving object and agent. We tune the distribution of the affine parameters to keep the range of the displacement fields within $[-6,+6]$, which is consistent with the V1Deform dataset. Together with the mask of the foreground object and the sampled transformation parameters, we render the image pair $({\bf I}_{t},{\bf I}_{t+1})$ and its displacement field $(\delta(x))$ for each pair of the background image and foreground image.

For the foreground objects, we obtain t he estimated masks from (tev 2006), resulting in 96 objects with $72$ views per object available. We generate $14,411$ synthetic image pairs with their corresponding displacement fields and further split $12,411$ pairs for training and $2,000$ pairs for testing. Compared with previous optical flow dataset like Flying Chairs (Dosovitskiy et al. 2015) and scene flow dataset like FlyingThings3D (Mayer et al. 2016), the proposed V1FlyingObjects dataset has various foreground objects with more realistic texture and smoother displacement fields, which simulates more realistic environments.

We shall release the two synthetic datasets, which are suitable for studying local motions and perceptions. Besides, we also use two public datasets:

MPI-Sintel. MPI-Sintel (Butler et al. 2012; Wulff et al. 2012) is a public dataset designed for the evaluation of optical flow derived from rendered artificial scenes, with special attention to realistic image properties. Since MPI-Sintel is relatively small, which contains around a thousand image pairs, we use it only for testing the learned models in the inference of the displacement field. We use the final version of MPI-Sintel and resize each frame into size $128\times 128$. We select frame pairs whose motions are within the range of $[-6,+6]$, resulting in $384$ frame pairs in total.

MUG Facial Expression. MUG Facial Expression dataset (Aifanti, Papachristou, and Delopoulos 2010) records natural facial expression videos of $86$ subjects sitting in front of one camera. This dataset has no ground truth of the displacement field, which we use for unsupervised learning. $200$ videos with $30$ frames are randomly selected for training, and anther $100$ videos are sampled for testing.

6.2 Learned Gabor-like units with quadrature phase relationship

In this subsection, we show and analyze the learned units. The size of the filter is $16\times 16$, with a sub-sampling rate of $8$ pixels. Fig. 2(a) displays the learned units, i.e., rows of $W$, on V1Deform dataset. The units are learned with non-parametric $M(\delta)$, i.e., we learn a separate $M(\delta)$ for each displacement and $\delta(x)$ is discretized with an interval of $0.5$. V1-like patterns emerge from the learned units. Moreover, within each sub-vector, the orientations and frequencies of learned units are similar, while the phases are different and approximately follow a quadrature relationship, consistent with the observation of biological V1 simple cells (Pollen and Ronner 1981; Emerson and Huang 1997). Similar patterns can be obtained by using a parametric version of $M(\delta)$. See Supplementary for more results of learned filters, including filters learned with different dimensions of sub-vectors using different datasets. It is worthwhile to mention that the dimension of sub-vectors is not constrained to be $2$. V1-like patterns also merge when the dimension is $4$ or $6$.

To further analyze the spatial profiles of the learned units, we fit every unit by a two dimensional Gabor function (Jones and Palmer 1987): $h(x^{\prime},y^{\prime})=A\exp(-(x^{\prime}/\sqrt{2}\sigma_{x^{\prime}})^{2}-(y^{\prime}/\sqrt{2}\sigma_{y^{\prime}}))\cos(2\pi fx^{\prime}+\phi)$, where $(x^{\prime},y^{\prime})$ is obtained by translating and rotating the original coordinate system $(x_{0},y_{0})$: $x^{\prime}=(x-x_{0})\cos\theta+(y-y_{0})\sin\theta,y^{\prime}=-(x-x_{0})\sin\theta+(y-y_{0})\cos\theta$. The fitted Gabor patterns are shown in Fig. 2(b), with the average fitting $r^{2}$ equal to 0.96 ($\text{std}=0.04$). The average spatial-frequency bandwidth is 1.13 octaves, with range of $0.12$ to $4.67$. Fig. 2(c) shows the distribution of the spatial-frequency bandwidth, where the majority falls within range of $0.5$ to $2.5$. The characteristics are reasonably similar to those of simple-cell receptive fields in the cat (Issa, Trepel, and Stryker 2000) (weighted mean $1.32$ octaves, range of $0.5$ to $2.5$) and the macaque monkey (Foster et al. 1985) (median $1.4$ octaves, range of $0.4$ to $2.6$). To analyze the distribution of the spatial phase $\phi$, we follow the method in (Ringach 2002) to transform the parameter $\phi$ into an effective range of [math] to $\pi/2$, and plot the histogram of the transformed $\phi$ in Fig. 2(c). The strong bimodal with phases clustering near [math] and $\pi/2$ is consistent with those of the macaque monkey (Ringach 2002).

In the above experiment, we fix the size of the convolutional filters ($16\times 16$ pixels). A more reasonable model is to have different sizes of convolutional filters, with small size filters capturing high-frequency content and large size filters capturing low-frequency content. For fixed-size filters, they should only account for the image content within a frequency band. To this end, we smooth every image by two Gaussian smoothing kernels (kernel size $8$, $\sigma=1,4$), and take the difference between the two smoothed images as the input image of the model. The effect of the two smoothing kernels is similar to a bandpass filter so that the input images are constrained within a certain range of frequencies. The learned filters on V1Deform dataset are shown in Fig 3(a). We also fit every unit by a two dimensional Gabor function, resulting in an average fitting $r^{2}=0.83$ ($\text{std}=0.12$). Following the analysis of (Ringach 2002; Rehn and Sommer 2007), a scatter plot of $n_{x}=\sigma_{x}f$ versus $n_{y}=\sigma_{y}f$ is constructed in Fig. 3(b) based on the fitted parameters, where $n_{x}$ and $n_{y}$ represent the width and length of the Gabor envelopes measured in periods of the cosine waves. Compared to the sparse coding model (a.k.a. Sparsenet) (Olshausen and Field 1996, 1997), the units learned by our model have more similar structure to the receptive fields of simples cells of Macaque monkey.

We also quantitatively compare the learned units within each sub-vector in Fig. 3(c). Within each sub-vector, the frequency $f$ and orientation $\theta$ of the paired units tend to be the same. More importantly, most of the paired units differ in phase $\phi$ by approximately $\pi/2$, consistent with the quadrature phase relationship between adjacent simple cells (Pollen and Ronner 1981; Emerson and Huang 1997).

6.3 Inference of displacement field

We then apply the learned representations to inferring the displacement field $(\delta(x))$ between pairs of frames $({\bf I}_{t},{\bf I}_{t+1})$. To get valid image patches for inference, we leave out those displacements at image border ($8$ pixels at each side).

We use non-parametric $M(\delta)$ and the local mixing motion model (Eq. (10)), where the local support $\mathcal{S}$ is in a range of $[-4,+4]$, and $\text{d}x$ is taken with a sub-sampling rate of $2$. After obtaining the inferred displacement field $(\delta(x))$ by the learned model, we also train a CNN model with ResNet blocks (He et al. 2016) to refine the inferred displacement field. In training this CNN, the input is the inferred displacement field, and the output is the ground truth displacement field, with least-squares regression loss. See Supplementary for the architecture details of the CNN. For biological interpretation, this refinement CNN is to approximate the processing in visual areas V2-V6 that integrates and refines the motion perception in V1 (Gazzaniga, Ivry, and Mangun 2002; Lyon and Kaas 2002; Moran and Desimone 1985; Born and Bradley 2005; Allman and Kass 1975). We learn the models from the training sets of V1Deform and V1FlyingObjects datasets respectively, and test on the corresponding test sets. We also test the model learned from V1FlyingObjects on MPI-Sintel Final, whose frames are resized to $128\times 128$.

Table 1 summarizes the average endpoint error (AEE) of the inferred results. We compare with several baseline methods, including FlowNet 2.0 and its variants (Dosovitskiy et al. 2015; Ilg et al. 2017). For baseline methods, we test the performance using two models: one retrained on our datasets (‘trained’) and the released model by the original authors which are pre-trained on large-scale datasets (‘pre-trained’). Note that for MPI-Sintel, a pre-trained baseline model gives better performance compared to the one trained on V1FlyingObjects, probably because these methods train deep and complicated neural networks with large amount of parameters to predict optical flows in supervised manners, which may require large scale data to fit and transfer to different domains. On the other hand, our model can be treated as a simple one-layer auto-encoder network, accompanied by weight matrices representing motions. As shown in Table 1, our model has about $88$ times fewer parameters than FlowNet 2.0 and $21$ times fewer parameters than the light FlowNet2-C model. We achieve competitive performance compared to these baseline methods. Fig. 4 displays several examples of the inferred displacement field. Inferred results from the FlowNet 2.0 models are shown as a qualitative comparison. For each dataset, we show the result of FlowNet 2.0 model with lower AEE between the pre-trained and trained ones.

6.4 Unsupervised learning

We further perform unsupervised learning of the proposed model, i.e., without knowing the displacement field of training pairs. For unsupervised learning, we scale the images to size $64\times 64$. The size of the filters is $8\times 8$, with a sub-sampling rate of $4$ pixels. Displacements at the image border ($4$ pixels at each side) are left out. We train the model on MUG Facial Expression and V1FlyingObjects datasets. Fig. 5 shows some examples of inferred displacement fields on the testing set of MUG Facial Expression. The inference results are reasonable, which capture the motions around eyes, eyebrows, chin, or mouth. For the model trained on V1FlyingObjects, we test on the testing set of V1FlyingObjects and MPI-Sintel. Table 2 summarizes the quantitative results. We include comparisons with several baseline methods for unsupervised optical flow estimation: Unsup (Jason, Harley, and Derpanis 2016) and UnFlow (Meister, Hur, and Roth 2018) and its variants, which are also trained on V1FlyingObjects. The proposed model achieves better performance compared to baseline methods. See Supplementary for qualitative comparisons and more inference results.

6.5 Multi-step frame animation

The learned model is also capable of multi-step frame animation. Specifically, given the starting frame ${\bf I}_{0}(x)$ and a sequence of displacement fields $\{\delta_{1}(x),...,\delta_{T}(x),\forall x\}$, we can animate the subsequent multiple frames $\{{\bf I}_{1}(x),...,{\bf I}_{T}(x)\}$ using the learned model. We use the model with local mixing. We introduce a re-encoding process when performing multi-step animation. At time $t$, after we get the next animated frame ${\bf I}_{t+1}$, we take it as the observed frame at time $t+1$, and re-encode it to obtain the latent vector $v_{t+1}$ at time $t+1$. Fig. 6 displays two examples of animation for 6 steps, learned with non-parametric version of $M$ on V1Deform and V1FlyingObjects. The animated frames match the ground truth frames well. See Supplementary for more results.

6.6 Frame interpolation

Inspired by the animation and inference results, we show that our model can also perform frame interpolation, by combining the animation and inference together. Specifically, given a pair of starting frame ${\bf I}_{0}$ and end frame ${\bf I}_{T}$, we want to derive a sequence of frames $({\bf I}_{0},{\bf I}_{1},...,{\bf I}_{T-1},{\bf I}_{T})$ that changes smoothly. Let $v_{0}(x)=W{\bf I}_{0}[x]$ and $v_{T}(x)=W{\bf I}_{T}[x]$ for each $x\in D$. At time step $t+1$, like the inference, we can infer displacement field $\delta_{t+1}(x)$ by steepest descent:

[TABLE]

Like the animation, we get the animated frame ${\bf I}_{t+1}$ by decoding $\hat{v}_{t+1}(x,\delta_{t+1}(x))$, and then re-encode it to obtain the latent vector $v_{t+1}(x)$.

The algorithm stops when ${\bf I}_{t}$ is close enough to ${\bf I}_{T}$ (mean pixel error $<10$). Fig. 7 shows four examples, learned with non-parametric $M$ on V1Deform and V1FlyingObjects. For $96.0\%$ of the testing pairs, the algorithm can accomplish the frame interpolation within $10$ steps. With this algorithm, we are also able to infer displacements larger than the acceptable range of $\delta$ by accumulating the displacements along the interpolation steps. See Supplementary for more results.

6.7 Ablation study

We perform ablation studies to analyze the effect of two components of the proposed model: (1) dimensionality of sub-vectors and (2) sub-sampling rate. Besides comparing the average endpoint error (AEE) of motion inference, we also test if the learned model can make accurate multi-step animation of image frames given the sequence of displacement fields. Per pixel mean squared error (MSE) of the predicted next five image frames is reported. Table 3 summarizes the results learned from V1Deform dataset. The dimensionality of sub-vectors controls the complexity of the motion matrices, which is set to a minimum of $2$ in the experiments. As the dimensionality of sub-vectors increases, the error rates of the two tasks decrease first and then increase. On the other hand, sub-sampling rate can be changed to make the adjacent image patches connect with each other more loosely or tightly. As shown in Table 3, a sub-sampling rate of $8$, which is half of the filter size, leads to the optimal performance.

7 Conclusion

This paper proposes a simple representational model that couples vector representations of local image contents and matrix representations of local motions, so that the vector representations are equivariant. Unlike existing models for V1 that focus on statistical properties of natural images or videos, our model serves a direct purpose of perception of local motions. Our model learns Gabor-like units with quadrature phases. We also give biological interpretations of the learned model and connect it to the spatiotemporal energy model. It is our hope that our model adds to our understanding of motion perception in V1.

Our motion model can be integrated with the sparse coding model. For sparse coding, we can keep the top-down decoder of the tight frame auto-encoder for each image frame, and impose sparsity on the number of sub-vectors that are active. For motion, we then assume that each sub-vectors is transformed by a sub-matrix that represents the local displacement.

This paper assumes linear decoder for image frames and matrix representation of local motion. We can generalize the linear decoder to a neural network and generalize the matrix representation to non-linear transformation modeled by non-linear recurrent network.

In our future work, we shall study the inference of ego-motion, object motions and 3D depth information by generalizing our model based on equivariant vector representations and their transformations. We shall also apply our model to stereo in binocular vision.

Acknowledgement

The work is supported by NSF DMS-2015577, ONR MURI project N00014-16-1-2007, DARPA XAI project N66001-17-2-4029, and XSEDE grant ASC170063. We thank Prof. Tai Sing Lee for sharing his knowledge and insights on V1.