Optical Flow on Evolving Sphere-Like Surfaces

In this work we consider optical flow on evolving Riemannian 2-manifolds which can be parametrised from the 2-sphere. Our main motivation is to estimate cell motion in time-lapse volumetric microscopy images depicting fluorescently labelled cells of a live zebrafish embryo. We exploit the fact that the recorded cells float on the surface of the embryo and allow for the extraction of an image sequence together with a sphere-like surface. We solve the resulting variational problem by means of a Galerkin method based on vector spherical harmonics and present numerical results computed from the aforementioned microscopy data.

1. Introduction. Motion estimation is a fundamental problem in image analysis and computer vision. An important task within is optical flow computation. It is concerned with the inference of a vector field which describes the displacements of brightness patterns, such as moving objects, in a sequence of images. Ever since the seminal work of Horn and Schunck [18] a variety of reliable and efficient methods have been proposed and successfully applied in a wide number of fields.
Primarily, optical flow is computed in the plane. However, it is readily generalised to non-Euclidean settings allowing, for instance, for cell motion analysis in timelapse microscopy data. It has been only recently that high-resolution observations of biological model organisms such as the zebrafish became possible. Despite its importance for tissue and organ formation, little is known about cell migration and proliferation patterns during the zebrafish's early embryonic development [1,34]. Fluorescence microscopy nowadays allows to record time-lapse images on the scale of single cells, see e.g. [20,28,34]. In view of increasing spatial as well as temporal resolutions producing vast amounts of data, manual analysis through visual inspection by humans is impracticable. Automated cell motion estimation therefore is key to large-scale analysis of such data. Optical flow computation delivers necessary quantitative methods and leads to insights into the underlying cellular mechanisms and the dynamic behaviour of cells. See, for example, [2,29,33,34] and the references therein.
The primary biological motivation for this work is the desire to analyse cell motion in a living zebrafish during early embryogenesis. The data at hand depict endodermal cells expressing a green fluorescent protein. By virtue of laser-scanning microscopy, (volumetric time-lapse) 4D images of these labelled cells can be recorded without capturing the background. It is known that endodermal cells float on a socalled monolayer during early embryonic development meaning that they do not stack on top of each other [39]. Figure 1 depicts two frames of the captured sequence, containing only the upper hemisphere of the animal embryo. Observe the salient formation of the cells and the noise present in the images. More precisely, one can see the nuclei of cells forming a round surface in a single layer. For more details on the microscopy data we refer to Sec. 5.1.
Recent efforts which aim at efficient motion estimation in such highly-sparse microscopy data regard this layer of cells as a two-dimensional surface and consider only the restriction-or a suitable projection-of the volumetric data to this prescribed geometry, see e.g. [23,24,34]. As a result, the spatial dimension of the data reduces by one, while simultaneously preserving the essential information.
In Schmid et al. [34], the geometry of choice is the static round sphere. Conceptually, we build upon this approach, however, address a concern raised in this article. In particular, in this article they acknowledge the need for more accurate shape representations, as the spherical model can suffer from inaccuracies. See, for instance, Fig. 1, the presumed surface in Fig. 2, and also the discussion in [34,Suppl. Inf.].
During early zebrafish development the initially spherical yolk undergoes deformation as endodermal cells converge towards the dorsal midline and the embryonic axis forms, see [21,Figs. 11 and 15] and Fig. 1. The major goal of this article is to address the following points. First, we allow for radial deviations from the perfect sphere and thereby effectively reduce errors from projections. Second, we model the layer of cells as an evolving surface, meaning that this surface reflects changes in the geometry as the embryo develops over time. Third, we estimate cell motion directly on the moving surface in contrast to computing it in-possibly distorting-map projections.
In what follows, we model the embryo of a zebrafish during its early development as a sphere-like surface. It is topologically diffeomorphic to the 2-sphere S 2 and most commonly defined as the set of points The functionρ : S 2 → (0, ∞) can be thought of as a radial deformation of S 2 and will have a dependence on time in the present paper. As a consequence, changes in the embryo's geometry are attributed accordingly, albeit valid only during early stages of its development as cells tend to cluster subsequently. The main intention of this work is to conceive cell motion only on this moving 2-dimensional manifold. Figure 2 depicts two frames of the surface together with images obtained by projecting the volumetric microscopy data in Fig. 1 onto the sphere-like geometry. The sequence contains a total number of 151 frames recorded at time intervals of 120 s. Fluorescence response is indicated by blue colour and is proportional to the observed intensity. The embryonic axis of the animal forms around the clearly visible dent. All dimensions are in micrometer (µm).
In this work we model the data as a time-dependent non-negative functionf . Its value directly corresponds to the fluorescence response of the observed cells. For a fixed time instant t ∈ [0, T ], the domain off is presumed to be a closed surface M t ⊂ R 3 . We assume that this surface can be parametrised by a smooth radial map from the 2-sphere. The temporal evolution of the dataf can then be tracked by solving an optical flow problem on this moving surface or, more conveniently, an equivalent problem on the round sphere.
Traditionally, the starting point for optical flow is the assumption of constant brightness: a point moving along a trajectory does not change its intensity over time. On a moving domain M = {M t } t one equivalently seeks, for every time t ∈ [0, T ], a tangent vector fieldv that solves a generalised optical flow equation at every point x ∈ M, wheref is the image sequence living on M. Here, for a fixed time t, ∇ M denotes the (spatial) surface gradient, dot the standard inner product, and dV tf an appropriate temporal derivative. The optical flow problem is ill-posed meaning that equation (1) is not uniquely solvable. A common approach to deal with non-uniqueness is Tikhonov regularisation, which consists of computing a minimiser of The first term of the sum is usually the squared L 2 norm of the left-hand side of (1) and, in the present article, the second term will be an H 1 Sobolev norm.
1.1. Contributions. The primary concern of this article is optical flow computation on evolving 2-dimensional Riemannian manifolds which can be parametrised from the sphere. Motivated by the aforementioned zebrafish microscopy data we consider closed surfaces for which the mapping is a diffeomorphism between the 2-sphere and M t for every time t ∈ [0, T ]. As a prototypical example we restrict ourselves to radially parametrised surfaces as they suit quite naturally to the given data. The contributions of this work are as follows. First, we give a variational formulation of optical flow on 2-dimensional closed Riemannian manifolds. We assume a dependence on time and speak of evolving surfaces. The main idea is to solve the problem by a Galerkin method in a finite-dimensional subspace of an appropriate (vectorial) Sobolev space. We take advantage of the fact that tangential vector spherical harmonics form a complete orthonormal system for L 2 (S 2 , T S 2 ). The sought vector field is thus uniquely determined when expanded in terms of the pushforward-by means of the differential of (2)-of these functions. From that we arrive at a minimisation problem over R n , where n is the dimension of the finitedimensional space, and state the optimality conditions. They can be written purely in terms of spherical quantities and solved on the 2-sphere. To this end, we use a standard polyhedral approximation and locally interpolate spherical functions by piecewise quadratic polynomials. For numerical integration we employ appropriate quadrature rules on the approximated sphere.
Second, to obtain the smooth sphere-like surface, which is described by the map (2), from the observed microscopy data, we formulate another variational problem on the sphere. The problem is essentially surface interpolation with H s Sobolev seminorm regularisation. Approximate cell centres serve as sample points of the surface. In particular, our microscopy data are supported only on the upper hemisphere, see Figs. 1 and 2. Scalar spherical harmonics are the appropriate choice for the numerical solution of the surface fitting problem, as they provide great flexibility with respect to the chosen space H s .
Finally, we present numerical experiments on the basis of the mentioned cell microscopy data of a live zebrafish. To this end we compute an approximation of the sphere-shaped embryo and obtain a sequence of images living on this moving surface. Eventually, we solve for the optical flow and present the results in a visually adequate manner.
1.2. Related Work. The first variational formulation of optical flow is commonly attributed to Horn and Schunck [18]. They attempted to compute a displacement field in R 2 by minimising a Tikhonov-regularised energy functional. It favours spatially regular vector fields by penalising the squared H 1 Sobolev seminorm. For introductory material on the subject we refer to [4,5] and to [40] for a survey on various optical flow functionals. Well-posedness of the aforementioned energy was first shown by Schnörr [35]. Moreover, there the problem was extended to irregular planar domains and solved by means of finite elements.
Weickert and Schnörr [42] considered a spatio-temporal model by extending the domain to R 2 × [0, T ]. It additionally favours temporal regularity of the solution by including first derivatives with respect to time. Such models are of particular interest whenever trajectories are to be computed from the optical flow field. A unifying framework including several spatial as well as temporal regularisers was proposed in [41]. For the purpose of evaluation and flow field visualisation a framework was created by Baker et al. [6].
Recently, generalisations to non-Euclidean domains have gained increasing attention. In [19] and [37] optical flow was considered in a spherical setting. Lefèvre and Baillet [27] adapted the Horn-Schunck functional to surfaces embedded in R 3 . Following Schnörr [35], they proved well-posedness of their formulation and employed a finite element method for solving the discrete problem on a triangle mesh. With an application to cell motion analysis, Kirisits et al. [22,24] recently considered optical flow on evolving surfaces with boundary. They generalised the spatio-temporal model in [42] to a non-Euclidean and dynamic setting. Eventually, the problem was tackled numerically by solving the corresponding Euler-Lagrange equations in the coordinate domain. Similarly, Bauer et al. [7] studied optical flow on time-varying domains, with and without spatial boundary. They proposed a treatment on surfaces parametrised by product manifolds, constructed an appropriate Riemannian metric, and proved well-posedness of their formulation.
In Kirisits et al. [23], the authors considered various decomposition models for optical flow on the 2-sphere. The proposed functionals were solved by means of projection to a finite-dimensional space spanned by vector spherical harmonics. Concerning projection methods, Schuster and Weickert [36] solved the optical flow problem in R 2 solely based on regularisation by discretisation.
Regarding sphere-like surfaces and spherical harmonics expansion of closed surfaces we refer to [32] and the references therein.
Finally, let us mention [2,29,33,34], where optical flow was employed for the analysis of cell motion in microscopy data. In particular, in Schmid et al. [34] the embryo of a zebrafish was modelled as a round sphere and motion of endodermal cells computed in map projections.
The remainder of this article is structured as follows. In Sec. 2, we formally introduce evolving sphere-like surfaces, recall the definition of vectorial Sobolev spaces on manifolds, and discuss both scalar and vector spherical harmonics on the 2-sphere. Section 3 is dedicated to optical flow on evolving surfaces and our variational formulation. In Sec. 4 we discuss the numerical solution. In particular, we propose to solve the resulting energy in a finite-dimensional subspace and rewrite the optimality conditions to be defined solely on the 2-sphere. Moreover, we show how to fit a sphere-like surface to the labelled cells in the microscopy data. Finally, in Sec. 5, we solve for the optical flow field and visualise the results. The appendix contains deferred material.

