# HG changeset patch
# User John Tsiombikas <nuclear@member.fsf.org>
# Date 1402121657 -10800
# Node ID 53a98c148bf8302bbc1932996976aded4d0e6982
# Parent  e78f68d03ae9269c14c5d4df4b93c10c57a8ae81
- introduced SurfaceGeometry to carry all the geometric information input to
  BRDF sampling and evaluation functions.
- made Reflectance keep an (optional) pointer to its material
- simplified PhongRefl::sample_dir, with the help of SurfaceGeometry
- worked around microsoft's broken std::thread implementation's deadlock on join

diff -r e78f68d03ae9 -r 53a98c148bf8 liberebus/src/brdf.cc
--- a/liberebus/src/brdf.cc	Sat Jun 07 06:10:21 2014 +0300
+++ b/liberebus/src/brdf.cc	Sat Jun 07 09:14:17 2014 +0300
@@ -1,16 +1,85 @@
+#include <math.h>
 #include <algorithm>
 #include "brdf.h"
+#include "material.h"
 #include "erebus_impl.h"
 
+// --- class SurfaceGeometry ---
+SurfaceGeometry::SurfaceGeometry(const Vector3 &norm, VecLength st)
+	: SurfaceGeometry(norm, Vector2(0, 0), st)
+{
+}
+
+SurfaceGeometry::SurfaceGeometry(const Vector3 &norm, const Vector2 &uvarg, VecLength st)
+	: normal(norm), uv(uvarg)
+{
+	if(st == unknown) {
+		normal.normalize();
+	}
+
+	tangent = Vector3(1.0f, 0.0f, 0.0f);
+	if(fabs(dot_product(normal, tangent)) - 1.0f < 1e-4f) {
+		tangent = Vector3(0.0f, 0.0f, 1.0f);
+	}
+	Vector3 bitan = cross_product(normal, tangent);
+	tangent = cross_product(bitan, normal);
+}
+
+SurfaceGeometry::SurfaceGeometry(const Vector3 &norm, const Vector3 &tang, VecLength st)
+	: SurfaceGeometry(norm, tang, Vector2(0, 0), st)
+{
+}
+
+SurfaceGeometry::SurfaceGeometry(const Vector3 &norm, const Vector3 &tang, const Vector2 &uvarg, VecLength st)
+	: normal(norm), tangent(tang), uv(uvarg)
+{
+	if(st == unknown) {
+		normal.normalize();
+		tangent.normalize();
+	}
+}
+
+Vector3 SurfaceGeometry::spherical_to_world(float theta, float phi) const
+{
+	float x = cos(theta) * sin(phi);
+	float y = sin(theta) * sin(phi);
+	float z = cos(phi);
+	return sample_to_world(Vector3(x, y, z));
+}
+
+Vector3 SurfaceGeometry::sample_to_world(const Vector3 &v) const
+{
+	Matrix3x3 xform;
+	xform.set_column_vector(tangent, 0);
+	xform.set_column_vector(cross_product(normal, tangent), 1);
+	xform.set_column_vector(normal, 2);
+	return v.transformed(xform);
+}
+
+Vector3 SurfaceGeometry::world_to_sample(const Vector3 &v) const
+{
+	Matrix3x3 xform;
+	xform.set_row_vector(tangent, 0);
+	xform.set_row_vector(cross_product(normal, tangent), 1);
+	xform.set_row_vector(normal, 2);
+	return v.transformed(xform);
+}
+
 // ---- class Reflectance ----
 Reflectance::Reflectance()
 {
+	mtl = 0;
 }
 
-float Reflectance::sample(const Vector3 &norm, const Vector3 &outdir, Vector3 *indir) const
+Reflectance::Reflectance(const Material *mtl)
 {
-	*indir = sample_dir(norm, outdir);
-	return eval(norm, outdir, *indir);
+	this->mtl = mtl;
+}
+
+float Reflectance::sample(const SurfaceGeometry &geom, const Vector3 &outdir, Vector3 *indir) const
+{
+	*indir = sample_dir(geom, outdir);
+	return eval(geom, outdir, *indir);
 }
 
 // --- class CompositeRefl ----
@@ -19,7 +88,13 @@
 	valid_checked = false;
 }
 
