An Implementation of a Convex Hull Algorithm

Version 1.0

Michael Müller Joachim Ziegler

MPI–I–94–105 February 1994 Acknowledgements

This work was supported in part by the ESPRIT Basic Research Actions Program of the EC under contract No. 7141(ALCOM II) and the BMFT (Förderungskennzeichen ITS 9103).

Abstract

We give an implementation of an incremental construction algorithm for convex hulls in IR^{d} using Literate Programming and LEDA in C++. We treat convex hulls in arbitrary dimensions without any non-degeneracy assumption. The main goal of this paper is to demonstrate the benefits of the literate programming approach. We find that the time we spent for the documentation parts is well invested. It leads to a much better understanding of the program and to much better code. Besides being easier to understand and thus being much easier to modify, it is first at all much more likely to be correct. In particular, a literate program takes much less time to debug. The difference between traditional straight forward programming and literate programming is somewhat like the difference between having the idea to a proof of some theorem in mind versus actually writing it down accurately (and thereby often recognizing that the proof is not as easy as one thought).

Keywords

Literate Programming, Convex Hull, WEB, LEDA, Object Oriented Programming, Computational Geometry, Incremental Construction Contents

1. INTRODUCTION 1

1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 6. The basic structure of the program . . . . . . . . . . . . . . . . . . . 4 11. The fundamental data structures . . . . . . . . . . . . . . . . . . . . 6 17. Computing plane equations . . . . . . . . . . . . . . . . . . . . . . . 12 24. The insert procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 31. Finding x-visible hull facets . . . . . . . . . . . . . . . . . . . . . . . 25 45. The Dimension Jump . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 51. Output Routines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 54. The main function . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 60. Solving a system of linear equations . . . . . . . . . . . . . . . . . . . 49 74. Dependency tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 76. Useful literature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 1. Introduction.

We give an implementation of an incremental construction algorithm for convex hulls in IR^{d} using Literate Programming (cf. [5]) and LEDA (cf. [6, 7]) in C++. The algorithm has been developed by Clarkson, Mehlhorn and Seidel (cf. [4]). In [2], a minor modification of this algorithm is described which maintains convex hulls in arbitrary dimensions without any non-degeneracy assumption.

2. Our main goal was to show how a complex algorithm can be implemented in a way such that everybody can easily understand the program. Therefore, we used literate programming. From LEDA (a Library of Efficient Data types and Algorithms) we took some useful and well known data structures. The reader not familiar with LEDA should not worry about lines of code like

list⟨vector⟩ L; because they all have their natural meaning: L is a list of vectors. All LEDAcommands are selfexplanatory.

We will first introduce the notation and describe the strategy of the algorithm. To do so, we will essentially cite parts of [4] and [2]. The citations appear in a smaller font and are terminated by a mention of the source (e.g., cited text

(cf. reference)).

3. The convex hull is constructed incrementally.

Let R = {x_{1}, …, x_{n}} be the multi–set of points whose convex hull has to be maintained and let π = x_{1} … x_{n} be the insertion order. Let π_{i} = x_{1} … x_{i}, R_{i} = {x_{1}, …, x_{i}} and let conv R_{i} be the convex hull of the points in R_{i}. Let d = dimR be the dimension of the convex hull of R and let DJ = {x_{j_{1}}, x_{j_{2}}, …, x_{j_{d+1}}} with 1 ≤ j_{1} ≤ … ≤ j_{d+1} ≤ n be the set of dimension jumps where x_{k} is called a dimension jump if dimR_{k−1} < dimR_{k}. Clearly, j_{1} = 1. In the incremental construction of conv R we maintain a triangulation Δ(π_{i}) of conv R_{i}: a simplicial complex whose union is conv R_{i} (a simplical complex is a collection of simplices such that the intersection of any two is a face of each^{1}). The vertices of the simplices in Δ(π_{i}) are points in R_{i}. The triangulation Δ(π_{i}) induces a triangulation CH(π_{i}) of the boundary of conv R_{i}: it consists of all facets of Δ(π_{i}) which are incident to only one simplex of Δ(π_{i}). If x ∈ aff R_{i} then a facet F of CH(π_{i}) is called visible from x or x–visible (we also say: x can see the facet) if x does not lie in the closed halfspace of aff R_{i} that is supported by F and contains conv R_{i}.

The triangulation Δ(π_{1}) consists of the single simplex {x_{1}}. For i ≥ 2, the triangulation Δ(π_{i}) is obtained from Δ(π_{i−1}) as follows. If x_{i} is a dimension jump, i.e., x_{i} ∉ aff R_{i−1}, then x_{i} is added to the vertex set of every simplex of Δ(π_{i−1}). If x_{i} is not a dimension jump then a simplex S(F ∪ {x_{i}}) = conv(F ∪ {x_{i}}) is added for every x_{i}–visible facet of CH(π_{i−1}). Figure 1 gives an example. For a simplex S let vert(S) denote the set of vertices that define this simplex. It is clear that Δ(π) contains a simplex whose vertex set is precisely the set of dimension jumps. We call this simplex the origin simplex of Δ(π). For every simplex S (besides the origin simplex) we call the vertex in vert(S) − DJ, that has been inserted last, the peak of S and the facet of S opposite to the peak the base facet of S. (cf. [2], p. 4)

1Note that the empty set is a facet of every simplex.

[figure: x_{11}; x_{9}; x_{3}; x_{6}; x_{7}; x_{2}; x_{10}; x_{1}; x_{8}; x_{5}; x_{4}]

Figure 1: A triangulation. The dimension jumps are the points x_{1}, x_{2}, and x_{5}.

4. It is convenient to extend Δ(π) to a triangulation \overline{Δ}(π) by also making the facets of CH(π) the base facet of some simplex: \overline{Δ}(π) is obtained from Δ(π) by adding a simplex S(F ∪\{\overline{O}\}) with base facet F and peak \overline{O} for every facet F of CH(π). Here \overline{O} is a fictitious point without geometric meaning. We propose to store the triangulation \overline{Δ}(π) as the set of its simplices together with some additional information: For each simplex S ∈ \overline{Δ}(π) we store its set of vertices, the equation of its base facet normalized such that the peak lies in the positive halfspace, and for each simplex S and vertex x ∈ vert S we store the other simplex ^{2}sharing the facet with vertex set vert(S) \ {x}. We also store a pointer to the origin simplex and a suitable representation of aff R, e.g., a maximal set of affinely independent points. The simplical complex \overline{Δ}(π_{1}) consists of two simplices: the bounded simplex S(\{x_{1}\}) and the unbounded simplex S(\{\overline{O}\}). (A simplex is called bounded if \overline{O} does not belong to its vertex set and it is called unbounded otherwise.) (cf. [2], p. 5)

\overline{O} is also called the anti-origin. Figure 2 shows the extended triangulation of the example of Figure 1. The points are numbered according to their insertion time. A base facet is indicated by an extra line. You can see that only the origin simplex has no base facet. All outer simplices have \overline{O} as peak. The point x_{12} sees the base facets conv(x_{2}, x_{7}) and conv(x_{7}, x_{9}). The simplex opposite to the vertex x_{1} with respect to the simplex conv(x_{1}, x_{3}, x_{5}) is the simplex conv(x_{3}, x_{5}, x_{8}). The vertex opposite to x_{1} in this simplex is the vertex x_{8}.

2They mean: a pointer to the other simplex

[figure: .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; x_{11}; x_{3}; x_{6}; x_{7}; x_{2}; x_{10}; x_{1}; x_{8}; x_{5}; x_{4}; x_{12}; x_{9}; \overline{O}; \overline{O}; \overline{O}; \overline{O}; \overline{O}; \overline{O}]

Figure 2: An extended triangulation

5. We give additional details of the insertion process. Consider the addition of the i–th point x = x_{i}, i ≥ 2. First decide whether x is a dimension jump (an O(d^{3}) test). If x is a dimension jump then add x_{i} to vert S for every simplex of \overline{Δ}(π_{i−1}) and add the simplex S(F ∪ \{\overline{O}\}) to \overline{Δ}(π_{i}) for every bounded simplex F of \overline{Δ}(π_{i−1}).

If x_{i} is not a dimension jump then we proceed as described in [4]. We first compute all x_{i}–visible facets F of CH(π_{i−1}) and then update the extended triangulation \overline{Δ} as follows: For each x_{i}–visible facet F of CH(π_{i−1}) (≡ x_{i}–visible base facet of an unbounded simplex in \overline{Δ}(π_{i−1})) we alter the simplex S(F ∪ \{\overline{O}\}) of \overline{Δ}(π_{i−1}) into S(F ∪ {x_{i}}). Moreover, for each new hull facet F ∈ CH(π_{i}) \ CH(π_{i−1}) we add the unbounded simplex S(F ∪\{\overline{O}\}). In other words, for each horizon ridge f of CH(π_{i−1}), i.e., ridge where exactly one of the incident facets is x_{i}–visible, we add the simplex S(f ∪ \{x_{i}, \overline{O}\}). The set of x_{i}–visible facets F of CH(π_{i−1}) can be found by visiting simplices according to the rule: Starting at the origin simplex visit any neighbor of a visited simplex that has an x_{i}–visible base facet. (cf. [2], pp. 5–6)

We call this search method the visibility search method.

Another search method is as follows. Walk through the simplices along a segment Ox_{i}^{− ; − ; − ; →} from a point O in the origin simplex to x_{i}. This leads us to the simplex containing x_{i}. If this simplex is unbounded, we have found an x_{i}-visible hull facet from which we can reach all other x_{i}-visible hull facets. This method is called the segment walking method (cf. [4], p. 11).

6. The basic structure of the program.

The fundamental data structures for the simplicial complex will be called class Simplex and class Triangulation. With these terms, we can now give a short overview of the program.

⟨ Header files to be included 7 ⟩

⟨ class Simplex 15 ⟩

⟨ class Triangulation 11 ⟩

⟨ Member functions of class Triangulation 18 ⟩ ⟨ Main program 54 ⟩

7. From LEDA we use the data types array, list, and we use streams for I/O. ⟨ Header files to be included 7 ⟩ ≡

#include <LEDA/array.h>
#include <LEDA/list.h>
#include <LEDA/stream.h>

See also sections 8 and 10. This code is used in section 6.

8. In order to show triangulations on the screen, we implement a function that draws the triangulation onto the screen in the two dimensional case using LEDA’s window type. Therefore, we have to include the appropriate LEDA header files. We are working with the X11R5 (xview) window environment. ⟨ Header files to be included 7 ⟩ +≡

#include <LEDA/window.h>
#include <LEDA/plane.h>

9. We have added to the LEDA-type matrix a function linear solver ( ), which solves a system of linear equations in the following sense: given a matrix A and a vector b, the function applies the Gaussian elimination algorithm for solving a linear system and returns a list of vectors which characterize the affine-linear space of the solution. When the list returned is empty, there is no solution at all. Otherwise, the list has the form (a, d_{1}, …, d_{r}), which means that d_{1}, …, d_{r} are the spanning vectors of the affine-space and a is a base point which lies in that space. Upon this function a function affine dependency ( ) is based, which determines for a given list of vectors whether they are affine-linearly dependent. Almost all of the calculations needed are carried out in these two functions. They are not yet part of LEDA, but soon will be. We give the code of these functions in Section 60. In the current version of LEDA, the entries of a vector and of a matrix are of type double, but in later versions the user will be able to choose between several numerical data types, for example double and rational. In the algorithm we use a macro number defined as double or rational, but it can be easily changed to any other type providing the operations =,+,–,*,/ and an absolute value function like fabs ( ). To choose rational arithmetics just define a macro RATIONAL. The macro number will be defined in linalg.h accordingly.

10. We include the appropriate header files depending on the kind of arithmetics we use. ⟨ Header files to be included 7 ⟩ +≡

#ifdef RATIONAL
#include "rat_matrix.h"

#else

#include <LEDA/matrix.h>

#endif #include "linalg.h" // see Sec. 64 11. The fundamental data structures.

Now we can begin to define our fundamental data structures (cf. Section 4, page 2). The whole simplicial complex will be managed by the class Triangulation. In this class, we store the coordinate vectors of the points given so far (list coordinates ), the dimension of the convex hull of these points (int dcur ), the dimension of the coordinate vectors of the input points (int dmax ) and a pointer to the origin simplex, from which we can reach all other simplices. also store the coordinates of a point which lies in the interior the origin simplex. When we compute the equation for the base facet of an unbounded simplex, it is useful to know a point which lies in the interior of the origin simplex (cf. Section 17) and we also need such a point as a starting point for the segment walking method. An appropriate point is the center point of the origin simplex

p = ∑_{i=0}^{dcur} \frac{v_{i}}{dcur + 1},

where v_{0}, ⋯, v_{dcur} are the coordinate vectors of the vertices of the origin simplex. To avoid the numerically problematic division, we store in the variable quasi center only the sum of the v_{i}’s and when we need p, we have to remind that p = quasi center /(dcur +1). Furthermore, we store a list of all constructed simplices (list all simplices ) which makes it easier to traverse all simplices (for instance in the destructor of the class or when displaying the simplicial complex). With this list, the interested reader may implement a copy constructor for the class.

During the insertion of some x_{i}, we have to find the x_{i}-visible facets of CH(π_{i−1}). For this purpose, we have implemented three search methods: the visibility search method and the segment walking method described in [4] and a modification of the visibility search method. For the selection of the search method, we introduce an enumeration type with the elements VISIBILITY,

MODIFIED_VISIBILITY and SEGMENT_WALK, respectively.

The public member method of Triangulation determines the search method to be used; it can be changed by the user at any time and its default value is SEGMENT_WALK. Each of these search methods stores its result (i.e., pointers to the unbounded simplices having the x_{i}-visible facets of CH(π_{i−1}) as base facets) in the list visible simplices .

As the representation of aff R, we use the (affine linearly independent) vertices of the origin simplex. ⟨ class Triangulation 11 ⟩ ≡

enum search method {

    VISIBILITY, SEGMENT_WALK, MODIFIED_VISIBILITY

};

class Triangulation { private: list⟨vector⟩ coordinates ; // the coordinate vectors of the x_{i} int dcur ; // dimension of the current convex hull int dmax ; // dimension of the coordinate vectors Simplex ∗origin simplex ; // pointer to the origin simplex vector quasi center ;

// sum of the coordinate vectors of the vertices of the origin simplex list⟨Simplex ∗⟩ all simplices ; // list of all simplices list⟨Simplex ∗⟩ visible simplices ; // result of search method ⟨ Further member declarations of class Triangulation 17 ⟩ void print (Simplex ∗); // writes some statistics about S to stdout

public:

search method method ; int searched simplices ; // used for statistical reasons int created simplices ( );

// returns the number of simplices that have been created void insert (const vector &x); // insertion routine void show (window &W);

// draws the triangulation onto the screen void print all ( ); // calls print() for all simplices Triangulation(int d, search method m = SEGMENT_WALK);

// constructor function with default argument ∼Triangulation( ); // destructor function

};