Sphere-Like Surfaces. Let
We denote by a smooth (local) parametrisation of S 2 mapping coordinates ξ = (ξ 1 , ξ 2 ) ∈ Ω to points x = (x 1 , x 2 , x 3 ) ∈ S 2 . Furthermore, let I := [0, T ] ⊂ R denote a time interval and let M = {M t } t∈I be a family of closed smooth 2-manifolds M t ⊂ R 3 . Each M t , t ∈ I, is assumed to be regular and oriented by the outward unit normal fieldN(t, x) ∈ R 3 , x ∈ M t . We assume that M (locally) admits a smooth parametrisation of the form whereρ : I × S 2 → (0, ∞), and call M an evolving sphere-like surface. We denote byf : M → R a smooth function on the moving surface. Its coordinate representation f : I × Ω → R and its corresponding spherical representatioñ f : I × S 2 → R are given by As a notational convention we indicate functions living on S 2 with a tilde and functions on M with a hat, respectively. Their corresponding coordinate version is treated without special indication. For convenience, we define a smooth (spatial) extension off to which is constant along radial lines. The extensionf (t, ·) is chosen such that it coincides withf (t, ·) on S 2 and withf (t, ·) on M t , respectively. For functions such asρ(t, ·) in (4) which are only introduced on S 2 , the extension is well-defined by considering M t = S 2 , that is,ρ(t, ·) ≡ 1 for all t ∈ I. Note that, whilef (t, ·) is Figure 3. Schematic illustration of a cut through the surfaces S 2 and M t intersecting the origin. In addition, we show a radial line along which the extensionf (t, ·) is constant. The surface normals are shown in grey.
constant in the direction of the surface normal of S 2 , it is in general not constant in the direction of the surface normal of M t . We point at Fig. 3 illustrating the setting.
Similarly, for a vector-valued functionv : M → R 3 the extension to R 3 \ {0} is defined component-wise and for all times t ∈ I, and is denoted byv. As a notational convention, boldface letters are used to denote vector fields. Moreover, we distinguish between lower and upper case boldface letters. The former identify tangent vector fields and their extensions to R 3 \ {0} whereas the latter indicate general vector fields in R 3 , with the exception of the parametrisations x and y. At this point let us refer to Tab. 1 at the end of this section containing a summary of the notation.
For a differentiable function f : I × Ω → R, we write ∂ i f as an abbreviation for the partial derivative of f with respect to ξ i . That is, The tangent plane at a point y(t, ξ) ∈ M t is denoted by T y(t,ξ) M t and the tangent bundle by T M t = {y(t, ξ)} × T y(t,ξ) M t : ξ ∈ Ω . The orthogonal projector onto the tangent plane T x M t at x ∈ M t , t ∈ I, is given by In particular if M t = S 2 , that isρ in (4) is identically one for all t ∈ I, the outward unit normal and the orthogonal projector are given byÑ and P S 2 , respectively.
In what follows, we define spatial differential operators. As they are identical to those on static surfaces we consider time t ∈ I arbitrary but fixed. Then, the surface gradient off , as given in (5), is defined by where ∇ R 3 denotes the usual gradient of the embedding space andf is the extension from (6). Let us stress that it is independent of the chosen extension, see e.g. [13, p. 389]. We emphasise that, in particular, if M t = S 2 for all t ∈ I it follows that

