Modelling Human Visual Motion Processing with Trainable Motion Energy Sensing and a Self-attention Network

Spatiotemporal separable Gabor filter: Since we aimed at an image-computable model, the input is a sequence of grayscale images $\mathbf{S}(\mathbf{p},\mathbf{t})$ for all spatial positions $\mathbf{p}=(x,y)$ within the image domain $\mathbf{\Omega}$ and for all times $\mathbf{t>0}$. The goal of the first-stage neuron is to capture local motion energy at a specific spatiotemporal frequency, which is associated with the function of a direction-selective neuron in the V1 cortex. The responses of neurons can be modeled as 3D Gabor filters [30, 31]. Gabor filters are known to be optimal in the sense that they achieve maximal resolution in both the spatiotemporal and associated frequency domains[32]. To save computational complexity, we decomposed 3D spatiotemporal filters into filters separable in space and time. The spatial component is described by 2D Gabor filters $\mathcal{G(\cdot)}$, and the temporal component $\mathcal{T(\cdot)}$ is described by a 1D sinusoidal function with exponential decay. Specifically, given $x^{\prime}=x\cos\theta+y\sin\theta$ and $y^{\prime}=-x\sin\theta+y\cos\theta$ within the receptive field, the impulse responses of spatial and temporal complex filters are defined as:

