Subversion Repositories gelsvn

Rev

Rev 405 | Rev 596 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
405 jab 1
\documentclass[a4paper]{article}
2
\usepackage{listings}
3
\usepackage[isolatin]{inputenc}
4
\usepackage{makeidx,times}
5
\usepackage{fancyhdr}
6
\usepackage{graphicx}
7
\usepackage{float}
8
\usepackage{alltt}
9
%%\usepackage{doxygen}
10
\usepackage[pdftex,
11
            pagebackref=true,
12
            colorlinks=true,
13
            linkcolor=blue
14
           ]{hyperref}
15
\makeindex
16
\lstset{tabsize=1,language=C++,basicstyle=\small}
17
\author{J. Andreas B{\ae}rentzen}
465 jab 18
\title{Introduction to GEL\\\vspace{0.25cm} \normalsize
19
\textit{An informal introduction to the GEL library with many examples of usage}.}
405 jab 20
\setcounter{tocdepth}{2}
21
\begin{document}
22
\maketitle
23
\tableofcontents
465 jab 24
\newpage\section{Introduction}
405 jab 25
 
26
GEL is an abbreviation of GEometry and Linear algebra. GEL is a C++ library with tools for manipulating digital representations of geometry such as polygonal meshes, triangle meshes, point clouds, and voxel grids. Since linear algebra tools are invariably needed for this, GEL also contains tools for linear algebra. Finally, GEL contains tools for drawing geometry using OpenGL and various utilities.
27
 
