A Numerical Framework for Efficient Motion Estimation on Evolving Sphere-Like Surfaces based on Brightness and Mass Conservation Laws

by   Lukas F. Lang, et al.

In this work we consider brightness and mass conservation laws for motion estimation on evolving Riemannian 2-manifolds that allow for a radial parametrisation from the 2-sphere. While conservation of brightness constitutes the foundation for optical flow methods and has been generalised to said scenario, we formulate in this article the principle of mass conservation for time-varying surfaces which are embedded in Euclidean 3-space and derive a generalised continuity equation. The main motivation for this work is efficient cell motion estimation in time-lapse (4D) volumetric fluorescence microscopy images of a living zebrafish embryo. Increasing spatial and temporal resolution of modern microscopes require efficient analysis of such data. With this application in mind we address this need and follow an emerging paradigm in this field: dimensional reduction. In light of the ill-posedness of considered conservation laws we employ Tikhonov regularisation and propose the use of spatially varying regularisation functionals that recover motion only in regions with cells. For the efficient numerical solution we devise a Galerkin method based on compactly supported (tangent) vectorial basis functions. Furthermore, for the fast and accurate estimation of the evolving sphere-like surface from scattered data we utilise surface interpolation with spatio-temporal regularisation. We present numerical results based on aforementioned zebrafish microscopy data featuring fluorescently labelled cells.



There are no comments yet.


page 1

page 2

page 3

page 4


Optical Flow on Evolving Sphere-Like Surfaces

In this work we consider optical flow on evolving Riemannian 2-manifolds...

Optical Flow on Evolving Surfaces with an Application to the Analysis of 4D Microscopy Data

We extend the concept of optical flow to a dynamic non-Euclidean setting...

Decomposition of Optical Flow on the Sphere

We propose a number of variational regularisation methods for the estima...

Residual viscosity stabilized RBF-FD methods for solving nonlinear conservation laws

We formulate an oversampled radial basis function generated finite diffe...

An a posteriori error analysis based on non-intrusive spectral projections for systems of random conservation laws

We present an a posteriori error analysis for one-dimensional random hyp...

Temporal Wasserstein non-negative matrix factorization for non-rigid motion segmentation and spatiotemporal deconvolution

Motion segmentation for natural images commonly relies on dense optic fl...

A task-driven implementation of a simple numerical solver for hyperbolic conservation laws

This article describes the implementation of an all-in-one numerical pro...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Recent advances in microscopy imaging techniques allow to study cellular dynamics of biological model organisms in more detail than ever before, see e.g. [30, 31, 36]. Time-lapse volumetric (4D) image sequences of the development of entire living animals can be captured in high resolution and on a sub-cellular scale. However, increasing spatial and temporal resolutions require additional efforts in dealing with the resulting large volumes of data. The need for efficient methods to analyse such data has already been acknowledged and is considered a major interdisciplinary challenge [29, 46].

One promising approach in dealing with image sequences of this type is dimensional reduction. A geometric model of the observed organism is introduced and the captured data is considered only with respect to this geometry, see [23, 47]. These efforts focus on the true shape—or an approximation—of the specimen and thereby reduce the spatial dimension of the data by considering only the restriction, or a suitable projection, to this geometry. Due to the spatial sparsity of the volumetric data the essential information is preserved.

A major gain of this approach is that it can also reduce the computational effort during analysis of the recorded material, see e.g. [33, 34, 35, 37, 47]. In addition, introducing a geometric representation of the specimen allows to compute accurate measurements, such as distances, on curved surfaces rather than in—possibly distorting—projections. For the quantitative analysis of cellular processes this leads to a considerable improvement, see [23].

The zebrafish is a popular and well-established animal research model that can be observed in vivo. Understanding its developmental process is of major interest. We refer to [32] for a detailed discussion and illustrations. Cellular dynamics of endodermal cells are crucial for organ and tissue formation during early development of the organism. Despite its importance, there is a lack of understanding of their migration and proliferation patterns [1, 47]. However, endodermal cells are known to form a so-called monolayer, meaning that they do not stack on top of each other but rather float side by side forming a contiguous single-cell layer [54]. For the purpose of observation, these cells can be fluorescently labelled and recorded separately from the background by means of confocal laser-scanning microscopy. Figure 1 illustrates a section of a captured image sequence containing only the upper hemisphere of the embryo. Shown are nuclei of endodermal cells during the gastrula period forming a round surface in a single-cell layer.

