Insect-Inspired Self-Motion Estimation with Dense Flow Fields—An Adaptive Matched Filter Approach

The control of self-motion is a basic, but complex task for both technical and biological systems. Various algorithms have been proposed that allow the estimation of self-motion from the optic flow on the eyes. We show that two apparently very different approaches to solve this task, one technically and one biologically inspired, can be transformed into each other under certain conditions. One estimator of self-motion is based on a matched filter approach; it has been developed to describe the function of motion sensitive cells in the fly brain. The other estimator, the Koenderink and van Doorn (KvD) algorithm, was derived analytically with a technical background. If the distances to the objects in the environment can be assumed to be known, the two estimators are linear and equivalent, but are expressed in different mathematical forms. However, for most situations it is unrealistic to assume that the distances are known. Therefore, the depth structure of the environment needs to be determined in parallel to the self-motion parameters and leads to a non-linear problem. It is shown that the standard least mean square approach that is used by the KvD algorithm leads to a biased estimator. We derive a modification of this algorithm in order to remove the bias and demonstrate its improved performance by means of numerical simulations. For self-motion estimation it is beneficial to have a spherical visual field, similar to many flying insects. We show that in this case the representation of the depth structure of the environment derived from the optic flow can be simplified. Based on this result, we develop an adaptive matched filter approach for systems with a nearly spherical visual field. Then only eight parameters about the environment have to be memorized and updated during self-motion.

available to animals which have to rely on other means to gain information about their position and self-motion. A direct method to measure self-motion for a walking artificial or biological agent is counting the steps or, in the case of a wheeled vehicle, to monitor the turns of the wheels. In contrast, most flying agents rely on their visual system to solve this task.
The visual system of an artificial or biological agent obtains information about self-motion from pixel shifts in the retinal image over time. These pixel shifts can be described by vectors, the optic flow vectors. The flow vectors depend on both the rotational and translational components of self-motion as well as on the viewing direction. Moreover, for the translational component it also depends on the distance to objects in the environment.
For small translations and rotations, the flow vector for viewing directiond i is given by (see [1] for derivation)p i ðt;r; m i Þ ¼ Àm i ðt À ðt Ád i Þd i Þ Àr Âd i ¼ Àm i ðd i Ât Âd i Þ Àr Âd i ; where μ i is the inverse distance ("nearness") to the object seen in directiond i ,t is the translation vector, andr is the rotation vector (defining a rotation of angle r ¼jr j around the axis given byr= jr j). According to Eq (1), the flow vectorp i is perpendicular to the corresponding viewing directiond i . (We use 3d vectors to represent optic flow vectors. Otherwise one would have to define a tangential plane for every viewing direction.) There are two principally different ways to use optic flow information for self-motion estimation. One way is to identify features in the retinal image at one point in time and find the same features at the next time point in order to compute their displacement. Several technical estimation methods for self-motion are based on these feature correspondences [2,3]. They rely on a small number of corresponding image points and have to concern about outliers. Such estimation methods are widely used in the technical literature to determine the movement and/or the calibration of a camera [4]. When the self-motion steps are small or the frame rates are high an alternative way to extract self-motion information is possible. Instead of extracting features in the image the local pixel shifts on the retina, called optical flow, produced by the self-motion of the agent is determined through spatiotemporal intensity correspondences in the pattern. This can be done by a gradient-based detector like the Lucas-Kanade detector [5], which compares spatial and temporal derivatives, or by a biologically inspired detector, like the elementary movement detector of the correlation type [6,7], which uses spatiotemporal auto-correlation signals.
Here we propose a new adaptive approach which combines the advantages of two methods for self-motion estimation based on optical flow, the matched filter approach (MFA) proposed by Franz and Krapp [8] and an algorithm proposed by Koenderink and van Doorn (KvD algorithm) [1]. The MFA estimates self-motion by using linear filters, so called matched filters. Matched filters have the structural form of the pattern they are meant to detect [9] and are the optimal detectors for patterns, which are disturbed by Gaussian errors. In this case the linear filters of the MFA resemble ideal flow fields. Franz and Krapp [8] introduced six filters of this type for the six self-motion components, three for translation and three for rotation. Each of these six filters was tuned to one of the flow fields generated by the six self-motion components, although in general the filters react also to flow generated by the other self-motion components. There is one exception: For a flow field which covers the whole viewing sphere and for isotropic distances, i.e. in the center of a sphere, Borst and Weber [10] showed that model neurons acting as such linear filters are not influenced by other flow fields. To eliminate the sensitive only to low order depths functions, which are the dipole and quadrupole moments of the depth structure.

Results
A major objective of this study is to show that two well-established self-motion estimators are mathematically equivalent: The MFA equals one iteration step of the KvD algorithm when the inverse distances to objects in the environment are assumed to be known. To achieve this, we derive the MFA in an alternative way.
We then show that the KvD approach leads to an biased estimator in the general case when distances are unknown and have to be estimated together with the self-motion parameters. We present a modification of the KVD iteration equation that removes the bias and derive an adaptive MFA from this corrected version, which includes a simple but, with respect to self-motion estimation, complete depth model.
Before dealing with these topics, both the basic equations underlying the MFA approach and the KvD algorithm need to be introduced.

