/**

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
@version 0.11
created 21 August 2002

*/

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

// 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, P1,P2,P3, N1,N2,N3, fx, fy, fz, smoothing, depth) 
   #if (depth = 0)
      #if (smoothing != 1)
         triangle {P1,P2,P3}
      #else
         smooth_triangle {P1,N1,P2,N2,P3,N3}
      #end
   #else
      #local uv12 = (uv1+uv2)/2;
      #local uv13 = (uv1+uv3)/2;
      #local uv23 = (uv2+uv3)/2;
      #local p12 = <fx(uv12.x,uv12.y),fy(uv12.x,uv12.y),fz(uv12.x,uv12.y)>;
      #local p13 = <fx(uv13.x,uv13.y),fy(uv13.x,uv13.y),fz(uv13.x,uv13.y)>;
      #local p23 = <fx(uv23.x,uv23.y),fy(uv23.x,uv23.y),fz(uv23.x,uv23.y)>;
      #local n12 = pnormal(fx,fy,fz,uv12);
      #local n13 = pnormal(fx,fy,fz,uv13);
      #local n23 = pnormal(fx,fy,fz,uv23);
      #local redtriangletemp = redtriangle;
      #local redtriangle = 1-redtriangletemp;
      bezier3triangles_rec(uv1,uv12,uv13,P1,p12,p13,N1,n12,n13,fx,fy,fz,smoothing,depth-1)
      bezier3triangles_rec(uv2,uv12,uv23,P2,p12,p23,N2,n12,n23,fx,fy,fz,smoothing,depth-1)
      bezier3triangles_rec(uv3,uv23,uv13,P3,p23,p13,N3,n23,n13,fx,fy,fz,smoothing,depth-1)
      #local redtriangle = redtriangletemp;
      bezier3triangles_rec(uv12,uv13,uv23,p12,p13,p23,n12,n13,n23,fx,fy,fz,smoothing,depth-1)
   #end
#end

