Artur

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