The matched filter approach
In the original MFA [8,11,33,34] the depth structure of the environment is not determined from the current flow field but described statistically with a fixed distribution that is assumed to be known. The first statistical parameter considered in [11] is the average inverse distance m i , which is measured in every viewing direction over a number of learning flights in different environments. The variability of the distances is given by the covariance matrix C μ . Secondly, the noise in the flow measurement is determined for each viewing direction n i where a zero mean is assumed. The noise values are combined in the covariance matrix C n . The third statistical parameter is the distribution of the translations t. It is assumed that the agent does not translate in every possible direction with the same probability. The corresponding statistical parameter is the covariance matrix C t .
An optic flow vectorp i has only two degrees of freedom because it is the projection of object motion on the retina and thus orthogonal to the corresponding viewing directiond i . To consider only these degrees of freedom Franz et al. [11] introduce a local two-dimensional vector space for each viewing directiond i which is orthogonal to the directiond i : whereũ andṽ are the basis vectors of the new vector space. The values x i and y i represent the two degrees of freedom ofp i . The measured vectorp i consists of the true optic flow vectorp 0 i and an additive noise n i . In [11] the weights W for the matched filters which are multiplied with the optic flow componentsx (wherex is a 2N dimensional vector containing all flow components x i and y i , i = 1, 2, . . ., N) to estimate the six self-motion componentsỹ est , are derived by a least-square principle: whereỹ are the true self-motion components. The weight matrix that minimizes the error e is: The covariance matrix C combines the covariance matrices C μ , C n and C t . The matrix F is given by where m i is the average or expected inverse distance for direction d i . The introduction of the matched filter approach is kept short because an alternative derivation of this approach is introduced in section 2.3, which is more suitable for the comparison with the KvD algorithm.

The Koenderink-van-Doorn (KvD) algorithm
As described by Koenderink and van Doorn [1], a straightforward approach for estimating the self-motion parameters is to find, in accordance with Eq (1), a translation vectort, a rotation vectorr and inverse distances {μ i } i = 1, 2, . . ., N that minimize the mean squared error between the theoretical optical flow vectors according to Eq (1), fp i ðt;r; m i Þg i¼1;2;...;N and the measured optical flow vectors fp i g i¼1;2;...;N :

Alternative derivation and properties of the coupling matrix in the MFA
To compare the two self-motion estimators one needs another mathematical form of the coupling matrix and the filters as given by the MFA (see 2.1). This can be achieved in two different ways. One can either transform the original Eqs (4), (6) and (7) of the MFA or derive the MFA in an alternative manner and show the equivalence to the original MFA afterward. Here, the second way is chosen. The matched filters will be derived according to the theory of optimal filters. The coupling of the estimated self-motion parameters will then be determined by inserting the filters in the flow Eq (1). The optimal weights in the MFA depend on the statistics of the distances, of the noise and of the preferred translations. Here we assume that nothing is known about these distributions and thus consider the simplest case: The noise values are set to the same value independent of the viewing direction. We assume no preference for specific translation directions. The average inverse distances m i are regarded as known. The theory of optimal filters states that for a uniform Gaussian noise the best linear filter is a matched filter which has the same form as the pattern it has to detect [9]. Therefore, the templates for estimating translational flow fields have the form of a translational flow field.
We have three translational templates one for each possible direction of translation represented by the three basis vectorsẽ a (a = 1,2,3). Similarly, we have three rotational templates The three translational and three rotational templatesT t a andT r a will be called "standard templates".
The scalar productTÁp D E of a flow fieldp and a templateT , where the brackets stand for the mean over all viewing directions, can be interpreted as the output of a specific model In general, the model neurons do not only react to the flow fields they are tuned to, but also to other flow fields. To solve this problem Franz et al. [11] introduced a matrix (F T C −1 F) −1 (see Eq (6)) which compensates for the coupling to other flow fields. The coupling between the self-motion estimates and therefore the coupling matrix are determined be inserting the templatesT t a andT r a in the flow Eq (1). For this, the translationt and the rotationr must be separated into their components,t ¼ t 1ẽ1 þ t 2ẽ2 þ t 3ẽ3 and r ¼ r 1ẽ1 þ r 2ẽ2 þ r 3ẽ3 , with the six self-motion parameters t 1 , t 2 , t 3 and r 1 , r 2 , r 3 .
Since the cross product is linear aã þ bb Âc ¼ aã Âc ð Þþbb Âc , the overall flow field is the sum of the six standard templatesT A (A = 1, 2, . . ., 6) weighted by the six selfmotion components θ A : Following our notation, the response of model neuron a A with corresponding templateT A to the flow fieldp (Eq (15)) is Combining the responses of all six model neurons in one equation gives where the vectorsã andỹ are considered as six dimensional vectors. Each entry of the 6 × 6 dimensional matrixM, the coupling matrix, can be seen as the generalized scalar product of two of the six standard templates:M We can also write the coupling matrix as where the indices t and r of the 3 × 3 sub-matrices indicate which templates are multiplied.
Using the inverse of the coupling matrix we can estimate the motion parametersỹ from the responses of the model neurons, The productT Áp D E and the multiplication withM À1 are linear transformations of the optical flow. One can therefore define new templates T 0 which include the linear transformation given by the matrixM À1 . Dahmen et al. [12] tested previously the self-motion estimation performance for different fields of view. The self-motion estimation performance even for error prone flow vectors is high, if the flow fields corresponding to different self-motion components differ essentially over the field of view. For a restricted field of view, for example a small region in front of the agent, upward translation cannot be distinguished from a pitch rotation of the agent. In this case the coupling matrix with constant distances becomes nearly singular and cannot be properly inverted.
In section 5.4 of the appendix it is shown by means of a coordinate transformation of two dimensional flow vectors (in tangent planes) into three dimensional ones (on the sphere) that Eq (20) is equivalent to Eq (4) with the weights given by Eq (6).