#macro bezier3triangles(A,B,C,D,E,F,G,H,I,J,K,L,depth,smoothing,rad)
   
   #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 fx2 = function(u,v,w) {
      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*u*w*w*(1-v*(1+u))*r201x +
      
      9*w*v*u*u*r112x +
      9*w*v*v*u*r121x +
      9*w*w*v*u*r211x
      
      -(
      pow(w,3)*(1-3*v*u) +
      pow(v,3)*(1-3*w*u) +
      pow(u,3)*(1-3*w*v) +
      
      3*w*w*v*(1-u*(1+v)) +
      3*w*v*v*(1-u*(1+w)) +
      3*v*v*u*(1-w*(1+u)) +
      3*v*u*u*(1-w*(1+v)) +
      3*u*u*w*(1-v*(1+w)) +
      3*u*w*w*(1-v*(1+u)) +
      
      9*w*v*u*u +
      9*w*v*v*u +
      9*w*w*v*u - 1
      ) *
      (r112x+r121x+r211x)/3
   }
   
   #local fy2 =  function(u,v,w) {
      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*u*w*w*(1-v*(1+u))*r201y +
      
      9*w*v*u*u*r112y +
      9*w*v*v*u*r121y +
      9*w*w*v*u*r211y
      
      -(
      pow(w,3)*(1-3*v*u) +
      pow(v,3)*(1-3*w*u) +
      pow(u,3)*(1-3*w*v) +
      
      3*w*w*v*(1-u*(1+v)) +
      3*w*v*v*(1-u*(1+w)) +
      3*v*v*u*(1-w*(1+u)) +
      3*v*u*u*(1-w*(1+v)) +
      3*u*u*w*(1-v*(1+w)) +
      3*u*w*w*(1-v*(1+u)) +
      
      9*w*v*u*u +
      9*w*v*v*u +
      9*w*w*v*u - 1
      ) *
      (r112y+r121y+r211y)/3
   }
   
   #local fz2 =  function(u,v,w) {
      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*u*w*w*(1-v*(1+u))*r201z +
      
      9*w*v*u*u*r112z +
      9*w*v*v*u*r121z +
      9*w*w*v*u*r211z
      
      -(
      pow(w,3)*(1-3*v*u) +
      pow(v,3)*(1-3*w*u) +
      pow(u,3)*(1-3*w*v) +
      
      3*w*w*v*(1-u*(1+v)) +
      3*w*v*v*(1-u*(1+w)) +
      3*v*v*u*(1-w*(1+u)) +
      3*v*u*u*(1-w*(1+v)) +
      3*u*u*w*(1-v*(1+w)) +
      3*u*w*w*(1-v*(1+u)) +
      
      9*w*v*u*u +
      9*w*v*v*u +
      9*w*w*v*u - 1
      ) *
      (r112z+r121z+r211z)/3
   }
   
   #macro fw()
      (1-fu-fv)/max(1-2*fu*fv,1)
   #end
   
   #declare Skew = function(v0,v1,v2) {
      v0*2/3+(1-v1-v2)/max(1-2*v1*v2,epsilon1)/3
   }
   
   #local fx = function(fu,fv) { fx2(Skew(fu,fv,fw()),Skew(fv,fw(),fu),Skew(fw(),fu,fv)) }
   #local fy = function(fu,fv) { fy2(Skew(fu,fv,fw()),Skew(fv,fw(),fu),Skew(fw(),fu,fv)) }
   #local fz = function(fu,fv) { fz2(Skew(fu,fv,fw()),Skew(fv,fw(),fu),Skew(fw(),fu,fv)) }
   
   union{
      #if (rad)
         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);
            #local p1 = <fx(0,0),fy(0,0),fz(0,0)>;
            #local p2 = <fx(0,1),fy(0,1),fz(0,1)>;
            #local p3 = <fx(1,0),fy(1,0),fz(1,0)>;
            #local n1 = pnormal(fx,fy,fz,<0,0>);
            #local n2 = pnormal(fx,fy,fz,<0,1>);
            #local n3 = pnormal(fx,fy,fz,<1,0>);				
            bezier3triangles_rec(<0,0>,<0,1>,<1,0>,p1,p2,p3,n1,n2,n3,fx,fy,fz,smoothing,depth)
         }
      #else
         #declare Points = pow(2,depth);
         #declare Vectors = array[Points+1][Points+1]
         #declare Normals = array[Points+1][Points+1]
         #declare C = 0;
         #debug "Creating arrays...\n"
         #while (C<=Points)
            #declare D = 0;
            #while (D<=Points)
               #if (C+D<=Points)
                  #declare Vectors[C][D] = <fx(C/Points,D/Points),fy(C/Points,D/Points),fz(C/Points,D/Points)>;
                  #declare Normals[C][D] = pnormal(fx,fy,fz,<C,D>/Points);
               #end
               #declare D = D+1;
            #end
            #declare C = C+1;
         #end
         #debug "Creating vertex_vectors...\n"
         mesh2 {
         vertex_vectors {
            (Points+1)*(Points+1)/2+(Points+1)/2
            #declare C = 0;
            #while (C<=Points)
               #declare D = 0;
               #while (D<=Points)
                  #if (C+D<=Points)
                     Vectors[C][D]
                  #end
                  #declare D = D+1;
               #end
               #declare C = C+1;
            #end
         }
         #if (smoothing)
         #debug "Creating normal_vectors...\n"
         normal_vectors {
            (Points+1)*(Points+1)/2+(Points+1)/2
            #declare C = 0;
            #while (C<=Points)
               #declare D = 0;
               #while (D<=Points)
                  #if (C+D<=Points)
                     Normals[C][D]
                  #end
                  #declare D = D+1;
               #end
               #declare C = C+1;
            #end
         }
         #end
         #debug "Creating face_indices...\n"
         face_indices {
            Points*Points
            #declare C = 0;
            #declare Nr = 0;
            #while (C<Points)
               #declare RowA = Nr;
               #declare RowB = Nr+Points-C+1;
               #declare D = 0;
               #while (D<Points)
                  #if (C+D<Points)
                     < D   +RowA,(D+1)+RowA, D   +RowB>
                     #if (C+D<Points-1)
                     <(D+1)+RowA, D   +RowB,(D+1)+RowB>
                     #end
                  #end
                  #declare D = D+1;
               #end
               #declare Nr = Nr+Points-C+1;
               #declare C = C+1;
            #end
         }
         }
         #debug "Finished!\n"
      #end
      
   }
   
#end


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

// SAMPLE SCENE

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


#include "colors.inc"

// Length of tangent vectors
#declare Pda = 2-1.5*sin((clock+0/6)*2*pi);
#declare Pdb = 2-1.5*sin((clock+2/6)*2*pi);
#declare Pdc = 2-1.5*sin((clock+2/6)*2*pi);
#declare Pdd = 2-1.5*sin((clock+4/6)*2*pi);
#declare Pde = 2-1.5*sin((clock+4/6)*2*pi);
#declare Pdf = 2-1.5*sin((clock+0/6)*2*pi);

#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,on,0.0) // depth, smoothing, radius of control points (0 = not enabled)
   }
   
   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}