Figure 1: Frames (left to right, top to bottom) of a volumetric zebrafish microscopy image sequence recorded during early embryogenesis. The sequence contains 151 frames recorded at intervals of . Blue colour indicates fluorescence response. As time evolves, the initially spherical yolk develops a clearly visible dent, which is where the embryonic axis forms and cells eventually converge to, see also [32, Figs. 11 and 15]. All dimensions are in micrometer (m).
Figure 2: Frame no. 150 of the zebrafish image sequence. The left image depicts the unprocessed volumetric microscopy data . The curved mesh in the center image illustrates a sphere-like surface fitted to approximate cell centres. The right image shows the surface data obtained by taking the radial maximum intensity projection of onto the surface within a narrow band. For details see Sec. 5. All dimensions are in micrometer (m).

The primary goal of this article is quantitative motion estimation of endodermal cells in fluorescence microscopy data of a living zebrafish embryo. Efficient motion estimation is crucial for the large-scale automated analysis of such datasets and can provide new insights into cellular mechanisms and the dynamic behaviour of cells. See e.g. [2, 10, 42, 45, 47].

We build upon previous work [37] where the deforming single-cell layer is modelled as a closed surface , , of the form

together with a time-dependent function that indicates fluorescence response and is assumed to be directly proportional to the observed intensity. Here, is a radial deformation of the 2-sphere . See Fig. 2 for the general idea and Fig. 3 for a sketch.

The main idea, which was developed in [33, 35], is to conceive the motion of a cell—as it migrates through Euclidean 3-space—only with respect to this moving surface. As a consequence, the velocity of a cell which always stays on this surface can be decomposed into the sum of a—prescribed and in general not tangential—surface velocity and a purely tangential velocity which is relative to . Here, denotes the tangent space at . See Fig. 4 for illustration. In further consequence, one can estimate from the data by solving a parametrised optical flow problem

on this evolving surface. Here, denotes a suitable temporal derivative, the (spatial) surface gradient, and dot the standard inner product. As a result, the velocity of a cell can be estimated as . While is relative to the chosen and should be interpreted with care, it is reasonable to assume that their sum is close to the true velocity of a cell. Integral curves then yield approximate cell trajectories.

In this model, is assumed to satisfy a brightness constancy assumption, which is typical for optical flow-based motion estimation: the intensity is conserved along trajectories of moving points. However, in many situations it is too restrictive and possibly violated, see e.g. the discussion in [14, Sec. 3].

In this article, we address this issue and assume that instead fulfils conservation of mass. We derive a suitable generalisation of the continuity equation to evolving surfaces which are embedded in and obtain the pointwise conservation law

Here, denotes the normal time derivative and the surface divergence, is related to surface curvature, and is the scalar normal velocity of the moving surface. The main advantage, compared to [33, 35, 37], is that one is able to directly infer the entire tangential velocity of cells from the data , where denotes the orthogonal projector onto the tangent space of . The normal component of is prescribed by the surface’s normal velocity and the total velocity can thus be estimated by adding the tangential part .

In view of the ill-posedness of above-mentioned conservation equations, we follow a variational approach and minimise a Tikhonov-type functional of the form

where is the squared norm of the left hand side of one of the above identities, is a regularisation functional, and a parameter balancing the two terms. One of the major advantages of applying this energy to mass preservation is that it favours regularity of rather than regularity of , which depends on the prescribed—and in practice often unknown—surface velocity.

We address in this work another major point. While dense motion estimation is often desired for complex natural scenes, it is redundant for aforementioned microscopy data. Such data are considerably simpler due to the characteristic shape of cell nuclei, the absence of occlusions, and their sparsity. In order to mitigate undesired fill-in effects of quadratic regularisation functionals to areas where is zero, we introduce novel regularisation functionals inspired by image segmentation models, see e.g. [7]. Given a segmentation of the cells, motion is only estimated in regions where data is present.

1.1 Contributions

