Objective Observer-Relative Flow Visualization in
Curved Spaces for Unsteady 2D Geophysical Flows
Peter Rautek, Matej Mlejnek, Johanna Beyer, Jakob Troidl,
Hanspeter Pfister, Thomas Theußl, and Markus Hadwiger

time

Fig. 1. Observer motion relative to the time evolution of features: Hurricane Isabel in a time-dependent global wind data set.
(Bottom left) The actual path of Isabel (from NHC/Wikipedia). (Top left) Our observer field u automatically follows the motion of Isabel
without explicit tracking of its path. The shown path is simply a path line of u. (Right) Feature-relative visualization, focused on Isabel inthe center, enabling analysis of its time evolution “in place.” The hurricane appears steady, with the Earth moving inversely underneath.
Abstract—Computing and visualizing features in fluid flow often depends on the observer, or reference frame, relative to which the
input velocity field is given. A desired property of feature detectors is therefore that they are objective, meaning independent of the input
reference frame. However, the standard definition of objectivity is only given for Euclidean domains and cannot be applied in curved
spaces. We build on methods from mathematical physics and Riemannian geometry to generalize objectivity to curved spaces, using
the powerful notion of symmetry groups as the basis for definition. From this, we develop a general mathematical framework for theobjective computation of observer fields for curved spaces, relative to which other computed measures become objective. An important
property of our framework is that it works intrinsically in 2D, instead of in the 3D ambient space. This enables a direct generalization of
the 2D computation via optimization of observer fields in flat space to curved domains, without having to perform optimization in 3D.
We specifically develop the case of unsteady 2D geophysical flows given on spheres, such as the Earth. Our observer fields in curved
spaces then enable objective feature computation as well as the visualization of the time evolution of scalar and vector fields, such thatthe automatically computed reference frames follow moving structures like vortices in a way that makes them appear to be steady.
Index Terms—Flow visualization, observer fields, frames of reference, objectivity, symmetry groups, intrinsic covariant derivatives

1

I NTRODUCTION

Flow visualization is concerned with velocity vector fields that describe
fluid motion. Since the definition of velocity is the infinitesimal changeof spatial position with respect to time, this description depends on• Peter Rautek, Matej Mlejnek, and Markus Hadwiger are with King Abdullah
University of Science and Technology (KAUST), Visual Computing Center,
Thuwal, 23955-6900, Saudi Arabia.
E-mail: {peter.rautek, matej.mlejnek, markus.hadwiger}@kaust.edu.sa.
• Johanna Beyer and Hanspeter Pfister are with Harvard University,
Cambridge, MA, USA. E-mail: {jbeyer, pfister}@seas.harvard.edu• Jakob Troidl is with King Abdullah University of Science and Technology(KAUST), Visual Computing Center, Thuwal, 23955-6900, Saudi Arabia,
and TU Wien, Vienna, Austria. E-mail: jakob.troidl@googlemail.com.
• Thomas Theußl is with King Abdullah University of Science and Technology(KAUST), Core Labs, Thuwal, 23955-6900, Saudi Arabia.
E-mail: thomas.theussl@kaust.edu.sa.
Manuscript received xx xxx. 201x; accepted xx xxx. 201x. Date of Publication
xx xxx. 201x; date of current version xx xxx. 201x. For information on
obtaining reprints of this article, please send e-mail to: reprints@ieee.org.
Digital Object Identifier: xx.xxxx/TVCG.201x.xxxxxxx

the notion of what is considered to be “the same position” over time.
Mathematically, this is the notion of a frame of reference, or an observer.
A frame of reference is a way of referencing positions in space overtime, i.e., a time-dependent transformation of space, with respect towhich derivatives are to be taken. (Not an observer as an actual person.)
For this reason, flow visualization as well as the computation offlow features, such as vortices, in general inherently depend on thechosen frame of reference. An important related concept therefore isinvariance with respect to certain changes of reference frame. Galilean
invariance [28] refers to frame changes with constant relative velocity,which is inherently related to inertial frames in physics [63]. Objectivity [31] considers frame changes with time-dependent translations androtations. The standard definition used in visualization and continuum
mechanics is the one of Truesdell and Noll [65], which has received
significant attention for flow fields since the work of Haller [31]. A
major basic limitation of this definition, however, is that it is given onlyfor Euclidean space and cannot be used in general curved spaces.
One area where objectivity is particularly relevant is that of geophysical flows, such as in the Earth’s atmosphere or oceans. The naturalreference frame, in this case, is rotating [51, p.14], i.e., non-inertial,
and Galilean invariance is not well-suited. As is common in flow visualization, we only consider motions by themselves (kinematics), notthe forces generating them (dynamics). We can therefore consider accelerating reference frames just as easily as inertial frames, because thedifferences in the forces acting are irrelevant. For example, in Fig. 2 weconsider an airplane as an observer measuring motion relative to its own
reference frame. Naturally, this does not change a phenomenon such asa hurricane, but it does change all relative velocity measurements.
Although the Earth’s surface is embedded in a flat 3D space (ignoring relativistic effects), it is beneficial to be able to view the domain
where a geophysical flow field is defined from an intrinsic, curved 2D
perspective. This not only allows for simplified computations in 2D
instead of in 3D, or considering different altitude or depth layers separately; it also corresponds perfectly to the reference frame of circularorbits, such as those of satellites or airplanes moving along great circles
around the Earth, as shown in Fig. 2. For visualization purposes, it isalso common to visualize spherical layers. Moreover, scientists haveeven discovered that some geophysical fluid systems, such as aspectsof our atmosphere or even that of Jupiter [29, 56], can sometimes bemodeled in a more realistic way in 2D than in 3D, because the atmosphere is so much thinner than wide. For this reason, it can be crucial
to perform computations in 2D, because the physical behavior of 2D
turbulence can be radically different from 3D turbulence [23, 35, 40].
With this motivation, we develop general techniques for flow visualization and computation of flow features in curved spaces. For practicalrelevance, we illustrate results for geophysical flows on spheres. We
generalize the notion of observer fields [30] from Euclidean space tocurved spaces. This enables “following” flow features such as vortices,and visualizing the corresponding time evolution of scalar and vector
attributes “in place.” For example, Fig. 1 shows hurricane Isabel (2003)
in a time-dependent global wind simulation. Our computed observer
field automatically follows the motion of the hurricane. Fig. 1 (right)
shows feature-relative visualization, where the observer field makes
the motion of Isabel appear to be steady. This enables “feature-centric”
analysis, without having to explicitly follow the motion of the hurricane.
1.1 Observer Fields and Observed Time Derivatives
An observer field u(x,t) is a time-dependent velocity field that, insteadof describing the motion of particles, describes the motion of a continuous field of observers [30]. The path lines of the observer fieldcorrespond to the world lines [3, p.8] of individual spatial locationsin the space of an observer that are viewed as being “the same point”
over time. If an observer field u describes a rigid motion, it is a Killing
vector field [53]. Figs. 2 (a,b,c) depict three time steps of observer fields
on a sphere. Fig. 3 shows one example input vector field v describing asimple rotating vortex, as it is seen by four different observer fields u.
The observed time derivative corresponding to a given u(x,t) measures the differential change of an arbitrary input field v(x,t) withrespect to the observer motion described by u. This concept was introduced in the context of Euclidean space [30]. We define this time
derivative, on any differential manifold M, as the differential operator
D∂:=+ Lu .
Dt∂t(1)

Here, Lu denotes the Lie derivative [22, Ch. 4] with respect to theobserver field u. In order to “steadify” an input vector field v, e.g., tobe able to compute and follow features, we compute observer fields u
by minimizing the derivative (D/Dt)(v − u). If Fig. 3 is interpreted as
showing four different input fields v, the observer-relative visualization
with respect to the optimized field u will always look like Fig. 3 (d).
1.2 Mathematical Framework for Curved Spaces
We develop a novel fully intrinsic mathematical framework for 2D flow
visualization and feature computation in curved spaces. Our major
contribution is the generalization of a prior Euclidean framework [30]for the objective computation of observer fields to curved spaces. We
compute these fields by solving an energy minimization problem, making unsteady input flows as steady as possible. We first introduce a
generalized notion of objectivity for curved spaces. While equivalent
to the standard one in Euclidean space [65], ours is fully intrinsic, using the concepts of continuous symmetry groups and group actions on

(a)

(b)

(c)

Fig. 2. Observers in curved spaces. An airplane (blue path) and asatellite (red path) are moving relative to the Earth. At each instant in
time, with respect to the Earth’s surface both fly in the direction of ageodesic (a great circle). Considering the observer (airplane or satellite)
as stationary, the instantaneous relative motion of the Earth is a Killing
vector field on the sphere, in the opposite direction (right). (a) The Killingfield for the satellite is a steady field (assuming the satellite ground trackis always the same geodesic); (b,c) The field for the airplane is unsteady:Each instant in time is determined by a different tangential velocity.

