/**

Macro for creating mesh with 3-sided bezier patches of degree 3.
Call with 

  bezier3triangles(A,B,C,D,E,F,G,H,I,J,K,L,depth,smooth,visible_points)
  
where A..L are control points:
- A..I points on the border in cyclic order. A,D,G are edges, B,C, 
  between A and D and so on.
- J,K,L are inner points, J closest to A, K closest to D, L closest 
  to G
  
and 
- depth: tesselation depth
- smooth: if set to 1 smooth_triangles are generated
- visible_points: if set to 1 control points are made visible

@author Micha Riser
created 21 August 2002

*/


// some small values bigger than zero
#declare epsilon1 = 0.001;
#declare epsilon2 = 0.001;

#macro pnormal(fx,fy,fz,uv)
	#local uu = uv.x; #local vv=uv.y;
	#local ru = -<fx(uu,vv),fy(uu,vv),fz(uu,vv)>;
	#local rv = ru;
	#local uu = uu+epsilon2;
	#local ru = ru + <fx(uu,vv),fy(uu,vv),fz(uu,vv)>;
	#local uu = uv.x;
	#local vv = vv+epsilon2;
	#local rv = rv + <fx(uu,vv),fy(uu,vv),fz(uu,vv)>;
	vcross(ru,rv)
#end

#macro bezier3triangles_rec(uv1, uv2, uv3, fx, fy, fz, smoothing, depth) 

	#if (depth =0)
	
 	   #if (smoothing != 1)
		triangle{
		   <fx(uv1.x,uv1.y),fy(uv1.x,uv1.y),fz(uv1.x,uv1.y)>,
		   <fx(uv2.x,uv2.y),fy(uv2.x,uv2.y),fz(uv2.x,uv2.y)>,
		   <fx(uv3.x,uv3.y),fy(uv3.x,uv3.y),fz(uv3.x,uv3.y)>}		 
	   #else	   
		   
		smooth_triangle{
		   <fx(uv1.x,uv1.y),fy(uv1.x,uv1.y),fz(uv1.x,uv1.y)>,pnormal(fx,fy,fz,uv1),
		   <fx(uv2.x,uv2.y),fy(uv2.x,uv2.y),fz(uv2.x,uv2.y)>,pnormal(fx,fy,fz,uv2),
		   <fx(uv3.x,uv3.y),fy(uv3.x,uv3.y),fz(uv3.x,uv3.y)>,pnormal(fx,fy,fz,uv3)}		 		   
	   #end		   
		   
	#else 
	
		#local uv12 = (uv1+uv2)/2;
		#local uv13 = (uv1+uv3)/2;
		#local uv23 = (uv2+uv3)/2;
		bezier3triangles_rec(uv1,uv12,uv13,fx,fy,fz,smoothing,depth-1)
		bezier3triangles_rec(uv2,uv12,uv23,fx,fy,fz,smoothing,depth-1)
		bezier3triangles_rec(uv3,uv23,uv13,fx,fy,fz,smoothing,depth-1)
		bezier3triangles_rec(uv12,uv13,uv23,fx,fy,fz,smoothing,depth-1)
	
	#end

#end

#macro w()
	(1-u-v)/max( 1-2*u*v,epsilon1 ) 
#end

#macro bezier3triangles(A,B,C,D,E,F,G,H,I,J,K,L,depth,smoothing,visible_points)

	#local r003x=A.x; #local r003y=A.y; #local r003z=A.z;
	#local r012x=B.x; #local r012y=B.y; #local r012z=B.z;
	#local r021x=C.x; #local r021y=C.y; #local r021z=C.z;
	#local r030x=D.x; #local r030y=D.y; #local r030z=D.z;
	#local r120x=E.x; #local r120y=E.y; #local r120z=E.z;
	#local r210x=F.x; #local r210y=F.y; #local r210z=F.z;
	#local r300x=G.x; #local r300y=G.y; #local r300z=G.z;
	#local r201x=H.x; #local r201y=H.y; #local r201z=H.z;	
	#local r102x=I.x; #local r102y=I.y; #local r102z=I.z;
	#local r112x=J.x; #local r112y=J.y; #local r112z=J.z;	
	#local r121x=K.x; #local r121y=K.y; #local r121z=K.z;	
	#local r211x=L.x; #local r211y=L.y; #local r211z=L.z;		


