//**************************************************************************
//
//      CSE 486 Project #4
//      By Anirudh Modi (anirudh@bart.aero.psu.edu) on 4/18/99-Sun
//	 & Jim Geis (geis@cse.psu.edu)
//	 & Howard Elikan (elikan@cse.psu.edu)
//
//**************************************************************************
// 	Have implemented the following ->
//		1. Image-Irradiance Equation
//**************************************************************************
#include <iostream.h>
#include <fstream.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include "image.h"
#include "matrix.h"
//**************************************************************************
#define FOCAL_LENGTH (5.0)
#define SPHERE_RADIUS (4.0)
#define BG_INTENSITY 64
#define IMAGE_DIM 256
#define XY_RANGE (1.0)
#define N_SAMPLE 512
//**************************************************************************
double ReflectanceSphere (double theta, double phi, double p_s, double q_s);
double Timer (void);
//**************************************************************************
int main (int argc, char *argv[])
{
if (argc != 4) 
	{
	 cerr << "Syntax: " << argv[0] << " <p_s> <q_s> <z0>" << endl;
	 exit(-1);
	}

double p_s = atof(argv[1]);
double q_s = atof(argv[2]);
double z0 = atof(argv[3]);

double r = SPHERE_RADIUS;
double lambda = FOCAL_LENGTH;
double theta_i = asin(r/(z0-lambda));
double ri = lambda*r*cos(theta_i)/(z0-lambda-r*sin(theta_i));

if (ri > XY_RANGE)
	{
	 double theta_i2 = acos(XY_RANGE/ri);
	 if (theta_i2 > theta_i)
		theta_i = theta_i2;
	}

Matrix T = TransformToImagePlane(0,0,0,0,0,1,0,0,0,lambda);

cout << "---------------------------------------" 
	<< "---------------------------------------" << endl;
cout << "Radius of sphere = " << r << endl;
cout << "Focal length of camera lens = " << lambda << endl;
cout << "Range of image plane = [" << -XY_RANGE << "," << XY_RANGE << "]x["
	<< -XY_RANGE << "," << XY_RANGE << "]" << endl;
cout << "---------------------------------------" 
	<< "---------------------------------------" << endl;
cout << "p_s = " << p_s << endl;
cout << "q_s = " << q_s << endl;
cout << "z0 = " << z0 << endl;
double s_mag = sqrt(p_s*p_s+q_s*q_s+1.0);
printf("Direction of point source = [%g,%g,%g] = [%.3f,%.3f,%.3f]\n",
	-p_s,-q_s,1.0,-p_s/s_mag,-q_s/s_mag,1.0/s_mag);
cout << "---------------------------------------" 
	<< "---------------------------------------" << endl;

if (z0 < (r+lambda))
	{
	 cerr << "Fatal error: Part of sphere lies between the"
		<< " lens and the image plane!" << endl;
	 cerr << "z0 should be >= " << r+lambda << endl;
	 exit(-1);
	}

cout << "Radius of sphere on image plane = " << ri << endl;

Matrix w(4,1), c_h(4,1);
w.p[3][0] = 1.0;

int N_phi, N_theta;
if (ri < 1.0)
	{
	 N_phi = (int) rint(N_SAMPLE*ri);
	 N_theta = (int) rint(N_SAMPLE*ri);
	}
else
	{
	 N_phi = N_SAMPLE;
	 N_theta = N_SAMPLE;
	}

cout << "Number of points sampled for image plane = " << N_phi
	<< "x" << N_theta << " = " << N_phi*N_theta << endl;
cout << endl;

int width = IMAGE_DIM;
int height = IMAGE_DIM;
GrayImage image(width, height, 255);
image.setValue(BG_INTENSITY);

int count;
double phi, theta, R, x, y, z, xi, yi;
int xx,yy,RR;

cout << "Generating shading information....." << endl;
double start_time = Timer();
fprintf(stderr, "%3d%% Completed",0);

for (int i = 0; i < N_phi; i++)
	for (int j = 0; j < N_theta; j++)
	{
		 phi = theta_i+(180.0*DEG2RAD-2*theta_i)*i/(N_phi-1);
		 theta = theta_i+(180.0*DEG2RAD-2*theta_i)*j/(N_theta-1);
		 R = ReflectanceSphere(theta, phi, p_s, q_s);
		 w.p[0][0] = x = - r * sin(theta) * cos(phi);;
		 w.p[1][0] = y = r * cos(theta);;
		 w.p[2][0] = z = z0 - r * sin(theta) * sin(phi);;
		 c_h = T * w;
		 if (fabs(c_h(3,0)) > 1.0e-4)
			{
		 	 xi = c_h(0,0)/c_h(3,0);
		 	 yi = c_h(1,0)/c_h(3,0);
			}
		 else
			{
		 	 xi = 1.0e+10;
		 	 yi = 1.0e+10;
			}
//		 PerspectiveTransform(T,x,y,z,&xi,&yi);
		 xx = (int) rint((xi+XY_RANGE)/(2*XY_RANGE)*(IMAGE_DIM-1));
		 yy = (int) rint((yi+XY_RANGE)/(2*XY_RANGE)*(IMAGE_DIM-1));
		 RR = (int) rint(R*255);
		 if (xx >= 0 && xx < IMAGE_DIM && 
		     yy >= 0 && yy < IMAGE_DIM &&
		     sqrt(xi*xi+yi*yi) < (ri-4.0/IMAGE_DIM))
			image.p[yy][IMAGE_DIM-1-xx] = RR;
		 count = i*N_theta+j+1;
		 if (count%1000 == 0)
			fprintf(stderr, "\r%3d%%",100*count/(N_theta*N_phi));
		}
fprintf(stderr, "\r%3d%% Completed!\n",100);
double end_time = Timer();
fprintf(stderr,"Consumed %.1f seconds.\n",end_time-start_time);
cout << "---------------------------------------" 
	<< "---------------------------------------" << endl;
		
//image.reverseThreshold(255);
image.WriteToFile("/tmp/sphere.ppm");
cout << "---------------------------------------" 
	<< "---------------------------------------" << endl;
system("xv /tmp/sphere.ppm");
system("rm -f /tmp/sphere.ppm");
}
//**************************************************************************
double ReflectanceSphere (double theta, double phi, double p_s, double q_s)
{
double p = -cos(phi)*sin(theta);
double q = cos(theta);
double r = sin(phi)*sin(theta);

double R = (p*p_s+q*q_s+r)/sqrt(1.0+p_s*p_s+q_s*q_s);
if (R < 0.0) 
	R = 0.0;

return R;
}
//**************************************************************************
double Timer (void)
{
time_t time_now;

time(&time_now);

return (double) (time_now);
}
//**************************************************************************