/**

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

*/

#declare Method = 0; // 0 = subdivision (red), 1 = double loop (blue)
#declare Depth = 3;

// Note that the triangles are not evenly distributed. The triangles at
// the upper right edge are larger than those at the two other edges.
// (The camera is right on, so it's not a trick of perspective, and the
// tangents for the three sides are the same, so that's not it either.)


// 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)>
		   pigment {rgb 1-<0.0,0.2,0.2>*redtriangle}
		}
	   #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)
		   pigment {rgb 1-<0.0,0.2,0.2>*redtriangle}
		}
	   #end		   
		   
	#else 
	
		#local uv12 = (uv1+uv2)/2;
		#local uv13 = (uv1+uv3)/2;
		#local uv23 = (uv2+uv3)/2;
		#local redtriangletemp = redtriangle;
		#local redtriangle = 1-redtriangletemp;
		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)
		#local redtriangle = redtriangletemp;
		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 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

   #if (Method=0)
		//mesh{
			#local redtriangle = mod(depth+1,2);
			bezier3triangles_rec(<0,0>,<0,1>,<1,0>,fx,fy,fz,smoothing,depth)
		//}
   #else
      #declare Points = pow(2,depth);
      #declare C = 0;
      #while (C<Points)
         #declare D = 0;
         #while (D<Points)
            #if (C+D<Points)
               #declare uv1 = <C  ,D  >/Points;
               #declare uv2 = <C+1,D  >/Points;
               #declare uv3 = <C  ,D+1>/Points;
               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)
                  pigment {rgb <0.8,0.8,1.0>}
               }
               #if (C+D<Points-1)
                  #declare uv1 = <C+1,D+1>/Points;
                  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)
                     pigment {rgb 1.0}
                  }
               #end
            #end
            #declare D = D+1;
         #end
         #declare C = C+1;
      #end
   #end
}		

#end


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

// SAMPLE SCENE

// Render as animation with these command line settings:
// +kc +kff50


#include "colors.inc"

// Length of tangent vectors
#declare Pda = 1.5;
#declare Pdb = Pda;
#declare Pdc = Pda;
#declare Pdd = Pda;
#declare Pde = Pda;
#declare Pdf = Pda;

#local rad=0.1;

#default {finish {phong 0.3}}
#declare Seed = seed(0);

union {
   
   object{ // rounded corner
      bezier3triangles(
         <2,0,2>,<2,0,2>-Pda*z,<2,2,0>-Pdb*y,
         <2,2,0>,<2,2,0>-Pdc*x,<0,2,2>-Pdd*z,
         <0,2,2>,<0,2,2>-Pde*y,<2,0,2>-Pdf*x,
         <2,0,2>-Pdf*x-Pda*z,<2,2,0>-Pdb*y-Pdc*x,<0,2,2>-Pdd*z-Pde*y, // inner points
      Depth,1,0) // depth, smoothing, enable visible control points
   }
   
   bicubic_patch {
      type 1 u_steps Depth v_steps Depth
      <2,0,2>,<2,0,2>-Pda*z,<2,2,0>-Pdb*y,<2,2,0>,
      <3,0,2>,<3,0,2>-Pda*z,<3,2,0>-Pdb*y,<3,2,0>,
      <4,0,2>,<4,0,2>-Pda*z,<4,2,0>-Pdb*y,<4,2,0>,
      <5,0,2>,<5,0,2>-Pda*z,<5,2,0>-Pdb*y,<5,2,0>
   }
   bicubic_patch {
      type 1 u_steps Depth v_steps Depth
      <2,0,2>,<2,0,2>-Pdc*z,<2,2,0>-Pdd*y,<2,2,0>,
      <3,0,2>,<3,0,2>-Pdc*z,<3,2,0>-Pdd*y,<3,2,0>,
      <4,0,2>,<4,0,2>-Pdc*z,<4,2,0>-Pdd*y,<4,2,0>,
      <5,0,2>,<5,0,2>-Pdc*z,<5,2,0>-Pdd*y,<5,2,0>
      matrix <0,1,0,0,0,1,1,0,0,0,0,0>
   }
   bicubic_patch {
      type 1 u_steps Depth v_steps Depth
      <2,0,2>,<2,0,2>-Pde*z,<2,2,0>-Pdf*y,<2,2,0>,
      <3,0,2>,<3,0,2>-Pde*z,<3,2,0>-Pdf*y,<3,2,0>,
      <4,0,2>,<4,0,2>-Pde*z,<4,2,0>-Pdf*y,<4,2,0>,
      <5,0,2>,<5,0,2>-Pde*z,<5,2,0>-Pdf*y,<5,2,0>
      matrix <0,0,1,1,0,0,0,1,0,0,0,0>
   }
   
   #declare Patch2 =
   bicubic_patch {
      type 1 u_steps Depth v_steps Depth
      <2,0,2>,<2,0,3>,<2,0,4>,<2,0,5>,
      <3,0,2>,<3,0,3>,<3,0,4>,<3,0,5>,
      <4,0,2>,<4,0,3>,<4,0,4>,<4,0,5>,
      <5,0,2>,<5,0,3>,<5,0,4>,<5,0,5>
   }
   object {Patch2}
   object {Patch2 matrix <0,1,0,0,0,1,1,0,0,0,0,0>}
   object {Patch2 matrix <0,0,1,1,0,0,0,1,0,0,0,0>}
   
   pigment {rgb 1}
   scale <1,-1,1>
}

camera {
   location <-3,3,-3>
   look_at <3,-3,3>
}

light_source {<+1,2,-3>*1000, color <1.0,1.0,0.8>}
light_source {<-3,2,+1>*1000, color <0.8,0.8,1.0>}

background {rgb 1}