Properties of the coupling matrix.
The equivalence between the MFA and one iteration of the KvD algorithm is shown in three steps. The first step was the alternative derivation of the MFA described above. The following second step is a simplification of the four submatrices of the coupling matrix.
Three cases have to be considered when calculating the entries of the matrix: the scalar product between two translational templates, the scalar product between two rotational templates and the scalar product between a translational and a rotational template.
The scalar product between two translational templates leads to the following expression: The scalar product between two rotational templates results in Finally, the scalar product between a translational and a rotational template leads to Eq (21) informs us that the estimates of the three translation parameters t 1 , t 2 and t 3 are coupled unless the termẽ a Ád Áẽ b Ád D E in M tt is proportional to the identity matrix. The same holds for the coupling between the three rotation parameters r 1 , r 2 and r 3 described by M rr , see Eq (22). Similarly, the term md D E in M tr , Eq (23), must be zero for the translation and rotation estimates to be uncoupled.

2.3.2
The case of constant distances and a spherical field of view. Borst and Weber [10] showed that for viewing directions homogeneously covering the whole sphere, and for identical distances (μ i = constant for all i) the model neurons respond only to the components of the flow field they are tuned to. This result can be easily verified within the conceptual framework provided here by replacing the sums in Eqs (21), (22) and (23)  ÞÁd Wφ sin WdWdφ which can be regarded as the product between a first-order dipole function and the zeroth order spherical harmonic function (which is a constant). Due to the orthogonality of the spherical harmonics this integral is zero.

The relationship between the MFA and the KvD algorithm
In a final step we will show that the MFA with the coupling matrix is equivalent to one iteration of the KvD algorithm [1] for known distances. Hence, the two approaches do not represent principally different methods, but rather one method with two different ways of dealing with the depth structure. An additional way to take the depth structure into account is given in section 2.6 where an adaptive MFA is proposed. Dahmen et al. [12] have already demonstrated that a MFA can be derived from the KvD algorithm, if certain terms are small and can be neglected. We will show that these terms are identical to entries of the coupling matrix, which was derived in the previous chapter.
The equivalence of the MFA and one iteration step of the KvD algorithm. Not only in the MFA, but also in the KvD approach coupling terms can be identified. In the translation Eq (11), the second term on the right side comprises the rotation. The rotation Eq (12) in turn contains a translational term. These terms were called "apparent rotation" and "apparent translation" by Dahmen et al. [12], and they couple translation and rotation. The third terms in both equations couple different components of the translation or the rotation whenever Ár contain off-diagonal components. The following derivation will show that these coupling terms are equivalent to the terms of the coupling matrix in the MFA. The derivation starts with the MFA including the coupling matrix. From Eq (17) one obtains together with the equations for the coupling matrix Eqs (21), (22) and (23): On the left side of Eq (24), the vectorst andr are multiplied with the entries of the coupling matrix. On the right side the optical flow vectorsp are multiplied with the templates. After some algebraic simplifications, which are given in appendix 5.1, and a rearrangement of the terms one obtains the known equations for translation Eq (11) and rotation Eq (12) of the KvD algorithm:t Hence, the MFA and the KvD algorithm are identical for one iteration step.

The bias of the KvD algorithm
The equivalence of the KvD algorithm and the MFA for known distances follows directly from the Gauss-Markov theorem [13], which states that an ordinary least-squares estimator is the best unbiased estimator for an estimation problem which is linear and has uncorrelated errors with equal variances. Both methods start with such a least-squares approach. The MFA minimizes the quadratic error between the six true self-motion values and the estimated values as can be seen in Eq (5), whereas the KvD algorithm minimizes the quadratic error between the measured optic flow and the theoretical optic flow as can be seen in Eq (8). For known distances, the optic flow and the self-motion values are connected through a linear transformation. Thus, the two least-square approaches lead to the same self-motion estimator. This estimator is the unique optimal estimator as stated by the Gauss-Markov theorem. The situation is different if the distances are not known and have to be estimated by the estimator together with the self-motion values. Then the problem is no longer linear and the Gauss-Markov theorem does not hold. Nonetheless, it seems likely that a least-square approach like in Eq (8) leads to an optimal estimator in the sense that the error in the estimated selfmotion components approaches zero with increasing number of measured optical flow values. However, as will be shown in this section, the KvD algorithm in general does not converge to the true self-motion values for an infinite number of flow vectors. It is still an open question from which minimization principle an optimal estimator can be derived. The increasing number of flow vectors raises the problem that for each flow vector which gives two additional error-prone values one additional value has to be estimated: the inverse distance in the respective viewing direction. Although the number of measured values increases towards infinity, the ratio between the number of estimated and measured values does not decrease to zero. Hence, even for an infinite number of flow vectors the estimated inverse distances are still afflicted with errors. However, it should still be possible to correctly estimate the fixed number of selfmotion values for an infinite number of flow vectors. In section 2.5.2 a modified KvD algorithm will be derived. The modification is tested numerically under two conditions where the original KvD algorithm turns out to be biased (section 2.5.3).
2.5.1 The non-vanishing error term. The KvD algorithm is an unbiased estimator only under certain conditions. To show this, the propagation of the error in the flow vectorsp i over the iterations will be analyzed. We model the measured flow vectorsp i as the sum of the true vectorp 0 i and a random error vector Dp i ,p i ¼p 0 i þ Dp i . Similar to vectorp i , the vectorsp 0 i and Dp i have only two degrees of freedom. It will be assumed that the random vectors Dp i are unbiased, i.e. the expectation values for all directions i are zero, E Dp i ð Þ ¼ 0. Two special conditions will be considered in the following: 1. The viewing directionsd i are equally distributed over the sphere.