See also sections 12, 13, and 14. This code is used in section 6.

12. At the end of the program we want to be able to print the number of simplices that have been created. If a simplex is the k-th simplex created, its component sim nr gets k (cf. Simplex ::Simplex( )). ⟨ class Triangulation 11 ⟩ +≡

int Triangulation ::created simplices ( ) {

Simplex Dummy (2); static dummys created = 0; dummys created {+ \atop +}; return Dummy .sim nr − dummys created ;

}

13. The constructor for class Triangulation is easy to implement. The default search method is segment walking.

⟨ class Triangulation 11 ⟩ +≡ Triangulation ::Triangulation(int d, search method m) {

dcur = −1; dmax = d; searched simplices = 0; origin simplex = nil ; method = m;

}

14. In the destructor for Triangulation, we have to release the storage which was allocated for the simplices. ⟨ class Triangulation 11 ⟩ +≡

Triangulation ::∼Triangulation( )

{ Simplex ∗S; forall (S, all simplices ) delete (S);

}

15. Now we define the class Simplex. We make class Triangulation a friend of class Simplex, so that it can access every private member of class Simplex. For each simplex, we store its vertices as an array vertices of pointers to the corresponding occurrences in the list coordinates of class Triangulation. For the anti-origin we store nil . The array has length dmax + 1 since a simplex has at most dmax + 1 vertices. When the current hull has dimension dcur , only the array elements 0 to dcur are used. Furthermore, we use the following convention:

the peak vertex of the simplex is always vertices [0].

In order to represent the neighborhood relation, we use a second array neighbors , such that neighbors [k] points to the simplex opposite to the vertex vertices [k].

Given a vertex v of a simplex V, let W be the neighbor of V opposite to v. It is also useful to find the vertex w opposite to v, i.e., the vertex w of W which is not a vertex of V. For this purpose, we use an array opposite vertices : if v is the k-th vertex of V, i.e., V vertices [k] ≡ v, and w is the l-th vertex of W, then V opposite vertices [k] ≡ l and vice versa W opposite vertices [l] ≡ k.

Figure 3 illustrates the connection between two adjacent simplices V and W. The numbers that stand outside the simplices are the numbers of the vertices of V, the others being the numbers of the vertices of W. In both simplices, the vertex with number 0 is the peak vertex. The connectivity of V and W is expressed as follows: we have V neighbors [2] ≡ W and W neighbors [0] ≡ V, indicated

[figure: W; V; 2; 0; w; v; 0; 1; 2; 1]

Figure 3: The connection of two simplices V and W

by the corresponding arrows. Furthermore, we have V opposite vertices [2] ≡ 0 and vice versa W opposite vertices [0] ≡ 2.

For the test whether a point sees a facet of a given simplex, we need the normal vector and the right side of the equation normal · x = alpha of the hyperplane which contains the facet. The normal vector must be oriented such that the vertex opposite to the face lies in the positive halfspace. When we need normal or alpha for a facet i of a simplex S (i.e, the facet opposite to the i-th vertex of S), we call the function normal (S, i) or alpha(S, i), respectively (cf. Section 17), which are members of class Triangulation. Once we have computed these values for a facet, we store them in the arrays normal values and alpha values of the corresponding simplex so that they not have to be computed again when they are used the next time. Unfortunately, after a dimension jump all entries of normal values and alpha values become invalid. Therefore, we store in an array valid equations the time, i.e. the current dimension dcur of the convex hull, when the values of the facet’s equation was computed. The functions normal (S, i) and alpha(S, i) check whether the respective values of the i-th facet of simplex S are still valid and if not they compute them. Then they return the valid values. The values are invalid iff valid equations [i] < dcur . Then the values are computed anew and valid equations [i] is set to dcur . Initially, they are set to 0. (In dimension 0 we never need the plane equations).

In the implementation of the deletion process we must not forget that we may have to set valid equations [i] to 0 again for all i and all simplices if we delete a vertex which was a dimension jump.

We also need a mark to indicate visited simplices when we traverse the triangulation (e.g., for the visibility search or for traversing all simplices when we process a dimension jump).

#define anti origin nil

⟨ class Simplex 15 ⟩ ≡ class Simplex { friend class Triangulation;

// Triangulation has unrestricted access private:

int sim nr ; // useful for debugging; unique number for each simplex

#ifdef USE_LEDA_ARRAYS

array⟨list item⟩ vertices ;

// pointers to the coordinate vectors of the vertices array⟨Simplex ∗⟩ neighbors ; array⟨int⟩ opposite vertices ; // indices of opposite vertices array⟨vector⟩ normal values ; // normal vectors of the facets array⟨number⟩ alpha values ;

// right side of the equation normal · x = alpha array⟨int⟩ valid equations ;

// dimension in which corresponding equation was computed #else

list item ∗vertices ;

// pointers to the coordinate vectors of the vertices Simplex ∗∗neighbors ; int ∗opposite vertices ; // indices of opposite vertices vector ∗normal values ; // normal vectors of the facets number ∗alpha values ;

// right side of the equation normal · x = alpha int ∗valid equations ;

// dimension in which corresponding equation was computed #endif

bool visited ;

// used to mark simplices when traversing the triangulation Simplex(int dmax ); // constructor function

#ifdef USE_LEDA_ARRAYS

∼Simplex( ) { } ; // destructor function #else

∼Simplex( ); // destructor function #endif

LEDA_MEMORY(Simplex);

};

See also section 16. This code is used in section 6.

16. The constructor for class Simplex sets the size of the arrays, allocates storage for the normal vectors and marks the simplex as not visited. ⟨ class Simplex 15 ⟩ +≡

#ifdef USE_LEDA_ARRAYS

Simplex ::Simplex(int dmax ): vertices (0, dmax ), neighbors (0, dmax ), opposite vertices (0, dmax ), normal values (0, dmax ), alpha values (0, dmax ), valid equations (0, dmax )

{

static int lfdnr = 0;

// each simplex gets a unique number (for debugging) sim nr = lfdnr {+ \atop +}; for (int i = 0; i ≤ dmax ; i{+ \atop +}) \{

normal values [i] = vector(dmax ); // initialize the normal vectors neighbors [i] = nil ; // to avoid illegal pointers when using print ( )

} visited = false;

} #else

Simplex ::Simplex(int dmax ) {

static int lfdnr = 0;

// each simplex gets a unique number (for debugging) vertices = new list item [dmax + 1]; neighbors = new Simplex ∗ [dmax + 1]; opposite vertices = new int [dmax + 1]; normal values = new vector [dmax + 1]; alpha values = new number [dmax + 1]; valid equations = new int [dmax + 1]; sim nr = lfdnr {+ \atop +}; for (int i = 0; i ≤ dmax ; i{+ \atop +}) \{

normal values [i] = vector(dmax ); // initialize the normal vectors neighbors [i] = nil ; // to avoid illegal pointers when using print ( ) valid equations [i] = 0;

} visited = false;

} Simplex ::∼Simplex( ) // destructor function {

delete vertices ; delete neighbors ; delete opposite vertices ; delete normal values ; delete alpha values ; delete valid equations ;

} #endif 17. Computing plane equations.

We treat now the functions for computing the values normal and alpha for the plane equation

normal · x = alpha

for the facet of a simplex. In Section 15, we mentioned that we call functions normal (S, i) and alpha(S, i), respectively, when we need one of these values for the i-th facet of the simplex S. These functions check whether the values have already been computed (and if they are still valid) and if this is not the case, they call the function plane equation( ) to compute them. ⟨ Further member declarations of class Triangulation 17 ⟩ ≡

void plane equation(Simplex ∗S, int dir );

// defined as a member of Triangulation vector &normal (Simplex ∗S, int i); number alpha(Simplex ∗S, int i);

See also sections 25, 31, 39, and 45. This code is used in section 11.

18.⟨ Member functions of class Triangulation 18 ⟩ ≡ vector &Triangulation ::normal (Simplex ∗S, int i) {

if (S valid equations [i] ≠ dcur ) plane equation(S, i); return S normal values [i];

}

number Triangulation ::alpha(Simplex ∗S, int i)

{ if (S valid equations [i] ≠ dcur ) plane equation(S, i); return S alpha values [i];

}

See also sections 19, 24, 26, 33, 34, 35, 36, 37, 42, 44, 47, 51, 52, and 53. This code is used in section 6.

19. The parameter dir of the function plane equation( ) determines the facet of the simplex S whose equation has to be computed; it is the facet opposite to the dir -th vertex of S. The result is stored in S normal values [dir ] and S alpha values [dir ], respectively. The equation is normalized, such that the vertex dir lies in the positive halfspace. We first check whether we have already computed the plane equation for the neighbor of S opposite to the dir -th vertex. If so, we can take the values from there and only have to reverse the signs.

⟨ Member functions of class Triangulation 18 ⟩ +≡ void Triangulation ::plane equation(Simplex ∗S, int dir ) {

if (S neighbors [dir ] valid equations [S opposite vertices [dir ]] ≡ dcur ) { // if the corresponding equation of the neighbor is valid

S normal values [dir ] =

−S neighbors [dir ] normal values [S opposite vertices [dir ]]; S alpha values [dir ] =

−S neighbors [dir ] alpha values [S opposite vertices [dir ]]; } else {

⟨ Compute the plane equation 20 ⟩ } S valid equations [dir ] = dcur ;

}

20. The system of equations that leads us to normal and alpha is a little bit complicated. By a hyperplane, we always mean an affine-linear subspace of the dmax -dimensional space which has dimension dmax − 1. Let p_{0}, …, p_{dcur} be the vertices of S. Choose one vertex p_{base} ≠ p_{dir} of S (we actually choose base = 0 if dir ≠ 0 and base = 1 otherwise). It is clear, that normal must be orthogonal to all connecting vectors p_{i} − p_{base}, 0 ≤ i ≤ dcur , i ≠ base, i ≠ dir . This gives us dcur − 1 equations with dmax variables. If dcur = dmax this system of equations determines the vector normal (up to a multiplication with a constant factor).

If dcur < dmax , we also require that the hyperplane is orthogonal to the affine subspace spanned by our current set of points, i.e., we postulate that normal is a linear combination of the direction vectors that span the plane in which the current original simplex lies, that is, of the current space we are working in. If we didn’t postulate this, it would be possible that the computed hyperplane contains the vertex with number dir . Let q_{0}, …, q_{dcur} denote the vertices of the origin simplex. Then the above condition is equivalent to:

There are λ_{1}, …, λ_{dcur} such that

normal + ∑_{i=1}^{dcur} λ_{i}(q_{i} − q_{0}) = 0.  (1)

Thus, we have another dmax equations and another dcur variables. The corresponding system of equations looks as follows: (for saving space we abbreviate dcur as c, dmax as m and dir as d).

M · ( normal_{0} ⋮ normal_{m} − 1 λ_{1} ⋮ λ_{c} ) = ( 0 ⋮ 0 0 ⋮ 0 )

with

M = ( \begin{matrix}p_{1,0} − p_{0,0} & … & p_{1,m−1} − p_{0,m−1} & 0 & … & 0 \\ ⋮ & & ⋮ & ⋮ & & ⋮ \\ p_{d−1,0} − p_{0,0} & … & p_{d−1,m−1} − p_{0,m−1} & 0 & … & 0 \\ p_{d+1,0} − p_{0,0} & … & p_{d+1,m−1} − p_{0,m−1} & 0 & … & 0 \\ ⋮ & & ⋮ & ⋮ & & ⋮ \\ p_{c,0} − p_{0,0} & … & p_{c,m−1} − p_{0,m−1} & 0 & … & 0 \\ 1 & … & 0 & q_{1,0} − q_{0,0} & … & q_{c,0} − q_{0,0} \\ ⋮ & ^{…} & ⋮ & ⋮ & & ⋮ \\ 0 & … & 1 & q_{1,m−1} − q_{0,m−1} & … & q_{c,m−1} − q_{0,m−1} \end{matrix} ).

Altogether, we have only dmax +dcur −1 equations but dmax +dcur variables. This means that we will get an affine space of dimension 1 for possible normal vectors and raw linear solver ( ) will compute a base point and one direction vector. Since the right side of all equations is 0, the zero vector is a solution of the system (the system is homogeneous). Therefore, the direction vector itself lies in the space of solution. Because it can’t be the zero vector (then the dimension of the solution space were 0, a contradiction), we can simply take it as a non-trivial solution from which we obtain normal .

To compute alpha, the right side of the equation, we see that a point x lies in the hyperplane iff normal · (x − p_{base}) = 0. This means, that we can take normal · p_{base} as alpha. Finally, we must test in which halfspace the vertex dir lies. If it does not lie in the positive halfspace defined by normal and alpha, we simply change the sign of normal and alpha, then it does. Note that we have to treat here another special case: if the vertex selected by dir is the antiorigin, the hyperplane must separate the anti-origin from any point inside the simplicial complex, so we simply test the point quasi center /(dcur + 1), which is an interior point of the origin simplex. After this test, normal will point to the outside of the simplicial complex.

All we have to do now is to fill a C-array appropriately with the matrix and the right side, to feed it to raw linear solver ( ), extract normal from the solution, compute alpha and adjust the sign of normal . ⟨ Compute the plane equation 20 ⟩ ≡

⟨ Set up the matrix 21 ⟩

/∗ now we can call raw linear solver ( ), which returns us a list L characterizing the space of the solution ∗/ list⟨vector⟩ L; L = raw linear solver (M, rows , cols ); ⟨ Extract normal from L 22 ⟩ ⟨ Compute the right hand side alpha and adjust the sign of normal 23 ⟩

This code is used in section 19.

21. To avoid many news and deletes we reallocate the memory for the matrix only if the size of the matrix has changed with respect to the last call. If dcur = dmax , condition (1) is obviously satisfied and we can omit the corresponding part of the system. ⟨ Set up the matrix 21 ⟩ ≡

static int rows , cols ; // size of the matrix static number ∗∗M;

// pointer to memory allocated for matrix and right side of equation int i, row , col ; // for stepping through the matrix /∗ we only reallocate memory if the sieze of the matrix has changed ∗/ if (rows ≠ (dcur ≡ dmax ? dmax − 1 : dmax + dcur − 1) ∨

