Perlin noise procedural derivative
The other day I was skimming through the forums and I stumbled across this post concerning reading material about noise functions. That subject being near and dear to my heart, and one of my projects being a library to modularize the generation of noise, I took a look at some of the posted links. A couple of them seemed interesting, notably:
The second of the above links details a couple noise variants that the author came up with, that he calls Swiss noise and Jordan noise (named for their apparent similarity to terrain found in Switzerland and Jordan, respectively). The variants utilize a technique that I first saw used by Inigo Quilez; that is, using the derivative of the noise of one layer/octave to modify the contributions of successive layers. It's a neat idea, but one that, unfortunately, isn't well-supported by the current versions of the noise generators in my Accidental Noise Library. The existing support for derivative noise uses a specific module that approximates the derivative using multi-sampling, and requires one module (with the attendant overhead of multi-sampling the source) for each axis of the noise function, so constructing a fractal based on derivative noise would be a bit clunky. I'm trying to change that, since including derivative calculations in the noise generation itself, while the noise value is being generated, opens up some pretty interesting effects.
The main reason it isn't currently supported is that the ANL library supports a number of different basis generators: classic value noise, Perlin noise, Perlin simplex noise, cellular noise, white noise, etc... Yet the library was designed such that the various variants can be mixed and matched, combined willy-nilly as the various octaves of a fractal, if so desired, and the various generators were not built to calculate derivative of the function alongside the value. In fact, figuring out how to calculate derivative alongside noise for some of the generators could prove quite tricky indeed.
A second reason lies in the process of calculating the derivative itself, and in how the noise is generated in ANL. You see, the typical Perlin classic noise process involves generating the noise values at the corners of a cube (or N-dimensional hypercube) and interpolating the values across the volume to obtain the final value. In ANL, I represented this as explicit interpolations; ie, 3D interpolated noise looks something like this:
It's not pretty, what with the exponentially exploding number of function arguments required as axes are added, but it has its own sort of elegance given the problem. The problem, you see, is that of exponential increase of the number of interpolations necessary to generate noise. In 2D noise, 4 points are interpolated using 3 interpolation operations. In 3D noise, 8 points are interpolated using 7 operations. Since ANL supports up to 6D noise (in order to perform distortionless seamless mapping in 3 dimensions as per the technique I detailed at http://www.gamedev.n...seamless-noise/ ), the 6D version requires interpolating 64 different points using 63 different interpolation operations.
Handling the 4D and 6D variants using the chained function call method is a bit less messy and more elegant than the alternatives, and I chose that approach to save myself some work.
However, in order to calculate the derivative of noise in-place, this chained function calling doesn't work that well. In Quilez's method, the interpolation operations are combined into a single large operation which is simplified/expanded to a polynomial from which coefficients are extracted, the coefficients being various combinations of the lattice sample point values. The polynomial has a number of terms equal to the number of sample points being interpolated. The coefficients of the terms are calculated from the sample points. Solving the polynomial for a given set of inputs and the enclosing sample points results in the noise value of the input point; and you can take the derivative of the polynomial with respect to a given input coordinate axis in order to calculate the derivative of the noise in that axis at the same time. Doing it this way is very simple and straightforward.
The problem, of course, is that while it is easy enough to construct the polynomial by hand for the 2D and 3D cases (and not extremely overwhelming even in the 4D case, though a bit of a pain), the 6D case is a mess. If you are going to try to simplify and expand a 64 term polynomial from a series of nested linear interpolation operations, you are going to need a rather large sheet of paper, and you are probably going to make mistakes. And especially considering how far in the past my last calculus course was, and how extremely rusty and shaky my skills are these days, doing that by hand just didn't seem appealing at all.
I don't have any kind of software to help me with these kinds of things. I'm sure such software exists, but if so, it's nothing I'm familiar with. So for a long time now, I've avoided trying to implement in-place derivative noise in ANL. Recently, though, I've had some free time at work to tinker with things, and I decided to try my hand at writing some Lua script to help me construct the 4D and 6D polynomials and extract the terms. It has required me to delve into old college textbooks on expression parsing, converting infix to prefix notation, etc... In other words, it's been fun.
I started with writing a function to convert an expression from infix to prefix. This, of course, is standard algorithm textbook stuff and only took a couple minutes to implement. I didn't want to write a general-case parser, so I imposed the requirement of the tokens in the input string being white-space delineated, so that I could extract the tokens using simple string.gmatch calls in Lua. The function outputs a string in prefix notation, whitespace-delineated.
From this string I build an expression tree. Each node of the tree can be one of three basic types: operator, unknown operand, or coefficient operand. The unknown operands correspond to the input coordinates (x, y, z, etc...), and the coefficients are derived from the sample points (A, B, C, etc...). I instituted this division between operand types to facilitate simplifying and expanding the terms, since I want to be able to isolate terms in regards to their unknowns, and group together the coefficients.
Once the expression tree is constructed, I recurse the tree and apply a few simple rules. The rules are as follows:
1) If a node is a "*" operation between two operands of type unknown (x, y, z, etc...), I concatenate the node into a single node of type unknown operand whose value is a string concatenation of the two initial unknowns. IE, if a node is represented as "x*y" then the result of this simplification will be a single token "xy".
2) If a node is a "*" operation between an unknown and another "*", and that right-hand side is a "*" with an unknown on the left, then the unknowns are concatenated. ie, a chain of nodes such as "x*(y*a)" reduces to "xy*a".
3) If a node is a "*" operation with the left operand being a coefficient and the right being an unknown, the two are swapped.
4) If a node is a "*" operation with the left side being an operand and the right side being a "+" or "-" operation, then the distributive property is applied and a new sub-tree constructed. ie, if the tree is equal to "a*(b+c)" the result will be a new tree "a*b + a*c".
The above rules are applied to the expression tree to construct the polynomial. However, the polynomial is still nested, requiring another pass to invert +/- signs of operations for subtraction operations. ie, a tree of the form "a-(b+c)" will be converted to "a-b-c". This involves recursing the tree and if a node is a subtraction operation, the right child tree will be recursed and every + changed to -, every - changed to +. At this point, the tree is actually technically "broken"; if you interpret the tree as-is, using tree traversal, the terms will be slightly wrong, if you assume that a sub-tree represents a sub-expression that can be parenthesized. However, if you traverse the tree and do a simple leftchild-root-rightchild traversal, outputting terms, then the result is our simple polynomial; or, rather, the polynomial before coefficients are grouped.
Here are the functions for the simplification operations:
The expansion/simplification is a bit hackish, but it works.
In order to construct the initial expressions that start the process, I wrote some helper functions which perform the chaining interpolation and generate an expression string:
For the purposes of illustration, here is a traversal of the expression tree generated by lerp6, before expansion into polynomial form.
Go ahead and try to reduce/expand/simplify that by hand into a polynomial. I dare you.
Expanding the tree results in the creation of terms, although it doesn't result in the final polynomial form. To do that, I had to come up with a way to combine terms that share unknown operand groups. If you are curious, here are the terms generated from lerp6, after expansion and simplification but before grouping coefficients:
If you were wondering, there are 728 separate terms there.
The next step is to combine the groups so that I can extract the coefficients of the terms. The way I elected to do this was to traverse the tree and build a string of the above terms (+/- unknown * coeff) delimited by some spacing character, "#", between the terms and whitespace between the sign, unknown, "*" and coeff portions. This string then is easy to parse into the individual terms.
From each term, the group is extracted and a bucket in a group list is created to hold that term if none exists. Then, a table is made from the sign and the coefficient, and this is added to the group list.
The result of this is a bucket full of groups, each group named for the unknowns (or interpolants) that are used in the term, and the elements of each group are the coefficients (grid samples) that contribute to the term. Here is the result of iterating the group list for lerp3:
y: +c-a x: +b-a zx: +f-e-b+a z: +e-a yx: +d-c-b+a zyx: +h-g-f+e-d+c+b-a zy: +g-e-c+a
From where I sit, those coefficients look correct, and are exactly the same as the coefficients in Quilez's code. I have verified the 2D and 4D cases as well, and they produce the expected results.
Once I have the terms of the polynomial, it is easy enough to build the actual code expressions to calculate the noise value, and construct expressions to calculate the derivative. The derivative, of course, is calculated with respect to a particular unknown, or axis, so the derivative of a noise function will be a vector of dimensionality equal to the dimensionality of the noise function. The final noise function can then be constructed. Here is the final 6D noise function:
The noise generation part works. I haven't yet gotten around to testing the derivative calculation. I'm sure it works, but the next step is to figure out how to re-design the basic modules in ANL to incorporate derivative calculation when (and only when) desired, in order to set up a proper test. That, of course, will take a bit more doing since I need to figure out how to calculate the derivative of simplex noise, cellular noise, etc..., as well as the derivatives of the combination/selection functions.