-int CompositeRefl::pick_brdf(const Vector3 &norm, const Vector3 &outdir) const
+CompositeRefl::CompositeRefl(const Material *mtl)
+	: Reflectance(mtl)
+{
+	valid_checked = false;
+}
+
+int CompositeRefl::pick_brdf(const SurfaceGeometry &geom, const Vector3 &outdir) const
 {
 	if(sub_brdf.empty()) {
 		return -1;
@@ -34,7 +109,7 @@
 	float *cdf = (float*)alloca(brdf_count * sizeof *cdf);
 	for(int i=0; i<brdf_count; i++) {
 		const SubRefl &sub = sub_brdf[i];
-		float w = sub.weight_func ? sub.weight_func(norm, outdir) : sub.weight;
+		float w = sub.weight_func ? sub.weight_func(geom, outdir) : sub.weight;
 		cdf[i] = i > 0 ? cdf[i - 1] + w : w;
 	}
 
@@ -83,93 +158,74 @@
 	"warning: don't call CompositeRefl's sample_dir and eval separately\n"
 	"         it'll pick different brdfs each time! use sample instead\n";
 
-Vector3 CompositeRefl::sample_dir(const Vector3 &norm, const Vector3 &outdir) const
+Vector3 CompositeRefl::sample_dir(const SurfaceGeometry &geom, const Vector3 &outdir) const
 {
 	if(!warned_composite_usage) {
 		fputs(composite_usage_warnstr, stderr);
 		warned_composite_usage = true;
 	}
-	int bidx = pick_brdf(norm, outdir);
-	return sub_brdf[bidx].brdf->sample_dir(norm, outdir);
+	int bidx = pick_brdf(geom, outdir);
+	return sub_brdf[bidx].brdf->sample_dir(geom, outdir);
 }
 
-float CompositeRefl::sample(const Vector3 &norm, const Vector3 &outdir, Vector3 *indir) const
+float CompositeRefl::sample(const SurfaceGeometry &geom, const Vector3 &outdir, Vector3 *indir) const
 {
-	int bidx = pick_brdf(norm, outdir);
-	return sub_brdf[bidx].brdf->sample(norm, outdir, indir);
+	int bidx = pick_brdf(geom, outdir);
+	return sub_brdf[bidx].brdf->sample(geom, outdir, indir);
 }
 
-float CompositeRefl::eval(const Vector3 &norm, const Vector3 &outdir, const Vector3 &indir) const
+float CompositeRefl::eval(const SurfaceGeometry &geom, const Vector3 &outdir, const Vector3 &indir) const
 {
 	if(!warned_composite_usage) {
 		fputs(composite_usage_warnstr, stderr);
 		warned_composite_usage = true;
 	}
-	int bidx = pick_brdf(norm, outdir);
-	return sub_brdf[bidx].brdf->eval(norm, outdir, indir);
+	int bidx = pick_brdf(geom, outdir);
+	return sub_brdf[bidx].brdf->eval(geom, outdir, indir);
 }
 
 // --- class LambertRefl ---
-Vector3 LambertRefl::sample_dir(const Vector3 &norm, const Vector3 &outdir) const
+Vector3 LambertRefl::sample_dir(const SurfaceGeometry &geom, const Vector3 &outdir) const
 {
 	Vector3 dir = Vector3{randf(-1, 1), randf(-1, 1), randf(-1, 1)}.normalized();
-	return dot_product(dir, norm) < 0.0 ? -dir : dir;
+	return dot_product(dir, geom.normal) < 0.0 ? -dir : dir;
 }
 
-float LambertRefl::eval(const Vector3 &norm, const Vector3 &outdir, const Vector3 &indir) const
+float LambertRefl::eval(const SurfaceGeometry &geom, const Vector3 &outdir, const Vector3 &indir) const
 {
-	return dot_product(norm, outdir);
+	return dot_product(geom.normal, outdir);
 }
 
 // --- class MirrorRefl ---
-Vector3 MirrorRefl::sample_dir(const Vector3 &norm, const Vector3 &outdir) const
+Vector3 MirrorRefl::sample_dir(const SurfaceGeometry &geom, const Vector3 &outdir) const
 {
-	return outdir.reflection(norm);
+	return outdir.reflection(geom.normal);
 }
 
-float MirrorRefl::eval(const Vector3 &norm, const Vector3 &outdir, const Vector3 &indir) const
+float MirrorRefl::eval(const SurfaceGeometry &geom, const Vector3 &outdir, const Vector3 &indir) const
 {
 	return 1.0f;
 }
 
 // --- class PhongRefl ---
-PhongRefl::PhongRefl()
+Vector3 PhongRefl::sample_dir(const SurfaceGeometry &geom, const Vector3 &outdir) const
 {
-	shininess = 42.0;
-}
+	Vector3 refl = outdir.reflection(geom.normal).normalized();
+	SurfaceGeometry refl_geom(refl, SurfaceGeometry::unit);
 
-Vector3 PhongRefl::sample_dir(const Vector3 &norm, const Vector3 &outdir) const
-{
-	// construct orthonormal basis with k being the reflection dir
-	Vector3 refl = outdir.reflection(norm).normalized();
-	Vector3 up = Vector3(0, 0, 1);
-	if(fabs(dot_product(refl, up)) - 1.0 < 0.001) {
-		up = Vector3(0, 1, 0);
-	}
-	Vector3 right = cross_product(up, refl);
-	up = cross_product(refl, right);
-
-	// construct the matrix transposed to move the sample from reflection
-	// space to world space
-	Matrix3x3 xform;
-	xform.set_column_vector(right, 0);
-	xform.set_column_vector(up, 1);
-	xform.set_column_vector(refl, 2);
+	float shininess = mtl ? mtl->get_attrib_value("shininess", geom.uv.x, geom.uv.y) : 42.0;
 
 	float phi = acos(pow(randf(), 1.0 / (shininess + 1)));
 	float theta = 2.0 * M_PI * randf();
 
-	Vector3 v;
-	v.x = cos(theta) * sin(phi);
-	v.y = sin(theta) * sin(phi);
-	v.z = cos(phi);
-	v.transform(xform);
-	return v;
+	return refl_geom.spherical_to_world(theta, phi);
 }
 
-float PhongRefl::eval(const Vector3 &norm, const Vector3 &outdir, const Vector3 &indir) const
+float PhongRefl::eval(const SurfaceGeometry &geom, const Vector3 &outdir, const Vector3 &indir) const
 {
-	Vector3 refl = outdir.reflection(norm);
+	float shininess = mtl ? mtl->get_attrib_value("shininess", geom.uv.x, geom.uv.y) : 42.0;
+
+	Vector3 refl = outdir.reflection(geom.normal);
 	float dot = std::max<float>(dot_product(indir, refl), 0.0f);
 	return pow(dot, shininess);
 }
diff -r e78f68d03ae9 -r 53a98c148bf8 liberebus/src/brdf.h
--- a/liberebus/src/brdf.h	Sat Jun 07 06:10:21 2014 +0300
+++ b/liberebus/src/brdf.h	Sat Jun 07 09:14:17 2014 +0300
@@ -4,18 +4,56 @@
 #include <vector>
 #include "texture.h"
 
+class Material;
+
+class SurfaceGeometry {
+public:
+	Vector3 normal, tangent;
+	Vector2 uv;
+
+	enum VecLength { unit, unknown };
+
+	explicit SurfaceGeometry(const Vector3 &norm, VecLength st = unknown);
+	SurfaceGeometry(const Vector3 &norm, const Vector2 &uv, VecLength st = unknown);
+	SurfaceGeometry(const Vector3 &norm, const Vector3 &tang, VecLength st = unknown);
+	SurfaceGeometry(const Vector3 &norm, const Vector3 &tang, const Vector2 &uv, VecLength st = unknown);
+
+	/** create a cartesian direction vector in sample space (zenith = +Z) and
+	 * transform it to world space.
+	 * \param theta the horizontal angle (azimuth, in radians) around the Z axis [0, 2pi]
+	 * \param phi the vertical angle (elevation, in radians) away from the Z axis [0, pi]
+	 */
+	Vector3 spherical_to_world(float theta, float phi) const;
+	/// transforms a direction vector from sample space (centered around Z) to world space
+	Vector3 sample_to_world(const Vector3 &v) const;
+	/// transforms a direction vector from world space to sample space (centered around Z)
+	Vector3 world_to_sample(const Vector3 &v) const;
+};
+
+/// abstract bidirection reflectance distribution function base class
 class Reflectance {
+protected:
+	const Material *mtl;	// pointer to the material we belong to
+
 public:
 	Reflectance();
+	explicit Reflectance(const Material *mtl);
 	virtual ~Reflectance() = default;
 
-	virtual Vector3 sample_dir(const Vector3 &norm, const Vector3 &outdir) const = 0;
-	virtual float sample(const Vector3 &norm, const Vector3 &outdir, Vector3 *indir) const;
-	virtual float eval(const Vector3 &norm, const Vector3 &outdir, const Vector3 &indir) const = 0;
+	void set_material(const Material *mtl);
+	const Material *get_material() const;
+
+	/// given an outgoing light direction generate an incidence direction to sample the BRDF
+	virtual Vector3 sample_dir(const SurfaceGeometry &geom, const Vector3 &outdir) const = 0;
+	/// given an outgoing direction, generate an incidence direction, and evaluate its probability
+	virtual float sample(const SurfaceGeometry &geom, const Vector3 &outdir, Vector3 *indir) const;
+	/// given an outgoing direction and an incidence direction, evaluate the probability of this path
+	virtual float eval(const SurfaceGeometry &geom, const Vector3 &outdir, const Vector3 &indir) const = 0;
 };
 
-typedef float (*CompReflWeightFunc)(const Vector3 &norm, const Vector3 &outdir);
+typedef float (*CompReflWeightFunc)(const SurfaceGeometry &geom, const Vector3 &outdir);
 
+/// composite BRDF, with multiple weighted sub-reflectances.
 class CompositeRefl : public Reflectance {
 private:
 	struct SubRefl {
@@ -26,41 +64,41 @@
 	std::vector<SubRefl> sub_brdf;
 	mutable bool valid_checked;
 
-	int pick_brdf(const Vector3 &norm, const Vector3 &outdir) const;
+	int pick_brdf(const SurfaceGeometry &geom, const Vector3 &outdir) const;
 
 public:
 	CompositeRefl();
+	explicit CompositeRefl(const Material *mtl);
 
 	virtual void add_brdf(Reflectance *brdf, float weight);
 	virtual void add_brdf(Reflectance *brdf, CompReflWeightFunc weight_func);
 
 	bool check_valid() const;
 
-	Vector3 sample_dir(const Vector3 &norm, const Vector3 &outdir) const override;
-	float sample(const Vector3 &norm, const Vector3 &outdir, Vector3 *indir) const override;
-	float eval(const Vector3 &norm, const Vector3 &outdir, const Vector3 &indir) const override;
+	Vector3 sample_dir(const SurfaceGeometry &geom, const Vector3 &outdir) const override;
+	float sample(const SurfaceGeometry &geom, const Vector3 &outdir, Vector3 *indir) const override;
+	float eval(const SurfaceGeometry &geom, const Vector3 &outdir, const Vector3 &indir) const override;
 };
 
+/// lambertian perfect diffuse reflectance
 class LambertRefl : public Reflectance {
 public:
-	Vector3 sample_dir(const Vector3 &norm, const Vector3 &outdir) const override;
-	float eval(const Vector3 &norm, const Vector3 &outdir, const Vector3 &indir) const override;
+	Vector3 sample_dir(const SurfaceGeometry &geom, const Vector3 &outdir) const override;
+	float eval(const SurfaceGeometry &geom, const Vector3 &outdir, const Vector3 &indir) const override;
 };
 
+/// perfect specular reflectance
 class MirrorRefl : public Reflectance {
 public:
-	Vector3 sample_dir(const Vector3 &norm, const Vector3 &outdir) const override;
-	float eval(const Vector3 &norm, const Vector3 &outdir, const Vector3 &indir) const override;
+	Vector3 sample_dir(const SurfaceGeometry &geom, const Vector3 &outdir) const override;
+	float eval(const SurfaceGeometry &geom, const Vector3 &outdir, const Vector3 &indir) const override;
 };
 
+/// glossy phong reflectance with lafortune sampling
 class PhongRefl : public Reflectance {
 public:
-	float shininess;
-
-	PhongRefl();
-
-	Vector3 sample_dir(const Vector3 &norm, const Vector3 &outdir) const override;
-	float eval(const Vector3 &norm, const Vector3 &outdir, const Vector3 &indir) const override;
+	Vector3 sample_dir(const SurfaceGeometry &geom, const Vector3 &outdir) const override;
+	float eval(const SurfaceGeometry &geom, const Vector3 &outdir, const Vector3 &indir) const override;
 };
 
 #endif	// BRDF_H_
diff -r e78f68d03ae9 -r 53a98c148bf8 liberebus/src/erebus.cc
--- a/liberebus/src/erebus.cc	Sat Jun 07 06:10:21 2014 +0300
+++ b/liberebus/src/erebus.cc	Sat Jun 07 09:14:17 2014 +0300
@@ -22,11 +22,13 @@
 
 struct erebus *erb_init(void)
 {
-	struct erebus *ctx;
+	struct erebus *ctx = 0;
 	try {
 		ctx = new struct erebus;
+		ctx->tpool = new ThreadPool;
 	}
 	catch(...) {
+		delete ctx;
 		return 0;
 	}
 
@@ -47,7 +49,11 @@
 
 void erb_destroy(struct erebus *ctx)
 {
-	delete ctx;
+	if(ctx) {
+		// make sure the threadpool stops BEFORE destroying the framebuffers etc in ctx
+		delete ctx->tpool;
+		delete ctx;
+	}
 }
 
 void erb_setopti(struct erebus *ctx, enum erb_option opt, int val)
@@ -154,9 +160,9 @@
 
 int erb_render_rect(struct erebus *ctx, int x, int y, int width, int height, long timeout)
 {
-	while(ctx->tpool.pending()) {
+	while(ctx->tpool->pending()) {
 		if(timeout > 0) {
-			long wait_interval = ctx->tpool.wait(timeout);
+			long wait_interval = ctx->tpool->wait(timeout);
 			timeout -= wait_interval;
 		} else {
 			return 1;
@@ -180,7 +186,7 @@
 			blk.sample = ctx->cur_sample;
 			blk.frame = ctx->cur_frame;
 
-			ctx->tpool.add_work(std::bind(render_block, ctx, blk));
+			ctx->tpool->add_work(std::bind(render_block, ctx, blk));
 
 			x += BLKSZ;
 		}
@@ -188,7 +194,7 @@
 	}
 
 	++ctx->cur_sample;
-	ctx->tpool.wait(timeout);	// wait for completion
+	ctx->tpool->wait(timeout);	// wait for completion
 	return ctx->cur_sample > erb_getopti(ctx, ERB_OPT_MAX_SAMPLES) ? 0 : 1;
 }
 
diff -r e78f68d03ae9 -r 53a98c148bf8 liberebus/src/erebus_impl.h
--- a/liberebus/src/erebus_impl.h	Sat Jun 07 06:10:21 2014 +0300
+++ b/liberebus/src/erebus_impl.h	Sat Jun 07 09:14:17 2014 +0300
@@ -28,7 +28,7 @@
 	Image<float> accum;	// sample accumulator per pixel
 	Option options[ERB_NUM_OPTIONS];
 
-	ThreadPool tpool;
+	ThreadPool *tpool;
 
 	// render state
 	float inv_gamma;
diff -r e78f68d03ae9 -r 53a98c148bf8 liberebus/src/rt.cc
--- a/liberebus/src/rt.cc	Sat Jun 07 06:10:21 2014 +0300
+++ b/liberebus/src/rt.cc	Sat Jun 07 09:14:17 2014 +0300
@@ -46,7 +46,7 @@
 	float shininess = mtl->get_attrib_value("shininess");
 
 	Vector3 sample_dir;
-	float prob = brdf->sample(norm, -hit.world_ray.dir, &sample_dir);
+	float prob = brdf->sample(SurfaceGeometry(norm), -hit.world_ray.dir, &sample_dir);
 	if(iter < max_iter && randf() <= prob && ray.energy * prob > 0.001) {
 		Ray sample_ray;
 		sample_ray.origin = ray.origin + ray.dir * hit.dist;
diff -r e78f68d03ae9 -r 53a98c148bf8 liberebus/src/threadpool.cc
--- a/liberebus/src/threadpool.cc	Sat Jun 07 06:10:21 2014 +0300
+++ b/liberebus/src/threadpool.cc	Sat Jun 07 09:14:17 2014 +0300
@@ -19,23 +19,52 @@
 	thread = new std::thread[num_threads];
 	for(int i=0; i<num_threads; i++) {
 		thread[i] = std::thread(&ThreadPool::thread_func, this);
+
+#ifdef _MSC_VER
+		/* detach the thread to avoid having to join them in the destructor, which
+		 * causes a deadlock in msvc implementation when called after main returns
+		 */
+		thread[i].detach();
+#endif
 	}
 	this->num_threads = num_threads;
 }
 
 ThreadPool::~ThreadPool()
 {
+#ifdef _MSC_VER
+	workq_mutex.lock();
+	workq.clear();
+	qsize = 0;
+	workq_mutex.unlock();
+#endif
+
 	quit = true;
 	workq_condvar.notify_all();
 
 	printf("ThreadPool: waiting for %d worker threads to stop ", num_threads);
 	fflush(stdout);
+#ifndef _MSC_VER
 	for(int i=0; i<num_threads; i++) {
 		thread[i].join();
 		putchar('.');
 		fflush(stdout);
 	}
+
+#else
+	// spin until all threads are done...
+	std::unique_lock<std::mutex> lock(workq_mutex);
+	while(nactive > 0) {
+		lock.unlock();
+		std::this_thread::sleep_for(std::chrono::milliseconds(128));
+		putchar('.');
+		fflush(stdout);
+		lock.lock();
+	}
+#endif	// _MSC_VER
+
 	putchar('\n');
+	delete [] thread;
 }
 
 void ThreadPool::add_work(std::function<void ()> func)