$$\left\{\begin{array}[]{lr}\mathcal{G}(x,y;{\color[rgb]{1,0,0}f_{s}},{\color[rgb]{1,0,0}\theta},{\color[rgb]{1,0,0}\sigma},{\color[rgb]{1,0,0}\gamma})=\exp\left(-\frac{x^{\prime 2}+\gamma^{2}y^{\prime 2}}{2\sigma^{2}}\right)\exp\left(i\left(2\pi{f_{s}}{x^{\prime}}\right)\right),s.t.~{}\{x,y~{}|(x^{2}+y^{2}\leq R^{2})\}&\\ \mathcal{T}\left(t;{\color[rgb]{1,0,0}f_{t}},{\color[rgb]{1,0,0}\tau}\right)=\exp{\left(-\frac{t}{\tau}\right)}\exp{(2\pi i\left(f_{t}t\right))},s.t.~{}\{t~{}|~{}0\leq t<T\}&\end{array}\right.$$

(1)

The parameters in red are designed to be trained to fit the dataset, where $f_{s}$ and $f_{t}$ denote the spatiotemporal frequency tunings of the filter with the preferred speed $v\propto\frac{f_{t}}{f_{s}}$, and $\theta$ determines the preferred moving orientation; $\sigma$ and $\gamma$ control the shape of the Gabor filter, and $\tau$ controls the degree of attenuation of the temporal impulse response. All parameters are subject to certain numerical constraints, such as $\theta$ being limited to [0,$2\pi$) to avoid redundancy; $f_{s}$ and $f_{t}$ are limited to less than 0.25 pixels per frame to avoid spectrum aliasing, etc. The response of a simple direction-selective cell $L_{n}$ to a video stimuli $\mathbf{S}(\mathbf{p},\mathbf{t})$ can be computed via separate convolutions:

$$L_{n}(x,y,t;\Theta)=(\mathbf{S}*\mathcal{G})*\mathcal{T}=\iiint\mathbf{S}(\mathcal{X},\mathcal{Y},\mathcal{T})\cdot\mathcal{G}_{n}(x-\mathcal{X},y-\mathcal{Y})\cdot\mathcal{T}_{n}(t-\mathcal{T})d(\mathcal{X},\mathcal{Y})d\mathcal{T}+{\color[rgb]{1,0,0}\alpha_{1}}$$

(2)

where $\alpha_{1}$ can be learned as spontaneous firing rates. Further, local motion energy is captured by a phase-insensitive complex cell in the V1 cortex, which computes the squared summation of the response from a pair of simple V1 cells with approximately orthogonal spatiotemporal receptive fields[24]. We denote the pair of orthogonal (even and odd) simple V1 cells as:

$$\left\{\begin{array}[]{lr}L_{n}^{o}(x,y,t;{\color[rgb]{1,0,0}\Theta})=(\mathbf{S}*{\rm Im}[\mathcal{G}])*{\rm Re}[\mathcal{T}]+(\mathbf{S}*{\rm Im}[\mathcal{G}])*{\rm Im}[\mathcal{T}]&\\ L_{n}^{e}(x,y,t;{\color[rgb]{1,0,0}\Theta})=(\mathbf{S}*{\rm Re}[\mathcal{G}])*{\rm Re}[\mathcal{T}]-(\mathbf{S}*{\rm Im}[\mathcal{G}])*{\rm Im}[\mathcal{T}],\end{array}\right.$$

(3)

where ${\rm Re}(\cdot)$ and ${\rm Im}(\cdot)$ extract the real and imaginary parts of a complex number; $*$ denotes convolution operations. Then, the response of a complex cell $L_{n}^{c}$ is obtained from a combination of the quadrature pair of the simple cells using the motion energy formulation:

$$L_{n}^{c}(x,y,t;{\color[rgb]{1,0,0}\Theta})=\left(L_{n}^{o}(x,y,t;{\color[rgb]{1,0,0}\Theta})\right)^{2}+(L_{n}^{e}(x,y,t;{\color[rgb]{1,0,0}\Theta}))^{2}$$

(4)

Multi-scale Wavelet Processing: The convolution kernel of our spatial filter has a fixed size of $15\times 15$, which imposes a physical limitation on the receptive field of each unit. To enhance the flexibility of the receptive field size, we employed a multi-scale processing strategy, as shown in Fig. 1 (A). Specifically, we constructed an image pyramid consisting of eight images scaled linearly from $H\times W$ to $\frac{H\times W}{16}$. In total, 256 independent complex cells are deployed across different scales, with the lower scale having a larger receptive field and a preference for a faster motion speed. This approach is computationally efficient for capturing large-scale displacements in engineering applications[33, 34], and enables the representation of different groups of cells sensitive to short- and long-distance motions [35]. The N complex cells $\{L_{n}^{c}\}_{i}^{N}$ capture motion energy on multiple scales. We applied energy normalization to each cell to ensure consistent energy levels:

$$\small\hat{L}_{n}^{c}(t)=\frac{{{\color[rgb]{1,0,0}K_{1}}}L_{n}^{c}(t)}{\sum_{\mathrm{i=1}}^{\mathrm{N}}L_{\mathrm{i}}^{c}(t)+{{\color[rgb]{1,0,0}\sigma_{1}}}},$$

(5)

where $\sigma_{1}$ is the semi-saturation constant of the normalization, and $K_{1}>0$ determines the maximum attainable response. We interpret the response, denoted by $\hat{L}_{n}(t)$, as the model equivalent of a post-stimulus time histogram (PSTH), which is a measure of the neuron’s firing rate. Physiologically, such responses could be computed via inhibitory feedback mechanisms[36, 37]. Considering the spatial arrangement of images, it is expected that motion energy responses should exist at each spatial location generated by the same complex cell groups. Bilinear interpolation is used to resize the multi-scale motion energy into the same spatial size = $\frac{H\times W}{8}$. In the context of DNNs, this is mainly to balance the trade-off between the spatial resolution and computational overhead, and the final output of the first stage is a 256-channel feature map $\mathbf{E}\in\mathbb{R}^{\frac{H}{8}\times\frac{W}{8}}$ that captures the underlying local motion energy, which partially characterizes the cellular patterns of the V1 cortex in a computational manner[24].

The receptive field size of first-stage neurons limits their ability to capture only local motion. Advanced spatial integration is an essential requirement for a motion perception system to solve the aperture problem[7]. Spatial integration could involve several mechanisms, such as object recognition and segmentation[38], depth estimation[39], contour structure inference[40], etc., which rely on substantial prior information that may go beyond the capabilities of classical modeling methods. DNNs, with their massive number of parameters and flexibility, offer a suitable solution. However, the spatial integration of local motion requires more flexible connectivity relations than what can be achieved by general convolutions, which are limited by their local receptive field[41]. To overcome this limitation, we propose a computational model based on the attention mechanism and recurrent processing to model the function of motion integration.

Constructing a Graph Structure: First, to establish more flexible connections between cells, we discard the Euclidean space structure within the image and construct topological spaces with an undirected weighted graph $\mathbf{G}=\{\mathbf{V},\mathbf{A}\}$, where $\mathbf{V}$ is the set of nodes; $\mathbf{A}$ represents the adjacent matrix based on the graph. Each spatial location $p(i,j)$ is treated as a node, and the feature of each node is derived from the whole set of local motion energies $\mathbf{E}(i,j)=\{\hat{L}_{n}^{c}(t)\}_{i=1}^{256}$. The connection between any pair of nodes is computed using a specific distance metric, and strong connection relationships are formed between nodes whose local motion energy patterns are similar.

Specifically, given a reshaped feature map $\mathbf{F}\in\mathbb{R}^{HW\times 256}$, the features at each of its locations are first locally embedded into the graph space, as $\varphi(\mathbf{F})=\operatorname{GELU}(\mathbf{F}\cdot\mathbf{W}_{\Theta 1})\cdot\mathbf{W}_{\Theta 2}$ , where $\mathbf{W}_{\Theta}\in\mathbb{R}^{256\times 256}$ is a group of trainable parameters. The distance between any pair of nodes $(i,j)$ is calculated by the cosine similarity, which is similar to the self-attention mechanism in current transformer structures[26, 42, 27]. We employ the adjacency matrix $\mathbf{A}\in\mathbb{R}^{HW\times HW}$ to represent the connectivity of the whole topological space, and $\mathbf{A}$ is a symmetric semi-positive definite matrix defined as:

$$\mathbf{A}(i,j)=\mathbf{A}(j,i)=\frac{\varphi(\mathbf{F})_{i}\cdot\varphi(\mathbf{F})_{j}}{\|\varphi(\mathbf{F})_{i}\|\|\varphi(\mathbf{F})_{j}\|}.$$

(6)

We perform exponential scaling of the connections between graphs using the matrix $\mathbf{A}$, given by $\exp(\mathbf{A}s)$, where $s$ is a learnable scalar restricted to the range (0,10) to avoid overflow. The smaller $s$, the smoother the connections across nodes and vice versa. Finally, a symmetric normalization operation is utilized to balance the energy, given by $\mathbf{A}:=\mathbf{D}^{-\frac{1}{2}}\exp(s\mathbf{A})\mathbf{D}^{-\frac{1}{2}}$, where $\mathbf{D}$ is the degree matrix with $\mathbf{D}=\operatorname{diag}\left(\left\{\sum_{j}\exp(s\mathbf{A}_{i,j})\right\}_{i-1}^{n}\right)$. As such, an energy-normalized undirected graph structure is constructed, as illustrated in the top side of Fig. 2 (B). Intuitively, this adjacency matrix represents the neuron’s affinity or connectivity within the space, with strong and global connections built across neurons with related motion responses.

Recurrent Integration Processing: Recurrent neural networks (RNNs) are often used to simulate neurons in the brain, as they are flexible in modelling temporal dependencies and feedback loops, which are fundamental aspects of neural processing in the brain[43]. We use a recurrent network, rather than multiple feedforward blocks, to simulate the process of local motion signals being gradually integrated into MT and eventually converging to a stable state.

As shown in Fig. 2 (A), the local motion energy from the first stage is divided into two recurrent streamlines. One is the motion energy $\small\mathbf{E}\in\mathbb{R}^{H\times W\times 256}$, which is continuously updated in the loop, while the other is embedded in the attention space to generate the graph adjacency matrix to control motion integration, denoted as $\small\mathbf{F}\in\mathbb{R}^{H\times W\times 256}$. In each iteration, the adjacency matrix is first constructed using $\mathbf{F}$. Subsequent motion integration is achieved through a simple matrix multiplication, which is computationally similar to the information propagation mechanism in transformers[26, 42] and can also be considered a simplified version of graph convolution[44]. The integrated motion information is passed through two independent Conv-GRU blocks to update the motion energy $\mathbf{E}$ and feature $\mathbf{F}$, respectively. The Conv-GRU represents a gated recurrent unit[45] implemented in a convolutional manner, and we adopt a spatio-temporal separable approach following RAFT[14]. The motion integration process approximates the ideal final convergence state of the motion energy $\mathbf{E}_{k}\rightarrow\mathbf{E}^{*}$ through recurrent iteration.

Decoding the Motion Flow: We adopt the same strategy of decoding the 2D optical flow in each iteration. Initially, the integrated motion is $\mathbf{E}$ projected back to the energy space with positive value using a square operation, followed by an energy normalization operation: $\hat{\mathbf{E}}(i,j)=\{\small{K_{2}}\mathbf{E}^{2}(i,j)\}/\{\small\sum_{{i,j}}^{\mathrm{H}\mathrm{W}}\mathbf{E}^{2}(i,j)+{\sigma_{2}}^{2}\}.$ The resulting response $\hat{\mathbf{E}}\in\mathbb{R}^{H\times W\times 256}$ is interpreted as a post-stimulus time histogram and the motion energy is constrained to the same energy space as the local motion energy from stage I, as we further design an identical flow decoder to project the energy of each spatial location into the optical flow. The structure of the flow decoder is a nonlinear mapping process consisting of multiple $1\times 1$ convolution blocks with residual connections, which are referred to in several current advanced optical flow estimation models[33, 46]. Additionally, a convex upsampling strategy[14] is employed to restore the optical flow’s original resolution. The entire architecture of stage II is illustrated in Fig. 2 (A, B).

Current methods for estimating optical flow using deep neural networks (DNNs) can be categorized as unsupervised/self-supervised and supervised learning approaches. While unsupervised learning methods are intuitively similar to creatures’ interaction with the world, most current methods based on differentiable image warping[47, 48, 49] still try to approximate physical motion ground truth. Therefore, we adopt a supervised learning approach in this work, which is more straightforward as recent research suggests that human perception of motion is reasonably similar to physical GT[16]. However, our primary goal lies in evaluating how well the model approximates human motion perception rather than its accuracy in predicting GT. To train and evaluate the model, we construct a dataset containing various natural and artificial motion scenes. Specifically, we incorporate the Sintel benchmark[50], the DAVIS[51] dataset with pseudo-labels generated by FlowFormer[52], as well as self-created multi-frame datasets with non-texture motions and drifting grating motion. Including simple motion stimuli and drifting gratings allows the model to generalize under different non-texture conditions while providing a potential slow-world Bayesian prior[9]. The model is first pre-trained with simple motion and subsequently fine-tuned on complex natural scenes to facilitate convergence[12]. More specific training details can be found in the supplementary materials.

Directional Tuning: Some V1 and MT neurons respond selectively to a specific range of motion directions. To investigate the directional tuning of the units in our model, we utilized in silico neurophysiology to measure the activation patterns of 256 units in response to drifting Gabor and plaid stimuli. A plaid consists of two superimposed drifting Gabors, as illustrated in Fig. 4 (C). Analysis of the data revealed three distinct groups of units based on their partial correlation to Gabor and plaid stimuli: 1) Component cells, which always respond to the direction of a Gabor component; 2) Pattern cells, which respond to the integrated (coherent) motion direction of plaid; and 3) Unclassified cells, which do not show a clear preference for either component or pattern motions, as illustrated in Fig. 4 (C).