#local rad=0.05;


#local fx =  function(u,v){   		  	
			pow(w(),3)*(1-3*v*u)*r300x +
			pow(v,3)*(1-3*w()*u)*r030x +
			pow(u,3)*(1-3*w()*v)*r003x +
			
			3*w()*w()*v*(1-u*(1+v))*r210x +
			3*w()*v*v*(1-u*(1+w()))*r120x +
			3*v*v*u*(1-w()*(1+u))*r021x +
			3*v*u*u*(1-w()*(1+v))*r012x +
			3*u*u*w()*(1-v*(1+w()))*r102x +
			3*w()*w()*u*(1-v*(1+u))*r201x +
			
			9*w()*v*u*u*r112x +
			9*w()*v*v*u*r121x +
			9*w()*w()*v*u*r211x
		  }
#local fy =  function(u,v){   		  	
			pow(w(),3)*(1-3*v*u)*r300y +
			pow(v,3)*(1-3*w()*u)*r030y +
			pow(u,3)*(1-3*w()*v)*r003y +
			
			3*w()*w()*v*(1-u*(1+v))*r210y +
			3*w()*v*v*(1-u*(1+w()))*r120y +
			3*v*v*u*(1-w()*(1+u))*r021y +
			3*v*u*u*(1-w()*(1+v))*r012y +
			3*u*u*w()*(1-v*(1+w()))*r102y +
			3*w()*w()*u*(1-v*(1+u))*r201y +
			
			9*w()*v*u*u*r112y +
			9*w()*v*v*u*r121y +
			9*w()*w()*v*u*r211y

		  }
#local fz =  function(u,v){   		  	
			pow(w(),3)*(1-3*v*u)*r300z +
			pow(v,3)*(1-3*w()*u)*r030z +
			pow(u,3)*(1-3*w()*v)*r003z +
			
			3*w()*w()*v*(1-u*(1+v))*r210z +
			3*w()*v*v*(1-u*(1+w()))*r120z +
			3*v*v*u*(1-w()*(1+u))*r021z +
			3*v*u*u*(1-w()*(1+v))*r012z +
			3*u*u*w()*(1-v*(1+w()))*r102z +
			3*w()*w()*u*(1-v*(1+u))*r201z +
			
			9*w()*v*u*u*r112z +
			9*w()*v*v*u*r121z +
			9*w()*w()*v*u*r211z

		  }

union{
#if (visible_points)
		sphere{A,rad pigment{Red}}
		sphere{D,rad pigment{Red}}
		sphere{G,rad pigment{Red}}
		sphere{B,rad pigment{Green}}
		sphere{C,rad pigment{Green}}
		sphere{E,rad pigment{Green}}
		sphere{F,rad pigment{Green}}		
		sphere{H,rad pigment{Green}}
		sphere{I,rad pigment{Green}}				
		sphere{J,rad pigment{Yellow}}
		sphere{K,rad pigment{Yellow}}		
		sphere{L,rad pigment{Yellow}}		
#end


		mesh{
			bezier3triangles_rec(<0,0>,<0,1>,<1,0>,fx,fy,fz,smoothing,depth)
		}
}		

#end


/////////////////////////////////////////////////////////////////////////////////////////

/* example scene */

#include "colors.inc"

object{ bezier3triangles( 
                 <0,0,0>,<0.3,0.3,-.2>,<0.6,0.3,0>,  // under line // first point is edge
                 <1,0,0>,<1,0.3,0.4>,<0.8,0.4,0.8>,  // right line // first point is edge
		 <0.6,0,1>,<0.4,.4,0.8>,<0.1,.1,0.5>, // left line // first point is edge
		 <0.3,1,0.3>,<0.7,.9,0.3>,<0.6,.8,0.6> // inner points
		 ,4,1,1) // depth, smoothing, enable visible control points
	pigment{rgb 1 }
	}


	
light_source{<0,5,0>,1}

camera{location <0.5,4,-4>
look_at 0.5
angle 30
}

plane{y,-0.01 pigment{rgb <0,0,1>}}