In the last two years and half, I have been working for Meteo Consult on an internal application running on a Mac, to create 3D meteorological maps, broadcast on the TV channel La Chaîne Météo.

A sphere, colour-mapped with NASA’s Blue Marble Next Generation. The Blue Marble measures 86400 by 43200 pixels, which is too big for most hardware, and too long to load anyway.

A detailed texture is needed

One major problem we have since the beginning is how to cover the Earth with a texture, since the texture has to be huge to be detailed enough. Currently we are stuck with a smaller texture, which presents two drawbacks:

The Earth is hardly detailed enough, so the minimum altitude of the camera is limited. For example, La Martinique or La Guadeloupe are only small blurry dots. We currently rely on 2D maps instead.
Even with its low resolution, the texture takes a lot of time to load on the GPU; about 3 s on my MacBook Pro 2013.

Hopefully, the application runs on a Mac Pro, which has a lot of GPU RAM; but even if we could load a big texture, GPU generally don’t handle textures larger than 16384 pixels, so we would be stuck anyway.

Tiles

Probably the solution was obvious to you: use tiles, Boy ! Of course, we thought of that since the beginning of the project, and I even tried to make something work, but to no avail. The major difficulty was to determine which tiles were visible. It’s rather easy on a flat 2D surface, but I could not find a reliable solution for the round 3D surface of the Earth.

Megatexture

Megatexture, also known as “Sparse Virtual Texture”, is a technique to compose a big virtual texture using tiles. The term was coined by John Carmack, who imagined this technique. I’ll stick with “Megatexture” since it sounds much cooler than “Virtual Texture”.

The virtual texture is made of tiles, which number varies with the Mipmap level.

Determining visible tiles

The great insight is how visible tiles are determined: the scene is rendered to an offscreen buffer, with a special fragment shader. In my case, the Megatexture is at most 256 by 256 tiles, and has a maximum of 8 mipmap levels, so the shader stores the tile’s x in the red component, the tile’s y in the green component, the mipmap level in the blue component, and the Texture ID in the alpha component. The scene is rendered to a RGBA8 offscreen buffer.

Texture ID

There may be several megatextures in a same scene. The texture ID permits to differentiate them in the Cache and in the Indirection Table later. Since objects which are not megatextured won’t be processed in the shader, you need to reserve a special Texture ID to mean “No texture”. It must corresponds to the buffer’s clear color, therefore I advise you choose the value 0x00, so it corresponds to a transparent color (since the texture Id is saved to the alpha channel).

Tiles determination shader

I’m sorry but I can’t provide my own code, so I’ll give you Sean Barrett’s instead, who was a pioneer in the technique, and made his code public:

    const float readback_reduction_shift = 2.0;
    const float vt_dimension_pages = 128.0;
    const float vt_dimension = 32768.0;
    uniform float mip_bias;
    // analytically calculates the mipmap level similar to what OpenGL does
    float mipmapLevel(vec2 uv, float textureSize)  {
        vec2 dx = dFdx(uv * textureSize);
        vec2 dy = dFdy(uv * textureSize);
        float d = max(dot(dx, dx), dot(dy, dy));
        return 0.5 * log2(d)   // explanation: 0.5*log(x) = log(sqrt(x)) + mip_bias - readback_reduction_shift;
    }

This first part determines the mipmap level. The formula is copied straight from OpenGL’s implementation, so everyone uses the same.

vt_dimensions_pages is the size of a Tile (what Barrett and a number of people call a “Page”, but which I find inappropriate).
vt_dimension is the size the megatexture at the most detailed level (mipmap 0).
you’ll see below that the CPU is going to read the pixels of the offscreen buffer. To save a lot of processing power, the scene is not rendered at full size. readback_reduction_shift is a power of two; since it equals to 2 here, the offscreen buffer is a quarter of the width and height of the final rendering. I personally set this value to 4, and set the width and height of the buffer to 1/16th of the size of my view.
I’m not sure what mip_bias is. I believe this is a way to make the shader less agressive in its changes of mipmap levels, at the cost of the texture being a little blurry at times. (I don’t use it my own implementation).

