/***********************************************************/
/**                                                       **/
/**                                                       **/
/**    Bezier subdivision based 3D -> Image rendering     **/
/**                                                       **/
/**    (C) Theo Verelst theover@tiscali.nl                **/
/**                     http://theover.tripod.com         **/
/**                                                       **/
/**    Non-commercial use of the code is permitted        **/
/**    under the condition that a pointer to me is        **/
/**    included in every copy or derivative.              **/
/**                                                       **/
/**    Commercial rights reserved, mail me.               **/
/**                                                       **/
/**    This is a prototype I made in one night,           **/
/**    so don't expect it to be all to tested or          **/
/**    refined, it is a brute force approach I redid      **/
/**    from software I made in the time of my masters'    **/
/**    thesis, which started with reading an article      **/
/**    Pulleyblank and Kapenga about the feasibility      **/
/**    for such an approach as a on-chip algorithm.       **/
/**                                                       **/
/**    For this version 0.1, you need a file test.cb      **/
/**    with the Bezier Surface definition(s)              **/
/**                                                       **/
/**                                                       **/
/***********************************************************/

#include <stdio.h>
#include <math.h>
#include <fcntl.h>
#include <string.h>
#include <malloc.h>

#include <sys/types.h>

unsigned char *im;
int *zbuf;
int w,h;

#define MAXBEZ (256*256*4*4*3)
int b1[MAXBEZ], b2[MAXBEZ], b3[MAXBEZ];
int b1max=0, b2max=0, b3max=0;

int premul=500;
unsigned char rc = 255;

#define PI 3.1415926535

unsigned char * createim(w,h)
int w,h;
{
   unsigned char *d;
   int x,y;
   d = (unsigned char *) malloc(3*w*h);
   for (x=0; x<w; x++)
      for (y=0; y<h; y++) {
         d[(y*w+x)*3+0] = y/40;
         d[(y*w+x)*3+1] = 0;
         d[(y*w+x)*3+2] = x/40;
      }
   return d;
}

int * createzbuf(w,h)
int w,h;
{
   int *d;
   int x,y;
   d = (int *) malloc(w*h*sizeof(int));
   for (x=0; x<w; x++)
      for (y=0; y<h; y++) {
         d[(y*w+x)] = 32000;           /* for small integers in fact */
      }                                /* +z is into the screen      */
   return d;
}

int writeppm(n,w,h)
char n[];
int w,h;
{
   FILE *fp;
   fp = fopen(n,"wb");
   if ((void *) fp == NULL ) return(-1);

   fprintf(fp,"P6\n%d %d\n255\n",w,h);
   fwrite(im, 3, w*h,fp);
   fclose(fp);
   return(0);
}

#define point(x,y) 

void drawline(x1,y1,x2,y2)
int x1,x2,y1,y2;
{
   int i,j;

   for (i=x1; i<x2; i++)
      point(i,y1);
}

void rectangle(x1,y1,x2,y2,g,z)
int x1,x2,y1,y2;
unsigned char g;
{
   int i,j;

   for (i=x1; i<x2; i++)
      for (j=y1; j<y2; j++)
        if (z<zbuf[j*w+i]) {
         zbuf[j*w+i] = z;
         im[(j*w+i)*3+0] = g;
         im[(j*w+i)*3+1] = g;
         im[(j*w+i)*3+2] = g;
      }
   rc -= 2;
}

main(argc,argv)
int argc;
char *argv[];
{
   int i,j;
   int *res;
   w = 1024; h = 1024;


   printf("starting\n"); fflush(stdout); 
   b1max = read_rbc("test.cb");
   printf("ints read: %d\n", b1max);
   if (b1max < 0) { printf("Error: Unable to read input file, aborting.\n"); }
   printf("read input:\n");
   print_rbc2(b1, b1max);

   printf("Creating Image and Z-buffer\n"); fflush(stdout); 
   im = createim(w,h);
   zbuf = createzbuf(w,h);
   printf("Transforming\n"); fflush(stdout); 
   transform(b1,b1max);
   printf("Subdividing\n"); fflush(stdout); 
   subdivn(b1,b3,&b1max,&b3max,3,&res);  
   printf("Rendering/subdividing\n"); fflush(stdout); 
   for (j=0; j<b3max; j++) {
      for (i=0; i<48; i++)
         b1[i] = b3[48*j+i];
      b1max = 1;
      subdivn(b1,b2,&b1max,&b2max,5,&res);  
                                /* mind res indicates result is in b1 or b2 */
      for (i=0; i<b2max; i++)
         flatproject_rbc(b2+48*i,im,w,h);
   }
                                             /* Axises in the middle at z=0 */
   rectangle(0,h/2,w,h/2+1,100,0); 
   rectangle(w/2,0,w/2+1,h,100,0);
   /* rectangle(10,100,100,110,100,0); */
   printf("Writing Image\n");
   fflush(stdout); 
   writeppm("im.ppm", w,h);
}

