1 Introduction
The visualization of vector field topology is a problem that arises naturally when studying the qualitative structure of flows that are tangential to some surface. As usual, we use the term surface for a real, smooth 2manifold (equipped with an atlas consisting of charts), see for example [O’N83] for an introduction to Riemannian Geometry. Having its roots in the theory of dynamical systems, the topological skeleton of a Hamiltonian flow on a surface with isolated critical points consists of these critical points and trajectories (streamlines) of the vector field that lie at the boundary of a hyperbolic sector and connect two of the critical points. Helman and Hesselink [HH90] introduced the concept of the topology of a planar vector field to the visualization community and proposed the following construction scheme: (1) critical points are located, (2) classified, and then (3) trajectories along hyperbolic sectors are traced and connected to their originating and terminating critical points or boundary points. Step (2)—the classification of a critical point—is usually based on the Jacobian of the vector field, see [Per91]
. Trajectories of step (3) are typically constructed by solving an ordinary differential equation for particle tracing.
Although a vast body of previous work in the field of flow visualization focuses on the problem of how to extend the method of Helman and Hesselink to vector fields on arbitrary surfaces as well as the second and third step of the above algorithm, not much attention has been paid to the first step. In this paper, we specifically address the identification and classification of critical points in parameter space.
As efficient computerbased visualization algorithms usually work with discrete parametrized versions of the surfaces involved — examples of popular discretization schemes are for example triangulated or quadrangulated versions of the surface — we will in this paper not focus on the wellresearched field of how to parametrize a given surface (see [FH05] for a recent survey) but assume that a surface always comes equipped with a globally continuous discrete parametrization that allows a cellwise (local) barycentric or bilinear interpolation scheme of a vector field tangential to the surface in parameter space.
While this task is rather easy for linear vector fields, the problem setting becomes more interesting for bilinearly interpolated fields. Bilinear interpolation is ubiquitous in scientific visualization because it is popular for widely used uniform or curvilinear grids representing planar or curved surfaces. Since bilinear interpolation is not linear, it can lead to higherorder critical points, which are neglected by oftenused linearization approaches.
(a)  (b)  (c) 
In this paper, we introduce a new method that locates and classifies all critical points within piecewise linearly (barycentrically) or bilinearly interpolated twodimensional grids. Our method determines the index of a critical point without the need to evaluate the Jacobian of the vector field in the critical point in order to determine its Poincaré index. For the case of bilinearly interpolated vector fields, our method is able to detect higherorder critical points and the presence of two firstorder critical points within one cell, which, to our knowledge, has not been achieved with the common methods [HH91] for the bilinear interpolation scheme before. Figure 1 shows a corresponding example and illustrates our classification method. Our approach is based on the idea that investigating the qualitative behavior of 0level sets of the components of the interpolated vector field provides information needed to compute the Poincaré index of a critical point. All qualitatively different possibilities of this behavior and the types of critical points yielded by each possibility are completely classified using the computational group theory tool GAP (Groups, Algorithms, and Programming) [GAP06]. See the enumeration of cases for marching cubes and generic substitope algorithms by Banks et al. [BLS04] for a previous example of an application of computational group theory in the field of scientific visualization. Furthermore, we discuss a cellbased topology simplification method as well as the question of numerical stability. Our approach results in an efficient, accurate, and robust cellbased algorithm for detecting all critical points of barycentrically or bilinearly interpolated 2D vector fields.
The paper is organized as follows. First, we will give a short review of the visualization literature dealing with vector field topology. Then, the theoretical foundations of vector field topology, namely the theory of the qualitative behavior of secondorder dynamical systems along with such fundamental notions as those of critical points, separatrices, and the Poincaré index of a critical point are reviewed. Following this, we will present our new approach—first the general framework will be discussed and then applied to two cases: barycentrically and bilinearly interpolated vector fields. Then, we will deal with open issues such as critical points on the boundary of cells and numeric stability followed by a more detailed description of the cellbased algorithm. We will conclude giving results and a short review of our method.
2 Previous Work
Topologybased methods for planar vector fields were first proposed to the visualization community by Helman and Hesselink [HH90], employing methods from the theory of dynamical systems [ALGM73] to planar linear (or linearized) vector fields in order to visualize flow characteristics. Their work has triggered a large body of further research on topologybased flow visualization, an overview of which is given by Post et al. [PVH03] and Scheuermann and Tricoche [ST05].
The case of a planar linear vector field is relatively simple to deal with, but it has a couple of drawbacks, most notably that only firstorder critical points can be detected. Much effort has been put into methods to overcome this drawback imposed by the interpolation scheme and to addresses the issue of detecting and processing higherorder critical points of interpolated vector fields, also of vector fields on arbitrary surfaces [LCJK09]. Such higherorder critical points can be found, for example, by using piecewise linear interpolation schemes in combination with a clustering of firstorder critical points according to Tricoche et al. [TSH00]. Another example is the work by Theisel [The02], who proposes a method for designing piecewise linear planar vector fields of arbitrary topology. Nonlinear interpolation schemes have been investigated by Scheuermann et al. [SHK97] and Zhang et al. [ZMT06]. Scheuermann et al. [SKMR98] propose a way to approximate higherorder critical points using Clifford algebra. Li et al. [LVRL06] use interpolation schemes based on a polar coordinate representation to detect vector field singularities on a surface.
In contrast to the global variational approach taken in [PP03] in which the authors construct a discrete Hodge decomposition in order to obtain location and type of critical points of vector fields on polyhedral surfaces, our approach is local and celloriented. It may thus be easier to implement when a discrete parametrization of the surface is already given and can be used in conjunction with other local, gridbased methods.
3 Theoretical Background
The methods used for extracting vector field topology are founded upon the theory of the qualitative behavior of dynamical systems. Most of the brief review of the essential theory in this section is taken from the books [DLA06, Per91, ALGM73].
3.1 Dynamical Systems
Definition 3.1
A dynamical system on , an open subset of , is a map , where with , , that satisfies