OPTICAL FLOW ON EVOLVING SPHERE-LIKE SURFACES
The last term of the sum on the right hand side is the normal derivative off , which according to the definition of the extension in (6) vanishes. Thus, For convenience let us observe that, by taking ∂ i f in (5), we arrive at due to the chain rule and the projection onto the tangent plane T x(ξ) S 2 .
Analogously to the surface gradient we define the spherical Laplace-Beltrami of f : where y is the parametrisation defined in (4), forms a basis of the tangent space T y(t,ξ) M t at y(t, ξ). Its elements form the gradient matrix Dy, which is derived as follows.
Letρ be the extension ofρ : I × S 2 → (0, ∞) according to (6). Then, y from (4) can be rewritten as By the chain rule, Using (8) and the fact thatρ equalsρ on S 2 gives where we have omitted the arguments (t, ξ) and (ξ) for better readability. Whenever convenient and no confusion will arise we will continue to do so. By applying (9) backwards and the fact thatρ(t, x(ξ)) = ρ(ξ) we have shown where Dx = (∂ 1 x, ∂ 2 x) is the gradient matrix associated with x.
Provided that x is a regular parametrisation it is straightforward to show that also y is regular. By computing Lagrange's identity one has ∂ 1 y × ∂ 2 y = 0 proving that the columns of Dy are linearly independent or, in other words, that y indeed is regular. We highlight that (19) below is in fact Lagrange's identity.
As a consequence, we can uniquely represent a tangent vectorv ∈ T y(t,ξ) [26,Prop. 3.15]. We call v i the components ofv.
In the sequel we will use Einstein summation convention. We sum over every index letter that appears exactly twice in an expression, once as a sub-and once as a superscript. For instance, we writev = v i ∂ i y for the sake of brevity.
We underline that the coordinate basis (11) is not orthogonal in general. We will, however, require an orthonormal frame {ê 1 (t, ξ),ê 2 (t, ξ)} of the tangent space T y(t,ξ) M t from this section on. In the coordinate basis it readŝ where α j i : I × Ω → R, i, j = {1, 2}, are functions obtained from the Gram-Schmidt process.
Combining (5) and (9) with the expressions derived for Dx and Dy we can conveniently state that Let us derive the following useful generalisation of (9). For a tangent vector where the third equality follows from the first equation in (14).
As soon as we have established the relation betweenv andṽ it will conveniently allow us to switch between (15) and (16). Moreover, the coordinate representation of the surface gradient (8) is derived as follows. Let us start out with the first equation in (14). By writing ∇ S 2f in the coordinate basis, that is ∇ S 2f = Dxv for some v, we obtain from (14) Furthermore, let t ∈ I be fixed and letf (t, ·) : where Jy is the Jacobian of y. According to Theorem 3 in [11, p. 88], it is given by and by using (12) yields Note that x·x = 1 and thus, terms of the form ∂ i x·x vanish. By the differentiability of x, ones has Therefore, ∂ i x · x = 0, meaning that tangential and normal vectors are orthogonal. We emphasise that Dy Dy is commonly referred to as Riemannian metric. It is positive definite and thus, (Jy(t, ξ)) 2 > 0 for all (t, ξ) ∈ I × Ω.
The parametrisations x and y defined in (3) and (4), respectively, suggest the straightforward construction of a smooth mapφ(t, ·) : S 2 → M t . It is given by the composition (y • x −1 )(t, ·), that is