manifolds from mathematical physics. From general symmetry groups,
we specialize to observer transformations modeled as elements of theisometry group of a manifold M. To simplify computations, we workwith infinitesimal isometries, given explicitly by vector fields on M.
To be able to achieve this, we derive the differential geometric
operators needed to compute observer fields in curved spaces: (1)Observed time derivatives; (2) Approximate Killing fields; and (3) The
velocity gradient tensor ∇v as a covariant derivative [22, Ch. 9], insteadof the Jacobian that is commonly used in flow visualization.
Methodology. Our generalization of objectivity and observer fields
to curved spaces requires mathematical machinery from Riemannian
geometry and mathematical physics that is not common in the flow
visualization literature. For this reason, we summarize background andour notation in the beginning, and provide extensive supplementary
appendixes. These techniques enable deriving a powerful generalizationof previous flow visualization techniques from flat space to curved
spaces in a clean and intrinsic manner. We believe that this powerful
methodology can be very useful for flow visualization in general, andthat it is worthwhile to introduce it to the flow visualization literature.
2 R ELATED W ORK
The framework presented here extends the work of Hadwiger et al. [30]from Euclidean space Rn to curved spaces, focusing on curved twomanifolds. We also model a collection of observers by an observer
velocity field describing their motion, which we also compute by solving a least-squares problem such that (1) it is objective, and (2) the input
flow field is as steady as possible relative to the observer field. However,
the generalization to curved spaces is conceptually and technically challenging. We propose a fully intrinsic generalization, working “inside”
the curved space. This results in natural and efficient generalizations,and facilitates computing 2D observer fields via 2D optimization.
Flow visualization on curved surfaces. All major flow visualization techniques have been extended from flat space to curved surfaces.
For example, LIC [14] and variants [43, 55, 60] have been extended
to curved surfaces, either for parametric surfaces [21], by texturing
triangle meshes [8, 59, 62], or by computing LIC in screen space [6, 36].
Image-based texture advection techniques [37, 67] have been extended
as well, such as image-based flow visualization on surfaces [39, 68].
Image space techniques typically project vector fields from 3D to 2D
screen space for visualization. We refer to the survey by Edmunds etal. [20]. Texture atlases have also been used for flow visualization, forexample flow charts [44]. We similarly use multiple coordinate charts.
Vector field design on surfaces. We extensively use differential
properties of manifolds and vector fields. These are particularly important for computations with vector fields on surfaces, such as vector
field design on surfaces [15, 16, 71]. We refer to the survey of Pengand Laramee [52]. We would like to particularly point out the use ofgeodesic polar maps and parallel transport along geodesics in the work

(a)

(b)

t=0

(c)

(d)

t=1

Fig. 3. A rotating vortex seen by four different observers. The observers rotate leftward around the sphere (here, around a vertical axis), withangular velocities ω = −2π, −π/4, −π/8, 0 [s−1 ] (a,b,c,d). All path lines are on the sphere (the depicted shadows), but we visualize the progression oftime as increasing radial distance and color (see inset). If these are seen as four different input fields v, minimizing the observed time derivative
D/Dt(v − u) gives observer fields u with rightward rotations ω = 2π, π/4, π/8, 0 [s−1 ]. (d) Observer-relative visualization is the same for all four cases.

of Palacios and Zhang [50]. The notion of covariant derivative, which
we employ in our work, is inherently related to parallel transport.
Covariant derivatives have been used for mesh processing [18] andcomputing vector [17] and direction [69] fields on surfaces. An important focus is often their discretization on triangle meshes [5], e.g., byusing discrete exterior calculus [18]. We discretize on triangle vertices
in 2D charts, and particularly focus on fully intrinsic 2D computations.
Interpolation. We use standard barycentric interpolation in triangles. However, spherical barycentric coordinates, i.e., mean value coordinates on the sphere, could also be used for better interpolation [38].
Killing vector fields on (Riemannian) manifolds [19] correspond totheir infinitesimal isometries. Their flows correspond to the intrinsic
isometries of a manifold. Basics of Killing fields for curved surfaces
are discussed by Ben-Chen et al. [9] and by Solomon [57]; a more advanced treatment is given by Petersen [53]. Approximate Killing fields
have been used to compute approximate intrinsic isometries of curved
surfaces [9], or to design approximate Killing fields on meshes [4, 5].
As-rigid-as-possible shape interpolation methods [1] work directly withisometries instead of with derivatives, and are computed extrinsically.
Observers. The concept of an observer as well as that of objectivity
in Euclidean space has been recognized to be of importance in fluid mechanics [31] as well as in flow visualization [28]. A standard referencein continuum mechanics is Truesdell and Noll [65], whose definition
for objectivity in Euclidean space is used often [7, 25, 31, 32, 34, 48].
Vortex detection and frame invariance. Vortices are an important
topic in fluid mechanics [54] and in flow visualization [28]. Many vortex detection methods were originally defined for steady flow, such as
Sujudi and Haimes [61], and were later extended to unsteady flow [70].
Galilean invariance of these critera is considered to be important [28].
Several known methods are Galilean-invariant, such as Weinkauf etal. [70], or the method of Bujack et al. [13]. The latter jointly considers
multiple observers for 2D flow fields. Specific vortex detectors havebeen developed to be objective by design [31, 32]. Günther et al. [7, 25]
present objective methods for Euclidean space that are generic by building on optimization, making non-objective vortex detectors objective.
Other kinds of invariance have been defined, particularly for flow fields,such as rotation invariance [26], or hyper-objective vortices [27].
Flow decomposition methods, such as variants of the HelmholtzHodge decomposition [10, 11], can also be used to remove backgroundflow, corresponding to the harmonic component [12]. These decompositions can be computed on curved surfaces [64]. However, they apply
only to time-independent fields, or alternatively they have to be applied
to each time step individually, without considering time derivatives.
3

Vector fields. A smooth vector field v on a manifold M is a smoothfunction giving a vector v(x) at every point x ∈ M on the manifold, i.e.,v : M → T M,
x 7→ v(x).

(2)

T M refers to the tangent bundle of M, the manifold of all tangent spacesof M, and a vector field is also referred to as a section of T M. Most
important, v(x) is an element of the tangent space at x, denoted by Tx M.
Where no confusion arises, we also denote a single vector v(x) by v.
Coordinate charts map open subsets U ⊂ M of an n-dimensional
manifold M to Rn . Coordinates on a curved manifold M are thus not
given by coordinate vectors, which is only possible in linear spaces,
but for each chart by n coordinate functions {xi }, with each xi : U → R.
For example, for n = 2, the coordinate functions x1 , x2 , or, sometimes,u, v. In contrast, basis vectors live in each tangent space Tx M, at x ∈ M.
Bases. Vectors are geometric objects independent of any chosen
basis. In order to refer a vector v to a basis, we expand it in components
as v = vi ei , where {ei } is a basis for Tx M, and, in general, the Tx M atdifferent x ∈ M have different bases. We employ the Einstein summation convention [22, p.59], implying summation over indices occurring
twice (once “upstairs” and “downstairs” each), e.g., vi ei := ∑i vi ei .
Dual bases. We will also need the concept of a dual basis {ω i },where ω i (e j ) = δ ij , with the Kronecker delta δ ij = 1 if i = j, andzero otherwise. Each ω i is a covector, or 1-form, which is a linear
function mapping a vector to a scalar. The dual basis {ω i } reads off
the components of a vector v referred to {ei }, such that v = ω i (v) ei .
Metrics. Lengths and angles of vectors are determined in an intrinsic
manner by defining a (Riemannian) metric g on M [19, p.35], which isa second-order tensor field defining an inner product hv, wi := g(v, w)
on each Tx M, varying smoothly over M. The vector norm is thendefined as kvk := hv, vi1/2 . Computations with g can use components
gi j := hei , e j i referred to {ei }, computing hv, wi as gi j vi w j . We also
extensively use covariant derivatives compatible with a given metric g.
Velocity gradients. An important derived quantity of a velocity vector field v is the velocity gradient tensor field ∇v. In flow visualization,∇v is often seen as a Jacobian ∂ j vi of partial derivatives in Cartesian coordinates. However, these coordinates do not exist in curved spaces, andin general ∇v 6= (∂ j vi ) ei ⊗ ω j , since ∂ j vi is not tensorial [22, p.241].

BACKGROUND AND N OTATION

Our intrinsic framework enables all computations, most importantly the
optimization of observer fields, to be 2D for curved surfaces, instead ofrequiring 3D ambient space computations. We perform all computations in an atlas of 2D coordinate charts (Fig. 4). This also simplifies
visualization computations on curved surfaces, by working directly incharts, for example for path line integration or Line Integral Convolution (LIC), as depicted in Fig. 3. However, already a surface as simple
as a sphere cannot be covered by a single chart. All computations incharts therefore need to be independent of the choice of coordinates.
We ensure this invariance by using general tensor methods [24], withthe corresponding transformation rules for transitions between charts.

