Prerequisites

This is not a Python nor a NumpPy beginner guide and you should have an intermediate level in both Python and NumPy. No prior knowledge of OpenGL is necessary because I'll explain everything.

Conventions

I will use usual naming conventions. If not stated explicitly, each script should import numpy, scipy and glumpy as:

import numpy as np

We'll use up-to-date versions (at the date of writing, i.e. August, 2017) of the different packages:

How to contribute

If you want to contribute to this book, you can:

• Report issues (https://github.com/rougier/python-opengl/issues)
• Suggest improvements (https://github.com/rougier/python-opengl/pulls)
• Correct English (https://github.com/rougier/python-opengl/issues)
• Star the project (https://github.com/rougier/python-opengl)
• Suggest a more responsive design for the HTML Book
• Spread the word about this book (Reddit, Hacker News, etc.)

Publishing

If you're an editor interested in publishing this book, you can contact me if you agree to have this version and all subsequent versions open access (i.e. online at this address), you know how to deal with restructured text (Word is not an option), you provide a real added-value as well as supporting services, and more importantly, you have a truly amazing latex book template (and be warned that I'm a bit picky about typography & design: Edward Tufte is my hero). Still here?

Book

This work is licensed under a Creative Commons Attribution-Non Commercial-Share Alike 4.0 International License. You are free to:

• Share — copy and redistribute the material in any medium or format
• Adapt — remix, transform, and build upon the material

The licensor cannot revoke these freedoms as long as you follow the license terms.

Under the following terms:

• Attribution — You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.
• NonCommercial — You may not use the material for commercial purposes.
• ShareAlike — If you remix, transform, or build upon the material, you must distribute your contributions under the same license as the original.

Code

The code is licensed under the OSI-approved BSD 2-Clause License.

The graphic pipeline

If you want to understand modern OpenGL, you have to understand the graphic pipeline and shaders. Shaders are pieces of program (using a C-like language) that are build onto the GPU and executed during the rendering pipeline. Depending on the nature of the shaders (there are many types depending on the version of OpenGL you're using), they will act at different stage of the rendering pipeline. To simplify this tutorial, we'll use only vertex and fragment shaders as shown below:

A vertex shader acts on vertices and is supposed to output the vertex position (gl_Position) on the viewport (i.e. screen). A fragment shader acts at the fragment level and is supposed to output the color (gl_FragColor) of the fragment. Hence, a minimal vertex shader is:

void main()
{
    gl_Position = vec4(0.0,0.0,0.0,1.0);
}

while a minimal fragment shader would be:

void main()
{
    gl_FragColor = vec4(0.0,0.0,0.0,1.0);
}

These two shaders are not very useful because the first shader will always output the null vertex (gl_Position is a special variable) while the second will only output the black color for any fragment (gl_FragColor is also a special variable). We'll see later how to make them to do more useful things.

One question remains: when are those shaders executed exactly ? The vertex shader is executed for each vertex that is given to the rendering pipeline (we'll see what does that mean exactly later) and the fragment shader is executed on each fragment (= pixel) that is generated after the vertex stage. For example, in the simple figure above, the vertex would be called 3 times, once for each vertex (1,2 and 3) while the fragment shader would be executed 21 times, once for each fragment.

Buffers

The next question is thus where do those vertices comes from ? The idea of modern GL is that vertices are stored on the CPU and need to be uploaded to the GPU before rendering. The way to do that is to build buffers onto the CPU and to send these buffers onto the GPU. If your data does not change, no need to upload them again. That is the big difference with the previous fixed pipeline where data were uploaded at each rendering call (only display lists were built into GPU memory).

But what is the structure of a vertex ? OpenGL does not assume anything about your vertex structure and you're free to use as many information you may need for each vertex. The only condition is that all vertices from a buffer have the same structure (possibly with different content). This again is a big difference with the fixed pipeline where OpenGL was doing a lot of complex rendering stuff for you (projections, lighting, normals, etc.) with an implicit fixed vertex structure. The good news is that you're now free to do anything you want, but the bad news is that you have to program just everything.

Let's take a simple example of a vertex structure where we want each vertex to hold a position and a color. The easiest way to do that in python is to use a structured array using numpy:

data = np.zeros(4, dtype = [ ("position", np.float32, 3),
                             ("color",    np.float32, 4)] )

We just created a CPU buffer with 4 vertices, each of them having a position (3 floats for x,y,z coordinates) and a color (4 floats for red, blue, green and alpha channels). Note that we explicitly chose to have 3 coordinates for position but we may have chosen to have only 2 if were to work in two-dimensions. Same holds true for color. We could have used only 3 channels (r,g,b) if we did not want to use transparency. This would save some bytes for each vertex. Of course, for 4 vertices, this does not really matter but you have to realize it will matter if your data size grows up to one or ten million vertices.

Variables

Now, we need to explain our shaders what to do with these buffers and how to connect them together. So, let's consider again a CPU buffer of 4 vertices using 2 floats for position and 4 floats for color:

data = np.zeros(4, dtype = [ ("position", np.float32, 2),
                             ("color",    np.float32, 4)] )

We need to tell the vertex shader that it will have to handle vertices where a position is a tuple of 2 floats and color is a tuple of 4 floats. This is precisely what attributes are meant for. Let us change slightly our previous vertex shader:

attribute vec2 position;
attribute vec4 color;
void main()
{
    gl_Position = vec4(position, 0.0, 1.0);
}

This vertex shader now expects a vertex to possess 2 attributes, one named position and one named color with specified types (vec2 means tuple of 2 floats and vec4 means tuple of 4 floats). It is important to note that even if we labeled the first attribute position, this attribute is not yet bound to the actual position in the numpy array. We'll need to do it explicitly at some point in our program and there is no magic that will bind the numpy array field to the right attribute, you'll have to do it yourself, but we'll see that later.

The second type of information we can feed the vertex shader is the uniform that may be considered as constant value (across all the vertices). Let's say for example we want to scale all the vertices by a constant factor scale, we would thus write:

uniform float scale;
attribute vec2 position;
attribute vec4 color;
void main()
{
    gl_Position = vec4(position*scale, 0.0, 1.0);
}

Last type is the varying type that is used to pass information between the vertex stage and the fragment stage. So let us suppose (again) we want to pass the vertex color to the fragment shader, we now write:

uniform float scale;
attribute vec2 position;
attribute vec4 color;
varying vec4 v_color;

void main()
{
    gl_Position = vec4(position*scale, 0.0, 1.0);
    v_color = color;
}

and then in the fragment shader, we write:

varying vec4 v_color;

void main()
{
    gl_FragColor = v_color;
}

The question is what is the value of v_color inside the fragment shader ? If you look at the figure that introduced the gl pipeline, we have 3 vertices and 21 fragments. What is the color of each individual fragment ?

The answer is the interpolation of all 3 vertices color. This interpolation is made using the distance of the fragment to each individual vertex. This is a very important concept to understand. Any varying value is interpolated between the vertices that compose the elementary item (mostly, line or triangle).

Ok, enough for now, we'll see an explicit example in the next chapter.

Bindings

• Pyglet is a pure python cross-platform application framework intended for game development. It supports windowing, user interface event handling, OpenGL graphics, loading images and videos and playing sounds and music. It works on Windows, OS X and Linux.
• PyOpenGL is the most common cross platform Python binding to OpenGL and related APIs. The binding is created using the standard ctypes library, and is provided under an extremely liberal BSD-style Open-Source license.
• ModernGL is a wrapper over OpenGL that simplifies the creation of simple graphics applications like scientific simulations, small games or user interfaces. Usually, acquiring in-depth knowledge of OpenGL requires a steep learning curve. In contrast, ModernGL is easy to learn and use, moreover it is capable of rendering with the same performance and quality, with less code written.
• Ctypes bindings can also be generated quite easily thanks to the gl.xml file provided by the Khronos group that defines the OpenGL and OpenGL API Registry. The number of functions and constants given in the table above have been computed using the code/chapter-02/registry.py program that parses the gl.xml file for each API and version and count the relevant features.

Engines

• The Visualization Toolkit (VTK) is an open-source, freely available software system for 3D computer graphics, image processing, and visualization. It consists of a C++ class library and several interpreted interface layers including Tcl/Tk, Java, and Python.
• Processing is a programming language, development environment, and online community. Since 2001, Processing has promoted software literacy within the visual arts and visual literacy within technology. Today, there are tens of thousands of students, artists, designers, researchers, and hobbyists who use Processing for learning, prototyping, and production.
• NodeBox for OpenGL is a free, cross-platform library for generating 2D animations with Python programming code. It is built on Pyglet and adopts the drawing API from NodeBox for Mac OS X. It has built-in support for paths, layers, motion tweening, hardware-accelerated image effects, simple physics and interactivity.
• Panda3D is a 3D engine: a library of subroutines for 3D rendering and game development. The library is C++ with a set of Python bindings. Game development with Panda3D usually consists of writing a Python or C++ program that controls the Panda3D library.
• VPython makes it easy to create navigable 3D displays and animations, even for those with limited programming experience. Because it is based on Python, it also has much to offer for experienced programmers and researchers.

Libraries

• Glumpy is a python library for scientific visualization that is both fast, scalable and beautiful. Glumpy leverages the computational power of modern Graphics Processing Units (GPUs) through the OpenGL library to display very large datasets and offers an intuitive interface between numpy and modern OpenGL. We'll use it extensively in this book.
• Vispy is the sister project of glumpy. It is a high-performance interactive 2D/3D data visualization library and offer a high-level interface for scientific visualization. The difference between glumpy and vispy is approximately the same as the difference between numpy and scipy even though vispy is independent of glumpy and vice-versa.

Normalize Device Coordinates

Figure

Normalized Device Coordinates (NDC) where bottom-left corner coordinate is (-1,-1) and top-right corner is (+1,+1).

Before even diving into actual code, it is important to understand first how OpenGL handles coordinates. More precisely, OpenGL considers only coordinates (x,y,z) that fall into the space where -1 ≤ x,y,z ≤ +1. Any coordinates that are outside this range will be discarded or clipped (i.e. won't be visible on screen). This is called Normalized Device Coordinates, or NDC for short. This is something you cannot change because it is part of the OpenGL API and implemented in your hardware (GPU). Consequently, even if you intend to render the whole universe, you'll have utlimately to fit it into this small volume.

The second important fact to know is that x coordinates increase from left to right and y coordinates increase from bottom to top. For this latter one, it is noticeably different from the usual convention and this might induce some problems, especially when you're dealing with the mouse pointer whose y coordinate goes the other way around.

Triangulation

Figure

Different triangulation of the same quad using from 2 to 5 triangles.

Triangulation of a surface means to find a set of triangles, which covers a given surface. This can be a tedious process but fortunately, there exist many different methods and algorithms to perform such triangulation automatically for any 2D or 3D surface. The quality of the triangulation is measured in terms of the closeness to the approximated surface, the number of triangles necessary (the smaller, the better) and the homogeneity of the triangles (we prefer to have triangles that have more or less the same size and to not have any degenerated triangle).

In our case, we want to render a square and we need to find the proper triangulation (which is not unique as illustrated on the figure). Since we want to minimize the number of triangles, we'll use the 2 triangles solution that requires only 4 (shared) vertices corresponding to the four corners of the quad. However, you can see from the figure that we could have used different triangulations using more vertices, and later in this book we will just do that (but for a reason).

Considering the NDC, our quad will thus be composed of two triangles:

• One triangle described by vertices (-1,+1), (+1,+1), (-1,-1)
• One triangle described by vertices (+1,+1), (-1,-1), (+1,-1)

Here we can see that vertices (-1,-1) and (+1,+1) are common to both triangles. So instead of using 6 vertices to describe the two triangles, we can re-use the common vertices to describe the whole quad. Let's name them:

• V₀: (-1,+1)
• V₁: (+1,+1)
• V₂: (-1,-1)
• V₃: (+1,-1)

Our quad can now be using triangle (V₀,V₁,V₂) and triangle (V₁,V₂,V₃). This is exactly what we need to tell OpenGL.

GL Primitives

Figure

Common OpenGL rendering primitives.

Ok, now things are getting serious because we need to actually tell OpenGL what to do with the vertices, i.e. how to render them? What do they describe in terms of geometrical primitives? This is quite an important topic since this will determine how fragments will actually be generated as illustrated on the image below:

Mostly, OpenGL knows how to draw (ugly) points, (ugly) lines and (ugly) triangles. For lines and triangles, there exist some variations depending if you want to specify very precisely what to draw or if you can take advantage of some implicit assumptions. Let's consider lines first for example. Given a set of four vertices (V₀,V₁,V₂,V₃), you migh want to draw segments (V₀,V₁)``(V₂,V₃) using GL_LINES or a broken line (V₀,V₁,V₂,V₃) using using GL_LINE_STRIP or a closed broken line (V₀,V₁,V₂,V₃,V₀,) using GL_LINE_LOOP. For triangles, you have the choices of specifying each triangle individually using GL_TRIANGLES or you can tell OpenGL that triangles follow an implicit structure using GL_TRIANGLE_STRIP. For example, considering a set of vertices (Vᵢ), GL_TRIANGLE_STRIP will produce triangles (Vᵢ,Vᵢ₊₁,Vᵢ₊₂). There exist other primitives but we won't used them in this book because they're mainly related to geometry shaders that are not introduced.

If you remember the previous section where we explained that our quad can be described using using triangle (V₀,V₁,V₂) and triangle (V₁,V₂,V₃), you can now realize that we can take advantage or the GL_TRIANGLE_STRIP primitive because we took care of describing the two triangles following this implicit structure.

Interpolation

Figure

The Barycentric interpolation f of a fragment p is given by f = 𝛌₁f₁ + 𝛌₂f₂ + 𝛌₃f₃

The choice of the triangle as the only surface primitive is not an arbitrary choice, because a triangle offers the possibility of having a nice and intuitive interpolation of any point that is inside the triangle. If you look back at the graphic pipeline as it has been introduced in the Modern OpenGL section, you can see that the rasterisation requires for OpenGL to generate fragments inside the triangle but also to interpolate values (colors on the figure). One of the legitimate questions to be solved is then: if I have a triangle (V₁,V₂,V₃), each summit vertex having (for example) a different color, what is the color of a fragment p inside the triangle? The answer is barycentric interpolation as illustrated on the figure on the right.

More precisely, for any point p inside a triangle A = (V₁,V₂,V₃), we consider triangles:

• A₁ = (P,V₂,V₃)
• A₂ = (P,V₁,V₃)
• A₃ = (P,V₁,V₂)

And we can define (using area of triangles):

• 𝛌₁ = A₁/A
• 𝛌₂ = A₂/A
• 𝛌₃ = A₃/A

Now, if we attach a value f₁ to vertex V₁, f₂ to vertex V₂ and f₃ to vertex V₃, the interpolated value f of p is given by: f = 𝛌₁f₁ + 𝛌₂f₂ + 𝛌₃f₃ You can check by yourself that if the point p is on a border of the triangle, the resulting interpolated value f is the linear interpolation of the two vertices defining the segment the point p belongs to.

This barycentric interpolation is important to understand even if it is done automatically by OpenGL (with some variation to take projection into account). We took the example of colors, but the same interpolation scheme holds true for any value you pass from the vertex shader to the fragment shader. And this property will be used and abused in this book.

Writing shaders

Now that your window has been created, we can start writing our program, that is, we need to write a vertex and a fragment shader. For the vertex shader, the code is very simple because we took care of using the normalized device coordinates to describe our quad in the previous section. This means vertices do not need to be transformed. Nonetheless, we have to take care of sending 4D coordinates even though we'll transmit only 2D coordinates (x,y) or the final result will be undefined. For coordinate z we'll just set it to 0.0 (but any value would do) and for coordinate w, we set it to 1.0 (see section Basic Mathematics for the explanation). Note also the (commented) alternative ways of writing the shader.

attribute vec2 position;
void main()
{
  gl_Position = vec4(position, 0.0, 1.0);

  // or gl_Position.xyzw = vec4(position, 0.0, 1.0);

  // or gl_Position.xy = position;
  //    gl_Position.zw = vec2(0.0, 1.0);

  // or gl_Position.x = position.x;
  //    gl_Position.y = position.y;
  //    gl_Position.z = 0.0;
  //    gl_Position.w = 1.0;
}

For the fragment shader, it is even simpler. We set the color to red which is described by the tuple (1.0, 0.0, 0.0, 1.0) in normalized RGBA notation. 1.0 for alpha channel means fully opaque.

void main()
{
  gl_FragColor = vec4(1.0, 0.0, 0.0, 1.0);

  // or gl_FragColor.rgba = vec4(1.0, 0.0, 0.0, 1.0);

  // or gl_FragColor.rgb = vec3(1.0, 0.0, 0.0);
  //    gl_FragColor.a = 1.0;
}

Compiling the program

We wrote our shader and we need now to build a program that will link the vertex and the fragment shader together. Building such program is relatively straightforward (provided we do not check for errors). First we need to request program and shader slots from the GPU:

program  = gl.glCreateProgram()
vertex   = gl.glCreateShader(gl.GL_VERTEX_SHADER)
fragment = gl.glCreateShader(gl.GL_FRAGMENT_SHADER)

We can now ask for the compilation of our shaders into GPU objects and we log for any error from the compiler (e.g. syntax error, undefined variables, etc):

vertex_code = """
  attribute vec2 position;
  void main() { gl_Position = vec4(position, 0.0, 1.0); } """

fragment_code = """
  void main() { gl_FragColor = vec4(1.0, 0.0, 0.0, 1.0); } """

# Set shaders source
gl.glShaderSource(vertex, vertex_code)
gl.glShaderSource(fragment, fragment_code)

# Compile shaders
gl.glCompileShader(vertex)
if not gl.glGetShaderiv(vertex, gl.GL_COMPILE_STATUS):
    error = gl.glGetShaderInfoLog(vertex).decode()
    print(error)
    raise RuntimeError("Vertex shader compilation error")

gl.glCompileShader(fragment)
if not gl.glGetShaderiv(fragment, gl.GL_COMPILE_STATUS):
    error = gl.glGetShaderInfoLog(fragment).decode()
    print(error)
    raise RuntimeError("Fragment shader compilation error")

Then we link our two objects in into a program and again, we check for errors during the process.

gl.glAttachShader(program, vertex)
gl.glAttachShader(program, fragment)
gl.glLinkProgram(program)

if not gl.glGetProgramiv(program, gl.GL_LINK_STATUS):
    print(gl.glGetProgramInfoLog(program))
    raise RuntimeError('Linking error')

and we can get rid of the shaders, they won't be used again (you can think of them as .o files in C).

gl.glDetachShader(program, vertex)
gl.glDetachShader(program, fragment)

Finally, we make program the default program to be ran. We can do it now because we'll use a single program in this example:

gl.glUseProgram(program)

Uploading data to the GPU

Next, we need to build CPU data and the corresponding GPU buffer that will hold a copy of the CPU data (GPU cannot access CPU memory). In Python, things are grealty facilitated by NumPy that allows to have a precise control over number representations. This is important because GLES 2.0 floats have to be exactly 32 bits long and a regular Python float would not work (they are actually equivalent to a C double). So let us specify a NumPy array holding 4×2 32-bits float that will correspond to our 4×(x,y) vertices:

# Build data
data = np.zeros((4,2), dtype=np.float32)

We then create a placeholder on the GPU without yet specifying the size:

# Request a buffer slot from GPU
buffer = gl.glGenBuffers(1)

# Make this buffer the default one
gl.glBindBuffer(gl.GL_ARRAY_BUFFER, buffer)

We now need to bind the buffer to the program, that is, for each attribute present in the vertex shader program, we need to tell OpenGL where to find the corresponding data (i.e. GPU buffer) and this requires some computations. More precisely, we need to tell the GPU how to read the buffer in order to bind each value to the relevant attribute. To do this, GPU needs to know what is the stride between 2 consecutive elements and what is the offset to read one attribute:

                 1ˢᵗ vertex  2ⁿᵈ vertex  3ʳᵈ vertex   …
                ┌───────────┬───────────┬───────────┬┄┄
                ┌─────┬─────┬─────┬─────┬─────┬─────┬─┄
                │  X  │  Y  │  X  │  Y  │  X  │  Y  │ …
                └─────┴─────┴─────┴─────┴─────┴─────┴─┄
     offset 0 → │ (x,y)     └───────────┘
                                stride

In our simple quad scenario, this is relatively easy to write because we have a single attribute ("position"). We first require the attribute location inside the program and then we bind the buffer with the relevant offset.

stride = data.strides[0]

offset = ctypes.c_void_p(0)
loc = gl.glGetAttribLocation(program, "position")
gl.glEnableVertexAttribArray(loc)
gl.glBindBuffer(gl.GL_ARRAY_BUFFER, buffer)
gl.glVertexAttribPointer(loc, 2, gl.GL_FLOAT, False, stride, offset)

We're basically telling the program how to bind data to the relevant attribute. This is made by providing the stride of the array (how many bytes between each record) and the offset of a given attribute.

Let's now fill our CPU data and upload it to the newly created GPU buffer:

# Assign CPU data
data[...] = (-1,+1), (+1,+1), (-1,-1), (+1,-1)

# Upload CPU data to GPU buffer
gl.glBufferData(gl.GL_ARRAY_BUFFER, data.nbytes, data, gl.GL_DYNAMIC_DRAW)

Rendering

We're done, we can now rewrite the display function:

def display():
    gl.glClear(gl.GL_COLOR_BUFFER_BIT)
    gl.glDrawArrays(gl.GL_TRIANGLE_STRIP, 0, 4)
    glut.glutSwapBuffers()

Figure

A red quad rendered using Python, raw OpenGL bindings and the venerable GLUT.

The 0,4 arguments in the glDrawArrays tells OpenGL we want to display 4 vertices from our current active buffer and we start at vertex 0. You should obtain the figure on the right with the same red (boring) color. The whole source ia available from code/chapter-03/glut-quad-solid.py.

All these operations are necessary for displaying a single colored quad on screen and complexity can escalate pretty badly if you add more objects, projections, lighting, texture, etc. This is the reason why we'll stop using the raw OpenGL interface in favor of a library. We'll use the glumpy library, mostly because I wrote it, but also because it offers a tight integration with numpy. Of course, you can design your own library to ease the writing of GL Python applications.

Uniform color

Figure

A blue quad rendered using a uniform variable specifying the color of the quad.

In the previous example, we hard-coded the red color inside the fragment shader source code. But what if we want to change the color from within the Python program? We could rebuild the program with the new color but that would not be very efficient. Fortunately there is a simple solution provided by OpenGL: uniform. Uniforms, unlike attributes, do not change from one vertex to the other and this is precisely what we need in our case. We thus need to slightly modify our fragment shader to use this uniform color:

uniform vec4 color;
void main()
{
  gl_FragColor = color;
}

Of course, we also need to upload a color to this new uniform location and this is easier than for attribute because the memory has already been allocated on the GPU (since the size is know and does not depend on the number of vertices).

loc = gl.glGetUniformLocation(program, "color")
gl.glUniform4f(loc, 0.0, 0.0, 1.0, 1.0)

If you run the new code/glut-quad-uniform-color.py example, you should obtain the blue quad as shown on the right.

Varying color

Figure

A colored quad using a per-vertex color.

Until now, we have been using a constant color for the four vertices of our quad and the result is (unsurprisingly) a boring uniform red or blue quad. We can make it a bit more interesting though by assigning different colors to each vertex and see how OpenGL will interpolate colors. Our new vertex shader would need to be rewritten as:

attribute vec2 position;
attribute vec4 color;
varying vec4 v_color;
void main()
{
  gl_Position = vec4(position, 0.0, 1.0);
  v_color= color;
}

We just added our new attribute color but we also added a new variable type: varying. This type is actually used to transmit a value from the vertex shader to the fragment shader. As you might have guessed, the varying type means this value won't be constant over the different fragments but will be interpolated depending on the relative position of the fragment in the triangle, as I explained in the Interpolation section. Note that we also have to rewrite our fragment shader accordingly, but now the v_color will be an input:

varying vec4 v_color;
void main()
{
  gl_FragColor = v_color;
}

We now need to upload vertex color to the GPU. We could create a new vertex dedicated buffer and bind it to the new color attribute, but there is a more interesting solution. We'll use instead a single numpy array and a single buffer, taking advantage of the NumPy structured array:

data = np.zeros(4, [("position", np.float32, 2),
                    ("color",    np.float32, 4)])
data['position'] = (-1,+1), (+1,+1), (-1,-1), (+1,-1)
data['color']    = (0,1,0,1), (1,1,0,1), (1,0,0,1), (0,0,1,1)

Our CPU data structure is thus:

           1ˢᵗ vertex              2ⁿᵈ vertex
   ┌───────────────────────┬───────────────────────┬┄
   ┌───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬─┄
   │ X │ Y │ R │ G │ B │ A │ X │ Y │ R │ G │ B │ A │ …
   └───┴───┴───┴───┴───┴───┴───┴───┴───┴───┴───┴───┴─┄
   ↑       ↑               └───────────────────────┘
 position  color                     stride
 offset    offset

Binding the buffer is now a bit more complicated but it is made relatively easy thanks to NumPy:

stride = data.strides[0]
offset = ctypes.c_void_p(0)
loc = gl.glGetAttribLocation(program, "position")
gl.glEnableVertexAttribArray(loc)
gl.glBindBuffer(gl.GL_ARRAY_BUFFER, buffer)
gl.glVertexAttribPointer(loc, 2, gl.GL_FLOAT, False, stride, offset)

offset = ctypes.c_void_p(data.dtype["position"].itemsize)
loc = gl.glGetAttribLocation(program, "color")
gl.glEnableVertexAttribArray(loc)
gl.glBindBuffer(gl.GL_ARRAY_BUFFER, buffer)
gl.glVertexAttribPointer(loc, 4, gl.GL_FLOAT, False, stride, offset)

Uniform color

Adding a uniform specified color is only a matter of modifying the fragment shader as in the previous section and directly assigning the color to the quad program (see code/chapter-03/glumpy-quad-uniform-color.py):

quad["color"] = 0,0,1,1

Varying color

Adding a per-vertex color is also only a matter of modifying the fragment shader as in the previous section and directly assigning the color to the quad program (see code/chapter-03/glumpy-quad-varying-color.py):

quad["color"] = (1,1,0,1), (1,0,0,1), (0,0,1,1), (0,1,0,1)

Homogeneous coordinates

Even though we're dealing with the three-dimensional Euclidean space, three dimensional coordinates are actually not the best representation we can use and this is the reason why we will use homogeneous coordinates that describe coordinates in a four-dimensional projective space (that includes the Euclidean space). We'll see in the next section that this allows us to express linear transformations (rotation, scaling), affine transformations (translations) and projection using 4×4 matrices. Homogeneous coordinatess are tightly linked with regular 3D coordinates with the noticeable difference that they require a fourth w coordinate that corresponds to the fourth dimension, let's call it the projective dimension. In order to explain it, we'll use a 1-dimensional space where point coordinates are single scalars indicating the position of the points onto the X-axis. This will make everything clearer hopefully.

Let us consider for example a simple set of points [-1.0, -0.5, 0.0, +0.5, +1.0] in this unidimensional space. We want to project onto another segment [-2,+2] that represents the screen (any point projected outside this segment is discared and won't be visible into the final projection). The question now is how do we project the points onto the screen?

Figure

5 sets of homogeneous coordinates, each of them corresponding to the set of unidimensional Cartesian coordinates [-1.0, -0.5, 0.0, +0.5, +1.0].

To answer this question, we need to know where is the camera (from where do we look at the scene) and where are the points positioned relatively to the screen. This is the reason why we introduce a supplementary w coordinate in order to indicate the distance to the screen. To go from our Euclidean representation to our new homogeneous representation, we'll use a conventional and default value of 1 for all the w such that our new point set is now [(-1.0,1.0), -(0.5,1.0), (0.0,1.0), (+0.5,1.0), (+1.0,1.0)]. Reciprocally, a point (x,w) in projective space corresponds to the point x/w (if w ≠ 0) in our unidimensional Euclidean space. From this conversion, we can see immediately that there exists actually an infinite set of homogenous coordinates that correspond to a single Cartesian coordinate as illustrated on the figure.

Figure

Different one dimensional projections using homogeneous coordinates.

Projections

We are now ready to project our point set onto the screen. As shown on the figure above, we can use an orthographic (all rays are parallel) or a linear projection (rays originate from the camera point and hit the screen, passing through points to be projected). For these two projections, results are similar but different. In the first case, distances have been exactly conserved while in the second case, the distance between projected points has increased, but projected points are still equidistant. The third projection is where homogenous coordinates make sense. For this (arbitrary) projection, we decided that the further the point is from the origin, the further away from the origin its projection will be. To do that, we measure the distance of the point to the origin and we add this distance to its w value before projecting it (this corresponds to the black circles on the figure) using the linear projection. It is to be noted that this new projection does not conserve the distance relationship and if we consider the set of projected points [P(-1.0), P(-0.5), P(0.0), P(+0.5), P(+1.0)], we have ║P(-1.0)-P(-0.5)]║ > ║P(-0.5)- P(0.0)║.

Back to our regular 3D-Euclidean space, the principle remains the same and we have the following relationship between Cartesian and homogeneous coordinates:

 (x,y,z,w) → (x/w, y/w, z/w) (for w ≠ 0)
Homogeneous    Cartesian

 (x,y,z)   → (x, y, z, 1)
Cartesian     Homogeneous

If you didn't understand everything, you can stick to the description provided by Sam Hocevar:

• If w = 1, then the vector (x,y,z,1) is a position in space
• If w = 0, then the vector (x,y,z,0) is a direction

Translation

Considering a vertex V = (x, y, z, 1) and a translation vector T = (tx, ty, tz, 0), the translation of V by T is (x+tx, y+ty, z+tz, 1). The corresponding matrix is given below:

┌          ┐   ┌   ┐   ┌                        ┐   ┌      ┐
│ 1 0 0 tx │ * │ x │ = │ 1*x + 0*y + 0*z + tx*1 │ = │ x+tx │
│ 0 1 0 ty │   │ y │   │ 0*x + 1*y + 0*z + ty*1 │   │ y+ty │
│ 0 0 1 tz │   │ z │   │ 0*x + 0*y + 1*z + tz*1 │   │ z+tz │
│ 0 0 0 1  │   │ 1 │   │ 0*x + 0*y + 0*z +  1*1 │   │ 1    │
└          ┘   └   ┘   └                        ┘   └      ┘

Scaling

Considering a vertex V = (x, y, z, 1) and a scaling vector S = (sx, sy, sz, 0), the scaling of V by S is (sx*x, sy*y, sz*z, 1). The corresponding matrix is given below:

┌            ┐   ┌   ┐   ┌                          ┐   ┌      ┐
│ sx 0  0  0 │ * │ x │ = │ sx*x +  0*y +  0*z + 0*1 │ = │ sx*x │
│ 0  sy 0  0 │   │ y │   │  0*x + sy*y +  0*z + 0*1 │   │ sy*y │
│ 0  0  sz 0 │   │ z │   │  0*x +  0*y + sz*z + 0*1 │   │ sz*z │
│ 0  0  0  1 │   │ 1 │   │  0*x +  0*y +  0*z + 1*1 │   │ 1    │
└            ┘   └   ┘   └                          ┘   └      ┘

Rotation

A rotation is defined by an axis of rotation A and an angle of rotation d. We defined below only the most common rotations, that is, around the X-axis, Y-axis and Z-axis.

┌                    ┐   ┌   ┐   ┌                                 ┐
│ 1   0       0    0 │ * │ x │ = │ 1*x      + 0*y      + 0*z + 0*1 │
│ 0 cos(d) -sin(d) 0 │   │ y │   │ 0*x + cos(d)*y - sin(d)*z + 0*1 │
│ 0 sin(d)  cos(d) 0 │   │ z │   │ 0*x + sin(d)*y + cos(d)*z + 0*1 │
│ 0   0       0    1 │   │ 1 │   │ 0*x      + 0*y +      0*z + 1*1 │
└                    ┘   └   ┘   └                                 ┘
                                 ┌                      ┐
                               = │ x                    │
                                 │ cos(d)*y - sin(d)*z  │
                                 │ sin(d)*y + cos(d)*z  │
                                 │ 1                    │
                                 └                      ┘

┌                    ┐   ┌   ┐   ┌                                  ┐
│  cos(d) 0 sin(d) 0 │ * │ x │ = │  cos(d)*x + 0*y + sin(d)*z + 0*1 │
│    0    1   0    0 │   │ y │   │       0*x + 1*y +      0*z + 0*1 │
│ -sin(d) 0 cos(d) 0 │   │ z │   │ -sin(d)*x + 0*y + cos(d)*z + 0*1 │
│    0    0    0   1 │   │ 1 │   │       0*x + 0*y      + 0*z + 1*1 │
└                    ┘   └   ┘   └                                  ┘
                                 ┌                      ┐
                               = │ cos(d)*x + sin(d)*z  │
                                 │ y                    │
                                 │ -sin(d)*x + cos(d)*z │
                                 │ 1                    │
                                 └                      ┘

┌                    ┐   ┌   ┐   ┌                                  ┐
│ cos(d) -sin(d) 0 0 │ * │ x │ = │  cos(d)*x - sin(d)*y + 0*z + 0*1 │
│ sin(d)  cos(d) 0 0 │   │ y │   │  sin(d)*x + cos(d)*y + 0*z + 0*1 │
│   0       0    1 0 │   │ z │   │       0*x +      0*y + 1*z + 0*1 │
│   0       0    0 1 │   │ 1 │   │       0*x +      0*y + 0*z + 1*1 │
└                    ┘   └   ┘   └                                  ┘
                                 ┌                      ┐
                               = │ cos(d)*x - sin(d)*y  │
                                 │ sin(d)*x + cos(d)*y  │
                                 │ z                    │
                                 │ 1                    │
                                 └                      ┘

A word of caution

OpenGL uses a column-major representation of matrices. This mean that when reading a set of 16 contiguous values in memory, relative to a 4×4 matrix, the first 4 values correspond to the first column while in Numpy (using C default layout), this would correspond to the first row. In order to stay consistent with most OpenGL tutorials, we'll use a column-major order in the rest of this book. This means that any glumpy transformations will appear to be transposed when displayed, but the underlying memory representation will still be consistent with OpenGL and GLSL. This is all you need to know at this stage.

Considering a set of 16 contiguous values in memory:

┌                                  ┐
│ a b c d e f g h i j k l  m n o p │
└                                  ┘

We get different representations depending on the order convention (column major or row major):

column-major                          row-major
  (OpenGL)                             (NumPy)
 ┌         ┐   ┌   ┐   ┌         ┐   ┌         ┐   ┌                   ┐
 │ a e i m │ × │ x │ = │ x y z w │ × │ a b c d │ = │ ax + ey + iz + mw │
 │ b f j n │   │ y │   └         ┘   │ e f g h │   │ bx + fy + jz + nw │
 │ c g k o │   │ z │                 │ i j k l │   │ cx + gy + kz + ow │
 │ d h l p │   │ w │                 │ m n o p │   │ dx + hy + lz + pw │
 └         ┘   └   ┘                 └         ┘   └                   ┘

For example, here is a translation matrix as returned by the glumpy.glm.translation function:

import glumpy
T = glumpy.glm.translation(1,2,3)
print(T)
[[ 1.  0.  0.  0.]
 [ 0.  1.  0.  0.]
 [ 0.  0.  1.  0.]
 [ 1.  2.  3.  1.]]
print(T.ravel())
[ 1.  0.  0.  0.  0.  1.  0.  0.  0.  0.  1.  0.  1.  2.  3.  1.]
                                                  ↑   ↑   ↑
                                                  13  14  15

So this means you would use this translation on the left when uploaded to the GPU, but you would use on the right with Python/NumPy:

T = glumpy.glm.translation(1,2,3)
V = [3,2,1,1]
print(np.dot(V, T))
[ 4.  4.  4.  1.]

Orthographic

┌                                         ┐ n: near
│ 2/(r-l)    0       0     -((r+l)/(r-l)) │ f: far
│   0     2/(t-b)    0     -((t+b)/(t-b)) │ t: top
│   0        0    -2/(f-n) -((f+n)/(f-n)) │ b: bottom
│   0        0      -1            0       │ l: left
└                                         ┘ r: right
          Orthographic projection

Perspective

┌                                               ┐ n: near
│ 2n/(r-l)    0       (r+l)/(r-l)       0       │ f: far
│    0     2n/(t-b)   (t+b)/(t-b)       0       │ t: top
│    0        0     -((f+n)/(f-n)) -(2nf/(f-n)) │ b: bottom
│    0        0           -1            0       │ l: left
└                                               ┘ r: right
            Perspective projection

At this point, it is not necessary to understand how these matrices were built. Suffice it to say they are standard matrices in the 3D world. Both assume the viewer (=camera) is located at position (0,0,0) and is looking in the direction (0,0,1).

There exists a second form of the perpective matrix that might be easier to manipulate. Instead of specifying the right/left/top/bottom planes, we'll use field of view in the horizontal and vertical direction:

┌                                     ┐ n: near
│ c/aspect  0       0          0      │ f: far
│    0      c       0          0      │ c: cotangent(fovy)
│    0      0  (f+n)/(n-f)  2nf/(n-f) │
│    0      0      -1          0      │
└                                     ┘
            Perspective projection

where fovy specifies the field of view angle in the y direction and aspect specifies the aspect ratio that determines the field of view in the x direction.

Model and view matrices

We are almost done with matrices. You may have guessed that the above matrices require the viewing volume to be in the z direction. We could design our 3D scene such that all objects are within this direction but it would not be very convenient. So instead, we use a view matrix that maps the world space to camera space. This is pretty much as if we were orienting the camera at a given position and look toward a given direction. In the meantime, we can further refine the whole pipeline by providing a model matrix that maps the object's local coordinate space into world space. For example, this is useful for rotating an object around its center. To sum up, we need:

• Model matrix maps from an object's local coordinate space into world space
• View matrix maps from world space to camera space
• Projection matrix maps from camera to screen space

This corresponds to the model-view-projection model. If you have read the whole chapter carefully, you may have guessed the corresponding GLSL shader:

uniform mat4 view;
uniform mat4 model;
uniform mat4 projection;
attribute vec3 P;
void main(void)
{
    gl_Position = projection*view*model*vec4(P, 1.0);
}

Colored cube

The previous cube is not very interesting because we used a single color for all the faces and this tends to hide the 3D structure. We can fix this by adding some colors and in the process, we'll discover why glumpy is so useful. To add color per vertex to the cube, we simply define the vertex structure as:

V = np.zeros(8, [("position", np.float32, 3),
                 ("color",    np.float32, 4)])
V["position"] = [[ 1, 1, 1], [-1, 1, 1], [-1,-1, 1], [ 1,-1, 1],
                 [ 1,-1,-1], [ 1, 1,-1], [-1, 1,-1], [-1,-1,-1]]
V["color"]    = [[0, 1, 1, 1], [0, 0, 1, 1], [0, 0, 0, 1], [0, 1, 0, 1],
                 [1, 1, 0, 1], [1, 1, 1, 1], [1, 0, 1, 1], [1, 0, 0, 1]]

And we're done ! Well, actually, we also need to slightly modify the vertex shader since color is now an attribute that needs to be passed to the fragment shader.

vertex = """
uniform mat4   model;         // Model matrix
uniform mat4   view;          // View matrix
uniform mat4   projection;    // Projection matrix
attribute vec4 color;         // Vertex color
attribute vec3 position;      // Vertex position
varying vec4   v_color;       // Interpolated fragment color (out)
void main()
{
    v_color = color;
    gl_Position = projection * view * model * vec4(position,1.0);
} """

fragment = """
varying vec4 v_color;         // Interpolated fragment color (in)
void main()
{
    gl_FragColor = v_color;
} """

Furthermore, since our vertex buffer fields corresponds exactly to program attributes, we can directly bind it:

cube = gloo.Program(vertex, fragment)
cube.bind(V)

But we could also have written

cube = gloo.Program(vertex, fragment)
cube["position"] = V["position"]
cube["color"] = V["color"]

Complete source code: code/chapter-05/color-cube.py

Outlined cube

We can make the cube a bit nicer by outlining it using black lines. To outline the cube, we need to draw lines between pairs of vertices on each face. 4 lines for the back and front face and 2 lines for the top and bottom faces. Why only 2 lines for top and bottom? Because lines are shared between the faces. So overall we need 12 lines and we need to compute the corresponding indices (I did it for you):

O = [0,1, 1,2, 2,3, 3,0,
     4,7, 7,6, 6,5, 5,4,
     0,5, 1,6, 2,7, 3,4 ]
O = O.view(gloo.IndexBuffer)

We then need to draw the cube twice. One time using triangles and the indices index buffer and one time using lines with the outline index buffer. We need also to add some OpenGL black magic to make things nice. It's not very important to understand it at this point but roughly the idea to make sure lines are drawn "above" the cube because we paint a line on a surface:

@window.event
def on_draw(dt):
    global phi, theta, duration

    window.clear()

    # Filled cube
    gl.glDisable(gl.GL_BLEND)
    gl.glEnable(gl.GL_DEPTH_TEST)
    gl.glEnable(gl.GL_POLYGON_OFFSET_FILL)
    cube['ucolor'] = .75, .75, .75, 1
    cube.draw(gl.GL_TRIANGLES, I)

    # Outlined cube
    gl.glDisable(gl.GL_POLYGON_OFFSET_FILL)
    gl.glEnable(gl.GL_BLEND)
    gl.glDepthMask(gl.GL_FALSE)
    cube['ucolor'] = 0, 0, 0, 1
    cube.draw(gl.GL_LINES, O)
    gl.glDepthMask(gl.GL_TRUE)

    # Rotate cube
    theta += 1.0 # degrees
    phi += 1.0 # degrees
    model = np.eye(4, dtype=np.float32)
    glm.rotate(model, theta, 0, 0, 1)
    glm.rotate(model, phi, 0, 1, 0)
    cube['model'] = model

Complete source code: code/chapter-05/outlined-cube.py

Textured cube

For making a textured cube, we need a texture (a.k.a. an image) and some coordinates to tell OpenGL how to map it to the cube faces. Texture coordinates are normalized and should be inside the [0,1] range (actually, texture coordinates can be pretty much anything but for the sake of simplicity, we'll stick to the [0,1] range). Since we are displaying a cube, we'll use one texture per side and the texture coordinates are quite easy to define: [0,0], [0,1], [1,0] and [1,1]. Of course, we have to take care of assigning the right texture coordinates to the right vertex or you texture will be messed up.

Furthemore, we'll need some extra work because we cannot share anymore our vertices between faces since they won't share their texture coordinates. We thus need to have a set of 24 vertices (6 faces × 4 vertices). We'll use the dedicated function below that will take care of generating the right texture coordinates.

def cube():
    vtype = [('position', np.float32, 3),
             ('texcoord', np.float32, 2)]
    itype = np.uint32

    # Vertices positions
    p = np.array([[1, 1, 1], [-1, 1, 1], [-1, -1, 1], [1, -1, 1],
                  [1, -1, -1], [1, 1, -1], [-1, 1, -1], [-1, -1, -1]],
                  dtype=float)

    # Texture coords
    t = np.array([[0, 0], [0, 1], [1, 1], [1, 0]])

    faces_p = [0, 1, 2, 3,  0, 3, 4, 5,   0, 5, 6, 1,
               1, 6, 7, 2,  7, 4, 3, 2,   4, 7, 6, 5]
    faces_t = [0, 1, 2, 3,  0, 1, 2, 3,   0, 1, 2, 3,
               3, 2, 1, 0,  0, 1, 2, 3,   0, 1, 2, 3]

    vertices = np.zeros(24, vtype)
    vertices['position'] = p[faces_p]
    vertices['texcoord'] = t[faces_t]

    filled = np.resize(
       np.array([0, 1, 2, 0, 2, 3], dtype=itype), 6 * (2 * 3))
    filled += np.repeat(4 * np.arange(6, dtype=itype), 6)

    vertices = vertices.view(gloo.VertexBuffer)
    filled   = filled.view(gloo.IndexBuffer)

    return vertices, filled

Now, inside the fragment shader, we have access to the texture:

vertex = """
uniform mat4   model;      // Model matrix
uniform mat4   view;       // View matrix
uniform mat4   projection; // Projection matrix
attribute vec3 position;   // Vertex position
attribute vec2 texcoord;   // Vertex texture coordinates
varying vec2   v_texcoord;   // Interpolated fragment texture coordinates (out)

void main()
{
    // Assign varying variables
    v_texcoord  = texcoord;

    // Final position
    gl_Position = projection * view * model * vec4(position,1.0);
} """


fragment = """
uniform sampler2D texture; // Texture
varying vec2 v_texcoord;   // Interpolated fragment texture coordinates (in)
void main()
{
    // Get texture color
    gl_FragColor = texture2D(texture, v_texcoord);
} """

Complete source code: code/chapter-05/textured-cube.py

Sample based methods

Figure

Rendering of a disc using multisample anti-aliasing. By mutisampling the same pixel, the edge can be rendered smoother. The thick black square pixel that was previously considered outside the shape can not be considered half-outside (2 samples outside) and half-inside (2 samples inside).

One of the simplest method to remove antialising consists in using several samples to estimate the final color of a fragment. Instead of only considering the center of the pixel, one case use several samples over the whole surface of a pixel in order to have a better estimate as shown on the figure on the right. A fragment that was previously considered outside, based on its center only, can now be considered half inside / half outside. This multi-sampling helps to attenuate the jagged edges we've seen in the previous section. On the figure, we used a very simple and straightforward multi-sample method, assigning fixed and equidistant locations for the four subsamples. There exist however better methods for multi-sampling as shown by the impressive list of sample based antialiasing techniques:

• SSAA: Supersampling antialiasing
• MSAA: Multisample antialiasing
• FSAA: Full screen anti-aliasing
• FXAA: Fast approximate antialiasing
• SMAA: Subpixel morphological antialiasing
• DLAA: Directionally localized antialiasing
• NFAA: Normal filter antialiasing
• HRAA: High-Resolution antialiasing
• TXAA: Temporal antialiasing
• EQAA: Enhanced quality at-prntialiasing
• CSAA: Coverage Sample antialiasing

Depending on the performance you need to achieve (in terms of rendering quality and speed) one method might be better than the other. However, you cannot expect to achieve top quality due to inherent limitations of all these methods. If they are great for real-time rendering such as video games (and some of them are really good), they are hardly sufficient for any scientific visualization as illustrated on the figure below.

Figure

Supersampling applied to a rasterized triangle, using various sub-pixel patterns. The more samples , the better the output but even using 64 samples, the rendering quality does not match the SDF rendering, especially if you consider triangle sharp vertices. Supersampled triangles have been rendered using a dedicated shader (see code/chapter-06/triangle-ssaa.py) and the SDF triangle has been rendered using a fake signed-distance triangle function (see below) and a stroke anti-alias function (see code/chapter-06/triangle-sdf.py)

This is the reason why we won't use them in the rest of this book. If you want more details on these techniques, you can have a look at this reddit discussion explaining antialiasing modes or this nice overview of MSAA

Coverage methods

Figure

Rendering of a disc using exact coverage anti-aliasing.

Another approach for anti-aliasing is to compute the exact coverage of a shape over each pixel surface as shown on the figure on the right. To do so, we need of course to know precisely the shape we want to display and where it is located in order to compute the coverage of the shape onto the pixel grid. On the image, this corresponds to the grey areas that give us direct access to the final color of the pixel (more precisely, the percentage of the color we have to mix with the background color or any other object in the vicinity). Unfortunately, such method is not possible to enforce in a full 3D scene because all the transformations and differen occlusuins would make the computation of the final shape too complex. In two dimensions however, this is probably the best method we can use and this is also the method that is used in the Anti-grain geometry library that constitutes the quality standard we aim at.

Figure

Actual versus approximated coverage.

But even in 2D, computing the exact coverage of the shape over the different pixels can become rapidly a complex and slow task. One way to greatly simplify the problem is to consider pixel to be round (instead of square or rectangle). With such asumption, we only need to compute the distance from the center of the pixel to the border of the shape (that is locally considered to be a line) to get a very accurate estimate of the coverage and this is exactly what we'll do in the next section.

If you wonder if our round pixel shape approximation makes any. sense at all, have a look at the subpixel zoo maintained by Ian Mallett and you'll understand our assumption is not so bad overall.

Signed distance function

From wikipedia (again):

A signed distance function (or oriented distance function) of a set Ω in a metric space determines the distance of a given point x from the boundary of Ω, with the sign determined by whether x is in Ω. The function has positive values at points x inside Ω, it decreases in value as x approaches the boundary of Ω where the signed distance function is zero, and it takes negative values outside of Ω.

Said differently: in order to render a shape, we need to find a function of x and y that returns a value that is the signed distance to the shape, that is, a signed distance to the border of the shape. Inside the shape, the value is positive, outside the shape the value is negative and on the border, the value is zero. Easy enough.

Of course, the question is now how do we find such function? Let's start with the most simple geometrical primitive: a circle centered on (xc,yc) with a radius r. For any point (x,y), we know the (non-negative) distance to the center is given by: d = sqrt((x-xc)*(x-xc)+(y-yc)*(y-yc)). To simplify computations, we'll consider the circle to centered on the origin, the distance now writes d = sqrt(x*x+y*y). This distance is not what we want since we target a signed distance to the border of the circle. However, this can be obtained very easily by subtracting the radius r from d(x,y). In the end, signed distance from a point (x,y) to a circle of radius r centered on the origin is given by:

d(x,y) = sqrt(x*x+y*y) - r

Figure

Signed distance to a circle. Inside is red, outside is blue, border is white.

See code/chapter-06/circle-sdf-distances.py

As an exercise, you can check that d(x,y) is zero if (x,y) is on the border, strictly negative if (x,y) is inside the circle and strictly positive outside the circle.

Now, let's check if OpenGL is consistent with our maths. We'll write a fragment shader that compute the color according to the distance to the shape. We'll use the blue color outside the circle, red color inside and white color on the border (with some tolerance or we won't see anything).

float distance(vec2 P, vec2 center, float radius)
{
    return length(P-center) - radius;
}

varying vec2 v_position;
void main()
{
    const float epsilon = 0.005;
    float d = distance(v_position, vec2(0.0), 0.5);
    if (d > +epsilon)
        gl_FragColor = vec4(abs(d), 0.0, 0.0, 1.0);
    else if (d < -epsilon)
        gl_FragColor = vec4(0.0, 0.0, abs(d), 1.0);
    else
        gl_FragColor = vec4(0.0, 0.0, 0.0, 1.0);
}

Geometrical primitives

We need now to define a few primitives usigned signed distance function. You'll understand in the next section why we only need a few primitives. In the meantime, we'll use a less boring palette than the one in the previous section. We'll use instead the palette that has become the standard for displaying SDF on Shadertoy (it has been designed by Íñigo Quílez to the best of my knowledge):

vec4 color(float d)
{
    vec3 white = vec3(1.0, 1.0, 1.0);
    vec3 blue  = vec3(0.1, 0.4, 0.7);
    vec3 color = white - sign(d)*blue;
    color *= (1.0 - exp(-4.0*abs(d))) * (0.8 + 0.2*cos(140.0*d));
    color = mix(color, white, 1.0-smoothstep(0.0,0.02,abs(d)) );
    return vec4(color, 1.0);
}

However, we don't want to copy this code in all the example. We can instead write a palette.glsl shader and include it in each of the example.

float SDF_circle(vec2 p, float radius)
{
    return length(p) - radius;
}

float SDF_plane(vec2 p, vec2 p0, vec2 p1)
{
  vec2 T = p1 - p0;
  vec2 O = normalize(vec2(T.y, -T.x));
  return dot(O, p0 - p);
}

// Code by Inigo Quilez
// See https://www.shadertoy.com/view/4llXD7
float SDF_box(vec2 p, vec2 size)
{
     vec2 d = abs(p) - size;
     return min(max(d.x,d.y),0.0) + length(max(d,0.0));
}

// Code derived from the true triangle code by Inigo Quilez
// See https://www.shadertoy.com/view/4llXD7
float SDF_round_box(vec2 p, vec2 size, float radius)
{
    return SDF_box(p, size) - radius;
}

float SDF_fake_box(vec2 p, vec2 size)
{
    return max(abs(p.x)-size.x, abs(p.y)-size.y);
}

// Code by Inigo Quilez
// See https://www.shadertoy.com/view/XsXSz4
float SDF_triangle(vec2 p, vec2 p0, vec2 p1, vec2 p2)
{
    vec2 e0 = p1 - p0;
    vec2 e1 = p2 - p1;
    vec2 e2 = p0 - p2;

    vec2 v0 = p - p0;
    vec2 v1 = p - p1;
    vec2 v2 = p - p2;

    vec2 pq0 = v0 - e0*clamp( dot(v0,e0)/dot(e0,e0), 0.0, 1.0 );
    vec2 pq1 = v1 - e1*clamp( dot(v1,e1)/dot(e1,e1), 0.0, 1.0 );
    vec2 pq2 = v2 - e2*clamp( dot(v2,e2)/dot(e2,e2), 0.0, 1.0 );

    float s = sign( e0.x*e2.y - e0.y*e2.x );
    vec2 d = min( min(
          vec2( dot( pq0, pq0 ), s*(v0.x*e0.y-v0.y*e0.x) ),
          vec2( dot( pq1, pq1 ), s*(v1.x*e1.y-v1.y*e1.x) )),
          vec2( dot( pq2, pq2 ), s*(v2.x*e2.y-v2.y*e2.x) ));
    return -sqrt(d.x)*sign(d.y);
}

// Code derived from the true triangle code by Inigo Quilez
// See https://www.shadertoy.com/view/XsXSz4
float SDF_round_triangle(vec2 p, vec2 p0, vec2 p1, vec2 p2, float radius)
{
    return SDF_triangle(p, p0, p1, p2) - radius;
}

float SDF_fake_triangle(vec2 p, vec2 p0, vec2 p1, vec2 p2)
{
    vec2 e0 = p1 - p0;
    vec2 e1 = p2 - p1;
    vec2 e2 = p0 - p2;

    vec2 v0 = p - p0;
    vec2 v1 = p - p1;
    vec2 v2 = p - p2;

    vec2 o0 = normalize(vec2(e0.y, -e0.x));
    vec2 o1 = normalize(vec2(e1.y, -e1.x));
    vec2 o2 = normalize(vec2(e2.y, -e2.x));

    return max(max(dot(o0,v0), dot(o1,v1)), dot(o2,v2));
}

// Code by Inigo Quilez
// See https://www.shadertoy.com/view/4sS3zz
float SDF_ellipse(vec2 p, vec2 ab)
{
    // The function does not like circles
    if (ab.x == ab.y) ab.x = ab.x*0.9999;

    p = abs( p ); if( p.x > p.y ){ p=p.yx; ab=ab.yx; }
    float l = ab.y*ab.y - ab.x*ab.x;
    float m = ab.x*p.x/l;
    float n = ab.y*p.y/l;
    float m2 = m*m;
    float n2 = n*n;
    float c = (m2 + n2 - 1.0)/3.0;
    float c3 = c*c*c;
    float q = c3 + m2*n2*2.0;
    float d = c3 + m2*n2;
    float g = m + m*n2;
    float co;

    if( d<0.0 ) {
        float p = acos(q/c3)/3.0;
        float s = cos(p);
        float t = sin(p)*sqrt(3.0);
        float rx = sqrt( -c*(s + t + 2.0) + m2 );
        float ry = sqrt( -c*(s - t + 2.0) + m2 );
        co = ( ry + sign(l)*rx + abs(g)/(rx*ry) - m)/2.0;
    } else {
        float h = 2.0*m*n*sqrt( d );
        float s = sign(q+h)*pow( abs(q+h), 1.0/3.0 );
        float u = sign(q-h)*pow( abs(q-h), 1.0/3.0 );
        float rx = -s - u - c*4.0 + 2.0*m2;
        float ry = (s - u)*sqrt(3.0);
        float rm = sqrt( rx*rx + ry*ry );
        float p = ry/sqrt(rm-rx);
        co = (p + 2.0*g/rm - m)/2.0;
    }
    float si = sqrt( 1.0 - co*co );
    vec2 r = vec2( ab.x*co, ab.y*si );
    return length(r - p ) * sign(p.y-r.y);
}

// Code by Inigo Quilez
// See https://www.shadertoy.com/view/MdfGWn
float SDF_fake_ellipse(vec2 p, vec2 size)
{
    float r = 0.2;
    float f = length( p*size );
    f = length(p*size);
    return f*(f-r)/length(p*size*size);
}

Raw points

Figure

"Point" as drawn by OpenGL (see points.py).

The most straightforward way to display points is to use the gl.GL_POINTS primitive that displays a quad that is always facing the camera (i.e. billboard). This is very convenient because a mathematical point has no dimension, even though we'll use this primitive to draw discs and circles in the next section. The size of the quad must be specified within the vertex shader using the gl_PointSize variable (note that the size is expressed in pixels). As shown on the figure, the result is quite ugly.

import numpy as np
from glumpy import app, gloo, gl

vertex = """
  attribute vec2 position;
  void main() {
      gl_PointSize = 5.0;
      gl_Position = vec4(position, 0.0, 1.0);
  } """

fragment = """
  void main() {
       gl_FragColor = vec4(vec3(0.0), 1.0);
  } """

window = app.Window(512, 512, color=(1,1,1,1))
points = gloo.Program(vertex, fragment, count=1000)
points["position"] = np.random.uniform(-1,1,(len(points),2))

@window.event
def on_draw(dt):
    window.clear()
    points.draw(gl.GL_POINTS)

app.run()

Antialiased points

For drawing antialiased point, the size of the quad must be slighlty larger than the actual diameter of the point because we need some extra space for the antialias area. Considering a point with a radius r, the size of the quad is thus 2+ceil(2*r) if we consider using 1 pixel for the antalias area. Finally, considering a point centered at center with radius radius, our vertex shader reads (see also previous chapter on signed distance field):

// Screen resolution as (width, height)
uniform vec2 resolution;

// Point center (in pixel coordinates)
attribute vec2 center;

// Point radius (in pixels)
attribute float radius;

varying vec2 v_center;
varying float v_radius;
void main()
{
    v_radius = radius;
    v_center = center;
    gl_PointSize = 2.0 + ceil(2.0*radius);
    gl_Position = vec4(2.0*center/resolution-1.0, 0.0, 1.0);
}

You may have noticed that we gave the window resolution to the shader using a uniform (that will be updated each time the window size has changed). The goal is to be able to use window coordinates (i.e. pixels) from within Python without taking care of the normalized device coordinate (this transformation has been done in the vertex shader above). We now have one problem to solve. A GL point is made from a single vertex and the apparent size of the resulting quad is controlled by the gl_PointSize variable resulting in several fragments. How things are interpolated between vertices knowing there is ony one vertex? The answer is that there is no interpolation. If we want to know the position of a fragment relatively to the center, we have to find it ourself. Luckily, there is one interesting variable gl_FragCoord that gives us the absolute coordinate of the fragment in window coordinates (bottom-left is (0,0)). Subtracting the center from this coordinate will give us the relative position of the fragment from which we can compute the distance to the outer border of the point. Finally, our fragment shader reads:

varying vec2 v_center;
varying float v_radius;
void main()
{
    vec2 p = gl_FragCoord.xy - v_center;
    float a = 1.0;
    float d = length(p) - v_radius;
    if(d > 0.0) a = exp(-d*d);
    gl_FragColor = vec4(vec3(0.0), a);
}

Last, we setup our python program to display some discs:

V = np.zeros(16, [("center", np.float32, 2),
                  ("radius", np.float32, 1)])
V["center"] = np.dstack([np.linspace(32, 512-32, len(V)),
                         np.linspace(25, 28, len(V))])
V["radius"] = 15

window = app.Window(512, 50, color=(1,1,1,1))
points = gloo.Program(vertex, fragment)
points.bind(V.view(gloo.VertexBuffer))

@window.event
def on_resize(width, height):
    points["resolution"] = width, height

@window.event
def on_draw(dt):
    window.clear()
    points.draw(gl.GL_POINTS)

app.run()

Figure

Discs positionned vertically with a 0.2 pixel increase. See discs-aligned.py

Figure

Circles positionned vertically with a 0.2 pixel increase. See circles-aligned.py

You can see the result on the image on the right. Not only the discs are properly antialiased, but they are also positionned at the subpixel level. In the image on the right, each disc is actually vertically shifted upward by 0.2 pixels compared to its left neightbour. However, you cannot see any artefacts (can you?): the discs are similar and properly aligned. For the disc outlines, we simply have to get the absolute distance instead of the signed distance.

varying vec2 v_center;
varying float v_radius;
void main()
{
    vec2 p = gl_FragCoord.xy - v_center;
    float a = 1.0;
    float d = length(p) - v_radius;
    if(abs(d) > 0.0) a = exp(-d*d);
    gl_FragColor = vec4(vec3(0.0), a);
}

Flat sphere

Figure

A lit sphere

If you look closely at a sphere, you'll see that that the projected shape on screen is actualy as disc as shown on the figure on the right. This is actually true independently of the viewpoint and we can take advantage of it. A long time ago (with the fixed pipeline), rendering a sphere meant tesselating the sphere with a large number of triangles. The larger the number of triangles, the higher the quality of the sphere and the slower the rendering. However, with the advent of shaders, things have changeg dramatically and we can use fake spheres, i.e. discs thar are painted such as to appear as spheres. This is known as "impostors". If you look again at the image, you might realize that the appeareance of the sphere is given by the shading that is not uniform and suggests instead a specific lighting that seems to come from the upper right corner. Let's see if we can reproduce this.

Figure

A black disc (sphere-1.py)

First thing first, Let's setup a scene in order to display a single and large disc. To do that, we simply test if a fragment is inside or outside the circle:

varying vec2 v_center;
varying float v_radius;
void main()
{
    vec2 p = gl_FragCoord.xy - v_center;
    float z = 1.0 - length(p)/v_radius;
    if (z < 0.0) discard;
    gl_FragColor = vec4(vec3(0.0), 1.0);
}

Figure

Sphere normals view on the xz plane.

To simulate lighting on the disc, we need to compute normal vectors over the surface of the sphere (i.e. disc). Luckily enough for us, computing the normal for a sphere is very easy. We can simply use the p=(x,y) coordinates inside the fragment shader and compute the z coordinate. How? you might ask yourself. This is actually correlated to the distance d to the center such that z = 1-d. If you want to convice yourself, just look at the figure on the right that show a side view of half a sphere on the xz plane. The z coordinate is maximal in the center and null on the border.

We're ready to simulate lighting on our disc using the Phong model. I won't give all the detail now because we'll see that later. However, as you can see on the source below, this is quite easy and the result is flawless.

Figure

A fake lit sphere (sphere-3.py)

varying vec2 v_center;
varying float v_radius;
void main()
{
    vec2 p = (gl_FragCoord.xy - v_center)/v_radius;
    float z = 1.0 - length(p);
    if (z < 0.0) discard;

    vec3 color = vec3(1.0, 0.0, 0.0);
    vec3 normal = normalize(vec3(p.xy, z));
    vec3 direction = normalize(vec3(1.0, 1.0, 1.0));
    float diffuse = max(0.0, dot(direction, normal));
    float specular = pow(diffuse, 24.0);
    gl_FragColor = vec4(max(diffuse*color, specular*vec3(1.0)), 1.0);
}

True sphere

Figure

A bunch of fake spheres.

We can use this technique to display several "spheres" having different sizes and positions as shown on the figure on the right. This can be used to represent molecules for examples. Howewer, we have a problem with sphere intersecting each other. If you look closely the figure, you might have notices that no sphere intersect any sphere. This is due to the depth testing of the unique vertex (remember gl.GL_POINTS) that is used to generate the quad fragments. Each of these fragments share the same z coordinate resulting in having sphre fully in front of another of fully behind another. For accurate rendering, we thus have to tell OpenGL what is the depth of each fragment using the gl_FragDepth variable (that must be between 0 and 1):

Figure

A bunch of fake spheres with correct intersections (spheres.py).

varying vec3 v_center;
varying float v_radius;
void main()
{
    vec2 p = (gl_FragCoord.xy - v_center.xy)/v_radius;
    float z = 1.0 - length(p);
    if (z < 0.0) discard;

    gl_FragDepth = 0.5*v_center.z + 0.5*(1.0 - z);

    vec3 color = vec3(1.0, 0.0, 0.0);
    vec3 normal = normalize(vec3(p.xy, z));
    vec3 direction = normalize(vec3(1.0, 1.0, 1.0));
    float diffuse = max(0.0, dot(direction, normal));
    float specular = pow(diffuse, 24.0);
    gl_FragColor = vec4(max(diffuse*color, specular*vec3(1.0)), 1.0);
}

You can see on the figures that the spheres now intersect each other correctly.

Tiny discs

Figure

Disc spiral

Adapting the shader from the "Dots, discs, circles" section, try to write a script to draw discs on a spiral as displayed on the figure on the right. Be careful with small discs, especially when the radius is less than one pixel. In such case, you'll have to find a convincing way to suggest the size of the disc...

Solution: spiral.py

Antialiased triangles

Try to adapt the code from the ellipses section to remake the animation on the right. Be careful with the computation of the bouding box.

Solution: triangles.py

GPU Voronoi

Figure

A voronoi diagram computed on the GPU.

We've seen when rendering sphere that the individual depth of each fragment can be controled withing the fragment shader and we computed this depth by taking the distance to the center of each disc/sphere. The goal of this exercise is thus to adapt this method to render a Voronoi diagram as shown on the right.

Solution: voronoi.py

Quiver plot

Figure

An dynamic quiver plot made of two triangles.

Now that we know how to draw arrows, we can make a quiver plot very easily. The obvious solution would be to draw n arrows using 2×n triangles (since one arrow is two triangle). However, if your arrows are evenly spaced as on the figure on the right, there is a smarter solution using only two triangles.

Solution: quiver.py

Light and shadows

As explained before, the shadertoy website is a great resource and you can learn a lot by reading the sources accompanying each demo. As an exercise, have a look at this wonderful demo by Marteen that shows two dimensional signed distance field functions with light and shadows.

Simply gorgeous...