int read_rbc(filename)
char *filename;
{
   int i,j;
   FILE *fp;

   b1max=0;
   i = 0;
   fp = fopen(filename,"r");
   if ((void *) fp == NULL ) return(-2);
   while (!feof(fp) && i<MAXBEZ) {
      if (fscanf(fp,"%d",&b1[i]) > 0) i++;
   }
   for (j=0; j<i; j++)          /* PREMULTIPLICATION */
      b1[j] *= premul;
   return i;
}

int print_rbc(b,m)
int *b,m;
{
   int i;

   for (i=0; i<m; i++) {
      if ( (i%3) == 0) printf("  ");
      printf(" %4d",b[i]);
      if ( (i%(3*4)) == (3*4-1)) printf("\n");
      if ( (i%(3*4*4)) == (3*4*4-1)) printf("\n");
   }
}

int print_rbc2(b,m)
int b[][4][4][3],m;
{
   int i;

   for (i=0; i<m; i++) {
      if ( (i%3) == 0) printf("  ");
      printf(" %4d",b[0][0][0][i]);
      if ( (i%(3*4)) == (3*4-1)) printf("\n");
      if ( (i%(3*4*4)) == (3*4*4-1)) printf("\n");
   }
}


int transform(b1,b1max)
int b1[][3], b1max;
{
   int i,j;
   int v[3];
   double w[3];
   double r[3][3], px,py,pz;

   /* Scale */
   w[0] = 2.3; w[1] = 2.3; w[2] = 2.3;
   for (i=0; i<b1max; i++)
      for (j=0; j<3; j++)
         b1[i][j] = (int) ((double) b1[i][j] * w[j]);
   /* Rotate */
   px=0.0; py=0.0; pz=1.0*PI/3.5;
   r[0][0] = cos(pz); r[1][0] = -sin(pz); r[2][0] = 0.0;
   r[0][1] = sin(pz); r[1][1] =  cos(pz); r[2][1] = 0.0;
   r[0][2] =     0.0; r[1][2] =      0.0; r[2][2] = 1.0;
   for (i=0; i<b1max; i++) {
      w[0] = (double) b1[i][0]; w[1] = (double) b1[i][1];
      for (j=0; j<2; j++)       /* only index 0 and 1 for rot around z axis */
         b1[i][j] = (int) (r[0][j]*w[0] + r[1][j]*w[1]);
   }
   /* Rotate */
   px=-PI/2.5; py=0.0; pz=0.0;
   r[0][0] =     1.0; r[1][0] =     0.0; r[2][0] =      0.0;
   r[0][1] =     0.0; r[1][1] = cos(px); r[2][1] = -sin(px);
   r[0][2] =     0.0; r[1][2] = sin(px); r[2][2] =  cos(px);
   for (i=0; i<b1max; i++) {
      w[1] = (double) b1[i][1]; w[2] = (double) b1[i][2];
      for (j=1; j<3; j++)       /* only index 1 and 2 for rot around x axis */
         b1[i][j] = (int) (r[1][j]*w[1] + r[2][j]*w[2]);
   }

   /* Translate */
   v[0] = 550*premul; v[1] = 300*premul; v[2] =   0*premul;
   for (i=0; i<b1max; i++)
      for (j=0; j<3; j++)
         b1[i][j] += v[j];
}

int xyzbbox(b,xmi,ymi,xma,yma,zmi,zma)
int *b, *xmi,*ymi,*xma,*yma, *zmi, *zma;
{
   int i;
   *xmi = *xma = b[0];
   *ymi = *yma = b[1];
   *zmi = *zma = b[2];
   for (i=3; i<4*4*3; i+=3) {
      if (b[i] < *xmi) {
         *xmi = b[i];
      } else {
         if (b[i] > *xma)
	    *xma = b[i];
      }
      if (b[i+1] < *ymi) {
         *ymi = b[i+1];
      } else {
         if (b[i+1] > *yma)
	    *yma = b[i+1];
      }
      if (b[i+2] < *zmi) {
         *zmi = b[i+2];
      } else {
         if (b[i+2] > *zma)
	    *zma = b[i+2];
      }
      
   }
}

