Artur
Romanesco broccoli pattern 2D shader
by
, 08-31-2012 at 07:39 AM (120059 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 ));
Next, Golden Angle for distribution pattern.Code:theta = atan( x, y) + PI; // Note that PI offset has been added for polar translation // equation to work (for secondary florets)
To draw a secondary florets or a circle on polar coordinates we can use polar translation equationsCode:golden = PI*(3-sqrt(5));
Code:r2 = sqrt(r*r-2*r*r0*cos( theta-th0 )+r0*r0);
Finally for an elliptical shape to compensate circle stretching as the cone height grows, following equation can be used. First in standard formCode:thera2 = atan( r*sin(theta)-r0*sin(th0), r*cos(theta)-r0*cos(th0));
But we need r. A bit of rearrangement is needed
Final displacement shaderCode:sqrt( r*r*(pow(sin(theta),2) + b*b*pow(cos(theta),2)) );
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.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 ); } }
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