cols ≠ (dcur ≡ dmax ? dmax : dmax + dcur )) { if (M) { // don’t free memory prior to first “re”-allocation

for (row = 0; row < rows ; row {+ \atop +}) delete M[row ];

delete M; } rows = (dcur ≡ dmax ? dmax − 1 : dmax + dcur − 1); cols = (dcur ≡ dmax ? dmax : dmax + dcur ); M = new number ∗ [rows ]; for (row = 0; row < rows ; row {+ \atop +}) \{

M[row ] = new number [cols + 1]; // store right side also in M M[row ][cols ] = 0;

} } /∗ Now we set up the equations for: normal is orthogonal to all connecting vectors p_{i} −p_{base}, 0 ≤ i ≤ dcur , i ≠ base, i ≠ dir . With row , we step through the rows 0, …, dcur − 2 of M. With col , we fill the columns with the vectors p_{i} − p_{base} ∗/ int base = (dir ≡ 0 ? 1 : 0); vector difference(dcur ); // gets the value of p_{i} − p_{base} for (i = 0, row = 0; row < dcur − 1; row {+ \atop +}, i{+ \atop +}) \{

while ((i ≡ dir ) ∨ (i ≡ base)) i{+ \atop +};

// do not insert p_{dir} − p_{base} and p_{base} − p_{base} difference = coordinates .contents (S vertices [i]) −

coordinates .contents (S vertices [base]); // p_{i} − p_{base} for (col = 0; col < dmax ; col {+ \atop +}) M[row ][col ] = difference[col ]; } if (dcur ≠ dmax ) {

// we need the equations corresponding to condition (1) for (row = 0; row < dcur − 1; row {+ \atop +}) \{

for (col = dmax ; col < cols ; col {+ \atop +}) M[row ][col ] = 0; // zero matrix in the upper right corner

}

/∗ the equations for: normal lies in the affine hull ∗/ for (row = dcur − 1; row < rows ; row {+ \atop +})

for (col = 0; col < dmax ; col {+ \atop +}) M[row ][col ] = (col ≡ row ? 1 : 0);

// unit matrix in the lower left corner for (col = dmax ; col < cols ; col {+ \atop +}) \{

// fill M column by column with q_{col−dmax+1} − q_{0} difference =

coordinates .contents (origin simplex vertices [col − dmax + 1]) − coordinates .contents (origin simplex vertices [0]);

for (row = dcur − 1; row < rows ; row {+ \atop +})

M[row ][col ] = difference[row ]; }

}

This code is used in section 20.

22. We mentioned already that the space of solution of our system M · x = b has dimension 1 and that the base point (i.e., the first element of L) is the zero vector. So we take the second element of L as a non-trivial solution of the system. normal consists of the first dmax entries of the solving vector of the system. ⟨ Extract normal from L 22 ⟩ ≡

vector help;

// we will throw away the last dcur components of this vector /∗ remove the base point from the list as described above ∗/ L.pop( ); /∗ the direction vector is a non-trivial solution ∗/ help = L.pop( ); /∗ compute normal from help copying the first dmax entries ∗/ for (int j = 0; j < dmax ; j{+ \atop +}) S normal values [dir ][j] = help[j];

This code is used in section 20.

23. We compute alpha and adjust the sign of normal .

⟨ Compute the right hand side alpha and adjust the sign of normal 23 ⟩ ≡

S alpha values [dir ] = S normal values [dir ] ∗

coordinates .contents (S vertices [base]); if (S vertices [dir ] ≡ anti origin) { // if vertex dir is the anti-origin

/∗ in the following expression we have to do an explicit cast since otherwise several overloaded operator ∗ ( )-functions would match this use. ∗/ if (S normal values [dir ] ∗ quasi center >

S alpha values [dir ] ∗ number(dcur + 1)) {

// does normal point to the interior of the triangulation? S normal values [dir ] = −S normal values [dir ]; S alpha values [dir ] = −S alpha values [dir ];

} } else {

if (S normal values [dir ] ∗ coordinates .contents (S vertices [dir ]) < S alpha values [dir ]) { // does p_{dir} lie in the negative halfspace ?

S normal values [dir ] = −S normal values [dir ]; S alpha values [dir ] = −S alpha values [dir ];

} }

This code is used in section 20.

24. The insert procedure.

We treat now the insertion procedure as described in Section 5. For the insertion of a point x, we distinguish three cases:

• x is the first point to be inserted.

• x is a dimension jump (and not the first point to be inserted).

• x is not a dimension jump.

⟨ Member functions of class Triangulation 18 ⟩ +≡ void Triangulation ::insert (const vector &x) { /∗ add x to the points already inserted and store its position in item x ∗/

list item item x = coordinates .append (x);

if (dcur ≡ −1) { // x is the first point to be inserted

⟨ Initialize the triangulation 27 ⟩ } else if ((dcur < dmax ) ∧ is dimension jump(x)) { // see Section 26

⟨ Dimension jump 46 ⟩ } else {

⟨ Non-dimension jump 28 ⟩ }

}

25. We need a function is dimension jump( ), which tells us whether x is a dimension jump or not. ⟨ Further member declarations of class Triangulation 17 ⟩ +≡

bool is dimension jump(const vector &x);

26. How can we test whether x is a dimension jump? x is a dimension jump iff x does not lie in the affine hull of the vertices of the origin simplex. Since all these vertices are affine-linearly independent by our construction, we only have to test whether x and all these vertices are affine-linearly dependent. We test this by using the function affine dependency ( ) (cf. Section 75), which gets as argument a list of all the vectors to test. ⟨ Member functions of class Triangulation 18 ⟩ +≡

bool Triangulation ::is dimension jump(const vector &x)

{ list⟨vector⟩ L; // the list for affine dependency ( ) /∗ we insert x and all vertices of the origin simplex into L ∗/ L.push(x); for (int i = 0; i ≤ dcur ; i{+ \atop +})

L.push(coordinates .contents (origin simplex vertices [i])); return (¬affine dependency (L));

}

27. When the first point x is inserted, we must initialize our triangulation, that means, we must build the first simplices by hand. This is easy to do. When we only have one point, the simplicial complex consists of two simplices: the origin simplex, containing x as peak, and an outer simplex outer simplex having the anti-origin as its peak. They both point to one another in a natural way. The origin simplex has no base facet by definition, and because dcur is 0 outer simplex has a (−1)-dimensional base facet, that means, it has no base facet either. The center point of the origin simplex is clearly x. ⟨ Initialize the triangulation 27 ⟩ ≡

Simplex ∗outer simplex ; // a pointer to the outer simplex

dcur = 0; // we jump from dimension -1 to dimension 0

all simplices .append (origin simplex = new Simplex (dmax ));

all simplices .append (outer simplex = new Simplex (dmax ));

origin simplex vertices [0] = item x ; // x is the only point and the peak origin simplex neighbors [0] = outer simplex ; origin simplex opposite vertices [0] = 0; outer simplex vertices [0] = anti origin; outer simplex neighbors [0] = origin simplex ; outer simplex opposite vertices [0] = 0; quasi center = x;

This code is used in section 24.

28. We discuss now how to handle insertions that are not dimension jumps. We first compute the set of all x-visible hull facets. This is described in detail in Section 31. As a result of this step, we get in visible simplices the list of unbounded simplices whose base facets see x. If there are none, then x lies within the current hull and we are done. Otherwise, we have to modify some simplices and to add some new ones as described in Section 5. Also the neighborhood information has to be updated. ⟨ Non-dimension jump 28 ⟩ ≡

⟨ Find x-visible hull facets 32 ⟩

if (¬visible simplices .empty ( )) {

list⟨Simplex ∗⟩ NewSimplices ;

// Simplices created to store horizon ridges Simplex ∗S; forall (S, visible simplices ) {

/∗ For each x-visible facet F of CH(π_{i−1}) alter the simplex S(F ∪ \{\overline{O}\}) of \overline{Δ}(π_{i−1}) into S(F ∪ \{x_{i}\}). Note that \overline{O} is the peak, i.e., S vertices [0]. ∗/ S vertices [0] = item x ; ⟨ For each horizon ridge add the new simplex 29 ⟩

}

visible simplices .clear ( );

⟨ Update the neighborhood relationship 30 ⟩ }

This code is used in section 24.

29. We now describe, how to update the neighborhood relationship and to compute the equations of the base facets of the new simplices.

At this point, we have found the current hull facets seeing x, in the form of the simplices whose base facets see x and with the anti-origin as their peak vertex. Let V be the set of such simplices. Now we update T by altering these simplices, and creating some others. The alteration is simply to replace the anti-origin with x in every simplex in V.

The new simplices correspond to new hull facets. Such facets are the hull of x and a horizon ridge f; a horizon ridge is a (d − 2)–dimensional face of conv R with the property that exactly one of the two incident hull facets sees x. Each horizon ridge f gives rise to a new simplex A_{f} with base facet conv(f ∪ \{x\}) and peak \overline{O}. For each horizon ridge of conv R there is a non-base facet G of a simplex in V such that x does not see the base facet of the other simplex incident to the facet G. Thus the set of horizon ridges is easily determined. (cf. [4], p. 11)

The figures 4 and 5 illustrate the situation. In figure 4, x sees the facets conv(f, g) and conv (g, h). There are two horizon ridges: the points f and h. The non-base facet G of the above text is the segment s which x does not see. In figure 5, x has been inserted. Two new unbounded simplices corresponding to the two horizon ridges have been added.

We find all horizon ridges incident to an updated simplex S with x-visible base facet by testing all its neighbors (except for the one opposite to its peak) whether their base facet is x-visible. If the base facet of a neighbor is not xvisible, we have found a horizon ridge f and have to create a new simplex T with base facet conv(f ∪ \{x\}) and peak \overline{O}. We collect all new simplices in the list NewSimplices .

We use the index k to run through the neighbors of S. When we have identified a horizon ridge, the vertices of the new simplex T are the vertices of S with the k-th vertex replaced by x. The peak of T is the anti-origin \overline{O}. We could therefore initialize the vertex set of the new simplex T by

T vertices = S vertices ;

T vertices [k] = item x ;

T vertices [0] = anti origin;

[figure: .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; f; g; h; x; s]

Figure 4: Before x is inserted

[figure: .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; f; g; h; x; s]

Figure 5: After x has been inserted In order to facilitate the update of the neighborhood relation, we proceed slightly differently: we make x the highest numbered vertex of T, i.e., we replace the second line by

T vertices [k] = S vertices [dcur ]; T vertices [dcur ] = item x ;

What are the neighbors of the new simplex T? The neighbor opposite to \overline{O} is S and the neighbor opposite to x is the neighbor of the old S (i.e., S before the replacement of its peak \overline{O} by x) incident to f ∪ \overline{O}. The neighbors opposite to the j-th vertex of T, with 1 ≤ j < dcur, are computed in the next section. ⟨ For each horizon ridge add the new simplex 29 ⟩ ≡

for (int k = 1; k ≤ dcur ; k{+ \atop +}) \{

if (normal (S neighbors [k], 0) ∗ x ≤ alpha(S neighbors [k], 0)) {

// x doesn’t see the base facet of the neighbor Simplex ∗T = new Simplex (dmax ); all simplices .append (T); NewSimplices .append (T); /∗ Take the vertices of S as the vertices of the new simplex, replacing the current vertex by the dcur -th, the first by x and the peak by \overline{O} ∗/ int ii ; for (ii = 0; ii ≤ dcur ; ii {+ \atop +}) T vertices [ii ] = S vertices [ii ]; T vertices [k] = S vertices [dcur ]; T vertices [dcur ] = item x ; T vertices [0] = anti origin; /∗ set the pointers to the two neighbors we already know and update the corresponding entries in the opposite vertices -arrays ∗/ T neighbors [dcur ] = S neighbors [k]; T opposite vertices [dcur ] = S opposite vertices [k]; T neighbors [0] = S; T opposite vertices [0] = k; /∗ Also set the reverse pointers from those two neighbors to the new simplex ∗/ S neighbors [k] neighbors [S opposite vertices [k]] = T; S neighbors [k] opposite vertices [S opposite vertices [k]] = dcur ; S neighbors [k] = T; S opposite vertices [k] = 0;

} }

This code is used in section 28.

30. We now complete the update of the neighborhood relation. How the neighborhood relationship has to be updated is described in [4] as follows.

24. THE INSERT PROCEDURE 23

It remains to update the neighbor relationship. Let A_{f} = S(conv(f ∪\{x\}), \overline{O}) be a new simplex corresponding to horizon ridge f. In the old triangulation (before adding x) there were two simplices \overline{V} and N incident to the facet conv(f ∪ \{\overline{O}\}); \overline{V} ∈ V ^{3} and N ∉ V. In the updated triangulation \overline{V} is replaced by a new simplex V that has the same base but peak x. The neighbor of A_{f} opposite to x is N and the neighbor opposite to \overline{O} is V. Now consider any vertex q ∈ f and let S = S_{f,q} be the set of simplices with peak x and including vertex(f) \ {q} ∪ {x} in their vertex set; for a face f we use vertex(f) to denote the set of vertices contained in f. We will show that the neighbor of A_{f} opposite to q can be determined by a simple walk through S. This walk amounts to a rotation about the (d−2)–face conv(vertex(f)\{q}∪{x}). Note first that V ∈ S. Consider next any simplex S = S(F, x) ∈ S. Then F = conv(f \ {q} ∪ {y_{1}, y_{2}}) for some vertices y_{1} and y_{2}. Thus S has at most two neighbors in S, namely the neighbors opposite to y_{1} and y_{2} respectively. Also, V has at most one neighbor in S, namely the neighbor opposite to q (Note that the neighbor opposite to y, where conv(f ∪ {y}) is the base facet of V, is the simplex A_{f} ∉ S.). The neighbor relation thus induces a path on the set S with V being one end of the path. Let V′ with base facet conv(f \ {q} ∪ {y_{1}, y_{2}}) be the other end of the path. Assume that the neighbor of V′ opposite to y_{1}, call it B, does not belong to S and that y_{1} = q if V = V′, i.e., the path has length zero. The simplex B includes vertex(f)\{q}∪{y_{2}, x} in its vertex set and does not have peak x. Thus B has peak \overline{O} and hence B is the neighbor of A_{f} opposite to q. This completes the description of the update step. (cf. [4], p. 11)

[figure: .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; f, q, y_{1}; N; \overline{O}; A_{f}; x; \overline{O}; B; y_{2}, y′_{1}; y′_{2}; \overline{V} , V; V′]

Figure 6: Updating the neighborhood relation

Figure 6 illustrates the situation described above in the two dimensional case. y′_{1} and y′_{2} are the new values of y_{1} and y_{2} after one rotation around x.

3V is the set of outer simplices which see x

This is the only rotation to be made. Then the neighbor of q with respect to A_{f} is found. It is B.