The second part determines the Tile’s x and y and renders them in the color buffer:

    void main() {
        // the tile x/y coordinates depend directly on the virtual texture coords
        gl_FragColor.rg = floor(gl_TexCoord[0].xy * vt_dimension_pages) / 255.0;
        gl_FragColor.b = mipmapLevel(gl_TexCoord[0].xy, vt_dimension) / 255.0;
        gl_FragColor.a = 1.0; // BGRA: mip, x, y, 255
    }

Note that there is a mistake here: the mipmap level must be floored! Otherwise there will be a discrepancy between the level determined here, and the one determined in the Texturing shader.

Result

A small image:

I changed the way colours are rendered so the image is visible, but the size is real. If the OpenGL view renders at 800 x 600, then the offset buffer is rendered at 1/16th of that, that is 50 x 37.

Loading tiles in the Cache

Reading back the offscreen buffer

I use glReadPixels to get the pixels. Every pixel corresponds to what I call a “Tile Id”: texture Id, x, y, mipmap level. Pixels with the “None” texture ID are discarded immediately. Others are converted to TileId objects, which are added to a NSMutableSet, in order to remove duplicates: since a same tile shows at several places, its TileId will appear several times.

It is not necessary to read the buffer, and therefore determine visible tiles, at every draw cycle. I do it only once every 4 frames (at 60fps = every 15th of second).

Determine the tiles to load

Now we have a list of visible tiles, we can compare them to the ones already in the Cache. The difference is the tiles to load.

In my implementation, tiles are loaded as background tasks, but textures are loaded in GPU memory on the main thread, because we have to with OpenGL. This textures obviously don’t have mip maps, but do use interpolation.

While the tile is loading, you will like to replace it with a “parent” tile— one with lower details — already in the cache. This is not too difficult, since the replacement only consists in a substitution in the Table of Indirection. Since the parent might not be in the Cache either, you should look for the grand-parent or grand-grand-parent, etc. I add the “base tile” (the lowest resolution one) to the set of visible tiles, so I’m always sure that at least the Base tile is in the Cache.

The Cache

The Cache itself is simply a texture (not mipmapped, but interpolated), which forms a grid of tiles. You need somewhere a table of correspondance between a position in the Cache and a TileId. I use a dictionary, indexed by the TileId. Use glTexSubImage2D() to replace only the part of the texture which contains the new tile.

When the Cache is full, some tiles must be dropped. People and I use a simple Least Recently Used mechanism to determine which ones. It’s simple, it works. I tried other heuristics, based on the mipmap level, to drop the least detailed tiles in last resort, but it did not work great, leading to load the most detailed tiles too frequently.

Dropping a Tile consists in marking its position as free in the table of correspondance. Since it does not perform OpenGL calls, it can be done at any time.

The Cache does not have to be huge: 16 x 16 tiles works. In my experience 8 x 8 tiles is not big enough on a Retina display: the program loads tiles and drops them continually. Make the Cache bigger if you want to remove some burden on the CPU, or have several Megatextures.

A 256 x 256 tile takes 250 KB of memory, so a 16 x 16 tiles cache takes 64 MB. That is very reasonable.

Table of Indirection

The Texturing Shader needs to know what are the coordinates of a Tile in the Cache texture. For that purpose, it is provided a Table of Indirection, which is a mip mapped texture.

A pixel of the texture contains the following information:

x position in the cache (stored in .r)
y position (in .g)
mipmap level (in .b).

For a particular mipmap level, this table has one pixel per tile. For instance, say that my megatexture measures 256 x 256 tiles at mipmap 0, then the texture measures 256 x 256 pixels at mipmap 0. There are only 128 x 128 tiles at mipmap 1, and hence the table measures 128 x 128 pixels at mipmap 1. There is a straight correspondance, so the Texturing Shader determines the tiles x and y according to the texture coordinate, and does a simple look up. (In other words, there is a table of indirection for every mipmap level. All these tables are combined in a single mipmapped texture).

