Small Blender Things: December 2013

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.