We implement the update of the neighborhood information as follows. For all new simplices corresponding to horizon ridges, the pointers to the neighbors opposite to x and \overline{O} are already set (cf. the previous section). It remains to do the following for every new simplex A_{f} corresponding to horizon ridge f:

For all vertices q of A_{f} except x and \overline{O} find the neighbor of A_{f} opposite to q and set the corresponding neighbor pointer.

Note that we do not need to set the pointer from the neighbor we have found to A_{f}, since the neighbor is also a new simplex and hence this pointer will be (or has been) set anyhow.

Determining the neighbor of A_{f} opposite to q is done as follows. We walk through the simplices T along the path through S starting at T = V = Af neighbors [0] as described in [4]. As long as T ∈ S (i.e., the peak of T is x) we go to the neighbor T′ of T opposite to y_{1} (for T = V we have y_{1} = q). The new y_{1} is the node of T′ equal to the vertex y_{2} of T. We store the indices of the vertices corresponding to y_{1} and y_{2} in two variables y1 and y2 respectively. In V, y_{2} is the vertex opposite to \overline{O} with respect to A_{f}. If T′ ∉ S (i.e., the peak of T′ is not x) we have found the neighbor B of A_{f} opposite to q. ⟨ Update the neighborhood relationship 30 ⟩ ≡

Simplex ∗Af ;

forall (Af , NewSimplices ) {

for (int k = 1; k < dcur ; k{+ \atop +}) \{

// for all vertices q of Af except x and \overline{O} find the opposite neighbor Simplex ∗T = Af neighbors [0]; for (int y1 = 0; T vertices [y1 ] ≠ Af vertices [k]; y1 {+ \atop +}) ; // exercise: show that we can also start with y1 = 1

int y2 = Af opposite vertices [0];

while (T vertices [0] ≡ item x ) { // while T ∈ S /∗ find new y_{1} ∗/ for (int new y1 = 0; T neighbors [y1 ] vertices [new y1 ] ≠

T vertices [y2 ]; new y1 {+ \atop +}) ; // exercise: show that we can also start with new y1 = 1

y2 = T opposite vertices [y1 ];

T = T neighbors [y1 ];

y1 = new y1 ; } Af neighbors [k] = T; // update the neighborhood relationship Af opposite vertices [k] = y1 ; // update the opposite neighbor

} }

This code is used in section 28.

31. FINDING X-VISIBLE HULL FACETS 25

31. Finding x-visible hull facets.

For finding the x-visible hull facets, we have implemented three search methods. The first method, the visibility search method, visits all simplices with x-visible base facet using depth first search starting in the origin simplex. It is implemented in the function visibility search( ).

The second method is a modification of the visibility search method. The difference is that if it has once reached an outer simplex, it restricts its search space to unbounded simplices. It uses the function search to outside( ), which is similar to visibility search( ) except that it stops when it has reached an unbounded simplex. It returns a pointer to the unbounded simplex that it has reached or nil if x lies in the interior of the hull. If it has reached an outer simplex, all unbounded x-visible simplices are collected using the function collect outer simplices ( ).

The third method is the segment walking method. This method walks through the simplices which are intersected by a ray Ox^{− ; − ; →} from a point O in the origin simplex to x. It returns a pointer to the simplex it has reached (even if this is a bounded simplex). The unbounded x-visible simplices are also collected using the function collect outer simplices ( ).

The visibility search method and the function collect outer simplices ( ) mark visited simplices as visited using the visited variable. We unmark them using the function clear visited marks ( ). ⟨ Further member declarations of class Triangulation 17 ⟩ +≡

void visibility search(Simplex ∗S, const vector &x); Simplex ∗search to outside(Simplex ∗S, const vector &x); Simplex ∗segment walk (const vector &x); void collect visible simplices (Simplex ∗S, const vector &x); void clear visited marks (Simplex ∗S);

32. The variable method which can be changed interactively by the user switches between the several search methods. ⟨ Find x-visible hull facets 32 ⟩ ≡

Simplex ∗last simplex ; // the simplex in which modified visibility search // and segment walking have stopped

switch (method ) {

case VISIBILITY: visibility search(origin simplex , x);

// generates list of unbounded simplices with x-visible base facet clear visited marks (origin simplex ); break;

  case MODIFIED_VISIBILITY:

last simplex = search to outside(origin simplex , x); if (last simplex ≠ nil ) // if x is not an interior point

collect visible simplices (last simplex , x);

// generates list of unbounded simplices with x-visible base facet clear visited marks (origin simplex ); break;

  case SEGMENT_WALK:

default: last simplex = segment walk (x);

if (last simplex vertices [0] ≡ anti origin) { // if x is not an interior point collect visible simplices (last simplex , x);

// generates list of unbounded simplices with x-visible base facet clear visited marks (last simplex );

}

break; }

This code is used in section 28.

33. How we can implement visibility search( ) is described in Section 5: starting at the origin simplex, we visit any unvisited neighbor of a visited simplex that has an x-visible base facet. Note that by this rule, we do not have to test the origin simplex (which by definition has indeed no base facet). The class Triangulation has a member list visible simplices , in which we store the outer simplices seeing x. The function visibility search( ) is recursive and gets as arguments a reference to the vector x and a pointer ∗S to the simplex to be visited. ⟨ Member functions of class Triangulation 18 ⟩ +≡

void Triangulation ::visibility search(Simplex ∗S, const vector &x) {

searched simplices {+ \atop +}; // only for statistical reasons

S visited = true; // we have visited S and never come back...

for (int i = 0; i ≤ dcur ; i{+ \atop +}) \{ Simplex ∗T = S neighbors [i]; // for all neighbors T of S if (¬T visited ) { // if the i-th neighbor has not been visited yet

if (normal (T , 0) ∗ x > alpha(T , 0)) {

// if x sees the base facet of the i-th neighbor if (T vertices [0] ≡ anti origin)

// if the i-th neighbor is an outer simplex visible simplices .push(T);

// we have found a visible simplex and store it visibility search(T , x); // do the recursive search

}

}

} }

34. Here is the first part of the possibly faster modified visibility search method: search from the origin simplex to the outside, then search on the outer facets recursively with depth first search. If x is an outer point, that means it is contained in one of the outer simplices, the function returns a pointer to the first outer simplex that is found. If x is an inner point, the function returns nil . When we say “possibly faster”, we have in mind that the searching to the outside (which is nothing but depth first search) will take exactly the same way as the normal visibility search if x is an interior point, so the time we have spent is unfortunately the same. It might only be faster if x lies not in the current hull. ⟨ Member functions of class Triangulation 18 ⟩ +≡

Simplex ∗Triangulation ::search to outside(Simplex ∗S, const vector

&x)

{ searched simplices {+ \atop +}; // only for statistical reasons S visited = true; // we have visited S and never come back... for (int i = 0; i ≤ dcur ; i{+ \atop +}) \{

Simplex ∗T = S neighbors [i]; // for all neighbors T of S if (¬T visited ) // if the i-th neighbor has not been visited yet

if (normal (T , 0) ∗ x > alpha(T , 0)) {

// if x sees the base facet of the i-th neighbor if (T vertices [0] ≡ anti origin)

// if the i-th neighbor is an outer simplex return T; // we have found to the outside

Simplex ∗result = search to outside(T , x); if (result ≠ nil ) return result ;

} } return nil ;

}

35. Now we collect all outer simplices which are visible from x. The collection process starts from an outer simplex S. ⟨ Member functions of class Triangulation 18 ⟩ +≡

void Triangulation ::collect visible simplices (Simplex ∗S, const vector

&x)

{ searched simplices {+ \atop +}; // only for statistical reasons S visited = true; // we have visited S and never come back... visible simplices .push(S); // store S as a visible simplex for (int i = 0; i ≤ dcur ; i{+ \atop +}) \{

Simplex ∗T = S neighbors [i]; // for all neighbors T of S if (¬T visited ∧ T vertices [0] ≡ anti origin)

// if the i-th neighbor has not been visited yet // and is an outer simplex if (normal (T , 0) ∗ x > alpha(T , 0))

// if x sees the base facet of the i-th neighbor collect visible simplices (T , x); // do the recursive collecting

} }

36. After a visibility search, we always must clear the visited -bits of the visited simplices. This is done by the recursive function clear visited marks ( ). It is very similar to the function visibility search( ). When we start this function, we also call it with the origin simplex as its argument. ⟨ Member functions of class Triangulation 18 ⟩ +≡

void Triangulation ::clear visited marks (Simplex ∗S)

{ S visited = false; // clear the visited -bit for (int i = 0; i ≤ dcur ; i{+ \atop +}) // for all neighbors of S

if (S neighbors [i] visited ) // if the i-th neighbor has been visited clear visited marks (S neighbors [i]); // clear its bit recursively

}

37. The following function implements the segment walk method to find the simplex containing the point x. Let O denote the “origin”, i.e., any point in the origin simplex; we can take quasi center /(dcur + 1) for O. The strategy is very simple: we start at the origin simplex and walk along the ray Ox^{− ; − ; →} through the simplices intersected by this ray until we reach the simplex containing x. There might be degenerate inputs, that means that the ray Ox^{− ; − ; →} passes through a vertex of the triangulation or that a segment of Ox^{− ; − ; →} lies whithin a facet of a traversed simplex. In order to deal with such degenerate inputs we will perturb O. The details of the perturbation will be described in Section 42. The perturbation scheme is similar to the well-known perturbation method for the simplex algorithm.

Assume now that on our walk we entered a simplex S through the facet in (i.e., the facet opposite to the in-th vertex of S) and we want to find out through which facet Ox^{− ; − ; →} leaves S. The points on the ray Ox^{− ; − ; →} satisfy the equation

r(λ) = O + λ(x − O) with 0 ≤ λ ≤ 1.

The hyperplane which contains the i-th facet of S is given by the equation

n^{(i)}p = α^{(i)},

where n^{(i)} = normal (S, i) is the normal vector and α^{(i)} = alpha(S, i) is the right side of the plane equation of the i-th facet of S. Thus substituting r(λ) for p, we get

λ^{(i)} = \frac{α^{(i)} − n^{(i)}O}{n^{(i)}(x − O)} ,

and this λ^{(i)} yields the intersection of the ray with the hyperplane. We do not need to worry about the case that the denominator in the above expression could be zero, because we never really compute any λ^{(i)} (this is due to the perturbation method described in Section 42). Since the vector x − O points from O to x, the values of λ increase along Ox^{− ; − ; →}. Hence the facet out through which we leave S is the facet for which λ^{(out)} is the smallest λ^{(i)} which is larger than λ^{(in)}. We really say larger and not larger or equal (which might be the case in the above equation), because by the perturbation, two different hyperplanes intersect Ox^{− ; − ; →} in two different points except if they both contain x, i.e., λ = 1. Since we only search for out if x ∉ S, we have λ^{(in)} < λ^{(out)} < 1 and thus λ^{(out)} ≠ λ^{(i)} for all i ≠ out . How we have already mentioned, we do not actually compute the values of the λ^{(i)}. We only need to know for some i and j whether λ^{(i)} < λ^{(j)} or not. This decision will be made by a function lambda cmp( ) (cf. Section 42). Only in this function, the perturbation of O plays a role.

The two variables in and out tell us the number of the facet through which we have entered and through which we will leave the current simplex, respectively. They change, when we walk from simplex to simplex. When we start our walk, there is no facet through which we have entered the current simplex (which is the origin simplex). This is indicated by setting the variable in to −1 at the beginning of segment walk ( ). In this case, we search for the facet of the origin simplex which is intersected by the “negative” part of Ox^{− ; − ; →}, i.e., we search the facet whose hyperplane intersects Ox^{− ; − ; →} with the largest possible negative value for λ and consider this as the entry facet. This gives us a starting value for in. We stop our walk, if we have found the simplex containing x (this might be an unbounded simplex, of course). We use two arrays nx and nO to store the values of the scalar products n^{(i)}x and n^{(i)}O which we will need several times.

When we have found the simplex S containing x, we stop and return it. ⟨ Member functions of class Triangulation 18 ⟩ +≡

Simplex ∗Triangulation ::segment walk (const vector &x)

{ Simplex ∗S = origin simplex ; // we start at the origin simplex bool x in S ;

// indicates whether we have found the simplex containing x int in = −1; // entry facet of the origin simplex

int i; // for treating every facet of S

#ifdef USE_LEDA_ARRAYS

array⟨number⟩ nx (0, dcur ); // scalar products n^{(i)}x array⟨number⟩ nO (0, dcur ); // scalar products n^{(i)}O

#else number ∗nx = new number [dcur + 1]; // scalar products n^{(i)}x number ∗nO = new number [dcur + 1]; // scalar products n^{(i)}O

#endif

while (true) { searched simplices {+ \atop +}; // only for statistical reasons ⟨ Compute the arrays nx and nO and test whether x ∈ S 38 ⟩ if (x in S ) {

#ifndef USE_LEDA_ARRAYS

delete nx ;

delete nO ;

#endif

return S; } /∗ We cannot set the next statement in front of this while-loop because we need the values of the arrays nx and nO to compute the initial value of in ∗/ if (in ≡ −1) { // if we are still in the origin simplex

⟨ Find the facet with largest λ < 0 40 ⟩ } ⟨ Go to the next Simplex on the ray Ox^{− ; − ; →} 41 ⟩

} }

38. The computation of the arrays is easily done by scalar multiplications. At the same time, we can test whether x lies in the current simplex. ⟨ Compute the arrays nx and nO and test whether x ∈ S 38 ⟩ ≡

x in S = true; // remains true until we find a facet which doesn’t see x if (S vertices [0] ≠ anti origin)

// otherwise we have reached the outside for (i = 0; i ≤ dcur ; i{+ \atop +}) \{

nx [i] = normal (S, i) ∗ x; nO [i] = normal (S, i) ∗ quasi center ; // this is n^{(i)}O(dcur + 1) if (nx [i] < alpha(S, i)) x in S = false;

}

This code is used in section 37.

39. When we start the segment walking method in the origin simplex we have to find the facet whose corresponding hyperplane intersects Ox^{− ; − ; →} with the largest λ < 0. For the comparison of two λ’s, we use the function lambda cmp( ) defined in Section 42. This function contains the actual perturbation method. A call lambda cmp(S, nx [i], nO [i], i, nx [j], nO [j], j) returns true iff λ^{(i)} < λ^{(j)}. We also need a function lambda negative( ) defined in Section 44, which decides whether λ^{(i)} is negative or not. ⟨ Further member declarations of class Triangulation 17 ⟩ +≡

bool lambda cmp(Simplex ∗S, number nix , number niO , int i, number

njx , number njO , int j); bool lambda negative(Simplex ∗S, number nx , number nO , int i);

