An Analysis of Minecraft-like Engines

Voxel engines are everywhere…

…and given the huge success of Minecraft, I think they are going to become a permanent fixture in the landscape of game engines.  Minecraft (and Infiniminer) live in a unique space somewhere between high-resolution voxels and static prefabricated environments.  Modifying a Minecraft world is both intuitive and elegant; you simply pick up and place blocks.  In many ways, it is the natural 3D generalization of older 2D tile based games.  Of course creating a Minecraft engine is much more challenging than making a 2D engine for a number of reasons.  First of all, it is in 3D and second it is generally expected that the environments must be much more dynamic, with interactive lighting, physics and so on.

When building a voxel game, it is important to choose a data structure for representing the world early on.  This decision more than any other has the greatest impact on the overall performance, flexibility, and scale of the game.  This post discusses some of the possible choices that can be made along these lines, and hopefully give a better explanation of what sort of technical considerations are involved.  In fact, I posit that some commonly held performance assumptions about Minecraft-type engines are probably wrong. In particular, the importance of random-access reads/writes for modifications is often overestimated, while the relative cost of iteration is greatly underestimated.  In concluding this article, I introduce a new data structure that exploits this observation, and ultimately gives much faster iteration (for typical Minecraft-esque maps) at the expense of slightly slower random accesses.

Conventional data structures

But before getting ahead of ourselves, let’s review some of the approaches were are currently being used to solve this problem.  We start with the absolute simplest choice, which is the “flat array”.  In the flat array model, you simply store all of your level data in one gigantic array:

As Minecraft and Infiniminer haver shown us this strategy can be viable — providing you limit the number of voxel types and make the world small enough.  Storing blocks in a big array has many advantages, like constant time random access reads and writes, but it is far from perfect. The most obvious downside is that the amount of memory consumed grows with the cube of the linear dimension, and that they are ultimately fixed in size.  These properties make them suitable for small, sand-box worlds (like the old-school Minecraft creative mode), but prevent them from being used in an infinite world (like Minecraft survival mode).

The most naive solution to the perceived disadvantages of flat arrays (and I say this because it is almost invariably the first thing that anyone suggests) is to try putting all the voxels in an octree.  An octree is a tree that recursively subdivides space into equal sized octants:

On the face of it, this does sound like a good idea.  Minecraft maps are mostly empty space (or at least piecewise constant regions), and so removing all those empty voxels ought to simplify things quite a bit.  Unfortunately, for those who have tried octrees, the results are typically not so impressive:  in many situations, they turn out to be orders of magnitude slower than flat arrays.  One commonly cited explanation is that access for neighborhood queries times in octrees are too slow, causing needless cache misses and that (re)allocating nodes leads to memory fragmentation.

However, to understand why octrees are slower for Minecraft games, it isn’t really necessary to invoke such exotic explanations.  The underlying reason for the worse performance is purely algorithmic: each time an arbitrary voxel is touched, either when iterating or performing a random access, it is necessary to traverse to the entire height of the octree.  Since octrees are not necessarily balanced, this can be as much as log(maximum size of game world)!  Assuming the coordinate system is say 32 bit, this brings an added overhead of 32 additional non-coherent memory accesses per each operation that touches the map (as opposed to the flat array, where everything is simply constant time).

A more successful method — and even less sophisticated — method is to use paging or “virtual memory” to expand the size of the flat array.  To do this one breaks down the flat 3D array into a collection of pages or “chunks”, and then uses a hash table to sparsely store only a subset of the pages which are currently required by the game:

Pages/chunks are typically assigned in the same way as any virtual memory system: by reserving the upper -bits of the coordinate for the chunk ID.  Using a virtual  memory system has the additional advantage that chunks can be lazily initialized, which is great for games with procedural content.  The same concept also allows chunks to be mapped back to disk when they are not needed (that is when there are no players nearby).  Using a hash map to store the chunks allows one to maintain constant time random access, while simultaneously taking advantage of sparsity as in an octree.

It seems like the current direction of Minecraft engines has begun to converge on the “virtualization” approach.  (And I base this claim on numerous blog posts).  But is this really the best solution?  After all, before we jump in and try to solve a problem we should at least try to quantitatively understand what is going on first; otherwise we might as well be doing literary criticism instead of computer science.  So I ask the following basic question:

What factors really determine the performance of a voxel engine?

There is a superificial answer to this question, which is “read/write access times”; after all every operation on a Minecraft map can be reduced to simple voxel queries.  If you accept this claim, then a virtualization is optimal — end of story.  However for reasons which I shall soon explain, this assertion is false.  In fact, I hope to persuade you that really it is the following that is true:

Claim:  The main bottleneck in voxel engines is iteration — not random access reads and writes!

The basic premise hinges on the fact that the vast majority of accesses to the voxel database come in the form of iterations over 3D ranges of voxels.  To demonstrate this point let us first consider some of the high level tasks involving the world map that a voxel engine needs to deal with:

• Placing and removing blocks
• Collision detection
• Picking/raytracing
• Lighting updates
• Physics updates
• Mesh generation
• Disk serialization
• Network serialization

The last two items here could be considered optional, but I shall leave them in for the sake of completeness (since most voxel games require some level of persistence, and many aspire at least to one day be multiplayer games).  We shall classify each task in terms of whatever operations are required to implement it within the game loop, which are either read/write and random access/iteration.  To measure cost, we count the number of raw voxel operations (read-writes) that occur over some arbitrary interval of gameplay.  The idea here is that we wait some fixed unit of time (t) and count up how many times we hit the voxels (ie the number of bytes we had to read/write).  Under the assumption that the world is much larger than cache (which is invariably true), these memory operations are effectively the limiting factor in any computation involving the game world.

