// -*- 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 spheres []Sphere;
var plane []Plane;
var primitives []Primitive;

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

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

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

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(isect.n);
		}
	}
}

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(n Vec) (basis [3]Vec) {
	basis[2].x = n.x;
	basis[2].y = n.y;
	basis[2].z = n.z;
	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] = normalize(cross(basis[1], basis[2]));
	basis[1] = normalize(cross(basis[2], basis[0]));

	return basis;
}

func ambient_occlusion(rnd *rand.Rand, isect Isect) (col Vec) {
	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;
	
	basis := orthoBasis(isect.n);

	var occlusion float64 = 0.0;
			
	var ray Ray;
	var occIsect Isect;

	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);

			// 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;
			
			ray.org = p;
			ray.dir.x = rx;
			ray.dir.y = ry;
			ray.dir.z = rz;
			
			occIsect.t = 1.0e+17;
			occIsect.hit = 0;

			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); 
			}
			 */

			if occIsect.hit > 0 {
				occlusion += 1.0;
			}
			
		}
	}
	occlusion = (float64(ntheta * nphi) - occlusion) / float64(ntheta * nphi);
	
	col.x = occlusion;
	col.y = occlusion;
	col.z = occlusion;

	return col;
}


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) {
	ch := make(chan int);
	for y := 0; y < h; y++ {
		go render_line(ch, img, y, w, h, nsubsamples);
	}

	for n := 0; n < h; n++ {
		<- ch;
	}
}

var freeRndList = make(chan *rand.Rand, 10)

func render_line(ch chan int, img []uint8, y int, w int, h int, nsubsamples int) {
	fimg := make([]float64, w * 3);

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

				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 {
					rnd, ok := <- freeRndList;
					if !ok {
						sec, nsec, _ := os.Time();
						rnd = rand.New(rand.NewSource(1e9 * sec + nsec));
					}
					col := ambient_occlusion(rnd, isect);
					_ = freeRndList <- rnd;			

					fimg[3 * x + 0] += col.x;
					fimg[3 * x + 1] += col.y;
					fimg[3 * x + 2] += col.z;
				}
				
			}
		}
		
		fimg[3 * x + 0] /= float64(nsubsamples * nsubsamples);
		fimg[3 * x + 1] /= float64(nsubsamples * nsubsamples);
		fimg[3 * x + 2] /= float64(nsubsamples * nsubsamples);
		
		img[3 * (y * w + x) + 0] = clamp(fimg[3 * x + 0]);
		img[3 * (y * w + x) + 1] = clamp(fimg[3 * x + 1]);
		img[3 * (y * w + x) + 2] = clamp(fimg[3 * x + 2]);
	}
	ch <- 0;
}

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];
	}
}

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() {
	img := make([]uint8, WIDTH * HEIGHT * 3);

	init_scene();

	render(img, WIDTH, HEIGHT, NSUBSAMPLES);

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