40. Now we can search for the “starting facet” in the origin simplex. ⟨ Find the facet with largest λ < 0 40 ⟩ ≡

for (i = 0; i ≤ dcur ; i{+ \atop +}) \{

if (lambda negative(S, nx [i], nO [i], i)) // λ^{(i)} < 0 ?

if (in ≡ −1 ∨ lambda cmp(S, nx [in], nO [in], in, nx [i], nO [i], i)) in = i; }

This code is used in section 37.

41. Now we describe how to find the facet with number out through which we will leave the current simplex. We have to find out such that for all i ≠out either λ^{(i)} < λ^{(in)} or λ^{(out)} < λ^{(i)}. Note that we have λ^{(i)} = λ^{(j)} for i ≠ j only if λ^{(i)} = λ^{(j)} = 1 since we perturb O (cf. Section 42). But λ^{(out)} < 1 always holds since x ∉ S. When we enter the new simplex we have to set in to the index of the facet over which we enter it. The index of a facet is the index of the vertex opposite to it. Hence we set in to S opposite vertices [out ].

⟨ Go to the next Simplex on the ray Ox^{− ; − ; →} 41 ⟩ ≡ int out = −1; // we have not yet found the desired facet for (i = 0; i ≤ dcur ; i{+ \atop +}) \{

if (i ≡ in) continue; // we surely do not go back on the ray...

if (lambda cmp(S, nx [in], nO [in], in, nx [i], nO [i], i)) // if λ^{(in)} < λ^{(i)} if (out ≡ −1 ∨ lambda cmp(S, nx [i], nO [i], i, nx [out ], nO [out ], out ))

// if we have not yet a candidate for out

// or if the current λ^{(i)} is better then the old one out = i;

}

in = S opposite vertices [out ];

S = S neighbors [out ];

This code is used in section 37.

42. It remains to describe how we decide whether λ^{(i)} < λ^{(j)}. As we have already mentioned we perturb^{4} O to get O^{E} := O+E, where E = (ϵ, ϵ^{2}, …, ϵ^{dcur}) for some sufficiently small ϵ > 0. Thus, if we consider ϵ small enough, the perturbated ray will not go through any vertex or any other intersection of the hyperplanes of the triangulation. Therefore, the facets which are intersected by this ray are totally linearly ordered. We have λ^{(i)} < λ^{(j)} iff

\frac{α^{(i)} − n^{(i)}O^{E}}{n^{(i)}(x − O^{E})} < \frac{α^{(j)} − n^{(j)}O^{E}}{n^{(j)}(x − O^{E})}  (2)

⟺

n^{(j)}(x − O^{E})(α^{(i)} − n^{(i)}O^{E}) <^{σ} n^{(i)}(x − O^{E})(α^{(j)} − n^{(j)}O^{E}),  (3)

where <^{σ} is defined to be < if the signs of the two denominators are equal and > otherwise. If n(x − O) ≠ 0, the sign of n(x − O^{E}) is equal to the sign of n(x − O), because we can choose ϵ sufficiently small. Otherwise, we have

n(x − O^{E}) = n(x − O) − nE = −nE,

and thus the sign of the expression is equal to the inverse sign of n_{k}, where k is the smallest index for which n_{k} ≠ 0 holds (because the term n_{k}ϵ^{k} is dominating).

After a horrendous computation, we find out that (3) is equivalent to

α^{(i)}n^{(j)}x − α^{(i)}n^{(j)}O − (n^{(j)}x)(n^{(i)}O) − (α^{(i)}n^{(j)} + (n^{(j)}x)n^{(i)})E
<^{σ} α^{(j)}n^{(i)}x − α^{(j)}n^{(i)}O − (n^{(i)}x)(n^{(j)}O) − (α^{(j)}n^{(i)} + (n^{(i)}x)n^{(j)})E

If the parts which do not depend on E are different, they can be used to decide whether λ^{(i)} < λ^{(j)} (because ϵ can be chosen arbitrarily small, so that the absolute value of the parts depending on E become smaller than any positive number). Otherwise, we take the minimal k such that

α^{(i)}n_{k}^{(j)} + (n^{(j)}x)n_{k}^{(i)} ≠ α^{(j)}n_{k}^{(i)} + (n^{(i)}x)n_{k}^{(j)}

and can also make our decision.

The computations for the comparison are carried out in the function lambda cmp( ), which is a member function of Triangulation and gets the current simplex S as an input parameter. It also gets the values of the scalar products n^{(i)}x, n^{(i)}O · (dcur + 1), n^{(j)}x and n^{(j)}O · (dcur + 1) in the form of the input parameters nix , niO , njx and njO , respectively. Furthermore, it gets the indices i and j of the hyperplanes to be compared in order to access alpha(S, i), alpha(S, j), normal (S, i) and normal (S, j).

We first have to decide upon <^{σ}. For this reason, we use a variable sigma which is true iff <^{σ} is >.

4This perturbation method was proposed by Kurt Mehlhorn

In the following, we have to do several explicit casts to number, since otherwise several overloaded operator∗ ( )-functions would match the corresponding expressions which leads to compiler errors.

We need sevaral times the expression njx ∗ number(dcur + 1) − njO and nix ∗ number(dcur + 1) − niO , respectively. Tests have shown that in large examples we spent a significiant amount of time evaluating these expressions, so we do it only once in each call of lambda cmp( ) and store its value in the variable njx njO and nix niO , respectively. ⟨ Member functions of class Triangulation 18 ⟩ +≡

bool Triangulation ::lambda cmp(Simplex ∗S, number nix , number niO , int i, number njx , number njO , int j)

{ bool sigma = false; // we assume <^{σ} to be < number njx njO = njx ∗ number(dcur + 1) − njO ; number nix niO = nix ∗ number(dcur + 1) − niO ; ⟨ Decide whether <^{σ} is < or > 43 ⟩ /∗ Now we test the parts which are not depending on E. Remember that niO = n^{(i)}O(dcur + 1), so we have to “adjust” nix and njx ∗/ number diff = njx njO ∗ alpha(S, i) − njx ∗ niO

− (nix niO ∗ alpha(S, j) − nix ∗ njO ); if (diff < 0) // λ^{(i)} <^{σ} λ^{(j)}

return ¬sigma; else if (diff > 0) return sigma; else { // the comparison depends on the factor of E

for (int k = 0; k < dcur ; k{+ \atop +}) \{

diff = −(alpha(S, i) ∗ normal (S, j)[k] + njx ∗ normal (S, i)[k]) + (alpha(S, j) ∗ normal (S, i)[k] + nix ∗ normal (S, j)[k]);

if (diff ≠ 0) // λ^{(i)} <^{σ} λ^{(j)} if diff < 0

return (diff < 0 ? true : false); } /∗ If the program reaches the next line of code, we have k ≡ dcur . This means that even the E- factors are equal. In this case, we have λ^{(i)} = λ^{(j)} and can return false ∗/ return false;

} }

43. Here we examine the direction of <^{σ}. ⟨ Decide whether <^{σ} is < or > 43 ⟩ ≡

/∗ We examine the denominators of inequality (2); we first look at the left denominator ∗/ if (nix niO < 0) // again, remember that niO = n^{(i)}O(dcur + 1)

sigma = ¬sigma; else if (nix niO ≡ 0) {

for (int k = 0; ; k{+ \atop +})

// search the smallest k with normal (S, i)[k] ≠ 0 if (normal (S, i)[k] ≠ 0) {

if (normal (S, i)[k] > 0) sigma = true; else sigma = false; break;

} } /∗ and now at the right one ∗/ if (njx njO < 0) // again, remember that njO = n^{(j)}O(dcur + 1)

sigma = ¬sigma; else if (njx njO ≡ 0) {

for (int k = 0; ; k{+ \atop +})

// search the smallest k with normal (S, j)[k] ≠ 0 if (normal (S, j)[k] ≠ 0) {

if (normal (S, j)[k] > 0) sigma = ¬sigma; break;

} }

This code is used in section 42.

44. To decide upon the sign of a λ^{(i)}, we use the function lambda negative( ) which is similar to lambda cmp( ) and returns true iff λ^{(i)} < 0. Since

λ^{(i)} = \frac{α^{(i)} − n^{(i)}O}{n^{(i)}(x − O)} ,

we have to check whether the signs of the numerator and the denominator are different. First note that we have always x ≠ O, because in this case we would have found out that x lies in the origin simplex and left segment walk ( ) earlier. Note also that always λ ≠ 0 since λ = 0 would imply that O lies on a facet of the origin simplex, a contradiction.

If by chance the denominator is zero, the ray Ox^{− ; − ; →} intersects facet i in an infinite point. We can return false, because such a facet is never a candidate for the intersection with the largest λ < 0. ⟨ Member functions of class Triangulation 18 ⟩ +≡

bool Triangulation ::lambda negative(Simplex ∗S, number nx , number

nO , int i)

{ bool negative = false; // we assume that λ is positive if (nx ∗ number(dcur + 1) − nO < 0) negative = true;

// if the denominator is less then zero, inverse the sign else if (nx ∗ number(dcur + 1) − nO ≡ 0) return false;

// like discussed above if (alpha(S, i) ∗ number(dcur + 1) − nO > 0) return (negative);

// the sign does not change else return (¬negative); // the sign changes

} 45. The Dimension Jump.

If the point being inserted is a dimension jump, we have to add it to the set of vertices of every simplex of the extended triangulation \overline{Δ}(π_{i−1}) and for every simplex F of Δ(π_{i−1}), we have to add a new simplex S(F ∪ \{\overline{O}\}) whose base facet is the corresponding simplex of the old triangulation and whose peak is the anti origin. To do so, we visit all simplices of the old triangulation starting at the origin simplex and visiting all neighbors of a visited simplex recursively. This is done by the function dimension jump( ). ⟨ Further member declarations of class Triangulation 17 ⟩ +≡

void dimension jump(Simplex ∗S, list item x);

46. Before we do a dimension jump, we compute the new center of the extended origin simplex. ⟨ Dimension jump 46 ⟩ ≡

dcur {+ \atop +};

quasi center += x;

dimension jump(origin simplex , item x );

clear visited marks (origin simplex );

This code is used in section 24.

47. In this section we describe the function dimension jump( ). Before we do this, we give an example of a dimension jump in the the two dimensional case. The following figure shows a typical constellation of vertices before a dimension jump.

[figure: .; .; .; .; .; .; .; .; .; .; .; .; .; x_{1}; x_{2}; x_{3}; \overline{O}; \overline{O}]

Figure 7: We are in dimension 1

The origin simplex is conv(x_{1}, x_{2}). The point x_{3} is not a dimension jump, because it lies on the line through x_{1} and x_{2}. At this point we have 4 simplices: two unbounded ones to the left and to the right with the anti-origin as their peak, the origin simplex and the simplex conv(x_{2}, x_{3}) with peak x_{3} and the base facet x_{2}.

Now we do a dimension jump by inserting a point x_{4} not colinear with the other ones.

We have jumped to dimension 2. Now we have 6 simplices. For the simplices conv(x_{1}, x_{2}) and conv(x_{2}, x_{3}) we have added two unbounded simplices below

[figure: .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; \overline{O}; \overline{O}; \overline{O}; \overline{O}; x_{1}; x_{2}; x_{4}; x_{3}]

Figure 8: A dimension jump

having \overline{O} as peak. The origin simplex is now conv(x_{1}, x_{2}, x_{4}). The simplex conv(x_{2}, x_{3}, x_{4}) has base facet conv(x_{2}, x_{3}) and peak x_{4}. It is the neighbor of the simplex with the vertices x_{2}, x_{3} and \overline{O} opposite to the vertex \overline{O}.

We can divide the simplices of \overline{Δ}(π_{i}) into three classes:

• Bounded extended simplices: they result from bounded simplices of \overline{Δ}(π_{i−1}) by adding x to the set of vertices.

• Unbounded extended simplices: they result from unbounded simplices of \overline{Δ}(π_{i−1}) by adding x to the set of vertices.

• New simplices: they result from bounded simplices of \overline{Δ}(π_{i−1}) by adding \overline{O} to the set of vertices.

In this and the subsequent sections we will use the following notation. For a simplex F of \overline{Δ}(π_{i−1}) let v_{0}, …, v_{dcur − 1} be its vertices (v_{0} is the peak). Let S = S(F ∪{x}) denote the simplex resulting from extending F. If F is bounded, let S new = S(F ∪ \{\overline{O}\}) be the new simplex constructed for F. The simplices of \overline{Δ}(π_{i}) look as follows:

• A (bounded or unbounded) extended simplex S has the vertices v_{0}, …, v_{dcur − 1}, x. Since the peak of a simplex is defined to be the vertex inserted last that was not a dimension jump, the peak of S is the same as the peak of F. Thus we append x to the list of vertices, i.e., we write the appropriate entries at position dcur into the arrays vertices , neighbors and opposite vertices .

• A new simplex S new has the vertices \overline{O}, v_{0}, …, v_{dcur − 1}, where \overline{O} is the peak.

In the following description we will continue to make the distinction between a simplex F of \overline{Δ}(π_{i−1}) and the extended simplex S resulting from it. In the implementation, both correspond to the parameter S of dimension jump( ). So every occurrence of F in the following corresponds to S in the program.

dimension jump( ) works as follows. Starting at the origin simplex it visits all simplices of \overline{Δ}(π_{i−1}) using depth-first-search. When a simplex F is visited it is declared visited, x is added to its set of vertices (this turns F into S = S(F ∪ {x})), and if the simplex is bounded then a new unbounded simplex S new = S(F ∪\{\overline{O}\}) is created. Then all neighbors of F in \overline{Δ}(π_{i−1}) are visited recursively. (Note that only the neighbors F neighbors [0], …, F neighbors [dcur − 1] are inspected). Once all neighbors are visited we update the neighbor relation. There we distinguish cases according to whether the simplex is bounded or not. ⟨ Member functions of class Triangulation 18 ⟩ +≡

void Triangulation ::dimension jump(Simplex ∗S, list item x)

{ Simplex ∗S new ; S visited = true; S vertices [dcur ] = x; if (S vertices [0] ≠ anti origin) \{ // S is bounded iff peak ≠ \overline{O}

⟨ Add a new unbounded simplex 48 ⟩ } /∗ The neighbor opposite to x might not yet exist. We make a call of dimension jump( ) for all unvisited neighbors of S. ∗/ for (int k = 0; k ≤ dcur − 1; k{+ \atop +}) \{ // for all neighbors of F if (¬S neighbors [k] visited ) dimension jump(S neighbors [k], x);

}

if (S vertices [0] ≡ anti origin) {

⟨ Complete neighborhood information if F is unbounded 49 ⟩ } else {

⟨ Complete neighborhood information if F is bounded 50 ⟩ }

}

