/*
0.642 sec on Core2 Quad Q9000 2.00GHz (VC++ 2010 Beta1, Windows 7)
based on "ao.c" @ AOBench
http://lucille.atso-net.jp/aobench/
*/

#include <windows.h>
#include <assert.h>
#include <ppl.h>

#include <sys/types.h>
#include <sys/timeb.h>

using namespace Concurrency;

static struct _timeb start, end;
static void reset_and_start_timer()
{
	_ftime(&start);
}
static double get_elapsed_seconds()
{
	_ftime(&end);
	return (1000 * (end.time - start.time) + (end.millitm - start.millitm)) / 1000.0;
}

static long long drand48_x = 0x1234ABCD330E;
static double drand48()
{
	drand48_x = drand48_x * 0x5DEECE66D + 0xB;
	return (drand48_x & 0xFFFFFFFFFFFF) * (1.0 / 281474976710656.0);
}
#define M_PI 3.14159265358979

#define WIDTH        256
#define HEIGHT       256
#define NSUBSAMPLES  2
#define NAO_SAMPLES  8

typedef struct _vec
{
	double x;
	double y;
	double z;
} vec;


typedef struct _Isect
{
	double t;
	vec    p;
	vec    n;
	int    hit; 
} Isect;

typedef struct _Sphere
{
	vec    center;
	double radius;

} Sphere;

typedef struct _Plane
{
	vec    p;
	vec    n;

} Plane;

typedef struct _Ray
{
	vec    org;
	vec    dir;
} Ray;

Sphere spheres[3];
Plane  plane;

static double vdot(vec v0, vec v1)
{
	return v0.x * v1.x + v0.y * v1.y + v0.z * v1.z;
}

static void vcross(vec *c, vec v0, vec v1)
{
	c->x = v0.y * v1.z - v0.z * v1.y;
	c->y = v0.z * v1.x - v0.x * v1.z;
	c->z = v0.x * v1.y - v0.y * v1.x;
}

static void vnormalize(vec *c)
{
	double length = sqrt(vdot((*c), (*c)));

	if (length > 1.0e-17) {
		c->x /= length;
		c->y /= length;
		c->z /= length;
	}
}

void
ray_sphere_intersect(Isect *isect, const Ray *ray, const Sphere *sphere)
{
	vec rs;
	double B, C, D;

	rs.x = ray->org.x - sphere->center.x;
	rs.y = ray->org.y - sphere->center.y;
	rs.z = ray->org.z - sphere->center.z;

	B = vdot(rs, ray->dir);
	C = vdot(rs, rs) - sphere->radius * sphere->radius;
	D = B * B - C;

	if (D > 0.0) {
		double t = -B - sqrt(D);

		if ((t > 0.0) && (t < isect->t)) {
			isect->t = t;
			isect->hit = 1;

			isect->p.x = ray->org.x + ray->dir.x * t;
			isect->p.y = ray->org.y + ray->dir.y * t;
			isect->p.z = ray->org.z + ray->dir.z * t;

			isect->n.x = isect->p.x - sphere->center.x;
			isect->n.y = isect->p.y - sphere->center.y;
			isect->n.z = isect->p.z - sphere->center.z;

			vnormalize(&(isect->n));
		}
	}
}

void
ray_plane_intersect(Isect *isect, const Ray *ray, const Plane *plane)
{
	double d, v, t;
	d = -vdot(plane->p, plane->n);
	v = vdot(ray->dir, plane->n);

	if (fabs(v) < 1.0e-17) return;

	t = -(vdot(ray->org, plane->n) + d) / v;

	if ((t > 0.0) && (t < isect->t)) {
		isect->t = t;
		isect->hit = 1;

		isect->p.x = ray->org.x + ray->dir.x * t;
		isect->p.y = ray->org.y + ray->dir.y * t;
		isect->p.z = ray->org.z + ray->dir.z * t;

		isect->n = plane->n;
	}
}

void
orthoBasis(vec *basis, vec n)
{
	basis[2] = n;
	basis[1].x = 0.0; basis[1].y = 0.0; basis[1].z = 0.0;

	if ((n.x < 0.6) && (n.x > -0.6)) {
		basis[1].x = 1.0;
	} else if ((n.y < 0.6) && (n.y > -0.6)) {
		basis[1].y = 1.0;
	} else if ((n.z < 0.6) && (n.z > -0.6)) {
		basis[1].z = 1.0;
	} else {
		basis[1].x = 1.0;
	}

	vcross(&basis[0], basis[1], basis[2]);
	vnormalize(&basis[0]);

	vcross(&basis[1], basis[2], basis[0]);
	vnormalize(&basis[1]);
}