The distribution of these cell types is not uniform, with component cells being more commonly found in the primary visual cortex (V1), while pattern cells are more often observed in the MT and MST regions[29]. In agreement with this, as illustrated in Fig. 4 (A) (B), our model demonstrates that in the first stage, most units tend to be component cells, whereas the number of pattern and unclassified cells increases in the second stage. In addition, we employed the maximizing activation method [56] to reversely render their cellular preferences, showing that the second-stage unclassified units respond to more complicated motion patterns consisting of both central and background motions, as presented at the bottom of Fig. 4 (B). This suggests that the classical classification of motion neurons into component and pattern cells might be insufficient to characterize the motion integration properties of these neurons.

Spectral Receptive Field and Speed Tuning: We test the spectral (spatial-frequency-vs-temporal-frequency) receptive field of the model’s units using a combination of drifting Gabors with different spatiotemporal frequencies. Two-dimensional oriented Gaussian contours are used to depict the receptive field of the 256 cells, fitted by minimizing the least square error. The results are shown in Fig 4 (E) . Visually, the distribution of the receptive field tilt angle spreads from horizontal/vertical directions to oblique directions (Fig. 4, E). This indicates that the units in the stage II have significant speed tuning compared to the stage I. Speed tuning is a characteristic of higher-order visual motion neurons[28] and is often found in the MT area[57]. This tendency can be seen from the distribution of the partial correlation between the actual receptive field and its speed prediction/independent prediction[55], as demonstrated in Fig. 4 (F). Our two-stages process shows a degree of consistency with the change in mammalian neural distribution from the V1 to MT area.