is a linear map and is given by
It follows from a direct calculation akin to the derivation of Dy in (12). Let us exhibit the action of Dφ(t, x), for x = x(ξ) and t ∈ I, onto a tangent In other words, the components (v 1 , v 2 ) are preserved whenever a tangent vector is mapped from S 2 to M t via the differential (20).
As a matter of fact, given a tangent vector fieldṽ = v i ∂ i x on S 2 , the differential Dφ gives rise to a unique tangent vector fieldv = v i ∂ i y on M t , see [26,Chapter 8]. Whenever we useṽ andv in the sequel we refer to their unique identification via the differential (20) and callv the pushforward ofṽ. At this point, the reader might find it helpful to have a look at Fig. 5.
With the above definitions at hand we are able to relate the surface integral (18) to an integral on S 2 via a change of variables. The key is to compute a meaningful surface element as |det(Dφ)| is the magnitude of the change of the volume element.
The following lemma provides the required form.
and letf : M → R andρ : I × S 2 → (0, ∞) be as above. Then, for t ∈ I, Proof. Let us denote byẽ 1 (ξ) andẽ 2 (ξ) the orthogonal unit vectors on S 2 in direction of ξ 1 and ξ 2 , respectively, which are obtained by normalising the coordinate basis Moreover, a straightforward calculation gives Dx Dx = 1 0 0 sin 2 ξ 1 and thus, the surface gradient ofρ in spherical coordinates (17) is given by where we have replaced the coordinate basis with the orthonormal basis.
Using Dx Dx in (19), the Jacobian Jy can be written as Here, we have omitted the argument (t, x(ξ)) of ∇ S 2ρ. Then, the integral turns out to be where the last equation follows from (18) if M t = S 2 , the fact that sin ξ 1 ≥ 0, and The concepts introduced above, and further properties thereof, may be found in any standard differential geometry book. For instance, in [9,10,25,26].

Vectorial Sobolev Spaces on Manifolds.
We briefly introduce the appropriate function spaces required for the variational optical flow formulation on Riemannian manifolds. Again, let us consider time t ∈ I arbitrary but fixed and recall that, for a tangent vector fieldv(t, ·) on M t , the component-wise radially constant extension ofv is denoted byv, cf. (6).
We denote by ∇ûv(t, x) the covariant derivative ofv at x ∈ M t along the direction of a tangent vectorû ∈ T x M t . It is defined as the tangential part of the usual directional derivative of the extensionv alongû in the embedding space, that is, and in matrix-vector form reads It is a linear operator ∇v(t, x) : T x M t → T x M t and its Hilbert-Schmidt norm is given by where {ê 1 ,ê 2 } denotes the orthonormal basis of the tangent space T x M t , cf. (13). We stress that (24) is invariant with respect to the chosen parametrisation.
For each t ∈ I, we define the Sobolev space H 1 (M t , T M t ) as the completion of the space of vector fields where the surface integral is defined in (18). Let us emphasise that (25) is indeed a norm whenever M t is diffeomorphic to the 2-sphere. The reason is that, by virtue of the Hairy Ball Theorem, there is no covariantly constant tangent vector field but v = 0, see e.g. [17, p. 125]. From now on we will omit the arguments (t, x) inside norms for the sake of simplicity whenever clear from the context. Alternatively, one can define Sobolev spaces of vector fields such that each component of a vector field originates from a scalar Sobolev space. See, for instance, Lefèvre and Baillet [27]. On the 2-sphere, however, they are typically introduced by means of the spherical Laplace-Beltrami operator, see e.g. [
We refer to Theorem 5.6 and Lemma 5.8 in [30, Sec. 5.1] for detailed proofs of the previous statements.
The set Ỹ n,j : n ∈ N 0 , j = 1, . . . , 2n + 1 is a complete orthonormal system of L 2 (S 2 ) with respect to the inner product ·, · L 2 (S 2 ) on S 2 . In further consequence, for a functionf ∈ L 2 (S 2 ), we have the Fourier series representatioñ Again, we refer to [30,Sec. 5.1] for the proofs, in particular to Theorem 5.25. In the present article we employ fully normalised spherical harmonics. For the explicit construction see [30,Sec. 5.2]. Moreover, the norm of L 2 (S 2 ) is readily stated in terms of the coefficients in the above expansion via Parseval's identity For an arbitrary real number s ∈ R, we define the Sobolev space H s (S 2 ) as the completion of all C ∞ (S 2 ) functions with respect to the norm We stress that, by (10), ∆ S 2f = −∆ R 3f and we have λ n ≥ 0 for all n ∈ N 0 yielding a sound definition. Accordingly, for s ∈ R, we define the H s seminorm of order s by OPTICAL FLOW ON EVOLVING SPHERE-LIKE SURFACES 13 Now that the space L 2 (S 2 ) is endowed with a basis, we can proceed to define an orthonormal system for square integrable tangent vector fields on the sphere. This will immediately allow us to treat vector-valued problems consistently.
By definition,ỹ n is a normal field whereasỹ (2) n andỹ (3) n are tangent vector fields. Consequently, the latter are called tangential vector spherical harmonics. Note that, by means of the Hairy-Ball Theorem, no tangential vector spherical harmonics of degree zero exist.
In further consequence, let us denote by L 2 (S 2 , T S 2 ) the space of square integrable tangent vector fields on S 2 equipped with the inner product Here, dS 2 denotes the usual spherical surface measure, see also Lemma 2.1.

Optical Flow on Evolving Surfaces.
3.1. Generalised Optical Flow. Optical flow models are typically based on the assumption of constant brightness. Given a sequence of (planar) images γ(t, ξ)) stays constant over time when moving along a trajectory γ(·, ξ) : I → Ω starting at ξ ∈ Ω. In other words, in the planar setting, we have n,j via the differential Dφ Table 1. Summary of notation used throughout the paper.
which is termed optical flow equation and must hold for all ξ ∈ Ω and all t ∈ I. For the sake of consistency, we denote by ∂ t the partial and by d/dt the total derivative with respect to time.
It is possible to generalise the idea to a non-Euclidean setting where the image lives on a, potentially moving, manifold. To this end, let us be given an evolving surface specified by a parametrisation y : I × Ω → R 3 as in (4) together with a functionf , its domain being M. For a time t ∈ I, is then an image on the surface. Adapting the above idea of constant brightness to the new setting requires that, along a smooth trajectory γ(·, x) : t → γ(t, x) ∈ M t that starts at x ∈ M 0 and always stays on the surface, we must havê However, in order to proceed as above one needs to define a meaningful derivative with respect to time.
One possibility, which is pursued in [22,24], is to consider derivatives along trajectories following the moving surface. Let y be as above and let ∂ t y =V be the surface velocity, its domain being t∈I ({t} × M t ) ⊂ R 4 . We emphasise that V is in general not tangent to M t , t ∈ I, and hence in our notation denoted by a boldface capital letter. Then, is the time derivative off at x 0 = y(t 0 , ξ) along the parametrisation y(·, ξ). As a consequence, one can deduce that is the time derivative off in normal direction. It is defined analogously to (31) albeit following a trajectory ψN through x 0 ∈ M t0 for which ∂ t ψN(t 0 , x 0 ) is orthogonal to T x0 M t0 . From that one can immediately formulate the above idea of constant brightness (30) along γ. To this end, we define byM := ∂ t γ the velocity of a point moving along the trajectory γ and demand that must hold. Equation (32) is a generalised optical flow equation. In Fig. 4 we sketch the various trajectories through the evolving surface and their corresponding velocities.
Since we are, however, interested in a coordinate representation of γ, we define a family of trajectories β : I × Ω → Ω such that γ(t, y(0, ξ)) = y(t, β(t, ξ)) holds for all t ∈ I and all ξ ∈ Ω. In other words, we want the composition of β with y, and γ to coincide. By taking the total derivative d/dt on both sides of the above equation we get ∂ t γ = ∂ t y + ∂ t β i ∂ i y.
Let us denotev := ∂ t β i ∂ i y and recall that ∂ t y =V is the surface velocity. The above relation states that the total velocityM = ∂ t γ along a level line of constant intensity isM =V +v (33) andv is a tangential velocity relative to the prescribed surface velocityV. Solving the generalised optical flow equation (32), however, is inconvenient as ψN and, in further consequence, dN t is unknown or hard to estimate. Nevertheless, one can relate (32) and (33), as shown in [24,Lemma 2], and arrive at the parametrised optical flow equation dV tf + ∇ Mf ·v = 0. (34) Solving for the optical flow then means finding a (time-varying) vector fieldv that is tangent to the surface at all times and satisfies the above equation at every point x ∈ M on the moving surface.
Let us conclude this subsection with the following remarks. In general, there exist infinitely many parametrisations y for a given evolving surface M. Nevertheless, the solutionv to the optical flow problem is independent of the choice of the parametrisation. See also [7,Sec. 3] regarding invariance with respect to reparametrisations.
The actual surface velocityV might be unknown or cannot be estimated from the data, as it is the case in this work. As a remedy we propose to decompose the total motionM into a presumed velocityV and a tangential componentv relative to it. Ideally,V is chosen so that it is consistent with M. We refer to [24] for a detailed discussion.
Moreover, we stress that the sought tangent vector fieldv depends on the chosen V and should be interpreted with care. However, it is reasonable to assume that the total motionM is close to the true velocity of a cell. The actual trajectories γ can be reconstructed by finding the integral curves of (33). For this precise approach we point the reader to [24].