void
ambient_occlusion(vec *col, const Isect *isect)
{
	int    i, j;
	int    ntheta = NAO_SAMPLES;
	int    nphi   = NAO_SAMPLES;
	double eps = 0.0001;

	vec p;
	vec basis[3];
	double occlusion = 0.0;

	p.x = isect->p.x + eps * isect->n.x;
	p.y = isect->p.y + eps * isect->n.y;
	p.z = isect->p.z + eps * isect->n.z;

	orthoBasis(basis, isect->n);

	for (j = 0; j < ntheta; j++) {
		for (i = 0; i < nphi; i++) {
			double theta = sqrt(drand48());
			double phi   = 2.0 * M_PI * drand48();

			double x = cos(phi) * theta;
			double y = sin(phi) * theta;
			double z = sqrt(1.0 - theta * theta);

			// local -> global
			double rx = x * basis[0].x + y * basis[1].x + z * basis[2].x;
			double ry = x * basis[0].y + y * basis[1].y + z * basis[2].y;
			double rz = x * basis[0].z + y * basis[1].z + z * basis[2].z;

			Ray ray;
			Isect occIsect;

			ray.org = p;
			ray.dir.x = rx;
			ray.dir.y = ry;
			ray.dir.z = rz;

			occIsect.t   = 1.0e+17;
			occIsect.hit = 0;

			ray_sphere_intersect(&occIsect, &ray, &spheres[0]); 
			ray_sphere_intersect(&occIsect, &ray, &spheres[1]); 
			ray_sphere_intersect(&occIsect, &ray, &spheres[2]); 
			ray_plane_intersect (&occIsect, &ray, &plane); 

			if (occIsect.hit) occlusion += 1.0;
		}
	}

	occlusion = (ntheta * nphi - occlusion) / (double)(ntheta * nphi);

	col->x = occlusion;
	col->y = occlusion;
	col->z = occlusion;
}

unsigned char
clamp(double f)
{
	int i = (int)(f * 255.5);

	if (i < 0) i = 0;
	if (i > 255) i = 255;

	return (unsigned char)i;
}

void
render(unsigned char *img, int w, int h, int nsubsamples)
{
//	int x, y;
//	int u, v;

	double *fimg = (double *)malloc(sizeof(double) * w * h * 3);
	memset((void *)fimg, 0, sizeof(double) * w * h * 3);
	parallel_for(0, h, [&](int y){
		parallel_for(0,w,[&](int x){
			for (int v = 0; v < nsubsamples; v++) {
				for (int u = 0; u < nsubsamples; u++) {
					double px = (x + (u / (double)nsubsamples) - (w / 2.0)) / (w / 2.0);
					double py = -(y + (v / (double)nsubsamples) - (h / 2.0)) / (h / 2.0);

					Ray ray;
					Isect isect;

					ray.org.x = 0.0;
					ray.org.y = 0.0;
					ray.org.z = 0.0;

					ray.dir.x = px;
					ray.dir.y = py;
					ray.dir.z = -1.0;
					vnormalize(&(ray.dir));

					isect.t   = 1.0e+17;
					isect.hit = 0;

					ray_sphere_intersect(&isect, &ray, &spheres[0]);
					ray_sphere_intersect(&isect, &ray, &spheres[1]);
					ray_sphere_intersect(&isect, &ray, &spheres[2]);
					ray_plane_intersect (&isect, &ray, &plane);

					if (isect.hit) {
						vec col;
						ambient_occlusion(&col, &isect);
						fimg[3 * (y * w + x) + 0] += col.x;
						fimg[3 * (y * w + x) + 1] += col.y;
						fimg[3 * (y * w + x) + 2] += col.z;
					}
				}
			}
			fimg[3 * (y * w + x) + 0] /= (double)(nsubsamples * nsubsamples);
			fimg[3 * (y * w + x) + 1] /= (double)(nsubsamples * nsubsamples);
			fimg[3 * (y * w + x) + 2] /= (double)(nsubsamples * nsubsamples);

			img[3 * (y * w + x) + 0] = clamp(fimg[3 *(y * w + x) + 0]);
			img[3 * (y * w + x) + 1] = clamp(fimg[3 *(y * w + x) + 1]);
			img[3 * (y * w + x) + 2] = clamp(fimg[3 *(y * w + x) + 2]);
		});
	});
}

void
init_scene()
{
	spheres[0].center.x = -2.0;
	spheres[0].center.y =  0.0;
	spheres[0].center.z = -3.5;
	spheres[0].radius = 0.5;

	spheres[1].center.x = -0.5;
	spheres[1].center.y =  0.0;
	spheres[1].center.z = -3.0;
	spheres[1].radius = 0.5;

	spheres[2].center.x =  1.0;
	spheres[2].center.y =  0.0;
	spheres[2].center.z = -2.2;
	spheres[2].radius = 0.5;

	plane.p.x = 0.0;
	plane.p.y = -0.5;
	plane.p.z = 0.0;

	plane.n.x = 0.0;
	plane.n.y = 1.0;
	plane.n.z = 0.0;

}

void
saveppm(const char *fname, int w, int h, unsigned char *img)
{
	FILE *fp;

	fp = fopen(fname, "wb");
	assert(fp);

	fprintf(fp, "P6\n");
	fprintf(fp, "%d %d\n", w, h);
	fprintf(fp, "255\n");
	fwrite(img, w * h * 3, 1, fp);
	fclose(fp);
}

int
main(int argc, char **argv)
{
	unsigned char *img = (unsigned char *)malloc(WIDTH * HEIGHT * 3);
	int i, N = 20;
	double avg = 0.0;
	for (i = 0; i < N; i++) {
		reset_and_start_timer();
		init_scene();
		render(img, WIDTH, HEIGHT, NSUBSAMPLES);//});
		avg += get_elapsed_seconds();
	}
	avg /= N;
	printf("Elapsed: %7.3f[sec] (N=%d)\n", avg, N);
	saveppm("ao.ppm", WIDTH, HEIGHT, img); 
	scanf_s(" ");
	return 0;
}