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 development often requires estimation and analysis of cell motion. However, the amount of data captured is tremendous and therefore manual analysis is not an option.
The specific biological motivation for this work is to understand the motion and division behaviour of fluorescently labelled endodermal cells of a zebrafish embryo. The marked 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 deforms over time; see Fig. 1.
We take these biological facts into account and restrict our attention to the analysis of cell motion on the yolk’s surface. With this approach it is possible to reduce the amount of data by one space dimension. The resulting problem consists in the estimation of motion of brightness patterns that are restricted to an itself moving surface. We approach this problem by adapting the classical concept of optical flow to the present setting, where the image domain is both non-Euclidean and dynamic. Note that due to the monolayer structure cell occlusions cannot occur. This makes the optical flow field a more reliable approximation to the true motion field.
Our contributions in the field of optical flow are as follows. First, we formulate the optical flow problem on an evolving two-dimensional manifold and give two equivalent ways of linearising the brightness constancy assumption (Secs. 2.1 and 2.2). One uses a parametrisation of the evolving surface, the other one is parameter-independent. Second, we use a generalisation of the Horn-Schunck model to regularise the optical flow field (Sec. 2.3). For a given global parametrisation of the evolving surface, we solve the associated Euler-Lagrange equations in the parameter domain with a finite difference scheme (Sec. 3). Finally, we apply this technique to obtain qualitative results from the afore-mentioned zebrafish data (Sec. 4). Our experiments show that the optical flow is an appropriate tool for analysing these data. It is capable of estimating global trends as well as individual cell movements and, in particular, it is able to indicate cell division events.
1.0.1 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 has been developed. See  for a recent survey.
Miura  observed that until 2005 optical flow has been mostly disregarded as a method for motion extraction in cell biological data. Since then, a few articles have explored this direction: Melani et al.  and Hubený et al.  extended variational optical flow methods to volumetric images to obtain 3D displacement fields. In the former article, the resulting algorithm is also applied to zebrafish microscopy data. Quelhas et al.  use optical flow to detect cell divisions in a live plant root. However, they work with 2D (plus time) data only. Therefore, their approach suffers from errors caused by 3D off-plane motion.
Clearly, certain 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. [7, 16] considered optical flow for spherical images. In a more general setting, Lefèvre and Baillet  extended the Horn-Schunck method to 2-Riemannian manifolds and showed well-posedness. They solve 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.
2 Optical Flow on Evolving Surfaces
2.1 Brightness Constancy
Let , , be a compact smooth two-dimensional manifold evolving smoothly over time. We assume the velocity to be unknown. Moreover, denote by a scalar time-dependent quantity defined on the surface
We begin with a Lagrangian specification of the optical flow field. That is, for every starting point we seek a trajectory where the data are conserved. More precisely, we want to find a function
for all , for all ,
is a diffeomorphism between and for all ,
is fulfilled and which satisfies a “brightness” constancy assumption (BCA)
In classical optical flow computations it is common practice to linearise the BCA by taking its time derivative and to solve the resulting equation for the Eulerian unknown .111To simplify expressions we use Newton’s notation for those time derivatives that correspond to actual velocities, for example . We also take this route, but differentiation of is more involved. Observe, for example, that for an arbitrary and the usual partial derivative
is not well-defined, simply because, in general, is not an element of for all .
Linearisation after pull-back.
Let be a compact domain and
be a parametrisation of the evolving surface. Denote by the coordinate representation of , that is,
be the coordinate counterpart of . This means, if we let , then gives the coordinates of in (see Fig. 2). In other words, we have the identity
which is a coordinate version of the BCA. After differentiation with respect to it becomes
where is the two-dimensional spatial gradient. Note that the last equation is nothing but the classical optical flow constraint (OFC) for Euclidean data and a displacement field .
We turn to our second derivation. While, as pointed out above, the partial derivative is undefined in general, it does make sense to differentiate following the surface movement. Let be a point on and an arbitrary smooth trajectory through the evolving surface satisfying . Now we can compute
to obtain a valid derivative of . Since this time derivative only depends on the vector
. Since this time derivative only depends on the vector, we denote it by . A natural candidate for a trajectory along which to differentiate is given by the parametrisation . Another possible choice would be a trajectory that is normal to . The resulting normal time derivative is accordingly denoted by .
Finally, we also need the surface gradient . If is a smooth extension of to an open neighbourhood of in , then the surface gradient of at is defined as the projection of the three-dimensional spatial gradient onto the tangent plane to
where is the unit normal to . The surface gradient only depends on the values of on the surface; see e.g. [4, p. 389]. Thus, is well-defined.
The spatial and temporal derivatives of introduced above are related in a simple way. As shown in , they satisfy the equality
where is the tangential surface velocity, that is, the projection of onto the tangent plane to . This decomposition of into normal and tangential components is clearly valid for any trajectory in place of , and therefore in particular for the unknown . This means we can use (5) in order to differentiate the BCA (1) with respect to . The resulting OFC reads
We conclude this section with a brief comparison of the two OFCs derived above. We start by showing how to obtain (4) from (6) and vice versa. To this end we again assume the existence of a global parametrisation and rewrite all quantities in (6) in terms of . First observe that, by (3), the velocity of equals the surface velocity plus a purely tangential component
where is the Jacobian matrix of with respect to . On the other hand, by (5), the normal time derivative is equal to the time derivative of following minus its tangential component
Using the last two equations to rewrite the left-hand side of (6) yields
which is already the left-hand side of (4) in terms of . It only remains to observe that and to replace the surface gradient by its coordinate expression , where is the coefficient matrix of the Riemannian metric; see e.g. .
We highlight the qualitative difference between the constraints (4) and (6). Note that in the former the unknown is , while in the latter it is . This means that (4) constrains the motion relative to the tangential surface velocity , while (6) constrains the absolute tangential motion.
From now on we fix an arbitrary and turn to the actual solution of the parametrised OFC for . Recall that with this notation contains the coefficients of the tangential vector field with respect to the tangential basis of . Note also that, by fixing , there is no more time-dependence in our problem which makes it effectively an optical flow problem on a static surface. Hence we omit any reference to from now on and write instead of .
The sought vector field is underdetermined by the OFC alone. We overcome this by minimising a functional that penalises violation of the OFC while imposing an additional smoothness restriction on . More precisely, we adopt a recent extension of the original quadratic Horn-Schunck regularisation to a Riemannian setting . Basically, they propose to minimise
Here, is the regularisation parameter and is the matrix containing the coefficient functions of the covariant derivatives
of . Using the Christoffel symbols (see Sec. 3) associated to the parametrisation the coefficients are given by
Rewriting (7) as an integral over the coordinate domain, we arrive at the functional
where is the Frobenius norm.
3 Numerical Solution
We solve the problem of minimising functional via its associated Euler-Lagrange equations. Regarding the integrand of as a function , they read
where subscripts of denote partial derivatives. The resulting pair of linear PDEs is of the form
The coefficient vectors are rather lengthy functions of the data and metric tensor
and metric tensor, which is why we do not write them out in full here. Letting for simplicity, the natural boundary conditions of the variational problem are
where . In case of a flat manifold, e.g. , the Euler-Lagrange equations (9) reduce to those of the original Horn-Schunck functional and the boundary conditions become the usual homogeneous Neumann ones. For more details on the calculus of variations we refer to .
Due to the nature of the microscopy data (see Sec. 4.1 and Fig. 1), the manifold modelling the deforming yolk is a surface with boundary that is most easily parametrised as the graph of a function . Hence, we set . Accordingly, for the metric we get
where is the identity matrix. The Christoffel symbols turn out to be
is the identity matrix. The Christoffel symbols turn out to be
Partial derivatives of and of the projected data were approximated by central differences. The system (9) with boundary conditions (10) was then solved with a standard finite difference scheme. In the following section numerical results are presented.
As mentioned before, the biological motivation for this work are cellular image data 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  for the imaging techniques. This type of image shows a high contrast at cell boundaries and a low signal-to-noise ratio in general. Our videos 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; see e.g.  for many illustrations and detailed explanations. For the biological methods such as the fluorescence marker and the embryos used in this work we refer to . The important aspect about endodermal cells is that they are known to form a monolayer during gastrulation , meaning that the radial extent is only a single cell. This crucial fact allows for the straightforward extraction of a surface together with a two-dimensional image of the stained cells. Since only a cuboid region of approximately of the pole region is captured by the microscope, this surface can easily be parametrised; cf. Sec. 3. The spatial resolution of the Gaussian filtered images is pixels and all intensities are given in the interval . Our sequence contains frames recorded in intervals of with clearly visible cellular movements and cell divisions.
4.2 Numerical results
In the following we present qualitative results and demonstrate the feasibility of our approach. For every subsequent pair of frames we minimised the functional (8) as outlined in Sec. 3. We chose grid size as well as temporal displacement as and the regularisation parameter was set to . For demonstration purpose we make use of the standard flow colour-coding , which maps (normalised) flow vectors to a colour space defined inside the unit circle. It is easy to see that the same colours are valid all over the manifold due to the parametrisation.
As representative candidates for this discussion we chose the displacement field between frames 57 and 58 for the following reasons. First, the surface is distinctly developed. Second, a considerable number of cells is present in the image, and third, the interval contains cell divisions. Figure 3, left, shows the colour-coded tangential vector field and the colour space whereas Fig. 3, right, displays the same motion field as computed in the parameter space.222Some figures may appear in colour only in the online version. A visual inspection of the dataset shows that cells tend to move towards the embryo’s body axis, which roughly runs along the main diagonal in Fig. 3, right. Clearly, the velocity field is sufficiently smooth and suggests this behaviour in an adequate manner on a large scale. The expected change in orientation along the body axis is well represented by the colour shift from orange-yellow below the main diagonal to purplish blue in the region above. On the contrary, the choice of the regularisation parameter ensures that individual movements are well preserved as can be observed from the image.
Figure 4 gives a detailed view of the section outlined by a (red) rectangle in Fig. 3, right. This section was chosen because it depicts a cell division. Figure 4, left, and Fig. 4, right, display the frames before and after the event, respectively. Moreover, in Fig. 4, left, the velocity field is shown. From the raw data we observed that when a cell actually splits, the two daughter cells drift apart in a angle with respect to the mother cell. The displacement field clearly shows the anticipated pattern caused by the diverging daughter cells. In Fig. 3, right, the event is point up by two areas which are coloured mutually opposite with respect to the colour space. Our results suggest that cell division can be indicated reasonably well by our model. Both implementation and data are available on our website.333http://www.csc.univie.ac.at
Detailed view of a cell division occurring between frames 57 (left) and 58 (right). All vectors are scaled and only every fourth vector is shown. Intensities are interpolated for smooth illustration
Aiming at efficient motion analysis of 4D cellular microscopy data, we generalised the Horn-Schunck method to videos defined on evolving surfaces. The biological fact that the observed cells move along an itself deforming surface allows for motion estimation in 2D (plus time). In the course of this work, we presented two ways to linearise the brightness constancy assumption and showed that one could be obtained from the other and vice versa. The resulting optical flow constraint was solved by means of quadratic regularisation and verified on the basis of the afore-mentioned data. Our qualitative results suggest that both global trends as well as individual movements including cell division are well shown in the surface velocity field. However, so far we only laid the basic groundwork in terms of a mathematical model.
We thank Pia Aanstad from the University of Innsbruck for sharing her biological insight and for kindly providing the microscopy data. This work has been supported by the Vienna Graduate School in Computational Science (IK I059-N) funded by the University of Vienna. In addition, we acknowledge the support by the Austrian Science Fund (FWF) within the national research networks “Photoacoustic Imaging in Biology and Medicine” (project S10505-N20, Reconstruction Algorithms for PAI) and “Geometry + Simulation” (project S11704, Variational Methods for Imaging on Manifolds).
-  S. Baker, D. Scharstein, J. P. Lewis, S. Roth, M. J. Black, and R. Szeliski. A Database and Evaluation Methodology for Optical Flow. Int. J. Comput. Vision, 92(1):1–31, November 2011.
-  P. Cermelli, E. Fried, and M. E. Gurtin. Transport relations for surface integrals arising in the formulation of balance laws for evolving fluid interfaces. J. Fluid Mech., 544:339–351, 2005.
-  R. Courant and D. Hilbert. Methods of mathematical physics. Vol. I. Interscience Publishers, Inc., New York, N.Y., 1953.
D. Gilbarg and N. Trudinger.
Elliptic Partial Differential Equations of Second Order. Classics in Mathematics. Springer Verlag, Berlin, 2001. Reprint of the 1998 edition.
-  B. K. P. Horn and B. G. Schunck. Determining optical flow. Artificial Intelligence, 17:185–203, 1981.
-  J. Hubený, V. Ulman, and P. Matula. Estimating large local motion in live-cell imaging using variational optical flow. In VISAPP: Proc. of the Second International Conference on Computer Vision Theory and Applications, pages 542–548. INSTICC, 2007.
-  A. Imiya, H. Sugaya, A. Torii, and Y. Mochizuki. Variational analysis of spherical images. In A. Gagalowicz and W. Philips, editors, Computer Analysis of Images and Patterns, volume 3691 of Lecture Notes in Computer Science, pages 104–111. Springer Berlin, Heidelberg, 2005.
-  C. B. Kimmel, W. W. Ballard, S. R. Kimmel, B. Ullmann, and T. F. Schilling. Stages of embryonic development of the zebrafish. Devel. Dyn., 203(3):253–310, 1995.
-  J. M. Lee. Riemannian manifolds, volume 176 of Graduate Texts in Mathematics. Springer-Verlag, New York, 1997. An introduction to curvature.
-  J. Lefèvre and S. Baillet. Optical flow and advection on 2-Riemannian manifolds: A common framework. IEEE Trans. Pattern Anal. Mach. Intell., 30(6):1081–1092, June 2008.
-  S. G. Megason and S. E. Fraser. Digitizing life at the level of the cell: high-performance laser-scanning microscopy and image analysis for in toto imaging of development. Mech. Dev., 120(11):1407–1420, 2003.
-  C. Melani, M. Campana, B. Lombardot, B. Rizzi, F. Veronesi, C. Zanella, P. Bourgine, K. Mikula, N. Peyriéras, and A. Sarti. Cells tracking in a live zebrafish embryo. In Proceedings of the 29th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBS 2007), pages 1631—1634, 2007.
-  K. Miura. Tracking Movement in Cell Biology. In J. Rietdorf, editor, Microscopy Techniques, volume 95 of Advances in Biochemical Engineering/Biotechnology, pages 267–295. Springer, 2005.
-  T. Mizoguchi, H. Verkade, J. K. Heath, A. Kuroiwa, and Y. Kikuchi. Sdf1/Cxcr4 signaling controls the dorsal migration of endodermal cells during zebrafish gastrulation. Development, 135(15):2521–2529, 2008.
-  P. Quelhas, A. M. Mendonça, and A. Campilho. Optical flow based arabidopsis thaliana root meristem cell division detection. In A. Campilho and M. Kamel, editors, Image Analysis and Recognition, volume 6112 of Lecture Notes in Computer Science, pages 217–226. Springer Berlin Heidelberg, 2010.
-  A. Torii, A. Imiya, H. Sugaya, and Y. Mochizuki. Optical Flow Computation for Compound Eyes: Variational Analysis of Omni-Directional Views. In M. De Gregorio, V. Di Maio, M. Frucci, and C. Musio, editors, Brain, Vision, and Artificial Intelligence, volume 3704 of Lecture Notes in Computer Science, pages 527–536. Springer Berlin, Heidelberg, 2005.
-  R. M. Warga and C. Nüsslein-Volhard. Origin and development of the zebrafish endoderm. Development, 126(4):827–838, February 1999.