Variational Formulation.
The optical flow equation (34) derived above is underdetermined and, in general, a unique solution is not ensured. A common technique to deal with non-uniqueness is Tikhonov regularisation, where one finds a minimiser of dV tf + ∇ Mf ·v 2 L 2 (M) + αR(v). Here, R(v) is a regularisation functional and α > 0 a regularisation parameter, balancing the two terms. The first term is typically referred to as data term whereas the second is called smoothness term. The latter enforces uniqueness and incorporates prior knowledge about favoured solutions.
A common choice for R(v) is the squared H 1 Sobolev seminorm, involving first derivatives with respect to space and time. It favours spatial as well as temporal regularity and is of particular interest when trying to estimate trajectories of objects, albeit computationally more demanding. See, for example [24,42] and [7].
Alternatively, one can omit temporal regularisation leading to a regulariser of the form which is equivalent to solving for each time instant separately. It resembles the original formulation in [18] and its extension to 2-manifolds [27]. In the present article we follow this approach and attempt in finding the unique minimiserv ∈ T M t of the energy Dy(t, ·) y(t, ·) y pỹpŷ p Figure 5. Commutative diagram relating spaces Ω, S 2 , and M t , and tangent vector fields. We highlight that y p is the coordinate representation, see Sec. 2.1, of a particular tangential vector spherical harmonicỹ p andŷ p is its uniquely identified tangent vector field on M t .

Finite-dimensional Projection.
For the subsequent discussion we let t ∈ I be arbitrary but fixed and assume to be given a parametrisation y(t, ·) : Ω → M t as defined in (4). We defer the question of how to find it to Sec. 4.3. Moreover, for notational convenience, we relabel the set of tangential vector spherical harmonics (28) using a single index letter p ∈ N. For instance, for the expansion of a tangent vector field on S 2 we simply writeũ = p u pỹp , where u p ∈ R are the coefficients. We intend to approximate the solution of the problem min v∈H 1 (Mt,T Mt) where E α is defined in (35). We define this space as where J U ⊂ N is a finite index set andŷ p is the pushforward of a particular vector spherical harmonicỹ p via the differential Dφ, see (20). Figure 5 gives a descriptive view of the relation between the introduced spaces and tangent vector fields. The sought vector field is then uniquely expanded aŝ with v p ∈ R, p ∈ J U , being the coefficients. Minimisation of functional (35) results in a finite-dimensional optimisation problem over R |J U | . Plugging ansatz (36) into (35) and writing out definition (25) of the Sobolev H 1 (M t , T M t ) norm gives By applying the definition of the Hilbert-Schmidt norm (24), using linearity of the covariant derivative ∇ûv with respect tov, and the definition of the norm of 18

OPTICAL FLOW ON EVOLVING SPHERE-LIKE SURFACES
for the regularisation term.
The optimality conditions for the discrete minimisation problem (37) are obtained by taking ∂E α /∂v p = 0, for all p ∈ J U , and are given by In matrix form they read