The contributions of this article are as follows. First, we discuss and introduce brightness and mass conservation laws on evolving surfaces, and relate these two concepts to each other. While conservation of brightness is the foundation for the optical flow equation and has been dealt with in [33, 35], we generalise in this article the principle of mass conservation to time-varying surfaces which are embedded in Euclidean 3-space and derive a generalised continuity equation. For numerical convenience we devise a parametrised version thereof.

Second, we propose new spatially varying regularisation functionals for motion estimation based on the discussed conservation laws. They are specially tailored to mentioned fluorescence microscopy data and indicate motion only in regions with cells present.

Third, for the numerical solution we propose a Galerkin method based on compactly supported (vectorial) basis functions. Resulting sparsity effects lead to vast improvements in performance compared to [37]

, which uses globally supported basis functions. Moreover, we provide a formula for the Hilbert-Schmidt norm of the covariant derivative of a vector field, which is commonly used for tangent vector field regularisation. As a result, the Gram-Schmidt orthonormalisation of the tangent basis is rendered redundant, yielding another major performance gain compared to

[8, 35, 37].

Fourth, for extracting a sphere-like surface together with surface image data from aforementioned microscopy image sequences we propose surface interpolation with spatio-temporal regularisation. Compared to [37], where only spatial regularisation is used, this leads to a more accurate estimation of the surface’s (normal) velocity and the surface data, which—in turn—should improve the accuracy of the computed cell velocities.

Fifth, we present numerical results based on aforementioned zebrafish microscopy data. We compute and compare cell motion estimated by imposing either of the two discussed conservation laws.

1.2 Related Work

Concerning dense motion estimation in , Horn and Schunck [27] were the first to propose a variational approach based on conservation of brightness. They suggested minimising a Tikhonov-type functional with Sobolev semi-norm regularisation, favouring spatially regular vector fields. For a general introduction to the topic see e.g. [4, 5] and for a survey on various optical flow functionals see [55]. Well-posedness of the Horn-Schunck functional was proved by Schnörr [48], where the problem was treated on irregular planar domains and solved by means of a finite element method.

Weickert and Schnörr [57] proposed an extension to the domain . The model includes spatial as well as temporal first-order regularisation and is particularly appealing whenever integral curves are to be recovered. A framework unifying various spatial and temporal regularisers was established by the same authors in [56]. For the comparison of different motion estimation methods an evaluation framework was developed in [6].

Only recently, generalisations to non-Euclidean and non-static domains have received increasing attention. For the purpose of robot vision, optical flow on the static round sphere was considered in [28, 52]. With an application to brain image analysis, Lefèvre and Baillet [40] generalised the Horn-Schunck functional to static surfaces which are embedded in and proved well-posedness. Numerically, the problem was solved on a triangle mesh with a finite element method.

With the aim of analysing cell motion in fluorescence microscopy data, Kirisits et al. [33, 35] considered a generalisation of the Horn-Schunck functional to evolving surfaces with boundary. In particular, in [35] the authors proposed a generalisation of the spatio-temporal model in [57]. Minimisation was performed by solving the associated Euler-Lagrange equations in the coordinate domain with a finite-difference scheme. In [34], they studied several decomposition models for optical flow on the static 2-sphere. The problems were solved by means of projection to finite-dimensional spaces spanned by tangent vector spherical harmonics.

In Bauer et al. [8], optical flow on moving manifolds with and without spatial boundary was investigated. The authors considered product manifolds for which an appropriate Riemannian metric was constructed and well-posedness of their formulation was shown.

In Lang and Scherzer [37], the embryo of a zebrafish was modelled as an evolving sphere-like surface. The generalised optical flow problem was rewritten as an equivalent problem on the 2-sphere and solved by means of a Galerkin method based on tangent vector spherical harmonics. In order to find the sphere-like surface from microscopy data, surface interpolation from approximate cell centres was proposed.

We also refer to [2, 42, 45], where the optical flow was computed to track cells in microscopy data, and to [10], where the optical flow was utilised to infer the motion of neural crest cells in zebrafish microscopy images. Moreover, in Schmid et al. [47], the sphere was used to model the embryo of a zebrafish and the motion of endodermal cells was computed in map projections by means of fluid image registration.

