// Ray class function declarations

#include "StdAfx.h"

// open file and initialize ray class
ray::ray(const wchar_t* fname)
{
	FILE *stream;
	char line[100];

	_wfopen_s( &stream, fname, L"r" );
	
	fgets(line, 100, stream);
	sscanf_s(line, "%d",&imgsz);
	fgets(line, 100, stream);
	sscanf_s(line, "%lf %lf",&amb, &diff);
	fgets(line, 100, stream);
	sscanf_s(line, "%lf %lf %lf",&light.x, &light.y, &light.z);
	fgets(line, 100, stream);
	sscanf_s(line, "%d", &kbr);
	fgets(line, 100, stream);
	sscanf_s(line, "%d", &kbs);
	fgets(line, 100, stream);
	sscanf_s(line, "%d", &knt);
	fgets(line, 100, stream);
	sscanf_s(line, "%d", &nballs);
	fgets(line, 100, stream);
	sscanf_s(line, "%lf", &iso);
#ifdef DEBUG	
	printf("%d %lf %lf %lf %lf %lf %d %d %d %d %lf\n", imgsz, amb, diff, light.x, light.y, light.z, kbr, kbs, knt, nballs, iso);
#endif
	ball=new metaball[nballs];	// create the metaballs
	for (int i=0; i<nballs; i++) {
		fgets(line, 100, stream);	// and initialize them
		sscanf_s(line, "%lf %lf %lf %lf %lf", &ball[i].r0.x, &ball[i].r0.y, &ball[i].r0.z, &ball[i].lambda, &ball[i].sigma);
#ifdef DEBUG
		printf("%lf %lf %lf %lf %lf \n", ball[i].r0.x, ball[i].r0.y, ball[i].r0.z, ball[i].lambda, ball[i].sigma);
#endif
	}
}

bool ray::findRoot(double x, double y, double &z ) {
	
	double a=-1.0;	// a,b endpoints
	double b=1.0;
	const double epsilon=2.0/(double)(imgsz-1); // root finding tolerance
	 	
	int i;
	
	// search for root (bracket the root)
	if (kbr) {
	
		double h=(b-a)/kbr;
		double y1=f(x,y,z=a);
		double y2;
		for ( i=0; i<kbr; i++) {
			//if (y1==0.0) {
			//	return true;	// exact root found
			//}
			y2=f(x,y,z+=h);
			//if (y2==0.0) {
			//	return true;	// exact root found
			//}
			if (y1*y2<0.0) {
				a=z-h;
				b=z;
#ifdef DEBUG 
				printf("a=%lf b=%lf\n",a,b);
#endif
				break; // root bracketed
			}
#ifdef DEBUG 
			printf("i=%d y1=%lf y2=%lf\n",i,y1,y2);
#endif 
			y1=y2;
		}
	}
	
	double u=f(x,y,a);
	double v=f(x,y,b);
	if (u * v > 0.0) return false; // unable to bracket root

		
	// bisecting (expects root is bracketed)

	if (kbs) {
		double w;
		double error=b-a;
		for (i=0; i<kbs; i++) {
			error/=2.0;
			z=a+error;
			w=f(x,y,z);
#ifdef DEBUG
			printf("i=%d\tz=%lg\tf(x)=%lg\te=%lg\n",i,z,w,error); 
#endif
			if (fabs(error)<epsilon) {
				return true;	 // root converged
			}
			if (w * u < 0.0) {
				b=z;
				v=w;
			}
			else {
				a=z;
				u=w;
			}
		}
	}
	// and finally, use sir Isac Newton method
	if (knt) {
		double d;
		for (i=0; i<knt; i++) {
			d=f(x,y,z)/fdz(x,y,z);
			z-=d;
#ifdef DEBUG
			printf("i=%d z=%lg\n",i,z);
#endif
			if (fabs(d)<epsilon) {
				return true;	// root converged
			}
		}
	}
	
	return false; // root did not converge
}

//The implicit function f
double ray::f(double x, double y, double z) {

	int i;
	double val=0.0;
	
	for (i=0; i<nballs; i++) {
			val += ball[i].gaussian(vector(x,y,z));
	}
	return (val - iso);
}
// df/dx
double ray::fdx(double x, double y, double z) {
	int i;
	double val=0.0;
	
	for (i=0; i<nballs; i++) {
			val += ball[i].dx_gaussian(vector(x,y,z));
	}
	return (val);
}


// df/dy
double ray::fdy(double x, double y, double z) {
	int i;
	double val=0.0;
	
	for (i=0; i<nballs; i++) {
			val += ball[i].dy_gaussian(vector(x,y,z));
	}
	return (val);
}

// df/dz
double ray::fdz(double x, double y, double z) {
	int i;
	double val=0.0;
	
	for (i=0; i<nballs; i++) {
			val += ball[i].dz_gaussian(vector(x,y,z));
	}
	return (val);
}

// numeric differentiation by central differences

double ray::nfdx(double x, double y, double z) {
	const double h=0.001;
	return ((f(x+h, y, z) - f(x-h, y, z)) / (2.0 * h)); 
}

double ray::nfdy(double x, double y, double z) {
	const double h=0.001;
	return ((f(x, y+h, z) - f(x, y-h, z)) / (2.0 * h)); 
}

double ray::nfdz(double x, double y, double z) {
	const double h=0.001;
	return ((f(x, y, z+h) - f(x, y, z-h)) / (2.0 * h)); 
}

double ray::illuminate(double x, double y, double z) {

	vector gradient(nfdx(x,y,z), nfdy(x,y,z), nfdz(x,y,z));
	gradient.norm();	// convert gradient to unit length
	double illum=gradient.dot(light);	// find cos of the angle
	if (illum<0.0) illum=0.0;	// normal points away
	illum=amb + diff * illum;	// add the constant ambient to the diffuse;
	//illum=0.5*(1.0-z);	// to see the root location;
	if (illum>1.0) illum=1.0;	// error chk
	if (illum<0.0) illum=0.0;	// error chk
	return illum;
}

// method for writing the image
void ray::writeImg(const wchar_t* fname, System::Windows::Forms::ProgressBar^ progressBar1) {

	const int nx=imgsz;
	const int ny=imgsz;
	double x,y,z;
	double dx=2.0/(double)(imgsz-1);
	double dy=dx;
	progressBar1->Maximum=imgsz;	// indicates current progress

	FILE *stream;
	_wfopen_s(&stream,fname,L"wb");
	fprintf(stream,"P6\n%d %d\n255\n",nx, ny);
	const unsigned char background[3] ={0,0,255};
	light.norm(); // normalize light vector
	y=-1.0;
	for(int iy = 0; iy < ny; iy++, y+=dy) {
		x=-1.0;
		progressBar1->Value=iy;	// update progress bar
		for (int ix = 0; ix < nx; ix++, x+=dx) {
			if (findRoot(x, y, z)) {	// root found ?
				double intensity = illuminate(x, y, z); // yes, calculate light itensity at root
				unsigned char gray[3];
				gray[0]=gray[1]=gray[2]=(int) (255.0 * intensity);
				fwrite(gray, 1, 3, stream);
			}
			else {
				fwrite(background, 1, 3, stream);	// root not found, use background color
			}
		}
	}
	fclose(stream);
}