chart
metricFig. 4. Our intrinsic framework employs a perspective from within acurved manifold M. Thus, all computations on curved surfaces are 2Dcomputations in 2D coordinate charts. Each point x ∈ M has intrinsic
differential properties: A (Riemannian) metric, and covariant derivative
operators. A vector field v is a vector v(x) in each tangent space Tx M.
Per chart, each Tx M has a coordinate basis {ei }, and a dual basis {ω i }.
The symmetric metric tensor field is visualized with glyphs on the right.

We note that even in flat R2 with polar coordinates, ∂ j vi is not a tensor.
This has important implications, such as that a tensor that vanishes in
one coordinate system must vanish in all coordinate systems [24, p.82].
This is not true for ∂ j vi , and thus the Jacobian behaves only as it shouldin special cases, such as by restriction to Cartesian tensors [2, Ch.2].
Our framework requires a more general definition of ∇v as thecovariant derivative ∇ j vi , which is standard in mathematical physics.
However, we have found no explicit use of this concept in the flow
visualization literature. We derive this notion for our purposes in Sec. 6.
Approximate observer Killing fields. In general, an observer field
is allowed to be given by an arbitrary vector field u. However, it makes
sense to restrict the motion described by u to be “as rigid as possible.”
We therefore optimize observer fields such that they correspond to
approximately isometric deformations, i.e., to “almost-rigid” motions.
In the same vein as Hadwiger et al. [30], we minimize the rate-ofdeformation tensor of u. However, in curved spaces M this computation
requires using the covariant derivative ∇u, and the metric g on M.
4 O BSERVERS IN C URVED S PACES
The first crucial notion for us to define is the meaning of an observer ina curved space. In analogy with the standard notion of observer transformations in Euclidean space [65], it is natural to consider observertransformations to in general be time-dependent distance-preserving
transformations of the underlying space, i.e., time-dependent isometries. In the general case, however, we have to consider the intrinsic
isometries of a smooth manifold M with metric g,1 corresponding togeodesic distances instead of Euclidean distances. Generally stated,
Remark. We employ a (non-relativistic) concept of frames of reference
postulating that all observers agree on the pairwise spatial geodesic
distances of events happening at the same time, and all observers also
agree on the time difference between events (i.e., time is absolute).
Observer transformations are therefore time-dependent intrinsic distance-preserving transformations. Mathematically, we model theseobserver transformations as paths t 7→ g(t) ∈ G, through a Lie group G,
chosen as the isometry group G := Isom(M) of the manifold M withmetric g. In our framework, we often also work directly with the derivatives of g(t), forming the corresponding path t 7→ X(t) ∈ g through the
Lie algebra g = isom(M), the Lie algebra of infinitesimal isometries.
4.1 Observers on the Sphere
Figs. 2 and 5 illustrate this idea for M = S2 , the standard two-sphere.
The motion of the Earth relative to a fixed observer, such as the airplanein Fig. 2, is given by a path of isometries t 7→ g(t) ∈ SO(3), which here
are different rotations of the Earth. At any instant in time t, the derivative of this path corresponds to a Lie algebra element X(t) ∈ so(3),which defines an infinitesimal isometry of M. However, X(t) is justan anti-symmetric matrix (Fig. 5, bottom). The specific corresponding infinitesimal isometry of the sphere is given by an isomorphism
between so(3) and the Lie algebra of Killing vector fields on the sphere(Fig. 5, top). This isomorphism is mathematically described by a Lie
algebra action (see App. J). In this way, the infinitesimal isometry on
the surface of the Earth is the Killing vector field x(t) corresponding toX(t). At each point x ∈ M, this gives a relative velocity vector x(x,t)(Fig. 5, top). We also know that x describes an infinitesimal isometryof M, because ∇x, as the covariant derivative of x, is anti-symmetric.2
In fact, any Killing field x is uniquely determined by its value x(x) andcovariant derivative (∇x)x at one point x ∈ M [53, Proposition 8.1.4].
4.2 Observer Fields
Generalizing further from the concept of one observer, as one timedependent isometric transformation of space, we want to be able todescribe an observer field—many observers—on a manifold M as avector field u on M, where u does not necessarily correspond to an
exact (infinitesimal) isometry. One of the most important properties of
1 A (Riemannian) metric g is required to define the meaning of isometry.
2 It is important to note, however, that checking ∇x for anti-symmetry needs

to be done with care, because it is a tensor of mixed type. See Sec. 7.1.

Fig. 5. Killing vector fields give the infinitesimal isometries of a manifold
M, here the sphere S2 . The Lie algebra elements X, Xi ∈ so(3) (bottom)
generate the Killing fields x, xi on M (top) through their Lie algebra action.
Like the matrices Xi , the Killing fields xi are linearly independent, forminga basis of the Lie algebra of Killing fields on S2 . That is, as we can
expand X = a X1 + b X2 + c X3 , we likewise get x = a x1 + b x2 + c x3 , where
the latter means point-wise addition of vectors in each Tx M at each x ∈ M.

such an observer field u is that we can compute it such that it minimizes
the observed time derivative (D/Dt)(v − u), given an arbitrary input
vector field v, whose time-dependent behavior we want to “follow.”
The second important property of an observer field u is that if u isindeed a Killing vector field, it in fact describes a single observer, as
defined above. Typically, we want to optimize the vector field u suchthat it describes observers that are as similar as possible. We do this byquantifying how much (or how “fast,” in terms of spatial distance) the
observers differ from one global isometry, by computing the Killingenergy of the vector field u over M. This is the integralZ

Eu dA.

(3)

M

The point-wise energy term Eu is given by Eq. 29. Eq. 3 is a scalarmeasure for how much u differs from an infinitesimal isometry on M.
See Hadwiger et al. [30] for a more detailed rationale in Euclidean
space. In this work, we now generalize this concept to curved spaces.
4.3 Feature-Relative Visualization
Fig. 1 and Fig. 6 illustrate how the computation of an observer field on
the surface of the Earth enables focusing a visualization on a featuresuch as a hurricane. In this example, the input field v(x,t) is given bya time-dependent global wind data set from the European Copernicus
project. The top row of Fig. 6 shows path lines following the motion of
the hurricane over time. In this frame of reference, it is hard to comparefeatures of the hurricane “in place,” due to its motion. In contrast,
the bottom row of Fig. 6 shows a visualization that is relative to an
automatically computed observer field u(x,t), following the motionof the hurricane. Since now a “feature-centric” observer is used, thehurricane has become steady. (Now the Earth is moving, relative to thenew observer, underneath the hurricane.) Path lines seeded at the same
positions as before now clearly highlight the vortex of the hurricane.
5 O BJECTIVITY IN C URVED S PACES
Given a particular observer transformation t 7→ g(t) ∈ G, it is natural todefine that a tensor field is objective if it simply undergoes the sametransformation as the observer. Because we model all possible observertransformations as paths through a continuous symmetry group G, i.e.,a Lie group, we can now rigorously formulate a generalized concept ofobjectivity for arbitrary manifolds M in a fully intrinsic manner, using
the concept of group action on M with a Lie group element g ∈ G.
We first summarize the definition of Truesdell and Noll [65], whichis standard in continuum mechanics [34, 48], but is restricted to Euclidean space. We then introduce our completely general definition that
includes curved spaces, but is equivalent in Euclidean space. App. D
derives the definition of Truesdell and Noll more generally, illuminating
the correspondence between their approach and ours in more detail.
5.1 Objectivity in Euclidean space
Truesdell and Noll [65, p.41] define objectivity with respect to a changeof reference frame in R3 via the distance-preserving transformationx∗ = c(t) + Q(t) x,t ∗ = t − a.

(4)To obtain a generalized definition of objectivity, a crucial propertyof the diffeomorphism φg is that it enables us to use the corresponding
differential, or pushforward. See Fig. 7. The pushforward is a map
dφg : T M → T M,
(a)where each (dφg )x at a point x ∈ M is a linear map(b)(dφg )x : Tx M → Tφg (x) M,
v 7→ (dφg )x (v).

(c)

(d)

Fig. 6. Feature-centric visualization of hurricane Isabel. (a,b) Pathlines in the original flow field. (c,d) Observer-relative pathlines; thehurricane appears steady. (a,c) First time step. (b,d) Last time step.
From (c) to (d), the Earth has moved underneath the steady hurricane.

