Decomposition of Optical Flow on the Sphere

12/16/2013 ∙ by Clemens Kirisits, et al. ∙ 0

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.



There are no comments yet.


page 6

page 18

Code Repositories


Optical Flow on Evolving Sphere-Like Surfaces

view repo


Brightness and Mass Conservation Laws on Evolving Sphere-Like Surfaces

view repo


Decomposition of Optical Flow on the Sphere

view repo
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

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

[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.

Figure 1: Frames no. 57 (left) and 58 (right) of the zebrafish microscopy image sequence. Blue colour indicates fluorescence response. The spherical image is obtained by a radial projection of the unfiltered data onto a fitted sphere. Top and bottom row differ by a rotation of degrees around the -axis. See Sec. 5 for more details on the data and preprocessing.

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.

Related Work

The first variational optical flow method is usually attributed to [8], where they used an seminorm for regularisation. In [25] it was shown that this particular choice leads to a well-posed problem. We refer to [3] for a gentle introduction to optical flow, to [31] for an overview of different optical flow functionals and to [4] 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 [15] 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 [8] 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 [15], however, the problem was solved by finite elements on a surface triangulation. Finally, we mention the work by Schuster and Weickert [26], 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 [19], 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 [1]. Hierarchical models, originally introduced in [27] 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 [10].

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 [24], 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.

2 Notation and Background

2.1 Scalar Spherical Harmonics


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 corresponding eigenspaces

have dimension and are mutually orthogonal in . Their direct sum equals

. Every eigenfunction

lies 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.

2.2 Vector Spherical Harmonics

Denote by the outward unit normal to and let

be the surface gradient of .

Definition 1.

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 .

Figure 2: Fully normalised scalar and corresponding vector spherical harmonics of degree . Scalar spherical harmonics are depicted using a blue to red colour bar. Type vector spherical harmonics are visualised with red arrows, type ones with blue arrows. Note that the length of the vectors has been scaled for better illustration.

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


for all and . In Fig. 2 a handful of the elements of both bases (3) and (5) are depicted. Every has the following Fourier series representation

In particular, we have Parseval’s identity

For a comprehensive and unified treatment of both scalar and vector spherical harmonics we refer to [6].

2.3 Sobolev Spaces on the Sphere

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 .

3 Decomposition Models for Optical Flow

3.1 Optical Flow on the Sphere

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


3.2 Regularisation

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 [8]. 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).

3.3 Optical Flow Decomposition

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 [1].


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 [29]. Generalising (13) to a decomposition into , instead of two, constituents is possible, but will not be considered here [7].

Hierarchical Models

Hierarchical image decomposition models have been introduced in [27]. 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 [27] 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.

Helmholtz Decomposition

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

4 Numerical Solution

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.

4.1 Optical Flow

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.

Hierarchical Decomposition

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

or simply

4.2 Decomposition

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 .

4.3 Evaluation of Integrals

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

, . 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

5 Experiments

Figure 3: Top view of frames no. 57 (left) and 58 (right) of the processed zebrafish images. The embryo’s body axis runs from bottom left to top right.

5.1 Description of Microscopy Data

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 [17] and [21], 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 [30]. 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 [11].

5.2 Acquisition of Spherical Data

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 [12], here we consider the unfiltered microscopy data for optical flow estimation.

The support of the obtained data is contained in the northern hemisphere. Thus, it suffices to consider only half of the triangulation leading to a total number of around 164000 faces. Figure 1 depicts a sample of two frames . In Fig. 3, a top view of the two frames is shown.

5.3 Visualisation of Tangent Vector Fields

Figure 4: The first image shows a tangential velocity field with the adjusted colour-coding. The second image depicts the colour space in the unit circle. The third and fourth images portray the vector field and the colour space, respectively, but in a top view.

In order to visualise our results we will apply the standard flow colour-coding [4] 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 [12]. 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
Table 1: Radii of the colour disks used in the different experiments below.

As a second way of illustrating steady velocity fields, we employ streamlines, see e.g. [33]. 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


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 [13], 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.

5.4 Experimental Results

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.222

5.4.1 Optical Flow

Figure 5: Minimiser of (top), with , , and increasing values , , and from left to right. The middle row depicts the curl-free component whereas the bottom row depicts the divergence-free component. The embryo’s body axis roughly runs from bottom left to top right in all images.
Figure 6: Streamlines illustrating the velocity fields from Fig. 5. With increasing parameter the streamlines change colours from yellow (bright) to green (dark). The top row shows the total motion whereas middle and bottom rows depict the curl-free and divergence-free parts of the Helmholtz decomposition, respectively. Images are arranged in accordance with Fig. 5.

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 .

Figure 7: Colour-coded velocity fields for sequence are shown in the top row. Parameters were chosen as and , , and and are arranged from left to right. The bottom row depicts the corresponding streamlines.

5.4.2 Decomposition

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.

Figure 8: Decomposition of the total motion (left) into structural parts (middle) and (right). Sequences were chosen as and and parameters were set to , , , and .
Figure 9: Decomposition of the total motion (left) into structural parts (middle) and (right). Sequences were chosen as and and parameters were set to , , , and .

Figure 10: Detailed view of a cell division. The same parameters as in Fig. 9 are used. The top left image depicts with superimposed. The top right image shows . The bottom left image illustrates whereas the bottom right image shows . For better illustration, and have been lightened up and vectors have been scaled.

5.4.3 Hierarchical Decomposition

Figure 11: Solutions of velocity fields obtained by the hierarchical decomposition. Images are aligned from left to right and top to bottom. At iteration the sequence was set to with parameters and .
Figure 12: Solutions of velocity fields obtained by the hierarchical decomposition. Images are aligned from left to right and top to bottom. At iteration the sequence was set to with parameters and .

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.

6 Conclusion

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).


  • [1] J. Abhau, Z. Belhachmi, and O. Scherzer. On a decomposition model for optical flow. In

    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.
  • [2] F. Amat, E. W. Myers, and P. J. Keller. Fast and robust optical flow for time-lapse microscopy using super-voxels. Bioinformatics, 29(3):373–380, 2013.
  • [3] G. Aubert and P. Kornprobst. Mathematical problems in image processing, volume 147 of Applied Mathematical Sciences. Springer, New York, second edition, 2006. Partial differential equations and the calculus of variations, With a foreword by Olivier Faugeras.
  • [4]