Optical Flow on Evolving Sphere-Like Surfaces
We propose a number of variational regularisation methods for the estimation and decomposition of motion fields on the 2-sphere. While motion estimation is based on the optical flow equation, the presented decomposition models are motivated by recent trends in image analysis. In particular we treat u+v decomposition as well as hierarchical decomposition. Helmholtz decomposition of motion fields is obtained as a natural by-product of the chosen numerical method based on vector spherical harmonics. All models are tested on time-lapse microscopy data depicting fluorescently labelled endodermal cells of a zebrafish embryo.READ FULL TEXT VIEW PDF
In this work we consider optical flow on evolving Riemannian 2-manifolds...
The aim of this paper is to discuss and evaluate total variation based
We extend the concept of optical flow to a dynamic non-Euclidean setting...
In this work we consider brightness and mass conservation laws for motio...
The proposed method uses live image footage which, based on calculations...
We address the problem of motion estimation in images operating in the
The goal of this paper is propose a mathematical framework for optical f...
Optical Flow on Evolving Sphere-Like Surfaces
Brightness and Mass Conservation Laws on Evolving Sphere-Like Surfaces
Decomposition of Optical Flow on the Sphere
Motion estimation is a fundamental task for the analysis of spatiotemporal data, the prototypical example of which are sequences of images taken by a camera. The term optical flow has been coined to designate the apparent motion in such data. Its accurate and efficient estimation has been a major topic in the fields of computer vision and image processing for more than 30 years. However, the applicability of optical flow algorithms is by no means limited to flat two-dimensional projections of real world scenes. The advance of microscopy techniques has led to a particularly promising application of optical flow: cell motion analysis. Reliable optical flow algorithms supplied with microscopy images of sufficiently high spatial and temporal resolution can obviously help understanding cellular dynamics in transparent organisms, see for example
has been coined to designate the apparent motion in such data. Its accurate and efficient estimation has been a major topic in the fields of computer vision and image processing for more than 30 years. However, the applicability of optical flow algorithms is by no means limited to flat two-dimensional projections of real world scenes. The advance of microscopy techniques has led to a particularly promising application of optical flow: cell motion analysis. Reliable optical flow algorithms supplied with microscopy images of sufficiently high spatial and temporal resolution can obviously help understanding cellular dynamics in transparent organisms, see for example[2, 18, 23, 24].
The particular dataset we are working with in this article depicts a living zebrafish embryo during early embryogenesis. Main feature of this dataset are the embryo’s endodermal cells, which have been labelled with a fluorescent protein and are known to develop on the surface of the zebrafish’s spherical yolk, see Fig. 1 and Sec. 5.1. The distribution of these cells can be modelled by a nonnegative function depending on time and position on the -sphere, such that the number is directly proportional to the fluorescence response of a point at time . The models we propose below are motivated, but surely not restricted, to this specific type of data. We argue that the problem of extracting and analysing motion from spherical data is sufficiently general so as to be of potential interest to a wider audience.
Our motion models are based on the optical flow equation
which every vector field describing the temporal evolution of should (approximately) satisfy for all and . Here denotes the surface gradient on the sphere. We derive this equation in more detail in Sec. 3.1. Directly solving the optical flow equation for is infeasible. We therefore use Tikhonov regularisation to compute an approximate solution to (1). Tikhonov regularisation consists in minimising a functional which is a weighted sum of two terms. The first one, usually called data or similarity term, is the squared norm of the left hand side of (1). The second term is a regularising functional , which in this article will always be a Sobolev norm (Sec. 3.2). These norms are introduced in Sec. 2.3.
Next, we extend the motion estimation method outlined above to two types of decomposition models that are also variational in nature (Sec. 3.3). While the input of those is again , their outputs are now two or more vector fields capturing different structural parts of the total motion. Both models are adaptations of recently proposed image decomposition techniques to the optical flow setting. The first one is a decomposition. Its idea is to replace in the data term with a sum , and then to add two different regularising functionals, one for and one for . The second one is a hierarchical model. Roughly speaking, one repeatedly minimises an optical flow functional, while the amount of regularisation is constantly decreased. In every iteration, the in the data term is replaced by a sum , where only is optimised and the are the results from previous iterations.
Finally, all optimisation problems are solved by projecting them onto finite dimensional spaces spanned by vector spherical harmonics (Sec. 4). One advantage of this method is that it automatically yields Helmholtz decompositions of all computed vector fields. Mathematical background on (vector) spherical harmonics is presented in Secs. 2.1 and 2.2. In Sec. 5 we provide details of our implementation, give a more detailed account of the used microscopy data, and show experimental results with these data.
To summarise, the main novelties presented in this article are decomposition models for optical flow on the sphere together with their application to microscopy data.
The first variational optical flow method is usually attributed to , where they used an seminorm for regularisation. In  it was shown that this particular choice leads to a well-posed problem. We refer to  for a gentle introduction to optical flow, to  for an overview of different optical flow functionals and to  for a recent survey and benchmark.
Optical flow algorithms have only recently been extended to data defined on non-Euclidean domains. In [9, 28] images defined on the sphere were treated, whereas in  the original functional by Horn and Schunck was generalised to -Riemannian manifolds and well-posedness was verified. Most recently, optical flow on evolving manifolds has been considered in [12, 13].
Horn and Schunck  numerically solved the variational problem by applying a finite difference scheme to the Euler-Lagrange equations. A similar approach was adopted in [12, 13] after parametrisation of the surface. In , however, the problem was solved by finite elements on a surface triangulation. Finally, we mention the work by Schuster and Weickert , where they used projection methods to solve the optical flow equation in the plane. Instead of Tikhonov-regularising their solution, they solely relied on regularisation by discretisation. The main reason for choosing a numerical method based on vector spherical harmonics in this article is that -type regularisers, for arbitrary real , are handled very easily in contrast to most other methods.
The aim of image decomposition models, as pioneered in , is to separate the cartoon and texture parts of images. While the cartoon component should capture large-scale structural components and should therefore be piecewise smooth, the texture component is supposed to consist of high-frequency oscillating patterns. Since the original model was promising but hard to implement, a large number of modifications and approximations have been proposed. In some of them the problematic -norm was approximated by an norm [22, 29]. Recently, models have been extended to the optical flow setting . Hierarchical models, originally introduced in  for image analysis, have not yet been tried in combination with optical flow. They have, however, the preferable property of producing arbitrarily fine multiscale descriptions of input data. As a concluding remark about vector field decompositions, let us remark that Helmholtz-Hodge decompositions of motion fields have enjoyed a certain degree of attention in recent years, not only in the plane [14, 34, 35], but also on surfaces .
Applying optical flow algorithms to cell microscopy data has become increasingly popular lately. See for example [2, 12, 13, 18, 24] and the references therein. We highlight the article , where also endodermal cells of zebrafish embryos have been analysed. There the authors point out that, although of immense importance for developmental biology, only little is known about the motion behaviour of this type of cells.
be the two-sphere embedded in . For functions we define the Laplace-Beltrami operator by
where is the usual Laplacian of and is the radially constant extension of to . The eigenvalues of
. The eigenvalues ofare
The corresponding eigenspaces . Every eigenfunction
The corresponding eigenspaceshave dimension and are mutually orthogonal in . Their direct sum equals
. Every eigenfunctionlies in and is called (scalar) spherical harmonic of degree . Their name derives from the equivalent characterisation of as the restriction to of the space of harmonic polynomials that are homogeneous of degree . From now on
always refers to a particular orthonormal basis of consisting of real-valued spherical harmonics. For numerical experiments we use the so-called fully normalised spherical harmonics, see [20, Sec. 5.2] for a detailed construction.
Denote by the outward unit normal to and let
be the surface gradient of .
Let and . Whenever a function that does not vanish identically admits one of the following three representations
then is called a vector spherical harmonic of degree and type . For obvious reasons we refer to types and as tangential vector spherical harmonics.
Note that there is no tangential spherical harmonic of degree .
We are mainly interested in the space of square integrable tangent vector fields on endowed with the inner product
where is the usual surface measure on the sphere. An orthonormal basis of this space is obtained from (3) by setting
In particular, we have Parseval’s identity
For a comprehensive and unified treatment of both scalar and vector spherical harmonics we refer to .
For an arbitrary real number , the space is commonly defined as the domain of . See [16, p. 37] or [20, Sec. 6.2] for example. In this section we introduce the spaces by means of the vectorial counterpart of .
For tangent vector fields we define the Laplace-Beltrami operator by
where application of to is understood componentwise and is the orthogonal projector onto the tangent plane , compare [6, Def. 5.26]. The tangential vector spherical harmonics introduced in Def. 1 are eigenfunctions of this operator to the same eigenvalues as their scalar counterparts: If we let
for every . The are as defined in (2), the only difference being that now and therefore the spectrum of is strictly positive. Applying functional calculus, we formally define the -th power of by
Finally, for every , set
Note that, in contrast to the scalar setting, this functional is an actual norm, rather than only a seminorm. We therefore define, for every real , as the space of all distributions for which the series in (6) is finite. Clearly, if is any sequence satisfying
for two positive constants , and for all , then replacing with in (6) leads to an equivalent norm and thus to the same space. For every sequence of positive numbers we denote the resulting norm simply by .
Let be a time interval. We assume to be given a scalar time-varying (brightness) function
The problem of estimating optical flow consists in tracking the temporal evolution of the data by means of a time-dependent vector field. Our optical flow model is based on the so-called brightness constancy assumption: We assume existence of a function satisfying
for all and . Intuitively this means that for every starting point on the sphere, the function remains constant along the trajectory . In addition we require that is a diffeomorphism of for every . The first equation in (8) implies that
This equation is typically written in terms of the vector field whose integral curves are the trajectories , which is defined by the equation . The resulting optical flow equation reads
Solving the optical flow equation directly is problematic. In general, a solution to (9) need not exist, and if it exists, it cannot be unique. The typical remedy is Tikhonov regularisation, where one minimises a functional of the form
with being a regularising functional that incorporates a-priori knowledge about desirable solutions. The parameter controls the amount of regularisation. In the context of optical flow one usually tries to enforce spatial (and temporal) smoothness on the solution. A natural candidate for would then be the squared Sobolev (semi-)norm, which penalises first derivatives in space and time equally, compare [13, 32]
Another, and in fact more popular, possibility is to drop time regularisation, in which case minimisation of (10) is equivalent to minimising
for each instant separately. This corresponds to the original approach of Horn and Schunck . From now on we denote the above data term by . Instead of (11) we consider the more general class of optical flow functionals
The regularisation parameter is omitted, as it can be considered a constant factor in the sequence . Functional (12) forms the basic optical flow setting of this article. All variational models considered here are extensions of (12).
The following two decomposition models are inspired by techniques that have recently been developed in the context of image analysis. The fact that motion estimation based on (9) can be viewed as denoising of vector-valued images suggests the translation of said image decomposition models to the optical flow setting .
The aim is now not to extract one, but two vector fields and in such a way that they capture different structural parts of the total motion of . The idea is to solve the following variational problem
|where the functional is defined as|
Choosing, for instance, the two regularisers to be and norms, respectively, would lead to a model which, in spirit, comes closest to the image decomposition model considered in . Generalising (13) to a decomposition into , instead of two, constituents is possible, but will not be considered here .
Hierarchical image decomposition models have been introduced in . There, an original image is decomposed by repeatedly applying denoising steps. The input of one such step is the residual of the previous one. In every step the degree of regularisation is decreased. In contrast to models, hierarchical decomposition models provide multiscale descriptions of the input data.
This iterative procedure can be transferred to the optical flow setting as follows. Let , , be a sequence of norms as defined in Section 2.3, such that
for all and . For every such sequence of sequences we propose the following iterative scheme,
The resulting sequence of accumulated solutions
provides a multiscale representation which, with an appropriate choice of sequences , can be made arbitrarily fine.
The hierarchical model as formulated above is a slight generalisation of the originally proposed one, in the sense that the sequences of regularising functionals considered in  are always of the form . That is, the regularising functional of step is the same as the one of previous steps save for a smaller regularisation parameter. Requirement (14) allows a more general setup.
The Helmholtz decomposition theorem states that every continuously differentiable tangent vector field on the sphere can be uniquely represented as the sum of a consoidal (curl-free) and a toroidal (divergence-free) vector field. More precisely, there exist two uniquely determined tangent vector fields and satisfying
See [6, Sec. 5.3], for example.
The projection method we use to numerically solve the variational problems presented above leads to solutions that are finite Fourier sums
where . Now, from the definition of basis (5), and the fact that for all sufficiently smooth functions , the Helmholtz decomposition of is obtained immediately
In the first subsection below we describe the numerical optimisation of the optical flow functional (12) based on the optical flow equation (9) and explain the modifications necessary for the hierarchical decomposition. Subsection 4.2 is dedicated to the decomposition model. Finally, we explain how the resulting spherical integrals are approximated (Sec. 4.3).
For convenience we relabel the orthonormal basis (5) using a single index and write, for example, from now on.
Let and let be a sequence comparable to in the sense of (7). Then, a minimiser of , if it exists, has to be in . We solve the problem of finding
by a projection method. That is, we let range only over a finite-dimensional subspace of , where
and is a finite index set. The unknown vector field can now be written as
and the problem of finding an optimal simplifies to a minimisation problem over . Plugging (16) into the optical flow functional gives
which is minimal, if the optimality conditions , for all , are satisfied. They read
or in matrix-vector form
where is the vector of unknown coefficients, the elements of matrix read
is a diagonal matrix that corresponds to the regularisation term and the right hand side is given by
With a slight abuse of notation we identified the set with in the definitions of . We continue to do so below.
The hierarchical model only needs a minor modification for the case . We can rewrite the data term from (15) as
where . Therefore, only the right hand side of the optimality system (18) has to be updated in every step. For simplicity we can assume that the approximation space is the same in every step, so that has the representation , where the coefficients are already known from previous steps. Letting denote the right hand side of the optimality system for step , we calculate
The projection approach explained above is easily adapted to the decomposition problem. Again, let be real numbers and choose two sequences , so that is a norm for and is one for . Now, we solve
are finite dimensional spaces. Proceeding as in the previous section, we obtain the following optimality conditions
where the coefficients and are as defined in (19) and (20), respectively. Concatenating the two coefficient vectors and into a single vector , so that the occupy the first entries while the occupy the last entries, the linear system reads
The matrix is given by
and is concatenated from two versions of in the same way as .
It remains to discuss the numerical evaluation of the integrals (19), (20). First, we approximate the -sphere with a polyhedron defined by a set of vertices and a set of triangular faces. Each triangle is associated with a tuple identifying the corresponding vertices . How the triangulated sphere is obtained in practice, is explained in Sec. 5.2.
Second, in every experiment data are given only at the vertices and for two time steps and . We set and and extend both functions to all of by linear interpolation on every triangle. Thus we obtain two continuous piecewise linear functions
by linear interpolation on every triangle. Thus we obtain two continuous piecewise linear functions, . The time derivative of is approximated by a simple forward difference
which is again piecewise linear on . The surface gradient is replaced by a vector field that is constant on every triangle. It is given by
where is the height vector of the triangle pointing from vertex to the opposite side, compare [5, Sec. 3.3.3].
Finally, we also replace the fully normalised scalar spherical harmonics by their piecewise linear approximations defined on . As before we obtain piecewise constant approximations of . The resulting approximated integrals read
where denotes the area of , and
The data which motivated the study of the proposed decomposition models are time-lapse volumetric (4D) images. The obtained sequence depicts a live zebrafish embryo during the gastrula period, approximately five to ten hours after fertilisation. With the help of confocal laser-scanning microscopy, endoderm cells expressing a green fluorescence protein were recorded separately from the background. For details on the imaging techniques and the fluorescence marker we refer to  and , respectively.
The sequence obtained by the microscope captures a cuboid region of approximately . The spatial resolution is voxels and the intensity is in the range . A total number of 77 images were taken, one every . In the following, the microscopic data will be denoted by
During this early stage, endodermal cell proliferation is known to take place on a so-called monolayer . In other words, cells move and divide without stacking and, as the yolk is ball-shaped, admit for the extraction of a spherical image sequence. For further explanations and numerous illustrations of the developmental process of zebrafish embryos we refer to .
We extracted spherical images from the dataset by first fitting a sphere to the approximate cell centres in each pair of consecutive frames. For simplicity we restrict our attention to one such pair of frames which we denote by . Cell centres are typically characterised by local maxima in intensity and can be found by applying a Gaussian filter and simple thresholding. Without loss of generality the radius of the fitted sphere is assumed to be . In a second step, we created a point grid starting from an icosahedron inscribed in the sphere. In each iteration every triangular face is split into four sub-triangles by connecting the edge midpoints with each other and projecting them onto the sphere. Thus, the total number of faces is , where is the number of refinements. In our experiments we found that iterations suffice.
In order to project the volumetric time-lapse data , onto the grid , we define
for the said pair of consecutive frames . Here is a piecewise linear extension to of and is sufficiently large. Deviations of the monolayer from a perfect sphere are thereby corrected. The obtained data are subsequently scaled to the range . Note that, in contrast to our previous work , here we consider the unfiltered microscopy data for optical flow estimation.
In order to visualise our results we will apply the standard flow colour-coding  using a colour disk. Figure 4 (rightmost image) depicts this colour space. Each vector is assigned a colour determined by its angle and length. However, this colour-coding is defined for planar vector fields only. As a possible remedy we suggest to first project tangential velocities to the plane and then correct the length. To this end, let be the orthogonal projector of onto the --plane. Accordingly, given a tangent vector field , the planar vector field which we visualise is
This construction is chosen so that it preserves the length of . This additional rescaling is different to . The resulting colour image is finally mapped back onto the hemisphere. Figure. 4 shows a tangent vector field visualised with the proposed approach.111Some figures may appear in colour only in the online version of this article. From now on we will visualise velocity fields only in top view, as in the right hand side of Fig. 4. In addition, for every figure the colour disk’s radius was chosen to be equal to the length of the longest vector under consideration. Specific values of are given in Table 1.
|Figure||5 left||5 mid||5 right||7 left||7 mid||7 right||8||9||11||12|
As a second way of illustrating steady velocity fields, we employ streamlines, see e.g. . In all our experiments we consider time as fixed and compute the optical flow for one pair of frames, cf. Sec. 3.2. Given a tangential velocity field and a starting point , a streamline on is the solution to the ordinary differential equation
is the solution to the ordinary differential equation
Numerically, we approximated (21) by solving
where is a step size, for a number of approximately 1300 initial points and iterations. The step size was chosen as . This use of integral curves is different from , where we computed approximate cell trajectories in a nonsteady velocity field. The visualisation by means of the colour coding is rich in detail and is even capable of indicating individual cell motion. Nevertheless, it fails to deliver intuition about the Helmholtz decomposition. Streamlines provide the anticipated effect.
We performed numerous experiments and minimised functionals (12), (13), and (15) as outlined in Sec. 4 for the two frames shown in Fig. 1 and Fig. 3. In all experiments the finite-dimensional spaces introduced in Sec. 4 were chosen as
All resulting linear systems were solved using the Generalized Minimal Residual Method (GMRES) on an Intel Xeon E5-1620 workstation with RAM. Solutions converged to a relative residual of 0.02 within 100 iterations. The runtime was governed by the evaluation of the integrals, cf. Sec. 4.3, and amounts to approximately five hours for the chosen bases and the chosen triangulation. Nevertheless, once the integrals are computed they can be used in all of the proposed models and the linear systems can be solved in a few seconds for different parameters and different norms. Our Matlab implementation and the data are available on our website.222http://www.csc.univie.ac.at
In the first experiment, we minimised functional as defined in (12) for and different values of and . Figure 5, top row, depicts the optical flow field for and values , , and . The presented results are in accordance with our findings in [12, 13]. As explained in Sec. 3.3, a Helmholtz decomposition is obtained immediately. Figure 5, middle row, shows whereas Fig. 5, bottom row, shows . Furthermore, in Fig. 6, streamlines for the same velocity fields are portrayed and the individual plots are arranged accordingly. In addition, Fig. 7 shows the velocity fields for parameters , , , and .
In a next experiment, we computed a minimiser for functional (13) in order to obtain a decomposition of the optical flow. The sequences and were chosen as and , respectively. In Figs. 8 and 9, the resulting decomposition is shown for two different parameter settings. The motion field in Fig. 8 was obtained with parameters , , , and . As anticipated, and capture different structural parts of the motion. is sufficiently smooth whereas contains spatial oscillations.
The result in Fig. 9 was computed by setting , , , and . Expectedly, the velocity field is, by choice of , smoother than in the previous setting. In addition, Fig. 10 illustrates the characteristics of the velocity fields during a cell division in more detail. While is smooth, clearly indicates the cell division.
As a final experiment, we computed two types of hierarchical decompositions as proposed in Sec. 3.3. First we chose such that is halved in every iteration. Solutions were obtained using parameters and . In Fig. 11, the subsequence is shown. As increases, the motion field expands on the details. In a second run, was set to , decreasing the exponent of by in every step. We iteratively computed solutions with parameters and , which were kept constant this time. Figure 12 depicts the subsequence . Note that in Figs. 11 and 12 the colour-coded visualisation is relative to the chosen subsequence. As an exception, here we allowed a maximum number of iterations instead of for the linear system solve to ensure a relative residual of 0.025 in the first step of the hierarchical decomposition.
We provided a set of variational methods for the analysis of motion fields. While their applicability is limited to data given on the sphere, the proposed models have great flexibility in terms of possible regularising functionals. In fact, the chosen numerical method based on tangential vector spherical harmonics allows for a straightforward usage of Sobolev norms for any real . Combined with both and hierarchical decomposition models, which we adapted to the spherical optical flow setting, this flexibility makes it possible to capture different motion characteristics with ease. Feasibility of the proposed models was verified on a microscopy dataset depicting endodermal cells of a zebrafish embryo.
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).
Energy Minimization Methods in Computer Vision and Pattern Recognition, volume 5681 of Lecture Notes in Computer Science, pages 126–139. Springer-Verlag, Berlin, Heidelberg, 2009.