Q(t) is a proper orthogonal tensor (a rotation), c(t) is a point (position
vector), and a ∈ R. This transformation assumes absolute time. It isthus sufficient to consider a = 0, disregarding time shifts, giving t ∗ = t.
With respect to this transformation, a scalar field is objective if it isunchanged; a vector field v is objective if it transforms according tov∗ = Q(t) v; a second-order tensor field S, as a linear transformation ofvectors, is objective if it transforms as S∗ = Q(t) S Q(t)T [65, p.42].
This entire definition depends on the domain being Euclidean: points
are position vectors; the difference between two points is a vector; alltangent spaces are copies of R3 with trivial parallel transport. Thisdefinition is therefore not valid for non-Euclidean (curved) manifolds.
5.2 Generalization of Objectivity
To generalize objectivity, we define this concept as a general notionof tensor fields being invariant with respect to a continuous symmetrygroup G, which is a Lie group. (Symmetry refers to a notion of being
the same.) For example, if the group G is chosen as the isometry groupof a (Riemannian) manifold, two tensor fields are “the same” if they
are isometric. Two fields being symmetries of each other then meansthat there exists a group element g ∈ G, such that the transformation
rules given below hold. Then, given any time-dependent observertransformation t 7→ g(t) ∈ G, a given tensor field is objective if, for each
fixed t, it simply follows the corresponding transformation g := g(t).
5.2.1 Symmetry groups and group actions
Our notion of symmetry corresponds to the transformation behavior
under a group action Φ, with a given Lie group element g ∈ G, where G
is the chosen symmetry group. An action Φ, specifically a smooth left
action, of a Lie group G on a manifold M, is a smooth map [33, p.209]
Φ : G × M → M,(g, x) 7→ Φ(g, x),(5)such that for every g ∈ G, the map
φg : M → M, with φg (x) := Φ(g, x), is a diffeomorphism.

(7)

(6)

Here, we focus on the general use of group actions Φ in our context,
and defer details to later sections. For now, it is sufficient to understandthat the diffeomorphisms φg will correspond to the flows of specificvector fields on M. These vector fields are generated by the action of
the Lie algebra g of the Lie group G on M. See App. J for details.
For example, if G is the group of all diffeomorphisms of M, these
vector fields are all possible (smooth) vector fields on M. The importantcase for our framework is choosing the group G as the isometry groupof M. The corresponding vector fields are then the Killing vector fields
on M, whose flows correspond to the isometries of M. See Sec. 7.

(8)

The notation (·)x means that the quantity in parentheses is located atx ∈ M, and Tx M denotes the tangent space at x. We can simply imagine
that the diffeomorphism φg transforms curves on M, and the differential
dφg transforms their tangent vectors accordingly. See also App. U.
In components, the map (dφg )x at any x ∈ M can be given by thecorresponding n × n matrix. See Fig. 7 for the case of a sphere (n = 2).
Euclidean space. When φg is an isometry of M = R3 , the pushforward dφg is a globally constant proper orthogonal (rotation) tensor Q,
i.e., (dφg )x = Q, with the same Q at all x ∈ M. See O’Neill [49, p.107].
Curved spaces. In general, however, the linear map (dφg )x is different for different points x ∈ M. In components, each (dφg )x can still begiven by a matrix, but it will be a different matrix for each point x ∈ M.
5.2.2

Objective scalar fields

Being objective should mean invariant under transformation, which for
scalar fields is trivial. We therefore define that a scalar field f : M → R
on a manifold M is objective when, under any diffeomorphism φg ,given by the group action Φ of a symmetry group G, it transforms asf ∗ φg (x) = f x .

(9)Abbreviated, we could write f ∗ = f , but it is crucial to note that f ∗ isevaluated at the point φg (x), whereas f is evaluated at the point x.
5.2.3

Objective vector fields

We now define that an arbitrary vector field v on a manifold M isobjective (with respect to a given symmetry group G), if, under thecorresponding group action Φ with any g ∈ G, it transforms as(v∗ )φg (x) = (dφg )x v .

(10)

We emphasize that v∗ is an element of the tangent space Tφg (x) M,
whereas v is an element of Tx M. Likewise, it is important to note thatthe differential (dφg )x is a linear map defined on Tx M. We can say
Remark. A vector field is objective, if it is simply pushed forward by
any diffeomorphism φg , defined according to the group action Φ. Thisdefinition of objectivity is valid for any smooth manifold where a notionof (smooth) symmetry is defined by a (smooth) symmetry group G.
Abbreviated transformation rule. For brevity, we can define theaction Φ, with g ∈ G, on any vector field v on M, by the differential in
Eq. 8, and abbreviate the objectivity criterion of Eq. 10 simply asv∗ = g v.

(11)

However, it is crucial that the meaning of the transformation represented
by g ∈ G in this shorthand notation is given by Eq. 10. In general, gcannot be mapped to the same globally defined matrix, corresponding to
the pushforward dφg , even though this is possible in the Euclidean case.
Nevertheless, this abbreviated form makes it easy to see the analogy
with the definition of Truesdell and Noll. In Euclidean space, the twoare equivalent. See App. D for more details. Our definition, however,gives a well-defined notion of objectivity for arbitrary manifolds M.

5.2.4 Objective tensor fields
The definition above for objective vectors generalizes to arbitrary tensor
fields T, in the way that is derived in full in App. A. Here, we derive thecommon special case of a second-order tensor field S, interpreted
 as alinear map of vectors, i.e., S : T M → T M, v 7→ S(v). Such a 11 tensor
field is objective, if under the group action Φ with g it transforms as(S∗ )φg (x) (v∗ ) = (dφg )x S (dφg−1 )φg (x) (v∗ ) .
(12)
The pushforward dφg−1 of the diffeomorphism φg−1 maps v∗ fromTφg (x) M to the corresponding v in Tx M according to Eq. 10 (inverted),
and then applies Eq. 10 again to S(v). φg−1 is the inverse of φg , meaning
φg−1 (φg (x)) = φg−1 g (x) = φe (x) = x. As in Eq. 11, we can abbreviate
S∗ = g S g−1 .

(13)

However, it is again crucial to note that the meaning of this rule is given
by Eq. 12. Fig. 7 illustrates how the pushforwards of φg and φg−1 map
tangent vectors between the tangent spaces Tφg (x) M and Tx M.
We note that to prove objectivity of the Killing operator K (Eq. 25),
which we use for computing observer fields, the generaldefinition of
objectivity from App. A is required, since Ku is a 02 tensor (App. F.1).
5.3 Objectivity for Intrinsic Isometries
To realize the notion of observer transformations as isometric transformations of the underlying space in the general mathematical framework
given above, we now consider the group G := Isom(M), i.e., the isometry group of the manifold M. The group action Φ is an isometry, withrespect to a given metric g, i.e., a smoothly varying inner product h·, ·i,on M, if for all g ∈ G, and for any two vector fields v, w on M, we have(dφg )x (v), (dφg )x (w) φ (x) = v, w x .
g(14)

As a direct consequence of Eq. 14 holding for the map φg , we get thatfor any vector field v that is objective according to Eq. 10, we havekv∗ kφg (x) = kvkx .

(15)

Infinitesimal isometries and Killing vector fields. Computationally,
instead of working with isometries directly, we work with their derivatives, which are infinitesimal isometries, corresponding to the Killingfields on M. See Sec. 7. To describe Killing fields efficiently, we need
the intrinsic velocity gradient tensor ∇v of a velocity field v on M.
6 G ENERAL V ELOCITY G RADIENT T ENSORS
For any smooth vector field v given on a manifold M, the velocity
gradient tensor ∇v is a second-order tensor field on M that allows
computing directional derivatives of the vector field v in any direction.
In Cartesian coordinates, ∇v can simply be defined by partial derivatives ∂ j vi . However, this definition cannot be used for intrinsic computations on curved manifolds, because ∂ j vi is in general not a tensor.
For a curved (or flat) manifold M, the tensor ∇v can, however, always
be computed intrinsically as the covariant derivative of v [22, p.241].

6.1

Intrinsic Covariant Derivatives

The covariant derivative ∇v of a vector field v on a manifold M is asecond-order tensor field on M, which, at any x ∈ M, maps an arbitrary
vector x ∈ Tx M to the derivative of the field v in the direction x, i.e.,(∇v)x : Tx M → Tx M,
x 7→ ∇x v = (∇v)x (x) .

This notation means that the covariant derivative of a vector field v,
evaluated in the direction x, is the vector ∇x v, which we get by computing ∇v (x). We can refer ∇v to a basis {ei ⊗ ω j }, with ⊗ the tensor
product, and write the components as n × n matrix ∇ j vi , with row andcolumn indices i, j, respectively, and n the dimension of M, i.e., fora curved surface n = 2. However, unless the manifold M is flat andCartesian or affine coordinates are used, ∇ j vi 6= ∂ j vi . The difference isa connection on M [22, p.242], corresponding to a notion of paralleltransport [19, p.52] of vectors. The connection can be given in components via Christoffel symbols Γi jk , and we get for ∇v in components,∇ j vi := ∂ j vi + Γi jk vk .

(17)

In flat space with Cartesian or affine coordinates, all Christoffel symbols
vanish [24, p.72] (Γi jk ≡ 0), but only in this special case is ∇ j vi = ∂ j vi .
The Γi jk determine the covariant derivatives of the basis vector fields{ei }, at each point x ∈ M, referred back to the same basis {ei } at x, i.e.,∇e j ek = Γi jk ei , and therefore Γi jk = ω i ∇e j ek .

(18)