Rewriting the Optimality Conditions.
Even tough directly solving the derived optimality conditions (38) is perfectly legitimate, we take a different approach. The goal of this section is to rewrite the optimality conditions in terms of quantities defined on the 2-sphere, thereby allowing a more general treatment. On the one hand, we want to deal with all surfaces M t for all t ∈ I in a unified manner and, on the other hand, we aim at evaluating (38) numerically on the (approximated) sphere without having to deal with multiple charts, see e.g. [9]. The following is a straight-forward generalisation of [24, Lemma 2].

Lemma 4.1. Consider time t ∈ I arbitrary but fixed. Letṽ
be two tangent vector fields on S 2 and M t , respectively, such that they are related via the differential (20). Then, the parametrised optical flow equation (34) is equivalent to Proof. According to the definitions (31) and (5), we have and it remains to show the identity where we have omitted the arguments (t, y(t, ξ)) on the left and (t, x(ξ)) on the right hand side, respectively. It follows directly from the coordinate representation of the directional derivatives (15) and (16).
In order to give coordinate expressions for the terms in (38) arising from the regularisation term we locally choose an orthonormal frame {ê 1 (t, ξ),ê 2 (t, ξ)} of the tangent space, see (13). As a consequence, the sought tangent vector fieldv can be written asv = w iê i (39) for some components (w 1 , w 2 ) . The reason for expressing the unknown in an orthonormal frame, rather than the coordinate frame, is to simplify matters with regard to the Hilbert-Schmidt norm (24) of the covariant derivative. However, the chosen Galerkin method expands the unknownv in terms of the pushfoward of vector fields which are defined on the 2-sphere, cf. (36). We necessarily need to establish the relation between the intended form (39) and the expression in terms of the coordinate frame.
(⇒) Supposev = Dφ(ũ). Let us take the inner product withê i on both sides. For the left hand side we havê v ·ê i = w jê j ·ê i = w j δ ji = w i . For the right hand side we first observe that, by inversion of the matrix α in (13), it holds that ∂ j y = (α −1 ) jê . Then, and we conclude that w i = u (α −1 ) i as required.
With the above relation at hand we obtain the following form.
Γ j ik denote the Christoffel symbols with regard to the orthonormal frame {ê 1 ,ê 2 } and are defined as We refer to [24,Lemma 3] for their derivation.
Proof. First let us show that, forû = w jê j as in (39), it holds that By the product rule for the covariant derivative (23), Consider the first term of the sum and letê i be represented in the coordinate basis as in (13). Then, Linearity of the lower argument of the covariant derivative with respect to C ∞ (M t ) functions, cf. (16), yields and by realising that ∇ ∂ k y w j is just the directional derivative (16) along ∂ k y we obtain Moreover, in the second term of the sum in (41) we use definition (40). Thus, by summing up all terms in (41) we obtain Applying the previous lemma gives coefficients D i u j and D i v j in the intended form. Finally, it remains to observe that since by definitionê i ·ê j = δ ij .
Finally, by combining Lemmas 2.1, 4.1, and 4.3 we are able to express the optimality conditions (38) in terms of integrals on the 2-sphere. Thus, we arrive at the optimality conditions whereρ ∇ S 2ρ 2 +ρ 2 arises from the Jacobian (19), see also Lemma 2.1.
The entries of the matrices A, D and of the vector b, respectively, are then given by and 4.3. Surface Parametrisation. In order to actually compute the above optimality conditions it remains to determine the radiusρ : I × S 2 → (0, ∞) in the presumed parametrisation (4). Again, we continue the discussion for one particular but fixed time t ∈ I and drop the argument whenever convenient. Estimatingρ(t, ·) : S 2 → (0, ∞) is closely related to surface interpolation from scattered data. Given noisy dataρ δ and a parameter β > 0, it amounts to finding the unique minimiser of the functional where s > 0 is a sufficiently large real number, cf. definition (27). The first term penalises deviation from the observed data whereas the second term enforces spatial regularity of the solution.
In practice, however, N > 0 evaluations {ρ δ (x i ) : x i ∈ S 2 } N i=1 are given at pairwise distinct points on the 2-sphere. In our particular application these correspond to taking the norm in R 3 of pairwise distinct sampling points lying on the sphere-like We again point the reader to Fig. 3. Furthermore, before turning to the numerical solution of (46), let us briefly discuss the regularity requirements. In [7], the authors demand twice continuous differentiability for both the manifold M t and the map y(t, ·) to obtain well-posedness of the optical flow problem. By definition of the parametrisation (4) we require that ρ(t, ·) ∈ C 2 (S 2 ). As a consequence of Theorem 2.7 in [15, Chapter 2.6] regarding Sobolev embeddings, the space H s (S 2 ) for s > 3 is the appropriate choice, i.e. H s (S 2 ) ⊂ C 2 (S 2 ).
Numerically, we approximate the solution of the problem miñ ρ∈H s (S 2 ) F β (ρ)

OPTICAL FLOW ON EVOLVING SPHERE-LIKE SURFACES
by considering a finite-dimensional subspace Q ⊂ H s (S 2 ) and point evaluations (47). In contrast to above, the space where J Q ⊂ N 0 again is an index set, is spanned by scalar spherical harmonics. The sought function is expanded asρ where the unknowns are the coefficients ρ p ∈ R, for p ∈ J Q . Plugging into (46), applying definition (27), and taking ∂F/∂ρ p , for all p ∈ J Q , gives the optimality conditions Denoting by = (ρ 1 , . . . , ρ |J Q | ) ∈ R |J Q | the vector of unknown coefficients, the equations (48) can be written in matrix-vector form as The entries of the matrix L = (l pq ) pq are the matrix M = diag(λ s 1 , . . . , λ s |J Q | ) is a diagonal matrix, and