According to [14], Schunck [50] was the first to propose motion estimation in image sequences based on the continuity equation. Since then, it has been utilised in numerous works, as mass preservation is a particularly appealing alternative for fluid motion estimation. For instance, in [9, 14, 60] it is was used to analyse meteorological satellite images. In [59], fluid flow was estimated from image sequences and in [51], the continuity equation was utilised to find the 3D deformation of a beating human heart in tomography images. Moreover, in [3] it was used to analyse blood flow and in [13] fluid flow was estimated by means of an integrated continuity equation paired with second-order regularisation. In [15], they proposed to use the continuity equation for cardiac motion correction of 3D images obtained by positron emission tomography. See also [24] for a survey of variational methods for fluid flow estimation. Finally, in [16], for the purpose of joint motion estimation and image reconstruction, both the optical flow and the continuity equation were used.

The remainder of this article is structured as follows. In Sec. 2, we introduce sphere-like evolving surfaces and their basic properties. Moreover, we discuss vectorial Sobolev spaces on manifolds, and introduce compactly supported (vectorial) basis functions and scalar spherical harmonics. In Sec. 3, we discuss brightness and mass conservation on evolving surfaces and introduce for each conservation principle a variational formulation. Section 4 is dedicated to their numerical solution. We derive necessary and sufficient conditions, which are evaluated and solved on the 2-sphere. In order to find a sphere-like surface from real microscopy data, we propose surface interpolation with spatio-temporal regularisation and discuss its numerical solution by means of scalar spherical harmonic expansion. In Sec. 5 we discuss, compare, and visualise numerical results based on aforementioned microscopy data of a zebrafish. Finally, Sec. 6 concludes the article.

2 Notation and Background

2.1 Sphere-Like Evolving Surfaces

Let us denote by the 2-sphere embedded in the 3-dimensional Euclidean space and let denote the norm of , . Moreover, we denote by

a regular parametrisation of mapping points in the coordinate domain to points on the sphere. The outward unit normal at is denoted by .

Let denote a time interval. We consider a family of closed smooth 2-manifolds and assume that each is regular and oriented by the outward unit normal field , . Furthermore, we assume that admits a smooth and smoothly evolving parametrisation of the form


with being a sufficiently smooth (radius) function. We refer to as evolving sphere-like surface.

Let us denote by a smooth function on , by its coordinate representation, and by its representation on , respectively. For and , they are related by


The partial derivative with respect to is abbreviated by . Accordingly, the partial derivative with respect to time is denoted by . At this point let us clarify further notational conventions. Functions or vector fields for which the domain is will be indicated with a tilde and functions for which the domain is will be indicated with a hat. Their corresponding coordinate representation is treated without special indication.

Moreover, we define a smooth (spatial) extension of to which is constant along radial lines. It is given by


Figure 3 illustrates the setting.

Figure 3: Left: Illustration of a cut through surfaces and intersecting the origin. The extension is constant along the shown radial line. Surface normals are depicted in grey. Right: Vectorial basis functions (red) and (white), and the corresponding zonal function centred at the north pole . The parameters of the corresponding function were chosen as and , cf. Sec 2.3.

In the following, let us consider time arbitrary but fixed. We denote the tangent plane at a point by and the tangent bundle by . The orthogonal projector onto , , is denoted by and is given by


The set


forms a basis for the tangent space at . As a consequence, a tangent vector can uniquely be represented as , with being its coordinate representation. The elements are called components of . For a tangent vector field we define its smooth extension to component-wise and analogous to (3).

In the following we will use Einstein summation convention and sum over each index letter appearing exactly twice in an expression, one as a sub- and once as a superscript. For instance, we will write for the sake of brevity. As a further notational convection, boldface letters are used to denote vector fields. In particular, lower case boldface letters refer to tangent vector fields, whereas upper case boldface letters refer to general vector fields in , with the exception of the parametrisations and . Moreover, we will drop arguments, such as or , whenever clear from the context.

The elements of (5) form the gradient matrix

with being the gradient matrix associated with the parametrisation . See [37, Sec. 2.1] for the derivation. The positive definite matrix is commonly referred to as Riemannian metric and is given by