We note that although Γi jk comprises n3 components, it is not a (thirdorder) tensor [24, p.82]. In general, a connection needs to be chosen [19,p.50]. We employ the Levi-Civita connection,3 which can be derived
intrinsically from the metric of M [22, p.229]. See App. O, Eq. O.10.
Even more easily, if we know an immersion of M into an ambient space
Rm , as for a curved surface immersed in R3 , we can obtain the Γi jk asΓi jk = ω̃ i ∂ j ẽk .

(19)

The tilde symbols denote ambient space (co)vectors, e.g., for M withn = 2 and ambient m = 3, {∂ j ẽk } comprises four 3D vectors, and {ω̃ i }comprises two 3D covectors. It is crucial that the {ω̃ i } are chosen suchthat they correspond to orthogonal projection from Rm onto the tangent
plane of the surface. See App. Q for a complete calculation for a simpleatlas of charts for the sphere. The Levi-Civita connection given by
this choice of Γi jk corresponds to the metric on M that is induced by
the immersion of M into the ambient space Rm [19, p.39]. Crucially,
however, the resulting Christoffel symbols are nevertheless intrinsic.
7

I SOMETRIES AND K ILLING F IELDS IN C URVED S PACES

Our framework requires an efficient mathematical description of the
intrinsic isometries of the manifold M. A powerful way to do this is
to work with the derivatives of diffeomorphisms that are isometries,
i.e., to consider infinitesimal isometries. These correspond to specificvector fields on the manifold M, which are called Killing vector fields.
A vector field u is a Killing vector field on M, if the flow φt (App. T)that it generates preserves the metric g on M. The most general way ofstating this mathematically is to require that the Lie derivative (App. G)of the metric g, with respect to u, must vanish everywhere on M. I.e,
Lu g = 0.

Fig. 7. The pushforward of a diffeomorphism φg is a linear map between
the tangent spaces Tx M and Tφg (x) M. That is, the pushforward (dφg )xmaps a tangent vector v ∈ Tx M to the vector (dφg )x (v) ∈ Tφg (x) M. The
pushforward (dφg−1 )φg (x) maps a tangent vector in the opposite direction.
In components, for a two-manifold M, each pushforward is a 2 × 2 matrix.

(16)

(20)

This means that for all x ∈ M, we have (Lu g)x = 0. See [46, p.250].
The set of all possible Killing vector fields on a given (Riemannian)manifold M constitutes a Lie algebra, which is, in particular, a vector
space. This enables us to consider the dimensionality of the set of
3 This is the unique metric-compatible connection, corresponding to a given(Riemannian) metric g on the manifold M, that has zero torsion [22, p.242].

(a)

(b)

(c)

(d)

(e)

Fig. 8. Rotating four centers relative to three different observers. Path lines seeded close to the four centers are shown in space-time as in Fig. 3.
(a) Observer rotating counter-clockwise around the center point between the four vortices. (b) The observer rotation in (a) is automatically removed
relative to an optimized observer field u, which corresponds to clockwise rotation: the four centers appear steady. (c) The field in (a) relative to an
observer rotating around the vertical axis. (d) The observer field u (path lines) optimized for the field in (c). (e) The field in (c) relative to u from (d).

Killing fields on M, and to write any Killing field as a linear combination of basis Killing fields. Fig. 5 shows an example on the sphere.
Apps. L and M describe the isometries of the two-sphere in detail, withfull calculations. App. K gives more details on Killing vector fields on
curved surfaces in general. However, for computing observer fields via
optimization (Sec. 9), as the most important part of our framework, we
directly evaluate Killing’s equation on the unknown vector field u.
7.1 Killing’s Equation
A vector field u = ui ei is a Killing field, if an equation known as
Killing’s equation holds at all x ∈ M [22, p.529]. However, in its
general form this equation is given for covectors, with components ui ,and not for vectors, with components ui . The ui can be computed asui = gi j u j [24, p.88], using a metric g, referred to components gi j .
Killing’s equation is then given as∇ j ui + ∇i u j = 0.

(21)This is the same as ∇u being anti-symmetric, which is the same ash∇u(x), xi = 0 for all x. However, in general it is not correct to simply
check whether a component matrix of ∇u is anti-symmetric, because∇u is a tensor of mixed type (contravariant and covariant indices).
Therefore, Eq. 21 is given using the covariant derivative of a covector
field [24, p.108]. In components, referred to a basis {ω i ⊗ ω j }, this is∇ j ui := ∂ j ui − Γki j uk .

(22)

We can introduce the vector components ui directly, and write Eq. 21 as∇ j gik uk + ∇i g jk uk = 0.
(23)
Applying the product rule for covariant derivatives [24, p.119], we
get ∇ j gik uk = (∇ j gik )uk + gik ∇ j uk = gik ∇ j uk , because the covariant
derivative of the metric vanishes.4 Therefore, Eq. 21 is equivalent togik ∇ j uk + g jk ∇i uk = 0.

(24)

In matrix notation, we can therefore rewrite Killing’s equation as thevanishing of the differential “Killing operator” K applied to the field u,T
Ku := g ∇u + g ∇u = 0.
(25)
The tensor Ku is still a 02 tensor, just like (∇ j ui + ∇i u j ). However,Eq. 25 implies that now we can simply check any matrix representationof g ∇u for anti-symmetry, to determine whether u is a Killing field.
7.2 Killing Energy
If the expression Ku is not exactly zero, in order to quantify how close
to being Killing an arbitrary vector field u is, we use the tensor normof the tensor Ku to define a scalar Killing energy (see also [9, Def. 5]).
As for vectors, the tensor norm of a higher-order tensor T is defined as
kTk2g := hT, Tig .

(26)

However, the tensor inner product required here, between two covariant
second-order tensors T, S, is given by the double contraction [45, p.82]
hT, Sig := gik g jl Tkl Si j = T i j Si j .
4 This (∇g = 0) is the definition of a metric connection [22, p.242].

(27)

This is similar to the Frobenius inner product tr(TT S) of matrices,which determines the Frobenius norm tr(TT T)1/2 . The crucial difference here is that because pairs of covariant and contravariant indices
are contracted, Eq. 27 is an invariant tensor expression, giving the same
result in all coordinate systems. hT, Tig = T i j Ti j thus is the square of
the norm. In matrix notation, we can write the inner product above ashT, Sig = tr g−1 TT g−1 S .

(28)

We define the Killing energy Eu, of the field u, at a point x ∈ M, as2

Eu := Ku g = Ku, Ku g .

(29)Sec. 9.1 defines the Killing energy of M as the integral of Eu over M.
Eq. 29 is then the density of the Killing energy of M per unit area.
8

O BSERVED T IME D ERIVATIVES IN C URVED S PACESThe observed time derivative defined in Eq. 1 measures the differential change of a tensor field with respect to the motion of a given
observer field u. We note that D/Dt is, in fact, the time-dependent
Lie derivative [45, p.95] with respect to the flow of the field u. The
crucial difference here is the semantic meaning of an observer field u,
as an intrinsic part of our operator. Our most important use of D/Dt isapplying it to the relative velocity field (v − u), which we minimize tocompute the unknown observer field u (Sec. 9). We get (Lu u = 0),D∂v ∂u(v − u) =−+ Lu v,
Dt
∂t
∂t
∂v ∂u
=
−
+ ∇v (u) − ∇u (v) .
∂t
∂t

(30)

Lie derivatives in curved spaces. A basic property of Lie derivatives
is that they are independent of the metric g defined on the manifold M,
and likewise independent of the connection on M [45, p.96]. That is,
Lu v = ∇v (u) − ∇u (v) ,= ∂ j vi u j − ∂ j ui v j ei .

(31)

The expansion in components in the second row above comprises only
partial derivatives. The full derivation is given in App. H. Thus, theobserved time derivative can be derived solely from partial derivatives,
even when the manifold M is curved, including the Lu term in Eq. 30.
9C OMPUTING O BSERVER F IELDS

We compute a global observer field u on any two-manifold M withmultiple charts using one global linear least-squares solve. The use of asingle global sparse linear system is enabled by “baking” all transition
maps (their Jacobian matrices) between charts into the system matrix.
9.1

Objective Function

Our optimization is conceptually very similar to the approach of Hadwiger et al. [30] for Euclidean space. However, all differential operators
have to be substituted by their generalized version for arbitrary manifolds, i.e., covariant derivatives, Lie derivatives, and tensor norms.

We formulate a global optimization problem to find an observer
field u, from the space of all possible fields in some function space Vof vector fields (e.g., Cn fields) on the manifold M, as the minimizer of1u∈V 2

Z tmaxZ

E +λ D+µ R

min

tmin




u, x,t dA dt.

(32)

M

Here, dA is an infinitesimal area element of M at a point x ∈ M, and theconstants λ , µ ∈ R are relative weights. The individual terms are
E (u, x,t) := Ku (x,t)

2g,D(v − u) (x,t) ,Dtg

R (u, x,t) := (v − u) (x,t)

2g.

