This page comes about because of two things: an assignment for my computer graphics course, and a bit of inspiration from a computer graphics YouTuber, acerola, who recently did some work with fractal motion. Though the page was actually a requirement for the course, so, you tell me what inspired this page more :P
The theme for our competition in our course was "Chaos in Harmony". If we can find a good, stable random function that we can efficiently raytrace over, that'd be really nice. The chaos would literally be in harmony.
In the video linked above, there was a claim along the lines of, "we can't draw a good bounding box for all the points, because there is no general solution for this". And I wondered: are there specific cases where we can? Could we use that to implement a raytracer?
A bit of background
I knew I wanted to do something with fractals for my project, so I turned to a book: "The Science of Fractal Images". It was in this book that I encountered a method for fractal Brownian motion (effectively, noise produced by random choices) that operated by splitting a line at its midpoint, over and over, and applying a decaying random displacement to the line at each of these points.
To illustrate this effect, you can see a simple implementation of this pattern below, iterated through its steps:
The implementation here is actually quite simple, and differs from the book in a number of ways. Namely, the random
displacement at the midpoint is determined as a random offset chosen below some noise N, minus N/2. At each
iteration, N is multiplied by some decay factor r, so the real displacement is in [-(r^I)N/2, (r^I)N/2] with I
iterations. This allows us to very easily find the actual bounds on the height of displacement after infinite iterations
by simply applying the corresponding geometric sum formula.
This was the first iteration of my implementation. We can easily generalise it to higher dimensions by identifying that
the number of points we need to compute is always the number of 2^D-1 where D are the number of controlling
dimensions (i.e. the dimensions we split across, not the final height dimension we produce by random generation).
Conveniently, all the midpoints we need to split are defined by the binary encodings of 1..=2^D, where an offset in
each dimension is either present as identified by 1 in binary or not present when identified by 0. This works because
this successfully identifies all combinations of dimensions we need to split, e.g., in 1 dimension, along the line; in
2 dimensions, along both axes and the diagonal away from the origin; in three dimensions, each axis and diagonal along
a cube; and so forth.
And, like before, we can visualise this in a gif:
Since my initial implementation was actually in Rust, we can very easily take this and actually just slam it into the site. You can play with the Brownian noise generator below:
Loading WebAssembly...
The corresponding source code is on GitHub and is generalized for an arbitrary number of dimensions, so e.g., you could have a 3 dimensional surface and 1 time dimension and slide through time to see a smoothly transitioning surface. I ended up not being able to implement a ray tracer for my original version in a suitable time to do this, though.
Okay, but why are we going through all this effort?
The original problem was that we want to generate lots of noise down to a very precise point efficiently and without precomputation of the surface. But above, you can see that we slow down a lot as we gain more iterations, but without much increase in visually perceivable data. What we would like, is to be able to perform this midpoint splitting only as necessary for actually intersecting the surface and compute the surface only to the degree which is necessary for doing so. After all, we can't see all the points if we intersect the surface somewhere earlier.
The trick here is that our noise generation is stable. For any point on the surface, we compute it by this random process of splitting at the midpoints and applying a random displacement. The key to having stable generation of this surface is to notice that, no matter what, each point always has the same sequence of midpoint splittings to compute it, so as long as our random generation is stable per midpoint then our whole surface is stably generated. To accomplish this, our random displacement function is defined as the hash of the endpoints of the line at each iteration and the seed selected at the beginning. This is actually exactly the method used above.
Then, since we know the general upper bound after a given number of iterations, we can simply treat the generation process itself as an acceleration structure; we know that, as a ray approaches the surface, it will pass through each theoretical upper bound as it approaches the surface. We can use this to sparsely compute the surface by only computing as much as is necessary to know that the ray definitely hasn't yet hit the surface, and then finally compute the set of triangles on the surface at the final iteration.
This is what ultimately led to my first renders, using two of these surfaces, a sphere, and I think this environment map (my dealer won't tell me where he gets it). I particularly like the ocean aesthetic.
And, this already looks surprisingly good, but it looks dull. The sun doesn't reflect in the way I'd like -- note no specular effects on the surface of the water anywhere but the points which are immediate reflections. This is kind of disappointing, and mostly is because we are currently using a ray-caster -- the effects from the light towards the eye are somewhat missed, and we need bidirectional ray-tracing for this. Beyond that, we would also like to have our image be a bit more striking.
Optimisation interlude
For a while here, I got a bit stuck in the weeds of optimisation. Just for personal edification, I added them here:
Binary nonsense
Because our iterations are encoded as midpoint splitting, our distance halved every time. At the first iteration, the whole surface can be thought of as from 0 to 2^32 in each direction. Then, in the next step, we halve each of these such that our midpoints fall on 2^31. Then, our midpoints fall along 2^31 plus-or-minus 2^30, and so on and so forth. When we cross a boundary between one region with a given upper bound and another region with another upper bound, we have to recompute the upper bound, and we don't want to recompute it all the way from the beginning iteration. We use the number of trailing zeros to identify how many iterations we've already computed as we enter the next region, and then continue, as the points which define the edges of the next zone internally contain the information as to how many iterations have already been computed. This saves about 40% of the compute time.
Visibility falloff
As you get further away from something, the less precise your view of that thing becomes. We can thus apply an optimisation to simply stop computing at a certain distance away from the camera, since we won't see it anyway, which conveniently can be computed by simply extracting the exponent out of the floating point value for the (squared) length of the vector from the camera to the point in space that we want to know the iterations-necessary-for-precise-viewing for, shifting it right (because division after logarithm is the same as a root before the logarithm), then subtracting the log of our intended falloff distance -- which in effect gives us a precision scaling method with zero floating point operations, since we can precompute the log of our intended falloff distance. This saves another 15-30% off of the compute time, especially when looking at the horizon.
"Just don't cache lol"
As we compute more rays out from the camera, one might be tempted to cache this information across rays. In practice, though, we end up using a couple of hundred GBs for our render with basically zero speedup, so we just don't cache our surface between intersection events at all.
Final image
Now, we've implemented bidirectional path tracing using SPPM, which ensures that our surface "lights up" even with a small number of light events actually hitting our surface where we are. I had originally implemented this to highlight caustics formed by the noise function, but ended up struggling to consistently get an image that had interesting caustic effects.
Instead, I went for a space scene. I found the lighting here to be particularly ominous, especially as the light blasts away much of the colour on the surface of our noise while lighting from random reflection allows us to continue to see the regions normally obscured from the light. In order to ensure that my scene received a majority of the light, I even implemented a spotlight which actually produces the light from the "sun" in the background. My final render is below:
You can also find a JPEG version and an EXR version.
What's particularly nice is, that we don't need a lot to specify our scene. The whole scene description, and all of its assets, are defined below (with comments):
<test type="image" id="brownian_space">
<integrator type="sppm" depth="4" photons="10000000" iterations="10" radius="16" updateInterval="1">
<scene>
<camera type="perspective" id="camera">
<integer name="width" value="2048"/>
<integer name="height" value="1024"/>
<string name="fovAxis" value="y"/>
<float name="fov" value="30"/>
<!-- We look from 2050,10050,10 towards the farthest corner of our surface. It's huge. -->
<transform>
<lookat origin="2050,10050,10" target="4194304,4194304,-229376" up="0,0,-1"/>
</transform>
</camera>
<!-- This is the "sun". You can find the sphere definition way below, but this allows us to constrain our lights. -->
<light type="spotlight" position="4194304,2097152,524288" direction="-4192254,-2087102,-524278" angle="0.1"
power="4194304"/>
<!-- This is the sand, which is moderately perturbed and lighter than the mountains behind. -->
<!-- Note the "interval", which defines the number of initial iterations which are forced to be zero. -->
<!-- Note the "vislimit", which defines the distance from the camera at which the surface detail halves. -->
<!-- Note the "lower" and "upper", which defines where the surface is clipped (in shape-space) to avoid
computing regions too far away that don't contribute to the image. -->
<instance>
<shape type="brownian" decay="0.5" interval="15" vislimit="65536" lower="2048000,10240000,-10240"
upper="6144000,12288000,256000"/>
<bsdf type="diffuse">
<texture name="albedo" type="constant" value="0.71,0.41,0.20"/>
</bsdf>
<transform>
<!-- The scale on the z axis makes our sand actually only very lightly perturbed. -->
<scale x="0.0009765625" y="0.0009765625" z="131072"/>
</transform>
</instance>
<!-- This is the mountain, which is strongly perturbed. -->
<instance>
<shape type="brownian" decay="0.5" interval="10" seed="15" vislimit="65536"
lower="2048000,10240000,-10240"
upper="6144000,12288000,256000"/>
<bsdf type="diffuse">
<texture name="albedo" type="constant" value="0.49,0.27,0.19"/>
</bsdf>
<transform>
<!-- The scale on the z axis makes our mountains jut out high above the surface. -->
<scale x="0.0009765625" y="0.0009765625" z="1048576"/>
</transform>
</instance>
<!-- Wait, the sun is just an emissive sphere? -->
<!-- Always was. -->
<instance>
<shape type="sphere"/>
<emission type="lambertian">
<texture name="emission" type="constant" value="0.8,0.8,0.8"/>
</emission>
<transform>
<scale value="524288"/>
<translate x="8388608" y="4194304" z="524288"/>
</transform>
</instance>
</scene>
<!-- Embrace the noise. It lends to our ominous aesthetic. -->
<sampler type="independent" count="16"/>
</integrator>
</test>
I ended up using a slightly wider 2:1 aspect ratio as I found it to be more ominous, and I'm pretty confident my SPPM implementation was bugged and led to things being washed out, but in this case it actually looked rather cool!
Last bits
I worked on this project alone. I used "royal we" above because, hell, I'm an academic, I'm used to saying "we" for research!
As for a title? How about, "Space Brownie"?