Fourier motion: In the "missing fundamental illusion"[60], as shown in Fig. 5 (C), when the first spatial harmonic is removed from a square-wave grating with a quarter-cycle shift, the perceived motion direction appears reversed. Our model, whose first stage estimates motion from the Fourier components[61], can predict this reversal. In contrast, computer vision models designed to infer optical flow based on structural correspondence do not exhibit this bias.

Motion Integration: The direction of 1D motion stimuli, such as drifting Gabors, is ambiguous due to the aperture problem. When presented alone, they are perceived to move in orthogonal directions of stripes. When superimposed with other 1D motion components in different orientations, 2D directions consistent with both components are perceived. Pattern motion neurons in MT may contribute to this phenomenon. The second stage of our model can explain this as well, as shown in Fig. 5. Psychophysics also demonstrates motion integration across space. Fig. 5 (A) shows a psychophysical stimulus consisting of drifting Gabors [58]. These Gabors have a variety of local directions and speeds, yet all of them are consistent with one global 2D motion (downward in this case). When viewed as a whole, humans do perceive coherent downward motion. Stage II of our model predicts both the perceived direction and speed of the global Gabor motion.

Fig. 5 (B) compares spatial motion integration between 1D Gabor motion (left) and 2D plaid motion (right). Humans are able to perceive global downward motion only in the former case: The heat maps depict how units with high activity establish long-distance connections to resolve the aperture problem when subjected to Gabor (ambiguous motion) stimuli. In contrast, the plaid stimuli (determined motion) suppress these long-distance connections. In the latter case, local integration of motion signals takes priority over global integration. Once the local ambiguity is resolved, the global integration process is suppressed. Our model can predict such adaptive motion pooling in human visual processing [58]. Furthermore, the barber pole illusion demonstrates how locally ambiguous 1D motion is affected by the shape of the moving area [62]. Specifically, as the height-width ratio of the visual area varies, human perception of direction shifts from oblique to vertical [59] (Fig. 5, D). Our model can predict the shift in perceived direction in stage II, showing its ability to integrate motion signals with boundary orientations. For more video demonstrations, such as the reverse phi illusion[63], please see the supplementary material.