The random vectors Dp i are uncorrelated and their variances constant, independent of the
The Gauss-Markov theorem assumes condition (2) to be fulfilled. We will also consider deviations from this condition, because the KvD algorithm and its modified version that will be derived in section 2.5.2 behave differently then. Condition (2) is violated if, for example, the error of the optic flow measurement depends on the length of the measured optic flow vectors. In subsection 5.2 of the appendix it is shown that, in general, the translation estimated by the KvD algorithm contains errors that do not vanish even for an increasing number of flow vectors. There are two error terms which are additive to the real translationt 0 .
The index 1 of the brackets hi stands for the limit of an infinite number of flow vectors. The discrete direction vectorsd i can be exchanged by their continuous counterpartsd Wφ , and the sum over i can be replaced by an integral over the field of view. The iteration equation forr, Eq (12), does not contain terms that lead to a bias in the estimated motion values. However, the estimated rotation will be affected indirectly by errors in the estimated translation.
First the integrals containing the denominators D 1 ¼ 1 Àt 0 Ád 2 and D 2 ¼ are analyzed, disregarding the numerators. The integrals over 1 D 1 and 1 D 2 are zero only if the first condition of equally distributed flow vectors is fulfilled. To avoid the singularity ford ¼t 0 a small constant ε was added to the denominators.
If conditions (1) and (2) are fulfilled, the terms a and b are zero and the KvD algorithm is an unbiased estimator.
If condition (1) holds but condition (2) does not as, for instance, in the case of realistic EMDs or gradient-based detectors, we have to integrate over a direction dependent function resulting from the direction dependent flow errors Dp. Hence the terms a and b converge to finite values.
The error terms a and b do not play a role if they are proportional to the identity matrix, because of the rescaling of the translation vector, which ensures ktk ¼ 1. The matrices Most interestingly, if condition (2) is fulfilled (a pre-condition of the Gauss-Markov theorem) but condition (1) is not, the terms a and b converge to finite values. In this case the integrals over the denominators D 1 , D 2 and the integrals over the numerators E 1 , E 2 have finite values. This means that the ordinary least-squares approach from Eq (8) leads to a biased selfmotion estimator.

Modification of the KvD iteration equations.
To improve the KvD algorithm a modified version of the iteration Eq (11) will be derived. The flow Eq (1) can be transformed: By taking the average and solving fort we obtain where ξ ensures thatt is normalized. Compared with the original equation for the translation Eq (11) (1) is not. The translation error of the original KvD algorithm is significantly larger and increasingly deviates from that of the modified version with increasing number of flow vectors. Due to the coupling of the iteration equations, the error of the rotation, in case of the original KvD algorithm, is affected by the translation error and, thus, also deviates from that of the modified version.
The right part of the figure shows results where the standard deviation of 4p i is proportional to the length of the local flow vectorp i . This time, the viewing directions cover the whole sphere homogeneously and, thus, condition (1) of section 2.5.1 holds while condition (2) does not. Again, the original KvD algorithm is biased (see also section 2.5.1). However, the error of the rotation in the original KvD algorithm is not influenced by the translation error as a consequence of the spherically distributed viewing directions.
The error of the modified KvD algorithm is in both analyzed cases inversely proportional to the square root of the number of flow vectors (see black line in Fig 1). The modified KvD algorithm shows therefore an error behavior as is characteristic of an unbiased linear estimator.

An adaptive MFA
The approach of Franz and Krapp [8] is based on the assumption that the statistics of the depth structure of the environment is fixed as well as the preferred translation directions of the agent are known. For the simplest statistical model, which assumes that the distance variability, the noise in the flow measurements and the preferred translation of the agent are independent from the viewing direction, the dedicated covariance matrices are proportional to the identity matrix. Most important, one has to specify the depth structure of the environment by defining the average inverse distances hμi i or, as in [8], the average distances hDi i . Franz and Krapp [8] modeled the distances hDi i by where D 0 denotes a typical distance in the upper hemisphere, ε i the elevation angle of viewing direction i and b ¼ h D 0 < 1 the ratio of the average flight altitude h and D 0 . This model takes into account that the distances are usually smaller for viewing directions below the horizon.
It is obvious that the performance of self-motion estimation is poor when the fixed depth model Eq (39) is not a good description of the depth structure at the current position of the agent. Therefore, we propose an adaptive MFA, which allows the agent to adapt its depth model to the current environment. This is only possible, if we can assume that the depth structure properties of the environment do not change abruptly from one time step to the next. A time constant describes the intervals in which the depth model will be updated. Between updates the depth model remains fixed and has to be memorized by the agent. The depth information is obtained from the optic flow, as well as the self-motion values. Since the original KvD algorithm is biased in this case, we use our modified KvD version as the starting point for deriving the adaptive model. doi:10.1371/journal.pone.0128413.g001

