Small Blender Things

Drawing cubic splines in OSL

Some time ago I wrote about a routine to calculate the closest distance to a quadratic bezier spline and recently I reused this routine again in a slightly optimized form to create wiggles.
Being able to calculate the closest distance to a quadratic spline is useful but I really wanted to use cubic splines because they offer more control, allowing for easier ways to create nice shapes.

update: made the code 10% faster by making it 2D only. The GitHub version relects that change but original code can be used by removing a single macro
Warning; long read ahead :-) The first section describes how to use the shader and where to download the code and is pretty short. The second section dives deep into the code and is definitely not an easy read, sorry for that but fortunately there is no need to go that deep if you know what a cubic spline is and just want to use the code.

If you would like to know more about programming OSL you might be interested in my book "Open Shading Language for Blender". More on the availability of this book and a sample can be found on this page.

The shader

The approach to finding the distance to a cubic spline is the same in principle as for a quadratic spline: find all values for t that satify the relation dB(t)/dt . (B(t) - P) = 0, or in English: those values of t whre the tangent to the spline is perpendicular to the vector pointing to the point being shaded (remember that the dot product of two vectors is zero if they are perpendicular).
For a quadratic bezier curve the equation to solve is cubic and is quite simple but for a cubic bezier curve this is a fifth order equation for which there is no analytical solution. And although cubic splines are used everwhere (for example in fonts) a numerical solution for this problem was surprisingly hard to find.
Fortunately after some serious searching the Graphic Gems series of books proved again to be incredible useful, specifically "A Bezier Curve-Based Root-Finder" by Philip J. Schneider.
The original code is in C so it neede quite some work to convert it to OSL. Normally converting C code to OSL is pretty straightforward but this code uses multidimensional arrays and recursion, both unknown to OSL. On the plus side, OSL's native point type and built-in dot() make for a lot less code.

Code availability and node setup

My implementation (with a sample shader in the same file) is available on GitHub. The node setup used to create the heart shape is shown below, it is made of two mirrored splines:

The values for the control points are:P0 = ( 0.5, 0.1, 0.0)
P1 = (-0.5, 0.7, 0.0)
P2 = ( 0.4, 1.2, 0.0)
P3 = ( 0.5, 0.8, 0.0)

Code explanation

I am not going to list all of the code because it's way to long and rather boring but I will tell some more about the FindRoots() which is nontrivial and needed to be reimplemented in a non-recursive form.
The general idea of the shader to cast the fifth order equation whose roots we'd like to find in a form that is itself a bezier curve (I.e. a bernstein polynomial) and exploit the specific properties of a bezier curve to efficiently find the roots.
The mathematical specifics are explained very well in Schneider's article so let us just look at the code:

#define POP wtop+=6;wn--