At a fairly coarse level, we can estimate this complexity in terms of  quantities in the game state, like the number of MOBs and players in the world, or in terms of the number of atomic events that occur in the interval, like elapsed frames, chunk updates etc.  This gives an asymptotic estimate for the cost of each task in terms of intelligible units:

Task	Operation	Cost (in voxel ops)
Block placement	Write	number of clicks
Collision check	Read	(number of MOBs) * (frame count)
Picking	Read*	(pick distance) * (frame count)
Lighting	Iteration (RW)	(torch radius)^3 * (number of torches) + (active world size) * (sun position changes)
Physics update	Iteration (RW)	(active world size) * (frame count)
Mesh generation	Iteration (R)	(chunk size) * (chunk updates)
Disk serialization	Iteration (R)	(active world size) * (save events)
Network IO	Iteration (R)	(active world size) * (number of players) + (chunk size) * (chunk updates)

* Assuming that picking is implemented by ray-marching.

Disclaimer:  The entries in this chart are estimates, and I do not claim that the dimensions for each task are 100% accurate.  If you disagree with anything here, please leave a comment.

Our goal is to pinpoint which operation(s) in this chart is the bottleneck, and to do so we need dimensional analysis to convert all our units into a common set of dimensions.  Somewhat arbitrarily, we choose frame count and number players, and so our units will be in (voxel ops * frames * players).  Going through the list of each item, we now try to estimate some plausible bounds for each dimension:

1. (number of clicks) is the amount of times the user clicked the mouse in the interval.  Since a user can click at most once per frame (and in practice is likely to click much much less), this value is bounded above by:
   • (number of clicks) < (frame count) * (number of players)
2. (number of MOBs)is the number of active entities, or movable objects in the game.  Assuming that they are not moving too fast (ie there is a maximum speed of a constant number of blocks per frame), then the number of reads due to a collision check is going to be proportional to the number of MOBs and the number of frames.  We assume that number of MOBs in the world is proportional to the number of players, according to some constant ,
   • (number of MOBs) = * (number of players)
3. (pick distance) is the maximum range a block can be selected by a player, and we shall assume that it is just a small constant.
4. Again, (torch radius) is a small number.  In Minecraft, this could theoretically be as high as 16, but we will simply leave it as a tiny constant for now.
5. (number of torches)is the number of torches placed / removed during the interval.  Torch placement is generally an irregular event, and so we estimate that a player will place a torch somewhat infrequently, at a rate goverened by some constant ,
   • (number of torches) = (number of players) * (frame count) *
6. (sun position changes)is the number of times that the sun changes.  This happens regularly every few frames, and so we estimate that,
   • (sun position changes) = (frame count) *
7. (active world size)is the size of the world.  It is basically proportional to the number of chunks visible by each player.  Assuming a Minecraft style world with a height cap, this would vary quadratically with the visible radius, , or cubically for a world with unlimited z-value.
   • (active world size) = (number of players) * (chunk size) *
8. (chunk updates) occur whenever a chunk changes.  Assuming this is random, but proportional to the size of the active world, we get
   • (chunk updates) = (active chunk size)  / (chunk size) *
9. (save events) occur at regular intervals and save the map to disk.  Again, because they are regular we can just bound them by a constant times the number of frames.
   • (save events) = (frame count) *

Plugging all that in, we get the following we get the following expression for the total amount of random access voxel ops:

(random voxel ops) = C * (number of players) * (frame_count) 

And for the sequential or iteration voxel ops, we get:

(sequential voxel ops) = C * (number of players) * (frame count) * (chunk size) * 

In general, the last quantity is much larger, by a factor of (chunk size) * , or (chunk size) * if you have unlimited height.  So, if you believe these estimates are reasonable, then you should be convinced that iteration by far dominates the performance of a Minecraft style game.  In fact, this estimate explains a few other things, like why visible radius is such an important performance factor in Minecraft.  Increasing the draw radius linearly, degrades the performance quadratically!

Can we make iteration faster?

The main conclusion to draw from the above estimate is that we should really focus our optimization efforts on iteration.  If we can bring that cost down measurably, then we can expect to see large improvements in the game’s performance.  But the question then becomes: is this even possible?  To show that at least in principle it is, we make use of a simple observation:  In any situation where we are iterating, we only need to consider local neighborhoods around each cell.  For example, if you are computing a physics update (in minecraft) you only consider the cells within your Moore neighborhood.  Moreover, these updates are deterministic:  If you have two voxels with the same Moore neighborhood, then the result of applying a physics update to each cell will be the same.  As a result, if we can iterate over the cells in groups which all have the same neighborhood, then we can apply the same update to all cells in that group at once!

This is the basic idea behind our strategy for speeding up iteration.  We will now try to explain this principle in more detail.  Let be the set of block types, and let be a map representing our world (ie an assignment of block types to integral coordinates on a regular grid).  A radius -Moore neighborhood is an element of the set .  The Moore neighborhood of a coordinate in is determined by a map,  such that,

Then an update rule is a map sending a Moore neighborhood to a block, , and the application of to is the rule sending,

Then the fact that we can group cells corresponds to the observation that:

So, if our world data allows us to group up similar Moore neighborhoods, then we can expect that on average we should be able to reduce the complexity of iteration by a substantial amount.

