Advances in laser-scanning microscopy and fluorescent protein technology have increased resolution of microscopy imaging up to a single cell level . They allow for four-dimensional (volumetric time-lapse) imaging of living organisms and shed light on cellular processes during early embryonic development. Understanding cellular processes often requires estimation and analysis of cell motion. However, the amount of data that is recorded is tremendous and therefore in many cases automated image analysis is necessary.
The specific biological motivation for this work is to understand the motion and division behaviour of fluorescently labelled endodermal cells of a zebrafish embryo. Although of considerable importance for developmental biology, knowledge about the migration patterns of these cells is scarce . The dataset under consideration consists of volumetric time-lapse images taken by a laser-scanning microscope. The recorded sequence depicts a cuboid section of said zebrafish embryo, whose endodermal cells express a fluorescent protein. We model this sequence by a scalar function
that assigns to every pair a nonnegative value proportional to the fluorescence response of point at time .
Optical flow methods are used regularly to estimate cellular motion, see Sec. 1.3. Applying them directly to our data to obtain a dense 3D velocity field
is possible but problematic from a computational point of view , even more so if temporal regularisation is to be included. We propose a solution to this by adapting our model according to biological facts about the nature of the marked cells.
Endodermal cells develop on the surface of the embryo’s yolk, where they form a non-contiguous monolayer . Loosely speaking, they only sit next to each other but not on top of each other. Moreover, the yolk’s shape is roughly spherical and deforms over time. This means that the yolk’s surface can be modelled by an embedded two-dimensional manifold , the subscript indicating dependence on time. In practice, can be approximated by fitting piecewise polynomials, for instance, to the cell centres.111Sometimes it is possible to already capture the yolk’s surface with the microscope in a second sequence of images. We do not, however, use such additional data in this article. Consequently it is possible to reduce the data dimension by only considering the restriction of to this moving surface; see Fig. 2. More details on the acquisition and preprocessing of the microscopy data are given in Sec. 5.2. This dimension reduction, in turn, necessitates the development of an optical flow model for data defined on an evolving surface, which is the main contribution of this article.
Let be a fixed instant of time and . Assume a cell located at , indicated by a relatively high value of , moves with velocity . On the other hand, suppose the yolk’s surface has velocity . The purely tangential vector
. The purely tangential vector
describes the cell’s velocity relative to . Put differently, the total observed velocity of a cell is the sum of the surface velocity and the cell’s tangential velocity . Compare Fig. 1. While the former is a quantity extrinsic to the surface the latter is intrinsic. A motion estimation method dealing with the full 4D dataset would directly try to calculate for all . The method proposed in this article, however, only computes the tangential field for a given surface velocity . The total velocity can then be recovered by adding the two vector fields.
In practice the true velocity of a moving surface might not be known and might even be impossible to determine from available data. This is also the case for the microscopy data considered in this paper. Our solution consists in picking one surface velocity that is consistent with , of which there are infinitely many in general, and to estimate the tangent field relative to this chosen surface velocity. While the resulting must be interpreted with care, it is reasonable to assume that the sum is close to the true total velocity . The selected surface velocity ideally strikes a balance between being easy to implement while being not too unnatural. While modelling the optical flow on an evolving surface is the main novelty of this article, from the viewpoint of our particular application, it can be regarded as a subproblem making the computation of 3D velocities feasible, namely by reducing the data dimension while keeping as much accuracy as possible.
The contributions of this article are as follows. First, we formulate the optical flow problem on an evolving two-dimensional manifold and derive a generalised optical flow constraint. Second, we translate the classical functional by Horn and Schunck  and its spatiotemporal extension by Weickert and Schnörr  to the setting of moving manifolds. The associated Euler-Lagrange equations are solved with a finite difference scheme requiring a global parametrisation of the moving manifold. Finally, we apply this technique to obtain qualitative results from the aforementioned zebrafish data. Our experiments show that the optical flow is an appropriate tool for analysing these data. It is capable of visualising global trends as well as individual cell movements. In particular, the computed flow field can indicate cell divisions, while its integral curves approximate cell trajectories.
Finally, we address a point raised in the recent publication by Schmid et al. , who also analysed endodermal cell dynamics in a zebrafish embryo. They approximated the surface by a sphere, used different map projections to reduce the amount of data by one dimension, and subsequently computed cell motion in the plane. They acknowledge, however, the need for more exact, and supposedly slower, imaging techniques that do not discard any 3D information. While our approach still requires the volume data to be projected onto a surface and thus is faster than comparable 3D approaches, it does not require the surface to be very simple — e.g. spherical or planar — or static.
This article is structured as follows. In the next subsection we review related literature. Section 2 is devoted to providing the necessary mathematical background, notations, and definitions. Sections 3 and 4 introduce our variational model of optical flow on evolving surfaces and contain the continuous and discretised optimality conditions, respectively. In Sec. 5 we explain our microscopy data and the necessary preprocessing steps, summarise our approach, and finally present numerical results.
1.3 Related work
Optical flow is the apparent motion in a sequence of images. Its estimation is a key problem in Computer Vision. Horn and Schunck  were the first to propose a variational approach assuming constant brightness of moving points and spatial smoothness of the velocity field. Since then, a vast number of modifications have been developed. See [3, 30] for recent surveys.
Using optical flow to extract motion information from cell biological data has gained popularity over the last decade. See, for example [1, 2, 5, 8, 14, 23, 24, 26, 27]. In these works displacement fields are computed either from 3D images or from 2D projections of the 3D data. While projections can suffer from inaccuracies [26, 27], the extraction of dense velocities from volumetric time-lapse data poses computational challenges . In the present article we avoid both of these problems.
Many natural scenarios are more accurately described by a velocity field on a non-flat surface rather than on a flat domain. With applications to robot vision, Imiya et al. [15, 28] considered optical flow for spherical images. Lefèvre and Baillet  extended the Horn-Schunck method to general 2-Riemannian manifolds, showed well-posedness, and applied it to brain imaging data. They solved the numerical problem with finite elements on a surface triangulation. In all of the above works the underlying imaging surface is fixed over time, while in this paper it is not.
A preliminary version of this paper appeared in . The main differences to the present article are as follows. First, our current implementation allows us to regularise spatiotemporally as well as only spatially. In  we only treated spatial regularisation. Second, the spatial regularisation functional has been improved in the sense that it is now parametrisation invariant. We have also conducted new experiments with the cell microscopy data and, in contrast to , computed approximate cell trajectories. Finally, we added some recent references.
2 Notation and Background
Whenever convenient we make use of the Einstein summation convention. Every index that appears exactly twice in an expression, once as a sub- and once as a superscript, is summed over.
2.1 Evolving Surfaces
Let be a family of compact smooth 2-manifolds indexed by a time interval . Each is assumed to be oriented by the unit normal field . For every and the orthogonal projector onto the tangent plane is given by
We call an evolving surface, if there is a smooth function
such that is a diffeomorphism between and for every , and is the identity on . Note that cannot be unique in general. With every there is associated a surface velocity. Denote the inverse of by . The surface velocity at a point is then defined by
In contrast to the domain of is not , but rather the 3-manifold
In other words, is a Eulerian specification of , while is a Lagrangian one. Even though different functions give rise to different velocities , the normal velocity of is independent of the choice of . That is, . We provide a short proof of this statement in Proposition 1 in the Appendix. Given a Eulerian specification of , we can obtain, at least locally, a Lagrangian one by solving the ordinary differential equation (
, we can obtain, at least locally, a Lagrangian one by solving the ordinary differential equation (3) for with initial condition . From now on we consider and fixed. See Sec. 5.2 for the specific and we use in the numerical computations.
Let be a parametrisation of mapping local coordinates to points of Euclidean space. By composing and we obtain a parametrisation of the evolving surface
With this convention we always have . Differentiation with respect to will be denoted by . The set forms a basis of . Note that this basis is not orthonormal in general. Using dot notation for the standard inner product of , the components of the first fundamental form are given by
The elements of its inverse are denoted by upper indices .
Let be a scalar function and its coordinate representation,222Distinguishing between a surface quantity and its coordinate representation is often avoided. We decided, however, to make this distinction for the data , and only for , as we found it helpful especially in Sec. 3. that is
The integral of over the evolving surface is then given by
where denotes the surface measure.
2.2 Derivatives on Evolving Surfaces
The spatial differential operators introduced below are not different from those on static manifolds. Therefore can be considered fixed in this paragraph.
The surface gradient of is the tangent vector field which points in the direction of greatest increase of . In local coordinates it is given by
where we omitted the arguments on the left and on the right hand side, respectively. The surface gradient is just the tangential part of the gradient. More precisely, if is a smooth extension of to an open neighbourhood of in , then
Note that the last expression does not depend on the choice of .
Similarly, for two tangent vector fields , on the covariant derivative of along is the tangential part of the conventional directional derivative of along . That is
where is an extension of as above and is the Jacobian of applied to . Let and be their representations in the coordinate basis. The covariant derivative then reads
The Christoffel symbols are defined by the action of on the coordinate basis
An explicit expression for the Christoffel symbols in terms of the first fundamental form is given by
Recall that the coordinate basis is in general not orthonormal. In Sec. 4, however, we want to rewrite the regularisation functional in terms of an orthonormal basis in order to simplify subsequent calculations. Therefore we now make the little extra effort of expressing the covariant derivative in terms of an arbitrary, possibly non-coordinate, frame . Writing and in this basis, the corresponding formula reads
For scalar functions like the covariant derivative is just the directional derivative along . It can be computed by using linearity of the covariant derivative with respect to its lower argument
The are the symbols associated to the new frame . In analogy to (8), they are defined by
For an orthonormal frame the following transformation law describes the relation between the two types of symbols
where is the -coordinate of , that is, and is the Kronecker delta. We give a short derivation of the equation above in Lemma 3 in the Appendix.
The covariant derivative of at a point is a linear operator on , mapping tangent vectors to tangent vectors . Its 2-norm (Frobenius norm) can be computed via
where now is an arbitrary orthonormal basis of the tangent space , that is, . Note that, if is a global parametrisation, then we can obtain a frame which is orthonormal everywhere by Gram-Schmidt orthonormalisation of the coordinate basis .
Let . Denote by a trajectory through with . We define the time derivative of following at as333Note that this composition of with is necessary, because the conventional partial derivative is meaningless in general.
There are a few special cases of this derivative that are worth mentioning. Let be a trajectory for which the vector is orthogonal to . The corresponding derivative is called normal time derivative and denoted by
Every Lagrangian coordinate system of engenders a time derivative like (13) in a natural way. For the time derivative of following is defined by
We choose the notation and , because the derivative (13) in fact only depends on the velocity of at , see Lemma 1. Finally, if is parametrised according to (4), which we assume from now on, then . For illustration see Fig. 3.
We stress that if is the physical surface velocity, then is the natural time derivative for functions defined on , since it measures the temporal change along trajectories of surface points. These trajectories are not cell trajectories in general. They coincide only if the cells do not move by themselves and all the motion is surface motion.
With the definitions from above, we have
The main idea in this derivation from  is to consider the normally constant extension of : Let be an open neighbourhood of . If is chosen sufficiently small, it is possible to define a function that is smooth, constant on normal lines through for every , and agrees with on . Therefore
The last equality holds because, by construction, equals and
Finally, note that by definition (3) we have . ∎
Note that, since is tangential, actually only depends on the tangential part of . Here is the orthogonal projector defined in (2).
Let be a tangent vector field on the evolving surface , that is, a function such that
where application of to is understood componentwise. Again we have . A normal time derivative for could be defined as well but will not be needed in the sequel. As in the scalar case, can be considered the natural time derivative for a tangent vector field , if is the physical surface velocity. By setting
we arrive at the following expression for in local coordinates
The new symbols have the explicit representation
which can be verified by taking inner products of both sides of (17) with the coordinate basis vectors.
Again, in order to simplify calculations later on, we want to express this derivative in terms of an orthonormal frame . We have
where the symbols are defined as before and satisfy an analogous transformation law
3 Model Statement
3.1 Generalised Optical Flow Equation
We assume to be given an evolving surface together with a known Lagrangian specification or, equivalently, a Eulerian one . In addition we are given scalar data on which we want to track over time.
Our optical flow model is based on the so-called brightness constancy assumption. For every we seek a trajectory along which is constant. In other words, we assume existence of a Lagrangian specification of such that
This implies that the time derivative of following has to vanish identically. We deduce from Lemma 1 that the following generalised optical flow equation has to hold
Let us continue the discussion of Sec. 1.1. According to our definition of , a cell located at moves with velocity
where is the inverse of , is the total observed velocity of a cell as introduced in Sec. 1.1 and is its velocity relative to . The second equality above is due to decomposition (1). According to our assumptions at the beginning of this section, we consider as given so that the actual unknown is .
The remaining part of this subsection is devoted to rewriting (22) in terms of local coordinates. First, we give an interpretation of the coordinates of with respect to the basis . Let be the coordinate counterpart of , defined by the equation
See also Fig. 4. Taking time derivatives on both sides and dropping arguments yields
since . We can conclude that , which means that is just the 2D velocity of the parametrised trajectory . It remains to rewrite (22) in terms of and .
The optical flow equation (22) is equivalent to
In local coordinates it reads
We prove the assertion in two steps. First we show that
and afterwards rewrite the right hand side in local coordinates.
By Lemma 1 the normal time derivative can be written as
The other summand of (22) rewrites as
Note that is not assumed to be normal to , so that the term does not vanish in general. However, it does appear twice with opposite signs. Finally recall that and by the definition of the first fundamental form
It is worth noting that the parametrised optical flow equation has precisely the same form as the classical 2D equation.
Directly solving the optical flow equation in the new setting is just as ill-posed as it is in the classical setting. We use variational regularisation to overcome this. In particular, we propose to minimise the following quadratic spatiotemporal functional to recover a vector field describing the tangential motion of data .
Here and are regularisation parameters. Recall from Sec. 2 that is temporally regularised according to the assumed surface motion . Functional (24) is a generalisation of the one presented in  for the Euclidean setting.
Moreover, if , minimisation of (24) is equivalent to minimising
for every instant separately. If for all , the functional reduces to that of . The spatial regularisation term as defined in (12) is independent of the chosen parametrisation. This is an improvement over the functional chosen in .
We end this section with a brief explanation, from an applied point of view, of why we regularise with covariant derivatives. Consider as a toy manifold the non-moving unit circle with parametrisation , and tangent basis . Consider the tangent vector vector field , where is a fixed number. This field would describe a uniform translation of data along the circle, and thus should not be penalised by a regularisation term that enforces spatial smoothness. But while conventional differentiation does not yield a vanishing vector field
covariant differentiation does
Here we used the fact that and .
An analogous argument explains our penalisation of of instead of the unprojected derivative .
4 Euler-Lagrange Equations
To simplify matters from now on we will assume having a global parametrisation of and thus a global parametrisation of the whole evolving surface, cf. (4). In addition, we express the functional (24) in an orthonormal non-coordinate basis with
This leads to wearisome calculations at first, which however pay off eventually when we compute the optimality conditions for the coordinates of with respect to this frame. Note that an orthonormal coordinate basis does not exist in general .
In this section we use the following notational convention. First, we identify with . In addition, Latin indices are always understood to run over the set , while Greek indices are reserved for .
4.1 Rewriting Functional (24)
If we set and , the coefficients of above can be rewritten using the unified notation
where and . Consequently, defining the operator , the integrand of the regularisation term becomes a weighted 2-norm of the matrix . The parametrised version of energy functional (24) now takes the following compact form
where and is the first fundamental form as introduced in (5). Observe that the simple form of the regulariser originates from representing and in an orthonormal basis. This also simplifies the computation of the optimality conditions.
4.2 Optimality System
Denote the interior of by . Functional (28) takes the general form
where the Lagrangian is a smooth function of all and . Denote partial derivatives of by subscripts. A minimiser of has to satisfy the following second-order elliptic system
for and where is the outward normal to . The derivatives of the Lagrangian read
4.3 Discretisation and Numerical Aspects
We solve the Euler-Lagrange equations (29) with a standard finite difference scheme. The spatiotemporal domain is assumed to be the unit cube and is approximated by a Cartesian grid with spacing of in the direction of , where . Grid points are denoted by . Thus, refers to the numerical approximation of at . Partial derivatives of the unknowns are approximated using central differences. They read
where the symbols and denote the neighbours of in the grid along coordinates and , respectively. From the choice of the discrete derivatives an eleven-point stencil is obtained; see Fig. 5. Derivatives of the data and the surface parametrisation are handled likewise, using central differences in the interior and inward differences at the boundaries.
However, the resulting (sparse) linear system is underdetermined from equations (29) alone, because the approximations used for the mixed derivatives of refer to points not occurring in any boundary condition. Thus, at every grid point with
additional boundary conditions are needed. At these points we set in the boundary condition (29), which is a vector pointing in the direction of the undetermined grid neighbour. This leads to expressions of the form , which, interpreted as a directional derivative, can be approximated by
5.1 Zebrafish Microscopy Data
As mentioned before, the biological motivation for this work are cellular image sequences of a zebrafish embryo. Endoderm cells expressing green fluorescent protein were recorded via confocal laser-scanning microscopy resulting in time-lapse volumetric (4D) images. See e.g.  for the imaging techniques.
The microscopy images were obtained during the gastrula period, which is an early stage in the animal’s developmental process and takes place approximately five to ten hours post fertilisation. In short, the fish forms on the surface of a spherical-shaped yolk, which itself deforms over time. Detailed explanations and numerous illustrations can be found in . For the biological methods such as the fluorescence marker and the embryos used in this work we refer to .
The captured area is approximately and shows the pole region of the yolk. Figure 2, left column, depict two frames of the raw data. The sequence contains frames recorded in intervals of with clearly visible cellular movements and cell divisions. The spatial resolution of the data is voxels. Intensities are in the range . In the following we denote by
the unprocessed microscopy data approximating from Sec. 1.1.
The important aspect about endodermal cells is that they are known to form a monolayer during gastrulation . In other words, the radial extent is only a single cell. This crucial fact allows for the straightforward extraction of a surface together with an image of the stained cells. Figure 2 illustrates the idea for two particular frames.
5.2 Preprocessing and Acquisition of Surface Data
In this section, we relate the mathematical concepts introduced in Sec. 2 to the 4D microscopic images. We give a concrete global parametrisation suitable for this type of data and discuss the necessary preprocessing steps leading to an approximation of the evolving surface together with an approximation of the scalar quantity .
The first step is to extract approximate cell centres from the raw microscopy data. As the positions of cells are characterised by local maxima in intensity they can be reliably obtained as follows. For every frame, a Gaussian filter is applied to the volumetric data . Then, local maxima with respect to intensity are computed and treated as cell centres whenever they exceed a certain threshold.
Next we fit a surface to the cell positions. For every frame, this is done by least squares fitting of a piecewise linear function combined with first-order regularisation. From that we get a height field which completely describes the discrete evolving surface. Finally, the numerical approximation of is calculated by linear interpolation of
is calculated by linear interpolation ofand evaluation at surface points determined by .
The combination of all processing steps described above turns the original 4D array into two three-dimensional arrays
Let us quickly relate to the quantities introduced in Sec. 2.1. The mapping
where ranges over a grid, is the discrete parametrisation. The corresponding is the function that identifies surface points with identical coordinates. Thus, the surface motion occurs only in direction of . However, we stress that this particular parametrisation was chosen due to the lack of knowledge about the true motion of material points on the surface.
5.3 Solving for the Velocity Fields
After the preprocessing of the microscopy data as explained above, the following steps lead to the desired solution:
From the parametrisation compute approximations of , , , , as explained in Sec. 2. Like all other quantities the are functions of space and time. They can be computed, for example, by Gram-Schmidt orthonormalisation of the coordinate basis .
Compute coefficients (34) of discretised optimality system from the quantities calculated in step 1.
Solve resulting linear system for unknowns , see Sec. 5.5.
Compute relative tangential velocity via (27) and recover total velocity .
In order to illustrate the computed tangential velocity fields we apply the standard flow colour-coding from .444Some figures may appear in colour only in the online version. This coding turns vector fields into colour images according to a particular 2D colour space.
However, we are interested in visualising tangential vector fields on an embedded manifold, which are functions with values in . To be able to apply the colour-coding mentioned above we turn the computed optical flow fields into vector fields in the following way
where is the orthogonal projection onto the - plane. If the scaling factor were omitted, the new vector field would simply be the original one as viewed from above. The reason for including this scaling are vectors having a large -component, which would otherwise seem unnaturally short. Finally, the image resulting from the colour-coding of vector field (31<