int FindRoots(point w0[6], float t[5]){
  point w[STACKSIZE];
  int wtop=0;
  int wend=0;
  int wn=0;
  for(int k=0;k<6 br="" k="" w0="" w="" wend="">  wn++;
  int nc=0;
  while(wn>0){
    //printf("[%d]",wn);
    int cc=CrossingCount(w,wtop); // w[:5]
    if(cc==0){
      POP;
      continue;
    }
    if(cc==1){
      if(wn>=MAXDEPTH){
        t[nc++]= (w[wtop][0]+w[wtop+5][0])/2;
        //printf(" MAX ");
        POP;
        continue;
      }
      if(ControlPolygonFlatEnough(w,wtop)){ // w[:5]
        t[nc++]= ComputeXIntercept(w,wtop);
        //printf(" FLAT (%.1f | %.1f)",w[wtop],w[wtop+5]);
        POP;
        continue;
      }
    }
    // not flat or more than one crossing
    DeCasteljau(w,0.5,wtop,wend);  //left=[wend:wend+5], right=[wend+6;wend+11]
    for(int k=0;k<6 br="" k="" w="" wend="" wtop="">    wn++; wend+=6;
  }
  return nc;
}

A fifth order bezier curve has 6 control points and at most 5 roots. The control points are passed as the w0 argument and the t array will receive the root values. The function's return value will be the number of roots found.
The central idea of this algorithm is to find a segment of the curve that crosses the x-axis just once and then determine the intercept, i.e. a root. To achieve this the spline is subdived into to halves, each again a fifth order spline, until such segment contains a single crossing.
Because OSL doesn't support recursion, we have to maintain a stack of control points ourselves and that's where the w array comes in.
Starting with the for-loop in line 8 we transfer the initial set of 6 control points to the stack and adjust wn to reflect that we now have one set on the stack. At this point wtop (the index of the first control point in the topmost set) is 0, and wend (the index of the first free slot in w) is 6. nc holds the number of roots found sofar.
Next we iterate as long as there is a set of control points on the stack (line 11). We caculate how many times the curve defined by this set crosses the x-axis. If there is no crossing there is no root in this curve so we pop the set of the stack and continue.
If there is a single root and we have recursed deep enough we simply approximate the root as lying halfway between the x-coordinates of the first and last control points and pop the set (line 20). If there is room to recurse we check whether the set of control points is "flat enough" (line 25) and if so compute the intercept with the x-axis of the line between the first and last control point. We store this root and pop the set of control points (see Schneider's article or the original source code to see what is meant by "flat enough").
If the control points are not "flat enough" or we have more than one root in this curve, we split the curve in two using DeCasteljau's method (line 33). The function will produce two sets of six control points at position wend in the w array. The save space we copy the second set to top of the stack (because the control points for the original curve are no longer needed), which means we only have to increase the stacksize wn by one.
Finally, if we run out of the while-loop we simply return the number of roots found.

Speed

Now you might think that all this recursion and the numerous calculations would be very slow but you will be surprised. The function seldom needs to recurse very deep because the accuracy used to determine if a curve is "flat enough" need not be that good (we can't see differences much less than a pixel). On my machine this means that I can tweak the node setup for the example image at interactive speeds if I set the number of preview samples to 1.
There is however plenty of room for optimization, especially in areas were we do calculations with 3D points. If we only want to work in 2D we might save up to 30% on unnecessary calculations.

Creating unique snowflakes with OSL

Falling snowflakes are easy enough to dimulate in Blender with a particle system but if quite a few of them might be seen fairly close up, the real challenge is in making sure each one is unique.

The shader presented in this article will create a unique snowflake for each value of its Seed. Using the object info node we may obtain a unique random number for each object that we may use for this purpose as we will see when we examine the node setup.

Snowflakes

The code below defines the following utility functions;

hex: this function will produce the x and y coordinates of a point as it is rotated through the six segments of a hexagon
in_hexagon: this function will return a non zero value if a point lies within a hexagon
pattern: this function does the real work. It returns non zero if a point is positioned within a number of hexagons that are randomly positioned along the x-axis. It also draws a connecting rod that extends as far as the farthest hexagon.

The body of the shader itself just copies the position being shaded along a hexagon and then calls pattern(). Note that because we expect uv-coordinates for a simple square (I.e in de range [0,1], we transform these coordinates so that the center is at (0.5, 0.5) [line 68].

#define SLOPE 1/sqrt(3)
#define SIDE sqrt(.75)
#define D60c cos(radians(60))
#define D60s sin(radians(60))
 
void hex(float x, float y, output float hx[6], output float hy[6]){
  hx[0]=x;
  hy[0]=y;
  hx[1]=x*D60c-y*D60s;
  hy[1]=y*D60c+x*D60s;
  hx[2]=x*D60c+y*D60s;;
  hy[2]=y*D60c-x*D60s;
  hx[3]=-hx[0];
  hy[3]=-hy[0];
  hx[4]=-hx[1];
  hy[4]=-hy[1];
  hx[5]=-hx[2];
  hy[5]=-hy[2];
}

int in_hexagon(
  float px, float py, 
  float cx, float cy, float r, 
  output float d
){
  d=hypot(px-cx,py-cy);
  if(d>r){ return 0; }
  float hx[6],hy[6];
  hex(px-cx, py-cy, hx, hy);
  for(int h=0; h < 6; h++){
    if((abs(hy[h]) < SLOPE*hx[h]) && (hx[h]< r*SIDE)){
      d=abs(hx[h]);
      return 1;
    }
  }
  return 0;
}

#define CELL noise("cell",seed++)

float pattern(float x, float y, int Seed, int Kernels){
  int seed=Seed;
  int n=(int)(1+Kernels*CELL);
  float hx=0, maxx=0;
  for(int f=0; f < n; f++){
    float hy=0;
    float r=0.2*CELL;
    float d;
    if(in_hexagon(x,y, hx,hy,r, d)){
      return d;
    }
    hx=SIDE*CELL;
    if(hx>maxx){maxx=hx;}
  }
  if(x < maxx && abs(y) < 0.01){ return 1; }
}

shader snowflake(
  point Pos=P,

  int Seed=0,
  int Kernels=15,
  
  output float Fac=0
){
  float hx[6],hy[6];

  hex(2*(Pos[0]-0.5), 2*(Pos[1]-0.5), hx, hy);

  for(int h=0; h<6 br="" fac="=0;" h="">    if(abs(hy[h]) < SLOPE*hx[h]){
      Fac=pattern(hx[h],hy[h],Seed, Kernels);
    }
  }
}

Example node setup

The node setup used to shade the particles in the image at the start of this article looks like this:

The random value from the object info node lies in the range [0,1], so we multiply it by a large value to get a unique integer.The output of the shader node is checked to see if it is non-zero. if so we use a group of shaders to mimic ice, otherwise a transparent shader. The Fac output may vary smoothly so we use a color ramp node with constant interpolation the create stepped values that we can use to drive an ice-like bump pattern. The Kernels input determines how many random hexagons are within each snow flake, so this essentially controles how dense the flakes look.

Code availability

You can copy the code from the listing above but you may also dowload it from GitHub.

Full book review: Blender 3D printing essentials

Last weekend I had the pleasure to read Blender 3D printing esentials by Gordon Fisher.

Gordon previously authored Blender 3D Basics and his experience and subject skill is clearly shown in this new short book 3D printing. The book is Blender specific of course but pretty much 3d-printer agnostic, so it's useful not only for printer owners but also for people who want to use online print shops.

Pros:

short
about 100 pages real content. Clear language as well, so you won't need much time to get up to speed

to the point
it won't tell you much about modelling as such but does show what elements of modelling need extra care when preparing a model for printing and how to fix problems

complete
all you need to know for most common printing materials and printing techniques but of course with new materials becoming available every day you'll have to check the material specifications for things like minimum wall thickness if you want to try a new material but that is not a critique on the book. On the contrary, it does an excellent job of describing things like what influences surface detail, why wall thickness is important and when and how to deal with overhang.

good description of essential tools
Blender's 3D printing toolkit and mesh analysis tools are covered well and the use of the protractor tool (to measure distance and angles in your model) is explained very clearly.

Cons:

uv-unwrapping and texture painting
(for printers that support multicolored models) would have benefited from more in depth attention as these subjects are generally considered difficult by beginners

the price
I personally feel that the €13,56 list price for a 100 page e-book is steep, not to mention the almost €26,- price tag for the print version, but of course online retailers may give you a better deal.

Conclusion

Good book, well written. Based on the content alone I would give it 4.5 out of 5 stars but the unrealistic price tags lowers that to 4.

Book Review: Blender 3D printing essentials

The next few days I will be reading Blender 3D printing esentials and I hope to write up a full review shortly.

wiggles / noodles shader for OSL

The basic idea of this shader is to scatter around wiggly lines that can be used as fibers, hairs or, like in the image below, noodles.

Each line consists of a number of segments connected end to end. Each segment is angled by a random amount relative to the previous one. Segments are not straight though, but curved by a certain amount. The parameters that control the shape of the line are illustrated below, with on the right an indication of what the basic pattern looks like:

The curved segments are implemented as quadratic splines, using routines discussed in an earlier article.

Code and node setup

The code is pretty straightforward. Apart from a rather long list of input parameters it is mainly concerned with calculating a list of segments for each line we want to draw (there may be more than one per cell). The main trick is in line 86, where we make certain that point p1, the point that is used to control the curvature of each segment, lies on the line line through the previous control point and the end point. This ensures that each segment joins the previous one smoothly.

#include "equations.h"

#define DOT(a,b) (a[0]*b[0]+a[1]*b[1])
#define SUB(a,b) vector(a[0]-b[0],a[1]-b[1],0)

// determine if point M is inside a rectangle with a margin
int in_rectangle(point M, point a, point b,
  vector u, float W, vector v, float linewidth){
  point A=a+linewidth*(-u-v);
  point B=b+linewidth*(u-v);
  point D=B+(W+2*linewidth)*v;
  vector AM=SUB(M,A);
  vector AD=SUB(D,A);
  vector AB=SUB(B,A);
  float dotamad=DOT(AM, AD);
  float dotadad=DOT(AD, AD);
  float dotamab=DOT(AM, AB);
  float dotabab=DOT(AB, AB);
  return (dotamad > 0 && dotamad < dotadad) && 
         (dotamab > 0 && dotamab < dotabab);
}

#define CELL noise("cell", cp, seed++)
#define CELL2 vector(CELL, CELL, 0)

shader wiggles(
  point Pos=P,
  float Scale=1,

  int Number=1,

  float Length=0.5,
  float LengthVar=0,
  float Kink=0,
  float Curl=0.2,
  float Wave=30,    // degrees
  int Steps=2,
  float StepsVar=0,

  float Width=0.02,
  float WidthVar=0,

  int Seed=0,
  
  output float Fac=0
){
  point p = Pos * Scale;
  p[2]=0;
  point ip= point(floor(p[0]),floor(p[1]),0);

  int nn=1+(int)ceil(Steps*Length);

  for(int xx=-nn; xx <= nn; xx++){
    for(int yy=-nn; yy <= nn; yy++){
      int seed=Seed;
      point cp = ip + vector(xx, yy, 0);
      for(int wiggle=0; wiggle < Number; wiggle++){
        vector start = cp + CELL2;
        start[2]=0;
        vector dir = CELL2 - 0.5;
        dir[2]=0;
        dir = normalize(dir);
          
        vector perp = vector(dir[1],-dir[0],0);
        float k=0.5 + Kink * (CELL-0.5);
        float c=Curl*(CELL-0.5);
        point p1=start+k*dir+c*perp;
        for(int step=0; step < Steps; step++){
          vector ldir = dir;
          ldir *= Length + LengthVar*CELL;
          point end=start+ldir;
          if(in_rectangle(p, start, end, dir, c/2, perp, Width+WidthVar)){
            float d,t;
            if(splinedist(start, p1, end, p, d, t)){
              float localwidth = Width+WidthVar*noise("uperlin",start,t);
              if(d < localwidth){
                Fac = (localwidth - d)/localwidth;
                return;
              }
            }
          }

          if(CELL < StepsVar){
            break;
          }else{
            p1 = end + (end - p1)*(1+noise("perlin",end)*Kink); 
            start = end;
            dir = rotate(dir, radians(Wave*noise("perlin", start)), vector(0,0,0), vector(0,0,1));
          }
        }
      }
    }
  }
}

The only other issue that needs attention is the generation of random numbers. In each cell we need a number of them and they need to be unique. We therefore add an extra seed argument to the call to noise. However, we must take care that all those numbers are generated in a repeatable way so we reset this seed for each cell to the value provided by the Seed input. This allows us to generate unique patterns for different objects sharing the same material.
The example image at the start was created with a node setup like this:

Note that the Fac output isn't simply 1 or 0 but contains the distance to the edge of the fiber and we use that to drive a bump node (through a square root math node (power 0.5) to give it a smoothly curved appearance). We use the object info node to generate a random number because the heap of noodles consists of three separate squashed half spheres. The shader expects an integer so we multiply the random value by a million to get a unique integer.
One final node of caution: this isn't a cheap shader because calculating the distance to a spline is rather expensive. OSL is quite good at optimizing expressions but still I did dpend quite some time on optimizing the splinedist() by hand. Did did indeed shave off some small percentage but the biggest win was the conversion of all calculations to two dimensions and the test to see if we are within the bound of the control rectangle before actually checking the distance to the spline (line 72 in the code)
A final thing is that most random vectors we generate don't need the z component but OSL has no notion of 2D vectors. I rewrote that in a way that doesn't waste a third of the random values calculated (the CELL2 macro). Even with these optimizations the image with the noodles took an hour to render (200 samples on a hexacore machine). That might be a bit too much but for adding realism to a sweater (in my case that means with cat hairs all over it ;-) this might be a interesting opton.

Code availability

The code is available on GitHub. For ease of use I inlined the necessary functions from equations.h so the shader can be used as is, without external dependencies.

Sprinkles revisited

A while ago when I started this blog I created a simple shader for chocolate sprinkles and while I was reviewing comments I noticed that the code could do with some extra love.

(Whether you think this is a tasty image probably depends on you being a five year old or not :-)
There are a couple of issues with the original code:

Unnecessary functions: When the code was written OSL was lacking an implementation of the function to calculate the distance to a line segment. This is no longer the case and the functionality is provided by the three argument version of the distance() function.
Unsophisticated use of cellnoise: In the original code we added all sorts of vectors to get more than one uncorrelated cell noise value but in fact all noise variants in OSL support 4 dimensions so we can use a vector plus an additional value to get different noise values for a given cell.
Use of a deprecated function: The cellnoise() function is deprecated and we should use noise("cell", ...) in stead.
Unnecessary include: stdosl.h is automatically included by Blender already so there is no need for us to include it again.

And as a commenter mentioned, the code could do with an extra output to identify individual sprinkles so we can color them.

Code and node setup

The cleaned up code is pretty straight forward (and quite a bit shorter):

shader sprinkles(
 point Pos = P,
 float Scale = 1,
 int Np = 1,
 int Seed = 42,
 float Radius = 0.05,
 float Size = 1,
 output float Fac = 0,
 output float Random = 0
){
 point p = Pos * Scale;
 point f = floor(p);
 
 int xx,yy,np;
 vector one = 1;
 
 for( xx=-1; xx<=1; xx++){
  for( yy=-1; yy<=1; yy++){
   point ff = f + vector(xx,yy,0);
   float u=Seed;
 
   for( np=0; np < Np; np++){
    vector pd1 = 2*noise("cell",ff,u)-one;
    vector pd2 = 2*noise("cell",ff,u+1)-one;
    
    point p1 = ff + pd1;
    point p2 = ff + pd2;
    
    p2 = (p2 - p1)*Size+p1;
    
    // reduce to 2D 
    p1[2]=0;
    p2[2]=0;
    p [2]=0;
    
    float r = distance(p1,p2,p);
    if ( r < Radius ) {
      Fac = 1 - r/Radius;
      Random = noise("cell",ff,u+2);
    }
    u+=3;
   }
  }
 }
}

The node setup used to make the image at the top looks like this:

Code availability

The shader is available on GitHub. If you would like to know more about programming OSL you might be interested in my book "Open Shading Language for Blender". More on the availability of this book and a sample can be found on this page.

Some compositing experiments, mimicing black and white print

The other day I was reading about old printing techniques and I wanted to create the look and feel of screen printing with half tones and screen tones or even recreate a wood block printed look.

The approach I came up with and want to share here is pretty simple and versatile and completely implemented in compositor nodes and in the section below I show the necessary node setup. Although done in the compositor I dearly missed the possibility to create OSL nodes in the compositor because it would be nice to have some more textures to choose from.

Node setup

The general idea is to take a black and white version of your rendered image and combine it with some pattern. If that pattern is a smoothly varying grey scale image then it is possible compare that pattern against the b&w image and produce just black and white pixels but more of each will be produced at darker and lighter areas respectively. If that sounds confusing, bear with me it will get clearer in a minute.

The image at the start of the article was created with the node setup shown above. The image was renderered in Cycles but composited in Blender internal to be able to use textures in the compositor. This a bit awkward but not all that difficult: just set up your scene in cycles and if it renders ok, switch to bkender internal, add a new material to an arbitrary object (but don't assign it to any actual geometry) and then add a wood texture. This wood texture will give us paralel lines we can use in the compositor even though it is not actually used in any material. If you enable compositing in the post processing panel and press F12 to render, the scene will render in Cycles even though the render engine us set to internal.

In the noodle above we convert the image to black and white by plugging it into a value socket. This happens to be a colorramp node which allows use to tweak the contrast etc. The texture node refers to the wood texture we created earlier. This texture is also scaled so it will give us the line spacing we want. The real work is done in the two group nodes. These take the image and our texture, rotate the texture and do the comparison. Multiplying two versions with the lines at 90 degrees to each other give us the point pattern in the opening image.

The node group not only rotates the texture but compares it to the image. We do not simply compare whether the black & white value of the image is larger tgan that of the texture because this would give rise to horrible aliasing effects. Therefore the light orange nodes implement a smoothstep() whose smoothness is controlled by the red value node (More on this here).

Variations

Of course, in stead of multiplying two perpendicular sets if regular lines you could use any suitable pattern. The first image uses a single pattern of horizontal lines and the second image uses a marble pattern to create a wood block like appearance.

Even though you cannot uze OSL in the compositor I think it would be possible to setup two scenes in cycles, one with your image and one with just an OSL. pattern on a plane and combine these in the ckmpositor.