48. For every bounded simplex F of \overline{Δ}(π_{i−1}) we add a new simplex S new = S(F ∪ \{\overline{O}\}) with peak \overline{O}. It is the neighbor of the bounded extended simplex S = S(F ∪\{x\}) opposite to x, and \overline{O} is the vertex opposite to x. For all vertices v of S different from x the neighbor of S opposite to v is the simplex S(F′ ∪{x}) where F′ is the neighbor of F opposite to v. Thus no action is required in the algorithm. ⟨ Add a new unbounded simplex 48 ⟩ ≡

S new = S neighbors [dcur ] = new Simplex (dmax );

all simplices .append (S new );

S opposite vertices [dcur ] = 0; S new vertices [0] = anti origin; for (int k = 1; k ≤ dcur ; k{+ \atop +}) S new vertices [k] = S vertices [k − 1];

This code is used in section 47.

49. We discuss how to compute the neighbors of unbounded extended simplices. The neighbor of an unbounded extended simplex S = S(F ∪ {x}) opposite to x is the simplex T with vert(F) ⊂ vert(T) and x ∉ vert(T). Consider the neighbor F′ ∈ \overline{Δ}(π_{i−1}) of F opposite to \overline{O}. F′ is bounded. Hence we constructed a simplex S new′ with vert(S new′) = vert(F′) ∪ \{\overline{O}\}. Since vert(F) \ \{\overline{O}\} ⊂ vert(F′) we have vert(F) ⊂ vert(S new′). Furthermore x ∉ vert(S new′). Thus T = S new′ is the neighbor of S opposite to x. We reach T from F (or S, respectively) by first going to the 0-th neighbor (that is F′ or S′, respectively) and then going to the dcur -th neighbor of S′ which is S new′ = S(F′ ∪ {x}) = T. The vertex opposite to x with respect to S is the vertex w opposite to \overline{O} with respect to F. Note that if w is the i-th vertex of F′ then it is the (i + 1)-st vertex of S new′ since we have inserted the anti-origin in vertices [0].

As in the previous section (bounded extended simplex), the neighborhood information for vertices v ≠ x of S is the same as for F and hence there is nothing to do for them. ⟨ Complete neighborhood information if F is unbounded 49 ⟩ ≡

S neighbors [dcur ] = S neighbors [0] neighbors [dcur ];

S opposite vertices [dcur ] = S opposite vertices [0] + 1;

This code is used in section 47.

50. Let F be a bounded simplex of \overline{Δ}(π_{i−1}). It gives rise to the extended simplex S = S(F ∪ \{x\}) and the new simplex S new = S(F ∪ \{\overline{O}\}). The neighbors of S were already computed in Section 48. We still need to determine the neighbors of S new . In order to create the neighborhood information for a new simplex S new , we step through the neighbors of F.

To find the neighbor of S new opposite to v ≠ \overline{O} consider the neighbor F′ ∈ \overline{Δ}(π_{i−1}) of F opposite to v. If F′ is unbounded, the neighbor of S new opposite to v is S′ and the vertex opposite to v is x. If F′ is bounded, the neighbor of S new opposite to v is the simplex S new′ constructed for F′ and the vertex opposite to v remains the same as in F. Note that a pointer to S new′ has been added to the neighbors array of F′ at position dcur during a recursive or a previous call of dimension jump( ).

The neighbor of a new simplex S new opposite to \overline{O} is S. The vertex opposite to \overline{O} is x. Recall that the k-th vertex of S is the k + 1-st vertex of S′.

⟨ Complete neighborhood information if F is bounded 50 ⟩ ≡

for (int k = 0; k < dcur ; k{+ \atop +}) \{

if (S neighbors [k] vertices [0] ≡ anti origin) { // if F′ is unbounded

S new neighbors [k + 1] = S neighbors [k];

// the neighbor of S new opposite to v is S′ S new opposite vertices [k + 1] = dcur ; // x stands in position dcur

}

else { // F′ is bounded

S new neighbors [k + 1] = S neighbors [k] neighbors [dcur ];

// neighbor of S new opposite to v is S new′ S new opposite vertices [k + 1] = S opposite vertices [k] + 1; // … vertex opposite to v remains the same … // again remember the ‘shifting’ of the vertices one step to the right

} } /∗ the simplex opposite to \overline{O} with respect to S new is S, and the vertex is x ∗/ S new neighbors [0] = S; S new opposite vertices [0] = dcur ;

This code is used in section 47.

51. Output Routines.

In order to demonstrate our program, we now add to Triangulation a function show ( ), which draws (in the special case dmax ≡ 2) the simplicial complex into a LEDA-window. Running through the list all simplices we draw each simplex. For each simplex, we draw its vertices and for each vertex of a simplex we draw the edges connecting it to the other vertices of the simplex. Clearly we do not draw the anti-origin and the edges incident to it. Thus the for-loop which steps through all vertices starts with v = 0 if S is bounded (i.e., S vertices [0] ≠ anti origin) and with v = 1 if S is unbounded (i.e., S vertices [0] ≡ anti origin). Furthermore, we draw every point that we have inserted so far onto the screen (there may be many points that are not vertices of any simplex). We do this by running through the list coordinates . ⟨ Member functions of class Triangulation 18 ⟩ +≡

void Triangulation ::show (window &W)

{ /∗ We first draw every simplex ∗/ Simplex ∗S; forall (S, all simplices ) {

for (int v = (S vertices [0] ≡ anti origin ? 1 : 0); v ≤ dcur ; v{+ \atop +}) \{

// for each vertex except the anti-origin vector x = coordinates .contents (S vertices [v]); point a(x[0], x[1]); for (int e = v + 1; e ≤ dcur ; e{+ \atop +}) \{

// draw undrawn edges incident to vertex vector y = coordinates .contents (S vertices [e]); point b(y[0], y[1]); /∗ draw the edges of unbounded simplices as thick lines ∗/ if (S vertices [0] ≡ anti origin) W .set line width(3); else W .set line width(1); W .draw segment (a, b);

} } } /∗ Now we draw every point ∗/ vector x; forall (x, coordinates ) { point a(x[0], x[1]); W .draw point (a);

} }

52. print all ( ) prints information about all simplices to stdout (similar to show). This was useful for debugging. The information of a single simplex is printed by the function print ( ). ⟨ Member functions of class Triangulation 18 ⟩ +≡

void Triangulation ::print all ( ) {

Simplex ∗S; forall (S, all simplices ) print (S);

}

53. Here is a short function that prints the data of a simplex. ⟨ Member functions of class Triangulation 18 ⟩ +≡

void Triangulation ::print (Simplex ∗S) {

cout ≪ "\n[" ≪ S sim nr ≪

  "]-----------------------------------------------------\n";

for (int i = 0; i ≤ dcur ; i{+ \atop +}) \{ if (S vertices [i] ≡ anti origin) cout ≪ "anti"; else cout ≪ coordinates .contents (S vertices [i]); cout ≪ "␣["; if (S neighbors [i]) cout ≪ S neighbors [i] sim nr ; else cout ≪ "*"; cout ≪ "]␣"; if ((S vertices [0] ≠ anti origin ∨ i ≡ 0) ∧ dcur > 0) {

// cout ≪ ";␣normal:␣" ≪ normal (S, i); cout ≪ ";␣normal:␣" ≪ S normal values [i];

// cout ≪ ";␣␣alpha:␣" ≪ alpha(S, i); cout ≪ ";␣␣alpha:␣" ≪ S alpha values [i]; cout ≪ ";␣valid_equations:␣" ≪ S valid equations [i];

} cout ≪ endl ;

} cout .flush( );

} 54. The main function.

There are three ways to feed the data into the program: we can take the input from the keyboard, from a file or via mouse input from a graphics window (only if we work in dimension 2). If the input is taken from the keyboard or from a file, the first number must be an integer specifying the dimension of the following coordinate vectors. If the input is taken from a file, the second number in the file is read but ignored by our program (in order to be able to use input files that are created by the program rbox which generates random input files; it is a tool of the QHULL–system (cf. [1])). The remaining numbers in the file are taken as the coordinates of the points. We can call the program from a shell with the following command line arguments in an arbitrary order:

• m: read input from mouse.

• k: read input from keys, first entering the dimension we will work in, then the coordinates of the points. The input process stops with an end-of-file

     (ctrl-D).

• f: read input from a file whose name must be given as the next argument in the command line.

• p: print informations about all simplices after each insertion.

• n: no display: when working in dimension 2 only draw the final result.

• s: suppress any display when working in dimension 2

• V: use the visibility search method.

• M: use the modified visibility search method.

• S: use the segment walking method.

We first give a function that tells the user the correct usage of the command line arguments of the program when he makes a mistake when invoking the program.

In the command line, the user can give any number of the above arguments, but only the last ones are valid. ⟨ Main program 54 ⟩ ≡

void tell usage(string prg name)

// prg name is the name of the executable program {

cout ≪ "Usage:␣" ≪ prg name ≪

        "␣[␣m␣|␣k␣|␣f␣filename␣|␣p␣|␣n␣|␣s␣|␣V␣|␣M␣|␣S]*" ≪ endl ;

// a regular expression exit (1);

}

See also section 55. This code is used in section 6.

55. The main program first reads in the command line setting the options, then it processes the data. ⟨ Main program 54 ⟩ +≡

enum input method {

  MOUSE, KEYS, INPUTFILE

}; main(int argc, char ∗∗argv ) {

⟨ Read the command line 56 ⟩ ⟨ Process the data 57 ⟩;

}

56. In the command line, every option consists of a single character. ⟨ Read the command line 56 ⟩ ≡

string istreamargs (argc, argv );

// create an input stream from the command line string prg name; // the name of the compiled, executable program args ≫ prg name; // get the name from the command line string option; // the options we will take from the command line string data file = "/dev/null";

// the name of the file that contains the data; /∗ data file is initialized to "/dev/null" to avoid complicated special treatment when no input file is specified ∗/ int dimension; // the dimension we will work in int number of points ;

// appears in input files generated by rbox, not used by our program input method read from = MOUSE; // default: read from mouse search method m = SEGMENT_WALK; // default: segment walking bool draw all = true; // draw every insert step (if dimension ≡ 2) bool suppress = false; // suppress any display (if dimension ≡ 2) bool print simplices = false;

// print information about all simplices after an insert while (true) {

args ≫ option;

if (args .eof ( )) break; // as long as we have command line arguments if (option.length( ) ≠ 1)

// if the current argument has more than one character tell usage(prg name); // tell the correct usage of the program

switch (option[0]) { // which option is to be processed? case 'm': read from = MOUSE;

break; case 'k': read from = KEYS;

break; case 'f':

/∗ this argument must be followed by another argument which is taken as the name of a file from which we read the data ∗/ args ≫ data file; // get the filename if (args .eof ( )) // print error message if no file is specified

tell usage(prg name); read from = INPUTFILE; // we read from a file break;

case 'p': print simplices = true;

break; case 'n': draw all = false;

break; case 's': suppress = true;

break;

case 'V': m = VISIBILITY;

break;

case 'M': m = MODIFIED_VISIBILITY;

break;

case 'S': m = SEGMENT_WALK;

break; default: tell usage(prg name);

break; }

}

This code is used in section 55.

57. Here is how we process the data. ⟨ Process the data 57 ⟩ ≡

/∗ if the input is not taken from the mouse, we need a file from which we read the data ∗/ file istream file in(data file); // file istream is a LEDA type if (¬file in) {

cout ≪ "unable␣to␣open␣file␣" ≪ data file ≪ endl ;

exit (2); } switch (read from) {

case MOUSE:

{

⟨ Input from mouse 58 ⟩ } exit (0); break;

case KEYS: cout ≪ "Dimension␣of␣coordinate␣vectors:␣";

cin ≫ dimension;

break; case INPUTFILE: file in ≫ dimension;

file in ≫ number of points ; // we do not use this value break;

} ⟨ Input from keyboard or file 59 ⟩

This code is used in section 55.

58. We use LEDA’s window type to implement a graphical input tool. We are working with the X11R5 (xview) window system.

By a click of the left mouse button, we can input a new two dimensional point into the whole complex. Then the triangulation will be drawn onto the screen. The convex hull is represented by thick lines, whereas the other lines of the triangulation are drawn as thin lines. A click of the right mouse button ends the program. The input points are automatically logged to the file chull.pts. ⟨ Input from mouse 58 ⟩ ≡

window W; Triangulation T (2, m);

// we are working in the plane with search method m double a, b; // coordinates of a point in the window file ostream protocol ("chull.pts"); int mouse = 0; // variable to indicate which mouse button was pressed protocol ≪ 2 ≪ endl ; // write the dimension to chull .pts while (mouse ≠ 3) { // while mouse click is not the right button

mouse = W .read mouse(a, b);

// read the window coordinates into a and b if (mouse ≡ 1) { // left button pressed

vector x(2); // get a two dimensional vector x[0] = a; x[1] = b; protocol ≪ x ≪ endl ; T .insert (x); W .clear ( ); T .show (W); if (print simplices ) T .print all ( );

} } cout ≪ endl ≪ "Searched␣Simplices:␣" ≪ T .searched simplices ≪ endl ;

// only for statistical reasons

This code is used in section 57.

59. If we take the input from the keyboard or from a file, we read for each point its dimension coordinates and insert it. If dimension ≡ 2 we also open a graphics window in order to display the triangulation. The window remains on the screen until the right mouse button is pressed in it. Keyboard input is terminated by an end-of-file (ctrl-D). ⟨ Input from keyboard or file 59 ⟩ ≡

Triangulation T (dimension, m); vector x(dimension); if (dimension ≡ 2 ∧ ¬suppress ) {

window W; while (¬(file in.eof ( ) ∨ cin.eof ( ))) {

if (read from ≡ KEYS) cin ≫ x; else file in ≫ x; T .insert (x); if (¬suppress ∧ draw all ) {

W .clear ( );

T .show (W); } if (print simplices ) T .print all ( );

} W .clear ( ); T .show (W);

cout ≪ "Press␣the␣right␣button␣in␣the␣drawing␣w\
  indow␣to␣terminate.\n";

cout .flush( ); while (W .read mouse( ) ≠ 3) ;

} else {

while (¬(file in.eof ( ) ∨ cin.eof ( ))) { if (read from ≡ KEYS) cin ≫ x; else file in ≫ x; T .insert (x); if (print simplices ) T .print all ( );

} } cout ≪ endl ≪ "Searched␣Simplices:␣" ≪ T .searched simplices ≪ endl ; cout ≪ "Simplices␣created:␣" ≪ T .created simplices ( ) ≪ endl ;

