// -*- mode: go; tab-width:4 -*-
package main

import (
	"math";
	"rand";
	"os";
	"fmt";
)


const (
	WIDTH = 256;
	HEIGHT = 256;
	NSUBSAMPLES = 2;
	NAO_SAMPLES = 8;
)

type Vec struct {
	x, y, z float64;
}

type Isect struct {
	t float64;
	p Vec;
	n Vec;
	hit int;
}

type Sphere struct {
	center Vec;
	radius float64;
}

type Plane struct {
	p Vec;
	n Vec;
}

type Ray struct {
	org Vec;
	dir Vec;
}

type Primitive interface {
	ray_intersect(isect *Isect, ray *Ray);
}


var primitives []Primitive;
var rnd *rand.Rand;


func (v0 *Vec) dot(v1 *Vec) float64 {
	return v0.x * v1.x + v0.y * v1.y + v0.z * v1.z;
}

func (c *Vec) cross(v0 *Vec, v1 *Vec) {
	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;
}

func (c *Vec) normalize() {
	var length float64 = math.Sqrt(c.dot(c));
	c.x /= length;
	c.y /= length;
	c.z /= length;
}

func (sphere Sphere) ray_intersect(isect *Isect, ray *Ray) {
	var rs Vec;
	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 := rs.dot(&ray.dir);
	C := rs.dot(&rs) - sphere.radius * sphere.radius;
	D := B * B - C;

	if D > 0 {
		t := -B - math.Sqrt(D);
		
		if t > 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;
			
			isect.n.normalize();
		}
	}
}

func (plane Plane) ray_intersect(isect *Isect, ray *Ray) {
	d := -plane.p.dot(&plane.n);
	v := ray.dir.dot(&plane.n);
	
	if math.Fabs(v) < 1.0e-17 {return}

	t := - (ray.org.dot(&plane.n) + d) / v;

	if t > 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;
	}
}


func orthoBasis(basis *[3]Vec, n Vec) {
	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;
	}
	
	basis[0].cross(&basis[1], &basis[2]);
	basis[0].normalize();
	
	basis[1].cross(&basis[2], &basis[0]);
	basis[1].normalize();
}




func ambient_occlusion(col *Vec, isect *Isect) {
	ntheta := NAO_SAMPLES;
	nphi := NAO_SAMPLES;
	const eps float64 = 0.0001;

	var p Vec;

	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;
	
	var basis [3]Vec;
	orthoBasis(&basis, isect.n);

	var occlusion float64 = 0.0;

	for j := 0; j < ntheta; j++ {
		for i := 0; i < nphi; i++ {
			theta := math.Sqrt(rnd.Float64());
			phi := 2.0 * math.Pi * rnd.Float64();
			
			x := math.Cos(phi) * theta;
			y := math.Sin(phi) * theta;
			z := math.Sqrt(1.0 - theta * theta);

			//fmt.Printf("%v %v %v\n", x, y, z);
			
			// local . global
			rx := x * basis[0].x + y * basis[1].x + z * basis[2].x;
			ry := x * basis[0].y + y * basis[1].y + z * basis[2].y;
			rz := x * basis[0].z + y * basis[1].z + z * basis[2].z;
			
			var ray Ray;
			
			ray.org = p;
			ray.dir.x = rx;
			ray.dir.y = ry;
			ray.dir.z = rz;
			
			var occIsect Isect;
			occIsect.t = 1.0e+17;
			occIsect.hit = 0;
			//fmt.Printf("%v\n", ray);
			
			/*
			 spheres[0].ray_intersect(&occIsect, &ray);
			 spheres[1].ray_intersect(&occIsect, &ray);
			 spheres[2].ray_intersect(&occIsect, &ray);
			 plane[0].ray_intersect (&occIsect, &ray);
			 */
			for n := 0; n < len(primitives); n++ {
				primitives[n].ray_intersect(&occIsect, &ray); 
			}
			
			//fmt.Printf("%v\n", occIsect);
			if occIsect.hit > 0 {
				occlusion += 1.0;
			}
			
		}
	}
	/*if occlusion > 0 {
	 fmt.Printf("%v\n", occlusion);
	 }*/
	occlusion = (float64(ntheta * nphi) - occlusion) / float64(ntheta * nphi);
	
	col.x = occlusion;
	col.y = occlusion;
	col.z = occlusion;
}


func clamp(f float64) uint8 {
	i := int(f * 255.5);

	if i < 0 {
		i = 0;
	}
	if i > 255 {
		i = 255;
	}
	
	return uint8(i);
}

func render(img []uint8, w int, h int, nsubsamples int) {
	fimg := make([]float64, w * h * 3);
	for i := 0; i < len(fimg); i++ {
		fimg[i] = 0;
	}

	for y := 0; y < h; y++ {
		for x := 0; x < w; x++ {

			for v := 0; v < nsubsamples; v++ {
				for u := 0; u < nsubsamples; u++ {
px := (float64(x) + (float64(u) / float64(nsubsamples)) - (float64(w) / 2.0)) / (float64(w) / 2.0);
					py := -(float64(y) + (float64(v) / float64(nsubsamples)) - (float64(h) / 2.0)) / (float64(h) / 2.0);

					var ray Ray;

					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;
					ray.dir.normalize();

					var isect Isect;
					isect.t = 1.0e+17;
					isect.hit = 0;
					/*
					 spheres[0].ray_intersect(&isect, &ray);
					 spheres[1].ray_intersect(&isect, &ray);
					 spheres[2].ray_intersect(&isect, &ray);
					 plane[0].ray_intersect (&isect, &ray);
					 */
					
					for n := 0; n < len(primitives); n++ {
						primitives[n].ray_intersect(&isect, &ray); 
					}

					if isect.hit > 0 {
						var col Vec;
						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] /= float64(nsubsamples * nsubsamples);
			fimg[3 * (y * w + x) + 1] /= float64(nsubsamples * nsubsamples);
			fimg[3 * (y * w + x) + 2] /= float64(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]);
		}
	}

}

func init_scene() {
	spheres := make([]Sphere, 3);
	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 := make([]Plane, 1);
	plane[0].p.x = 0.0;
	plane[0].p.y = -0.5;
	plane[0].p.z = 0.0;

	plane[0].n.x = 0.0;
	plane[0].n.y = 1.0;
	plane[0].n.z = 0.0;

	primitives = make([]Primitive, len(spheres) + len(plane));

	for i := 0; i < len(spheres); i++ {
		primitives[i] = spheres[i];
	}
	for i := 0; i < len(plane); i++ {
		primitives[i + len(spheres)] = plane[i];
	}
	//fmt.Printf("%v %v\n", primitives, len(primitives));
}

func saveppm(fname string, w int, h int, img []uint8) os.Error {
	file, err := os.Open(fname, os.O_CREAT | os.O_WRONLY, 0644);

	fmt.Fprintf(file, "P6\n");
	fmt.Fprintf(file, "%d %d\n", w, h);
	fmt.Fprintf(file, "255\n");
	file.Write(img);
	file.Close();

	return err;
}

func main() {
	sec, nsec, err := os.Time();
	if err != nil {
		fmt.Print(err);
		return;
	}
	rnd = rand.New(rand.NewSource(1e9 * sec + nsec));

	img := make([]uint8, WIDTH * HEIGHT * 3);

	init_scene();

	render(img, WIDTH, HEIGHT, NSUBSAMPLES);

	saveppm("ao.ppm", WIDTH, HEIGHT, img); 
}