Self-Motion Estimation with Adaptive Filters
As will be shown in the following, it is not necessary to represent the full depth information, because it is sufficient to memorize just eight distance dependent parameters. Three parameters are contained in the distance-dependent term of the rotation Eq (12) The first three orders of spherical harmonic functions, the zeroth, first and second order, comprise nine parameters, but only eight parameters can be determined from the optic flow. The zeroth order remains undefined, because one cannot distinguish between the situation, where all objects have half the distance to the agent, and the situation, where the agent translates with double speed: The optical flow will be the same. Hence, from optic flow fields only the direction of the translation can be identified.
Between updates of the depth model, rotations, in particular, impair its validity, but can be easily compensated by counter-rotating the depth model. This is achieved by multiplying the depth dependent coupling matrix with a rotation matrix obtained from the rotation parameters of the current self-motion estimate.
The adaptive model will be derived first for the general case with an arbitrary field of view. Then a more biologically plausible adaptive model with a spherical field of view will be presented. A spherical field of view facilitates an intuitive interpretation of the depth model. In this specific case, the nine depth-dependent parameters are exactly the coefficients of the first three orders of the spherical harmonics expansion, i.e the dipole and quadrupole moments of the depth structure μ ϑφ .
2.6.1 Motivation for an adaptive MFA. In Fig 2 the initialization phase of the adaptive model is shown. The agent flies inside a sphere in such a way, that the depth model is the same for every trajectory step (see Figs 2 and 3B and methods 4.4 for details). The error in the estimated self-motion parameters decreases exponentially. Hence, the error is corrected to a large extent in the first few iteration steps.
In natural environments one can detect subspaces where the distances to a moving agent do not change over a certain time [35]. Hence, one can assume that the overall depth model changes only slightly from one step to the next in a given environment and can expect good self-motion estimates even for only a single iteration step of the KvD algorithm based on the old depth model. Furthermore, the depth model can be used even for a longer time interval. In this case the depth model is not updated instantaneously after receiving new optic flow information, but less frequently after several optic flow processing steps. Because the self-motion parameters and the depth model are formulated in the body coordinate system of the agent, the depth model has to be rotated together with the agent.
2.6.2 Matched filters and depth-dependent coupling matrix for the adaptive MFA. The equations of the modified KvD algorithm form the basis of the adaptive MFA (a small constant ε is inserted in Eq (33), the equation of the inverse distance, to avoid the singularity in cased i ¼t):  r ¼p Âd As shown in section 2.4 these iteration equations can be decomposed into the product of the current optical flow and a standard template, on the one side, and a coupling matrix, on the other side (see following equation). The coupling matrix is the part of the adaptive model that contains the depth values μ i : M tr ¼ À1dÂ where the matrix1dÂ is defined by1 dÂṽ ¼d Âṽ; i.e. the cross product of the two vectorsd andṽ can be expressed as multiplication of matrix ½dÂ and vectorṽ. 2.6.3 The case of spherically distributed flow vectors. For the special case of spherically distributed flow vectors the depth-dependent coupling matrix (see subsection (2.3.1) is given explicitly. In the simplest case of an agent being in the center of a sphere, i.e. if all inverse distances have the same value, the depth-dependent coupling matrix is proportional to the identity matrix (see end of subsection (2.3.1). In the general case, i.e. when the inverse distances can have arbitrary values, the entries of the depth-dependent coupling matrix are the expansion coefficients of the first three orders of the real valued spherical harmonic functions.
Then the environmental depth structure μ i can be described by a real-valued function μ ϑφ parameterized by the azimuth angle φ and the elevation angle ϑ (in a spherical coordinate system). Such functions can be described by an expansion of real-valued spherical harmonic functions [36] (see appendix 5.3 and Fig 4). Lower orders of these functions contain less details than higher order functions.
Given a function μ(ϑ, φ) depending on azimuth angle φ and elevation angle ϑ the expansion coefficients a ln are R ln (ϑ, φ) represent, for example, a dipole function for specific l and n. The corresponding coefficient a ln provides information about how pronounced the specific dipole part is in the expanded function μ (ϑ, φ).
With the help of all coefficients a ln the expanded function μ (ϑ, φ) is given by the reverse transformation Due to the linearity of an integral expression and the orthogonality of the spherical harmonic functions [36],  doi:10.1371/journal.pone.0128413.g004

Self-Motion Estimation with Adaptive Filters
replaced by the coefficients of spherical harmonic functions: where a is the coefficient of the zeroth spherical harmonic function, which cannot be determined by the algorithm as explained earlier. The three b's are the three coefficients of the firstorder spherical harmonic functions, the dipole functions. The five c's are the five coefficients of the second-order spherical harmonic functions, the quadrupole functions. The only non-constant parameters in the depth-dependent matrix are the nine coefficients of the expansion by spherical harmonics. The nine parameters the agent has to memorize during flight have, in the case regarded here, a physical interpretation: They are the dipole and quadrupole parts of the depth distribution of the environment. The non-existence of higherorder coefficients in the depth-dependent matrix indicates that these orders contain no information for solving the self-motion problem. If the distances are constant, the matrix M is the identity matrix (except for a normalization constant) as mentioned before.
The computation of the nine coefficients might not seem to be biologically plausible at first sight. However, the computation of one coefficient a ln of the expansion corresponds to a weighted wide-field integration and is reminiscent of the function of LPTCs in flies [16,17,25,26]. One could imagine a neuron for each of the nine parameters a ln , which represents a specific global property of the depth structure.
With regard to the computational effort a full inversion of a 6 × 6 matrix is not required. The submatrix M tr is zero in the spherical case, hence only an inversion of the submatrix M tt is required. If the quadrupoles in this submatrix are sufficiently small, the inverse matrix can be linearly approximated by a Neumann series [37], (I−A) −1 = I+A+A 2 +A 3 +. . ..
2.6.4 Test of the adaptive MFA. In this section we compare quantitatively the adaptive and the original MFA. We present a simulation in a very simple environment where the agent moves inside a sphere (see Fig 3 and methods 4.4). Nothing is known in advance about the flight directions and whereabouts of the agent inside the sphere. Hence, the covariance matrices of the original MFA are set proportional to the unit matrix. The chosen trajectory in this setting is a sinusoidal curve. On this trajectory the current depth distribution differs essentially from that in the center of the sphere (Fig 3A).
The amount of translation and rotation in each step varies along the trajectory. The steps are chosen so that the maximum rotation angle is about 8 degrees per step ensuring that the approximation of the KvD algorithm, which is valid for small rotation angles, still holds. The maximum rotation angles between the processed images is in the same range as the maximum rotation angles during saccades of flying insects [38,39].
On the whole, the adaptive MFA performs better than the original non-adaptive MFA ( Fig  5). If an Gaussian error is added to the optical flow with a standard deviation of ten percent of the averaged overall flow value (Fig 5), one could expect that a translational or rotational flow component, which is overlayed hundredfold by the other flow component disappears totally in this flow error, because the error is tenfold higher as the flow component in this situation. But the estimators use hundreds or a few thousands of flow vectors to estimate self-motion. Due to the large number of flow vectors (insects usually have between a few hundred (e.g. fruit fly) and a few thousand ommatidia (bee, dragonfly) per eye), the self-motion can be still estimated in this case within a useful error range. The results are shown for about 5000 flow vectors ( Fig  5). The error increases for both the translational and the rotational self-motion estimate to a value of about 10 degrees. This error is additive to the error described in the upper panels of Albeit the simplicity of the simulation it shows some basic features of the compared algorithms. The simulation does not generate any outliers due to moving objects or depth discontinuities. A small number of outliers will not affect the MFA as a consequence of the linear summation over thousands of optic flow vectors. More complex simulations in virtual environments with rendered images and EMDs is left to future work (Strübbe et al., in prep.).

Discussion
The aim of this study is to develop an adaptive matched filter approach to self-motion estimation which could be in principle the underlying concept of self-motion estimation in flying insects. As a novel characteristic, this approach assumes an adaptation to the depth structure in the insect visual motion pathway, an assumption that is supported by recent experimental evidence [31,40,41]. Our approach starts from a theoretical point of view by analysing and unifying the non-adaptive matched filter approach (MFA) to self-motion estimation together with the Koenderink van Doorn (KVD) algorithm which incorporates an estimation of the depth values.
To take advantage of both algorithms, some mathematical problems had to be solved. First, it was shown that the two algorithms are equivalent in case the distances to objects in the Self-Motion Estimation with Adaptive Filters environment are assumed to be known. Secondly, a bias in the KvD algorithm was removed by a small correction of the iteration equations. And last but not least, an analysis of the specific case of a spherical field of view, reminiscent of that of flying insects, shows that the depth structure can be represented by only eight parameters without losing relevant information for selfmotion estimation and that these eight parameters are the dipole and quadrupole moments of the spherical harmonics. Technical and biological systems have different origins and often operate under different conditions. Biological systems arise through evolutionary adaptation. They usually have to operate in a great variety of environments. Hence, the neural computations underlying the animal's behavior need to be particularly robust. In addition, the animal has restrictions with respect to its computational hardware. Neuronal circuits can perform linear transformations in parallel by a dendritic weighted summation of the inputs of a neuron and non-linear operations through the non-linear response behavior of a neuron to its overall input. Nonetheless, a non-linear operation, such as computing the inverse of a matrix, with changing entries, is not easy to implement by neuronal hardware. The bio-inspired computational model analyzed here is the MFA of self-motion estimation. It is a linear model for a fixed depth assumption, derived under the side condition of maximal robustness against errors in the measured optical flow field (see Eq (5) from Franz et al. [11]).
The KvD algorithm of self-motion estimation which is compared with the MFA model was derived analytically in a technical framework on the basis of a minimization principle. The resulting iteration equations represent a gradient descent where the current self-motion parameters are used to determine a better depth model and the new depth model is used to determine a better estimate of self-motion in the next iteration step. If one considers only one iteration step, where the depth model is seen as fixed, the self-motion estimation is linear. It is not only linear, but also equivalent to the biologically inspired MFA which uses a fixed depth distribution.
The equivalence of both models becomes evident within the framework of linear estimator theory. There exists an unique optimal estimator for a linear estimation problem with errorprone inputs. The Gauss-Markov theorem describes this estimator. Both compared methods represent this optimal solution.
However, some differences exists. In the MFA Franz et al. [11] weighted the filters by matrices that represent additional assumptions about the situation under which the self-motion is estimated. If these assumptions are correct, the weighting of the filters improves self-motion estimation; however, if these assumptions are incorrect in the current situation, the estimator gets worse. Hence, the additional matrices make the estimator more specific. These matrices can also be implemented in the KvD algorithm by a modification of the minimization principle. When, for example, it is known, that the optical flow can be measured more accurately below the horizon, because the objects are generally closer there, this knowledge can be taken into account by introducing weights in the initial equation. We argue that it is not always useful to take knowledge about the preferred self-motion directions into account. Even when the moving agent solely translates in the forward direction a disturbance can lead to a passive translation also in other directions.
From a mathematical point of view the bias of the KvD algorithm is remarkable. When the depth distribution of the environment has to be determined together with the self-motion parameters, the estimation problem is no longer linear. The standard procedure for estimating parameters from inputs, which are disturbed by Gaussian errors, is the minimization of the mean squared error. It might be counter-intuitive that the true self-motion values are not even a local minimum. The standard approach fails, because the standard condition assumes that an increasing number of measured values, here the flow vectors, are accompanied by a fixed number of estimation parameters. However, for the non-linear estimation problem the number of distance values increases together with the number of flow vectors. Only the number of selfmotion parameters remains constant. With every additional flow vector additional information about the self-motion parameters is obtained, because one gets two additional independent values from the flow vector, but only one additional parameter has to be estimated (the distance corresponding to the flow vector). However, the standard approach does not use this additional information in an optimal way.
Here we derived a modified version of the KvD algorithm that is not derived from a minimization principle. Hence it is not clear whether it leads to the best estimate for a given finite number of flow vectors. Rather, the numerical simulation indicates that the modified version has the desired property that the algorithm converges to the real self-motion values for an infinite number of flow vectors. It is left to further mathematical work to analyze optimal criteria for the non-linear estimation problem in case of a finite number of flow vectors.
Based on this modified version of the KvD algorithm an adaptive MFA was derived. It was shown that it is a critical issue to correctly determine the dipole components. If a small rotation is superimposed by a large translation the non-adaptive MFA cannot provide useful rotation estimates, whereas the adaptive MFA is accurate up to a few degrees. A situation where a relatively large translation encounters a relatively small rotation is given in the inter-saccadic phases of insect flight [38,39,42]. In these phases the insect tries to avoid any rotation. If the insect stabilizes its flight with the help of the visual system, the non-adaptive MFA cannot be the underlying concept to detect small rotations in these phases in environments it is not tuned to. To estimate rotations, which are superimposed by a large translation, one has to determine the current dipole components, as is done by the adaptive MFA.
The adaptive MFA was inspired by the finding that for a spherical field of view the depth structure of the environment can be represented by only eight parameters without losing relevant information for self-motion estimation and by the fact that the visual system of insects has an almost spherical field of view. The spherical field of view is also a desirable property of technical systems which are designed to estimate their self-motion on the basis of optical flow fields. Such systems can be realized by panoramic cameras [43,44].
Adaptation to the depth structure of the environment means that the adaptation takes place on another time scale than the image processing itself. Hence, some information about the depth structure has to be memorized by the system. The result that exact self-motion can only be estimated for a spherical field of view, if eight parameters about the depth structure of the environment are known, is therefore in accordance with the limited computational resources of insects.
Motion adaptation was analyzed in the insect visual pathway and found to depend on the overall velocity in the visual field [27,31,32,40,41,45]. Since, at least during translational motion, the overall retinal velocity depends on the depth distribution of the environment. The experimentally characterized processes of motion adaptation may well play a role in an adaptive mechanism of self-motion estimation as proposed in the present study. Here we give a short analysis from a theoretical point of view which components are needed for the adaptive MFA. Minimalistically, one needs eight model neurons for the eight depth parameters. The weighted summation over the inputs of one of these model neurons corresponds to one of the eight integrals over the depth distribution, where the spherical harmonic functions play the role of the weighting parameters. Examples of neurons performing such an integration are the LPTC neurons of flies, the neuronal candidates for the six model neurons, which represent the matched filters for the six self-motion components. Given the properties of LPTCs [16,17,[24][25][26], it is likely that one hypothetical model neuron for depth representation does not cover the whole sphere. Due to the linearity of self-motion estimation LPTCs can be combined to represent one self-motion component. On this basis, it might be possible that one LPTC codes information for both translation and rotation when the corresponding flow fields resemble each other within the receptive field of the neuron. Hence, the hypothetical depth neurons could be realized by a network of neurons with each neuron receiving input from only part of the visual field.
The hypothetical neurons representing the depth structure need some pre-and post-processing. In the adaptive MFA, only the pre-processing contains non-linear operations, namely the transformation of the optical flow into local depth values that are integrated afterwards by the hypothetical model neurons. One can assume that the determination of the depth values is simplified during nearly pure translation as is characteristic of the insect saccadic flight strategy [19]. In these phases the depth structure can be determined more easily, because the optic flow is not superimposed by a rotational component [46,47].
The post-processing concerns the determination of the depth-dependent matrix which corrects the outputs of the six model neurons, corresponding to the motion sensitive LPTC cells. From a mathematical point of view, two subsequent linear transformations can be combined in a single linear transformation. In the adaptive MFA we have two subsequent linear transformations: the fixed linear transformation by the six model neurons that receive direct optical flow input and the adaptive linear transformation by the depth-dependent matrix, the entries of which are the responses of the eight depth neurons. There are two options where the adaptation could take place: The two linear transformations could be merged into one linear transformation, which means that the adaptation takes place at an early stage of optic flow processing. Alternatively, one could assume that the two linear transformations are spatially separated, and the depth-dependent matrix is realized by an adapting linear circuitry, which wires the early stage neurons.
The linear transformation given by the depth dependent matrix can be obtained without a matrix inversion by applying an appropriate linearization of the inverse depth dependent matrix (see section 2.6.3). With this simplification and the above simplification of depth capturing the adaptive MFA can be realized by relatively simple neuronal circuits.

Numerical test and simulation
We used a numerical test (Fig 1) and a simulation (Fig 5) to show the performance of the considered self-motion estimators based on optical flow fields. The optical flow fields used here are computed in all cases from the flow Eq (1) where the distances, the viewing directions and the self-motion parameters are given. In some cases a Gaussian error was added to the optical flow values. The task of the self-motion estimators is to determine the self-motion parameters from these flow fields.
Whereas the distances in the numerical test are obtained by a random process, the distances in the simulation are defined by the environment. In the numerical test the individual estimates are independent of each other. In the simulation a trajectory through the environment is constructed to enable an agent to use the depth model of the environment for a series of estimates at subsequent trajectory points.
For programming and testing the algorithm the programming language Matlab was used.

Construction of a spherical field of view
For the simulation and the numerical test a spherical visual field of view is needed. It is not a trivial task to arrange a number of viewing directions on a sphere in a way that the density is equally distributed. Here we used an iterative solution. The iteration starts with a Platonic solid, an octahedron. An octahedron has eight faces: Eight equilateral triangles with the same side lengths, which surround a symmetric solid.
In each iteration step every triangle is replaced by four new triangles with one triangle placed in the middle of the old triangle and the other three are placed in corners. The mid-points of the sides of the old triangle coincide with the corners of the new triangles. These midpoints are projected radially on the surface of the sphere which surround the object so that all corners of the new triangles lie on this sphere.
In each iteration step the sphere is covered with nearly identical triangles. After the last iteration step the mid-points of each triangle give a viewing directions. With an arbitrary number of iterations n the number of viewing directions are 8 Á 4 n .

Numerical test of the bias of the KvD algorithm
In the numerical test of the bias of the KvD algorithm two configurations are tested, one with a spherical field of view and one with a non-spherical field. The non-spherical field of view is realized by the same procedure as in the spherical case (see methods 4.2), but two starting triangles are left out. The omitted triangles are opposite triangles both placed in the upper part of the sphere. This configuration has some obvious symmetry, but avoids the case that each omitted viewing direction has an omitted counterpart on the other side of the sphere.
The number of viewing directions is increased by increasing the number of iterations starting from an octahedron. Each step increases the number of viewing directions by four. For each number of viewing directions 40 self-motion estimation trials are tested.

Simulation of the adaptive MFA
The environment used for the simulation is a unit sphere in which the agent flies on a sinusoidal trajectory. The effect of adaptation should be larger if, for example, the agent flies outdoors and then enters a tunnel which is a common experimental set-up for navigation studies in flying insects [48][49][50][51]. However, the test of the adaptive MFA in this kind of environment needs a 3-d engine for rendering images, on which the flow vectors are estimated with motion detectors like the Reichardt detector or the Lucas-Kanade detector. This issue will be addressed in another study (Strübbe et al., in prep.).
For the trajectory a sinusoidal curve was chosen (see Fig 3 right). Although it is typical that flying insects use a saccadic flight strategy to separate translation and rotation, the self-motion estimators in this study were analyzed under the condition of a combined translation and rotation. In each run the agent flies along the trajectory and the angle error between the true selfmotion axes and the estimated axes of the tested estimators are shown in Fig 5.

Derivation of the equivalence of the MFA and KvD algorithm
From Eq (17) one obtains together with the equations for the coupling matrix Eqs (21), (22) and (23) and the definitions of the templates Eqs (13) and (14): Self-Motion Estimation with Adaptive Filters Multiplyingt andr with the matrices andp with the templates leads to the two equations: The term To show that the KvD algorithm does not converge to the real values for an infinite number of flow vectors, we start by assuming the opposite. If the KvD algorithm had a minimum at the true valuest 0 ,r 0 and m 0 i we could insert these values in the iteration Eqs (10), (11) and (12) and would get the same values back. Substituting the true values for translation and rotation in the equation for the nearness Eq (10), we get the nearness error 4μ i which directly depends on the flow vector errors: The estimated translation after an infinite number of iterationstð1Þ ¼ lim n ! 1t ðnÞ shows the following equation. Please note that we consider two limits: the limit of infinite iteration steps, lim n ! 1 , and the limit of an infinite number of flow vectors, indicated by an infinity symbol as index of the brackets which stands for the mean over all viewing directions.
Iftð1Þ is a stable minimum 4t must vanish: For the analysis of the conditions under which the term 4t vanishes, see the text (subsection 2.5.1).

Expression ofd andd d by spherical harmonics
The vectord and the dyadic productd d are expressed by linear combinations of real-valued spherical harmonics, which are itself linear combinations of complex spherical harmonicsY lm , The real-valued spherical harmonics from zero order to second-order are: Zeroth-Order First-order  1, 2, 3). This can be seen when the components ofd are written in a spherical coordinate system: The components d 1 , d 2 and d 3 are equal in this arrangement to the first-order functions f 1 , f 2 and f 3 except of a normalization factor.

Expression ofd d through spherical harmonics.
First the dyadic productd d is formulated in a spherical coordinate system The following analysis shows that the off-diagonal elements of this matrix can each be r g 0

The weight matrix of the original MFA
To see that Eq (20) corresponds to Eq (4) as derived by Franz et al. in 2004 [11] for the original MFA, a change of the basis of the vector space is applied. In [11] the coordinates of each flow vectorp i were given with respect to different basis vectorsũ i ,ṽ i , spanning the tangential plane on the sphere for viewing directiond i . In the following we derive the transformation matrix to transform the coordinates of a vector given in the standard Euclidean basis vectorsẽ 1 ,ẽ 2 andẽ 3 into the basis defined by ðũ i ;ṽ i ;d i Þ i¼1;2;...;N . With N optical flow vectors we define a 3 × N dimensional vector space, which is the N-fold Cartesian product of the three basis vectorsẽ 1 ,ẽ 2 andẽ 3 . All flow vectorsp i are represented now by a single stacked column vectorp, which has the dimension of 3 × N. In the same way the templates are written as where the index A stands for one of the six standard templates. The six templatesT A are combined to matrixT with six columns and 3 × n rows: T ¼T The responses of the six model neuronsã in this notation arẽ a ¼T T Áp; ð82Þ whereT T stands for the transpose ofT . Together with the coupling matrix one obtains for the self-motion components (see Eq (20)): y est ¼M À1 Áã ¼M À1 ÁT T Áp: ð83Þ The complete vector basis of the above templates and flow vectors is described by two indices j = 1,2,3 and i = 1, 2, . . ., N and has the form B ¼ẽ 1;1 ;ẽ 2;1 ;ẽ 3;1 ;ẽ 1;2 ;ẽ 2;2 ;ẽ 3;2 ; . . . ;ẽ 1;N ;ẽ 2;N ;ẽ 3;N À Á ; ð84Þ where theẽ 1;i ,ẽ 2;i andẽ 3;i are the same for each i. The basis vectors in the notation of Franz et al. [11] are B 0 ¼ũ 1 ;ṽ 1 ;d 1 ;ũ 2 ;ṽ 2 ;d 2 ; . . . ;ũ n ;ṽ n ;d n ; ð85Þ where the viewing directionsd i supplemented the local vector spaces to a three dimensional space.
The new basis vectors B 0 can be expressed by the old B, The matrixV transforms the coordinates of a vector given in basis B, Eq (84), into basis B 0 , Eq (85). Eq (83) becomes The termV Áp is equivalent tox in Eq (4) from the original MFA of Franz et al. [11]. The term x combines all optic flow components of the local two dimensional vector spaces as spanned by the two local vectorsũ andṽ. The third component is zero becausep andd are orthogonal.

! T
the first row ofV ÁT is scrutinized.

Acknowledgments
We are grateful to Jens Peter Lindemann for programming support. This study was supported by the Deutsche Forschungsgemeinschaft (DFG).