You might wonder why the mipmap level is stored in the table, since it can be determined by the shader. Actually, this is what allows to substitute a parent tile; in that case, the pixel contains the x, y, and mipmap of the parent — not the child. The mipmap level of the parent tile is needed to determine correctly the position within the parent tile.

Texturing Shader

Finally, we arrive at the end of the chain !

const float phys_tex_dimension_pages = 32.0;
const float page_dimension = 256.0;
const float page_dimension_log2 = 8.0;
uniform float mip_bias;
uniform sampler2D pageTableTexture;
uniform sampler2D physicalTexture;
// converts the virtual texture coordinates to the physical texcoords
vec2 calulateVirtualTextureCoordinates(vec2 coord) {
        float bias = page_dimension_log2 - 0.5 + mip_bias;
        vec4 pageTableEntry = texture2D(pageTableTexture, coord, bias) * 255.0;
        float mipExp = exp2(pageTableEntry.a); // alpha channel has mipmap-level
        vec2 pageCoord = pageTableEntry.bg; // blue-green has x-y coordinates
        vec2 withinPageCoord = fract(coord * mipExp);
        return ((pageCoord + withinPageCoord) / phys_tex_dimension_pages);
    }
    void main(void) {
        vec2 coord = calulateVirtualTextureCoordinates(gl_TexCoord[0].xy);
        vec4 vtex = texture2D(physicalTexture, coord);
        gl_FragColor = vtex;
    }

Let’s first see the main() function: you’ll notice that texturing of the object is done as usual, using the boring texture2D() function.

phys_tex_dimension_page is the size of the megatexture at the maximum mipmap level, expressed in Tiles
page_dimension is the size of a Tile, it is not used directly because…
… page_dimension_log2 is the same size, but expressed as a power of two (2^8 = 256)
pageTableTexture is the Indirection Table texture
physicalTexture is the Cache texture

Now, I would like to focus on this:

float bias = page_dimension_log2 - 0.5 + mip_bias;
vec4 pageTableEntry = texture2D(pageTableTexture, coord, bias) * 255.0;

Older versions of OpenGL (like the one I’m constrained to use, because of Apple), did not allow to sample the texture for a particular mip map level. However, texture2D() may take a third parameter, which is a value added to the implicit mipmap level (the one computed by OpenGL). I don’t know why 0.5 is substracted, but it works better this way.

I must say that I had a lot of problems with this principle because it assumes that:

Tiles are square
The megatexture is square

Since I had to cover the Earth, my megatexture was not square, but had a 2:1 ratio instead. And my tiles were 512 x 256 pixels. If this is not the case, you will run into troubles, because the computation of the mipmap level is right vertically, but not horizontally, and you will have visual artifacts, since the wrong mipmap level is sampled from the indirection texture.So, don’t do that: make your tiles square and stretch your megatexture if needed. It will save you a lot of pain.

(With a more recent version of GLSL, you might use texture2DLod(), and compute the mipmap level like in the Determination Shader, and not have this problem).

The MegaTexture is precise enough so the Richat Structure (“The Eye of Sahara”) can be seen.

Generating Tiles

We’re not done yet! Remember that the megatexture is a huge image that must be cut into tiles. I personally wrote a Python program that uses Image Magick to cut tiles and resize them. I won’t go into details here, but you must know that Image Magick is slow, and not very user friendly (and that by default, rescaling is proportional). You may do it otherwise, maybe using a Photoshop script or whatever.

Seams

There is one final problem with the principle of the megatexture itself. Because tiles are all stored in the Cache texture in a random order, a tile is unrelated with its neighbours. This causes visual problems because of the linear interpolation of tiles, which will cause half a pixel of neighbour textures to show:

The solution is well known: leave a margin of 1 pixel on each side, and sample at this size. Hence the actual usable size is 254 x 254 pixels on my 256 x 256 tiles.