In fact, this same idea also applies to meshing:  If you have a sequence of adjacent cubes that all have the same 3x3x3 neighborhood, then you can fuse their faces together to make a smaller mesh.  For example, instead of generating a mesh for each cube one at a time:

We can fuse cells with the same Moore neighborhood together:

So if we can iterate over the cells in batches, we should expect that not only will our physics loop get faster, but we should also expect that mesh quality should improve as well!

Run length encoding and interval trees

At this point, we now have some pretty good intuition that we can do iteration faster, so the obvious next step is to figure out exactly how this works.  Here is one strategy that I have used in a simple javascript based Minecraft game (which you can get here).  The idea follows from the simple observation that Minecraft style worlds can be easily run-length-encoded.  In run length encoding, one first flattens the 3D array describing a chunk down to a 1D string.  Then, substrings of repeated characters or “runs” are grouped together.  For example, the string:

     aaaaabbbbacaa

Would become:

     5a4b1a1c2a

Which is a size reduction of about 25%.  You can imagine replacing the characters in the string with voxel types, and then compressing them down this way.  Now, applying run-length-encoding to chunks is by no means a new idea.  In fact, it is even reasonable to do so as a preprocessing step before g-zipping the chunks and sending them over the network/spooling to disk.  The fact that run-length encoding does compress chunks works to our advantage:  We can use this fact to iterate over the voxel set in units of runs instead of units of pixels, and can easily be extended to walk runs of Moore neighborhoods as well.

This sounds great for iterations, but what about doing random-access reads and writes?  After all, we still need to handle collision detection, block placement and raytracing somehow.  If we were dumb about it, just accessing a random block from a run-length-encoded chunk could take up to linear time.  It is obvious that we shouldn’t be happy with this performance, and indeed there is a simple solution.  The key observation is that a run-length encoding of a string is formally equivalent to an interval tree representation.

An interval tree is just an ordinary binary tree, where the key of each node is the start of a run and the value is the coordinate of the run.  Finding a greatest lower bound for a coordinate is equivalent to finding the interval containing a node.  To do insertion is a bit trickier, but not by much.  It does involve working through a few special cases by hand, but it is nothing that cannot be done without taking a bit of time and a pad of paper.  If we implement this data structure using some type of self-balancing binary tree (for example a red-black tree or a splay tree), then we can perform random reads and writes in logarithmic time on the number of runs.  Best of all: we also get improved mesh generation and in-memory compression for free!

Now in practice, you would want to combine this idea with virtualization.  Basically, you would store each chunk as an interval tree, then represent the entire world as a hash map.  The reason for doing this would be to integrate paging and chunk generation more seamlessly.  In fact, this is the method I took in the last two minecraft type games that I wrote, and you can find some example code illustrating this data structure right here.

Comparisons and conclusion

Ok, now that we’ve gone through the details, let’s break down our options by data structure and time complexity:

Data structure	Random access time	Iteration time	Size	Notes
Flat array
Virtualized array				is the number of occupied (virtual) chunks
Octree				is the height of the octree
Interval tree + Hash table				is the number of runs in a region, is the number of runs a single chunk.

Under our previous estimates, we can suppose quite reasonably that an interval tree should outperform a virtualized array for typical maps.  (Of course this claim goes out the window if your level geometry is really chaotic and not amenable to run-length encoding).  To test this hypothesis, I coded up a pretty simple benchmark.  At the outset the map is initialized to a 256x256x256 array of volume, with 5 different layers of blocks.  Randomly, 2^16 blocks from 32 different types are sprinkled around the domain.  To assess performance, we measure both how long an average random read takes and how long a sequential iteration requires (for Moore radii of 0 and 1).  The results of the benchmark are as follows:

Data structure	Avg. random read	Avg. sequential read (Moore radius = 0)	Avg. sequential read (Moore radius = 1)
Flat array	0.224 μs/read	0.178 μs/read	0.060 μs/read
Virtual array	0.278 μs/read	0.210 μs/read	0.107 μs/read
Octree	2.05 μs/read	0.981 μs/read	0.933 μs/read
Interval tree + hashing	0.571 μs/read	0.003 μs/read	0.006 μs/read

In each column, the best result is underlined in bold.  There are a few things to take home from this.  First, for random access reads and writes, nothing beats a flat array.  Second, octrees are not a very good idea for voxel worlds.  On all accounts, a virtual array is far superior.  Finally, while the access time for random reads is noticeably slower in an interval tree, the speed benefits of improved sequential iteration make them more than worthwhile in practice.  The code for the benchmark can be found here.

Of course, interval trees will not perform well in worlds which are more-or-less random, or if there are not many runs of similar blocks that can be grouped together.  But when your world can be described in this way (which is certainly true of a Minecraft like game), switching to intervals can offer a huge improvement in performance.  I also firmly believe that there is still quite a bit of room for improvement.  If you have any suggestions or ideas, please leave a comment!

Goofing around in node.js

In my spare time I’ve been messing around a bit with node.js, and have been trying out building a browser based multiplayer online game.  Right now I am using sprites/minecraft style levels since they are easy to make and low bandwidth; but ultimately I would like to try out doing some more stuff with meshes.  You can give it a try if you have a browser which supports WebGL.   Here is a screen shot:

The URL for the game is http://mmotest.nodester.com .  Right now I am using the freebie hosting coupon for nodester, but I have no idea how much longer it will last.  The database hosting is from mongolabs, and so again is fairly limited.  However, this does show that you can build something like a minecraft game in the browser pretty easily.  Of course the latency with a shared hosting provider like nodester starts to become a bigger problem as you get more players in the game.  I will try to post the source code once I get a chance to clean some things up.  I am having some issues configuring the nodester deployment so that I don’t have to hardcode in various strings like the database password.  Unfortunately, for the life of me I can’t figure out how to pass those sort of config things in on the command line via nodesters REST API, and so for now they are just stored in version control.  Once I can figure out how to remove this information, I will make the github repository public so that others can try it out.  On the other hand, if I can’t solve the problem I may try just do something irresponsible like release a repository with all the passwords in it and let everyone go crazy hacking my demo accounts.

Ludum Dare 21 Results

Well the results for Ludum Dare 21 are in!  Here is the final rank for my entry:

Ratings

#4	Innovation(Jam)	3.90
#29	Overall(Jam)	2.95
#30	Humor(Jam)	2.33
#31	Graphics(Jam)	3.30
#32	Fun(Jam)	2.55
#45	Audio(Jam)	2.12
#60	Coolness	8%
#63	Theme(Jam)	2.40
#104	Community	3.13

I have to say that I am a bit disappointed about only getting 4th place for innovation .  (At least I got in the top 10 in that category, which is somewhat respectable given the huge number of entries!)  I think that the controls definitely could have been done a bit better.  In retrospect it was a mistake to use a separate parallel transport frame for moving the camera and for applying forces.  This caused the camera to drift slightly which made the controls more confusing than they needed to be (plus it would have been a really easy fix had I thought of it during the competition).  I’m not sure if it is worth developing this concept into a full game.  Judging by the rating, it didn’t seem like it was that popular, but perhaps I am reading too much into the scores.  Also I think that the screen shot I picked for the entry was not nearly exciting enough (the wood level would have been much better), and the some of the levels definitely needed more polish.

At least there is a mini-LD coming up next weekend, so hopefully that will go a bit better!

Polynomials in C++

Last time, we showed how one can use symmetric tensors to conveniently represent homogeneous polynomials and Taylor series.  Today, I am going to talk about how to actually implement a generic homogeneous polynomial/symmetric tensor class in C++.  The goal of this implementation (for the moment) is not efficiency, but rather generality and correctness.  If there is enough interest we can investigate some optimizations perhaps in future posts.

How to deal with polynomials?

There are quite a few ways to represent polynomials.  For the moment, we aren’t going to assume any type of specialized sparse structure and will instead simply suppose that our polynomials are drawn at random from some uniform distribution on their coefficients.  In some situations, for example in the solution of various differential equations, this is not too far from the truth, and while there may be situations where this assumption is too pessimistic, there are at least equally many situations where it is a good approximation of reality.

In 1D, the most direct method to represent a dense, inhomogeneous polynomial is to just store a big array of coefficients, where the entry corresponds to the coefficient for the degree i monomial.  Naively generalizing a bit, the same idea works for multivariate polynomials as well, where you simply make the array multidimensional.  For example, suppose you had the 2D coefficient matrix:

Which we could then expand out to give a polynomial like:

Though it is a bit simplistic, working in this format does have a few advantages.  For example:

1. Evaluating a polynomial is easy.  Just compute the vectors and the then plug them into the coefficient tensor using ordinary Einstein summation.
2. Multiplying can also be pretty fast.  Observe that coefficient matrix of the product of two polynomial matrices is just the convolution of their two matrices.  Using fast convolution algorithms (for example an FFT based method would work), this can be computed with only a mere log factor overhead!

The amount of storage for representing a polynomial on variables with a maximum per-coefficient degree of using this scheme requires exactly scalars worth of memory. However, this representation scheme does have some major problems:

1. It doesn’t capture all the terms above degree r properly.  For example, if you take the term from the above example, it is actually second order, and so is the polynomial that it represents.  However, we are missing the extra terms for .
2. Related to the previous point, it is not closed under linear coordinate transformations.  If you were to say do a rotation in , that term could change into something like , and thus we would need to resize the matrix to keep up with the extra stuff.
3. If you only care about storing terms with degree at most , then the above representation is wasteful by at least a factor of about .
4. And if you want to deal with terms that are exactly degree d, then it is even worse with an overhead of .

The underlying cause for each of these problems is similar, and ultimately stems from the fact that the function is not multilinear.  As we saw last time it can be useful to deal with multivariate polynomials using tensors.  Similarly, it is also useful to split them up into homogeneous parts (that is grouping all the terms which have the same degree together).  We also showed how this approach greatly simplified the expression of a general multidimensional Taylor series expansion.  So today, we are going to investigate a different, and more “tensor-flavored” method for representing homogeneous polynomials.

Quick Review

Recall that a homogeneous polynomial on d variables, , with degree r is a polynomial,

We also showed that this could be written equivalently using tensor notation:

Which allows us to reinterpret p as a “r-multilinear function” on x.  We say that a function is multilinear if for any of its arguments , it is a linear function, that is:

While it may seem a bit restrictive to work in homogeneous polynomials only, in practice it doesn’t make much difference.  Here is a simple trick that ought to be familiar to anyone in computer graphics who knows a bit about projective matrices.  Let us consider the following 1D example (which trivially generalizes):  given a polynomial in one variable,

Let us create a new variable , and construct the new polynomial:

Now we observe that:

1. is a homogeneous polynomial of degree 2.
2. If we pick , then .