,

.
Dynamical systems are closely related to autonomous systems of differential equations. On the one hand, let be open and let be Lipschitz continuous on . Then the initial value problem of the autonomous system of differential equations
(3.1) 
with for any has a unique solution defined for all by virtue of the Picard–Lindelöf theorem. For each initial value this induces a mapping referred to as a trajectory of the system (3.1). The mapping then lies in and thus is a dynamical system in the sense of Definition 3.1. It is called the dynamical system induced by the system of differential equations (3.1). On the other hand, if is a dynamical system on , then
defines a vector field on and for each , is the solution to the initial value problem with of (3.1).
3.2 Critical Points
An important concept in the field of dynamical systems is the notion of a critical point:
Definition 3.2
An equilibrium or critical point of a dynamical system is a point where . If the dynamical system is induced by a system of differential equations (3.1), then a critical point of is a point where . If the Jacobian
has only eigenvalues with nonvanishing real part,
is called a hyperbolic critical point. If , then is called nondegenerate or firstorder critical point. Otherwise it is called degenerate or higherorder critical point.A system of differential equations can be approximated by its linearization around a critical point without changing its qualitative behavior if is a hyperbolic critical point of that system (HartmanGrobman theorem [Per91]). For planar linear systems, only certain types of critical points can occur and these can be classified in terms of the eigenvalues of the Jacobian as shown in Fig. 2 (here only nondegenerate cases with are considered).
3.3 Poincaré Index of a Critical Point
In order to classify critical points of vector fields one can use the notion of the Poincaré index (or index for short) of a critical point.
Definition 3.3
Let be a vector field on some open . If is an isolated critical point of and a Jordan curve such that is the only critical point of in its interior, then the Poincaré index of (or index for short) is
with .
It can be shown that isolated firstorder critical points (i.e. isolated critical points for which the Jacobian of the vector field in the critical point has no eigenvalue of 0) have a Poincaré index of 1 and that a saddle has a Poincaré index of 1, whereas nonsaddles have a Poincaré index of 1 (see Fig. 2).
3.4 Topological Equivalence and Sectors
Let us now establish the fundamental notion of topological equivalence of vector fields:
Definition 3.4
Suppose that and with open sets . The two autonomous systems of differential equations and and thier induced vector fields are said to be topologically equivalent if there exists an orientation preserving homeomorphism that maps trajectories of the first system onto trajectories of the seconds system.
Markus [Mar54] showed that for planar systems of differential equations the condition of being topologically equivalent is equivalent to the systems having the same separatrix configurations, where a separatrix of a system (3.1) is a trajectory of (3.1) which is either a critical point, a limit cycle, or a trajectory lying on the boundary of a hyperbolic sector as defined below. This justifies the use of the term vector field topology for the topological skeleton of a vector field consisting of separatrices of that field.
The notion of sectors was first introduced by Poincaré [Poi93] to investigate higherorder critical points of planar systems, and later extended by Bendixon [Ben01] and Andronov [ALGM73]. The idea is that one can describe the qualitative behavior of a planar vector field in a suitable neighborhood of an isolated critical point of in terms of connected regions, so called sectors, which form a partition of . Within each sector the trajectories of exhibit a behavior that is characteristic for this type of sector. It can be shown that there exist three topologically different types of sectors:
Definition 3.5
A sector of a critical point can be classified as a hyperbolic, parabolic, or elliptic sector according to its topological structure as shown in Fig. 3.
Tricoche et al. [TSH00] also used the idea of sectors to model higherorder critical points with a piecewise linear interpolation scheme.
3.5 Bifurcation Theory
Bifurcation theory is based on the notion of structural stability of a vector field due to Andronov and Pontryagin [AP37]: if the qualitative behavior of a dynamical system (3.1) does not change for small perturbations of the vector field , then that vector field is called structurally stable. If is not structurally stable, the topological skeleton of the vector field changes even under small perturbations of the vector field and is said to be structurally unstable or to lie within the bifurcation set.
The perturbation of the vector field is usually modeled by an additional parameter :
(3.2) 
A value of for which the system (3.2) lies in the bifurcation set is called a bifurcation point of (3.2) and is then called bifurcation value of (3.2).
Bifurcation theory has been studied extensively, see for example the book by Guckenheimer and Holmes [GH90]. It also explains the splitting of higherorder critical points into several nearby firstorder critical points: it can be shown that, if a vector field has an isolated critical point of higher order, there exists a perturbation of such that splits into several isolated firstorder critical points nearby.
4 Classifying Critical Points Without Using the Jacobian
In this section, we show how the Poincaré index of an isolated firstorder critical point can be computed by evaluating the sign configuration of the vector field’s components on a finite set of sample points in a neighborhood of the critical point, reminiscent of the marchingcubes classification applied to isosurfaces in scalar fields.
4.1 Setting
From now on let be a vector field defined on an open such that is Lipschitz continuous on and only has nondegenerate firstorder isolated critical points.
4.2 level sets, Areas of Characteristic Behavior
We introduce the notion of areas of characteristic behavior that will enable us to calculate the index of a critical point.
Lemma 4.1
Let be a vector field like in 4.1
and an isolated firstorder critical point of
, i.e.
. Then, lies in
the intersection of the 0level sets and of and
. Furthermore, there exists an such that and
are infinitesimally straight lines (i.e. lie infinitesimally
close to straight lines) in an neighborhood of .
Proof.
By definition, a critical point of has to lie in the intersection set of the 0level sets of the components of . Since is linearizable around , Taylor’s theorem leads to with in the ball . Since has full rank, the 0level sets of are straight lines intersecting in .
Definition 4.2
Let be a vector field like in 4.1 and an isolated firstorder critical point of . Then for an the 0level sets of and of partition the ball of into four disjoint open subsets called areas of characteristic behavior or areas for short.
In each area, the signs of and do not change, i.e. for arbitrary and the following holds:
The set of areas around a critical point can be seen as an ordered, cyclic sequence of areas according to the order in which they intersect with the boundary of , as illustrated in Fig. 4. We consider clockwise traversal direction unless otherwise noted.
Since each of the two vector field components can either be positive or negative in one area, there exist four different types of characteristic areas as shown in Fig. 5(a). In the following, areas are written as ordered 2tuples over the set or equivalently as elements of the set , where , , , . Area sequences can then be written as ordered 4tuples over the set of areas, i.e. as ordered 4tuples over the set . Two areas that lie next to each other are called adjacent areas. As and are adjacent, the indices of the areas in the area sequence are cyclic, i.e. the area sequence is glued together at and .
Remark 4.3
For two adjacent areas of an area sequence with , , either
holds, i.e. exactly one component flips its sign for two adjacent areas but not both.
4.3 Area Sequence and Types of Critical Points
Definition 4.4
Let be an isolated firstorder critical point of defined like
above with an area sequence
. Then for a pair
of adjacent characteristic areas of the area
sequence of (see Fig. 5(a)), a
turning in the characteristic behavior of the vector field
can be defined. This turning can either be a clockwise
turning or a counterclockwise turning as defined in
Figs. 5(b) and 5(c),
respectively.
Since Figs. 5(b) and 5(c)
contain all 8 possible configurations of pairs of characteristic
areas, only a clockwise or counterclockwise turning behavior can
occur and the list in Figs. 5(b) and
5(c) is complete.
Remark 4.5
One adjacent pair of characteristic areas already determines the whole area sequence in terms of its turning behavior. This is the case because the vector field components whose signs flip are alternating in the area sequence. Thus, an area sequence can either show a clockwise turning behavior or a counterclockwise turning behavior, but not a mixture of these and we will also speak of a clockwise or counterclockwise turning behavior of the area sequence as a whole—a clockwise or counterclockwise area sequence for short.
We can use the turning behavior of area sequences around a critical point to determine its Poincaré index:
Theorem 4.6
Let be a vector field like in 4.1 and let be an isolated critical point of of first order, such that has full rank. If the area sequence of is counterclockwise, then the Poincaré index of is and is a topological saddle of . If the area sequence of is clockwise, then and is a nonsaddle firstorder critical point of .
Proof.
The area sequence can also be interpreted as a sequence of four qualitative samples of the vector field values lying on a piecewise linear closed curve that contains (qualitative in the sense that just the signs of the vector field components are sampled). We will now prove that summing up the angle change of the vector field between these discrete samples and performing this for all four samples yields the same result as evaluating the continuous integral of the change of angle of over a continuous Jordan curve around , only containing one critical point . Therefore, the Poincaré index of can be computed by just identifying the sequence of characteristic areas around .
Let be the linearized system of defined on an ball around . Since for two adjacent areas only the sign of one component of changes, the vector at the first sampling point and the vector at the second sampling point lie in two different but adjacent quadrants around (see Fig. 6). As the value of the linear approximation changes linearly on between two sample points and is continuous, the turning behavior is uniquely determined by the two sample points; and since the two sample points lie in adjacent quadrants defined by the coordinate axes, the vector rotates about an angle with when walking on from one sample point to the next. Thus, the whole area sequence (for four pairs of areas) yields the total change in angle , with . As the area sequence wraps around the critical point, the first and the last vector of our sample sequence are the same and thus the change in angle of the vector field has to be with , see Fig. 5(b) and Fig. 5(c). As , for some also for , i.e. the vector field is not constant in a neighborhood of and one has . It is because . Thus, only is possible and the critical point of the linearized field is of Poincaré index .
Since for linear systems all critical points of index have been classified (see Sections 3.2 and 3.3 as well as Fig. 2), a critical point can either be a saddle (for ), if the area sequence shows a counterclockwise turning behavior, or a nonsaddle (i.e. a source or sink for ), if the area sequence shows a clockwise turning behavior. According to the HartmanGrobman theorem, the Poincaré indices of firstorder isolated critical points are invariant under linearization and has a saddle in if and only if has a topological saddle in , which completes the proof of the theorem.
4.4 Invariant Operations on the Area Sequence
There are 8 possible area sequences as there are four different areas and each area sequence can be clockwise or counterclockwise. Since these sequences yield only two types of critical points distinguishable by their Poincaré index, it is desirable to build equivalence classes of area sequences that yield critical points of the same Poincaré index. These equivalence classes can be directly constructed using elementary group theory. We refer to the book [Hup67] for a comprehensive overview of the theory of finite groups.
First it is obvious but nonetheless interesting to observe that the turning behavior of the area sequence is invariant under rotation and also invariant under a simultaneous sign flip of both vector field components. Figure 5 illustrates this: for example, the first and the last area pair of Fig. 5(b) and Fig. 5(c) can be related through a sign flip, likewise the second and third area pair.
Rotations and sign flips can be modeled as group operations. Given a group and a set , a (left) group action of on is a mapping denoted such that for all (here denotes the neutral element in ) and for all and all . The orbit of an element under a group action of on is the set . Orbits of a group action on a set define equivalence classes on , and the set of all orbits is a partition of .
The rotation of the area sequence is identical to the operation of the cyclic group on the indices of the area sequence. The sign flip can be modeled as the operation of a group isomorphic to the symmetric group on the signs of the components. From now on, the first group is referred to as the shape group and the second group as the (sign) flip group . The direct product of and is called the coloring group .
Now we define the group action of on the set of area sequences. Let be the projection of onto and the projection of onto . Further, let be an element of the coloring group with its projections and onto the shape group and the flip group, respectively. Then, acts on an area sequence as defined below:
where is a permutation of the indices and a selfinverse permutation on the set of areas, . can be interpreted as an element that flips the signs of all components of the vector field: areas of type are mapped to areas of type and vice versa; an area of type 2 is interchanged with an area of type 3.
The orbits of this group action on the set of all possible area sequences yield equivalence classes of area sequences such that each equivalence class contains all area sequences that can be mapped onto each other by rotations and sign flips.
4.5 Interpolated Vector Fields
Typical vector field data is given on a grid. Therefore, the vector field needs to be interpolated to obtain values at nongrid points. Local (i.e. cellwise) interpolation schemes are often chosen for the sake of simplicity and speed. In the following two sections, we will have a closer look at two interpolation schemes: barycentric interpolation on triangles and bilinear interpolation on rectangles.
Both interpolation schemes share some important properties. On the one hand, inside cells, both interpolation schemes are of class and thus linearizable. On the other hand, they are defined locally or in a piecewise way: a different interpolant is used for each cell. The interpolation is only of class across cell boundaries and it is linear and continuous along cell boundaries. Basic concepts and tools are developed for the simpler case of the barycentric interpolant. Then, these tools are adapted and extended for bilinear interpolation.
Note that the case of linear interpolation on triangular meshes can equally efficiently be solved by calculating the rotation along each triangle directly. As there can either be no or exactly one critical point of Poincaré index inside each cell, a cell with zero rotation has no critical point on its inside, whereas a nonzero rotation implies that there is a critical point inside the cell and already determines its Poincaré index. None the less we will describe the barycentric setting in the following as the notions introduced there will also be used in the more subtle case of the bilinear interpolant.
5 Barycentric Case
A simplex on vertices offers a natural way for the linear interpolation via barycentric coordinates. For the following, we will again restrict the dimension to , where a simplex is identical with a triangle.
5.1 Cells
From now, a cell denotes a triangle with vertices , , and with one real 2D vector (from the vector field) attached to each vertex such that a cell becomes a 3tuple . Most of the time, we are only interested in the vector field values and not the position of the vertices, just writing . The vector field inside a cell will be written , when is the set of convex combinations of the vectors . Further, let . We restrict the vector field values on the vertices to be nonzero such that the interpolated vector field has isolated singularities only.
5.2 Sector Sequences and level sets of the Interpolant
The vector field is interpolated componentwise inside the cell using barycentric coordinates. The interpolation is linear, and the interpolant is continuous in the cell and on its boundaries.
Definition 5.1
Given the cell , let us look only at one component of , say (). If two adjacent vertices (i.e. two endpoints of an edge of the cell) have values of () such that one value is smaller than and the other one bigger than , then by virtue of continuity and linearity of the interpolant on cell edges there exists exactly one point on the edge with an value (value) of . Such an edge is called an active edge for (). If the two values are both smaller or bigger than , there exists no such point and the edge is called an inactive edge for ().
Remark 5.2
We can now observe the following:
1) For one cell, there exist at most two active edges for each component and, thus, at most two points on the cell boundary with value for each component if all vertex values of each component differ from .
2) Since a critical point of the interpolated vector field can only occur in a crossing point of the 0level sets of the components of the vector field, a cell with an interior critical point requires two active edges per component. Furthermore, there can at most be one critical point inside a cell because the 0level sets of the interpolant are straight lines.
3) Since there can only be one critical point inside a cell, the sector sequence around a critical point can be found by looking for intersections of the 0level sets of the components with the cell boundary. Each active edge yields exactly one 0value point on a boundary edge, the position of which can be obtained through linear interpolation between adjacent vertices.
If one collects the sequence of 0value points of the components while traversing the boundary of the cell (as shown in Fig. 7), one can not only construct the area sequence but also determine whether the 0level sets cross and thus whether there is a critical point of the interpolated vector field. Let denote a 0value point of the first component and let denote a 0value point of the second component of the vector field. Then, for the sequences , , , or , the cell does not contain a critical point, whereas for the sequences and , there is a critical point. These intersection sequences are all possible sequences for the barycentric interpolant because there can be at most two 0value points for each component on the boundary of the cell.
The sequence of the crossings of the 0level sets of the components together with the information of how the signs of the components change defines the sequence of characteristic areas around a critical point and thus its Poincaré index as we have shown in Section 4.3. We will make use of this in the following section.
5.3 An EdgeColoring Problem
Before we can classify critical points, a suitable data structure for describing the area sequence around a critical point is needed. Just looking at the sign configuration of the components in the vertices (which can be seen as a coloring of the vertices) does not allow us to distinguish the types of critical points because the area sequence is not fully determined just by the signs of the components in the vertices (see Fig. 8).
(a)  (b) 
We use an edgebased data structure to uniquely describe the sequence of 0values of the vector field components on the boundary of the cell and the area sequence: each cell edge is given two slots that can be filled with 0value points of the components (see Fig. 7). The slots represent the order of the two components’ zero values while walking around the boundary of a cell. The slots also indicate whether a component changes from positive to negative or from negative to positive values as it passes through a 0value point. There is no more than one 0value for each component along a boundary edge of a cell because of the linearity of the interpolation along cell edges. Hence each edge can be given two slots and 13 types of boundary edges can occur, as listed in Table 1. The notation denotes that changes from a positive value to a negative value, changes from a negative value to a positive value, and the sign change of occurs before the sign change of as one traverses the edge on the cell boundary in clockwise direction. If a component is not listed, it indicates that this component has no sign change on that edge and thus no 0value.
Edge type  0values of components 