The elements of its inverse are denoted by

. Both are tensors and obey a transformation law when changing from one basis to another. To this end, let

be an arbitrary basis of at such that


Moreover, let be the inverse of the matrix . Then, for a -tensor of order which is defined in the basis , its representation in the basis is given by


See e.g. [38] for details.

The surface gradient of a function , as given in (2), is defined by

where is the usual gradient of the embedding space and is the extension defined in (3). In particular, for we have . Moreover, for a tangent vector , , we have


see [37, Sec. 2.1].

For a function we define the spherical Laplace-Beltrami in accordance to the surface gradient as


where denotes the Laplacian of and is the extension defined in (3).

For an arbitrary surface embedded in 3-space, the total curvature, which is twice the mean curvature, is defined as


A numerically convenient representation is

where denotes the trace of a matrix. See [38, Chap. 8] for details.

Naturally, the chosen parametrisation of admits a smooth map of the form

The differential of is a linear map given by


see [37, Sec. 2.1] for its derivation. It provides a unique identification of a tangent vector field on with a tangent vector field on and for we have . In other words, the differential acts solely on the tangent basis, cf. [37, Sec. 2.1].

For , the surface integral of a function , as defined in (2), is given by


where is the Jacobian of . See Thm. 3 in [17, p. 88]. Moreover, by [37, Lemma 2.1], we have


For further details on the concepts discussed above we refer the reader to standard differential geometry books, such as [18, 19, 38, 39].

2.2 Vectorial Sobolev Spaces on Manifolds

In the following, we consider and arbitrary but fixed. Recall that denotes the component-wise extension (3) of a tangent vector field on . We define the covariant derivative of at a point along a tangent vector as

In particular, for being an element of the coordinate basis (5) and , it reads in terms of coordinates

see e.g. [38, Lemma 4.3]. Here, denote the Christoffel symbols with regard to the coordinate basis, that is, . Let us denote the above coefficients by


They obey the usual tensorial transformation law (7), see. e.g. [38, Lemma 4.7].

Furthermore, the covariant derivative is a linear operator and its Hilbert-Schmidt norm is given by


where is an arbitrary orthonormal basis of the tangent space . We highlight that (15) is invariant with regard to the chosen parametrisation . The following lemma provides a convenient way for its computation:

Lemma 1.

Let and be arbitrary for some . Then, for , it holds that


where we have omitted the arguments on the left-hand and on the right-hand side.


First, let us show that the right-hand side of (16) is parametrisation independent. To this end, let and be arbitrary bases for such that its relation is given by (6). Then, by [38, Lemma 4.7] and (7), we have the transformation law

for the components (14) of the covariant derivative. Moreover, and transform as and , respectively. Recall that, by definition, and . As a consequence,

Suppose now that is orthonormal so that and . Then,

where the second equality follows from the fact that

and the last equality is by definition (15). Finally, the claim follows from combining these equations. ∎

We define for each the Sobolev space as the completion of tangent vector fields with respect to


Let us add that (17) is a norm whenever is diffeomorphic to the 2-sphere since, by virtue of the Hairy Ball Theorem, no covariantly constant tangent vector field but exists, see e.g. [26, p. 125]. We refer to [22, 53] for more details on Sobolev spaces on Riemannian manifolds.

For a tangent vector field , its surface divergence is defined as


where is an orthonormal basis of the tangent space and is defined analogous to (14), see [38] for details.

2.3 Compactly Supported Basis Functions

Let and let . Then, we define the one-dimensional piecewise polynomial function as

The parameter controls the support and is its degree. For a point we define the -zonal function


which, as a consequence, is compactly supported on . See [21, 49] for further details. Moreover, we define the tangent vector fields


where is the outward unit normal of . See Fig. 3 for illustration.

2.4 Scalar Spherical Harmonics

Let us consider the space of homogeneous harmonic polynomials in which are of degree . We restrict their domain to the sphere and denote this space by . Then, by Thm. 5.6 in [43, Sec. 5.1], we have .

For , an element

is an infinitely often differentiable eigenfunction of the Laplace-Beltrami operator

, as defined in (9), and is referred to as a (scalar) spherical harmonic

