/* draw_rgb.c 
  *
  *  Jiwon Hahn
  *  
  *  This program reads a material file, calculate XYZ and rgb values for each pixel, 
  *  and create corresponding r,g,b image files.
  *  
  *  Calling convention: prog materal_file 
  *  Output files are out.b out.g out.r 
  */


#include <stdio.h>
#include <stdlib.h>
#include <fcntl.h>
#include <math.h>
//choose the parameters 
#define M 0.5
#define Lx 0.578
#define Ly 0.578
#define Lz 0.578 
#define Z0 0 
#define R 100 
#define C 1
#define I_SIZE 250 

//Calculation of Fresnel Reflectance
float fresnel(float n,float k, int th)
{
    float a, b, Rn, Rp;
    a=sqrt(0.5)*sqrt(sqrt((n*n-k*k-sin(th)*sin(th))*(n*n-k*k-sin(th)*sin(th))+4*n*n*k*k)+(n*n-k*k-sin(th)*sin(th)));
    b=sqrt(0.5)*sqrt(sqrt((n*n-k*k-sin(th)*sin(th))*(n*n-k*k-sin(th)*sin(th))+4*n*n*k*k)-(n*n-k*k-sin(th)*sin(th)));
    Rn=(a*a+b*b-2*a*cos(th)+cos(th)*cos(th))/(a*a+b*b+2*a*cos(th)+cos(th)*cos(th)); 
    Rp=Rn*(a*a+b*b-2*a*sin(th)*tan(th)+sin(th)*sin(th)*tan(th)*tan(th)) / (a*a+b*b+2*a*sin(th)*tan(th)+sin(th)*sin(th)*tan(th)*tan(th));
    return 0.5*(Rn+Rp);
}

//Interpolation - makes a continuous function and returns the value corresponding to given lambda
float getvalue(float l, float *x, float *f, int N)
{
    int i;
    for(i=0; i<N-1; i++)
        if(x[i]<=l && l<=x[i+1])
            return ((f[i+1]-f[i])/(x[i+1]-x[i])*(l-x[i])+f[i]);
    return 0.0;  // in case the given lambda is out of range
}

void print()
{	
    printf("******The parameters used in this Program******\n");
    printf("      Light source vector=(%.2f,%.2f,%.2f)\n",Lx,Ly,Lz);
    printf("      Surface roughness: m=%.1f\n",M);
}
    