And so we have faithfully converted our inhomogeneous degree 2 polynomial in one variable into a new homogeneous degree 2 polynomial in two variables!

In general, converting from an inhomogeneous polynomial to a homogeneous polynomial is a simple matter of adding an extra variable.  The fully generalized version works like this:  Suppose that we have a polynomial of degree r in d unknowns.  Then we add a variable, w, and construct the new homogeneous degree r polynomial,

And again, .

Indexing, iteration and combinatorics

So with that stuff above all out of the way, we shall now settle on an implementation: we’re going to use symmetric tensors to represent homogeneous polynomials.  Let us start off by listing a few basic properties about symmetric tensors:

1. A tensor is called symmetric if for all permutations of an index, and so on, .
2. As a result, all vectors (aka a rank-1 tensors) are trivially symmetric.
3. Similarly, a matrix is symmetric if and only if .
4. A rank 3 tensor is symmetric if (and so on).
5. Consequently is the dimension of each of the indices in a symmetric tensor must be equal, and so we shall sometimes refer to this quantity as the “dimension of the tensor”, or when we need to be more specific we shall call it the index dimension.
6. The dimension of the space of all symmetric tensors with rank r and (index) dimension d is .

And so the first problem that we must solve is to figure out how to pack the coordinates of a general rank r, dimension d symmetric tensor into an array of size .  Closely related to this issue is the problem of indexing, and we are going to need methods for both randomly accessing elements of the tensor and for performing sequential iteration over the unique entries.  In fact, it is really this last issue which is the most essential, since if we can define an efficient method for enumerating the unique elements of the tensor, then we can simply use the order in which we enumerate the entries of the tensor as the indexes into a flat array.

By definition, the entries of a symmetric tensor are invariant under index permutation.  This leads to two distinct ways to index into a symmetric rank r dimension d tensor:

1. Tensor style indexing – This is the ordinary way that we write the index for a generic tensor using subscripts.  Entries of the tensor are parameterized by a length r vector of integers, where for all .
2. Polynomial style (degree) indexing – This is the style of indexing is unique to symmetric tensors and comes from their close association with polynomials.  Here, we take the index to be the degree of one of the monomial terms in the corresponding homogeneous polynomial   These indices are enumerated by length d vectors of integers, where for all k, and .

We can convert between these two types of indexes up to a point.  Converting from tensor indexing to degree indexing is straightforward, where we just pick

the number of elements in that are equal to k

Going the other way is a bit harder, since there are multiple tensor indices which map to the same degree index.  To resolve the ambiguity, we can simply take the lexicographically first item.  This same idea also turns out to be useful for enumerating the entries within our tensor.  For example, suppose that we wanted to enumerate the components of a rank 3, degree 4 tensor.  Then, we could use the following ordering,

Pos.   Tensor   Degree
-------------------------
0      (0,0,0)  (3,0,0,0)
1      (1,0,0)  (2,1,0,0)
2      (2,0,0)  (2,0,1,0)
3      (3,0,0)  (2,0,0,1)
4      (1,1,0)  (1,2,0,0)
5      (2,1,0)  (1,1,1,0)
6      (3,1,0)  (1,1,0,1)
7      (2,2,0)  (1,0,2,0)
8      (3,2,0)  (1,0,1,1)
9      (3,3,0)  (1,0,0,2)
10     (1,1,1)  (0,3,0,0)
11     (2,1,1)  (0,2,1,0)
12     (3,1,1)  (0,2,0,1)
13     (2,2,1)  (0,1,2,0)
14     (3,2,1)  (0,1,1,1)
15     (3,3,1)  (0,1,0,2)
16     (2,2,2)  (0,0,3,0)
17     (3,2,2)  (0,0,2,1)
18     (3,3,2)  (0,0,1,2)
19     (3,3,3)  (0,0,0,3)

Here, the left hand column represents the order in which the entries of the tensor were visited, while the middle column is the set of unique tensor indices listed lexicographically.  The right hand side is the corresponding degree index, which is remarkably in colexicographic order!  This suggests other interesting patterns as well, such as the following little theorem:

Proposition: Colexicographically sort all the degree indices in a rank r, d-dimensional symmetric tensor.  Then the position of the degree index, , is

(That is supposed to be a binomial coefficient, but the parenthesis seem to be a bit messed up in wordpress.)  I will leave this proof as an exercise to the reader, but offer the hint that it can be proven directly by induction.  It is also interesting to note that the formula for the position does not depend on either the rank of the tensor or the first coefficient.

Putting it all together

Now let us try to work out a method for iterating the tensor index incrementally.  A simple strategy is to use a greedy method:  We just scan over the index until we find a component that can be incremented, then scan backwards and reset all the previous indices.  In psuedo-C++, the code looks something like this:

bool next() {
  for(int i=0; i<rank; ++i) {
    if(tensor_index[i] < dimension - 1) {
      ++tensor_index[i];
      for(int j=i-1; j>=0; --j) {
          tensor_index[j] = tensor_index[i];
      }
      return true;
    }
  }
  return false;
}

The tensor_index array contains rank elements, and is initialized to all 0′s before this code is called.  Each call to next() advances tensor_index one position.

Superficially, there are some similarities between this algorithm and the method for incrementing a binary counter.  As a result, I propose the following conjecture (which I have not investigated carefully enough to say is true confidently):

Conjecture: Amortized over course of a full traversal, the cost  of next is constant.