no 0value, neither for nor for  
For a triangle cell, each of the three edges can now be of types –, which represent a coloring of the triangle edges with 13 colors. Therefore, such an edgecoloring can be seen as an ordered 3tuple over the set of the 13 colors. Note that not every possible 3tuple over the set of colors is a valid coloring as there exist the following types of invalid colorings, i.e. colorings that are not (or not uniquely) mappable to a sign configuration of the vector field at the cell vertices:
1) Colorings that have an invalid number of active edges:
For each component, the cell has to have at least two active edges and the number of active edges has to be even.
Example: the coloring
would have one active edge for
both and and would thus
be an invalid coloring.
2) Colorings that are not mappable to a sign configuration of the vector
field in the cell vertices at all.
Example: a component, say , cannot
change from to twice in a row. Thus, for example,
a coloring would be invalid.
Each valid coloring of a triangle can by construction be uniquely mapped to a sign configuration of the values of the vector field components on the vertices (see Fig. 11). However, the edgecoloring of a triangle carries more information than just the sign configuration of the vector field components on the vertices (namely a unique description of the intersection topology of the 0level sets with the cell boundaries) such that several colorings of the triangle may be mapped to the same sign configuration.
5.4 Classification
As discussed in Section 4.4, there are certain operations that leave the clockwise and counterclockwise turning behavior of the area sequence invariant. These operations can be modeled as a group action of a certain group on the set of area sequences.
Since there exists a bijective map from the set of edgecolorings of a cell to the set of area sequences (each coloring describes exactly one area sequence), one can make use of the ideas presented in Section 4.4 to construct equivalence classes of cell colorings yielding the same area sequence. Then, all possible edge colorings of a cell can be constructed and classified by using a group action that builds equivalence classes of equivalent colorings. Colorings related by rotation or color flip thus yield the same type of critical point.
Theorem 5.3
Every configuration of a cell (as a coloring of the cell) that results in an intersection of the level sets of the two components inside of and thus contains a critical point is topologically equivalent to one of 8 representative configurations of which four contain a saddle firstorder critical point (Poincaré index ) and four contain a nonsaddle firstorder critical point (Poincaré index ). Representatives for the 8 orbits are listed in Table 2 and visualized (Online Resource 1). Colorings that are not included in this list do not contain a critical point.
Class  Cell coloring  Index 