main(int argc, char **argv)
{
   int i, j, Nxyz, NF; //number of xyz and fresnel values
   float *l, *x, *y, *z, *La, *F; //arrays for storing lambda,x,y,z and fresnel values
   float n,k; //temporary storages for scanning material file and calculating fresnel values
   float X=0.0,Y=0.0,Z=0.0;  //XYZ values
   float r,g,b;
   float max=0.0;
   float Ti[3][3]={{3.194,-1.259,0.260},{-1.848,2.230,-0.193},{-0.344,0.044,0.937}};
   FILE *fp, *fxyz;
   float alpha, D;
   float N[3], L[3], V[3]={0,0,1}, H[3];
   float NH, NV, VH, NL, LH;
   float G, Gb, Gc;
   float c; //counter
   int thi; //theta i
   int cx, cy, cz; //x,y,z coordinate
   float BRDFco=0.0; //brdf coefficients, which is brdf divided by fresnel reflectance 
   float rim[I_SIZE][I_SIZE];
   float gim[I_SIZE][I_SIZE];
   float bim[I_SIZE][I_SIZE];
   unsigned char rimage[I_SIZE][I_SIZE];
   unsigned char gimage[I_SIZE][I_SIZE];
   unsigned char bimage[I_SIZE][I_SIZE];
   FILE *fr, *fg, *fb; //file descriptors

    print();
    L[0]=Lx; L[1]=Ly; L[2]=Lz; //given Light source direction
    //V[0]=0; V[1]=0; V[2]=1;   viewer is always (0,0,1)
    H[0]=Lx/sqrt(Lx*Lx+Ly*Ly+(Lz+1)*(Lz+1)); H[1]=Ly/sqrt(Lx*Lx+Ly*Ly+(Lz+1)*(Lz+1)); H[2]=(Lz+1)/sqrt(Lx*Lx+Ly*Ly+(Lz+1)*(Lz+1));  //H bisects the angle between viewer and light
    VH=H[2];

//Read xyz file
   if((fxyz=fopen("xyz","r"))==NULL){
      fprintf(stderr, "ERROR: could not open the xyz file");
      exit(1);
    }
    fscanf(fxyz,"%d\n",&Nxyz); //the number of xyz
    l=(float *)malloc(sizeof(float)*Nxyz);
    x=(float *)malloc(sizeof(float)*Nxyz);
    y=(float *)malloc(sizeof(float)*Nxyz);
    z=(float *)malloc(sizeof(float)*Nxyz);
    for(i=0; i<Nxyz; i++)
        fscanf(fxyz,"%f%f%f%f\n",l+i,x+i,y+i,z+i); //store lambda, x,y,z in arrays.
    fclose(fxyz);
    
//Read material file
    if((fp=fopen(argv[1],"r"))==NULL){
       fprintf(stderr, "ERROR: could not open the output file");
       exit(1);
    }
    fscanf(fp,"%d\n",&NF);
    La=(float *)malloc(sizeof(float)*NF);
    F=(float *)malloc(sizeof(float)*NF);
    LH=L[0]*H[0]+L[1]*H[1]+L[2]*H[2];
    thi=acos(LH); 
    for(i=0;i<NF;i++){
        fscanf(fp,"%f%f%f\n",La+i,&n,&k); //store lambda in array La
        *(F+i)=fresnel(n,k,thi); //store fresnel values in F     
    }
    fclose(fp);

// This is the main part of the program
// Rough homogeneous surface model is used to generate the pixel values
    for(i=0; i<I_SIZE; i++){
        for(j=0; j<I_SIZE; j++){
            cx=j-I_SIZE/2;
            cy=i-I_SIZE/2;
            cz=Z0+sqrt(R*R-(cx*cx+cy*cy));

           if(cx*cx+cy*cy<R*R){ //inside the circle
            N[0]=cx/sqrt(cx*cx+cy*cy+cz*cz); N[1]=cy/sqrt(cx*cx+cy*cy+cz*cz); N[2]=cz/sqrt(cx*cx+cy*cy+cz*cz);
            NH=N[0]*H[0]+N[1]*H[1]+N[2]*H[2];
            NV=N[2];
            NL=N[0]*L[0]+N[1]*L[1]+N[2]*L[2];
            alpha=acos(NH);
            D=exp(-alpha*alpha/(M*M));
            //calculating G
            Gb=2*NH*NV/VH; Gc=2*NH*NL/VH;
            if(Gb>=1 && Gc>=1) G=1;
            else if(Gb<Gc) G=Gb;
            else G=Gc;
            if(4*NL*NV!=0)   BRDFco=D*G/(4*NL*NV);

            //Integration for getting X, Y, Z
            for(c=380.0; c<=780.0; c+=1.0){
                X+=getvalue(c,l,x,Nxyz)*getvalue(c,La,F,NF);
                Y+=getvalue(c,l,y,Nxyz)*getvalue(c,La,F,NF);
                Z+=getvalue(c,l,z,Nxyz)*getvalue(c,La,F,NF);
            }
            X=X*BRDFco;
            Y=Y*BRDFco;
            Z=Z*BRDFco;
		printf("i=%d,j=%d, X=%.3f,Y=%.3f,Z=%.3f\n",i,j,X,Y,Z);
            //Transformation to r,g,b values
            r=Ti[0][0]*X+Ti[1][0]*Y+Ti[2][0]*Z;
            g=Ti[0][1]*X+Ti[1][1]*Y+Ti[2][1]*Z;
            b=Ti[0][2]*X+Ti[1][2]*Y+Ti[2][2]*Z;
	   
	    //max rgb determinatioN
	    if(max<r) max=r;
	    if(max<g) max=g;
	    if(max<b) max=b;
	  
	    //store the r,g,b values in arrays
	    rim[i][j]=r;
            gim[i][j]=g;
            bim[i][j]=b;        
          }    
        }
     }
     printf("\nmaximum rgb value is %.2f\n",max);
    //free dynamic array allocations
    free(l); free(x); free(y); free(z);
    free(La); free(F);

    //scale the r,g,b values to maximum 255.
    for(i=0; i<I_SIZE; i++)
	for(j=0;j<I_SIZE;j++){
	  rimage[i][j]=rim[i][j]*255/max;
          bimage[i][j]=bim[i][j]*255/max;
          gimage[i][j]=gim[i][j]*255/max;
	}
   // generate raw image files 
    fr = fopen("out.r", "w");
    fg = fopen("out.g", "w");
    fb = fopen("out.b", "w");
    for(i=0;i<I_SIZE;i++)
	for(j=0; j<I_SIZE; j++)
           fputc(rimage[i][j],fr);
    for(i=0;i<I_SIZE;i++)
	for(j=0; j<I_SIZE; j++)
           fputc(gimage[i][j],fg);
    for(i=0;i<I_SIZE;i++)
	for(j=0; j<I_SIZE; j++)
           fputc(bimage[i][j],fb);
    fclose(fr); fclose(fg); fclose(fb);
}