Artur

Romanesco broccoli pattern 2D shader by , 08-31-2012 at 06:39 AM (118604 Views) 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)) );
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)
{
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]