An OSL wood knot shader

Back in january I presented a wood shader that produces decent results but one of things that was missing was a way to add realistic knots to planks. In this article a present a shader that is meant as a first step in producing knots. It is not a finished shader yet, but it is fully functional as a texture warping tool and in a future article I might detail how to combine it with a wood shader and a shader to texture the inside of the knot properly. Because that will take some time so I thought it better to present it as a WIP as I am currently unable to spend more than a few minutes behind my desk (and there's no Blender for Android alas).

Warping space

The idea is quite simple and might be adaptable to more than just producing knots: we distribute points randomly and warp the position coordinates around these points. These warped coordinates are then used as the basis for a texture. In the example image we use the warped coordinates as the input to a bands texture:

As you can see the bands of the texture seem to flow around the knots as if they were repelled by them and in fact that is pretty much how the algorithm works.

Implementation


int bend(vector p, vector k, float r, float a, float m, output vector B){
vector D = k - p;
float L = length(D);
if( L < r ){
float c = L/r;
float d = m * pow( 1 - c , a);
if( d < L ){
B = d * normalize(D);
return 1;
}else{
B = D;
return 2;
}
}
return 0;
}

shader knot(
vector Pos = P,
float Scale = 5,

float R = 0.8,
float Falloff = 1,
float Strength = 0.9,
float Knots=0.5,

output vector Vec = P,
output float Fac = 0
){
vector p = Scale * Pos;
vector sdp = 0;

float TR = ceil(R);
for(float dx=-TR; dx <= TR; dx++){
for(float dy=-TR; dy <= TR; dy++){
for(float dz=-TR; dz <= TR; dz++){
vector ip = floor(p)+vector(dx,dy,dz);
for(int ik=0; ik < (int)Knots; ik++){
vector k = noise("cell",ip,ik);
vector dp= 0;
int ret = bend(p,ip+k,R,Falloff,Strength,dp);
if(ret != 0){
Fac=max(Fac,ret==2);
sdp+=dp;
}
}
if( noise("cell",ip,-1) < mod(Knots,1.0) ){
vector k = noise("cell",ip,-2);
vector dp= 0;
int ret = bend(p,ip+k,R,Falloff,Strength,dp);
if(ret != 0){
Fac=max(Fac,ret==2);
sdp+=dp;
}
}
}
}
}
if( Fac < 1 ){
Vec = p + sdp;
}else{
Vec = sdp;
}
}
The magic is in the bend() function. It calculates the difference vector from the point being shaded to the center of the repulsion. If the distance is short enough, a translation vector is returned that is a distance dependent fraction of the difference vector. It takes some thinking to see that in order to create the illusion of repulsion we actually have to replace the point being shaded towards the center of repulsion: remember that at the point we are shading we want to see the lines closer to the center.

When we calculate the translation towards the center and we are close enough to it, we might overshoot the center point. In which case we return a value of 2, signalling we are inside the knot and return not the translated shading point but the vector pointing to the center (which might be useful to render some concentric pattern in the knot).

The shader itself is nothing more than checking if we might be in range of one of the randomly distributed points and calling the bend() function to do the actual work. We do allow for a fractional number of knots per unit cell, in which case the fraction acts as the probability that a knot might occur.

Room for improvement

Obviously knots are not spherical marbles embedded randomly in some wood but as the remnants of branches they are more reminiscent of stubby cylinders which might be better modeled by determining the distance to a randomly oriented line segment instead, something I intend to implement in the near future.

Another area that needs attention is to way the inside of the knot appears. In the example setup we produced some concentric circles but a real knot is a bit more complex than that.

OSL debugging tips

The other day I was debugging a shader that didn't produce the results I expected, in fact it didn't produce any output at all. In its basic form it resembled the following shader:


shader S(
vector position = P,

output color Color = 0
){
... many calculations ...
if( ... some condition ...){ Color = 1; }
}

When I hooked this shader up in some node setup, I'd expected it to produce a pattern but all I got was a uniform black. Obviously the condition was never true and Color was never assigned a value of 1. So I had to determine what the values of some intermediate values were to get an idea why the if statement didn't fire.

So I inserted a printf() call:


shader S(
vector position = P,

output color Color = 0
){
... many calculations ...
printf("%.2f\n",position);
if( ... some condition ...){ Color = 1; }
}

I clicked the recompile button, looked in the console and ... nothing was printed! To be 100% sure that the shader was evaluated I changed it to:


shader S(
vector position = P,

output color Color = 0
){
... many calculations ...
printf("%.2f\n",position);
if( ... some condition ...){ Color = 1; }
Color = 0.5;
}

Now the node setup produced the expected shade of gray, so it was hooked up correctly and was indeed being executed but more surprisingly, the console showed the printf output.

This baffled me: if I commented out the final assignment to Color, the printf() wasn't executed. I thought I found a bug in OSL and reported it. But was it?

Constant folding

It turned out that this is intended behavior. The OSL compiler is a clever piece of software and does its best to reduce unnecessary calculations. One of the tricks it employs is called constant folding, a method in common use in many compilers. Basically it tries to calculate at compile time anything that can be compiled, for example:

twopi = 2 * 3.14;fourpi = 2 * twopi;

becomes after compilation:

twopi = 6.28;fourpi = 12.56;

This keeps your code readable but reduces actual calculation at runtime. And this can lead to significant speedups when you are shading milions of samples.

OSL has a lot more of these tricks up its sleeve and can reduce quite complex expressions, even removing branches from if/else statements. But it goes a step further still. If the compiler detects that none of the output values is changed from its default by a shader, this shader is not executed at all.

For my example shader this implies that because the condition in the if statement is apparently based on things that can be calculated at compile time and never true, the if statement is optimized away. The result of this is that no output value will ever be changed by this shader and so the whole shader is optimized away, even if we put printf(), warning() or error() calls in this shader, because statements that produce side effects only and do not contribute to the output of the shader are not by themselves considered enough to keep an if/else branch or even a shader from being optimized away.

Debugging OSL code

Now if we want to debug some code we might have a problem because simply inserting a printf() call is no guarantee that the information we are interested in will be printed. Also, assigning a value to an output socket might not be enough:


shader S(
vector position=P,

output color Color=0
){
...
printf(.....);
Color = 1;
Color = position[0];
}

Neither of the last two statements will ensure that the printf() call is executed. The first one us an assigment of a constant and can be calculated at compile time. The last one only depends on the input and can be 'hardwired', that is, there is no need to execute the shader; the compiler can simply connect input and output.

So if we want to make sure the printf() executes, we need an assignment to Color that cannot be calculated at compile time, for example one of:

Color = position[0]>0;
Color = rand();

Another option is to create an extra output socket just for debugging and assign the variable that you are interested in to this socket. If you connect it to a diffuse shader there is a lot you can tell based on the color, but of course not all variables are suitable for this approach (strings to name one).

New debugging options

From release r57660 (download it from the daily builds, I tested the 64bit windows version) you can disable OSL optimizations by assigning a value to the OSL_OPTIONS environment variable. If you use Windows this can be done like this:

  • open a command prompt
  • set the variable (don't use quotes)
    set OSL_OPTIONS=optimize=0
  • run Blender from the command line
    C:\path\to\blender.exe
  • don't forget to recompile your OSL shader
On unix-like systems you start from an xterm or equivalent and type export OSL_OPTIONS="optimize=0"
path/to/blender
(if you do not use bash but tcsh use setenv instead of export). Setting these variables is local to the command prompt/xterm you are working in so when you exit Blender and close the command prompt the optimization override will be gone. Changing the enviroment on a more global level is not recommended. Also when you're done it might be a good idea to remove all printf calls from you production code because printing to a console is extremely slow. You might still want to print information in exceptional situations but you might consider using error() for that purpose.

Conclusion

In many cases using printf() calls to debug your OSL shaders is still the method of choice but if nothing shows up on the console make sure that these calls are not optimized away.

Working with the incidence vector in OSL

The fact that in OSL we have the incidence vector I available has proven its utility already in the many cool toon and carpaint shaders people have implemented. In this article we take this concept a bit further and create a shader that does not change it lighting model as a function of the angle between incidence vector and normal but changes its complete behavior depending on the angle between incidence vector an object orientation, allowing us the model a 3D-object (Venetian blinds) with just a 2D-tecture.

Lenticular printing

The idea of providing a different image based on the viewing angle is not new and is fact readily available in traditional media under the name Lenticular printing.But in OSL we are able to act on any incidence vector, not just the one from the camera. This means for example that shadow rays act as they should and that we can make our material transparent in places to let it cast correct shadows.
In the sample image and in the video we used this trick to render venetian blinds as seen from different position with just s single plane.
The concept is simple: the horizontal slats of the blinds appear taller when seen from an angle so we calculate the angle between the incidence vector and the x-axis (because unfortunately there is no attribute object:rotation to determine the up axis of an object) and divide that by maximum angle that would let us see something through the blinds. That fraction is used to paint part of a horizontal strip representing a slat. Finally we replace the shading normal by one that is perpendicular to the slat allowing for s bit of a curve across the slat:
surface venetianblinds(
  point Pos = P,
  float Spacing = 0.1,
  float Width = 0.1,
  float Curve = 0.2,
  
  output float Fac = 0,
  output float Top = 0,
  output normal Normal = N
){
  Top = dot(I,vector(0,0,1))>0;
  
  float maxangle = atan(Spacing/Width);
  float angle = abs(M_PI_2-acos(dot(normalize(-I),vector(0,0,1))));

  float height = Spacing;
  // TODO we should look at the sign here to prevent 'flip'
  if (angle < maxangle){ height *= angle/maxangle; }

  float z = mod(Pos[2],Spacing);
  Fac = (z < height);
  
  if(z< height){
    Normal = normalize(vector(0,0,1)+N*Curve*(height/2-z)/height);
  }
  
}

Example node setup

The node setup used to render both the still and the video looks like this: This material is applied to a simple plane. The lighting consists of a single mesh light on the left.

Rayleigh scattering in Cycles

Prompted by a question in a comment by Joakim Poromaa Helger I came up with this approach to emulate Rayleigh scattering (the effect that turns the sky blue and the sun red at the same time) in Cycles: by combining a transparent and a glass shader.


In the picture above there is a bright white sun on the right. In a material that exhibits Rayleigh scattering (like air) blue light is scattered more than red light. That means that light traveling from a lightsource through such material will be tinted red because blue light is scattered to all sides. In a direction away from the direct path of the light we see just this scattered blue light and not (or much less) the direct light.

Now if we combine a red(ish) transparent shader with a blue(ish) glass shader, the transparent shader will cause some of the light to be passed through the object without regard for index of refraction and give it a redish hue, while the glass shader will pass some (refracted) blue light through and reflect some. (In the sample picture I have turned off caustics because that takes forever to give a good quality render but also because in this setup caustics exhibit a blue color while it should be redish as well). The node setup looks like this:

With some extra trickery we can correct the color for any caustics as well:

(Note that this image took 10,000 samples and the caustics are still very noisy. Clearly caustics is not Cycles strong point. The light/shadow inside the orb are due to the fact I placed it slightly below the ground but am to lazy to fix it).

We achive the caustic coloring by arranging that the camera sees the blue glass shader while any object receiving diffuse light (like the floor) sees a red glass shader, a distinction we can make with help of the light path node:
Now obiously the effect here is exagerated for demonstration purposes and arguably it is not a real simulation of Rayleigh scattering but for artistic purposes this might be a good start.



An OSL starfield shader for Blender Cycles

Inspired by a BlenderHD video I implemented a simple starfield shader in OSL. The distribution of stars can be controlled by a density map and the color of the stars follows (very roughly) the probability of colors for visible stars. As I will show in an example, adding dust clouds and twinkling can be done with procedural textures.

Starfield example

The image below is a still from a short animation.

(It is pretty dark, for cinematic effects you probably would choose a far brighter setup)
The stars are generated with the starfield shader and their distribution is controlled by a noise texture. The faint green glow of stellar dust and the faint orange color on the horizon are generated with a separate noise and gradient texture respectively. (example node setup at the end)
The animation features additional twinkling that is not visible in the still image. This is effected by blending time animated 4d simplex noise.
(In the animation the contrast/lightness was autocorrected for a better view)

The code

The code below contains two functions: one that calculates a probability dependent voronoi distribution and the shader itself. The code is pretty straight forward but note that the Falloff value should not be set too high to prevent artifacts in the star shapes when seen in close-up. (These are due to the fact that we divide space in cells and the fact that we only considered the nearest neighboring cells for our voronoi distribution. Taking into account further neighbors would solve this but impact the performance of the shader in a rather drastic manner).
The outputs of the shader consist of a color and a factor. The color is the emmisive power in W/m2. This could be used as an input to an emmision closure but in fact it might be simpler just to normalize this and use it as a color directly (as is done in all examples). The factor is just a value which is 1.0 for the center of the stars aand decreases exponentially when further away.
point voronoi3dp(point p, float density, output float d2)
{
 int xx, yy, zz, xi, yi, zi;

 xi = (int)floor(p[0]);
 yi = (int)floor(p[1]);
 zi = (int)floor(p[2]);
 
 float dbest = 1e10;
 point pbest = 1e10;
 vector dz = vector(7,111,19);
 for (xx = xi - 1; xx <= xi + 1; xx++) {
  for (yy = yi - 1; yy <= yi + 1; yy++) {
   for (zz = zi - 1; zz <= zi + 1; zz++) {
    vector ip = vector(xx, yy, zz);
    if(cellnoise(ip)< density){
     point  vp = ip + cellnoise(ip+dz);
     vector dp = p-vp;
     float  d  = dot(dp,dp);
     if (d < dbest) {
      dbest = d;
      pbest = vp;
     }
    }
   }
  }
 }
 d2=dbest;
 return pbest;
}

shader stars(
 point Pos = P,
 float Stardensity = 0.001,
 float Scale = 1000.0,
 float Falloff = 10,
 
 output color col = 0.0,
 output float fac = 0.0
){
 point p = Pos * Scale;
 
 float d;
 voronoi3dp(p,Stardensity,d);
 fac = exp(-d*Falloff);
 if(d<100 br="">  col = blackbody(3500.0+6000.0*cellnoise(p+vector(17,17,17)));
 }
}

Example node setup

The node setup used to generate both the still image and the animation is shown below. It's fairly large but not really complicated.

(click to enlarge)
Each texture contribution is ordered in a horizontal row of nodes. From top to bottom: the starfield (dark blue, note the normalization of the color), the twinkling (light blue), the dust (green) and the horizon glow (yellow). The glare was added in the compositor along with some contrast enhancement (for the video sample I added extra contrast enhancement in an external program). An example .blend file bundled with the stars.osl script and the simplexnoise.osl script is available as a .zip file on GitHub. (click View Raw button to download)

100k visits, hurrah!



This blog has reached the 100k visits mark just a couple of hours ago. Not bad for something that started just a few months ago :-) Apparently you all find it worthwhile reading material so I will definitely continue writing. (and to be honest, I feel flattered by all the nice comments I get, the Blender community is a really friendly bunch of people :-) So a big thank you to you and onward to the next milestone!

A new tree add-on part VI: speed improvements for complex trees

Trees that contain many branchpoints take a progressively longer time to generate as they become more complex. This is due to the time it takes to determine which branch segment is closest to any endpoint. The code up to version 0.0.7 did a linear scan through all the branch segments for each endpoint and while this was still fast enough for moderately sized trees (the time it takes doesn't explode exponentially because each iteration removes endpoints within kill distance of a branch) it took many seconds for a tree with thousand endpoints, short internode lengths and a short kill distance. In other words it got slow when the trees became interesting.

Version 0.0.8 (available on Github, click view raw to download) adds a significant speed up for larger trees (3-6 times depending on parameter settings). Quite complex tree skeletons can now be generated in a few seconds instead of half a minute or more (depending on the power of your pc of course. And remember, this is a pure python implementation so having multiple cores makes no difference, it's the raw speed of a single core that counts). Please note that because endpoints within kill distance are removed in a different order, trees generated in version 0.0.8 are not identical to trees generated in previous versions but their general appearance stays the same.

For the technically inclined: this speed up is realized by storing all branch segments in a kd-tree, so looking for a nearest neighbor now scales as the log of the number of branch segments instead of taking linear time. This is a very noticeable difference once your tree grows to hundreds or even thousands of branch segments. The kd-tree is a pure python implementation in its own, separate module kdtree.py that might be interesting for other applications as well.