// only for statistical reasons

This code is used in section 57.

60. Solving a system of linear equations.

In this part, we add to the LEDA-type matrix a function which solves a system of linear equations in the following sense: given a matrix A and a vector b, the function applies the Gaussian elimination algorithm for solving a linear system and returns a list of vectors which characterize the affine-linear space of the solution. When the list returned is empty, there is no solution at all. Otherwise, the list has the form (a, d_{1}, …, d_{r}), which means that d_{1}, …, d_{r} are the spanning vectors of the affine-space and a is a base point which lies in that space.

In the subsequent part we give two applications of the algorithm which we need: determine whether a set of vectors is linearly dependent or affine-linearly dependent.

61. The program has the following structure:

⟨ linalg.c 61 ⟩ ≡

⟨ Header files for linalg 63 ⟩ ⟨ The function 65 ⟩ ⟨ Dependency tests 74 ⟩

62. In the actual version of LEDA, the entries of a matrix are of type double, but in later versions the user will be able to choose between several kinds of arithmetical operation, for example rational arithmetics. Therefore, we define here a macro number and use it in the algorithm defined as double or rational, but it can be easily changed to any other type providing the operations =,+,–,*,/ and an absolute value function like fabs ( ) for C++’s double type. (If you want to compute in a prime field for example, you can forget about the fabs ( ) function. Alternatively to fabs ( ) you can also use the <-operator.) The only lines of code that must be changed are the ones in which we update the pivot-element (see below) and test it against zero. Because of the low stability of C++’s type double, we never test a double-variable to be zero. We concern it as zero, if its absolute value is less than a constant EPSILON.

We have also implemented some routines which compute with rational arithmetics and thus guarantee to be exact in every case. To choose this kind of arithmetics one simply has to define the RATIONAL macro. The macro definitions depend on whether we compute with rational arithmetics or not. ⟨ Macro definitions 62 ⟩ ≡

#ifdef RATIONAL

#define number rational #define vector rat vector #define matrix rat matrix

#define EPSILON 0

#define number abs abs

// for rational, the function is called “abs ( )” #else #define number double #define EPSILON 1 · 10^{−10}

#ifdef USE_FABS

#define number abs fabs #else #define number abs(a) ((a) < 0 ? −(a) : (a)) #endif #endif

This code is used in section 64.

63. We include the following LEDA header files, depending on the kind of arithmetics we use: ⟨ Header files for linalg 63 ⟩ ≡

#include <LEDA/list.h>
#include <LEDA/array.h>
#include <stream.h>
#include <math.h>
#ifdef RATIONAL
#include "rat_matrix.h"

#else

#include <LEDA/matrix.h>

#endif

#include "linalg.h"

This code is used in section 61.

64. We create the following header file for the functions we will write. The main work of linear solver ( ) is done in the function raw linear solver ( ) which gets a pointer C to a 2-dimensional C-array in which the calculation is done Furthermore, it gets as parameters rows and cols the dimensions of the matrix A of the underlaying system of equations. (I.e., since C also contains the right side b, C is an array of dimension rows × (cols + 1).). We can save some LEDA overhead when we call raw linear solver ( ) directly, avoiding the initialization of the LEDA-matrix A. We make use of this in plane equation( ).

⟨ linalg.h 64 ⟩ ≡

⟨ Macro definitions 62 ⟩

list⟨vector⟩ linear solver (matrix &A, const vector &b); list⟨vector⟩ raw linear solver (number ∗∗C, int rows , int cols ); bool linear dependency (list⟨vector⟩ L); bool affine dependency (list⟨vector⟩ L); 65. In the function linear solver ( ) we only have to initialize the matrix in which raw linear solver ( ) calculates and call raw linear solver ( ). We only want to allocate new memory for C if the dimensions have changed with respect to the last call. Thus we make C, rows and cols static. ⟨ The function 65 ⟩ ≡

list⟨vector⟩ linear solver (matrix &A, const vector &b) {

static number ∗∗C;

// the matrix in which we will calculate (C = (A|b)) static int rows ; // the number of rows of A static int cols ; // the number of columns of A int i, j; /∗ reallocate memory for C if necessary ∗/ if (rows ≠ A.dim1 ( ) ∨ cols ≠ A.dim2 ( )) {

if (C) { for (i = 0; i < rows ; i{+ \atop +}) delete C[i]; // delete each row delete C;

}

rows = A.dim1 ( );

cols = A.dim2 ( );

C = new number ∗ [rows ];

for (i = 0; i < rows ; i{+ \atop +}) C[i] = new number [cols + 1]; } /∗ copy A and b into C ∗/ for (i = 0; i < rows ; i{+ \atop +}) \{

for (j = 0; j < cols ; j{+ \atop +}) C[i][j] = A(i, j);

C[i][cols ] = b[i]; } return raw linear solver (C, rows , cols );

}

See also section 66. This code is used in section 61.

66. We can now give an overview of raw linear solver ( ). ⟨ The function 65 ⟩ +≡

list⟨vector⟩ raw linear solver (number ∗∗C, int rows , int cols ) {

list⟨vector⟩ L; // the list that is returned int i, j, k; // indices to step through the matrix ⟨ Make upper diagonal form 67 ⟩ ⟨ Make unity matrix on the left side 71 ⟩ ⟨ Compute the list L which contains the solution 72 ⟩

⟨ Deallocate memory and return the list L 73 ⟩ }

67. This part of the function changes C into upper diagonal form. Here is the description of the first step. We choose the pivot-element p (the element C_{i,j} with which we make appropriate columns to zero) the element of highest absolute value in C. Note that the last column of C is b and that we don’t search for p in this column. Note also that this rule guarantees maximal arithmetical stability. Having found p in row i and column j, we interchange row i with row 1 and column j with column 1, so that p stands now in C_{0,0}. In the array var , we store the indices of the variables for every column. Then we use the pivot-element p to set C_{0,1}, C_{0,2}, … = 0. The second step is like the first one, but it only works on the submatrix starting at C_{1,1}. If in this process p is found zero, we have upper diagonal form. ⟨ Make upper diagonal form 67 ⟩ ≡

#ifdef USE_LEDA_ARRAYS

array⟨int⟩ var (0, cols − 1); #else

int ∗var = new int [cols ]; #endif

number p;

for (i = 0; i < cols ; i{+ \atop +}) var [i] = i;

// at the beginning, variable x_{i} stands in column i /∗ here comes the main loop ∗/ for (i = 0; i < cols ; i{+ \atop +}) {

p = −1; // initialize p to a negative value

⟨ Search the pivot-element, interchange rows and columns 68 ⟩ ⟨ Do the pivoting 69 ⟩

}

⟨ Test whether the system has a solution 70 ⟩

This code is used in section 66.

68. We find the pivot-element by searching through the actual submatrix. Since the interior of this loop is executed very often, we use some dirty C constructions. ⟨ Search the pivot-element, interchange rows and columns 68 ⟩ ≡

int p row , p col ; // position of p

number ∗N; // to step through the elements of a row of C for (j = i; j < rows ; j{+ \atop +}) \{ // step through rows i to rows − 1

N = &(C[j][cols ]); // in each row start with N at the right end for (k = cols − i + 1; k^{−}; ) {

// step through columns cols − i + 1 + i − 1 to i − 1 + 1 /∗ N still points to the element right of the current column ∗/ if (∗^{−}N > p ∨ ∗N < −p) { // new pivot element found

p = number abs(∗N);

p row = j; // store the position of the new p p col = k + i − 1;

} }

} /∗ if p is zero, we can stop. The test whether p is zero actually is a test whether p is less then EPSILON ∗/ if (p ≤ number(EPSILON)) break; // exit the main loop /∗ We interchange rows i and p row by interchanging the pointers ∗/ number ∗help = C[i]; C[i] = C[p row ]; C[p row ] = help; /∗ We interchange columns i and p col by copying ∗/ if (p col > i)

for (j = 0; j < rows ; j{+ \atop +}) \{ number swap = C[j][i]; C[j][i] = C[j][p col ]; C[j][p col ] = swap;

}

/∗ We store the interchanging of the variables in var ∗/ int dummy = var [i]; var [i] = var [p col ]; var [p col ] = dummy ;

This code is used in section 67.

69. Now we are ready to do the pivot-step with the element C_{i,i}. Again we do a bit dirty C stuff to save a little time. ⟨ Do the pivoting 69 ⟩ ≡

p = C[i][i]; // note that p was always positive ! C[i][i] = 1; // C_{i,i} becomes 1 for (k = i + 1; k ≤ cols ; k{+ \atop +}) // treat row i

C[i][k] = C[i][k]/p; number ∗Cji ; // to step through C[j][i] number ∗Cii = &(C[i][i]); // will remain uncanged for (j = i + 1; j < rows ; j{+ \atop +}) \{ // for each row below row i

Cji = &(C[j][i]);

number factor = ∗Cji ;

k = cols − i;

while (k^{−}) // for each C[j][i + k] with i + k ∈ [i … cols − 1] Cji [k] −= factor ∗ Cii [k]; // note that Cji [k] ≡ C[j][i + k]

}

This code is used in section 67.

70. We introduce a new variable max row , which holds the highest index of a non-zero row created by the elimination algorithm. If we have stopped with p = 0, we must test whether the system has a solution at all: if there is now in the rightmost column for b in a row whose index is greater then max row a value not equal to zero, then the system has no solution at all. ⟨ Test whether the system has a solution 70 ⟩ ≡

int max row = i − 1; // we have gone one row too far in the loop if (p ≤ EPSILON) // we have stopped with p = 0

for (j = max row + 1; j < rows ; j{+ \atop +})

if (number abs(C[j][cols ]) > EPSILON) return L;

// we return the empty list L

This code is used in section 67.

71. We now have the diagonal elements C_{0,0}, …, C_{max row ,max row} = 1. All elements below them are zero. We make all elements above them zero, too. Then we have a unity matrix E_{max row+1} in the upper left corner. ⟨ Make unity matrix on the left side 71 ⟩ ≡

for (i = max row ; i ≥ 0; i^{−}) { // go up in the matrix number ∗Cji ; number ∗Cii = &(C[i][i]); // will remain uncanged for (j = i − 1; j ≥ 0; j^{−}) { // for each row above row i

Cji = &(C[j][i]);

number factor = ∗Cji ;

k = cols − i;

while (k^{−}) { // for each C[j][i + k] with i + k ∈ [i … cols − 1] Cji [k] −= factor ∗ Cii [k]; // note that Cji [k] ≡ C[j][i + k]

} }

}

This code is used in section 66.

72. We now have E_{max row +1} in the upper left corner, which means that we have max row + 1 variables, which depend on the cols − (max row + 1) free variables. The vector var tells us, which variables are free and which are dependent. The affine space of the solution is thus generated by cols − (max row + 1) spanning vectors and has a base point base. For every free variable, we create a vector. Then we give for every dependent vector its depending parts to the appropriate spanning vectors and to base. The dimension of the spanned space is stored in the variable dim of span. ⟨ Compute the list L which contains the solution 72 ⟩ ≡

int dim of span = cols − (max row + 1); /∗ dim of span > 0 means that the solution contains more than one point ∗/

#ifdef USE_LEDA_ARRAYS

array⟨vector⟩ span(0, (dim of span > 0) ? dim of span − 1 : 0); #else

vector ∗span = new vector [(dim of span > 0) ? dim of span : 1]; #endif

vector base(cols ); // base is initialized to be the 0-vector for (i = 0; i < dim of span; i{+ \atop +}) span[i] = base;

// initialize span[i] to be the 0-vector /∗ Next, we set the components of those spanning vectors to 1, which belong to the free variables ∗/ for (i = 0; i < dim of span; i{+ \atop +}) span[i][var [max row + 1 + i]] = 1; /∗ Now we treat every dependent variable ∗/ for (i = 0; i ≤ max row ; i{+ \atop +}) \{

base[var [i]] = C[i][cols ];

for (k = max row + 1; k < cols ; k{+ \atop +})

span[k − max row − 1][var [i]] = −C[i][k]; }

This code is used in section 66.

73. After we have deallocated the space for C, we can return the list L. ⟨ Deallocate memory and return the list L 73 ⟩ ≡

for (i = 0; i ≤ cols − max row − 2; i{+ \atop +}) L.push(span[i]); /∗ O.K., let’s insert base at the beginning of L and we are done ! ∗/ L.push(base);

#ifndef USE_LEDA_ARRAYS

delete var ; delete span;

#endif

return L;

This code is used in section 66.

74. Dependency tests.

We first give an algorithm which determines whether a given list of vectors v_{1}, …, v_{r} are linearly dependent. To test them, we ask if there are α_{1}, …, α_{r}, not all equal to zero, such that ∑_{i=1}^{r} α_{i}v_{i} = 0. ⟨ Dependency tests 74 ⟩ ≡

bool linear dependency (list⟨vector⟩ L) /∗ returns true iff the vectors in L are linearly dependent. ∗/ { /∗ we simply create the appropriate matrix and give it to linear solver ( ), the left side b being zero ∗/

int cols = L.size( ); // the number of vectors int rows = L.head ( ).dim( ); // the dimension we are working in if (cols > rows ) return true;

// more then rows vectors are always dependent matrix A(rows , cols ); vector b(rows ); // note that b is by default the zero vector vector v; int row ; int col = 0; forall (v, L) { // fill the columns of A with the vectors in L

for (row = 0; row < rows ; row {+ \atop +}) A(row , col ) = v[row ]; col {+ \atop +};

} list⟨vector⟩ R; // result list for the call of linear solver ( ) R = linear solver (A, b); if (R.size( ) > 1) return true;

// they are surely dependent because we have infinitely many // (α_{1}, …, αcols ) that solve the system above

/∗ Otherwise, we check whether any α_{i} is not equal to 0 ∗/ vector test = R.pop( ); for (row = 0; row < cols ; row {+ \atop +}) \{

if (number abs(test [row ]) > EPSILON) return true; } /∗ The only case left now is that only the trivial linear combination of the v_{i} is zero, that is, they are linearly independent ∗/ return false;

}

See also section 75. This code is used in section 61.

75. The next function tests whether the vectors v_{1}, …, v_{r} are affine-linearly dependent. Note that this is the case iff v_{2} − v_{1}, …, v_{r} − v_{1} are linearly dependent. Therefore, we can simply use the function linear dependency ( ) for the test. ⟨ Dependency tests 74 ⟩ +≡

bool affine dependency (list⟨vector⟩ L) /∗ returns true iff the vectors in L are affine-linearly dependent. ∗/ {

list⟨vector⟩ Help; vector v1 ; v1 = L.pop( ); while (¬L.empty ( )) Help.push(L.pop( ) − v1 ); return linear dependency (Help);

} 76. Useful literature.