Using all this stuff, I’ve put together a little C++ class which iterates over the coordinates of a symmetric tensor.  I’ve also uploaded it to github, where you can go download and play around with it (though it isn’t much).  I hope to keep up this trend with future updates.  For those who are curious, here is the URL:

https://github.com/mikolalysenko/0fpsBlog

Next time

In the next episode, I am going to continue on this journey and build an actual symmetric tensor class with a few useful operations in it.  More speculatively, we shall then move on to rational mappings, algebraic surfaces, transfinite interpolation and deformable objects!

Multidimensional Taylor Series and Symmetric Tensors

As a warm up post, I’m going to talk about an important generalization of something that should be familiar to anyone’s who has made through a semester of calculus: Taylor series!  (And if you haven’t seen these guys before, or are perhaps feeling a bit rusty, then by all means please head on over to Khan academy to quickly brush up.  Go right ahead, I’ll still be here when you get back.)

Now this probably sounds like an awfully scary title for the first article in this miniseries about mathematics in graphics programming.  Perhaps you’re right!  But I’d like to think that graphics programmers are a tough bunch, and that using this kind of a name may indeed have quite the opposite effect of emboldening those cocky individuals who think they can press on and figure things out in spite of any frightening mathematical jargon.

Taylor Series

At the very heart of this discussion we are going to deal with two of the most important tasks any graphics programmer needs to worry about:  approximation and book keeping.  Taylor series are of course one of the oldest and best known methods for approximating functions.  You can think of Taylor series in a couple of ways.  One possibility is to imagine that they are successively approximating a given input function by adding additional polynomial terms to it.  So the idea is that you can start with some arbitrary 1D function, let’s call it , and suppose that you are allowed the following two operations:

1. You can evaluate a function at 0.
2. You can take a derivative,

Then, we can compute the Taylor series expansion of f about 0 in the usual way,

and so on

If we’re really slick, we can save the first coefficients for these polynomials in a vector, call them say ., and then we can evaluate some approximation of f by summing up the first k terms in the above approximation:

Here is an animated gif showing the convergence of the Taylor series for the exponential function that I shamelessly ripped off from wikipedia:

Taylor series approximation for the function about the point , copyright the Wikimedia foundation. Credit Oleg Alexandrov

Higher Dimensional Taylor Series

It is easy to adapt Taylor series to deal with vector valued functions in a single variable, you just treat each component separately.  But what if the domain f was not one dimensional?  This could easily happen for example, if you were trying to approximate something like a surface, a height field, or an image filter locally.  Well, in this post I will show you a slick way to deal with Taylor series of all sizes, shapes and dimensions!  But before I do that, let me show you what the ugly and naive approach to this problem looks like:

Suppose that is a 2D scalar-valued function. Then, let us look for a best quadratic (aka degree 2) approximation to f within the region near .  By analogy to the 1D case, we want to find some 2nd order polynomial $p(x,y)$ in two variables such that:

And:

Phew!  That’s a lot more equations and coefficients than we got in the 1D case for a second order approximation.  Let’s work through solving for the coefficient $a_{20}$:  To do this, we take p and plug it into the appropriate equation:

If we generalize this idea a bit, then we can see that for an arbitrary Taylor series approximation in 2D, the coefficient is determined by the equation,

All this is well and good, but it has a few problems. First, the summation for p is quite disorganized.  How are we supposed to keep track of which coefficients go where?  If we want to store the coefficients of p on a computer, then we need some less ad-hoc method for naming them and packing them into memory.  Second, it seems like this is going to be a headache to deal with 3 or more dimensions, since we’ll need to basically repeat the same sort of reasoning over and over again.  As a result, if we want to start playing around with higher dimensional Taylor series, we’re going to need a more principled notation for dealing with higher order multivariate polynomials.

Tensors to the Rescue!

One such system is tensor notation!  Though often maligned by mathematical purists. algebraic geometers, and those more modern algebraically savvy category theorists, tensors are a simple and effective notation that dramatically simplifies calculations.  From a very simplistic point of view, a tensor is just a multidimensional array.  The rank of the tensor is the number of different indices, each of which has a distinct dimension.  In C-style syntax, you could declare a tensor in the following way:

  float vector[10];                          //A rank 1 tensor, with dimension (10)
  float matrix[5][5];                        //A rank 2 tensor, with dimensions (5, 5)
  float spreadsheet[ROWS][COLUMNS][PAGES];   //A rank 3 tensor, with dimensions (ROWS, COLUMNS, PAGES)
  float crazy_thing[10][16][3][8];           //A rank 4 tensor, with dimension (10, 16, 3, 8 )

We will usually write tensors as upper case letters.  To reference an individual entry in this array, we will use an ordered subscript, like so:

Which we can take to be equivalent to the C-code:

  M[i][j]

For the rest of the article, we are going to stick with this point of view that tensors are just big arrays of numbers.  We aren’t going to bother worrying about things like covariance/contravariance (and if you don’t know what those words are, forget I said anything), nor are we going to mess around too much with operators like tensor products.  There is nothing wrong with doing this, though it can be a bit narrow minded and it does somewhat limit the applications to which tensors may be applied.  If it bothers you to think about tensors in this way, here is a more algebraic/operational picture of what a tensor does: much as how a row vector can represent a scalar-valued linear function, (via the dot-product), a tensor can represents a multi-linear function, : that is, it takes in several vectors and spits out a scalar which varies linearly in each of its arguments.