σ (i, j) := n · (i − 1) + j, with i, j ∈ [1, n].

(39)

We can use this to rearrange the tensor inner product, to obtainT̄t := Tkl ,with t = σ (k, l),S̄s := Si j ,with s = σ (i, j),T̄ s = ḡst T̄t ,and therefore,

(40)

(33)

T̄ s S̄s = ḡst T̄t S̄s = gik g jl Tkl Si j = T i j Si j .

(34)Thus, the tensor inner product hT, Sig (Eq. 27) can simply be computed
as the contraction T̄ s S̄s . Using this approach, the tensor norm of theKilling energy in the least-squares problem given above is computed byinserting the components ḡst (Eq. 38) of ḡ−1 into the matrix Ḡ above.

2

D (u, x,t) :=The indices i, j, k, l are given by the inverse of the indexing function(35)
R

The integral of E over the domain gives the Killing energy M Eu dA of
the observer field u on M. || · ||g denotes the tensor and vector norms,
respectively, depending on the argument, with respect to the metric g.
Discretization. To solve Eq. 32, we discretize the manifold M witha triangle mesh. Vector components are stored at the vertices, and welabel the corresponding points on the manifold M with x1 , x2 , . . . , xN .
The Killing energy term E (Eq. 33) is computed by discretizing theKilling operator K (Eq. 25), to compute the covariant derivative termsof ∇u in Eq. 25. The weights for all differential operators to evaluate
Ku are stored in a matrix K. We compute partial derivatives via triangle
vertex 1-ring neighborhood stencils, computing weighted averages offunction values at the mesh vertices. The weights are independent of
the actual function values, and can therefore by pre-computed for a
fixed mesh geometry and be used for any function on M. See App. R.
The observed time derivative term D (Eq. 34) is computed frompartial derivatives alone (Eqs. 30 and 31). We compute time derivatives via central differences (in t) on each vertex. The weights for all
differential operators to compute D are stored in a matrix D. The regularization term R (Eq. 35) is represented by a matrix N, which is simply
an identity matrix times −1, to compute −u. All terms involving the(known) input field v are pre-computed and stored in the corresponding
locations of the matrices K, D, and the right-hand side b (see below).
9.2 Global Least-Squares Solution
To solve Eq. 32, we are looking for the minimizer of the least-squares
problem minu 21 ||Au − b||2g , with matrix A and vector b as given below.
The least-squares solution, weighted with the metric g, is the solutionof the normal equations of the weighted least-squares problem givenby the square, positive-definite system AT W A u = AT W b, with
 


 
0
K
Ḡ
0
0
A :=D, W := M · 0 λ · G
0 , b := − ∂∂tv  . (36)
N
0
0
µ ·G
v
The matrix M is diagonal, with weights giving the area on M assigned
to each vertex element xi , discretizing integration over area elements dA.
The other blocks of the matrix W encode g on M, given by the blocks
 −1

g
ḡx1
0
...
0
0 ...
0 
x1
−1

 0
0 
ḡx2
...
0 
 0 gx2 . . .
Ḡ:= .
..
.. .. , G:=
..
..
..
 ..
. (37)
 ..

.
.
.
.
.
.
.
0
0 . . . gxN00. . . ḡ−1xNHere, gxi refers to the n × n component matrix gi j of the metric g at thevertex xi ∈ M. For efficient computation of the tensor norm kTkg , weuse the following approach to convert the tensor inner product in Eq. 27
to a simple contraction of two first-order tensors. We introduce newn̄-dimensional first-order tensors T̄, S̄, and a symmetric n̄-dimensional
second-order tensor ḡ−1 , where n̄ := n2 . That is, n̄ = 4 for n = 2. We
denote the n̄ × n̄ components of ḡ−1 by ḡst , computed from gi j asḡst := gik g jl ,s,t ∈ [1, n̄], (i, j) = σ −1 (s), (k, l) = σ −1 (t). (38)

9.3 Multi-Chart Optimization
At each vertex xi ∈ M, there is only one vector for the field u. We
store vector components with respect to one selected chart for each xi .
However, finite differences to estimate derivatives (App. R) have to becomputed in the same chart, even when adjacent vertices are in different
charts. We enable this by “baking” all required chart transition maps,
i.e., the n × n Jacobian matrices transforming from the coordinate basisof one chart to that of another chart, into the matrices K and D. The
elements of these matrices are inserted such that all finite differences
are computed using vector components referencing the same chart.
More details for implementing this approach are explained in App. S.
9.4 Objectivity of Intrinsic Observer Field Computation
Remark. The computation of the observer field u as the minimizerof Eq. 32 is objective with respect to any observer transformation
t 7→ g(t) ∈ G, with G = Isom(M) the isometry group of the manifold M.
We give proofs in App. F. Except point (4) below, the reasons for thisare quite intuitive and can be summarized as follows: (1) The norms
used in Eq. 32 are invariant under isometries, because by definitionany
R isometry does not change them (Eq. 14). (2) The Killing energy
M Eu dA is time-independent, and measures how much the field u differs from an infinitesimal isometry. Since observer transformations areexact isometries of M, they do not change the Killing energy. (3) A timedependent observer transformation t 7→ g(t), which with G = Isom(M)is a time-dependent isometry t 7→ φg(t) of M, changes any velocity fieldon M by point-wise addition of the Killing field corresponding to thederivative X(t) of the isometry. However, the Killing field that is added
due to the observer transformation cancels out in the relative velocity
field (v − u). (4) The objectivity of the observed time derivative (Eq. 1)
is far from obvious. However, Marsden and Hughes [45, p.101] provethat the time-dependent Lie derivative of an objective tensor field withrespect to the flow of any vector field is objective. See App. F.3.
That is, the norm of all terms in Eq. 32 does not change under theisometry t 7→ φg(t) of any observer transformation t 7→ g(t), and because
the minimizer of Eq. 32 is unique, as in [30], its computation is thus
objective with respect to the group action Φ of the isometry group of M.
10 R ESULTS
We have implemented a general framework for flow visualization on
curved surfaces, given as triangle meshes, in C++, and OpenGL for
rendering. For the sphere, we use six very simple charts, as describedin detail in App. Q. We emphasize that the only 3D computations that
we use are for initially obtaining the metric of the sphere as a 2D tensor
field in each of the charts, and for visualization purposes in OpenGL.
All our vector fields are given via 2D components, and all computations
required for solving Eq. 32 are fully intrinsic in 2D. We solve theleast-squares system given in Sec. 9.2 by using the conjugate gradientsolvers provided by the C++ Eigen library, obtaining the minimizer u
in multiple iterations. We have also tried direct methods using matrix
factorization of the symmetric positive definite matrix AT W A. Due tothe size of our grids, conjugate gradients obtain the minimum faster.

Table 1. Data sets with grid/mesh resolution.

Table 2. Computation times for the observer field u.

Data set# Vert.

# Tri.

Time

Time Steps

Data setλµ# CG Iter.

Comp.-time

SYNTHETIC VORTEX

10,242
10,24240,962
163,842
32,40062,412

20,48020,48081,920327,680
64,796
124,820

0 . . . 1s
0 . . . 2π s
0 . . . π2 s
8 days8 days8 days

100
100
300
32
32
32

SYNTHETIC VORTEX

0.1
1
103
1013
1013

0
0
0
0
0

1,000
10,0005, 00010,00010,000

6s
78s
12min
3min
6min

SYNTHETIC FOUR CENTERSJUPITER VORTEX STREETEARTH FLOW

·EARTH FLOW ( SUBDOMAIN )·EARTH FLOW ( ADAPT. RES .)Data sets. Table 1 gives an overview of the data sets that we haveused in this paper and the video. We ran experiments with synthetic data
sets and with simulated data on the sphere. The optimizations are global.
For the earth flow data set we have, however, experimented with two
approaches for region of interest optimization to reduce computation
times. Table 2 gives computation times for the conjugate gradientsolves to optimize observer fields u. All computations were run on twoIntel Xeon E5-2600 processors with 16 cores running at 2.3GHz.
Visualization. We have adapted the observer-relative visualization
methods described by Hadwiger et al. [30] from Euclidean space tocurved spaces. Since we need coordinate charts to compute the optimization, we can also use them as texture parametrization for LIC
computations in texture space on the GPU. Corresponding to a given
observer field u, we compute LIC images of observed stream lines [30].
Likewise, we also compute the geometry of integral curves for observed
stream lines as well as for observed path lines. To visualize the timeparameter of path lines, we increase the radial distance from the sphereaccording to the parameter t, scaled by a user-specified factor. This is
helpful to see the time evolution of path lines. Additionally, we draw
path line shadows on the sphere, to ensure that their position in spacecan be correctly interpreted. Moreover, we visualize time-dependentscalar fields relative to the observer field via color coding, and in thesame way show texture patterns, e.g., Earth or Jupiter textures. We
can dynamically switch the observer for any flow field, and directly
visualize the corresponding observed fields. Figs. 1, 6, and 9 show
observer-relative scalar fields (planet textures, vorticity magnitude), as
well as observed LIC. Figs. 3, 6, and 8 show observed path lines.
10.1 Synthetic Data Sets
Synthetic Vortex. We transform a steady input field of a simple vortex
into different time-dependent fields by rotating the input around thevertical axis of a sphere at different angular velocities (Fig. 3). The
latter correspond to different scalings of the Killing field x2 in Fig. 5.
From each transformed, time-dependent input field, we computethe corresponding observer field u via optimization. In all cases, theminimizer u of Eq. 32 is a steady Killing field. The observer fields thenenable recovering the original steady vortex in Fig. 3 (d). Relative toits observer field, each unsteady field becomes the original steady field.
Synthetic Four Centers. Fig. 8(a) shows four vortices that rotate
around each other clockwise, around a common point in the centerof the sphere. We show that the observer field u that is computed
automatically via optimization recovers a steady field for the four
vortices, as depicted in Fig. 8(b). In Fig. 8(c), we then observe the input
field in (a) relative to a leftward rotation around the vertical axis of the