28
\subsection{GEL Homepage}
29
The online home of GEL is here:
30
\href{http://www.imm.dtu.dk/GEL/}{http://www.imm.dtu.dk/GEL/}
31
 
32
\subsection{Purpose of this document}
465 jab 33
GEL is not sufficiently documented, but a fair amount of Doxygen documentation is provided for those parts of GEL which are used a lot. The goal of this document is to provide more detail and to help people getting started. In particular, this document focuses on functionality in the CGLA and HMesh namespaces, but we shall also take a look at some other parts of GEL. Unfortunately, it is rather difficult to be complete: There is functionality in GEL not covered here. Consequently, if you are looking for something particular, we \textbf{entreat }you to look at the Doxygen web (start at the GEL homepage). Most functions and classes are in fact described in the auto-generated documentation.
405 jab 34
 
35
\subsection{Overview of GEL}
36
 
37
GEL contains many tools for geometry processing and a number of data structures. Perhaps the most important is the Manifold class which resides in the HMesh namespace. A Manifold is a halfedge based representation of a polygonal mesh. There is also a simpler triangle mesh data structure, TriMesh, based on an indexed face set. There are several voxel grid data structures, a fairly efficient k-D tree for storing points, bounding box hierarchies etc.
38
 
465 jab 39
A number of algorithms exist for manipulating these representations of geometric data. For instance, Manifold has a function for mesh simplification, several functions for smoothing (including anisotropic), mesh optimization, and mesh clean up.
405 jab 40
Another strong point of GEL is the fact that it contains good support for converting between representations; especially between meshes and voxel grids. There are several tools for polygonizing isosurfaces and also converting meshes to distance fields.  The tools pertaining to the Manifold data structure are in the HMesh namespace. The other geometry related tools are in the Geometry namespace.
41
 
42
There are two packages for linear algebra: CGLA is strictly for small vectors and matrices (dimensions 2,3, and 4). LinAlg is a Lapack wrapper for slightly larger problems. At this point, it is fairly simple but includes a number of functions for solving linear systems, factoring matrices and finding eigensolutions for symmetric matrices.
43
 
44
GLGraphics contains facilities for drawing entities from other parts of GEL via OpenGL. Especially TriMesh and Manifold have good support. For instance, there is a nice wireframe drawing class. There are also two virtual trackballs for interactive programs and some tools for viewing in interactive programs and SOIL a small open source library for image loading by Jonathan Dummer.
45
 
46
Finally there are some utilities in Util: Getting command line arguments, hash table classes, a 2D grid class, a resource manager, an XML parser by Jeppe Frisvad and other tools.
47
 
48
Apart from the namespaces (sublibraries) we have a number of demo applications. Some of these might be useful in themselves. One example is a tool for converting a mesh to a signed distance field and an OBJ viewer  which is also able to load X3D and PLY meshes. We are working on the MeshEdit application which is already a highly useful swiss knife for geometry processing.
49
 
50
\section{Compilation and Installation}
51
 
52
Life would be easier if one's library did not rely on dependencies. We try hard to have only a minimal set, but the following libraries \textit{are} required in order to build GEL. 
53
\begin{itemize}
54
\item OpenGL
55
\item GLUT
56
\item GLEW
57
\item Lapack (and its depencies)
58
\item GLConsole (only for MeshEdit).
59
\end{itemize}
60
OpenGL is typically installed on any platform. GLUT is typically available for any platform. A replacement such as FreeGLUT should be an option. Lapack is typically also available everywhere. The only problem is GLConsole. As of writing only MeshEdit depends on GLConsole which can be obtained from SourceForge:
61
 
62
\href{http://sourceforge.net/projects/glconsole/}{http://sourceforge.net/projects/glconsole/}
63
 
64
GLConsole is extremely convenient but a well kept, little documented secret. Please do not try to install it unless you need MeshEdit. We will update this text if we either include it in GEL (assuming the authors allow it) or change to some other library.
65
 
66
GEL is compiled using CMake. Currently CMake is the official way of producing Visual Studio Solution files. A separate Mac OSX XCode project file is also maintained but may be removed if CMake makes it redundant. A set of Makefiles for various unices are also found in the package but no longer maintained.
67
 
68
To install GEL, put the libraries and header files in the logical place. CMake should instruct your build environment on how to do this. Then you can use GEL simply by including the header files and linking against the GEL libraries. On windows using Visual Studio it is important that you define SECURE\_SCL to be 0. This is because Visual Studio is a bit pedantic about not allowing access to the raw memory of STL containers.
69
 
70
\section{CGLA}
71
 
72
CGLA is a set of numerical C++ vector and matrix classes and class
73
templates designed with computer graphics in mind. CGLA stands for
465 jab 74
``Computer Graphics Linear Algebra''. Let us get right down to the obvious question: Why create another
405 jab 75
linear algebra package?
76
Well, CGLA evolved from a few matrix and vector classes because I
77
didn't have anything better. Also, I created CGLA to experiment with
78
some template programming techniques. This led to the most important
79
feature of CGLA, namely the fact that all vector types are derived
80
from the same template. 
81
 
82
This makes it easy to ensure identical semantics: Since all vectors
83
have inherited, say, the * operator from a common ancestor, it works
84
the same for all of them.
85
 
86
It is important to note that CGLA was designed for Computer Graphics 
87
(not numerical computations) and this had a number of
88
implications. Since, in computer graphics we mainly need small vectors
89
of dimension 2,3, or 4 CGLA was designed for vectors of low
90
dimensionality. Moreover, the amount of memory allocated for a vector
91
is decided by its type at compile time. CGLA does not use dynamic
92
memory. CGLA also does not use virtual functions, and most functions
93
are inline. These features all help making CGLA relatively fast. 
94
 
95
Of course, other libraries of vector templates for computer graphics
96
exist, but to my knowledge none where the fundamental templates are
97
parametrized w.r.t. dimension as well as type. In other words, we have
98
a template (\texttt{ArithVec})that gets both type
99
(e.g. \texttt{float}) and dimension 
100
(e.g. 3) as arguments. the intended use of this template is as
101
ancestor of concrete types such as \texttt{Vec3f} - a 3D floating
102
point type. 
103
 
104
The use of just one template as basis is very important, I believe,
105
since it makes it extremely simple to add new types of
106
vectors. Another very generic template is \texttt{ArithMat} which is a
107
template for matrix classes. (and not necessarily NxN matrices). 
108
 
109
From a users perspective CGLA contains a number of vector and matrix
110
classes, a quaternion and some utility classes. In summary, the most
111
important features are
112
\begin{itemize}
113
\item A number of 2, 3 and 4 d vector classes.
114
\item A number of Matrix classes.
115
\item A Quaternion class.
116
\item Some test programs.
117
\item Works well with OpenGL.
118
\end{itemize}
119
 
120
 
121
\subsection{Naming conventions}
122
 
123
Vectors in CGLA are named \texttt{VecDT} where \texttt{D} stands for
124
dimension at  \texttt{T} 
125
for type. For instance a 3D floating point vector is named
126
\texttt{Vec3f}. Other types are d (double), s (short), i (int), uc
127
(unsigned char), and usi (unsigned shourt int).
128
 
129
Matrices are similiarly named \texttt{MatDxDT}. For instance a 4D
130
\texttt{double} matrix is called \texttt{Mat4x4d}.
131
 
132
\subsection{How to use CGLA}
133
 
134
If you need a given CGLA class you can find the header file that
135
contains it in this document. Simply include the header file and use
136
the class. Remember also that all CGLA functions and classes live in
137
the CGLA namespace! Lastly, look at the example programs that came
138
with the code.
139
 
140
An important point is that you should never use the Arith... classes
141
directly. Classes whose names begin with Arith are templates used for
142
deriving concrete types. It is simpler, cleaner, and the intended
143
thing to do to only use the derived types.
144
 
145
In some cases, you may find that you need a vector or matrix class
146
that I haven't defined. If so, it is fortunate that CGLA is easy to
147
extend. Just look at, say, \texttt{Vec4f} if you need a \texttt{Vec5d}
148
class.
149
 
150
For some more specific help look at the next section where some of the
151
commonly used operations are shown.
152
 
153
 
154
\subsection{CGLA: Examples of Usage}
155
 
156
The examples below are by no means complete. Many things are
157
possible but not covered below. However, most of the common usage is
158
shown, so this should be enough to get you started. Note that in the
159
following we assume that you are \texttt{using namespace CGLA} and
160
hence don't prefix with \texttt{CGLA::}.
161
 
162
In short, to compile the examples below you would need the following
163
in the top of your file
164
 
165
\begin{verbatim}
166
#include <iostream> // For input output
167
#include "CGLA/Vec3f.h"
168
#include "CGLA/Quaternion.h"
169
#include "CGLA/Mat4x4f.h"
170
 
171
using namespace std; // For input output
172
using namespace CGLA;
173
\end{verbatim}
174
 
175
To begin with let us create 3 3D vectors. This is done as follows:
176
\begin{verbatim}
177
  Vec3f p0(10,10,10);
178
  Vec3f p1(20,10,10);
179
  Vec3f p2(10,20,10);
180
\end{verbatim}
181
if we had left out the arguments the three vectors would have been
182
uninitialized, at least in release mode. If we compile in debug mode
183
they would have been initialized to ``Not a Number'' to aid
465 jab 184
debugging. However, the $[10 \;10 \;10]^T$ vector could also have been
405 jab 185
created differently, using
186
\begin{verbatim}
187
  Vec3f p0(10);
188
\end{verbatim}
189
 
190
A very common operation is to compute the normal of a triangle from
191
the position of its vertices. Assuming the three vectors represent the
192
vertices of a triangle, we can compute the normal by finding the
193
vector from the first vertex to the two other vertices, taking the
194
cross product and normalizing the result. This is a one-liner:
195
 
196
\begin{verbatim}
197
  Vec3f n = normalize(cross(p1-p0,p2-p0));
198
\end{verbatim}
199
 
200
Quite a lot goes on in the snippet above. Observe that the - operator
201
also works for vectors. In fact almost all the arithmetic operators
202
work for vectors. You can also use assignment operators (i.e
203
\texttt{+=}) which is often faster. Then there is the function
204
\texttt{cross} which simply computes the cross product of its
205
arguments. Another frequently used function is \texttt{dot} which
206
takes the dot product. Finally the vector is normalized using the
207
function normalize.
208
 
209
Of course, we can print all or at least most CGLA entities. For
210
example 
211
 
212
\begin{verbatim}
213
  cout << n << endl;
214
\end{verbatim}
215
 
216
will print the normal vector just computed. We can also treat a vector
217
as an array as shown below
218
 
219
\begin{verbatim}
220
  float x = n[0];
221
\end{verbatim}
222
 
223
here, of course, we just extracted the first coordinate of the
224
vector. 
225
 
226
CGLA contains a number of features that are not used very frequently,
227
but which are used frequently enough to warrent inclusion. A good
228
example is assigning to a vector using spherical coordinates:
229
 
230
\begin{verbatim}
231
  Vec3f p;
232
  p.set_spherical(0.955317, 3.1415926f/4.0f , 1);
233
\end{verbatim}
234
 
235
CGLA also includes a quaternion class. Here it is used to construct a
236
quaternion which will rotate the x axis into the y axis.
237
 
238
\begin{verbatim}
239
  Quaternion q;
240
  q.make_rot(Vec3f(1,0,0), Vec3f(0,1,0));
241
\end{verbatim}
242
 
243
Next, we construct a $4\times 4$ matrix \texttt{m} and assign a translation
244
matrix to the newly constructed matrix. After that we ask the
245
quaternion to return a $4\times 4$ matrix corresponding to its
246
rotation. This rotation matrix is then multiplied onto \texttt{m}.
247
 
248
\begin{verbatim} 
249
  Mat4x4f m = translation_Mat4x4f(Vec3f(1,2,3));
250
  m *= q.get_mat4x4f();
251
\end{verbatim}
252
 
253
Just like for vectors, the subscript operator works on
254
matrices. However, in this case there are two indices. Just using one
255
index will return the ith row as a vector as shown on the first line
256
below. On the second line we see that using two indices will get us an
257
element in the matrix. 
258
 
259
\begin{verbatim} 
260
  Vec4f v4 = m[0];
261
  float c = m[0][3];
262
\end{verbatim}
263
 
264
There is a number of constructors for vectors. The default constructor
265
will create a null vector as we have already seen. We can also specify
266
all the coordinates. Finally, we can pass just a single number
267
$a$. This will create the $[a\:a\:a]^T$ vector. For instance, below we
268
create the $[1\: 1\: 1]^T$ vector. Subsequently, this vector is
269
multiplied onto \texttt{m}. 
270
 
271
\begin{verbatim}
272
  Vec3f p(1);
273
  Vec3f p2 = m.mul_3D_point(p);
274
\end{verbatim}
275
 
276
Note though that \texttt{m} is a $4\times
277
4$ matrix so ... how is that possible? Well, we use the function
278
\texttt{mul\_3D\_point} which, effectively, adds a $w=1$ coordinate to p
279
making it a 4D vector. This $w$ coordinate is stripped afterwards. In
280
practice, this means that the translation part of the matrix is also
281
applied. There is a similar function \texttt{mul\_3D\_vector} if we want
282
to transform vectors without having the translation. This function,
283
effectively, sets $w=0$.
284
 
285
Finally, CGLA is often used together with OpenGL although there is no
286
explicit tie to the OpenGL library. However, we can call the
287
\texttt{get} function of most CGLA entities to get a pointer to the
288
contents. E.g. \texttt{p.get()} will return a pointer to the first
289
float in the 3D vector \texttt{p}. This can be used with OpenGL's
290
``v'' functions as shown below. 
291
 
292
\begin{verbatim}
293
  glVertex3fv(p.get());
294
\end{verbatim}
295
 
296
\section{HMesh}
297
 
298
HMesh's \texttt{Manifold} class is a halfedge based mesh data structure. This type of data structure restricts you to manifolds, but allows for general N-gons for arbitrary N and makes manipulation somewhat simpler than other types of data structures for meshes.  Note that this particular implementation does not allow more than one loop per face. 
299
 
300
The central element of the data structure is the halfedge. It is not the goal here to discuss the halfedge data structure, but it is probably helpful that you see the syntax for e.g. going from a halfedge to the next. It may also be useful to understand the basics of the data structure. The goal is to explain the basics and show how to navigate around the mesh.
301
 
302
First of all, we only access elements via iterators. An iterator can be thought of as a pointer, but it is in fact an iterator to an STL list. Given a 
303
\begin{verbatim}
304
Manifold m;
305
\end{verbatim}
306
we can loop over all its faces using
307
\begin{verbatim}
308
for(FaceIter f = m.faces_begin(); f != m.faces_end(); ++f)
309
  f->touched = -1;
310
\end{verbatim}
311
 
312
The loop above visited each face f of m and assigned -1 to a variable called touched. The touched variables are simply integers, and a touched variable is associated with each face.
313
 
314
Note that we only refer to faces via their iterators. The same is true of vertices and halfedges. Vertices and halfedges also have touched variables, so the above code could be repeated exactly for vertices and halfedges. Here goes:
315
\begin{verbatim}
316
for(VertexIter v = m.vertices_begin(); 
317
     v != m.vertices_end(); ++v)
318
  v->touched = -1;
319
 
320
for(HalfEdgeIter h = m.halfedges_begin(); 
321
     h != m.halfedges_end(); ++h)
322
  h->touched = -1;
323
\end{verbatim}
324
 
325
What would be the use of setting the touched values to -1? Well, if we want to keep track of whether we have visited a given halfedge, we can initially set the touched values of all halfedges to -1 and then as we visit a halfedge set its value to something different. The same idea can be applied to faces and vertices.
326
 
327
Another common use of the touched values is as indices into an array. In this way, you can store auxiliary data with the mesh. In fact, there is a function\\ \texttt{enumerate\_vertices()} which enumerates all vertices (so each has a unique positive (or zero) integer. Similar functions exist for faces and halfedges.
328
 
329
So far, we have only discussed touched, but halfedges, faces, and vertices contain other data members. Most of these are iterators pointing to other edges, faces, or vertices. In the following, we briefly mention each member of the datastructures and show how to get to it. To get the opposite halfedge, use
330
\begin{verbatim}
331
HalfEdgeIter h2 = h1->opp;
332
\end{verbatim}
333
 
334
To get the next halfedge (in counterclockwise cycle around the vertex), 
335
\begin{verbatim}
336
HalfEdgeIter h2 = h1->next;
337
\end{verbatim}
338
the previous halfedge, 
339
\begin{verbatim}
340
HalfEdgeIter h2 = h1->prev;
341
\end{verbatim}
342
the vertex they point at, 
343
\begin{verbatim}
344
VertexIter v = h->vert
345
\end{verbatim}
346
and the face whose counter clockwise halfedge loop they are one segment of.
347
\begin{verbatim}
348
FaceIter f = h->face;
349
\end{verbatim}
350
A vertex knows its position which is a Vec3f. To get the position, use
351
\begin{verbatim}
352
Vec3f pos = v->pos;
353
\end{verbatim}
354
A vertex also knows one outgoing halfedge:
355
\begin{verbatim}
356
HalfEdgeIter h = v->he;
357
\end{verbatim}
358
A face knows one halfedge in its counterclockwise cycle
359
\begin{verbatim}
360
HalfEdgeIter h = face->last;
361
\end{verbatim}
362
 
363
 
364
To move around a face, we can repeatedly ask for the next halfedge. Likewise, to move around a vertex, we can repeatedly ask for \texttt{h->opp->next} but it is not elegant, and in both cases we need a stop condition. Circulators encapsulate the task of visiting all edges delimiting a face or radiating from a vertex in a nice(r) way.
365
 
366
Given 
367
\begin{verbatim}
368
VertexIter v
369
\end{verbatim}
370
we can go:
371
\begin{verbatim}
372
VertexCirculator vc(v);
373
Vec3f p(0.0);
374
for( ; !vc.end(); ++vc)
375
  p += vc.get_vertex()->pos;
376
p/= vc.no_steps();
377
\end{verbatim}
378
 
379
So what does the code above do? It visits each neighboring vertex to a given vertex and computes the average position which is useful in a number of situations.
380
 
381
We can do the same thing with a FaceCirculator
382
\begin{verbatim}
383
FaceCirculator fc(f);
384
Vec3f p(0.0);
385
for( ; !fc.end(); ++fc)
386
  p += fc.get_vertex()->pos;
387
p/= fc.no_steps();
388
\end{verbatim}
389
 
390
Internally, both the vertex and face circulators simply keep track of a single halfedge. So \texttt{get\_vertex()} returns the vertex pointed to by that halfedge. There are similar functions for getting the halfedge itself, the face it points to or its opposite halfedge.
391
 
465 jab 392
Use circulators whenever you need to go around a face or vertex, but not when you need more complex navigation which is not easily expressed as circulation. In that case it is probably easier to just use the next, prev,  and opp to move around. To work with the functions and classes just described, you will need the following header files:
393
\begin{verbatim}
394
#include <HMesh/Manifold.h>
395
#include <HMesh/FaceCirculator.h>
396
#include <HMesh/VertexCirculator.h>
397
\end{verbatim}
398
\subsection{Extended Example: Edge Flipping}
399
The following example demonstrates how to create a Manifold and add polygons (in this case triangles) and finally how to flip edges of a manifold. First, let us define some vertices
400
\begin{verbatim}
401
	vector<Vec3f> vertices(3);
402
	vertices[0] = p1;
403
	vertices[1] = p2;
404
	vertices[2] = p3;
405
\end{verbatim}
406
and then a vector of faces. The faces vector contains the number of vertices of each face in the mesh we want tor produce. Initially, we want to make a mesh with just one single triangle, so we produce a vector of one element and set that element to 3.
407
\begin{verbatim}
408
vector<int> faces(1);
409
faces[0] = 3;
410
\end{verbatim}
405 jab 411
 
465 jab 412
Next, we need the index vector. This vector is a long list of indices to vertices. For each face, we store the indices to its vertices in this vector. Since we just have a triangle, we store three vertices in the vector.
413
\begin{verbatim}
414
vector<int> indices(3);
415
indices[0]=0;
416
indices[1]=1;
417
indices[2]=2;
418
\end{verbatim}
419
So, if we had had two triangles, we would have stored six indices. Finally, we call \texttt{build\_manifold} with the indexed face set data, we have compiled. 
420
\begin{verbatim}
421
build_manifold(mani,        	// The triangle mesh.
422
			 3,           		// Number of vertices.
423
			 &vertices[0],	// Pointer to vertices.
424
			 1,           		// Number of faces.
425
			 &faces[0],  	// Pointer to faces.
426
			 &indices[0]);	// Pointer to indices.
427
\end{verbatim}
428
The result of the above is a mesh with a single triangle. Note that \texttt{build\_manifold} is in
429
\begin{verbatim}
430
#include <HMesh/build_manifold.h>
431
\end{verbatim}
432
Next, we try to split an edge. We grab the first halfedge and split it: Now we have two triangles. Splitting the halfedge changes the containing triangle to a quadrilateral, and we call triangulate to divide it into two triangls again.
433
\begin{verbatim}
434
	HalfEdgeIter h = m.halfedges_begin();
435
	m.split_edge(h);
436
	m.triangulate(h->face); 
437
\end{verbatim}
438
Now, we obtain an iterator to the first of our (two) triangles.
439
\begin{verbatim}
440
	FaceIter f1 =m.faces_begin();
441
\end{verbatim}
442
We create a vertex, \texttt{v}, at centre of \texttt{f1} and insert it in that face:
443
\begin{verbatim}
444
	VertexIter v=m.create_vertex(centre(f1));
445
	m.face_insert_point(f1,v);
446
\end{verbatim}
447
The code below, marks one of a pair of halfedges (with the value 1).
448
\begin{verbatim}
449
for(HalfEdgeIter h = m.halfedges_begin(); h!=m.halfedges_end(); ++h)
450
  h->touched =0;
451
for(HalfEdgeIter h = m.halfedges_begin(); h!=m.halfedges_end(); ++h)
452
  if(h->opp->touched == 0)
453
    h->touched = 1;
454
\end{verbatim}
455
Next, we set the \texttt{flipper} variable to point to the first of these halfedges:
456
\begin{verbatim}
457
	flipper = m.halfedges_begin();
458
\end{verbatim}
459
The long piece of code below examines a halfedge pointed to by \texttt{flipper}, and if it is not a boundary halfedge it tries to flip it. Then it increments the \texttt{flipper} variable until it lands on a halfedge marked with 1. If \texttt{flipper} reaches the end of the list of halfedges, we start over.
460
\begin{verbatim}
461
if(is_boundary(flipper))
462
  cout << "boundary edge: could not flip" << endl;
463
else
464
  if(!m.flip(flipper))
465
    cout << "could not flip" << endl;
466
do
467
{
468
  ++flipper;
469
  if(flipper==m.halfedges_end())
470
    flipper = m.halfedges_begin();
471
}
472
while(flipper->touched == 0); 
473
\end{verbatim}
474
\subsection{Extended Example: Implicit Surface Polygonization}
475
This example demonstrates how we can produce an HMesh Manifold from an implicitly defined object. Say, we have some function, \texttt{f}, which represents an object by what we sometimes call a characteristic function (\texttt{f<0} inside the object, \texttt{f>0} outside, \texttt{f=0}  on the surface). We wish to extract the 0 set as a Manifold.
476
 
477
First, we create a 3D grid of floating point values as shown below.
478
\begin{verbatim}
479
RGrid<float> grid(Vec3i(N,N,N)); // Voxel grid on which to sample 
480
\end{verbatim}
481
We sample \texttt{f} on that grid
482
\begin{verbatim}
483
for(int i=0;i<N;++i)
484
  for(int j=0;j<N;++j)
485
    for(int k=0;k<N;++k)
486
      grid[Vec3i(i,j,k)] = f(i,j,k);
487
\end{verbatim}
488
Next, we call a function which polygonizes the 0 set. The first two arguments are the grid and the manifold, the third argument is set to 0, since we want the zero level set, and the final argument tells the polygonizer to push the vertices onto the zero set (as opposed to leaving them in their initial position which is close to but not actually on the zero set).
489
\begin{verbatim}
490
cuberille_polygonize(grid, mani, 0.0, true);
491
\end{verbatim}
492
In post-processing, we triangulate \texttt{mani} and remove ugly triangles. Caps are triangles with one big angle (close to 180 degrees). Needles are triangles with one very short edge.
493
\begin{verbatim}
494
triangulate(mani);
495
remove_caps_from_trimesh(mani, M_PI * 0.85);
496
remove_needles_from_trimesh(mani, 1e-2);
497
\end{verbatim}
498
Note that you need the following header files:
499
\begin{verbatim}
500
#include <Geometry/RGrid.h>
501
#include <HMesh/volume_polygonize.h>
502
\end{verbatim}
503
to use the grid and the polygonizer, respectively.
504
 
505
 
506
\section{GLGraphics}
507
 
508
GLGraphics contains a set of tools needed for visualization of geometrical objects. In particular, this namespace contains mesh rendering facilities, and in the following, a very simple example is given. This section is basically a walk through of the simplest GEL program that draws a mesh. Apart from GEL, we also use OpenGL and that requires a toolkit for connecting with the window system. With GEL programs, one generally uses GLUT, and this example is not an exception.
509
 
510
For starters, we need to include some files mostly from GLGraphics but also the header file for the mesh load function of HMesh and the GLEW header. Why incude glew.h rather than just gl.h? Well, HMesh/draw.h contains some draw functions which rely on shaders, and since gl.h is normally not up to date, it makes sense to use glew.h instead of gl.h even though it inflicts a dependency. We include GLGraphics/gel\_glut.h rather than the normal glut.h. That is because the GEL version works the same on Windows, OSX, Linux and probably other platforms.
511
 
512
We also open the namespaces we will need and declare two globals: The view controller \texttt{view\_ctrl} and the mesh \texttt{mani}. The view controller class is rather complex. It more or less insulates you from having to deal directly with the projections and transformations of OpenGL: It sets up a perspective projection and also the modelview transformation needed to view the object. From a user interface perspective, it defines a virtual trackball which allows you to rotate the object.
513
\begin{verbatim}
514
#include <GL/glew.h>
515
#include <GLGraphics/gel_glut.h>
516
#include <GLGraphics/draw.h>
517
#include <GLGraphics/GLViewController.h>
518
#include <HMesh/load.h>
519
 
520
using namespace std;
521
using namespace HMesh;
522
using namespace CGLA;
523
using namespace GLGraphics;
524
 
525
GLViewController* view_ctrl;
526
Manifold mani;
527
\end{verbatim}
528
 
529
Below is the display function which is called from the GLUT event loop whenever an event that requires redrawing is received. All drawing takes place inside this function. In this case, it is simple since all we need to do is clear the screen (and depth buffer), set up a new modelview transformation with the view controller and then call the draw function followed by a swap of buffers.
530
 
531
The draw function sends the polygons to the graphics card. It is very inefficient and uses the fixed function pipeline. At a minimum one should draw the polygons to a display list, but for a simple example, this works. The buffer swap is because we use double buffering; in other words, we draw to the back buffer.
532
\begin{verbatim}
533
void display() {
534
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
535
    view_ctrl->set_gl_modelview();
536
    draw(mani);
537
    glutSwapBuffers();
538
}
539
 
540
\end{verbatim}
541
The mouse callback below is called when a mouse event has occurred. This is where user input  for the virtual trackball is registered. We check whether a mouse is pressed or released. If it is pressed, we detect which button and then ``grab'' the virtual trackball in order to rotate, zoom, or pan depending on which button is pressed. If the button is released, we inform the view controller of this.
542
\begin{verbatim}
543
void mouse(int button, int state, int x, int y)  {
544
    Vec2i pos(x,y);
545
    if (state==GLUT_DOWN) 
546
    {
547
        if (button==GLUT_LEFT_BUTTON) 
548
            view_ctrl->grab_ball(ROTATE_ACTION,pos);
549
        else if (button==GLUT_MIDDLE_BUTTON) 
550
            view_ctrl->grab_ball(ZOOM_ACTION,pos);
551
        else if (button==GLUT_RIGHT_BUTTON) 
552
            view_ctrl->grab_ball(PAN_ACTION,pos);
553
    }
554
    else if (state==GLUT_UP)
555
        view_ctrl->release_ball();
556
}
557
\end{verbatim}
558
If the mouse moves with a button depressed, we register motion in the function below. The new mouse position is sent to the view controller
559
which then, e.g., rotates the trackball accordingly. Note that this function ends by posting a redisplay event, i.e. informs GLUT that display should be called at the earliest convenience.
560
\begin{verbatim}
561
void motion(int x, int y) {
562
    Vec2i pos(x,y);
563
    view_ctrl->roll_ball(Vec2i(x,y));
564
    glutPostRedisplay();
565
}
566
\end{verbatim}
567
 
568
Finally, in the main function, we first load the mesh and then setup glut. This mostly involves defining the callback function just described. We set up the view controller passing it window dimensions and the size and position of the bounding sphere of the object just loaded. Finally, we do some minimal OpenGL set up: Essentially clearing the depth buffer and enabling lights. Finally, we pass control to the GLUT event loop.
569
\begin{verbatim}
570
int main(int argc, char** argv) {
571
    string file = "bunny-little.x3d";
572
    load(file, mani);
573
 
574
    glutInitDisplayMode(GLUT_RGBA|GLUT_DOUBLE|GLUT_DEPTH);
575
    glutInitWindowSize(800, 800);
576
    glutInit(&argc, argv);
577
    glutCreateWindow("GEL Example Program");
578
    glutDisplayFunc(display);
579
    glutMouseFunc(mouse);
580
    glutMotionFunc(motion);
581
 
582
    Vec3f bsphere_center(0,0,0);   
583
    float bsphere_radius = 3;
584
    mani.get_bsphere(bsphere_center,bsphere_radius);
585
    view_ctrl = new GLViewController(800,800, 
586
        bsphere_center,bsphere_radius*2);
587
 
588
    glClearColor(0.0f, 0.0f, .5f, 0.f);
589
    glEnable(GL_DEPTH_TEST);
590
    glEnable(GL_LIGHTING);
591
    glEnable(GL_LIGHT0);
592
 
593
    glutMainLoop();
594
}
595
\end{verbatim}
596
Note that the mesh can be either an X3D, OBJ, PLY, or OFF files. However, only the geometry is loaded not the textures. In the case of OBJ files, we also convert all polygons to triangles, because the TriMesh loader (see the geometry class) is in fact used.
597
 
598
Apart from the drawing functions and the viewcontroller, GLGraphics also contains  other useful pieces: A few functions which greatly simplify shader loading and SOIL. SOIL is a small library for loading images. It is designed for OpenGL and makes it easy to e.g. save screen shots or load an image for use as a texture. While the need for image loading is fairly obvious, it might be more difficult to see the need for shader loading functions. However, it \textit{is} a little tricky to make shader programs in C++ using GLSL but the problems are almost all silly little things. For instance, how do you robustly read from a text file? How do you compile a shader and check for errors? In OpenGL What is the difference between a GLSL program and a shader anyway?
599
 
600
The short answer to the last question is this: A "shader" can be either a vertex shader, a geometry shader (in OpenGL 2.0) or a fragment shader. A "program" is a linked combination of these three types of shaders. Note that you don't have to have a geometry shader.
601
 
602
These functions attempt to obviate the need for an answer to the first two questions by providing an API for loading and compiling shaders and checking for errors. However, I don't include functions for creating the program, attaching shaders and linking. Why not? Because the need for flexibility means that the API would be just as complex as just using the OpenGL functions directly! However, loading shaders and checking for errors in compiled shaders is different. It makes sense to wrap that. It also makes sense to wrap the error checking for programs that are linked, so there is a function for that too.
603
 
604
Since shader loading and error check sandwhich the two calls needed for compilation, the most important function in this API, \texttt{create\_glsl\_shader}, loads a shader, compiles it, checks for errors and returns the shader handle. There is also a version which creates a shader from a string.
605
 
606
There is some code below to illustrate usage.	
607
\begin{verbatim}
608
GLuint vs = create_glsl_shader(GL_VERTEX_SHADER, 
609
  shader_path, "tri.vert");
610
GLuint gs = create_glsl_shader(GL_GEOMETRY_SHADER_EXT, 
611
  shader_path, "tri.geom");
612
GLuint fs = create_glsl_shader(GL_FRAGMENT_SHADER, 
613
  shader_path, "tri.frag");
614
 
615
prog_P0 = glCreateProgram();
616
 
617
if(vs) glAttachShader(prog_P0, vs);
618
if(gs) glAttachShader(prog_P0, gs);
619
if(fs) glAttachShader(prog_P0, fs);
620
 
621
// Specify input and output for the geometry shader.
622
// Note that this must be done before linking the program.
623
glProgramParameteriEXT(prog_P0,GL_GEOMETRY_INPUT_TYPE_EXT,
624
  GL_TRIANGLES);
625
glProgramParameteriEXT(prog_P0,GL_GEOMETRY_VERTICES_OUT_EXT,3);
626
glProgramParameteriEXT(prog_P0,GL_GEOMETRY_OUTPUT_TYPE_EXT,
627
  GL_TRIANGLE_STRIP);
628
 
629
// Link the program object and print out the info log
630
glLinkProgram(prog_P0);
631
print_glsl_program_log(prog_P0);
632
 
633
// Install program object as part of current state
634
glUseProgram(prog_P0);
635
 
636
// Set the value of a uniform
637
glUniform2f(glGetUniformLocation(prog_P0,"WIN_SCALE"), 
638
  win_size_x/2.0, win_size_y/2.0);
639
\end{verbatim}
640
 
641
\section{LinAlg}
642
Sometimes the simple 2,3,4 dimensional vectors from CGLA just don't cut it. We often need to solve large linear systems. The LinAlg namespace contains some vector and matrix classes and an interface to Lapack. This provides fairly easy way to do many numerical computations.
643
 
644
In the example below, we  find coefficients for $a x^2 + b y^2 +cxy + dx + ey + f$ such that the surface contains six specific points. The coefficients are found by solving the linear system $\mathbf A \mathbf x = \mathbf b$ where the rows of the $\mathbf A$ matrix contains $x^2 \; y^2 \; xy \; x \; y \; 1$ computed from the $xy$ positions of each point, $\mathbf x$ contains the unknown coefficients, and $\mathbf b$ on the right hand side contains the $z$ values of the six points. Using LinAlg, we implement it as follows: 
645
\begin{verbatim}
646
    CMatrix A(6,6);
647
    CVector b(6);
648
 
649
    for(int i=0;i<6;++i)
650
    {
651
        A.set(i,0,pts[i][0]*pts[i][0]);
652
        A.set(i,1,pts[i][1]*pts[i][1]);
653
        A.set(i,2,pts[i][0]*pts[i][1]);
654
        A.set(i,3,pts[i][0]);
655
        A.set(i,4,pts[i][1]);
656
        A.set(i,5,1);
657
        b[i] = pts[i][2];
658
    }
659
    CVector x;
660
    LinearSolve(A,b,x);
661
\end{verbatim}
662
While this example shows only the function for solving a system which should have a solution, there are also functions for solving overdetermined systems in the least square sense and finding the minimum norm solution to underdetermined systems using singular value decomposition. There is also a function for finding eigenvalues and eigenvectors of symmetric matrices and more.
663
\section{Geometry}
664
The Geometry namespace contains many different classes and functions. Roughly, we can divide the contents into
665
\begin{itemize}
666
\item Spatial data structures such as kd tress, bounding box hierarchies, and BSP trees.
667
\item A Triangle mesh data structure called TriMesh. TriMesh is different from HMesh in that it only contains triangles. It is much more an API for real time rendering which also facilitates material properties and not just geometry.
668
\item Classes and functions for dealing with volume data. HGrid and RGrid are, respectively, a class for hierarchically stored voxel grids (sparse grids) and regular grids. There is a number of functions for traversing voxel grids systematically and along rays.
669
\item Surprisingly, you will also find my C++ port of Jules Bloomenthals isosurface polygonizer and an implementation of Luiz Velho's parametric surface tiler in this namespace.
670
\end{itemize}
671
One of the most frequently used tools in Geometry is the kD tree class which merits an example. In the code below, we create a kD tree which had \texttt{Vec3f}s as keys and integers as data. Note that since the kD tree is a template, we could use any other type as both key and data, however, only CGLA vectors have been tested as key types, but there are no restrictions on data types.
672
\begin{verbatim}
673
    KDTree<Vec3f,int> tree;
674
    for(int i=0;i<N;++i)
675
        {
676
            Vec3f p0;
677
            make_ran_point(p0);
678
            tree.insert(p0, i);
679
        }
680
    tree.build();
681
\end{verbatim}
682
Note also, that before we use the tree it needs to be built. This is because the internal representation is a balanced binary heap which it is easier to build at the end rather than maintain during insertion of new points. To look for the points near a given point, we can call
683
\begin{verbatim}
684
    Vec3f p0;
685
    // ...
686
    std::vector<Vec3f> keys;
687
    std::vector<int> vals;
688
    float radius = 1.95f;
689
    int N = tree.in_sphere(p0, radius , keys, vals);
690
\end{verbatim}
691
This will find all points in a radius of 1.95 units from \texttt{p0}. The keys and values are stored in \texttt{keys} and \texttt{vals} when the functions returns. The return value of the function is the number of points found within the radius. 
692
 
693
Sometimes, we simply want the point closest to a given point \texttt{p0} which is done using \texttt{tree.closest\_point(p0, d, pos, x)} where the two first arguments indicate the point and the maximum search radius. The two last arguments are set to the key and data of the point found to be closest to \texttt{p0}. The function returns true if a closest point was found.
694
 
695
 
696
 
697
 
698
\section{Util}
699
The Util namespace contains a rather mixed bag of utilities many of which are quite useful and very diverse. 
700
\begin{trivlist}
701
\item The XML parser for instance is written from scratch and the basis for our X3D loader although it will load any XML file. 
702
\item \texttt{Grid2D} is often useful when we want to deal with e.g. simulation on a 2D grid. 
703
\item The \texttt{Timer} class is nearly indispensable when we need to time something, and it is very easy to use. 
704
\item \texttt{ArgExtracter} is a class for getting the command line arguments to a program. 
705
\item There is more - for instance a resource manager template, some string utility templates and a hash table implementation.
706
\end{trivlist}
707
Arguably, it would be more pretty if these pieces had been more thematically linked, but many are used a lot, and it seemed reasonable to have somewhere to stick various functions and classes that did not naturally belong elsewhere.
708
 
709
 
710
\section{Authors and License}  
405 jab 711
Many other people contribute, but the core of GEL was written (mostly) by Andreas B{\ae}entzen (jab{@}imm.dtu.dk), and any bug fixes, contributions, or questions should be addressed to me, unless you know precisely who could take care of the problem.
712
 
713
I was considering putting GEL under the LGPL. But it is a long complex
714
text. The longer any kind of document, the more chances for loopholes
465 jab 715
in my opinion. Instead, I list a few rules below. The most
405 jab 716
important one is that if you want to use GEL for some purpose, and it
717
is not crystal clear whether it is against the rules, contact me. As
718
for the rules:
719
 
720
You are allowed to use GEL for academic or commercial purposes. In
721
particular, you are welcome to give GEL to students as a basis for 
722
academic work or to use it for your own applications.
723
 
724
The biggest limitation that I want to impose is that you cannot repackage
725
GEL in any way. You are not allowed to distribute another library which
726
contains GEL or parts of GEL unless you make it clear that this other
727
library contains GEL. You are also not allowed to redistribute GEL in a
728
changed form. If you want changes to be made, contact me.
729
 
730
Of course, neither I nor my employer will give you any money or be
731
held responsible in any way, under any circumstances what so ever - no
732
matter what sort of damage GEL might inflict upon you. Not that I can
733
foresee GEL causing you any damage :-) 
734
 
735
If anything is unclear, please contact me. In fact, if you want to use
736
GEL for a bigger project, I'd appreciate an email to jab@imm.dtu.dk
737
 
738
In a project such as this, it is almost impossible to completely avoid
739
building upon fragments of other peoples source code. GEL includes an
740
obj loader based on work by Nate Robins. Some pieces of source code
465 jab 741
from Graphics Gems have also been used. Moreover, I have included
405 jab 742
rply by Diego Nehab, Princeton University, and the Simple OpenGL Image Loader by
743
Jonathan Dummer.  All of this amounts to only a fraction of the GEL source code and it should not be in violation of any license. In particular, SOIL and RPLY are under the MIT license and it is acceptable to include these packages as long as the copyright notice is retained (which it is.)
744
 
745
\end{document}