Numerical Approximation.
Let us finally discuss the numerical solution of the optimality conditions (42). In particular, one needs to (approximately) evaluate the integrals (43), (44), and (45). Even though integrals on the 2-sphere can be computed exactly and quadrature rules exist up to a certain degree, see e.g. [3,16], we instead prefer to use a triangulation together with an appropriate quadrature. The reason is that numerical quadrature on the sphere would have to be of rather high degree to reproduce small details and features of the data, contrary to the chosen quadrature, which can easily be refined up to the desired precision. Finally let us mention that, for a more accurate evaluation of the integrals, one can introduce an intermediate (radial) map from the polyhedron to geodesic triangles. See e.g. [16,Sec. 7.2]. We use a polyhedral approximation S 2 h = (V, T ) of the 2-sphere S 2 . It is defined by a set V = {v 1 , . . . , v n } ⊂ S 2 of vertices and a set T = {T 1 , . . . , T m } ⊂ V × V × V of triangular faces. Each triangle is most easily parametrised using barycentric coordinates, see e.g. [8,Chapter 5]. We associate with each triangle T i ∈ T a tuple (i 1 , i 2 , i 3 ) identifying the corresponding vertices (v i1 , v i2 , v i3 ), which are arranged in clockwise order. The parametrisation (3) then reads which is referred to as the reference triangle. The gradient matrix of T i is then simply The surface normal is constant on T i and is denoted byÑ i . We approximate all functions on S 2 by corresponding functions on the polyhedron S 2 h . A continuous functionf : S 2 → R is replaced by its piecewise polynomial interpolationf h : Here, {φ j } are N h = 6 quadratic shape functions forming a nodal basis together with nodal points {v j } ⊂ S 2 h andf is the usual radially constant extension, cf. (6) in Sec. 2.1. In other words,f h is both a radial projection from the 2-sphere to the polyhedron S 2 h and to piecewise quadratic functions. Note that the shape functions are defined on the triangular faces T i . Whenever a functionf has a dependence on time we simply compute its approximationf h separately for all times t ∈ I. We point the reader to Fig. 6 for a figurative illustration.
In further consequence, the fully normalised scalar spherical harmonics, which were introduced in (26), are substituted with their corresponding approximations on S 2 h . ForỸ ∈ Harm n , n ∈ N 0 , we havẽ We chose piecewise quadratic approximations forỸ so that we can adequately apply ∇ S 2 h and obtain piecewise linear vector fields. Accordingly, we define approximations of the vector spherical harmonics, introduced in (29), as follows. Proposition 1. LetỸ ∈ Harm n , n ∈ N. The piecewise linear interpolations of the corresponding tangential vector spherical harmonics on a triangular face T i ∈ T arẽ y (2) Their derivation is deferred to the appendix. Without loss of generality, letf h (0, ·) andf h (1, ·) be the approximations of the dataf at two subsequent frames. We define the derivative with respect to time by the forward difference ∂ tfh (·) :=f h (1, ·) −f h (0, ·). Moreover, we replace the surface gradient ∇ S 2f of a function on S 2 with its counterpart ∇ S 2 hf h on S 2 h , which is computed according to (17). The functionρ is obtained by solving (48) and, for numerical computations, is further replaced with its piecewise quadratic interpolationρ h as in (49). Coefficients α j i are computed by the Gram-Schmidt process at the nodal points. For numerical computations piecewise quadratic approximations, as defined in (49), are used.
Finally, for the calculation of the integrals we employ the standard quadrature on triangulated spheres, see e.g. [3,16]. Let ξ c = (1/3, 1/3) be the centroid of the reference triangle Ω. Then, we approximate the spherical integral over a functioñ f : S 2 → R on the 2-sphere by S 2f

Microscopy Data.
The present data consist of volumetric time-lapse (4dimensional) images of a live zebrafish embryo during the gastrula period. These videos were recorded approximately five to ten hours after fertilisation by means of confocal laser-scanning microscopy and feature endodermal cells expressing a green fluorescence protein. As a consequence, these labelled cells are recorded without background and allow for a separate treatment. We refer the reader to [21] for many illustrations and a detailed discussion of the zebrafish's developmental process. Regarding the imaging techniques used during data acquisition we refer to [28] and for the treatment of the specimen we point the reader to [31]. The crucial feature of endodermal cells is the fact that they form a so-called monolayer during early morphogenesis, see [39]. Essentially, it means that the labelled cells do not sit on top of each other but float side by side forming an artificial sphere-shaped layer. It can be regarded as a surface and allows for the straightforward extraction of an image sequence. Clearly, this surface is subject to geometric approximations. For instance, in [23,34] it is assumed an ideal sphere, whereas in [7] and [24] only a fraction of the data is considered and modelled as a moving manifold and a height field, respectively, both possessing a boundary.

Preprocessing and Surface Data Acquisition.
Let us briefly discuss the preprocessing steps required to obtain an image sequence together with the evolving surface. We limit our consideration to two consecutive frames and denote the respective volumetric data by f δ 0 and f δ 1 . For each frame, the approximate surface is found by minimising the functional (46) with approximate cell centres acting as sample points. They appear as local maxima in image intensity and are readily located by Gaussian filtering followed by plain thresholding. However, beforehand the points are centred around the origin by fitting a sphere to the cell centres of the first frame and subsequently subtracting the spherical centre.
The triangle mesh S 2 h is obtained by iterative refinement of an icosahedron that is inscribed in the 2-sphere, see e.g. [8,Chapter 1.3.3]. Every refinement step halves the edge lengths by connecting the edge midpoints and projecting them to the unit sphere. Consequentially, every triangular face is split into four smaller triangles and the total number of faces after k ∈ N 0 subdivisions is 20 · 4 k . In our case, k = 7 refinements are required to resolve the data adequately.
It remains to discuss the acquisition of the approximationsf h (0, ·) andf h (1, ·) on the polyhedron. For a frame j ∈ {0, 1}, we define the value at a nodal point v i ∈ S 2 h in (49) via the projection where ε > 0 is chosen sufficiently large.f δ j denotes the piecewise linear extension of f δ j to R 3 , which is necessary for gridded data. The above projection within the narrow band corrects for small deviations of the cells from the fitted surface. Again, we refer to Fig. 6 for illustration. Finally, all intensities are scaled to the interval [0, 1]. Figure 2 shows two frames of the extracted image sequence defined on the spherelike evolving surface. Figure 7 depicts the same matter but in a top view. For better illustration we have added an artificial mesh. Its radius has been widened by one percent.