References

[1] C. Bradford Barber, David P. Dobkin, Hannu Huhdanpaa “The quickhull Algorithm for Convex Hull” Geometry Center Techical Report GCG53, University of Minnesota, available via anonymous ftp from

geom.umn.edu:pub/qhull.tar.Z

[2] Christoph Burnikel, Kurt Mehlhorn, Stefan Schirra “On Degeneracy in Geometric Computations”

[3] K.E. Clarkson, K. Mehlhorn, Raimund Seidel “Four Results on Randomized Incremental Constructions” Technical Report MPI-I-92-112, March 1992

[4] K.E. Clarkson, K. Mehlhorn, Raimund Seidel “Four Results on Randomized Incremental Constructions” revised version of February 15, 1993

[5] Donald Knuth, Silvio Levy “The CWEB System of Structured Documentation” available via anonymous ftp from labrea.stanford.edu:

/pub/cweb/cweb.tar.gz  or  ftp.th-darmstadt.de:
/pub/programming/literate-programming/c.c++/cweb.tar.gz

[6] Kurt Mehlhorn, Stefan Näher: “LEDA, a Library of Efficient Data Types and Algorithms”, Proceedings of the 14th Symposium on Mathematical Foundations of Computer Science, LNCS Vol. 379, 88-106, 1989, (to appear in Communications of the ACM)

[7] Stefan Näher: “LEDA Manual Version 3.0”, Technical Report, MPI-I-93-106, Max-Planck-Institut für Informatik, Saarbrücken, 1993.

Index

A: \underline{9}, \underline{60}, \underline{64}, \underline{65}, \underline{74}. a: \underline{51}, \underline{58}. abs : \underline{62}. Af : \underline{30}. affine dependency : 9, 26, \underline{64}, \underline{75}. all simplices : \underline{11}, 14, 27, 29, 48,

51, 52. alpha: 15, \underline{17}, 18, 20, 23, 29, 33,

34, 35, 37, 38, 42, 44, 53. alpha values : \underline{15}, 16, 18, 19,

23, 53. anti-origin: 4. anti origin: \underline{15}, 23, 27, 29, 32,

33, 34, 35, 38, 45, 47, 48, 50, 51, 53.

append : 24, 27, 29, 48. argc: \underline{55}, 56. args : 56. argv : \underline{55}, 56. b: \underline{9}, \underline{51}, \underline{58}, \underline{60}, \underline{64}, \underline{65}, \underline{74}. base: 20, \underline{21}, 23, \underline{72}, 73. base facet: 3. bool: \underline{26}, \underline{42}, \underline{44}. bounded: 4. C: \underline{64}, \underline{65}, \underline{66}. chull : 58. Cii : \underline{69}, \underline{71}. cin: 57, 59. Cji : \underline{69}, \underline{71}. clear : 28, 58, 59. clear visited marks : \underline{31}, 32, 36, 46. col : \underline{21}, \underline{74}. collect outer simplices : 31. collect visible simplices : \underline{31}, 32,

35. cols : 20, \underline{21}, \underline{64}, \underline{65}, \underline{66}, 67, 68,

69, 70, 71, 72, 73, \underline{74}. contents : 21, 23, 26, 51, 53. coordinates : \underline{11}, 15, 21, 23, 24,

26, 51, 53. cout : 53, 54, 57, 58, 59. created simplices : \underline{11}, 12, 59. d: \underline{11}, \underline{13}. data file: \underline{56}, 57. dcur : \underline{11}, 13, 15, 18, 19, 20, 21, 22,

23, 24, 26, 27, 29, 30, 33, 34, 35, 36, 37, 38, 40, 41, 42, 43, 44, 46, 47, 48, 49, 50, 51, 53.

diff : \underline{42}. difference: \underline{21}. dim: 74. dim of span: \underline{72}. dimension: \underline{56}, 57, 59. dimension jump: 3. dimension jump: \underline{45}, 46, 47,

48, 50. dim1 : 65. dim2 : 65. dir : \underline{17}, \underline{19}, 20, 21, 22, 23. dmax : \underline{11}, 13, \underline{15}, \underline{16}, 20, 21, 22,

24, 27, 29, 48, 51. draw all : \underline{56}, 59. draw point : 51. draw segment : 51. dummy : \underline{68}. Dummy : \underline{12}. dummys created : \underline{12}. e: \underline{51}. empty : 28, 75. endl : 53, 54, 57, 58, 59. eof : 56, 59. EPSILON: \underline{62}, 68, 70, 74. exit : 54, 57. fabs : 9, \underline{62}. factor : \underline{69}, \underline{71}. false: 16, 36, 38, 42, 43, 44, 56, 74. file in: \underline{57}, 59. flush: 53, 59. head : 74. Help: \underline{75}. help: \underline{22}, \underline{68}. horizon ridge: 5. i: \underline{16}, \underline{17}, \underline{18}, \underline{21}, \underline{26}, \underline{33}, \underline{34}, \underline{35},

\underline{36}, \underline{37}, \underline{39}, \underline{42}, \underline{44}, \underline{53}, \underline{65}, \underline{66}. ii : \underline{29}. in: \underline{37}, 40, 41. input method: \underline{55}. INPUTFILE: 55, 56, 57. insert : \underline{11}, 24, 58, 59. int: \underline{12}. is dimension jump: 24, \underline{25}, 26. item x : \underline{24}, 27, 28, 29, 30, 46. j: \underline{22}, \underline{39}, \underline{42}, \underline{65}, \underline{66}. k: \underline{29}, \underline{30}, \underline{42}, \underline{43}, \underline{47}, \underline{48}, \underline{50}, \underline{66}. KEYS: 55, 56, 57, 59. L: \underline{2}, \underline{20}, \underline{26}, \underline{64}, \underline{66}, \underline{74}, \underline{75}. lambda cmp: 37, \underline{39}, 40, 41, 42, 44. lambda negative: \underline{39}, 40, 44. last simplex : \underline{32}.

LEDA_MEMORY: 15.

length: 56. lfdnr : \underline{16}. linear dependency : \underline{64}, \underline{74}, 75. linear solver : 9, \underline{64}, \underline{65}, 74. M: \underline{21}. m: \underline{11}, \underline{13}, \underline{56}. main: \underline{55}. matrix: \underline{62}. max row : \underline{70}, 71, 72, 73. method : \underline{11}, 13, 32.

MODIFIED_VISIBILITY: 11, 32,

56. mouse: \underline{58}. MOUSE: 55, 56, 57. N: \underline{68}. negative: \underline{44}. neighbors : \underline{15}, 16, 19, 27, 29, 30,

33, 34, 35, 36, 41, 47, 48, 49, 50, 53.

new y1 : \underline{30}. NewSimplices : \underline{28}, 29, 30. nil : 13, 15, 16, 31, 32, 34. niO : \underline{39}, \underline{42}, 43. nix : \underline{39}, \underline{42}.

nix niO : \underline{42}, 43.
njO : \underline{39}, \underline{42}, 43.
njx : \underline{39}, \underline{42}.
njx njO : \underline{42}, 43.
nO : \underline{37}, 38, \underline{39}, 40, 41, \underline{44}.
normal : 15, \underline{17}, 18, 20, 21, 22,

23, 29, 33, 34, 35, 37, 38, 42, 43, 53.

normal values : \underline{15}, 16, 18, 19,

22, 23, 53. number: \underline{18}, \underline{62}.

number abs: \underline{62}.
number of points : \underline{56}, 57.
nx : \underline{37}, 38, \underline{39}, 40, 41, \underline{44}.

opposite vertices : \underline{15}, 16, 19, 27, 29, 30, 41, 47, 48, 49, 50.

option: \underline{56}.

origin simplex: 3.

origin simplex : \underline{11}, 13, 21, 26,

27, 32, 37, 46. out : 37, \underline{41}. outer simplex : \underline{27}. p: \underline{67}. p col : \underline{68}. p row : \underline{68}. peak: 3. plane equation: \underline{17}, 18, 19, 64. pop: 22, 74, 75. prg name: \underline{54}, \underline{56}. print : \underline{11}, 16, 52, 53. print all : \underline{11}, 52, 58, 59. print simplices : \underline{56}, 58, 59. protocol : \underline{58}. pts : 58. push: 26, 33, 35, 73, 75. quasi center : \underline{11}, 20, 23, 27, 37,

38, 46. R: \underline{74}. RATIONAL: 9, 10, 62, 63. raw linear solver : 20, \underline{64}, 65, \underline{66}. read from: \underline{56}, 57, 59. read mouse: 58, 59. result : \underline{34}. row : \underline{21}, \underline{74}. rows : 20, \underline{21}, \underline{64}, \underline{65}, \underline{66}, 68,

69, 70, \underline{74}. S: \underline{14}, \underline{17}, \underline{18}, \underline{19}, \underline{28}, \underline{31}, \underline{33}, \underline{34},

\underline{35}, \underline{36}, \underline{37}, \underline{39}, \underline{42}, \underline{44}, \underline{45}, \underline{47}, \underline{51}, \underline{52}, \underline{53}.

S new : \underline{47}, 48, 49, 50. search method: \underline{11}, 13, 56. search to outside: \underline{31}, 32, 34. searched simplices : \underline{11}, 13, 33, 34,

35, 37, 58, 59.

SEGMENT_WALK: 11, 32, 56.

segment walk : \underline{31}, 32, 37, 44. set line width: 51. show : \underline{11}, 51, 58, 59. sigma: \underline{42}, 43. sim nr : 12, \underline{15}, 16, 53. Simplex: \underline{15}, \underline{16}. size: 74. span: \underline{72}, 73. stdout : 11, 52. string istream: 56. suppress : \underline{56}, 59. swap: \underline{68}. T: \underline{29}, \underline{30}, \underline{33}, \underline{34}, \underline{35}, \underline{58}, \underline{59}. tell usage: \underline{54}, 56. test : \underline{74}. Triangulation: \underline{11}, \underline{13}, \underline{14}, \underline{15},

\underline{18}, \underline{34}, \underline{37}. true: 33, 34, 35, 37, 38, 42, 43,

44, 47, 48, 56, 74. unbounded: 4.

USE_FABS: 62.
USE_LEDA_ARRAYS: 15, 16, 37,

67, 72, 73. v: \underline{51}, \underline{74}. valid equations : \underline{15}, 16, 18, 19, 53. var : \underline{67}, 68, 72, 73. vector: \underline{62}. vertices : \underline{15}, \underline{16}, 21, 22, 23, 26, 27,

28, 29, 30, 32, 33, 34, 35, 38, 47, 48, 49, 50, 51, 53.

VISIBILITY: 11, 32, 56.

visibility search: \underline{31}, 32, 33, 36. visible from x: 3.

visible simplices : \underline{11}, 28, 33, 35.
visited : \underline{15}, 16, 31, 33, 34, 35,

36, 47. void: \underline{19}, \underline{24}, \underline{33}, \underline{35}, \underline{36}, \underline{47},

\underline{51}, \underline{52}, \underline{53}.

v1 : \underline{75}.
W: \underline{11}, \underline{51}, \underline{58}, \underline{59}.

x: \underline{11}, \underline{24}, \underline{25}, \underline{26}, \underline{31}, \underline{33}, \underline{34}, \underline{35}, \underline{37}, \underline{45}, \underline{47}, \underline{51}, \underline{58}, \underline{59}.

x-visible: 3.
x in S : \underline{37}, 38.
y: \underline{51}.
y1 : \underline{30}.
y2 : \underline{30}.

List of Refinements

⟨ Add a new unbounded simplex 48 ⟩ Used in section 47. ⟨ Complete neighborhood information if F is bounded 50 ⟩ Used in section 47. ⟨ Complete neighborhood information if F is unbounded 49 ⟩ Used in section 47. ⟨ Compute the arrays nx and nO and test whether x ∈ S 38 ⟩ Used in section 37. ⟨ Compute the list L which contains the solution 72 ⟩ Used in section 66. ⟨ Compute the plane equation 20 ⟩ Used in section 19. ⟨ Compute the right hand side alpha and adjust the sign of normal 23 ⟩ Used

in section 20.

⟨ Deallocate memory and return the list L 73 ⟩ Used in section 66. ⟨ Decide whether <^{σ} is < or > 43 ⟩ Used in section 42.

⟨ Dependency tests 74, 75 ⟩ Used in section 61. ⟨ Dimension jump 46 ⟩ Used in section 24. ⟨ Do the pivoting 69 ⟩ Used in section 67. ⟨ Extract normal from L 22 ⟩ Used in section 20.

⟨ Find x-visible hull facets 32 ⟩ Used in section 28. ⟨ Find the facet with largest λ < 0 40 ⟩ Used in section 37. ⟨ For each horizon ridge add the new simplex 29 ⟩ Used in section 28. ⟨ Further member declarations of class Triangulation 17, 25, 31, 39, 45 ⟩ Used

in section 11.

⟨ Go to the next Simplex on the ray Ox^{− ; − ; →} 41 ⟩ Used in section 37. ⟨ Header files for linalg 63 ⟩ Used in section 61. ⟨ Header files to be included 7, 8, 10 ⟩ Used in section 6. ⟨ Initialize the triangulation 27 ⟩ Used in section 24. ⟨ Input from keyboard or file 59 ⟩ Used in section 57.

⟨ Input from mouse 58 ⟩ Used in section 57. ⟨ Macro definitions 62 ⟩ Used in section 64. ⟨ Main program 54, 55 ⟩ Used in section 6.

⟨ Make unity matrix on the left side 71 ⟩ Used in section 66. ⟨ Make upper diagonal form 67 ⟩ Used in section 66. ⟨ Member functions of class Triangulation 18, 19, 24, 26, 33, 34, 35, 36, 37, 42, 44, 47,

51, 52, 53 ⟩ Used in section 6.

⟨ Non-dimension jump 28 ⟩ Used in section 24.

⟨ Process the data 57 ⟩ Used in section 55.

⟨ Read the command line 56 ⟩ Used in section 55. ⟨ Search the pivot-element, interchange rows and columns 68 ⟩ Used in section 67.

⟨ Set up the matrix 21 ⟩ Used in section 20.

⟨ Test whether the system has a solution 70 ⟩ Used in section 67.

⟨ The function 65, 66 ⟩ Used in section 61.

⟨ Update the neighborhood relationship 30 ⟩ Used in section 28.

⟨ class Simplex 15, 16 ⟩ Used in section 6. ⟨ class Triangulation 11, 12, 13, 14 ⟩ Used in section 6.

⟨ linalg.c 61 ⟩

LIST OF REFINEMENTS 63

⟨ linalg.h 64 ⟩