Comparison of Complex Natural Scenes: The proposed model can effectively handle natural scenes, as demonstrated in Fig. 6 (A). Natural stimuli contain diverse spatiotemporal frequency components, leading to complex activation patterns in Stage I. From the decoded flow field, the motion of stage I is limited to localized areas. For example, local motion cannot be found in untextured road areas. This situation necessitates long-range interactions with the surrounding spatial context, which our recurrent integration process in stage II effectively accomplishes. It is evident from the affinity heat map (right side of subfigure A) that object and background areas are clearly segregated in stage II. This suggests that the integration mechanism based on attention has the potential to combine motion integration and object segmentation into a single framework. These two processes are considered highly relevant in the human visual system [64].

We used psychophysically measured optical flow from the Sintel dataset[16] as a benchmark for naturalistic scene flow perceived by humans. Our model was compared to several optical flow estimation methods used in computer vision, including classical algorithms like Farneback, as well as SOTA DNN-based models. We evaluated officially released DNN models that utilize a wide range of inference structures, such as multi-scale inference[11, 12], spatial recurrent models[14], graph reasoning[65], and vision transformers[15].

As shown in Table 1, we computed both the Pearson correlation coefficients and vector endpoint errors (EPEs) between the model prediction and human response, or ground truth (GT). Additionally, we examined the partial correlation between humans and models while controlling the impact of GT:

$$\rho_{\text{model }}=r_{\text{resp model }\cdot GT}=\frac{r_{\text{resp model }}-r_{\text{resp }GT}\cdot r_{\text{model }GT}}{\sqrt{1-r_{\text{resp }GT}^{2}}\sqrt{1-r_{\text{model }GT}^{2}}},$$

(7)

where $r$ is the Pearson correlation. This measure is critical in validating the models’ ability to capture the characteristics of the human response, as any model could appear to have a high correlation with the human response by just approximating the GT due to a high correlation between the human response and GT[16].

Quantitatively, the proposed model outperformed all compared models in terms of partial correlation. The RAFT-val in the table is the RAFT framework trained on our dataset as a validation, and indicates that our mixed training set also improves the explanatory power of the human response. Fig. 6 (C) shows that our model significantly improves the partial correlation across iterations in Stage II, indicating that the proposed recurrent motion integration architecture can generate more human-like deviations from GT, which is a trend not present in a similar recurrent network, RAFT. In Fig. 6 (D), one can directly see that our model prediction is more similar to human-perceived flow than the GT flow.

Table 1 also shows the results of three biologically-inspired models. FFV1MT[67] is a model capable of computing dense optical flow using direct decoding of the Simoncelli & Heeger V1-MT mechanism. The other two models are modified versions of MotionNet[20] and DorsalNet[23], recently proposed DNN-based models for the explanation of neural responses to visual motion stimuli. Since the original models are designed to recognize global motion only, we tested general multi-layer 3D CNN (a core component of MotionNet and DorsalNet) with residual connections, and a pre-trained DorsalNet with frozen parameters, with a linear flow decoder to compute dense flow, trained on natural dense optical flow datasets. Our model outperforms these biologically-inspired models in predicting human responses. The low performances of FFV1MT model and the modified DorsalNet also suggest that accurately estimating dense optical flow is a challenging task, requiring specific design considerations to address complex and long-range spatial interactions, large jumps, and boundary effects, the complexities of which are not adequately captured by simple mechanisms.