That said, even if we just think about tensors as arrays, there’s still a number of useful, purely formal, operations that one can perform on tensors.  For example, you can add them up just like vectors.  If A and B are two tensors with the same dimensions, then we can define their sum componentwise as follows:

Where we take the symbol i to be a generic index here.  The other important operation that we will need to use is called tensor contraction.  You can think of tensor contraction as being something like a generalization of both the dot product for vectors, and matrix multiplication .  Here is the idea:  Suppose that you have two tensors with dimensions   and ranks n,m respectively and some index between them with a common dimension.  Then we can form a new tensor with rank n + m – 1 by summing over that index in A and B simultaneously:

Writing all this out is pretty awkward, so mathematicians use a slick short hand called Einstein summation convention to save space.  Here is the idea:  Any time that you see a repeated index in a product of tensor coefficients that are written next to each other, you sum over that index.  For example, you can use Einstein notation to write the dot product of two vectors, , as follows:

Similarly, suppose that you have a matrix M, then you can write the linear transformation of a vector x by M in the same shorthand,

Which beats having to remember the rows-by-columns rule for multiplying vectors!  Similarly, you can multiply two matrices using the same type of notation,

It even works for computing things like the trace of a matrix:

We can also use tensor notation to deal with things like computing a gradient.  For example, we write the derivative of a function with respect to the component as .  Combined with Einstein’s notation, we can also perform operations such as taking a gradient of a function along a certain direction.  If is some direction vector, then the derivative of along evaluated at the point is,

Symmetric Tensors and Homogeneous Polynomials

Going back to multidimensional Taylor series, how can we use all this notation to help us deal with polynomials?  Well, let us define a vector whose components are just the usual coordinate variables, and pick some rank 1 tensor A with appropriate dimension.  If we just plug in x, then we get the following expression:

Which is of course just a linear function on x! What if we wanted to make a quadratic function?  Again, using tensor contraction this is no big deal:

Neat!  This gives us a quadratic multivariate polynomial on x.  Moreover, it is also homogeneous, which means that it doesn’t have any degree 1 or lower terms sitting around. In fact, by generalizing from this pattern if we wanted to say store an arbitrary degree n polynomial, we could pack all of its coefficients into a rank n tensor, and evaluate by contracting:

But there is some redundancy here.  Notice in the degree 2 case, we are assuming that the components of x commute, or in other words the terms are really the same and so we really don’t need to store both the coefficients and .  We could make our lives a bit easier if we just assumed that they were equal.  In fact, it would be really nice if whenever we took any index, like say and then permuted it to something arbitrary, for example , it always gave us back the same thing.  To see why this is, let us try to work out what the coefficient for  the monomial in the degree n polynomial given by .  Directly expanding using Einstein summation, we get a sum over all permutations of the indices :

If we assume that all those coefficients are identical, then that above nasty summation has a pretty simple form:  namely it is a multinomial coefficient!  As a result, if we wanted to say find the coefficient of in $A_{i j k} x_i x_j x_k$, then we could just use the following simple formula:

A tensor which has the property that its coefficients are invariant under permutation of its indices is called a symmetric tensor, and as we have just seen symmetric tensors provide a particularly efficient method for representing homogeneous polynomials.  But wait!  There’s more…

If we pick , then the above formula is almost exactly the right formula for the Taylor series coefficient of in the expansion of .  The only thing is that we are off by a factor of , but no matter, we can just divide that off.  Taking this into account, we get the following striking formula for the Taylor series expansion of a function about the origin,

Which is remarkably close to the formula for 1D Taylor series expansions!

Next Time

I will show how to actually implement symmetric tensors in C++11, and give some example applications of multidimensional Taylor series in implicit function modeling and non-linear deformations!  I’ll also show how the above expansion can be simplified even further by working in projective coordinates.

A call for more technical blogs

There is a trend these days to avoid writing about technical things in programming — and in particular game development — blogs.  Just go to places like r/programming or altdevblogaday, and so on and you find plenty of articles giving you great advice on everything EXCEPT math and programming!  What gives?