int subdiv_rbc(in,out)
int in[][4][4][3], out[][4][4][3];
{
   int i,j,k;
   int t[2][4][4][3];
   int a , b , c;
   int   d , e  ;
   int     f    ;

   for (i=0; i<4; i++) {
      for (j=0; j<3; j++) {
         a = (in[0][i][0][j] + in[0][i][1][j])/2 ;
         b = (in[0][i][1][j] + in[0][i][2][j])/2 ;
         c = (in[0][i][2][j] + in[0][i][3][j])/2 ;
	 d = (a+b)/2;
	 e = (b+c)/2;
	 f = (d+e)/2;

	 t[0][i][0][j] = in[0][i][0][j];
	 t[0][i][1][j] = a;
	 t[0][i][2][j] = d;
	 t[0][i][3][j] = f;

	 t[1][i][0][j] = f;
	 t[1][i][1][j] = e;
	 t[1][i][2][j] = c;
	 t[1][i][3][j] = in[0][i][3][j];
      }
   }

   for (i=0; i<4; i++) {
      for (j=0; j<3; j++) {
         a = (t[0][0][i][j] + t[0][1][i][j])/2 ;
         b = (t[0][1][i][j] + t[0][2][i][j])/2 ;
         c = (t[0][2][i][j] + t[0][3][i][j])/2 ;
	 d = (a+b)/2;
	 e = (b+c)/2;
	 f = (d+e)/2;

	 out[0][0][i][j] = t[0][0][i][j];
	 out[0][1][i][j] = a;
	 out[0][2][i][j] = d;
	 out[0][3][i][j] = f;

	 out[1][0][i][j] = f;
	 out[1][1][i][j] = e;
	 out[1][2][i][j] = c;
	 out[1][3][i][j] = t[0][3][i][j];
      }
   }

   for (i=0; i<4; i++) {
      for (j=0; j<3; j++) {
         a = (t[1][0][i][j] + t[1][1][i][j])/2 ;
         b = (t[1][1][i][j] + t[1][2][i][j])/2 ;
         c = (t[1][2][i][j] + t[1][3][i][j])/2 ;
	 d = (a+b)/2;
	 e = (b+c)/2;
	 f = (d+e)/2;

	 out[2][0][i][j] = t[1][0][i][j];
	 out[2][1][i][j] = a;
	 out[2][2][i][j] = d;
	 out[2][3][i][j] = f;

	 out[3][0][i][j] = f;
	 out[3][1][i][j] = e;
	 out[3][2][i][j] = c;
	 out[3][3][i][j] = t[1][3][i][j];
      }
   }

/*   print_rbc2(out,(4*48));
   flatproject_rbc(out[0],im,w,h);
   flatproject_rbc(out[1],im,w,h);
   flatproject_rbc(out[2],im,w,h);
   flatproject_rbc(out[3],im,w,h);
*/
}

subdivn(b1,b2,b1m,b2m,n,res)
int b1[][48], b2[][48], *b1m, *b2m, n,*res;
{
   int i,j;

   for (i=0; i<n/2; i++) {
      for (j=0; j<*b1m; j++)
         subdiv_rbc(b1[j],b2[4*j]);
      *b2m = 4* *b1m;
      for (j=0; j<*b2m; j++)
         subdiv_rbc(b2[j],b1[4*j]);
      *b1m = 4* *b2m;
      
   }
   if ((n%2) == 1) {
      for (j=0; j<*b1m; j++)
         subdiv_rbc(b1[j],b2[4*j]);
      *b2m = 4* *b1m;
      *res = 1;
   } else {
      *res = 0;
   }

}

int flatproject_rbc(r,im,width,height)
int *r;
unsigned char *im;
int width, height;
{
   int x1,x2,y1,y2,z1,z2;
   int intens;
   int t[48];
   int i;

/*   for (i=0; i<48; i+=3) {
      t[i] = r[i]/100;
      t[i+1] = r[i+1]/100;
      t[i+2] = r[i+2];
   } */
   xyzbbox(r,&x1,&y1,&x2,&y2,&z1,&z2);
   intens = (400 * abs(z2 - z1)) / ((x2-x1+y2-y1)/1+1);
   if (intens > 250) intens = 250;
   intens = 250-intens;
   x1/=premul;
   x2/=premul;
   y1/=premul;
   y2/=premul;
   if (x1<0) { x1=0; } else if (x1>(width-1)) x1=width-1;
   if (x2<0) { x2=0; } else if (x2>(width-1)) x2=width-1;
   if (y1<0) { y1=0; } else if (y1>(height-1)) y1=height-1;
   if (y2<0) { y2=0; } else if (y2>(height-1)) y2=height-1;
   rectangle(x1,y1,x2,y2,(unsigned char) intens, z1/premul);
}