. Its corresponding eigenvalue is

, see Lemma 5.8 in [43] for a proof. Moreover, it holds that


where , cf. Thm. 5.9 in [43].

The set is a complete orthonormal system in with respect to . As a consequence, every function can be uniquely expanded in its Fourier series representation as

See Thm. 5.25 in [43] for the details. In this article, we will assume that denote fully normalised spherical harmonics, see [43, Sec. 5.2] for their construction. By Parseval’s identity, we furthermore have

Again, see Thm. 5.25 in [43].

We define the Sobolev space for arbitrary by means of the completion of all functions with respect to the norm

For , we define the seminorm of order by


3 Problem Formulation

Let us consider an evolving sphere-like surface


which is specified in terms of a parametrisation as in (1). Every choice of gives rise to a surface velocity


where . We stress that the velocity depends on the chosen parametrisation of which, in general, infinitely many exist. However, its (scalar) normal component, given by

is intrinsic and thus independent of the choice of , see e.g. [35, Prop. 1]. As a consequence, (24) can be represented as


where is the normal velocity and is a vector field tangent to , .

In the following, we consider smooth trajectories of moving particles (or cells) which always stay on the evolving surface. To this end, we assume the existence of a Lagrangian specification


of the path of a particle which starts at and always stays on the surface. Expressing (26) with the help of a coordinate representation requires that


holds for all . As a consequence of (27) and with the help of (24) we find that


where is a purely tangential velocity. Therefore, the velocity of a particle moving along (26) can be decomposed into the surface velocity , which is prescribed by the chosen parametrisation , and a tangential part relative to it. See Fig. 4 for a sketch.

As a consequence of (26) and (27) we infer that the normal part of the velocity of a particle following equals the normal velocity of the surface movement. In other words,

Suppose now that the evolving surface (23) is embedded in a fluid which moves with a velocity , . For and , we denote the restriction of to the surface by . We stress that this fluid velocity is in general different from the surface velocity , defined in (24).

In the following we assume that a particle of interest following (26) convects with this fluid. In other words, for and we require that


From (28) and (25) we find that


Therefore, the surface (23) must evolve with (scalar) normal velocity . Since the fluid velocity can uniquely be decomposed into a normal and a tangential part, we conclude that the latter is given by


The primary goal of this article is to estimate the motion of cells as they move along trajectories (26) through Euclidean 3-space. The main assumption is that they form a surface structure which is deforming over time and can be estimated from image data . Hence, we focus on estimating rather than and utilise the fact that the unknown can be decomposed as in (30).

In the following we discuss two conceptually different ways of estimating the tangential part of the particle motion. One is based on conservation of the data along paths (26) and leads to a generalised optical flow equation. Given and a surface velocity , one tries to compute a tangential vector field relative to it. This precise approach has been pursued already in [33, 35, 37].

The other idea is based on conservation of mass and leads to a suitable generalisation of the continuity equation to evolving surfaces. Given and only the normal component of the surface velocity, one directly tries to infer the entire tangential part of the particle motion.

The main differences are as follows. First, they differ in the assumptions imposed on . One assumes conservation of brightness whereas the other assumes conservation of mass. Second, in the former the unknown is , whereas in the latter the unknown is . Third, as we employ a variational approach, they differ in their regularity assumptions. The first approach desires regularity of , which depends on the tangential part of the imposed surface velocity , whereas the second enforces regularity of the tangential part of the desired motion.

3.1 Conservation of Brightness

Let us be given a function such that, for time ,

is an image on the surface . In this section we assume that, along a smooth trajectory (26), this data satisfies


for all and all . Typically, this constraint is termed brightness constancy assumption and is the basis for many motion estimation methods.

In order to linearise (32) by differentiation with respect to time, one may consider temporal derivatives along trajectories, see [33, 35]. To this end, we define the time derivative of along a trajectory at as


In further consequence, the time derivative of at along the parametrisation is defined analogously as


For a trajectory that passes through at time and for which is orthogonal to , the so-called normal time derivative of is defined as


The relation between (33) and (35) is given by


See [12, Sec. 3.3] for the details. Figure 4 shows a sketch of the different trajectories introduced above and their velocities.