There’s just not enough articles any more that get down to the nuts and bolts of algorithms and data structures, or at an even more basic level actually bother to try explaining some interesting theoretical concept.  Once venerable clearing houses like flipcode or gamedev.net have either shutdown or become diluted into shallow social-networking-zombies.  Perhaps this is a symptom of our ever decreasing attention spans, or even more pessimistically a sign that indie devs have simply given up on trying to push the technological envelope out of fear of competing with big studios.  It seems that trying to innovate technically is now viewed as `engine-development’ and derided as unproductive low-level tinkering.  Wherever it comes from, I am convinced that these insubstantial discussions are making us all stupider, and that the we need to start writing and paying attention to hard technical articles.

So, rather than sit back and complain, I’ve decided to do something proactive and try to update this blog more often with useful — or at least interesting and substantial — technical content.  But before that, I will start by listing some of the things I am *not* going to talk about:

• Industry news
• Business advice
• Coding “best practices”
• Social networking
• Marketing
• Project management

As I’ve just described, there’s already an abundance of literature on these subjects, possibly because they are trivial to write about, and it is easy to have an opinion about them. As a result, I don’t really feel particularly compelled, or even much qualified as an industry-outsider-academic-type, to discuss any of those things. More pragmatically, I also find that discussing these sorts of “soft”, subjective issues leads to either pointless back patting or unproductive bickering.

Instead, I’m going to use this blog to talk about the “harder” foundational stuff.  When possible, I will try to provide real code here — though I will probably switch languages a lot — but my real focus is going to be on explaining “why” things work the way they do. As a result, there will be math, and this is unavoidable.  I’m also going to make an effort to provide relevant references and links when possible, or at least to the extent of my knowledge of the prior art.  However, if I do miss a citation, please feel free to chime in and add something.

I don’t have a set schedule for topics that I want to cover, but I do have some general ideas.  Here is a short and incomplete, list of things that I would like to talk about:

• Procedural generation and implicit functions
• Physical modeling using Lagrangians
• Mesh data structures
• Spatial indexing
• Non-linear deformable objects
• Collision detection and Minkowski operations
• Fourier transforms, spherical harmonics and representation theory
• Applications of group theory in computer graphics
• Audio processing

I’m also open to requests at this stage.  If there is a topic that more people are interested in, I can spend more time focusing on that, so please leave a request in the comments.

New paper out!

Haven’t updated this blog in a long time, but I figure that since it is summer now maybe I will finally be able to keep this thing regularly updated (see how long that lasts…)  Anyway, some good news is that my latest paper has been accepted to the SIAM/SIGGRAPH conference on solid and physical modeling.  Here is a link to the download:

http://sal-cnc.me.wisc.edu/index.php?option=com_remository&Itemid=143&func=fileinfo&id=183

There’s a bunch of neat ideas in here, but the big idea here is the Minkowski product.  The motivation for this comes from the basic Minkowski sum.  If we recall, for two sets , their Minkowski sum is defined as:

Or alternatively,

The idea behind the Minkowski product is to replace with an arbitrary group, .  If we do this, then we can define the Minkowski product over in the following way:

Much like how the Minkowski sum is useful for collision detection of translating objects, the generalized Minkowski product can be used for collision detection between translating and rotating bodies.  Similarly, we can also define a Minkowski quotient which is analogous to the Minkowski difference:

And even better yet, we show how to compute these quantities using convolution algebras!  These operators turn out to be very useful in understanding things like partial symmetries of solids, mechanism workspaces and robotics.  All of the gory details and more are in the paper!

BSP Trees and Collision Detection

Recently I published a paper in CAD on BSP tree merging via linear programming.  At a very high level, the idea behind this paper was to replace the tree partitioning step with a linear programming feasibility test.  Doing this leads to a bunch of nice things, such as the removal of the partition function and all 9 of its special cases.  It also has the benefit that that for some suitable linear programming algorithms, the cost of evaluating feasibility can be amortized over the tree traversal giving an optimal linear-time output sensitive algorithm.

Of course as it always turns out with these things, I recently discovered that I was not the first to think of using linear programming for solving this problem.  Just last week Prof. Alberto Paoluzzi visited our lab, and I was lucky enough to get a chance to talk to him about BSP trees.  As it just so happens, he proposed using linear programming for tree partitioning in his book Geometric Programming for Computer-Aided Design (which is actually a very good read.)  Of course, this shouldn’t diminish the contribution of our paper, as we do present a much simpler merge algorithm along with an analysis of both time complexity and robustness.  In any future works, I will definitely add a reference to his book.

Anyway, I wanted to talk about one idea which didn’t make it into that paper, but I still happen to think is rather interesting.  As an accidental side effect of our merging algorithm, we discovered a straightforward method for solving exact continuous collision detection between arbitrary moving non-convex objects.  The key idea to this process was initially suggested to me by my colleague Nicholas Smolinske (who sadly is no longer in graduate school.)  Here’s how it goes:

Suppose you have two convex shapes moving at a constant velocity.  For each shape, the object Minkowski summed with its trajectory traces out a 4D convex set space and time.  If you take these two sets for both shapes and intersect them, then the result is also a convex shape in 4D.  Moreover, the lowest point in time within the overlap region must be the point at which the two objects initially made contact.  This idea is illustrated in the following diagram:

Because all of the sets we are dealing with are convex, it is possible to exactly solve for the point of intersection using linear programming with an objective function .  Now, the next question is how to do this with BSP trees?  Well, it turns out that one side effect from the algorithm proposed in our paper is that it ends up evaluating linear programming on each convex cell contained in the final BSP tree.  We can make use of this by simply picking the same objective function as above, and storing the minimal point over all convex cells as our result.  Alternatively, if we only want to get the point of intersection, then we don’t even need to compute the result of the merge so we may instead use it as a kind of driver for solving linear programming over the cell-complex hierarchically.

The idea of BSP trees representing shapes in space/time is not new, and to the best of my knowledge was first explored by Agarwal et al. (they give some novel tree construction strategies, but don’t really deal with merging or related issues).  Another interesting way to look at this idea is that one could understand the BSP tree as generalizing some of the ideas of polytope hierarchies (such as Dobkin-Kirkpatrick) to non-convex sets.  Of course this still doesn’t answer the big question which is how to go about building decent BSP trees (or hierarchies for that matter) and also how to pick which order to merge them.  Someday I might go back and try refining these ideas to the point where they could make for a nice short conference paper, but for now it is just something fun to think about.

Hello World

Hi.  My name is Mikola Lysenko and this is my first post on 0FPS. I am currently a first year graduate student at the University of Wisconsin-Madison, working in the spatial automation lab.  I used to run the website assertfalse.com, but I ended up taking it down to save money on bandwidth expenses.  Anyway, I figure that lately I’ve been too lazy in my writing process and that a blog is a pretty good way to bring discipline back into my life.  I intend to focus on research topics in solid modeling and computer graphics, as well as other ideas I find interesting.