1  
2  
3  
4  
5  
6  
7  
8 
Proof.
This theorem is proven using a GAP program, see (Online Resource 3).
The basic idea is as follows. Using a combinatoric description, all possible colorings of a cell are constructed: the set of all ordered 3tuples over the set of all possible colors. Then, a group action on that set is considered to build equivalence classes of colorings such that all pairs of elements of an equivalence class can be mapped onto each other by means of operations leaving the type of critical point invariant, as described in Section 4.4.
Let be the set of colors as described above and let the set of 3tuples over , be the set of colored cells or tiles. For the coloring group , which is the direct product of the shape group and the flip group (see Section 4.4), an action of on the set of all colored tiles can be defined. Here, the shape group is the cyclic group as the subgroup of rotations of the symmetry group of the combinatoric triangle (which is given by the symmetric group ). The color flip group is chosen to be , which is isomorphic to ; the numbers in the cycles correspond to the color indices of edge colors. The choice of generator for is obvious and unique since simultaneous flipping of signs of the two components results in an interchange of edges with color and , and , etc.
Using our GAP program, first all possible colorings of a cell are generated and for each orbit a representative is checked for the validity of the coloring. Elimination of invalidly colored cells can be done on a representative level: if the representative of an orbit is an invalid coloring, all other elements in the orbit are invalid because they are equivalent colorings. Thus, one can talk of invalidly colored orbits.
After discarding all invalidly colored orbits, for the remaining orbits the sequence of 0value points of the components on the cell boundary is extracted from the coloring. Since the 0level sets of the components are straight lines, the sequence of 0value points determine whether the 0level sets intersect within the cell (thus yielding a critical point) or not (no critical point). If the 0level sets intersect within the cell, the area sequence and its turning behavior are extracted to classify the type of critical point as saddle (Poincaré index 1) or nonsaddle (Poincaré index 1). As all critical points are of firstorder, this classification by the Poincaré index covers all possible types of critical points. Since all possible colorings have been considered the list of equivalence classes is complete. See (Online Resource 1) for a list of all equivalence classes, including sample visualizations.
6 Bilinear Case
In this section, the concepts developed in the previous section are transferred to the bilinear case. Most elements can be immediately adopted, but some caution has to be taken because the interpolant is no longer linear.
6.1 Interpolation Scheme and Cells
Bilinear interpolation is an extension of linear interpolation for interpolating functions of two variables. This interpolation scheme is simple, fast to implement, and widely used in many visualization algorithms working in a cellwise fashion on uniform, rectilinear, or curvilinear 2D grids.
A cell is represented by an ordered 4tuple on the points with one realvalued 2D vector (the vector field values) attached to each vertex such that a cell becomes a 4tuple . Mostly, we are only interested in the vector field values, just writing .
Any uniform, rectilinear, or convex curvilinear cell can be transformed to the unit square in terms of a diffeomorphism, yielding local Euclidean coordinates within the cell. Without loss of generality we therefore only consider the case of Euclidean coordinates in the following, i.e. the setting when the interpolation scheme is employed in parameter space. We use the following notation for bilinear interpolation:
Let , , be bilinear functions for the two vector components (corresponding to and ), defined by
where are local coordinates within one cell and ,…, are the values of the vector field components at the four vertices . Then is called the bilinear interpolant of .
Remark 6.1
The bilinear functions can be rewritten as
(6.1) 
with , ,
and .
The following properties of the interpolant will be of importance throughout the rest of the section:

The interpolation is linear in each dimension and continuous. Furthermore, the interpolation along edges of the cell is linear: analogously to the barycentric case in Section 5, one can define active and inactive cell edges. There exist four points on a cell boundary with value and thus most four active edges, if all vertex values differ from .

In the case , is a linear function. If , is nonlinear and when looking at the partial derivatives and , one sees that has an extremal point at . This extremal point is a saddle because the Jacobian is singular at and for all the Hessian of has a negative determinant (and thus its eigenvalues are of opposite sign; see Figs. 9 and 12 for examples).
6.2 Bilinear Interpolation of Scalar Fields
Let us look at the qualitative behavior of a bilinearly interpolated scalar fields before extending this to vector fields. We only consider one here, say , as a placeholder for a scalar field. The qualitative behavior of the interpolant within the cell can be described in terms of its vertex configuration:
Theorem 6.2
Given and a cell with vertex values , , and , for , four cases for the qualitative behavior of level sets of the interpolant within a cell can occur: Fig. 9 summarizes these cases that depend on the classification of vertex values being greater or less than . Here, active and inactive edges for a cell are defined in the same way as for the barycentric interpolant in Definition 5.1. Please note that we use the terminology saddle cell (see Fig. 9) to describe a configuration of vertex values , which is different from the classification of a critical point as a saddle point. In the following, we denote the first always by the compound term saddle cell.
All of the statements follow immediately from the continuity of the interpolant and its linearity along cell boundaries.
Remark 6.3
For an easier notation of the different cell types in terms of the definitions of Theorem 6.2, we write a cell as a 4tuple over the set where the th entry of that tuple is set to if and to otherwise. Table 3 summarizes this tuple notation, referred to as vertex sign configuration, as in the barycentric case.
Cell type  Vertex sign configuration 