Visualisation of Results.
We employ the standard flow colour-coding [6] for the visualisation of the computed vector fields. Its purpose is to create a colour image by assigning every vector a colour from a pre-defined colour disk. The colour associated is determined by a vector's angle and its length. However, it was originally defined for planar vector fields and requires adaptation to our particular purpose of tangent vector field visualisation. To this end, we follow the idea developed in [23] by first projecting each vector to the plane and then rescaling its length. Let us denote by P x 3 : (x 1 , x 2 , x 3 ) → (x 1 , x 2 , 0) the orthogonal projector of R 3 onto the x 1 -x 2 -plane. For a tangent vector fieldv we apply the colour-coding to the planar vector field v P x 3v P x 3v.
It is constructed so that the length of individual vectors is preserved. Subsequently, the obtained colour image is mapped back onto the surface. Clearly, in the above construction, one has to distinguish the cases where x 3 ≥ 0 and x 3 < 0. Moreover, P x 3 is required to be injective in either case. The radius R of the colour disk is chosen to be equal to the longest vector in the respective vector field we attempt to visualise. Table 2 lists all values of R for the different figures in this section. In Fig. 9 we show a colour-coded tangent vector field together with the colour disk.
For simplicity reasons, for image functions as well as surfaces we plot their piecewise linear approximations. Moreover, the visualised vector fields are evaluated at the centroids and result in piecewise constant colour-coded images.

Results.
We performed several experiments on said zebrafish microscopy data. In order to obtain an approximation of the evolving surface, we minimised functional (46) by solving the optimality conditions (48). As mentioned in Sec. 5.2, approximate cell centres serve as input. The parameter of the Sobolev space H s (S) was chosen as s = 3 + , where = 2.2204 · 10 −16 is the machine precision, cf. also the discussion regarding theoretical requirements in Sec. 4.3. The regularisation parameter was set to β = 10 −4 and the finite-dimensional subspace was chosen as Q = span Ỹ n,j : n = 0, . . . , 30, j = 1, . . . , 2n + 1 .  Table 2. Radii R of the colour disks used for colour-coded visualisation of tangent vector fields.
In the second step, we computed a minimiser of functional (35) as outlined in Sec. 4.4. Here, the finite-dimensional subspace was chosen as U = span ŷ (i) n,j : n = 1, . . . , 50, j = 1, . . . , 2n + 1, i = 2, 3 . The linear systems resulting from optimality conditions (42) and (48) were solved by means of the General Minimal Residual Method (GMRES) using an Intel Xeon E5-1620 3.6 GHz workstation equipped with 128 GB RAM. Solutions to (42) and (48) converged within 1000 and 100 iterations, respectively, to a relative residual of 10 −2 . The overall runtime was dominated by the evaluation of the integrals (43), (44), and (45). In our Matlab implementation it amounts to several hours. However, the resulting linear system can typically be solved within seconds. Both implementation and data are available on our website. 1 Figures 8 portrays a minimising function of F β for frames 140 and 141 of the image sequence. The resulting surface is depicted in Fig. 2 Fig. 7. Clearly, it reflects the geometry appropriately and contains the desired cell features, cf. also the unprocessed microscopy data in Fig. 1.

and in
In a second step we solved for minimisers of E α for different values of the regularisation parameter α. The tangent vector field is visualised as discussed in Sec. 5.3. Note that in all figures the colour disk has been scaled for better illustration. In Fig. 9, we illustrate tangent vector fields by means of the colour-coding obtained for α = 10 −2 , α = 10 −1 , α = 1, and α = 10. In Fig. 10 the same results are shown but in a top view.
As the central theme of this article is cell motion estimation of fluorescently labelled cells, we also computed their total motion. Recall from (33) in Sec. 3.1 that it is given byM =V +v, whereV denotes the surface velocity andv the optical flow, respectively. By our choice of the parametrisation (4) of the evolving surface its velocity is given byV = ∂ t y = ∂ tρ x. Figure 11 depicts the surface velocityV in the usual colour-coding and the signed norm, i.e. sign(∂ tρ ) V , indicating radial expansions and contractions. Moreover, we visualise the optical flow for one particular value of α and the total motionM.
For comparison, we chose to compute the optical flow on the static round sphere. As this geometry is just a special case of the proposed framework it amounts to findingρ constant in time and space by choosing Q = span Ỹ n,j : n = 0, j = 1 .
Note that for static surfaces the surface velocityV vanishes, that is,M =v. Moreover, we highlight that for computing the optical flow on the static sphere a more efficient approach was proposed in [23]. However, since in this article a norm equivalent to the Hilbert-Schmidt norm (24) is used, we prefer to compare with the same definition and the same regularisation parameters. Figure 12 visualises the optical flow computed for the static spherical geometry with the same regularisation parameters as in Fig. 10. 6. Conclusion. With the goal of efficient cell motion analysis we considered optical flow on evolving surfaces. As a prototypical example we restricted ourselves to surfaces parametrised from the round sphere and showed that 4D microscopy data of a living zebrafish embryo can be faithfully represented in this way. In contrast to previous works, where only a section of the embryo or a spherical approximation was considered, our approach fully attributes the geometry and models the embryo as a closed surface of genus zero. The resulting energy functional was solved by means of a Galerkin method based on vector spherical harmonics. Moreover, the parametrisation of the moving sphere-like surface was obtained from the data by solving a surface interpolation problem. Scalar spherical harmonics expansion allows to easily meet the smoothness requirements of the surface. Finally, we conducted several experiments based on said microscopy data. Our results show that cell motion can be indicated reasonably well by the proposed approach.
Appendix. It remains to give the calculations regarding the piecewise linear approximations of vector spherical harmonics on S 2 h . Both equations in Prop. 1 follow directly by expanding the definitions of the tangential vector spherical harmonics (29). For the fist identity, that is (51), we havẽ y (2) The second identity, that is (52), follows from the fact that from the definition of the interpolation (50) ofỸ h , and the directional derivative (9) on the triangular face, that is