(a)

(b)

(c)

Fig. 9. Planet-scale flow patterns observed using LIC. Vorticity iscolor-coded on a red-blue CCW/CW scale. LIC is used to visualize
instantaneous velocities. (a) A vortex street at the equator. LIC cannot
show the vortices. (b) The globally optimized observer field that transforms the flow field (a) into a near-steady flow field (c), which shows thevortices in place while the planet is moving underneath instead. Note how
LIC of the instantaneous velocity field now clearly shows the vortices.

SYNTHETIC FOUR CENTERSJUPITER VORTEX STREETEARTH FLOW ( SUBDOMAIN )EARTH FLOW ( ADAPT. RES .)sphere. Fig. 8(e) again depicts the field that is recovered automatically,
this time relative to an observer field u, shown in Fig. 8(d), that rotatesclockwise as well as rightward around the vertical axis of the sphere.
10.2 Global Planet-Scale Optimization
Jupiter Vortex Street. Fig. 9 depicts a 2D von Kármán vortex street
centered around the equator of Jupiter. We have generated this data
set for testing purposes by mapping a flat vortex street to a sphere bydirectly mapping a Cartesian x-y grid to a latitude-longitude grid.
Earth Flow. For real planet-scale experiments we use open source
flow data from the EU Copernicus project. We use global wind data toshow that large phenomena like hurricane Isabel (2003) are automatically tracked by the solution of our optimization problem, although no
explicit feature tracking is performed. Figs. 1 and 6 show the input flow
data and the observed flow fields. We are able to compute an observed
flow field that is nearly steady. To get an undistorted visualization thatis centered on the hurricane, we can pick one time-dependent Killingfield as the observer that travels with the hurricane. The resulting
visualization applies the inverse observer motion to the Earth texture.
10.3 Region of Interest Optimization
Our approach optimizes for the best observer field u with a global
optimization. However, often an optimization of the whole domain at
the same resolution level is not necessary, since the user might only
be interested in a smaller subdomain. We have tested two methods forthis case. We either select a subdomain and run the global optimization
only in this subdomain, or we adaptively remesh to obtain different,
adaptive mesh resolutions for different regions on the sphere, and thenrun the optimization over the whole domain. The computation times
can be greatly reduced from a full-resolution grid either way (Table 2).
11 D ISCUSSION AND C ONCLUSIONS
We have presented the first objective, optimization-based framework forcomputing time-dependent reference frames in curved spaces with thegoal of turning unsteady input flow fields into (nearly) steady fields relative to a global observer field. A major motivation for this approach is toenable feature-relative visualization and analysis of the time evolutionof features, such as scalar and vector field attributes of hurricane simulations, “in place.” A powerful example of this is shown in Fig. 1 (topleft), where the path of hurricane Isabel is simply a path line of theautomatically computed observer field u. We have seeded this path line
manually, but this could also be done algorithmically. Moreover, pathlines of u anywhere inside the core region of the hurricane follow thegeneral time-dependent diffeomorphism that constitutes the warping
domain of the hurricane. We believe that our mathematical frameworkhas great potential in applications such as studying features of interestin climatology and oceanography in a reference frame-independentway. Our method is objective, meaning that it is invariant with respect
to any isometric observer transformation of the input flow field. This
enables it to automatically follow the motion of features of interest,and makes it independent of super-imposed motions, including thoseof moving measurement devices capturing real-world geophysical flow
phenomena. Finally, we believe that the fully intrinsic approach that
we have taken, and the methodology from Riemannian geometry that
we have used in our framework, could provide significant advantages
and insight to other areas and applications of flow visualization as well.
ACKNOWLEDGMENTS
We thank Anna Frühstück for the illustrations and for help with the figures and the video.
Hurricane Isabel data courtesy of EU Copernicus project, path from National Hurricane
Center/Wikipedia. This work was supported by King Abdullah University of Scienceand Technology (KAUST), and the KAUST Office of Sponsored Research (OSR) award
OSR-2015-CCF-2533-01. This research used resources of the Core Labs of KAUST.

R EFERENCES[1] M. Alexa, D. Cohen-Or, and D. Levin. As-rigid-as-possible shape interpolation. In Proceedings of SIGGRAPH 2000, pp. 157–164, 2000.
[2] R. Aris. Vectors, Tensors and the Basic Equations of Fluid Mechanics.
Dover Publications, Inc., 1990.
[3] V. I. Arnold. Mathematical Methods of Classical Mechanics. SpringerVerlag, 2nd ed., 1989.
[4] O. Azencot, M. Ben-Chen, F. Chazal, and M. Ovsjanikov. An operator
approach to tangent vector field processing. Computer Graphics Forum,32(5):73–82, 2013.
[5] O. Azencot, M. Ovsjanikov, F. Chazal, and M. Ben-Chen. Discrete derivatives of vector fields on surfaces–an operator approach. ACM Transactions
on Graphics, 34(3):Article No. 29, 2015.
[6] S. Bachthaler, M. Strengert, D. Weiskopf, and T. Ertl. Parallel texturebased vector field visualization on curved surfaces using GPU cluster computers. In Proceedings of Eurographics Symposium on Parallel Graphics
and Visualization (EGPGV), pp. 75–82, 2006.
[7] I. Baeza Rojo and T. Günther. Vector field topology of time-dependent
flows in a steady reference frame. IEEE Transactions on Visualization andComputer Graphics, 26(1):280–290, 2020.
[8] H. Battke, D. Stalling, and H.-C. Hege. Fast line integral convolution for
arbitrary surfaces in 3d. In Visualization and Mathematics: Experiments,
Simulations, and Environments, pp. 181–195. Springer, 1997.
[9] M. Ben-Chen, A. Butscher, J. Solomon, and L. Guibas. On discrete Killing
vector fields and patterns on surfaces. In Proceedings of Eurographics
Symposium on Geometry Processing, pp. 1701–1711, 2010.
[10] H. Bhatia, G. Norgard, V. Pascucci, and P.-T. Bremer. The HelmholtzHodge decomposition–a survey. IEEE Transactions on Visualization andComputer Graphics, 19(8):1386–1404, 2013.
[11] H. Bhatia, V. Pascucci, and P.-T. Bremer. The natural Helmholtz-Hodge
decomposition for open-boundary flow analysis. IEEE Transactions on
Visualization and Computer Graphics, 20(11):1566–1578, 2014.
[12] H. Bhatia, V. Pascucci, R. M. Kirby, and P.-T. Bremer. Extracting features from time-dependent vector fields using internal reference frames.
Computer Graphics Forum, 33(3):21–30, 2014.
[13] R. Bujack, M. Hlawitschka, and K. I. Joy. Topology-inspired Galilean
invariant vector field analysis. In Proceedings of IEEE Pacific Visualization
2016, pp. 72–79, 2016.
[14] B. Cabral and L. C. Leedom. Imaging vector fields using line integral
convolution. In Proceedings of SIGGRAPH ’93, pp. 263–270, 1993.
[15] G. Chen, V. Kwatra, L.-Y. Wei, C. D. Hansen, and E. Zhang. Design of
2d time-varying vector fields. IEEE Transactions on Visualization andComputer Graphics, 18(10):1717–1730, 2012.
[16] G. Chen, K. Mischaikow, R. S. Laramee, P. Pilarczyk, and E. Zhang. Vector
field editing and periodic orbit extraction using morse decomposition.
IEEE Transactions on Visualization and Computer Graphics, 13(4):769–
785, 2007.
[17] F. de Goes, M. Desbrun, and Y. Tong. Vector Field Processing on Triangle
Meshes. In SIGGRAPH Asia 2015 Courses, pp. 17:1–17:48, 2015.
[18] F. de Goes, B. Liu, M. Budninskiy, Y. Tong, and M. Desbrun. Discrete
2-tensor fields on triangulations. Computer Graphics Forum, 33(5):13–24,
2014.
[19] M. P. do Carmo. Riemannian Geometry. Birkhäuser, 1992.
[20] M. Edmunds, R. S. Laramee, G. Chen, N. Max, E. Zhang, and C. Ware.
Surface-based flow visualization. Computers & Graphics, 36(8):974–990,
2012.
[21] L. K. Forssell and S. D. Cohen. Using line integral convolution forflow visualization: curvilinear grids, variable-speed animation, and unsteady flows. IEEE Transactions on Visualization and Computer Graphics,
1(2):133–141, 1995.
[22] T. Frankel. The Geometry of Physics: An Introduction. Cambridge
University Press, 3rd ed., 2011.
[23] A. Frishman and C. Herbert. Turbulence statistics in a two-dimensionalvortex condensate. Phys. Rev. Lett., 120(20):204505, 2018.
[24] P. Grinfeld. Introduction to Tensor Analysis and the Calculus of Moving
Surfaces. Springer-Verlag, 2013.
[25] T. Günther, M. Gross, and H. Theisel. Generic objective vortices for flow
visualization. ACM Transactions on Graphics, 36(4):Article No. 141,2017.
[26] T. Günther, M. Schulze, and H. Theisel. Rotation invariant vortices forflow visualization. IEEE Transactions on Visualization and Computer
Graphics, 22(1):817–826, 2016.