inactive  
single active  
double active  
saddle cell 
Lemma 6.4 (Boundingbox lemma)
Let be a bilinear interpolant on a rectangular grid cell and let the cell be active. For each curve of the level set connecting two value points , on the boundary of the cell, the Cartesian product
defines a bounding box: the curve of the level set of is inside except at the value points at the border of the cell. See Fig. 9 for an illustration of the set for some sample configurations.
Proof.
We only give a sketch of the proof. The basic idea is the following.
Let the bounding box be given as the convex hull of the four points
Note that . Then, the bilinear interpolant restricted to the bounding box can be expressed in terms of another bilinear interpolant defined on and interpolating the values of the points . To leave , a 0isoline of has to cross the boundary of . Thus, a crossing 0isoline can only occur where the value of on the boundary of is 0. One can then show that apart from the two 0values in the corner points and of a , there exist no other boundary points of with 0value. Therefore, such a crossing cannot occur and the 0level set of completely lies inside and cannot leave .
Since the interpolation function along the boundary of is continuous and linear, it is sufficient to show that only two points in the set have 0value. It then follows immediately that these are the only points on the boundary with zero value. This can be shown easily for a representative of each of the cases listed in Theorem 6.2 and Fig. 9.
6.3 Analytical Computation of the Position of Critical Points
The position of critical points of the bilinear interpolant can be computed as intersection points of 0level sets of and by solving
(6.2) 
Since both and are linear or quadratic equations in , the system can either be combined into a single linear or quadratic equation in or . We will refer to this equation as the intersection equation. Since the linear case occurring here can be treated in exactly the same way as the barycentric case described in Section 5, this degenerate case will be excluded in the following and we will only speak of the (quadratic) intersection equation with discriminant . The equation can either have no real solution (), two different real solutions (), or a double real solution (), and the number of solutions of the intersection equation gives the number of critical points, see Fig. 10.
Since the sign of the discriminant determines whether the vector field has no, one, or two critical points, it can be interpreted as a bifurcation parameter of a dynamical system, where is the bifurcation value. Also can be interpreted as a measure for the stability of the critical point—critical points with are generally not stable, whereas the case with two () or no intersection points () can be more or less stable according to the magnitude of .
Given a rectangular cell with the vector field values
one can write the components of the bilinear interpolant
in normal form like done in (6.2) (see Remark 6.1).
The intersection points of the 0level sets of the components are obtained by solving (6.2). First, we determine which of the interpolants is linear in and (i.e. for which ); if both are linear, (6.2) becomes a linear system and can be solved separately yielding one or no solution (and not infinitely many solutions as the singularities are restricted to be isolated). If one of the interpolants is linear and the other nonlinear, the linear one can be solved for or and this can be plugged into the other to yield a quadratic equation in or . If both interpolants are nonlinear, one can be solved for or and plugged into the other.
All these equations have at most two real solutions yielding one coordinate of the 0value points of the interpolated vector field. The other coordinate can be obtained through the linear or quadratic equation between and from the first step. A solution within the cell is found by restricting the coordinates to be in . Note that the vector field values for each component are always restricted to be nonzero at cell vertices and that the two components are assumed not to be exactly the same.
6.4 Level Sets and Sector Sequences
Each solution of (6.2) lies in the intersection set of the 0level sets of and . The following cases can occur when looking at the intersections of level sets of and :
Theorem 6.5
Proof.
Clearly one has: and are disjoint which proves 1). The same holds for 2). Now it remains to show 2), i.e. that and cannot touch twice and that the intersection equation has a double solution iff and touch in one point.
Suppose that and touch in one point as shown in Fig. 10(b). A small perturbation of the vector field values now either results in a case where the hyperbola branches do not intersect (i.e. ) as shown in Fig. 10(a), or it yields a case where and intersect in two points, i.e. where the intersection equation has two different solutions as shown in Fig. 10(c). As each pair of branches can be perturbed separately a double touching case cannot exist as this could be perturbed to yield a triple intersection of the branches which is a contradiction to the maximum number of real solutions of the intersection equation. This proves 3).
6.5 Critical Points and Area Sequences
Let us now look at the types of critical points yielded for the different intersection cases for the 0level sets of the components as investigated in Theorem 6.5.
1) The case , where the 0level sets do not intersect, is trivial.
2) For the case , when the 0level sets of the components intersect twice, yielding two firstorder isolated critical points of the interpolated vector field of first order as none of the directional derivatives is in the intersection points. For each critical point separately, the same things hold as for critical points in the barycentric case; in particular we can look at area sequences, turning behavior of those, etc. for each critical point separately.
Remark 6.6
The area sequences of the two critical points always share three characteristic areas and these are traversed in opposite directions for each of the points (i.e. clockwise turning for an area pair becomes a counterclockwise turning and vice versa). Therefore, one of the two points has a Poincaré index of and the other one has a Poincaré index of . Another way to see this is the fact that to first observe that there are at most two critical points inside a cell and that the Poincaré index of a cell can either be , or . Now, by virtue of the summation theorem for the Poincaré index, two critical points inside a cell of the same index would force the cell to have an index of , a contradiction.
3) For the touching case () of the 0level sets, the tools developed in the previous sections cannot be applied: when linearizing the field around the touching point, the 0level sets of the two components are tangential to each other and the Jacobian has at least one zero eigenvalue. Thus, the critical point is not of firstorder and a welldefined area sequence does not exist for this case. However, a perturbation of the field yields two critical points, one of index and one of index . Then, Poincaré index theory states that the one critical point yielded by the touching of the 0level sets has an index equal to the sum of the indices of the two critical points it splits into, i.e. the secondorder critical point from a touching of the 0level sets has index 0.
6.6 Edge Coloring and Invariant Operations on the Colorings
As the bilinear interpolant is linear along the cell edges, one can employ the same ideas as in the barycentric case in order to determine existence and types of critical points inside a cell, again using an edgebased “slot” data structure. Using a data structure with two slots for each edge of a cell, one can see this as a cell coloring problem with a coloring of a cell represented as a 4tuple over the set of 13 colors like described in Section 5.3.
Investigating the intersection topology of the 0level sets becomes more complicated because the 0level sets are no longer straight lines. Using the same notation as for the barycentric case, a sequence of 0value points (where again stands for a 0value point of the first component of the vector field and for a 0value of the second component) not necessarily has no intersections. Here, it can either yield a case where the isolines have no point in common, touch in one point, or intersect twice as shown in Fig. 1.
The question is whether for a sequence of 0value points of the form or the subsequence or can yield a touching or intersection of the 0level sets or not. If not, the part or of the sequence is irrelevant for the intersection topology and can be ”collapsed”, e.g. a sequence for which cannot yield a touching of the isolines could be reduced to