Artur

Romanesco broccoli pattern 2D shader

Rate this Entry

photo: Romanesco broccoli

I started playing with fractal shaders and in particularity I wanted to reproduce romanesco broccoli in SL. At first it was just a simple experiment but after each obstacle I got more persistent so eventually I discovered the world of calculus.
I'm going to skip the theory part about phyllotaxis, Fibonacci etc and get straight to the shader math. Good site that explains it all is: http://www.math.smith.edu/phyllo/

Most of the math was found on Wolfram and WikiPedia pages

First off we need to convert our UVs from cartesian to polar coordinate system


Code:
 r = sqrt( pow( x, 2 ) + pow( y, 2 ));

Code:
theta = atan( x, y) + PI;
// Note that PI offset has been added for polar translation 
// equation to work (for secondary florets)
Next, Golden Angle for distribution pattern.


Code:
golden = PI*(3-sqrt(5));
To draw a secondary florets or a circle on polar coordinates we can use polar translation equations


Code:
r2 = sqrt(r*r-2*r*r0*cos( theta-th0 )+r0*r0);

Code:
thera2 = atan( r*sin(theta)-r0*sin(th0), r*cos(theta)-r0*cos(th0));
Finally for an elliptical shape to compensate circle stretching as the cone height grows, following equation can be used. First in standard form



But we need r. A bit of rearrangement is needed


Code:
sqrt( r*r*(pow(sin(theta),2) + b*b*pow(cos(theta),2)) );
Final displacement shader
Code:
displacement brocdisp()
{
  float spread		= .01;		// spiral spreading factor
  float n	 	= 48; 		// number of florets
  float fS 		= .0015;		// floret scale
  float fE		= 1;		// floret scale exponent
  float b		= 1.7;		// floret ellipse eccentricity
  float rPxS		= 0.6;		// phyllotaxis scale within a floret
  uniform float octv	= 3;		// number of octaves

  float x = s - 0.5;
  float y = t - 0.5;
  
  //cartesian to polar coordinate system
  float r = sqrt( pow( x, 2 ) + pow( y, 2 ));	// pow for elliptical shape
  float theta = atan( x, y) + PI;		// offseting with PI so that secondary theta can be calculated

  float golden = PI*(3-sqrt(5));		// golden angle 2.39996 radians

  float disp = 0, rPx = rPxS, thPx = 0, fcell = 1, fcellP = 1, fR, i, p, f ,r0, th0, rF, thF, rFloret;
  float dispL[] = 0; resize( dispL, octv + 1 );

  P += -1 * (1-r) * normalize(N); // Initial cone shape displacement
  N = calculatenormal( P ); // Calculate new normals

  for(i=0; i<= octv; i+=1)
  {
  // Reset variables
  rPx = rPxS; thPx = 0; disp = 0; fcell = 1;
  // To mask individual floret cells, we need first pass of p to build disp phyllotaxis
  // and then it can be used to comp a mask: step(disp, fR - rFloret)
    for(p=0; p<= 1; p+=1)
    {
        // Florets on fermats spiral
	for(f=0; f<= n; f+=1)
	{
	  r0 =  f * spread;					// spreading radius. sqrt(f) for linear spread
	  th0 =  mod( golden * f, PI*2 );			// rotate golden angle
	  fR = sqrt( pow( f, fE ) * fS );			// floret Radius and its exponential growth. sqrt for linear shape

		rF = sqrt(r*r-2*r*r0*cos( theta-th0 )+r0*r0);				// polar r translation equation
		thF = atan( r*sin(theta)-r0*sin(th0), r*cos(theta)-r0*cos(th0));	// polar theta translation equation
 
		thF = mod( thF - th0, 2*PI );						// Rotate theta for ellipse orientation
		rFloret = sqrt( rF*rF*(pow(sin(thF),2) + b*b*pow(cos(thF),2)) );	// Elliptical floret for cone side
		rFloret = rFloret * (1-(r*0.1));					// Normalize florets
		rFloret = pow( rFloret, 0.8);						// compensate mushroom shape

	  if( p != 1 ) disp = max( disp, fR - rFloret );				// build Phyllotaxis displacement
	  if( p != 1 ) fcell = mix( fcell, fR, step(disp, fR - rFloret));

	  if( p == 1 ) rPx = mix( rPx, rPxS/fR * rFloret, step(disp, fR - rFloret)); 	// build phyllotaxis r for next iteration
	  if( p == 1 ) thPx = mix( thPx, thF + PI, step(disp, fR - rFloret));		// build phyllotaxis theta for next iteration
    	}
    }

    r = rPx; // store r for next iteration
    theta = thPx; // store theta for next iteration

    dispL[i] = disp * fcellP; // displacement layers mult by cell size from previous iteration
    fcellP *= fcell; // comp cell sizes for the next iteration

    P += -0.2 * dispL[i] * normalize(N);
    N = calculatenormal( P );
  }
}
Now remember, it's not a 3D representation but 2D pattern. Although I implemented it as a displacement shader I used very little displacement. Stronger displacement causes micropolygons stretch too much and it all gets just messy. But persistence is a bitch so I'm already thinking of 3D version as procedural primitives. This means 3D math, C and Python.

Broccoli displacement shader on a single polygon


RIB:
Code:
Option "searchpath" "shader" [ "&:.:${RIBDIR}" ] 
Display "ip" "it" "rgb"
Format 1024 1024 1
Projection  "perspective" "fov" [ 20 ]
Translate  0 -0.3 4
Rotate 50 1 0 0 Rotate 180 0 1 0

Translate  0 0 0
Rotate -90 1 0 0
Rotate 0   0 1 0
Scale 1 1 -1
WorldBegin
    LightSource "pointlight" 1 "intensity" 20.0 "lightcolor" [0.85 0.9 1] "from" [2 -2 -4]
    LightSource "pointlight" 1 "intensity" 30.0 "lightcolor" [1 1 0.8] "from" [-4 -8 2]
    Surface "plastic"
    Displacement "brocdisp"
    Attribute "bound" "displacement" [0.9]
    Attribute "shade" "relativeshadingrate" [0.331]
    AttributeBegin
        Scale 2 2 2
        Polygon "P" [-0.5 0 -0.5  -0.5 0 0.5  0.5 0 0.5  0.5 0 -0.5]
                "st" [0 0  0 1  1 1  1 0]
    AttributeEnd
WorldEnd

Submit "Romanesco broccoli pattern 2D shader" to Digg Submit "Romanesco broccoli pattern 2D shader" to del.icio.us Submit "Romanesco broccoli pattern 2D shader" to StumbleUpon Submit "Romanesco broccoli pattern 2D shader" to Google

Updated 08-26-2013 at 05:13 AM by Artur

Categories
Technical

Comments