[27] T. Günther and H. Theisel. Hyper-objective vortices. IEEE Transactions
on Visualization and Computer Graphics, 25(1):1–1, 2018.
[28] T. Günther and H. Theisel. The State of the Art in Vortex Extraction.
Computer Graphics Forum, 37(6):149–173, 2018.
[29] A. Hadjighasem and G. Haller. Geodesic transport barriers in Jupiter’s
atmosphere: A video-based analysis. SIAM Review, 58(1):69–89, 2016.
[30] M. Hadwiger, M. Mlejnek, T. Theußl, and P. Rautek. Time-dependent
flow seen through approximate observer Killing fields. IEEE Transactions
on Visualization and Computer Graphics, 25(1):1257–1266, 2019.
[31] G. Haller. An objective definition of a vortex. Journal of Fluid Mechanics,
525:1–26, 2005.
[32] G. Haller, A. Hadjighasem, M. Farazmand, and F. Huhn. Defining coherent vortices objectively from the vorticity. Journal of Fluid Mechanics,
795:136–173, 2016.
[33] D. D. Holm, T. Schmah, and C. Stoica. Geometric Mechanics and Symmetry: From Finite to Infinite Dimensions. Oxford University Press, 2009.
[34] G. A. Holzapfel. Nonlinear Solid Mechanics: A Continuum Approach forEngineering. Wiley, 2000.
[35] B. L. Hua and P. Klein. An exact criterion for the stirring properties ofnearly two-dimensional turbulence. Physica D: Nonlinear Phenomena,
113(1):98–110, 1998.
[36] J. Huang, W. Pei, C. Wen, G. Chen, W. Chen, and H. Bao. Output-coherent
image-space LIC for surface flow visualization. In Proceedings of IEEE
Pacific Visualization, pp. 137–144, 2012.
[37] B. Jobard, G. Erlebacher, and M. Y. Hussaini. Lagrangian-Eulerian advection of noise and dye textures for unsteady flow visualization. IEEE
Transactions on Visualization and Computer Graphics, 8(3):211–222,
2002.
[38] T. Langer, A. Belyaev, and H.-P. Seidel. Spherical barycentric coordinates.
In Proceedings of Eurographics Symposium on Geometry Processing, pp.
81–88, 2006.
[39] R. S. Laramee, B. Jobard, and H. Hauser. Image space based visualizationof unsteady flow on surfaces. In Proceedings of IEEE Visualization 2003,
pp. 123–130, 2003.
[40] J. Laurie, G. Boffetta, G. Falkovich, I. Kolokolov, and V. Lebedev. Universal profile of the vortex condensate in two-dimensional turbulence. Phys.
Rev. Lett., 113(25):254503, 2014.
[41] J. M. Lee. Introduction to Smooth Manifolds. Springer-Verlag, 2nd ed.,2012.
[42] J. M. Lee. Introduction to Riemannian Manifolds. Springer-Verlag, 2nded., 2018.
[43] G.-S. Li, X. Tricoche, and C. Hansen. GPUFLIC: Interactive and accurate
dense visualization of unsteady flows. In Proceedings of Eurovis 2006, pp.
29–34, 2006.
[44] G.-S. Li, X. Tricoche, D. Weiskopf, and C. Hansen. Flow charts: Visualization of vector fields on arbitrary surfaces. IEEE Transactions on
Visualization and Computer Graphics, 14(5):1067–1080, 2008.
[45] J. E. Marsden and T. J. Hughes. Mathematical Foundations of Elasticity.
Dover Publications, Inc., 1994.
[46] A. McInerney. First Steps in Differential Geometry. Springer-Verlag,
2013.
[47] S. B. Myers. Isometries of 2-dimensional riemannian manifolds into
themselves. Proc. National Academy of Sciences of the USA, 22(5):297–
300, 1936.
[48] R. W. Ogden. Non-Linear Elastic Deformations. Dover Publications, Inc.,
1997.
[49] B. O’Neill. Elementary Differential Geometry. Academic Press, revised
2nd ed., 2006.
[50] J. Palacios and E. Zhang. Rotational symmetry field design on surfaces.
ACM Transactions on Graphics, 26(3):Article No. 55, 2007.
[51] J. Pedlosky. Geophysical Fluid Dynamics. Springer-Verlag, 2nd ed., 1987.
[52] Z. Peng and R. S. Laramee. Higher dimensional vector field visualization:
A survey. In Proceedings of EG UK Theory and Practice of Computer
Graphics, pp. 149–163, 2009.
[53] P. Petersen. Riemannian Geometry. Springer-Verlag, 3rd ed., 2016.
[54] P. G. Saffman. Vortex Dynamics. Cambridge University Press, 1995.
[55] H.-W. Shen and D. L. Kao. UFLIC: A line integral convolution algorithm
for visualizing unsteady flows. In Proceedings of IEEE Visualization ’97,
pp. 317–322, 1997.
[56] J. Sokol. Flattened fluids help scientists understand oceans andatmospheres, July 2018.
[Online:www.wired.com/story/
flattened-fluids-help-scientists-understand-oceans-\
and-atmospheres; posted 27-June-2018].

[57] J. Solomon. Discrete killing fields for pattern synthesis and symmetry
detection. Undergraduate honors thesis, Stanford University, Palo Alto,CA, 2010.
[58] M. Spivak. A Comprehensive Introduction to Differential Geometry (5volumes). Publish or Perish Press, 3rd ed., 1999.
[59] D. Stalling. LIC on surfaces. Technical report, Konrad Zuse Institute
Berlin (ZIB), Berlin, Germany, 1997.
[60] D. Stalling and H.-C. Hege. Fast and resolution independent line integral
convolution. In Proceedings of SIGGRAPH ’95, pp. 249–256, 1995.
[61] D. Sujudi and R. Haimes. Identification of swirling flow in 3-d vector fields.
In Proceedings of the 12th Computational Fluid Dynamics Conference,
pp. 792–799, 1995.
[62] C. Teitzel, R. Grosso, and T. Ertl. Line integral convolution on triangulated
surfaces. In Proceedings of WSCG, pp. 572–581, 1997.
[63] K. S. Thorne and R. D. Blandford. Modern Classical Physics: Optics,
Fluids, Plasmas, Elasticity, Relativity, and Statistical Physics. Princeton
University Press, 2017.
[64] Y. Tong, S. Lombeyda, A. N. Hirani, and M. Desbrun. Discrete multiscale
vector field decomposition. In Proceedings of SIGGRAPH, pp. 445–452,
2003.
[65] C. Truesdell and W. Noll. The Nonlinear Field Theories of Mechanics.
Springer-Verlag, 1965.
[66] L. W. Tu. Differential Geometry: Connections, Curvature, and Characteristic Classes. Springer-Verlag, 2017.
[67] J. J. van Wijk. Image based flow visualization. ACM Transactions on
Graphics, 21(3):745–754, 2002.
[68] J. J. van Wijk. Image based flow visualization for curved surfaces. In
Proceedings of IEEE Visualization 2003, pp. 123–130, 2003.
[69] A. Vaxman, M. Campen, O. Diamanti, D. Panozzo, D. Bommes, K. Hildebrandt, and M. Ben-Chen. Directional Field Synthesis, Design, and Processing. Computer Graphics Forum, 35(2):545–572, 2016.
[70] T. Weinkauf, J. Sahner, H. Theisel, and H.-C. Hege. Cores of swirling
particle motion in unsteady flows. IEEE Transactions on Visualization
and Computer Graphics, 13(6):1759–1766, 2007.
[71] E. Zhang, K. Mischaikow, and G. Turk. Vector field design on surfaces.
ACM Transactions on Graphics, 25(4):1294–1326, 2006.