Added xatlas as alternative to thekla, forced it on.

Did some hacks to it to avoid it from failing on bad geometry.

1 files changed, 7384 insertions, 0 deletions

diff --git a/thirdparty/xatlas/xatlas.cpp b/thirdparty/xatlas/xatlas.cpp
new file mode 100644
index 0000000000..f6a9ce64dc
--- /dev/null
+++ b/thirdparty/xatlas/xatlas.cpp
@@ -0,0 +1,7384 @@
+// This code is in the public domain -- castanyo@yahoo.es
+#include "xatlas.h"
+#include <assert.h>
+#include <float.h>
+#include <math.h>
+#include <stdarg.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <string.h>
+#include <time.h>
+#include <algorithm>
+#include <cmath>
+#include <memory>
+#include <unordered_map>
+#include <vector>
+
+#undef min
+#undef max
+
+#ifndef xaAssert
+#define xaAssert(exp)                                    \
+	if (!(exp)) {                                        \
+		xaPrint("%s %s %s\n", #exp, __FILE__, __LINE__); \
+	}
+#endif
+#ifndef xaDebugAssert
+#define xaDebugAssert(exp) assert(exp)
+#endif
+#ifndef xaPrint
+#define xaPrint(...)                            \
+	if (xatlas::internal::s_print) {            \
+		xatlas::internal::s_print(__VA_ARGS__); \
+	}
+#endif
+
+#ifdef _MSC_VER
+// Ignore gcc attributes.
+#define __attribute__(X)
+#endif
+
+#ifdef _MSC_VER
+#define restrict
+#define NV_FORCEINLINE __forceinline
+#else
+#define restrict __restrict__
+#define NV_FORCEINLINE __attribute__((always_inline)) inline
+#endif
+
+#define NV_UINT32_MAX 0xffffffff
+#define NV_FLOAT_MAX 3.402823466e+38F
+
+#ifndef PI
+#define PI float(3.1415926535897932384626433833)
+#endif
+
+#define NV_EPSILON (0.0001f)
+#define NV_NORMAL_EPSILON (0.001f)
+
+namespace xatlas {
+namespace internal {
+
+static PrintFunc s_print = NULL;
+
+static int align(int x, int a) {
+	return (x + a - 1) & ~(a - 1);
+}
+
+static bool isAligned(int x, int a) {
+	return (x & (a - 1)) == 0;
+}
+
+/// Return the maximum of the three arguments.
+template <typename T>
+static T max3(const T &a, const T &b, const T &c) {
+	return std::max(a, std::max(b, c));
+}
+
+/// Return the maximum of the three arguments.
+template <typename T>
+static T min3(const T &a, const T &b, const T &c) {
+	return std::min(a, std::min(b, c));
+}
+
+/// Clamp between two values.
+template <typename T>
+static T clamp(const T &x, const T &a, const T &b) {
+	return std::min(std::max(x, a), b);
+}
+
+static float saturate(float f) {
+	return clamp(f, 0.0f, 1.0f);
+}
+
+// Robust floating point comparisons:
+// http://realtimecollisiondetection.net/blog/?p=89
+static bool equal(const float f0, const float f1, const float epsilon = NV_EPSILON) {
+	//return fabs(f0-f1) <= epsilon;
+	return fabs(f0 - f1) <= epsilon * max3(1.0f, fabsf(f0), fabsf(f1));
+}
+
+NV_FORCEINLINE static int ftoi_floor(float val) {
+	return (int)val;
+}
+
+NV_FORCEINLINE static int ftoi_ceil(float val) {
+	return (int)ceilf(val);
+}
+
+NV_FORCEINLINE static int ftoi_round(float f) {
+	return int(floorf(f + 0.5f));
+}
+
+static bool isZero(const float f, const float epsilon = NV_EPSILON) {
+	return fabs(f) <= epsilon;
+}
+
+static float lerp(float f0, float f1, float t) {
+	const float s = 1.0f - t;
+	return f0 * s + f1 * t;
+}
+
+static float square(float f) {
+	return f * f;
+}
+
+static int square(int i) {
+	return i * i;
+}
+
+/** Return the next power of two.
+* @see http://graphics.stanford.edu/~seander/bithacks.html
+* @warning Behaviour for 0 is undefined.
+* @note isPowerOfTwo(x) == true -> nextPowerOfTwo(x) == x
+* @note nextPowerOfTwo(x) = 2 << log2(x-1)
+*/
+static uint32_t nextPowerOfTwo(uint32_t x) {
+	xaDebugAssert(x != 0);
+	// On modern CPUs this is supposed to be as fast as using the bsr instruction.
+	x--;
+	x |= x >> 1;
+	x |= x >> 2;
+	x |= x >> 4;
+	x |= x >> 8;
+	x |= x >> 16;
+	return x + 1;
+}
+
+static uint64_t nextPowerOfTwo(uint64_t x) {
+	xaDebugAssert(x != 0);
+	uint32_t p = 1;
+	while (x > p) {
+		p += p;
+	}
+	return p;
+}
+
+static uint32_t sdbmHash(const void *data_in, uint32_t size, uint32_t h = 5381) {
+	const uint8_t *data = (const uint8_t *)data_in;
+	uint32_t i = 0;
+	while (i < size) {
+		h = (h << 16) + (h << 6) - h + (uint32_t)data[i++];
+	}
+	return h;
+}
+
+// Note that this hash does not handle NaN properly.
+static uint32_t sdbmFloatHash(const float *f, uint32_t count, uint32_t h = 5381) {
+	for (uint32_t i = 0; i < count; i++) {
+		union {
+			float f;
+			uint32_t i;
+		} x = { f[i] };
+		if (x.i == 0x80000000) x.i = 0;
+		h = sdbmHash(&x, 4, h);
+	}
+	return h;
+}
+
+template <typename T>
+static uint32_t hash(const T &t, uint32_t h = 5381) {
+	return sdbmHash(&t, sizeof(T), h);
+}
+
+static uint32_t hash(const float &f, uint32_t h) {
+	return sdbmFloatHash(&f, 1, h);
+}
+
+// Functors for hash table:
+template <typename Key>
+struct Hash {
+	uint32_t operator()(const Key &k) const { return hash(k); }
+};
+
+template <typename Key>
+struct Equal {
+	bool operator()(const Key &k0, const Key &k1) const { return k0 == k1; }
+};
+
+class Vector2 {
+public:
+	typedef Vector2 const &Arg;
+
+	Vector2() {}
+	explicit Vector2(float f) :
+			x(f),
+			y(f) {}
+	Vector2(float x, float y) :
+			x(x),
+			y(y) {}
+	Vector2(Vector2::Arg v) :
+			x(v.x),
+			y(v.y) {}
+
+	const Vector2 &operator=(Vector2::Arg v) {
+		x = v.x;
+		y = v.y;
+		return *this;
+	}
+	const float *ptr() const { return &x; }
+
+	void set(float _x, float _y) {
+		x = _x;
+		y = _y;
+	}
+
+	Vector2 operator-() const {
+		return Vector2(-x, -y);
+	}
+
+	void operator+=(Vector2::Arg v) {
+		x += v.x;
+		y += v.y;
+	}
+
+	void operator-=(Vector2::Arg v) {
+		x -= v.x;
+		y -= v.y;
+	}
+
+	void operator*=(float s) {
+		x *= s;
+		y *= s;
+	}
+
+	void operator*=(Vector2::Arg v) {
+		x *= v.x;
+		y *= v.y;
+	}
+
+	friend bool operator==(Vector2::Arg a, Vector2::Arg b) {
+		return a.x == b.x && a.y == b.y;
+	}
+
+	friend bool operator!=(Vector2::Arg a, Vector2::Arg b) {
+		return a.x != b.x || a.y != b.y;
+	}
+
+	union {
+#ifdef _MSC_VER
+#pragma warning(push)
+#pragma warning(disable : 4201)
+#endif
+		struct
+		{
+			float x, y;
+		};
+#ifdef _MSC_VER
+#pragma warning(pop)
+#endif
+
+		float component[2];
+	};
+};
+
+Vector2 operator+(Vector2::Arg a, Vector2::Arg b) {
+	return Vector2(a.x + b.x, a.y + b.y);
+}
+
+Vector2 operator-(Vector2::Arg a, Vector2::Arg b) {
+	return Vector2(a.x - b.x, a.y - b.y);
+}
+
+Vector2 operator*(Vector2::Arg v, float s) {
+	return Vector2(v.x * s, v.y * s);
+}
+
+Vector2 operator*(Vector2::Arg v1, Vector2::Arg v2) {
+	return Vector2(v1.x * v2.x, v1.y * v2.y);
+}
+
+Vector2 operator/(Vector2::Arg v, float s) {
+	return Vector2(v.x / s, v.y / s);
+}
+
+Vector2 lerp(Vector2::Arg v1, Vector2::Arg v2, float t) {
+	const float s = 1.0f - t;
+	return Vector2(v1.x * s + t * v2.x, v1.y * s + t * v2.y);
+}
+
+float dot(Vector2::Arg a, Vector2::Arg b) {
+	return a.x * b.x + a.y * b.y;
+}
+
+float lengthSquared(Vector2::Arg v) {
+	return v.x * v.x + v.y * v.y;
+}
+
+float length(Vector2::Arg v) {
+	return sqrtf(lengthSquared(v));
+}
+
+float distance(Vector2::Arg a, Vector2::Arg b) {
+	return length(a - b);
+}
+
+bool isNormalized(Vector2::Arg v, float epsilon = NV_NORMAL_EPSILON) {
+	return equal(length(v), 1, epsilon);
+}
+
+Vector2 normalize(Vector2::Arg v, float epsilon = NV_EPSILON) {
+	float l = length(v);
+	xaDebugAssert(!isZero(l, epsilon));
+#ifdef NDEBUG
+	epsilon = 0; // silence unused parameter warning
+#endif
+	Vector2 n = v * (1.0f / l);
+	xaDebugAssert(isNormalized(n));
+	return n;
+}
+
+Vector2 normalizeSafe(Vector2::Arg v, Vector2::Arg fallback, float epsilon = NV_EPSILON) {
+	float l = length(v);
+	if (isZero(l, epsilon)) {
+		return fallback;
+	}
+	return v * (1.0f / l);
+}
+
+bool equal(Vector2::Arg v1, Vector2::Arg v2, float epsilon = NV_EPSILON) {
+	return equal(v1.x, v2.x, epsilon) && equal(v1.y, v2.y, epsilon);
+}
+
+Vector2 max(Vector2::Arg a, Vector2::Arg b) {
+	return Vector2(std::max(a.x, b.x), std::max(a.y, b.y));
+}
+
+bool isFinite(Vector2::Arg v) {
+	return std::isfinite(v.x) && std::isfinite(v.y);
+}
+
+// Note, this is the area scaled by 2!
+float triangleArea(Vector2::Arg v0, Vector2::Arg v1) {
+	return (v0.x * v1.y - v0.y * v1.x); // * 0.5f;
+}
+float triangleArea(Vector2::Arg a, Vector2::Arg b, Vector2::Arg c) {
+	// IC: While it may be appealing to use the following expression:
+	//return (c.x * a.y + a.x * b.y + b.x * c.y - b.x * a.y - c.x * b.y - a.x * c.y); // * 0.5f;
+	// That's actually a terrible idea. Small triangles far from the origin can end up producing fairly large floating point
+	// numbers and the results becomes very unstable and dependent on the order of the factors.
+	// Instead, it's preferable to subtract the vertices first, and multiply the resulting small values together. The result
+	// in this case is always much more accurate (as long as the triangle is small) and less dependent of the location of
+	// the triangle.
+	//return ((a.x - c.x) * (b.y - c.y) - (a.y - c.y) * (b.x - c.x)); // * 0.5f;
+	return triangleArea(a - c, b - c);
+}
+
+float triangleArea2(Vector2::Arg v1, Vector2::Arg v2, Vector2::Arg v3) {
+	return 0.5f * (v3.x * v1.y + v1.x * v2.y + v2.x * v3.y - v2.x * v1.y - v3.x * v2.y - v1.x * v3.y);
+}
+
+static uint32_t hash(const Vector2 &v, uint32_t h) {
+	return sdbmFloatHash(v.component, 2, h);
+}
+
+class Vector3 {
+public:
+	typedef Vector3 const &Arg;
+
+	Vector3() {}
+	explicit Vector3(float f) :
+			x(f),
+			y(f),
+			z(f) {}
+	Vector3(float x, float y, float z) :
+			x(x),
+			y(y),
+			z(z) {}
+	Vector3(Vector2::Arg v, float z) :
+			x(v.x),
+			y(v.y),
+			z(z) {}
+	Vector3(Vector3::Arg v) :
+			x(v.x),
+			y(v.y),
+			z(v.z) {}
+
+	const Vector3 &operator=(Vector3::Arg v) {
+		x = v.x;
+		y = v.y;
+		z = v.z;
+		return *this;
+	}
+
+	Vector2 xy() const {
+		return Vector2(x, y);
+	}
+
+	const float *ptr() const { return &x; }
+
+	void set(float _x, float _y, float _z) {
+		x = _x;
+		y = _y;
+		z = _z;
+	}
+
+	Vector3 operator-() const {
+		return Vector3(-x, -y, -z);
+	}
+
+	void operator+=(Vector3::Arg v) {
+		x += v.x;
+		y += v.y;
+		z += v.z;
+	}
+
+	void operator-=(Vector3::Arg v) {
+		x -= v.x;
+		y -= v.y;
+		z -= v.z;
+	}
+
+	void operator*=(float s) {
+		x *= s;
+		y *= s;
+		z *= s;
+	}
+
+	void operator/=(float s) {
+		float is = 1.0f / s;
+		x *= is;
+		y *= is;
+		z *= is;
+	}
+
+	void operator*=(Vector3::Arg v) {
+		x *= v.x;
+		y *= v.y;
+		z *= v.z;
+	}
+
+	void operator/=(Vector3::Arg v) {
+		x /= v.x;
+		y /= v.y;
+		z /= v.z;
+	}
+
+	friend bool operator==(Vector3::Arg a, Vector3::Arg b) {
+		return a.x == b.x && a.y == b.y && a.z == b.z;
+	}
+
+	friend bool operator!=(Vector3::Arg a, Vector3::Arg b) {
+		return a.x != b.x || a.y != b.y || a.z != b.z;
+	}
+
+	union {
+#ifdef _MSC_VER
+#pragma warning(push)
+#pragma warning(disable : 4201)
+#endif
+		struct
+		{
+			float x, y, z;
+		};
+#ifdef _MSC_VER
+#pragma warning(pop)
+#endif
+
+		float component[3];
+	};
+};
+
+Vector3 add(Vector3::Arg a, Vector3::Arg b) {
+	return Vector3(a.x + b.x, a.y + b.y, a.z + b.z);
+}
+Vector3 add(Vector3::Arg a, float b) {
+	return Vector3(a.x + b, a.y + b, a.z + b);
+}
+Vector3 operator+(Vector3::Arg a, Vector3::Arg b) {
+	return add(a, b);
+}
+Vector3 operator+(Vector3::Arg a, float b) {
+	return add(a, b);
+}
+
+Vector3 sub(Vector3::Arg a, Vector3::Arg b) {
+	return Vector3(a.x - b.x, a.y - b.y, a.z - b.z);
+}
+
+Vector3 sub(Vector3::Arg a, float b) {
+	return Vector3(a.x - b, a.y - b, a.z - b);
+}
+
+Vector3 operator-(Vector3::Arg a, Vector3::Arg b) {
+	return sub(a, b);
+}
+
+Vector3 operator-(Vector3::Arg a, float b) {
+	return sub(a, b);
+}
+
+Vector3 cross(Vector3::Arg a, Vector3::Arg b) {
+	return Vector3(a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x);
+}
+
+Vector3 operator*(Vector3::Arg v, float s) {
+	return Vector3(v.x * s, v.y * s, v.z * s);
+}
+
+Vector3 operator*(float s, Vector3::Arg v) {
+	return Vector3(v.x * s, v.y * s, v.z * s);
+}
+
+Vector3 operator*(Vector3::Arg v, Vector3::Arg s) {
+	return Vector3(v.x * s.x, v.y * s.y, v.z * s.z);
+}
+
+Vector3 operator/(Vector3::Arg v, float s) {
+	return v * (1.0f / s);
+}
+
+Vector3 lerp(Vector3::Arg v1, Vector3::Arg v2, float t) {
+	const float s = 1.0f - t;
+	return Vector3(v1.x * s + t * v2.x, v1.y * s + t * v2.y, v1.z * s + t * v2.z);
+}
+
+float dot(Vector3::Arg a, Vector3::Arg b) {
+	return a.x * b.x + a.y * b.y + a.z * b.z;
+}
+
+float lengthSquared(Vector3::Arg v) {
+	return v.x * v.x + v.y * v.y + v.z * v.z;
+}
+
+float length(Vector3::Arg v) {
+	return sqrtf(lengthSquared(v));
+}
+
+float distance(Vector3::Arg a, Vector3::Arg b) {
+	return length(a - b);
+}
+
+float distanceSquared(Vector3::Arg a, Vector3::Arg b) {
+	return lengthSquared(a - b);
+}
+
+bool isNormalized(Vector3::Arg v, float epsilon = NV_NORMAL_EPSILON) {
+	return equal(length(v), 1, epsilon);
+}
+
+Vector3 normalize(Vector3::Arg v, float epsilon = NV_EPSILON) {
+	float l = length(v);
+	xaDebugAssert(!isZero(l, epsilon));
+#ifdef NDEBUG
+	epsilon = 0; // silence unused parameter warning
+#endif
+	Vector3 n = v * (1.0f / l);
+	xaDebugAssert(isNormalized(n));
+	return n;
+}
+
+Vector3 normalizeSafe(Vector3::Arg v, Vector3::Arg fallback, float epsilon = NV_EPSILON) {
+	float l = length(v);
+	if (isZero(l, epsilon)) {
+		return fallback;
+	}
+	return v * (1.0f / l);
+}
+
+bool equal(Vector3::Arg v1, Vector3::Arg v2, float epsilon = NV_EPSILON) {
+	return equal(v1.x, v2.x, epsilon) && equal(v1.y, v2.y, epsilon) && equal(v1.z, v2.z, epsilon);
+}
+
+Vector3 min(Vector3::Arg a, Vector3::Arg b) {
+	return Vector3(std::min(a.x, b.x), std::min(a.y, b.y), std::min(a.z, b.z));
+}
+
+Vector3 max(Vector3::Arg a, Vector3::Arg b) {
+	return Vector3(std::max(a.x, b.x), std::max(a.y, b.y), std::max(a.z, b.z));
+}
+
+Vector3 clamp(Vector3::Arg v, float min, float max) {
+	return Vector3(clamp(v.x, min, max), clamp(v.y, min, max), clamp(v.z, min, max));
+}
+
+Vector3 saturate(Vector3::Arg v) {
+	return Vector3(saturate(v.x), saturate(v.y), saturate(v.z));
+}
+
+Vector3 floor(Vector3::Arg v) {
+	return Vector3(floorf(v.x), floorf(v.y), floorf(v.z));
+}
+
+bool isFinite(Vector3::Arg v) {
+	return std::isfinite(v.x) && std::isfinite(v.y) && std::isfinite(v.z);
+}
+
+static uint32_t hash(const Vector3 &v, uint32_t h) {
+	return sdbmFloatHash(v.component, 3, h);
+}
+
+/// Basis class to compute tangent space basis, ortogonalizations and to
+/// transform vectors from one space to another.
+class Basis {
+public:
+	/// Create a null basis.
+	Basis() :
+			tangent(0, 0, 0),
+			bitangent(0, 0, 0),
+			normal(0, 0, 0) {}
+
+	void buildFrameForDirection(Vector3::Arg d, float angle = 0) {
+		xaAssert(isNormalized(d));
+		normal = d;
+		// Choose minimum axis.
+		if (fabsf(normal.x) < fabsf(normal.y) && fabsf(normal.x) < fabsf(normal.z)) {
+			tangent = Vector3(1, 0, 0);
+		} else if (fabsf(normal.y) < fabsf(normal.z)) {
+			tangent = Vector3(0, 1, 0);
+		} else {
+			tangent = Vector3(0, 0, 1);
+		}
+		// Ortogonalize
+		tangent -= normal * dot(normal, tangent);
+		tangent = normalize(tangent);
+		bitangent = cross(normal, tangent);
+		// Rotate frame around normal according to angle.
+		if (angle != 0.0f) {
+			float c = cosf(angle);
+			float s = sinf(angle);
+			Vector3 tmp = c * tangent - s * bitangent;
+			bitangent = s * tangent + c * bitangent;
+			tangent = tmp;
+		}
+	}
+
+	Vector3 tangent;
+	Vector3 bitangent;
+	Vector3 normal;
+};
+
+// Simple bit array.
+class BitArray {
+public:
+	BitArray() :
+			m_size(0) {}
+	BitArray(uint32_t sz) {
+		resize(sz);
+	}
+
+	uint32_t size() const {
+		return m_size;
+	}
+
+	void clear() {
+		resize(0);
+	}
+
+	void resize(uint32_t new_size) {
+		m_size = new_size;
+		m_wordArray.resize((m_size + 31) >> 5);
+	}
+
+	/// Get bit.
+	bool bitAt(uint32_t b) const {
+		xaDebugAssert(b < m_size);
+		return (m_wordArray[b >> 5] & (1 << (b & 31))) != 0;
+	}
+
+	// Set a bit.
+	void setBitAt(uint32_t idx) {
+		xaDebugAssert(idx < m_size);
+		m_wordArray[idx >> 5] |= (1 << (idx & 31));
+	}
+
+	// Toggle a bit.
+	void toggleBitAt(uint32_t idx) {
+		xaDebugAssert(idx < m_size);
+		m_wordArray[idx >> 5] ^= (1 << (idx & 31));
+	}
+
+	// Set a bit to the given value. @@ Rename modifyBitAt?
+	void setBitAt(uint32_t idx, bool b) {
+		xaDebugAssert(idx < m_size);
+		m_wordArray[idx >> 5] = setBits(m_wordArray[idx >> 5], 1 << (idx & 31), b);
+		xaDebugAssert(bitAt(idx) == b);
+	}
+
+	// Clear all the bits.
+	void clearAll() {
+		memset(m_wordArray.data(), 0, m_wordArray.size() * sizeof(uint32_t));
+	}
+
+	// Set all the bits.
+	void setAll() {
+		memset(m_wordArray.data(), 0xFF, m_wordArray.size() * sizeof(uint32_t));
+	}
+
+private:
+	// See "Conditionally set or clear bits without branching" at http://graphics.stanford.edu/~seander/bithacks.html
+	uint32_t setBits(uint32_t w, uint32_t m, bool b) {
+		return (w & ~m) | (-int(b) & m);
+	}
+
+	// Number of bits stored.
+	uint32_t m_size;
+
+	// Array of bits.
+	std::vector<uint32_t> m_wordArray;
+};
+
+/// Bit map. This should probably be called BitImage.
+class BitMap {
+public:
+	BitMap() :
+			m_width(0),
+			m_height(0) {}
+	BitMap(uint32_t w, uint32_t h) :
+			m_width(w),
+			m_height(h),
+			m_bitArray(w * h) {}
+
+	uint32_t width() const {
+		return m_width;
+	}
+	uint32_t height() const {
+		return m_height;
+	}
+
+	void resize(uint32_t w, uint32_t h, bool initValue) {
+		BitArray tmp(w * h);
+		if (initValue)
+			tmp.setAll();
+		else
+			tmp.clearAll();
+		// @@ Copying one bit at a time. This could be much faster.
+		for (uint32_t y = 0; y < m_height; y++) {
+			for (uint32_t x = 0; x < m_width; x++) {
+				//tmp.setBitAt(y*w + x, bitAt(x, y));
+				if (bitAt(x, y) != initValue) tmp.toggleBitAt(y * w + x);
+			}
+		}
+		std::swap(m_bitArray, tmp);
+		m_width = w;
+		m_height = h;
+	}
+
+	bool bitAt(uint32_t x, uint32_t y) const {
+		xaDebugAssert(x < m_width && y < m_height);
+		return m_bitArray.bitAt(y * m_width + x);
+	}
+
+	void setBitAt(uint32_t x, uint32_t y) {
+		xaDebugAssert(x < m_width && y < m_height);
+		m_bitArray.setBitAt(y * m_width + x);
+	}
+
+	void clearAll() {
+		m_bitArray.clearAll();
+	}
+
+private:
+	uint32_t m_width;
+	uint32_t m_height;
+	BitArray m_bitArray;
+};
+
+// Axis Aligned Bounding Box.
+class Box {
+public:
+	Box() {}
+	Box(const Box &b) :
+			minCorner(b.minCorner),
+			maxCorner(b.maxCorner) {}
+	Box(const Vector3 &mins, const Vector3 &maxs) :
+			minCorner(mins),
+			maxCorner(maxs) {}
+
+	operator const float *() const {
+		return reinterpret_cast<const float *>(this);
+	}
+
+	// Clear the bounds.
+	void clearBounds() {
+		minCorner.set(FLT_MAX, FLT_MAX, FLT_MAX);
+		maxCorner.set(-FLT_MAX, -FLT_MAX, -FLT_MAX);
+	}
+
+	// Return extents of the box.
+	Vector3 extents() const {
+		return (maxCorner - minCorner) * 0.5f;
+	}
+
+	// Add a point to this box.
+	void addPointToBounds(const Vector3 &p) {
+		minCorner = min(minCorner, p);
+		maxCorner = max(maxCorner, p);
+	}
+
+	// Get the volume of the box.
+	float volume() const {
+		Vector3 d = extents();
+		return 8.0f * (d.x * d.y * d.z);
+	}
+
+	Vector3 minCorner;
+	Vector3 maxCorner;
+};
+
+class Fit {
+public:
+	static Vector3 computeCentroid(int n, const Vector3 *__restrict points) {
+		Vector3 centroid(0.0f);
+		for (int i = 0; i < n; i++) {
+			centroid += points[i];
+		}
+		centroid /= float(n);
+		return centroid;
+	}
+
+	static Vector3 computeCovariance(int n, const Vector3 *__restrict points, float *__restrict covariance) {
+		// compute the centroid
+		Vector3 centroid = computeCentroid(n, points);
+		// compute covariance matrix
+		for (int i = 0; i < 6; i++) {
+			covariance[i] = 0.0f;
+		}
+		for (int i = 0; i < n; i++) {
+			Vector3 v = points[i] - centroid;
+			covariance[0] += v.x * v.x;
+			covariance[1] += v.x * v.y;
+			covariance[2] += v.x * v.z;
+			covariance[3] += v.y * v.y;
+			covariance[4] += v.y * v.z;
+			covariance[5] += v.z * v.z;
+		}
+		return centroid;
+	}
+
+	static bool isPlanar(int n, const Vector3 *points, float epsilon = NV_EPSILON) {
+		// compute the centroid and covariance
+		float matrix[6];
+		computeCovariance(n, points, matrix);
+		float eigenValues[3];
+		Vector3 eigenVectors[3];
+		if (!eigenSolveSymmetric3(matrix, eigenValues, eigenVectors)) {
+			return false;
+		}
+		return eigenValues[2] < epsilon;
+	}
+
+	// Tridiagonal solver from Charles Bloom.
+	// Householder transforms followed by QL decomposition.
+	// Seems to be based on the code from Numerical Recipes in C.
+	static bool eigenSolveSymmetric3(const float matrix[6], float eigenValues[3], Vector3 eigenVectors[3]) {
+		xaDebugAssert(matrix != NULL && eigenValues != NULL && eigenVectors != NULL);
+		float subd[3];
+		float diag[3];
+		float work[3][3];
+		work[0][0] = matrix[0];
+		work[0][1] = work[1][0] = matrix[1];
+		work[0][2] = work[2][0] = matrix[2];
+		work[1][1] = matrix[3];
+		work[1][2] = work[2][1] = matrix[4];
+		work[2][2] = matrix[5];
+		EigenSolver3_Tridiagonal(work, diag, subd);
+		if (!EigenSolver3_QLAlgorithm(work, diag, subd)) {
+			for (int i = 0; i < 3; i++) {
+				eigenValues[i] = 0;
+				eigenVectors[i] = Vector3(0);
+			}
+			return false;
+		}
+		for (int i = 0; i < 3; i++) {
+			eigenValues[i] = (float)diag[i];
+		}
+		// eigenvectors are the columns; make them the rows :
+		for (int i = 0; i < 3; i++) {
+			for (int j = 0; j < 3; j++) {
+				eigenVectors[j].component[i] = (float)work[i][j];
+			}
+		}
+		// shuffle to sort by singular value :
+		if (eigenValues[2] > eigenValues[0] && eigenValues[2] > eigenValues[1]) {
+			std::swap(eigenValues[0], eigenValues[2]);
+			std::swap(eigenVectors[0], eigenVectors[2]);
+		}
+		if (eigenValues[1] > eigenValues[0]) {
+			std::swap(eigenValues[0], eigenValues[1]);
+			std::swap(eigenVectors[0], eigenVectors[1]);
+		}
+		if (eigenValues[2] > eigenValues[1]) {
+			std::swap(eigenValues[1], eigenValues[2]);
+			std::swap(eigenVectors[1], eigenVectors[2]);
+		}
+		xaDebugAssert(eigenValues[0] >= eigenValues[1] && eigenValues[0] >= eigenValues[2]);
+		xaDebugAssert(eigenValues[1] >= eigenValues[2]);
+		return true;
+	}
+
+private:
+	static void EigenSolver3_Tridiagonal(float mat[3][3], float *diag, float *subd) {
+		// Householder reduction T = Q^t M Q
+		//   Input:
+		//     mat, symmetric 3x3 matrix M
+		//   Output:
+		//     mat, orthogonal matrix Q
+		//     diag, diagonal entries of T
+		//     subd, subdiagonal entries of T (T is symmetric)
+		const float epsilon = 1e-08f;
+		float a = mat[0][0];
+		float b = mat[0][1];
+		float c = mat[0][2];
+		float d = mat[1][1];
+		float e = mat[1][2];
+		float f = mat[2][2];
+		diag[0] = a;
+		subd[2] = 0.f;
+		if (fabsf(c) >= epsilon) {
+			const float ell = sqrtf(b * b + c * c);
+			b /= ell;
+			c /= ell;
+			const float q = 2 * b * e + c * (f - d);
+			diag[1] = d + c * q;
+			diag[2] = f - c * q;
+			subd[0] = ell;
+			subd[1] = e - b * q;
+			mat[0][0] = 1;
+			mat[0][1] = 0;
+			mat[0][2] = 0;
+			mat[1][0] = 0;
+			mat[1][1] = b;
+			mat[1][2] = c;
+			mat[2][0] = 0;
+			mat[2][1] = c;
+			mat[2][2] = -b;
+		} else {
+			diag[1] = d;
+			diag[2] = f;
+			subd[0] = b;
+			subd[1] = e;
+			mat[0][0] = 1;
+			mat[0][1] = 0;
+			mat[0][2] = 0;
+			mat[1][0] = 0;
+			mat[1][1] = 1;
+			mat[1][2] = 0;
+			mat[2][0] = 0;
+			mat[2][1] = 0;
+			mat[2][2] = 1;
+		}
+	}
+
+	static bool EigenSolver3_QLAlgorithm(float mat[3][3], float *diag, float *subd) {
+		// QL iteration with implicit shifting to reduce matrix from tridiagonal
+		// to diagonal
+		const int maxiter = 32;
+		for (int ell = 0; ell < 3; ell++) {
+			int iter;
+			for (iter = 0; iter < maxiter; iter++) {
+				int m;
+				for (m = ell; m <= 1; m++) {
+					float dd = fabsf(diag[m]) + fabsf(diag[m + 1]);
+					if (fabsf(subd[m]) + dd == dd)
+						break;
+				}
+				if (m == ell)
+					break;
+				float g = (diag[ell + 1] - diag[ell]) / (2 * subd[ell]);
+				float r = sqrtf(g * g + 1);
+				if (g < 0)
+					g = diag[m] - diag[ell] + subd[ell] / (g - r);
+				else
+					g = diag[m] - diag[ell] + subd[ell] / (g + r);
+				float s = 1, c = 1, p = 0;
+				for (int i = m - 1; i >= ell; i--) {
+					float f = s * subd[i], b = c * subd[i];
+					if (fabsf(f) >= fabsf(g)) {
+						c = g / f;
+						r = sqrtf(c * c + 1);
+						subd[i + 1] = f * r;
+						c *= (s = 1 / r);
+					} else {
+						s = f / g;
+						r = sqrtf(s * s + 1);
+						subd[i + 1] = g * r;
+						s *= (c = 1 / r);
+					}
+					g = diag[i + 1] - p;
+					r = (diag[i] - g) * s + 2 * b * c;
+					p = s * r;
+					diag[i + 1] = g + p;
+					g = c * r - b;
+					for (int k = 0; k < 3; k++) {
+						f = mat[k][i + 1];
+						mat[k][i + 1] = s * mat[k][i] + c * f;
+						mat[k][i] = c * mat[k][i] - s * f;
+					}
+				}
+				diag[ell] -= p;
+				subd[ell] = g;
+				subd[m] = 0;
+			}
+			if (iter == maxiter)
+				// should not get here under normal circumstances
+				return false;
+		}
+		return true;
+	}
+};
+
+/// Fixed size vector class.
+class FullVector {
+public:
+	FullVector(uint32_t dim) { m_array.resize(dim); }
+	FullVector(const FullVector &v) :
+			m_array(v.m_array) {}
+
+	const FullVector &operator=(const FullVector &v) {
+		xaAssert(dimension() == v.dimension());
+		m_array = v.m_array;
+		return *this;
+	}
+
+	uint32_t dimension() const { return m_array.size(); }
+	const float &operator[](uint32_t index) const { return m_array[index]; }
+	float &operator[](uint32_t index) { return m_array[index]; }
+
+	void fill(float f) {
+		const uint32_t dim = dimension();
+		for (uint32_t i = 0; i < dim; i++) {
+			m_array[i] = f;
+		}
+	}
+
+	void operator+=(const FullVector &v) {
+		xaDebugAssert(dimension() == v.dimension());
+		const uint32_t dim = dimension();
+		for (uint32_t i = 0; i < dim; i++) {
+			m_array[i] += v.m_array[i];
+		}
+	}
+
+	void operator-=(const FullVector &v) {
+		xaDebugAssert(dimension() == v.dimension());
+		const uint32_t dim = dimension();
+		for (uint32_t i = 0; i < dim; i++) {
+			m_array[i] -= v.m_array[i];
+		}
+	}
+
+	void operator*=(const FullVector &v) {
+		xaDebugAssert(dimension() == v.dimension());
+		const uint32_t dim = dimension();
+		for (uint32_t i = 0; i < dim; i++) {
+			m_array[i] *= v.m_array[i];
+		}
+	}
+
+	void operator+=(float f) {
+		const uint32_t dim = dimension();
+		for (uint32_t i = 0; i < dim; i++) {
+			m_array[i] += f;
+		}
+	}
+
+	void operator-=(float f) {
+		const uint32_t dim = dimension();
+		for (uint32_t i = 0; i < dim; i++) {
+			m_array[i] -= f;
+		}
+	}
+
+	void operator*=(float f) {
+		const uint32_t dim = dimension();
+		for (uint32_t i = 0; i < dim; i++) {
+			m_array[i] *= f;
+		}
+	}
+
+private:
+	std::vector<float> m_array;
+};
+
+namespace halfedge {
+class Face;
+class Vertex;
+
+class Edge {
+public:
+	uint32_t id;
+	Edge *next;
+	Edge *prev; // This is not strictly half-edge, but makes algorithms easier and faster.
+	Edge *pair;
+	Vertex *vertex;
+	Face *face;
+
+	// Default constructor.
+	Edge(uint32_t id) :
+			id(id),
+			next(NULL),
+			prev(NULL),
+			pair(NULL),
+			vertex(NULL),
+			face(NULL) {}
+
+	// Vertex queries.
+	const Vertex *from() const {
+		return vertex;
+	}
+
+	Vertex *from() {
+		return vertex;
+	}
+
+	const Vertex *to() const {
+		return pair->vertex; // This used to be 'next->vertex', but that changed often when the connectivity of the mesh changes.
+	}
+
+	Vertex *to() {
+		return pair->vertex;
+	}
+
+	// Edge queries.
+	void setNext(Edge *e) {
+		next = e;
+		if (e != NULL) e->prev = this;
+	}
+	void setPrev(Edge *e) {
+		prev = e;
+		if (e != NULL) e->next = this;
+	}
+
+	// @@ It would be more simple to only check m_pair == NULL
+	// Face queries.
+	bool isBoundary() const {
+		return !(face && pair->face);
+	}
+
+	// @@ This is not exactly accurate, we should compare the texture coordinates...
+	bool isSeam() const {
+		return vertex != pair->next->vertex || next->vertex != pair->vertex;
+	}
+
+	bool isNormalSeam() const;
+	bool isTextureSeam() const;
+
+	bool isValid() const {
+		// null face is OK.
+		if (next == NULL || prev == NULL || pair == NULL || vertex == NULL) return false;
+		if (next->prev != this) return false;
+		if (prev->next != this) return false;
+		if (pair->pair != this) return false;
+		return true;
+	}
+
+	float length() const;
+
+	// Return angle between this edge and the previous one.
+	float angle() const;
+};
+
+class Vertex {
+public:
+	uint32_t id;
+	uint32_t original_id;
+	Edge *edge;
+	Vertex *next;
+	Vertex *prev;
+	Vector3 pos;
+	Vector3 nor;
+	Vector2 tex;
+
+	Vertex(uint32_t id) :
+			id(id),
+			original_id(id),
+			edge(NULL),
+			pos(0.0f),
+			nor(0.0f),
+			tex(0.0f) {
+		next = this;
+		prev = this;
+	}
+
+	// Set first edge of all colocals.
+	void setEdge(Edge *e) {
+		for (VertexIterator it(colocals()); !it.isDone(); it.advance()) {
+			it.current()->edge = e;
+		}
+	}
+
+	// Update position of all colocals.
+	void setPos(const Vector3 &p) {
+		for (VertexIterator it(colocals()); !it.isDone(); it.advance()) {
+			it.current()->pos = p;
+		}
+	}
+
+	bool isFirstColocal() const {
+		return firstColocal() == this;
+	}
+
+	const Vertex *firstColocal() const {
+		uint32_t firstId = id;
+		const Vertex *vertex = this;
+		for (ConstVertexIterator it(colocals()); !it.isDone(); it.advance()) {
+			if (it.current()->id < firstId) {
+				firstId = vertex->id;
+				vertex = it.current();
+			}
+		}
+		return vertex;
+	}
+
+	Vertex *firstColocal() {
+		Vertex *vertex = this;
+		uint32_t firstId = id;
+		for (VertexIterator it(colocals()); !it.isDone(); it.advance()) {
+			if (it.current()->id < firstId) {
+				firstId = vertex->id;
+				vertex = it.current();
+			}
+		}
+		return vertex;
+	}
+
+	bool isColocal(const Vertex *v) const {
+		if (this == v) return true;
+		if (pos != v->pos) return false;
+		for (ConstVertexIterator it(colocals()); !it.isDone(); it.advance()) {
+			if (v == it.current()) {
+				return true;
+			}
+		}
+		return false;
+	}
+
+	void linkColocal(Vertex *v) {
+		next->prev = v;
+		v->next = next;
+		next = v;
+		v->prev = this;
+	}
+	void unlinkColocal() {
+		next->prev = prev;
+		prev->next = next;
+		next = this;
+		prev = this;
+	}
+
+	// @@ Note: This only works if linkBoundary has been called.
+	bool isBoundary() const {
+		return (edge && !edge->face);
+	}
+
+	// Iterator that visits the edges around this vertex in counterclockwise order.
+	class EdgeIterator //: public Iterator<Edge *>
+	{
+	public:
+		EdgeIterator(Edge *e) :
+				m_end(NULL),
+				m_current(e) {}
+
+		virtual void advance() {
+			if (m_end == NULL) m_end = m_current;
+			m_current = m_current->pair->next;
+			//m_current = m_current->prev->pair;
+		}
+
+		virtual bool isDone() const {
+			return m_end == m_current;
+		}
+		virtual Edge *current() const {
+			return m_current;
+		}
+		Vertex *vertex() const {
+			return m_current->vertex;
+		}
+
+	private:
+		Edge *m_end;
+		Edge *m_current;
+	};
+
+	EdgeIterator edges() {
+		return EdgeIterator(edge);
+	}
+	EdgeIterator edges(Edge *e) {
+		return EdgeIterator(e);
+	}
+
+	// Iterator that visits the edges around this vertex in counterclockwise order.
+	class ConstEdgeIterator //: public Iterator<Edge *>
+	{
+	public:
+		ConstEdgeIterator(const Edge *e) :
+				m_end(NULL),
+				m_current(e) {}
+		ConstEdgeIterator(EdgeIterator it) :
+				m_end(NULL),
+				m_current(it.current()) {}
+
+		virtual void advance() {
+			if (m_end == NULL) m_end = m_current;
+			m_current = m_current->pair->next;
+			//m_current = m_current->prev->pair;
+		}
+
+		virtual bool isDone() const {
+			return m_end == m_current;
+		}
+		virtual const Edge *current() const {
+			return m_current;
+		}
+		const Vertex *vertex() const {
+			return m_current->to();
+		}
+
+	private:
+		const Edge *m_end;
+		const Edge *m_current;
+	};
+
+	ConstEdgeIterator edges() const {
+		return ConstEdgeIterator(edge);
+	}
+	ConstEdgeIterator edges(const Edge *e) const {
+		return ConstEdgeIterator(e);
+	}
+
+	// Iterator that visits all the colocal vertices.
+	class VertexIterator //: public Iterator<Edge *>
+	{
+	public:
+		VertexIterator(Vertex *v) :
+				m_end(NULL),
+				m_current(v) {}
+
+		virtual void advance() {
+			if (m_end == NULL) m_end = m_current;
+			m_current = m_current->next;
+		}
+
+		virtual bool isDone() const {
+			return m_end == m_current;
+		}
+		virtual Vertex *current() const {
+			return m_current;
+		}
+
+	private:
+		Vertex *m_end;
+		Vertex *m_current;
+	};
+
+	VertexIterator colocals() {
+		return VertexIterator(this);
+	}
+
+	// Iterator that visits all the colocal vertices.
+	class ConstVertexIterator //: public Iterator<Edge *>
+	{
+	public:
+		ConstVertexIterator(const Vertex *v) :
+				m_end(NULL),
+				m_current(v) {}
+
+		virtual void advance() {
+			if (m_end == NULL) m_end = m_current;
+			m_current = m_current->next;
+		}
+
+		virtual bool isDone() const {
+			return m_end == m_current;
+		}
+		virtual const Vertex *current() const {
+			return m_current;
+		}
+
+	private:
+		const Vertex *m_end;
+		const Vertex *m_current;
+	};
+
+	ConstVertexIterator colocals() const {
+		return ConstVertexIterator(this);
+	}
+};
+
+bool Edge::isNormalSeam() const {
+	return (vertex->nor != pair->next->vertex->nor || next->vertex->nor != pair->vertex->nor);
+}
+
+bool Edge::isTextureSeam() const {
+	return (vertex->tex != pair->next->vertex->tex || next->vertex->tex != pair->vertex->tex);
+}
+
+float Edge::length() const {
+	return internal::length(to()->pos - from()->pos);
+}
+
+float Edge::angle() const {
+	Vector3 p = vertex->pos;
+	Vector3 a = prev->vertex->pos;
+	Vector3 b = next->vertex->pos;
+	Vector3 v0 = a - p;
+	Vector3 v1 = b - p;
+	return acosf(dot(v0, v1) / (internal::length(v0) * internal::length(v1)));
+}
+
+class Face {
+public:
+	uint32_t id;
+	uint16_t group;
+	uint16_t material;
+	Edge *edge;
+
+	Face(uint32_t id) :
+			id(id),
+			group(uint16_t(~0)),
+			material(uint16_t(~0)),
+			edge(NULL) {}
+
+	float area() const {
+		float area = 0;
+		const Vector3 &v0 = edge->from()->pos;
+		for (ConstEdgeIterator it(edges(edge->next)); it.current() != edge->prev; it.advance()) {
+			const Edge *e = it.current();
+			const Vector3 &v1 = e->vertex->pos;
+			const Vector3 &v2 = e->next->vertex->pos;
+			area += length(cross(v1 - v0, v2 - v0));
+		}
+		return area * 0.5f;
+	}
+
+	float parametricArea() const {
+		float area = 0;
+		const Vector2 &v0 = edge->from()->tex;
+		for (ConstEdgeIterator it(edges(edge->next)); it.current() != edge->prev; it.advance()) {
+			const Edge *e = it.current();
+			const Vector2 &v1 = e->vertex->tex;
+			const Vector2 &v2 = e->next->vertex->tex;
+			area += triangleArea(v0, v1, v2);
+		}
+		return area * 0.5f;
+	}
+
+	Vector3 normal() const {
+		Vector3 n(0);
+		const Vertex *vertex0 = NULL;
+		for (ConstEdgeIterator it(edges()); !it.isDone(); it.advance()) {
+			const Edge *e = it.current();
+			xaAssert(e != NULL);
+			if (vertex0 == NULL) {
+				vertex0 = e->vertex;
+			} else if (e->next->vertex != vertex0) {
+				const halfedge::Vertex *vertex1 = e->from();
+				const halfedge::Vertex *vertex2 = e->to();
+				const Vector3 &p0 = vertex0->pos;
+				const Vector3 &p1 = vertex1->pos;
+				const Vector3 &p2 = vertex2->pos;
+				Vector3 v10 = p1 - p0;
+				Vector3 v20 = p2 - p0;
+				n += cross(v10, v20);
+			}
+		}
+		return normalizeSafe(n, Vector3(0, 0, 1), 0.0f);
+	}
+
+	Vector3 centroid() const {
+		Vector3 sum(0.0f);
+		uint32_t count = 0;
+		for (ConstEdgeIterator it(edges()); !it.isDone(); it.advance()) {
+			const Edge *e = it.current();
+			sum += e->from()->pos;
+			count++;
+		}
+		return sum / float(count);
+	}
+
+	// Unnormalized face normal assuming it's a triangle.
+	Vector3 triangleNormal() const {
+		Vector3 p0 = edge->vertex->pos;
+		Vector3 p1 = edge->next->vertex->pos;
+		Vector3 p2 = edge->next->next->vertex->pos;
+		Vector3 e0 = p2 - p0;
+		Vector3 e1 = p1 - p0;
+		return normalizeSafe(cross(e0, e1), Vector3(0), 0.0f);
+	}
+
+	Vector3 triangleNormalAreaScaled() const {
+		Vector3 p0 = edge->vertex->pos;
+		Vector3 p1 = edge->next->vertex->pos;
+		Vector3 p2 = edge->next->next->vertex->pos;
+		Vector3 e0 = p2 - p0;
+		Vector3 e1 = p1 - p0;
+		return cross(e0, e1);
+	}
+
+	// Average of the edge midpoints weighted by the edge length.
+	// I want a point inside the triangle, but closer to the cirumcenter.
+	Vector3 triangleCenter() const {
+		Vector3 p0 = edge->vertex->pos;
+		Vector3 p1 = edge->next->vertex->pos;
+		Vector3 p2 = edge->next->next->vertex->pos;
+		float l0 = length(p1 - p0);
+		float l1 = length(p2 - p1);
+		float l2 = length(p0 - p2);
+		Vector3 m0 = (p0 + p1) * l0 / (l0 + l1 + l2);
+		Vector3 m1 = (p1 + p2) * l1 / (l0 + l1 + l2);
+		Vector3 m2 = (p2 + p0) * l2 / (l0 + l1 + l2);
+		return m0 + m1 + m2;
+	}
+
+	bool isValid() const {
+		uint32_t count = 0;
+		for (ConstEdgeIterator it(edges()); !it.isDone(); it.advance()) {
+			const Edge *e = it.current();
+			if (e->face != this) return false;
+			if (!e->isValid()) return false;
+			if (!e->pair->isValid()) return false;
+			count++;
+		}
+		if (count < 3) return false;
+		return true;
+	}
+
+	bool contains(const Edge *e) const {
+		for (ConstEdgeIterator it(edges()); !it.isDone(); it.advance()) {
+			if (it.current() == e) return true;
+		}
+		return false;
+	}
+
+	uint32_t edgeCount() const {
+		uint32_t count = 0;
+		for (ConstEdgeIterator it(edges()); !it.isDone(); it.advance()) {
+			++count;
+		}
+		return count;
+	}
+
+	// The iterator that visits the edges of this face in clockwise order.
+	class EdgeIterator //: public Iterator<Edge *>
+	{
+	public:
+		EdgeIterator(Edge *e) :
+				m_end(NULL),
+				m_current(e) {}
+
+		virtual void advance() {
+			if (m_end == NULL) m_end = m_current;
+			m_current = m_current->next;
+		}
+
+		virtual bool isDone() const {
+			return m_end == m_current;
+		}
+		virtual Edge *current() const {
+			return m_current;
+		}
+		Vertex *vertex() const {
+			return m_current->vertex;
+		}
+
+	private:
+		Edge *m_end;
+		Edge *m_current;
+	};
+
+	EdgeIterator edges() {
+		return EdgeIterator(edge);
+	}
+	EdgeIterator edges(Edge *e) {
+		xaDebugAssert(contains(e));
+		return EdgeIterator(e);
+	}
+
+	// The iterator that visits the edges of this face in clockwise order.
+	class ConstEdgeIterator //: public Iterator<const Edge *>
+	{
+	public:
+		ConstEdgeIterator(const Edge *e) :
+				m_end(NULL),
+				m_current(e) {}
+		ConstEdgeIterator(const EdgeIterator &it) :
+				m_end(NULL),
+				m_current(it.current()) {}
+
+		virtual void advance() {
+			if (m_end == NULL) m_end = m_current;
+			m_current = m_current->next;
+		}
+
+		virtual bool isDone() const {
+			return m_end == m_current;
+		}
+		virtual const Edge *current() const {
+			return m_current;
+		}
+		const Vertex *vertex() const {
+			return m_current->vertex;
+		}
+
+	private:
+		const Edge *m_end;
+		const Edge *m_current;
+	};
+
+	ConstEdgeIterator edges() const {
+		return ConstEdgeIterator(edge);
+	}
+	ConstEdgeIterator edges(const Edge *e) const {
+		xaDebugAssert(contains(e));
+		return ConstEdgeIterator(e);
+	}
+};
+
+/// Simple half edge mesh designed for dynamic mesh manipulation.
+class Mesh {
+public:
+	Mesh() :
+			m_colocalVertexCount(0) {}
+
+	Mesh(const Mesh *mesh) {
+		// Copy mesh vertices.
+		const uint32_t vertexCount = mesh->vertexCount();
+		m_vertexArray.resize(vertexCount);
+		for (uint32_t v = 0; v < vertexCount; v++) {
+			const Vertex *vertex = mesh->vertexAt(v);
+			xaDebugAssert(vertex->id == v);
+			m_vertexArray[v] = new Vertex(v);
+			m_vertexArray[v]->pos = vertex->pos;
+			m_vertexArray[v]->nor = vertex->nor;
+			m_vertexArray[v]->tex = vertex->tex;
+		}
+		m_colocalVertexCount = vertexCount;
+		// Copy mesh faces.
+		const uint32_t faceCount = mesh->faceCount();
+		std::vector<uint32_t> indexArray;
+		indexArray.reserve(3);
+		for (uint32_t f = 0; f < faceCount; f++) {
+			const Face *face = mesh->faceAt(f);
+			for (Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) {
+				const Vertex *vertex = it.current()->from();
+				indexArray.push_back(vertex->id);
+			}
+			addFace(indexArray);
+			indexArray.clear();
+		}
+	}
+
+	~Mesh() {
+		clear();
+	}
+
+	void clear() {
+		for (size_t i = 0; i < m_vertexArray.size(); i++)
+			delete m_vertexArray[i];
+		m_vertexArray.clear();
+		for (auto it = m_edgeMap.begin(); it != m_edgeMap.end(); it++)
+			delete it->second;
+		m_edgeArray.clear();
+		m_edgeMap.clear();
+		for (size_t i = 0; i < m_faceArray.size(); i++)
+			delete m_faceArray[i];
+		m_faceArray.clear();
+	}
+
+	Vertex *addVertex(const Vector3 &pos) {
+		xaDebugAssert(isFinite(pos));
+		Vertex *v = new Vertex(m_vertexArray.size());
+		v->pos = pos;
+		m_vertexArray.push_back(v);
+		return v;
+	}
+
+	/// Link colocal vertices based on geometric location only.
+	void linkColocals() {
+		xaPrint("--- Linking colocals:\n");
+		const uint32_t vertexCount = this->vertexCount();
+		std::unordered_map<Vector3, Vertex *, Hash<Vector3>, Equal<Vector3> > vertexMap;
+		vertexMap.reserve(vertexCount);
+		for (uint32_t v = 0; v < vertexCount; v++) {
+			Vertex *vertex = vertexAt(v);
+			Vertex *colocal = vertexMap[vertex->pos];
+			if (colocal) {
+				colocal->linkColocal(vertex);
+			} else {
+				vertexMap[vertex->pos] = vertex;
+			}
+		}
+		m_colocalVertexCount = vertexMap.size();
+		xaPrint("---   %d vertex positions.\n", m_colocalVertexCount);
+		// @@ Remove duplicated vertices? or just leave them as colocals?
+	}
+
+	void linkColocalsWithCanonicalMap(const std::vector<uint32_t> &canonicalMap) {
+		xaPrint("--- Linking colocals:\n");
+		uint32_t vertexMapSize = 0;
+		for (uint32_t i = 0; i < canonicalMap.size(); i++) {
+			vertexMapSize = std::max(vertexMapSize, canonicalMap[i] + 1);
+		}
+		std::vector<Vertex *> vertexMap;
+		vertexMap.resize(vertexMapSize, NULL);
+		m_colocalVertexCount = 0;
+		const uint32_t vertexCount = this->vertexCount();
+		for (uint32_t v = 0; v < vertexCount; v++) {
+			Vertex *vertex = vertexAt(v);
+			Vertex *colocal = vertexMap[canonicalMap[v]];
+			if (colocal != NULL) {
+				xaDebugAssert(vertex->pos == colocal->pos);
+				colocal->linkColocal(vertex);
+			} else {
+				vertexMap[canonicalMap[v]] = vertex;
+				m_colocalVertexCount++;
+			}
+		}
+		xaPrint("---   %d vertex positions.\n", m_colocalVertexCount);
+	}
+
+	Face *addFace() {
+		Face *f = new Face(m_faceArray.size());
+		m_faceArray.push_back(f);
+		return f;
+	}
+
+	Face *addFace(uint32_t v0, uint32_t v1, uint32_t v2) {
+		uint32_t indexArray[3];
+		indexArray[0] = v0;
+		indexArray[1] = v1;
+		indexArray[2] = v2;
+		return addFace(indexArray, 3, 0, 3);
+	}
+
+	Face *addUniqueFace(uint32_t v0, uint32_t v1, uint32_t v2) {
+
+		int base_vertex = m_vertexArray.size();
+
+		uint32_t ids[3] = { v0, v1, v2 };
+
+		Vector3 base[3] = {
+			m_vertexArray[v0]->pos,
+			m_vertexArray[v1]->pos,
+			m_vertexArray[v2]->pos,
+		};
+
+		//make sure its not a degenerate
+		bool degenerate = distanceSquared(base[0], base[1]) < NV_EPSILON || distanceSquared(base[0], base[2]) < NV_EPSILON || distanceSquared(base[1], base[2]) < NV_EPSILON;
+		xaDebugAssert(!degenerate);
+
+		float min_x = 0;
+
+		for (int i = 0; i < 3; i++) {
+			if (i == 0 || m_vertexArray[v0]->pos.x < min_x) {
+				min_x = m_vertexArray[v0]->pos.x;
+			}
+		}
+
+		float max_x = 0;
+
+		for (int j = 0; j < m_vertexArray.size(); j++) {
+			if (j == 0 || m_vertexArray[j]->pos.x > max_x) { //vertex already exists
+				max_x = m_vertexArray[j]->pos.x;
+			}
+		}
+
+		//separate from everything else, in x axis
+		for (int i = 0; i < 3; i++) {
+
+			base[i].x -= min_x;
+			base[i].x += max_x + 10.0;
+		}
+
+		for (int i = 0; i < 3; i++) {
+			Vertex *v = new Vertex(m_vertexArray.size());
+			v->pos = base[i];
+			v->nor = m_vertexArray[ids[i]]->nor,
+			v->tex = m_vertexArray[ids[i]]->tex,
+
+			v->original_id = ids[i];
+			m_vertexArray.push_back(v);
+		}
+
+		uint32_t indexArray[3];
+		indexArray[0] = base_vertex + 0;
+		indexArray[1] = base_vertex + 1;
+		indexArray[2] = base_vertex + 2;
+		return addFace(indexArray, 3, 0, 3);
+	}
+
+	Face *addFace(uint32_t v0, uint32_t v1, uint32_t v2, uint32_t v3) {
+		uint32_t indexArray[4];
+		indexArray[0] = v0;
+		indexArray[1] = v1;
+		indexArray[2] = v2;
+		indexArray[3] = v3;
+		return addFace(indexArray, 4, 0, 4);
+	}
+
+	Face *addFace(const std::vector<uint32_t> &indexArray) {
+		return addFace(indexArray, 0, indexArray.size());
+	}
+
+	Face *addFace(const std::vector<uint32_t> &indexArray, uint32_t first, uint32_t num) {
+		return addFace(indexArray.data(), (uint32_t)indexArray.size(), first, num);
+	}
+
+	Face *addFace(const uint32_t *indexArray, uint32_t indexCount, uint32_t first, uint32_t num) {
+		xaDebugAssert(first < indexCount);
+		xaDebugAssert(num <= indexCount - first);
+		xaDebugAssert(num > 2);
+		if (!canAddFace(indexArray, first, num)) {
+			return NULL;
+		}
+		Face *f = new Face(m_faceArray.size());
+		Edge *firstEdge = NULL;
+		Edge *last = NULL;
+		Edge *current = NULL;
+		for (uint32_t i = 0; i < num - 1; i++) {
+			current = addEdge(indexArray[first + i], indexArray[first + i + 1]);
+			xaAssert(current != NULL && current->face == NULL);
+			current->face = f;
+			if (last != NULL)
+				last->setNext(current);
+			else
+				firstEdge = current;
+			last = current;
+		}
+		current = addEdge(indexArray[first + num - 1], indexArray[first]);
+		xaAssert(current != NULL && current->face == NULL);
+		current->face = f;
+		last->setNext(current);
+		current->setNext(firstEdge);
+		f->edge = firstEdge;
+		m_faceArray.push_back(f);
+		return f;
+	}
+
+	// These functions disconnect the given element from the mesh and delete it.
+
+	// @@ We must always disconnect edge pairs simultaneously.
+	void disconnect(Edge *edge) {
+		xaDebugAssert(edge != NULL);
+		// Remove from edge list.
+		if ((edge->id & 1) == 0) {
+			xaDebugAssert(m_edgeArray[edge->id / 2] == edge);
+			m_edgeArray[edge->id / 2] = NULL;
+		}
+		// Remove edge from map. @@ Store map key inside edge?
+		xaDebugAssert(edge->from() != NULL && edge->to() != NULL);
+		size_t removed = m_edgeMap.erase(Key(edge->from()->id, edge->to()->id));
+		xaDebugAssert(removed == 1);
+#ifdef NDEBUG
+		removed = 0; // silence unused parameter warning
+#endif
+		// Disconnect from vertex.
+		if (edge->vertex != NULL) {
+			if (edge->vertex->edge == edge) {
+				if (edge->prev && edge->prev->pair) {
+					edge->vertex->edge = edge->prev->pair;
+				} else if (edge->pair && edge->pair->next) {
+					edge->vertex->edge = edge->pair->next;
+				} else {
+					edge->vertex->edge = NULL;
+					// @@ Remove disconnected vertex?
+				}
+			}
+		}
+		// Disconnect from face.
+		if (edge->face != NULL) {
+			if (edge->face->edge == edge) {
+				if (edge->next != NULL && edge->next != edge) {
+					edge->face->edge = edge->next;
+				} else if (edge->prev != NULL && edge->prev != edge) {
+					edge->face->edge = edge->prev;
+				} else {
+					edge->face->edge = NULL;
+					// @@ Remove disconnected face?
+				}
+			}
+		}
+		// Disconnect from previous.
+		if (edge->prev) {
+			if (edge->prev->next == edge) {
+				edge->prev->setNext(NULL);
+			}
+			//edge->setPrev(NULL);
+		}
+		// Disconnect from next.
+		if (edge->next) {
+			if (edge->next->prev == edge) {
+				edge->next->setPrev(NULL);
+			}
+			//edge->setNext(NULL);
+		}
+	}
+
+	void remove(Edge *edge) {
+		xaDebugAssert(edge != NULL);
+		disconnect(edge);
+		delete edge;
+	}
+
+	void remove(Vertex *vertex) {
+		xaDebugAssert(vertex != NULL);
+		// Remove from vertex list.
+		m_vertexArray[vertex->id] = NULL;
+		// Disconnect from colocals.
+		vertex->unlinkColocal();
+		// Disconnect from edges.
+		if (vertex->edge != NULL) {
+			// @@ Removing a connected vertex is asking for trouble...
+			if (vertex->edge->vertex == vertex) {
+				// @@ Connect edge to a colocal?
+				vertex->edge->vertex = NULL;
+			}
+			vertex->setEdge(NULL);
+		}
+		delete vertex;
+	}
+
+	void remove(Face *face) {
+		xaDebugAssert(face != NULL);
+		// Remove from face list.
+		m_faceArray[face->id] = NULL;
+		// Disconnect from edges.
+		if (face->edge != NULL) {
+			xaDebugAssert(face->edge->face == face);
+			face->edge->face = NULL;
+			face->edge = NULL;
+		}
+		delete face;
+	}
+
+	// Triangulate in place.
+	void triangulate() {
+		bool all_triangles = true;
+		const uint32_t faceCount = m_faceArray.size();
+		for (uint32_t f = 0; f < faceCount; f++) {
+			Face *face = m_faceArray[f];
+			if (face->edgeCount() != 3) {
+				all_triangles = false;
+				break;
+			}
+		}
+		if (all_triangles) {
+			return;
+		}
+		// Do not touch vertices, but rebuild edges and faces.
+		std::vector<Edge *> edgeArray;
+		std::vector<Face *> faceArray;
+		std::swap(edgeArray, m_edgeArray);
+		std::swap(faceArray, m_faceArray);
+		m_edgeMap.clear();
+		for (uint32_t f = 0; f < faceCount; f++) {
+			Face *face = faceArray[f];
+			// Trivial fan-like triangulation.
+			const uint32_t v0 = face->edge->vertex->id;
+			uint32_t v2, v1 = (uint32_t)-1;
+			for (Face::EdgeIterator it(face->edges()); !it.isDone(); it.advance()) {
+				Edge *edge = it.current();
+				v2 = edge->to()->id;
+				if (v2 == v0) break;
+				if (v1 != -1) addFace(v0, v1, v2);
+				v1 = v2;
+			}
+		}
+		xaDebugAssert(m_faceArray.size() > faceCount); // triangle count > face count
+		linkBoundary();
+		for (size_t i = 0; i < edgeArray.size(); i++)
+			delete edgeArray[i];
+		for (size_t i = 0; i < faceArray.size(); i++)
+			delete faceArray[i];
+	}
+
+	/// Link boundary edges once the mesh has been created.
+	void linkBoundary() {
+		xaPrint("--- Linking boundaries:\n");
+		int num = 0;
+		// Create boundary edges.
+		uint32_t edgeCount = this->edgeCount();
+		for (uint32_t e = 0; e < edgeCount; e++) {
+			Edge *edge = edgeAt(e);
+			if (edge != NULL && edge->pair == NULL) {
+				Edge *pair = new Edge(edge->id + 1);
+				uint32_t i = edge->from()->id;
+				uint32_t j = edge->next->from()->id;
+				Key key(j, i);
+				xaAssert(m_edgeMap.find(key) == m_edgeMap.end());
+				pair->vertex = m_vertexArray[j];
+				m_edgeMap[key] = pair;
+				edge->pair = pair;
+				pair->pair = edge;
+				num++;
+			}
+		}
+		// Link boundary edges.
+		for (uint32_t e = 0; e < edgeCount; e++) {
+			Edge *edge = edgeAt(e);
+			if (edge != NULL && edge->pair->face == NULL) {
+				linkBoundaryEdge(edge->pair);
+			}
+		}
+		xaPrint("---   %d boundary edges.\n", num);
+	}
+
+	/*
+	Fixing T-junctions.
+
+	- Find T-junctions. Find  vertices that are on an edge.
+		- This test is approximate.
+		- Insert edges on a spatial index to speedup queries.
+		- Consider only open edges, that is edges that have no pairs.
+		- Consider only vertices on boundaries.
+	- Close T-junction.
+		- Split edge.
+
+	*/
+	bool splitBoundaryEdges() // Returns true if any split was made.
+	{
+		std::vector<Vertex *> boundaryVertices;
+		for (uint32_t i = 0; i < m_vertexArray.size(); i++) {
+			Vertex *v = m_vertexArray[i];
+			if (v->isBoundary()) {
+				boundaryVertices.push_back(v);
+			}
+		}
+		xaPrint("Fixing T-junctions:\n");
+		int splitCount = 0;
+		for (uint32_t v = 0; v < boundaryVertices.size(); v++) {
+			Vertex *vertex = boundaryVertices[v];
+			Vector3 x0 = vertex->pos;
+			// Find edges that this vertex overlaps with.
+			for (uint32_t e = 0; e < m_edgeArray.size(); e++) {
+				Edge *edge = m_edgeArray[e];
+				if (edge != NULL && edge->isBoundary()) {
+					if (edge->from() == vertex || edge->to() == vertex) {
+						continue;
+					}
+					Vector3 x1 = edge->from()->pos;
+					Vector3 x2 = edge->to()->pos;
+					Vector3 v01 = x0 - x1;
+					Vector3 v21 = x2 - x1;
+					float l = length(v21);
+					float d = length(cross(v01, v21)) / l;
+					if (isZero(d)) {
+						float t = dot(v01, v21) / (l * l);
+						if (t > 0.0f + NV_EPSILON && t < 1.0f - NV_EPSILON) {
+							xaDebugAssert(equal(lerp(x1, x2, t), x0));
+							Vertex *splitVertex = splitBoundaryEdge(edge, t, x0);
+							vertex->linkColocal(splitVertex); // @@ Should we do this here?
+							splitCount++;
+						}
+					}
+				}
+			}
+		}
+		xaPrint(" - %d edges split.\n", splitCount);
+		xaDebugAssert(isValid());
+		return splitCount != 0;
+	}
+
+	// Vertices
+	uint32_t vertexCount() const {
+		return m_vertexArray.size();
+	}
+	const Vertex *vertexAt(int i) const {
+		return m_vertexArray[i];
+	}
+	Vertex *vertexAt(int i) {
+		return m_vertexArray[i];
+	}
+
+	uint32_t colocalVertexCount() const {
+		return m_colocalVertexCount;
+	}
+
+	// Faces
+	uint32_t faceCount() const {
+		return m_faceArray.size();
+	}
+	const Face *faceAt(int i) const {
+		return m_faceArray[i];
+	}
+	Face *faceAt(int i) {
+		return m_faceArray[i];
+	}
+
+	// Edges
+	uint32_t edgeCount() const {
+		return m_edgeArray.size();
+	}
+	const Edge *edgeAt(int i) const {
+		return m_edgeArray[i];
+	}
+	Edge *edgeAt(int i) {
+		return m_edgeArray[i];
+	}
+
+	class ConstVertexIterator;
+
+	class VertexIterator {
+		friend class ConstVertexIterator;
+
+	public:
+		VertexIterator(Mesh *mesh) :
+				m_mesh(mesh),
+				m_current(0) {}
+
+		virtual void advance() {
+			m_current++;
+		}
+		virtual bool isDone() const {
+			return m_current == m_mesh->vertexCount();
+		}
+		virtual Vertex *current() const {
+			return m_mesh->vertexAt(m_current);
+		}
+
+	private:
+		halfedge::Mesh *m_mesh;
+		uint32_t m_current;
+	};
+	VertexIterator vertices() {
+		return VertexIterator(this);
+	}
+
+	class ConstVertexIterator {
+	public:
+		ConstVertexIterator(const Mesh *mesh) :
+				m_mesh(mesh),
+				m_current(0) {}
+		ConstVertexIterator(class VertexIterator &it) :
+				m_mesh(it.m_mesh),
+				m_current(it.m_current) {}
+
+		virtual void advance() {
+			m_current++;
+		}
+		virtual bool isDone() const {
+			return m_current == m_mesh->vertexCount();
+		}
+		virtual const Vertex *current() const {
+			return m_mesh->vertexAt(m_current);
+		}
+
+	private:
+		const halfedge::Mesh *m_mesh;
+		uint32_t m_current;
+	};
+	ConstVertexIterator vertices() const {
+		return ConstVertexIterator(this);
+	}
+
+	class ConstFaceIterator;
+
+	class FaceIterator {
+		friend class ConstFaceIterator;
+
+	public:
+		FaceIterator(Mesh *mesh) :
+				m_mesh(mesh),
+				m_current(0) {}
+
+		virtual void advance() {
+			m_current++;
+		}
+		virtual bool isDone() const {
+			return m_current == m_mesh->faceCount();
+		}
+		virtual Face *current() const {
+			return m_mesh->faceAt(m_current);
+		}
+
+	private:
+		halfedge::Mesh *m_mesh;
+		uint32_t m_current;
+	};
+	FaceIterator faces() {
+		return FaceIterator(this);
+	}
+
+	class ConstFaceIterator {
+	public:
+		ConstFaceIterator(const Mesh *mesh) :
+				m_mesh(mesh),
+				m_current(0) {}
+		ConstFaceIterator(const FaceIterator &it) :
+				m_mesh(it.m_mesh),
+				m_current(it.m_current) {}
+
+		virtual void advance() {
+			m_current++;
+		}
+		virtual bool isDone() const {
+			return m_current == m_mesh->faceCount();
+		}
+		virtual const Face *current() const {
+			return m_mesh->faceAt(m_current);
+		}
+
+	private:
+		const halfedge::Mesh *m_mesh;
+		uint32_t m_current;
+	};
+	ConstFaceIterator faces() const {
+		return ConstFaceIterator(this);
+	}
+
+	class ConstEdgeIterator;
+
+	class EdgeIterator {
+		friend class ConstEdgeIterator;
+
+	public:
+		EdgeIterator(Mesh *mesh) :
+				m_mesh(mesh),
+				m_current(0) {}
+
+		virtual void advance() {
+			m_current++;
+		}
+		virtual bool isDone() const {
+			return m_current == m_mesh->edgeCount();
+		}
+		virtual Edge *current() const {
+			return m_mesh->edgeAt(m_current);
+		}
+
+	private:
+		halfedge::Mesh *m_mesh;
+		uint32_t m_current;
+	};
+	EdgeIterator edges() {
+		return EdgeIterator(this);
+	}
+
+	class ConstEdgeIterator {
+	public:
+		ConstEdgeIterator(const Mesh *mesh) :
+				m_mesh(mesh),
+				m_current(0) {}
+		ConstEdgeIterator(const EdgeIterator &it) :
+				m_mesh(it.m_mesh),
+				m_current(it.m_current) {}
+
+		virtual void advance() {
+			m_current++;
+		}
+		virtual bool isDone() const {
+			return m_current == m_mesh->edgeCount();
+		}
+		virtual const Edge *current() const {
+			return m_mesh->edgeAt(m_current);
+		}
+
+	private:
+		const halfedge::Mesh *m_mesh;
+		uint32_t m_current;
+	};
+	ConstEdgeIterator edges() const {
+		return ConstEdgeIterator(this);
+	}
+
+	// @@ Add half-edge iterator.
+
+	bool isValid() const {
+		// Make sure all edges are valid.
+		const uint32_t edgeCount = m_edgeArray.size();
+		for (uint32_t e = 0; e < edgeCount; e++) {
+			Edge *edge = m_edgeArray[e];
+			if (edge != NULL) {
+				if (edge->id != 2 * e) {
+					return false;
+				}
+				if (!edge->isValid()) {
+					return false;
+				}
+				if (edge->pair->id != 2 * e + 1) {
+					return false;
+				}
+				if (!edge->pair->isValid()) {
+					return false;
+				}
+			}
+		}
+		// @@ Make sure all faces are valid.
+		// @@ Make sure all vertices are valid.
+		return true;
+	}
+
+	// Error status:
+
+	struct ErrorCode {
+		enum Enum {
+			AlreadyAddedEdge,
+			DegenerateColocalEdge,
+			DegenerateEdge,
+			DuplicateEdge
+		};
+	};
+
+	mutable ErrorCode::Enum errorCode;
+	mutable uint32_t errorIndex0;
+	mutable uint32_t errorIndex1;
+
+private:
+	// Return true if the face can be added to the manifold mesh.
+	bool canAddFace(const std::vector<uint32_t> &indexArray, uint32_t first, uint32_t num) const {
+		return canAddFace(indexArray.data(), first, num);
+	}
+
+	bool canAddFace(const uint32_t *indexArray, uint32_t first, uint32_t num) const {
+		for (uint32_t j = num - 1, i = 0; i < num; j = i++) {
+			if (!canAddEdge(indexArray[first + j], indexArray[first + i])) {
+				errorIndex0 = indexArray[first + j];
+				errorIndex1 = indexArray[first + i];
+				return false;
+			}
+		}
+		// We also have to make sure the face does not have any duplicate edge!
+		for (uint32_t i = 0; i < num; i++) {
+			int i0 = indexArray[first + i + 0];
+			int i1 = indexArray[first + (i + 1) % num];
+			for (uint32_t j = i + 1; j < num; j++) {
+				int j0 = indexArray[first + j + 0];
+				int j1 = indexArray[first + (j + 1) % num];
+				if (i0 == j0 && i1 == j1) {
+					errorCode = ErrorCode::DuplicateEdge;
+					errorIndex0 = i0;
+					errorIndex1 = i1;
+					return false;
+				}
+			}
+		}
+		return true;
+	}
+
+	// Return true if the edge doesn't exist or doesn't have any adjacent face.
+	bool canAddEdge(uint32_t i, uint32_t j) const {
+		if (i == j) {
+			// Skip degenerate edges.
+			errorCode = ErrorCode::DegenerateEdge;
+			return false;
+		}
+		// Same check, but taking into account colocal vertices.
+		const Vertex *v0 = vertexAt(i);
+		const Vertex *v1 = vertexAt(j);
+		for (Vertex::ConstVertexIterator it(v0->colocals()); !it.isDone(); it.advance()) {
+			if (it.current() == v1) {
+				// Skip degenerate edges.
+				errorCode = ErrorCode::DegenerateColocalEdge;
+				return false;
+			}
+		}
+		// Make sure edge has not been added yet.
+		Edge *edge = findEdge(i, j);
+		// We ignore edges that don't have an adjacent face yet, since this face could become the edge's face.
+		if (!(edge == NULL || edge->face == NULL)) {
+			errorCode = ErrorCode::AlreadyAddedEdge;
+			return false;
+		}
+		return true;
+	}
+
+	Edge *addEdge(uint32_t i, uint32_t j) {
+		xaAssert(i != j);
+		Edge *edge = findEdge(i, j);
+		if (edge != NULL) {
+			// Edge may already exist, but its face must not be set.
+			xaDebugAssert(edge->face == NULL);
+			// Nothing else to do!
+		} else {
+			// Add new edge.
+			// Lookup pair.
+			Edge *pair = findEdge(j, i);
+			if (pair != NULL) {
+				// Create edge with same id.
+				edge = new Edge(pair->id + 1);
+				// Link edge pairs.
+				edge->pair = pair;
+				pair->pair = edge;
+				// @@ I'm not sure this is necessary!
+				pair->vertex->setEdge(pair);
+			} else {
+				// Create edge.
+				edge = new Edge(2 * m_edgeArray.size());
+				// Add only unpaired edges.
+				m_edgeArray.push_back(edge);
+			}
+			edge->vertex = m_vertexArray[i];
+			m_edgeMap[Key(i, j)] = edge;
+		}
+		// Face and Next are set by addFace.
+		return edge;
+	}
+
+	/// Find edge, test all colocals.
+	Edge *findEdge(uint32_t i, uint32_t j) const {
+		Edge *edge = NULL;
+		const Vertex *v0 = vertexAt(i);
+		const Vertex *v1 = vertexAt(j);
+		// Test all colocal pairs.
+		for (Vertex::ConstVertexIterator it0(v0->colocals()); !it0.isDone(); it0.advance()) {
+			for (Vertex::ConstVertexIterator it1(v1->colocals()); !it1.isDone(); it1.advance()) {
+				Key key(it0.current()->id, it1.current()->id);
+				if (edge == NULL) {
+					auto edgeIt = m_edgeMap.find(key);
+					if (edgeIt != m_edgeMap.end())
+						edge = (*edgeIt).second;
+#if !defined(_DEBUG)
+					if (edge != NULL) return edge;
+#endif
+				} else {
+					// Make sure that only one edge is found.
+					xaDebugAssert(m_edgeMap.find(key) == m_edgeMap.end());
+				}
+			}
+		}
+		return edge;
+	}
+
+	/// Link this boundary edge.
+	void linkBoundaryEdge(Edge *edge) {
+		xaAssert(edge->face == NULL);
+		// Make sure next pointer has not been set. @@ We want to be able to relink boundary edges after mesh changes.
+		Edge *next = edge;
+		while (next->pair->face != NULL) {
+			// Get pair prev
+			Edge *e = next->pair->next;
+			while (e->next != next->pair) {
+				e = e->next;
+			}
+			next = e;
+		}
+		edge->setNext(next->pair);
+		// Adjust vertex edge, so that it's the boundary edge. (required for isBoundary())
+		if (edge->vertex->edge != edge) {
+			// Multiple boundaries in the same edge.
+			edge->vertex->edge = edge;
+		}
+	}
+
+	Vertex *splitBoundaryEdge(Edge *edge, float t, const Vector3 &pos) {
+		/*
+		  We want to go from this configuration:
+
+				+   +
+				|   ^
+		   edge |<->|  pair
+				v   |
+				+   +
+
+		  To this one:
+
+				+   +
+				|   ^
+			 e0 |<->| p0
+				v   |
+		 vertex +   +
+				|   ^
+			 e1 |<->| p1
+				v   |
+				+   +
+
+		*/
+		Edge *pair = edge->pair;
+		// Make sure boundaries are linked.
+		xaDebugAssert(pair != NULL);
+		// Make sure edge is a boundary edge.
+		xaDebugAssert(pair->face == NULL);
+		// Add new vertex.
+		Vertex *vertex = addVertex(pos);
+		vertex->nor = lerp(edge->from()->nor, edge->to()->nor, t);
+		vertex->tex = lerp(edge->from()->tex, edge->to()->tex, t);
+		disconnect(edge);
+		disconnect(pair);
+		// Add edges.
+		Edge *e0 = addEdge(edge->from()->id, vertex->id);
+		Edge *p0 = addEdge(vertex->id, pair->to()->id);
+		Edge *e1 = addEdge(vertex->id, edge->to()->id);
+		Edge *p1 = addEdge(pair->from()->id, vertex->id);
+		// Link edges.
+		e0->setNext(e1);
+		p1->setNext(p0);
+		e0->setPrev(edge->prev);
+		e1->setNext(edge->next);
+		p1->setPrev(pair->prev);
+		p0->setNext(pair->next);
+		xaDebugAssert(e0->next == e1);
+		xaDebugAssert(e1->prev == e0);
+		xaDebugAssert(p1->next == p0);
+		xaDebugAssert(p0->prev == p1);
+		xaDebugAssert(p0->pair == e0);
+		xaDebugAssert(e0->pair == p0);
+		xaDebugAssert(p1->pair == e1);
+		xaDebugAssert(e1->pair == p1);
+		// Link faces.
+		e0->face = edge->face;
+		e1->face = edge->face;
+		// Link vertices.
+		edge->from()->setEdge(e0);
+		vertex->setEdge(e1);
+		delete edge;
+		delete pair;
+		return vertex;
+	}
+
+private:
+	std::vector<Vertex *> m_vertexArray;
+	std::vector<Edge *> m_edgeArray;
+	std::vector<Face *> m_faceArray;
+
+	struct Key {
+		Key() {}
+		Key(const Key &k) :
+				p0(k.p0),
+				p1(k.p1) {}
+		Key(uint32_t v0, uint32_t v1) :
+				p0(v0),
+				p1(v1) {}
+		void operator=(const Key &k) {
+			p0 = k.p0;
+			p1 = k.p1;
+		}
+		bool operator==(const Key &k) const {
+			return p0 == k.p0 && p1 == k.p1;
+		}
+
+		uint32_t p0;
+		uint32_t p1;
+	};
+
+	friend struct Hash<Mesh::Key>;
+	std::unordered_map<Key, Edge *, Hash<Key>, Equal<Key> > m_edgeMap;
+	uint32_t m_colocalVertexCount;
+};
+
+class MeshTopology {
+public:
+	MeshTopology(const Mesh *mesh) {
+		buildTopologyInfo(mesh);
+	}
+
+	/// Determine if the mesh is connected.
+	bool isConnected() const {
+		return m_connectedCount == 1;
+	}
+
+	/// Determine if the mesh is closed. (Each edge is shared by two faces)
+	bool isClosed() const {
+		return m_boundaryCount == 0;
+	}
+
+	/// Return true if the mesh has the topology of a disk.
+	bool isDisk() const {
+		return isConnected() && m_boundaryCount == 1 /* && m_eulerNumber == 1*/;
+	}
+
+private:
+	void buildTopologyInfo(const Mesh *mesh) {
+		const uint32_t vertexCount = mesh->colocalVertexCount();
+		const uint32_t faceCount = mesh->faceCount();
+		const uint32_t edgeCount = mesh->edgeCount();
+		xaPrint("--- Building mesh topology:\n");
+		std::vector<uint32_t> stack(faceCount);
+		BitArray bitFlags(faceCount);
+		bitFlags.clearAll();
+		// Compute connectivity.
+		xaPrint("---   Computing connectivity.\n");
+		m_connectedCount = 0;
+		for (uint32_t f = 0; f < faceCount; f++) {
+			if (bitFlags.bitAt(f) == false) {
+				m_connectedCount++;
+				stack.push_back(f);
+				while (!stack.empty()) {
+					const uint32_t top = stack.back();
+					xaAssert(top != uint32_t(~0));
+					stack.pop_back();
+					if (bitFlags.bitAt(top) == false) {
+						bitFlags.setBitAt(top);
+						const Face *face = mesh->faceAt(top);
+						const Edge *firstEdge = face->edge;
+						const Edge *edge = firstEdge;
+						do {
+							const Face *neighborFace = edge->pair->face;
+							if (neighborFace != NULL) {
+								stack.push_back(neighborFace->id);
+							}
+							edge = edge->next;
+						} while (edge != firstEdge);
+					}
+				}
+			}
+		}
+		xaAssert(stack.empty());
+		xaPrint("---   %d connected components.\n", m_connectedCount);
+		// Count boundary loops.
+		xaPrint("---   Counting boundary loops.\n");
+		m_boundaryCount = 0;
+		bitFlags.resize(edgeCount);
+		bitFlags.clearAll();
+		// Don't forget to link the boundary otherwise this won't work.
+		for (uint32_t e = 0; e < edgeCount; e++) {
+			const Edge *startEdge = mesh->edgeAt(e);
+			if (startEdge != NULL && startEdge->isBoundary() && bitFlags.bitAt(e) == false) {
+				xaDebugAssert(startEdge->face != NULL);
+				xaDebugAssert(startEdge->pair->face == NULL);
+				startEdge = startEdge->pair;
+				m_boundaryCount++;
+				const Edge *edge = startEdge;
+				do {
+					bitFlags.setBitAt(edge->id / 2);
+					edge = edge->next;
+				} while (startEdge != edge);
+			}
+		}
+		xaPrint("---   %d boundary loops found.\n", m_boundaryCount);
+		// Compute euler number.
+		m_eulerNumber = vertexCount - edgeCount + faceCount;
+		xaPrint("---   Euler number: %d.\n", m_eulerNumber);
+		// Compute genus. (only valid on closed connected surfaces)
+		m_genus = -1;
+		if (isClosed() && isConnected()) {
+			m_genus = (2 - m_eulerNumber) / 2;
+			xaPrint("---   Genus: %d.\n", m_genus);
+		}
+	}
+
+private:
+	///< Number of boundary loops.
+	int m_boundaryCount;
+
+	///< Number of connected components.
+	int m_connectedCount;
+
+	///< Euler number.
+	int m_eulerNumber;
+
+	/// Mesh genus.
+	int m_genus;
+};
+
+float computeSurfaceArea(const halfedge::Mesh *mesh) {
+	float area = 0;
+	for (halfedge::Mesh::ConstFaceIterator it(mesh->faces()); !it.isDone(); it.advance()) {
+		const halfedge::Face *face = it.current();
+		area += face->area();
+	}
+	xaDebugAssert(area >= 0);
+	return area;
+}
+
+float computeParametricArea(const halfedge::Mesh *mesh) {
+	float area = 0;
+	for (halfedge::Mesh::ConstFaceIterator it(mesh->faces()); !it.isDone(); it.advance()) {
+		const halfedge::Face *face = it.current();
+		area += face->parametricArea();
+	}
+	return area;
+}
+
+uint32_t countMeshTriangles(const Mesh *mesh) {
+	const uint32_t faceCount = mesh->faceCount();
+	uint32_t triangleCount = 0;
+	for (uint32_t f = 0; f < faceCount; f++) {
+		const Face *face = mesh->faceAt(f);
+		uint32_t edgeCount = face->edgeCount();
+		xaDebugAssert(edgeCount > 2);
+		triangleCount += edgeCount - 2;
+	}
+	return triangleCount;
+}
+
+Mesh *unifyVertices(const Mesh *inputMesh) {
+	Mesh *mesh = new Mesh;
+	// Only add the first colocal.
+	const uint32_t vertexCount = inputMesh->vertexCount();
+	for (uint32_t v = 0; v < vertexCount; v++) {
+		const Vertex *vertex = inputMesh->vertexAt(v);
+		if (vertex->isFirstColocal()) {
+			mesh->addVertex(vertex->pos);
+		}
+	}
+	std::vector<uint32_t> indexArray;
+	// Add new faces pointing to first colocals.
+	uint32_t faceCount = inputMesh->faceCount();
+	for (uint32_t f = 0; f < faceCount; f++) {
+		const Face *face = inputMesh->faceAt(f);
+		indexArray.clear();
+		for (Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) {
+			const Edge *edge = it.current();
+			const Vertex *vertex = edge->vertex->firstColocal();
+			indexArray.push_back(vertex->id);
+		}
+		mesh->addFace(indexArray);
+	}
+	mesh->linkBoundary();
+	return mesh;
+}
+
+static bool pointInTriangle(const Vector2 &p, const Vector2 &a, const Vector2 &b, const Vector2 &c) {
+	return triangleArea(a, b, p) >= 0.00001f &&
+		   triangleArea(b, c, p) >= 0.00001f &&
+		   triangleArea(c, a, p) >= 0.00001f;
+}
+
+// This is doing a simple ear-clipping algorithm that skips invalid triangles. Ideally, we should
+// also sort the ears by angle, start with the ones that have the smallest angle and proceed in order.
+Mesh *triangulate(const Mesh *inputMesh) {
+	Mesh *mesh = new Mesh;
+	// Add all vertices.
+	const uint32_t vertexCount = inputMesh->vertexCount();
+	for (uint32_t v = 0; v < vertexCount; v++) {
+		const Vertex *vertex = inputMesh->vertexAt(v);
+		mesh->addVertex(vertex->pos);
+	}
+	std::vector<int> polygonVertices;
+	std::vector<float> polygonAngles;
+	std::vector<Vector2> polygonPoints;
+	const uint32_t faceCount = inputMesh->faceCount();
+	for (uint32_t f = 0; f < faceCount; f++) {
+		const Face *face = inputMesh->faceAt(f);
+		xaDebugAssert(face != NULL);
+		const uint32_t edgeCount = face->edgeCount();
+		xaDebugAssert(edgeCount >= 3);
+		polygonVertices.clear();
+		polygonVertices.reserve(edgeCount);
+		if (edgeCount == 3) {
+			// Simple case for triangles.
+			for (Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) {
+				const Edge *edge = it.current();
+				const Vertex *vertex = edge->vertex;
+				polygonVertices.push_back(vertex->id);
+			}
+			int v0 = polygonVertices[0];
+			int v1 = polygonVertices[1];
+			int v2 = polygonVertices[2];
+			mesh->addFace(v0, v1, v2);
+		} else {
+			// Build 2D polygon projecting vertices onto normal plane.
+			// Faces are not necesarily planar, this is for example the case, when the face comes from filling a hole. In such cases
+			// it's much better to use the best fit plane.
+			const Vector3 fn = face->normal();
+			Basis basis;
+			basis.buildFrameForDirection(fn);
+			polygonPoints.clear();
+			polygonPoints.reserve(edgeCount);
+			polygonAngles.clear();
+			polygonAngles.reserve(edgeCount);
+			for (Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) {
+				const Edge *edge = it.current();
+				const Vertex *vertex = edge->vertex;
+				polygonVertices.push_back(vertex->id);
+				Vector2 p;
+				p.x = dot(basis.tangent, vertex->pos);
+				p.y = dot(basis.bitangent, vertex->pos);
+				polygonPoints.push_back(p);
+			}
+			polygonAngles.resize(edgeCount);
+			while (polygonVertices.size() > 2) {
+				uint32_t size = polygonVertices.size();
+				// Update polygon angles. @@ Update only those that have changed.
+				float minAngle = 2 * PI;
+				uint32_t bestEar = 0; // Use first one if none of them is valid.
+				bool bestIsValid = false;
+				for (uint32_t i = 0; i < size; i++) {
+					uint32_t i0 = i;
+					uint32_t i1 = (i + 1) % size; // Use Sean's polygon interation trick.
+					uint32_t i2 = (i + 2) % size;
+					Vector2 p0 = polygonPoints[i0];
+					Vector2 p1 = polygonPoints[i1];
+					Vector2 p2 = polygonPoints[i2];
+
+					bool degenerate = distance(p0, p1) < NV_EPSILON || distance(p0, p2) < NV_EPSILON || distance(p1, p2) < NV_EPSILON;
+					if (degenerate) {
+						continue;
+					}
+
+					float d = clamp(dot(p0 - p1, p2 - p1) / (length(p0 - p1) * length(p2 - p1)), -1.0f, 1.0f);
+					float angle = acosf(d);
+					float area = triangleArea(p0, p1, p2);
+					if (area < 0.0f) angle = 2.0f * PI - angle;
+					polygonAngles[i1] = angle;
+					if (angle < minAngle || !bestIsValid) {
+						// Make sure this is a valid ear, if not, skip this point.
+						bool valid = true;
+						for (uint32_t j = 0; j < size; j++) {
+							if (j == i0 || j == i1 || j == i2) continue;
+							Vector2 p = polygonPoints[j];
+							if (pointInTriangle(p, p0, p1, p2)) {
+								valid = false;
+								break;
+							}
+						}
+						if (valid || !bestIsValid) {
+							minAngle = angle;
+							bestEar = i1;
+							bestIsValid = valid;
+						}
+					}
+				}
+				if (!bestIsValid)
+					break;
+
+				xaDebugAssert(minAngle <= 2 * PI);
+				// Clip best ear:
+				uint32_t i0 = (bestEar + size - 1) % size;
+				uint32_t i1 = (bestEar + 0) % size;
+				uint32_t i2 = (bestEar + 1) % size;
+				int v0 = polygonVertices[i0];
+				int v1 = polygonVertices[i1];
+				int v2 = polygonVertices[i2];
+				mesh->addFace(v0, v1, v2);
+				polygonVertices.erase(polygonVertices.begin() + i1);
+				polygonPoints.erase(polygonPoints.begin() + i1);
+				polygonAngles.erase(polygonAngles.begin() + i1);
+			}
+		}
+	}
+	mesh->linkBoundary();
+	return mesh;
+}
+
+} //  namespace halfedge
+
+/// Mersenne twister random number generator.
+class MTRand {
+public:
+	enum time_e { Time };
+	enum { N = 624 }; // length of state vector
+	enum { M = 397 };
+
+	/// Constructor that uses the current time as the seed.
+	MTRand(time_e) {
+		seed((uint32_t)time(NULL));
+	}
+
+	/// Constructor that uses the given seed.
+	MTRand(uint32_t s = 0) {
+		seed(s);
+	}
+
+	/// Provide a new seed.
+	void seed(uint32_t s) {
+		initialize(s);
+		reload();
+	}
+
+	/// Get a random number between 0 - 65536.
+	uint32_t get() {
+		// Pull a 32-bit integer from the generator state
+		// Every other access function simply transforms the numbers extracted here
+		if (left == 0) {
+			reload();
+		}
+		left--;
+		uint32_t s1;
+		s1 = *next++;
+		s1 ^= (s1 >> 11);
+		s1 ^= (s1 << 7) & 0x9d2c5680U;
+		s1 ^= (s1 << 15) & 0xefc60000U;
+		return (s1 ^ (s1 >> 18));
+	};
+
+	/// Get a random number on [0, max] interval.
+	uint32_t getRange(uint32_t max) {
+		if (max == 0) return 0;
+		if (max == NV_UINT32_MAX) return get();
+		const uint32_t np2 = nextPowerOfTwo(max + 1); // @@ This fails if max == NV_UINT32_MAX
+		const uint32_t mask = np2 - 1;
+		uint32_t n;
+		do {
+			n = get() & mask;
+		} while (n > max);
+		return n;
+	}
+
+private:
+	void initialize(uint32_t seed) {
+		// Initialize generator state with seed
+		// See Knuth TAOCP Vol 2, 3rd Ed, p.106 for multiplier.
+		// In previous versions, most significant bits (MSBs) of the seed affect
+		// only MSBs of the state array.  Modified 9 Jan 2002 by Makoto Matsumoto.
+		uint32_t *s = state;
+		uint32_t *r = state;
+		int i = 1;
+		*s++ = seed & 0xffffffffUL;
+		for (; i < N; ++i) {
+			*s++ = (1812433253UL * (*r ^ (*r >> 30)) + i) & 0xffffffffUL;
+			r++;
+		}
+	}
+
+	void reload() {
+		// Generate N new values in state
+		// Made clearer and faster by Matthew Bellew (matthew.bellew@home.com)
+		uint32_t *p = state;
+		int i;
+		for (i = N - M; i--; ++p)
+			*p = twist(p[M], p[0], p[1]);
+		for (i = M; --i; ++p)
+			*p = twist(p[M - N], p[0], p[1]);
+		*p = twist(p[M - N], p[0], state[0]);
+		left = N, next = state;
+	}
+
+	uint32_t hiBit(uint32_t u) const {
+		return u & 0x80000000U;
+	}
+	uint32_t loBit(uint32_t u) const {
+		return u & 0x00000001U;
+	}
+	uint32_t loBits(uint32_t u) const {
+		return u & 0x7fffffffU;
+	}
+	uint32_t mixBits(uint32_t u, uint32_t v) const {
+		return hiBit(u) | loBits(v);
+	}
+	uint32_t twist(uint32_t m, uint32_t s0, uint32_t s1) const {
+		return m ^ (mixBits(s0, s1) >> 1) ^ ((~loBit(s1) + 1) & 0x9908b0dfU);
+	}
+
+	uint32_t state[N]; // internal state
+	uint32_t *next; // next value to get from state
+	int left; // number of values left before reload needed
+};
+
+namespace morton {
+// Code from ryg:
+// http://fgiesen.wordpress.com/2009/12/13/decoding-morton-codes/
+
+// Inverse of part1By1 - "delete" all odd-indexed bits
+uint32_t compact1By1(uint32_t x) {
+	x &= 0x55555555; // x = -f-e -d-c -b-a -9-8 -7-6 -5-4 -3-2 -1-0
+	x = (x ^ (x >> 1)) & 0x33333333; // x = --fe --dc --ba --98 --76 --54 --32 --10
+	x = (x ^ (x >> 2)) & 0x0f0f0f0f; // x = ---- fedc ---- ba98 ---- 7654 ---- 3210
+	x = (x ^ (x >> 4)) & 0x00ff00ff; // x = ---- ---- fedc ba98 ---- ---- 7654 3210
+	x = (x ^ (x >> 8)) & 0x0000ffff; // x = ---- ---- ---- ---- fedc ba98 7654 3210
+	return x;
+}
+
+// Inverse of part1By2 - "delete" all bits not at positions divisible by 3
+uint32_t compact1By2(uint32_t x) {
+	x &= 0x09249249; // x = ---- 9--8 --7- -6-- 5--4 --3- -2-- 1--0
+	x = (x ^ (x >> 2)) & 0x030c30c3; // x = ---- --98 ---- 76-- --54 ---- 32-- --10
+	x = (x ^ (x >> 4)) & 0x0300f00f; // x = ---- --98 ---- ---- 7654 ---- ---- 3210
+	x = (x ^ (x >> 8)) & 0xff0000ff; // x = ---- --98 ---- ---- ---- ---- 7654 3210
+	x = (x ^ (x >> 16)) & 0x000003ff; // x = ---- ---- ---- ---- ---- --98 7654 3210
+	return x;
+}
+
+uint32_t decodeMorton2X(uint32_t code) {
+	return compact1By1(code >> 0);
+}
+
+uint32_t decodeMorton2Y(uint32_t code) {
+	return compact1By1(code >> 1);
+}
+
+uint32_t decodeMorton3X(uint32_t code) {
+	return compact1By2(code >> 0);
+}
+
+uint32_t decodeMorton3Y(uint32_t code) {
+	return compact1By2(code >> 1);
+}
+
+uint32_t decodeMorton3Z(uint32_t code) {
+	return compact1By2(code >> 2);
+}
+} // namespace morton
+
+// A simple, dynamic proximity grid based on Jon's code.
+// Instead of storing pointers here I store indices.
+struct ProximityGrid {
+	void init(const Box &box, uint32_t count) {
+		cellArray.clear();
+		// Determine grid size.
+		float cellWidth;
+		Vector3 diagonal = box.extents() * 2.f;
+		float volume = box.volume();
+		if (equal(volume, 0)) {
+			// Degenerate box, treat like a quad.
+			Vector2 quad;
+			if (diagonal.x < diagonal.y && diagonal.x < diagonal.z) {
+				quad.x = diagonal.y;
+				quad.y = diagonal.z;
+			} else if (diagonal.y < diagonal.x && diagonal.y < diagonal.z) {
+				quad.x = diagonal.x;
+				quad.y = diagonal.z;
+			} else {
+				quad.x = diagonal.x;
+				quad.y = diagonal.y;
+			}
+			float cellArea = quad.x * quad.y / count;
+			cellWidth = sqrtf(cellArea); // pow(cellArea, 1.0f / 2.0f);
+		} else {
+			// Ideally we want one cell per point.
+			float cellVolume = volume / count;
+			cellWidth = powf(cellVolume, 1.0f / 3.0f);
+		}
+		xaDebugAssert(cellWidth != 0);
+		sx = std::max(1, ftoi_ceil(diagonal.x / cellWidth));
+		sy = std::max(1, ftoi_ceil(diagonal.y / cellWidth));
+		sz = std::max(1, ftoi_ceil(diagonal.z / cellWidth));
+		invCellSize.x = float(sx) / diagonal.x;
+		invCellSize.y = float(sy) / diagonal.y;
+		invCellSize.z = float(sz) / diagonal.z;
+		cellArray.resize(sx * sy * sz);
+		corner = box.minCorner; // @@ Align grid better?
+	}
+
+	int index_x(float x) const {
+		return clamp(ftoi_floor((x - corner.x) * invCellSize.x), 0, sx - 1);
+	}
+
+	int index_y(float y) const {
+		return clamp(ftoi_floor((y - corner.y) * invCellSize.y), 0, sy - 1);
+	}
+
+	int index_z(float z) const {
+		return clamp(ftoi_floor((z - corner.z) * invCellSize.z), 0, sz - 1);
+	}
+
+	int index(int x, int y, int z) const {
+		xaDebugAssert(x >= 0 && x < sx);
+		xaDebugAssert(y >= 0 && y < sy);
+		xaDebugAssert(z >= 0 && z < sz);
+		int idx = (z * sy + y) * sx + x;
+		xaDebugAssert(idx >= 0 && uint32_t(idx) < cellArray.size());
+		return idx;
+	}
+
+	uint32_t mortonCount() const {
+		uint64_t s = uint64_t(max3(sx, sy, sz));
+		s = nextPowerOfTwo(s);
+		if (s > 1024) {
+			return uint32_t(s * s * min3(sx, sy, sz));
+		}
+		return uint32_t(s * s * s);
+	}
+
+	int mortonIndex(uint32_t code) const {
+		uint32_t x, y, z;
+		uint32_t s = uint32_t(max3(sx, sy, sz));
+		if (s > 1024) {
+			// Use layered two-dimensional morton order.
+			s = nextPowerOfTwo(s);
+			uint32_t layer = code / (s * s);
+			code = code % (s * s);
+			uint32_t layer_count = uint32_t(min3(sx, sy, sz));
+			if (sx == (int)layer_count) {
+				x = layer;
+				y = morton::decodeMorton2X(code);
+				z = morton::decodeMorton2Y(code);
+			} else if (sy == (int)layer_count) {
+				x = morton::decodeMorton2Y(code);
+				y = layer;
+				z = morton::decodeMorton2X(code);
+			} else { /*if (sz == layer_count)*/
+				x = morton::decodeMorton2X(code);
+				y = morton::decodeMorton2Y(code);
+				z = layer;
+			}
+		} else {
+			x = morton::decodeMorton3X(code);
+			y = morton::decodeMorton3Y(code);
+			z = morton::decodeMorton3Z(code);
+		}
+		if (x >= uint32_t(sx) || y >= uint32_t(sy) || z >= uint32_t(sz)) {
+			return -1;
+		}
+		return index(x, y, z);
+	}
+
+	void add(const Vector3 &pos, uint32_t key) {
+		int x = index_x(pos.x);
+		int y = index_y(pos.y);
+		int z = index_z(pos.z);
+		uint32_t idx = index(x, y, z);
+		cellArray[idx].indexArray.push_back(key);
+	}
+
+	// Gather all points inside the given sphere.
+	// Radius is assumed to be small, so we don't bother culling the cells.
+	void gather(const Vector3 &position, float radius, std::vector<uint32_t> &indexArray) {
+		int x0 = index_x(position.x - radius);
+		int x1 = index_x(position.x + radius);
+		int y0 = index_y(position.y - radius);
+		int y1 = index_y(position.y + radius);
+		int z0 = index_z(position.z - radius);
+		int z1 = index_z(position.z + radius);
+		for (int z = z0; z <= z1; z++) {
+			for (int y = y0; y <= y1; y++) {
+				for (int x = x0; x <= x1; x++) {
+					int idx = index(x, y, z);
+					indexArray.insert(indexArray.begin(), cellArray[idx].indexArray.begin(), cellArray[idx].indexArray.end());
+				}
+			}
+		}
+	}
+
+	struct Cell {
+		std::vector<uint32_t> indexArray;
+	};
+
+	std::vector<Cell> cellArray;
+
+	Vector3 corner;
+	Vector3 invCellSize;
+	int sx, sy, sz;
+};
+
+// Based on Pierre Terdiman's and Michael Herf's source code.
+// http://www.codercorner.com/RadixSortRevisited.htm
+// http://www.stereopsis.com/radix.html
+class RadixSort {
+public:
+	RadixSort() :
+			m_size(0),
+			m_ranks(NULL),
+			m_ranks2(NULL),
+			m_validRanks(false) {}
+	~RadixSort() {
+		// Release everything
+		free(m_ranks2);
+		free(m_ranks);
+	}
+
+	RadixSort &sort(const float *input, uint32_t count) {
+		if (input == NULL || count == 0) return *this;
+		// Resize lists if needed
+		if (count != m_size) {
+			if (count > m_size) {
+				m_ranks2 = (uint32_t *)realloc(m_ranks2, sizeof(uint32_t) * count);
+				m_ranks = (uint32_t *)realloc(m_ranks, sizeof(uint32_t) * count);
+			}
+			m_size = count;
+			m_validRanks = false;
+		}
+		if (count < 32) {
+			insertionSort(input, count);
+		} else {
+			// @@ Avoid touching the input multiple times.
+			for (uint32_t i = 0; i < count; i++) {
+				FloatFlip((uint32_t &)input[i]);
+			}
+			radixSort<uint32_t>((const uint32_t *)input, count);
+			for (uint32_t i = 0; i < count; i++) {
+				IFloatFlip((uint32_t &)input[i]);
+			}
+		}
+		return *this;
+	}
+
+	RadixSort &sort(const std::vector<float> &input) {
+		return sort(input.data(), input.size());
+	}
+
+	// Access to results. m_ranks is a list of indices in sorted order, i.e. in the order you may further process your data
+	const uint32_t *ranks() const {
+		xaDebugAssert(m_validRanks);
+		return m_ranks;
+	}
+	uint32_t *ranks() {
+		xaDebugAssert(m_validRanks);
+		return m_ranks;
+	}
+
+private:
+	uint32_t m_size;
+	uint32_t *m_ranks;
+	uint32_t *m_ranks2;
+	bool m_validRanks;
+
+	void FloatFlip(uint32_t &f) {
+		int32_t mask = (int32_t(f) >> 31) | 0x80000000; // Warren Hunt, Manchor Ko.
+		f ^= mask;
+	}
+
+	void IFloatFlip(uint32_t &f) {
+		uint32_t mask = ((f >> 31) - 1) | 0x80000000; // Michael Herf.
+		f ^= mask;
+	}
+
+	template <typename T>
+	void createHistograms(const T *buffer, uint32_t count, uint32_t *histogram) {
+		const uint32_t bucketCount = sizeof(T); // (8 * sizeof(T)) / log2(radix)
+		// Init bucket pointers.
+		uint32_t *h[bucketCount];
+		for (uint32_t i = 0; i < bucketCount; i++) {
+			h[i] = histogram + 256 * i;
+		}
+		// Clear histograms.
+		memset(histogram, 0, 256 * bucketCount * sizeof(uint32_t));
+		// @@ Add support for signed integers.
+		// Build histograms.
+		const uint8_t *p = (const uint8_t *)buffer; // @@ Does this break aliasing rules?
+		const uint8_t *pe = p + count * sizeof(T);
+		while (p != pe) {
+			h[0][*p++]++, h[1][*p++]++, h[2][*p++]++, h[3][*p++]++;
+#ifdef _MSC_VER
+#pragma warning(push)
+#pragma warning(disable : 4127)
+#endif
+			if (bucketCount == 8) h[4][*p++]++, h[5][*p++]++, h[6][*p++]++, h[7][*p++]++;
+#ifdef _MSC_VER
+#pragma warning(pop)
+#endif
+		}
+	}
+
+	template <typename T>
+	void insertionSort(const T *input, uint32_t count) {
+		if (!m_validRanks) {
+			m_ranks[0] = 0;
+			for (uint32_t i = 1; i != count; ++i) {
+				int rank = m_ranks[i] = i;
+				uint32_t j = i;
+				while (j != 0 && input[rank] < input[m_ranks[j - 1]]) {
+					m_ranks[j] = m_ranks[j - 1];
+					--j;
+				}
+				if (i != j) {
+					m_ranks[j] = rank;
+				}
+			}
+			m_validRanks = true;
+		} else {
+			for (uint32_t i = 1; i != count; ++i) {
+				int rank = m_ranks[i];
+				uint32_t j = i;
+				while (j != 0 && input[rank] < input[m_ranks[j - 1]]) {
+					m_ranks[j] = m_ranks[j - 1];
+					--j;
+				}
+				if (i != j) {
+					m_ranks[j] = rank;
+				}
+			}
+		}
+	}
+
+	template <typename T>
+	void radixSort(const T *input, uint32_t count) {
+		const uint32_t P = sizeof(T); // pass count
+		// Allocate histograms & offsets on the stack
+		uint32_t histogram[256 * P];
+		uint32_t *link[256];
+		createHistograms(input, count, histogram);
+		// Radix sort, j is the pass number (0=LSB, P=MSB)
+		for (uint32_t j = 0; j < P; j++) {
+			// Pointer to this bucket.
+			const uint32_t *h = &histogram[j * 256];
+			const uint8_t *inputBytes = (const uint8_t *)input; // @@ Is this aliasing legal?
+			inputBytes += j;
+			if (h[inputBytes[0]] == count) {
+				// Skip this pass, all values are the same.
+				continue;
+			}
+			// Create offsets
+			link[0] = m_ranks2;
+			for (uint32_t i = 1; i < 256; i++)
+				link[i] = link[i - 1] + h[i - 1];
+			// Perform Radix Sort
+			if (!m_validRanks) {
+				for (uint32_t i = 0; i < count; i++) {
+					*link[inputBytes[i * P]]++ = i;
+				}
+				m_validRanks = true;
+			} else {
+				for (uint32_t i = 0; i < count; i++) {
+					const uint32_t idx = m_ranks[i];
+					*link[inputBytes[idx * P]]++ = idx;
+				}
+			}
+			// Swap pointers for next pass. Valid indices - the most recent ones - are in m_ranks after the swap.
+			std::swap(m_ranks, m_ranks2);
+		}
+		// All values were equal, generate linear ranks.
+		if (!m_validRanks) {
+			for (uint32_t i = 0; i < count; i++) {
+				m_ranks[i] = i;
+			}
+			m_validRanks = true;
+		}
+	}
+};
+
+namespace raster {
+class ClippedTriangle {
+public:
+	ClippedTriangle(Vector2::Arg a, Vector2::Arg b, Vector2::Arg c) {
+		m_numVertices = 3;
+		m_activeVertexBuffer = 0;
+		m_verticesA[0] = a;
+		m_verticesA[1] = b;
+		m_verticesA[2] = c;
+		m_vertexBuffers[0] = m_verticesA;
+		m_vertexBuffers[1] = m_verticesB;
+	}
+
+	uint32_t vertexCount() {
+		return m_numVertices;
+	}
+
+	const Vector2 *vertices() {
+		return m_vertexBuffers[m_activeVertexBuffer];
+	}
+
+	void clipHorizontalPlane(float offset, float clipdirection) {
+		Vector2 *v = m_vertexBuffers[m_activeVertexBuffer];
+		m_activeVertexBuffer ^= 1;
+		Vector2 *v2 = m_vertexBuffers[m_activeVertexBuffer];
+		v[m_numVertices] = v[0];
+		float dy2, dy1 = offset - v[0].y;
+		int dy2in, dy1in = clipdirection * dy1 >= 0;
+		uint32_t p = 0;
+		for (uint32_t k = 0; k < m_numVertices; k++) {
+			dy2 = offset - v[k + 1].y;
+			dy2in = clipdirection * dy2 >= 0;
+			if (dy1in) v2[p++] = v[k];
+			if (dy1in + dy2in == 1) { // not both in/out
+				float dx = v[k + 1].x - v[k].x;
+				float dy = v[k + 1].y - v[k].y;
+				v2[p++] = Vector2(v[k].x + dy1 * (dx / dy), offset);
+			}
+			dy1 = dy2;
+			dy1in = dy2in;
+		}
+		m_numVertices = p;
+		//for (uint32_t k=0; k<m_numVertices; k++) printf("(%f, %f)\n", v2[k].x, v2[k].y); printf("\n");
+	}
+
+	void clipVerticalPlane(float offset, float clipdirection) {
+		Vector2 *v = m_vertexBuffers[m_activeVertexBuffer];
+		m_activeVertexBuffer ^= 1;
+		Vector2 *v2 = m_vertexBuffers[m_activeVertexBuffer];
+		v[m_numVertices] = v[0];
+		float dx2, dx1 = offset - v[0].x;
+		int dx2in, dx1in = clipdirection * dx1 >= 0;
+		uint32_t p = 0;
+		for (uint32_t k = 0; k < m_numVertices; k++) {
+			dx2 = offset - v[k + 1].x;
+			dx2in = clipdirection * dx2 >= 0;
+			if (dx1in) v2[p++] = v[k];
+			if (dx1in + dx2in == 1) { // not both in/out
+				float dx = v[k + 1].x - v[k].x;
+				float dy = v[k + 1].y - v[k].y;
+				v2[p++] = Vector2(offset, v[k].y + dx1 * (dy / dx));
+			}
+			dx1 = dx2;
+			dx1in = dx2in;
+		}
+		m_numVertices = p;
+	}
+
+	void computeAreaCentroid() {
+		Vector2 *v = m_vertexBuffers[m_activeVertexBuffer];
+		v[m_numVertices] = v[0];
+		m_area = 0;
+		float centroidx = 0, centroidy = 0;
+		for (uint32_t k = 0; k < m_numVertices; k++) {
+			// http://local.wasp.uwa.edu.au/~pbourke/geometry/polyarea/
+			float f = v[k].x * v[k + 1].y - v[k + 1].x * v[k].y;
+			m_area += f;
+			centroidx += f * (v[k].x + v[k + 1].x);
+			centroidy += f * (v[k].y + v[k + 1].y);
+		}
+		m_area = 0.5f * fabsf(m_area);
+		if (m_area == 0) {
+			m_centroid = Vector2(0.0f);
+		} else {
+			m_centroid = Vector2(centroidx / (6 * m_area), centroidy / (6 * m_area));
+		}
+	}
+
+	void clipAABox(float x0, float y0, float x1, float y1) {
+		clipVerticalPlane(x0, -1);
+		clipHorizontalPlane(y0, -1);
+		clipVerticalPlane(x1, 1);
+		clipHorizontalPlane(y1, 1);
+		computeAreaCentroid();
+	}
+
+	Vector2 centroid() {
+		return m_centroid;
+	}
+
+	float area() {
+		return m_area;
+	}
+
+private:
+	Vector2 m_verticesA[7 + 1];
+	Vector2 m_verticesB[7 + 1];
+	Vector2 *m_vertexBuffers[2];
+	uint32_t m_numVertices;
+	uint32_t m_activeVertexBuffer;
+	float m_area;
+	Vector2 m_centroid;
+};
+
+/// A callback to sample the environment. Return false to terminate rasterization.
+typedef bool (*SamplingCallback)(void *param, int x, int y, Vector3::Arg bar, Vector3::Arg dx, Vector3::Arg dy, float coverage);
+
+/// A triangle for rasterization.
+struct Triangle {
+	Triangle(Vector2::Arg v0, Vector2::Arg v1, Vector2::Arg v2, Vector3::Arg t0, Vector3::Arg t1, Vector3::Arg t2) {
+		// Init vertices.
+		this->v1 = v0;
+		this->v2 = v2;
+		this->v3 = v1;
+		// Set barycentric coordinates.
+		this->t1 = t0;
+		this->t2 = t2;
+		this->t3 = t1;
+		// make sure every triangle is front facing.
+		flipBackface();
+		// Compute deltas.
+		valid = computeDeltas();
+		computeUnitInwardNormals();
+	}
+
+	/// Compute texture space deltas.
+	/// This method takes two edge vectors that form a basis, determines the
+	/// coordinates of the canonic vectors in that basis, and computes the
+	/// texture gradient that corresponds to those vectors.
+	bool computeDeltas() {
+		Vector2 e0 = v3 - v1;
+		Vector2 e1 = v2 - v1;
+		Vector3 de0 = t3 - t1;
+		Vector3 de1 = t2 - t1;
+		float denom = 1.0f / (e0.y * e1.x - e1.y * e0.x);
+		if (!std::isfinite(denom)) {
+			return false;
+		}
+		float lambda1 = -e1.y * denom;
+		float lambda2 = e0.y * denom;
+		float lambda3 = e1.x * denom;
+		float lambda4 = -e0.x * denom;
+		dx = de0 * lambda1 + de1 * lambda2;
+		dy = de0 * lambda3 + de1 * lambda4;
+		return true;
+	}
+
+	bool draw(const Vector2 &extents, bool enableScissors, SamplingCallback cb, void *param) {
+		// 28.4 fixed-point coordinates
+		const int Y1 = ftoi_round(16.0f * v1.y);
+		const int Y2 = ftoi_round(16.0f * v2.y);
+		const int Y3 = ftoi_round(16.0f * v3.y);
+		const int X1 = ftoi_round(16.0f * v1.x);
+		const int X2 = ftoi_round(16.0f * v2.x);
+		const int X3 = ftoi_round(16.0f * v3.x);
+		// Deltas
+		const int DX12 = X1 - X2;
+		const int DX23 = X2 - X3;
+		const int DX31 = X3 - X1;
+		const int DY12 = Y1 - Y2;
+		const int DY23 = Y2 - Y3;
+		const int DY31 = Y3 - Y1;
+		// Fixed-point deltas
+		const int FDX12 = DX12 << 4;
+		const int FDX23 = DX23 << 4;
+		const int FDX31 = DX31 << 4;
+		const int FDY12 = DY12 << 4;
+		const int FDY23 = DY23 << 4;
+		const int FDY31 = DY31 << 4;
+		int minx, miny, maxx, maxy;
+		if (enableScissors) {
+			int frustumX0 = 0 << 4;
+			int frustumY0 = 0 << 4;
+			int frustumX1 = (int)extents.x << 4;
+			int frustumY1 = (int)extents.y << 4;
+			// Bounding rectangle
+			minx = (std::max(min3(X1, X2, X3), frustumX0) + 0xF) >> 4;
+			miny = (std::max(min3(Y1, Y2, Y3), frustumY0) + 0xF) >> 4;
+			maxx = (std::min(max3(X1, X2, X3), frustumX1) + 0xF) >> 4;
+			maxy = (std::min(max3(Y1, Y2, Y3), frustumY1) + 0xF) >> 4;
+		} else {
+			// Bounding rectangle
+			minx = (min3(X1, X2, X3) + 0xF) >> 4;
+			miny = (min3(Y1, Y2, Y3) + 0xF) >> 4;
+			maxx = (max3(X1, X2, X3) + 0xF) >> 4;
+			maxy = (max3(Y1, Y2, Y3) + 0xF) >> 4;
+		}
+		// Block size, standard 8x8 (must be power of two)
+		const int q = 8;
+		// @@ This won't work when minx,miny are negative. This code path is not used. Leaving as is for now.
+		xaAssert(minx >= 0);
+		xaAssert(miny >= 0);
+		// Start in corner of 8x8 block
+		minx &= ~(q - 1);
+		miny &= ~(q - 1);
+		// Half-edge constants
+		int C1 = DY12 * X1 - DX12 * Y1;
+		int C2 = DY23 * X2 - DX23 * Y2;
+		int C3 = DY31 * X3 - DX31 * Y3;
+		// Correct for fill convention
+		if (DY12 < 0 || (DY12 == 0 && DX12 > 0)) C1++;
+		if (DY23 < 0 || (DY23 == 0 && DX23 > 0)) C2++;
+		if (DY31 < 0 || (DY31 == 0 && DX31 > 0)) C3++;
+		// Loop through blocks
+		for (int y = miny; y < maxy; y += q) {
+			for (int x = minx; x < maxx; x += q) {
+				// Corners of block
+				int x0 = x << 4;
+				int x1 = (x + q - 1) << 4;
+				int y0 = y << 4;
+				int y1 = (y + q - 1) << 4;
+				// Evaluate half-space functions
+				bool a00 = C1 + DX12 * y0 - DY12 * x0 > 0;
+				bool a10 = C1 + DX12 * y0 - DY12 * x1 > 0;
+				bool a01 = C1 + DX12 * y1 - DY12 * x0 > 0;
+				bool a11 = C1 + DX12 * y1 - DY12 * x1 > 0;
+				int a = (a00 << 0) | (a10 << 1) | (a01 << 2) | (a11 << 3);
+				bool b00 = C2 + DX23 * y0 - DY23 * x0 > 0;
+				bool b10 = C2 + DX23 * y0 - DY23 * x1 > 0;
+				bool b01 = C2 + DX23 * y1 - DY23 * x0 > 0;
+				bool b11 = C2 + DX23 * y1 - DY23 * x1 > 0;
+				int b = (b00 << 0) | (b10 << 1) | (b01 << 2) | (b11 << 3);
+				bool c00 = C3 + DX31 * y0 - DY31 * x0 > 0;
+				bool c10 = C3 + DX31 * y0 - DY31 * x1 > 0;
+				bool c01 = C3 + DX31 * y1 - DY31 * x0 > 0;
+				bool c11 = C3 + DX31 * y1 - DY31 * x1 > 0;
+				int c = (c00 << 0) | (c10 << 1) | (c01 << 2) | (c11 << 3);
+				// Skip block when outside an edge
+				if (a == 0x0 || b == 0x0 || c == 0x0) continue;
+				// Accept whole block when totally covered
+				if (a == 0xF && b == 0xF && c == 0xF) {
+					Vector3 texRow = t1 + dy * (y0 - v1.y) + dx * (x0 - v1.x);
+					for (int iy = y; iy < y + q; iy++) {
+						Vector3 tex = texRow;
+						for (int ix = x; ix < x + q; ix++) {
+							//Vector3 tex = t1 + dx * (ix - v1.x) + dy * (iy - v1.y);
+							if (!cb(param, ix, iy, tex, dx, dy, 1.0)) {
+								// early out.
+								return false;
+							}
+							tex += dx;
+						}
+						texRow += dy;
+					}
+				} else { // Partially covered block
+					int CY1 = C1 + DX12 * y0 - DY12 * x0;
+					int CY2 = C2 + DX23 * y0 - DY23 * x0;
+					int CY3 = C3 + DX31 * y0 - DY31 * x0;
+					Vector3 texRow = t1 + dy * (y0 - v1.y) + dx * (x0 - v1.x);
+					for (int iy = y; iy < y + q; iy++) {
+						int CX1 = CY1;
+						int CX2 = CY2;
+						int CX3 = CY3;
+						Vector3 tex = texRow;
+						for (int ix = x; ix < x + q; ix++) {
+							if (CX1 > 0 && CX2 > 0 && CX3 > 0) {
+								if (!cb(param, ix, iy, tex, dx, dy, 1.0)) {
+									// early out.
+									return false;
+								}
+							}
+							CX1 -= FDY12;
+							CX2 -= FDY23;
+							CX3 -= FDY31;
+							tex += dx;
+						}
+						CY1 += FDX12;
+						CY2 += FDX23;
+						CY3 += FDX31;
+						texRow += dy;
+					}
+				}
+			}
+		}
+		return true;
+	}
+
+	// extents has to be multiple of BK_SIZE!!
+	bool drawAA(const Vector2 &extents, bool enableScissors, SamplingCallback cb, void *param) {
+		const float PX_INSIDE = 1.0f / sqrt(2.0f);
+		const float PX_OUTSIDE = -1.0f / sqrt(2.0f);
+		const float BK_SIZE = 8;
+		const float BK_INSIDE = sqrt(BK_SIZE * BK_SIZE / 2.0f);
+		const float BK_OUTSIDE = -sqrt(BK_SIZE * BK_SIZE / 2.0f);
+
+		float minx, miny, maxx, maxy;
+		if (enableScissors) {
+			// Bounding rectangle
+			minx = floorf(std::max(min3(v1.x, v2.x, v3.x), 0.0f));
+			miny = floorf(std::max(min3(v1.y, v2.y, v3.y), 0.0f));
+			maxx = ceilf(std::min(max3(v1.x, v2.x, v3.x), extents.x - 1.0f));
+			maxy = ceilf(std::min(max3(v1.y, v2.y, v3.y), extents.y - 1.0f));
+		} else {
+			// Bounding rectangle
+			minx = floorf(min3(v1.x, v2.x, v3.x));
+			miny = floorf(min3(v1.y, v2.y, v3.y));
+			maxx = ceilf(max3(v1.x, v2.x, v3.x));
+			maxy = ceilf(max3(v1.y, v2.y, v3.y));
+		}
+		// There's no reason to align the blocks to the viewport, instead we align them to the origin of the triangle bounds.
+		minx = floorf(minx);
+		miny = floorf(miny);
+		//minx = (float)(((int)minx) & (~((int)BK_SIZE - 1))); // align to blocksize (we don't need to worry about blocks partially out of viewport)
+		//miny = (float)(((int)miny) & (~((int)BK_SIZE - 1)));
+		minx += 0.5;
+		miny += 0.5; // sampling at texel centers!
+		maxx += 0.5;
+		maxy += 0.5;
+		// Half-edge constants
+		float C1 = n1.x * (-v1.x) + n1.y * (-v1.y);
+		float C2 = n2.x * (-v2.x) + n2.y * (-v2.y);
+		float C3 = n3.x * (-v3.x) + n3.y * (-v3.y);
+		// Loop through blocks
+		for (float y0 = miny; y0 <= maxy; y0 += BK_SIZE) {
+			for (float x0 = minx; x0 <= maxx; x0 += BK_SIZE) {
+				// Corners of block
+				float xc = (x0 + (BK_SIZE - 1) / 2.0f);
+				float yc = (y0 + (BK_SIZE - 1) / 2.0f);
+				// Evaluate half-space functions
+				float aC = C1 + n1.x * xc + n1.y * yc;
+				float bC = C2 + n2.x * xc + n2.y * yc;
+				float cC = C3 + n3.x * xc + n3.y * yc;
+				// Skip block when outside an edge
+				if ((aC <= BK_OUTSIDE) || (bC <= BK_OUTSIDE) || (cC <= BK_OUTSIDE)) continue;
+				// Accept whole block when totally covered
+				if ((aC >= BK_INSIDE) && (bC >= BK_INSIDE) && (cC >= BK_INSIDE)) {
+					Vector3 texRow = t1 + dy * (y0 - v1.y) + dx * (x0 - v1.x);
+					for (float y = y0; y < y0 + BK_SIZE; y++) {
+						Vector3 tex = texRow;
+						for (float x = x0; x < x0 + BK_SIZE; x++) {
+							if (!cb(param, (int)x, (int)y, tex, dx, dy, 1.0f)) {
+								return false;
+							}
+							tex += dx;
+						}
+						texRow += dy;
+					}
+				} else { // Partially covered block
+					float CY1 = C1 + n1.x * x0 + n1.y * y0;
+					float CY2 = C2 + n2.x * x0 + n2.y * y0;
+					float CY3 = C3 + n3.x * x0 + n3.y * y0;
+					Vector3 texRow = t1 + dy * (y0 - v1.y) + dx * (x0 - v1.x);
+					for (float y = y0; y < y0 + BK_SIZE; y++) { // @@ This is not clipping to scissor rectangle correctly.
+						float CX1 = CY1;
+						float CX2 = CY2;
+						float CX3 = CY3;
+						Vector3 tex = texRow;
+						for (float x = x0; x < x0 + BK_SIZE; x++) { // @@ This is not clipping to scissor rectangle correctly.
+							if (CX1 >= PX_INSIDE && CX2 >= PX_INSIDE && CX3 >= PX_INSIDE) {
+								// pixel completely covered
+								Vector3 tex2 = t1 + dx * (x - v1.x) + dy * (y - v1.y);
+								if (!cb(param, (int)x, (int)y, tex2, dx, dy, 1.0f)) {
+									return false;
+								}
+							} else if ((CX1 >= PX_OUTSIDE) && (CX2 >= PX_OUTSIDE) && (CX3 >= PX_OUTSIDE)) {
+								// triangle partially covers pixel. do clipping.
+								ClippedTriangle ct(v1 - Vector2(x, y), v2 - Vector2(x, y), v3 - Vector2(x, y));
+								ct.clipAABox(-0.5, -0.5, 0.5, 0.5);
+								Vector2 centroid = ct.centroid();
+								float area = ct.area();
+								if (area > 0.0f) {
+									Vector3 texCent = tex - dx * centroid.x - dy * centroid.y;
+									//xaAssert(texCent.x >= -0.1f && texCent.x <= 1.1f); // @@ Centroid is not very exact...
+									//xaAssert(texCent.y >= -0.1f && texCent.y <= 1.1f);
+									//xaAssert(texCent.z >= -0.1f && texCent.z <= 1.1f);
+									//Vector3 texCent2 = t1 + dx * (x - v1.x) + dy * (y - v1.y);
+									if (!cb(param, (int)x, (int)y, texCent, dx, dy, area)) {
+										return false;
+									}
+								}
+							}
+							CX1 += n1.x;
+							CX2 += n2.x;
+							CX3 += n3.x;
+							tex += dx;
+						}
+						CY1 += n1.y;
+						CY2 += n2.y;
+						CY3 += n3.y;
+						texRow += dy;
+					}
+				}
+			}
+		}
+		return true;
+	}
+
+	void flipBackface() {
+		// check if triangle is backfacing, if so, swap two vertices
+		if (((v3.x - v1.x) * (v2.y - v1.y) - (v3.y - v1.y) * (v2.x - v1.x)) < 0) {
+			Vector2 hv = v1;
+			v1 = v2;
+			v2 = hv; // swap pos
+			Vector3 ht = t1;
+			t1 = t2;
+			t2 = ht; // swap tex
+		}
+	}
+
+	// compute unit inward normals for each edge.
+	void computeUnitInwardNormals() {
+		n1 = v1 - v2;
+		n1 = Vector2(-n1.y, n1.x);
+		n1 = n1 * (1.0f / sqrtf(n1.x * n1.x + n1.y * n1.y));
+		n2 = v2 - v3;
+		n2 = Vector2(-n2.y, n2.x);
+		n2 = n2 * (1.0f / sqrtf(n2.x * n2.x + n2.y * n2.y));
+		n3 = v3 - v1;
+		n3 = Vector2(-n3.y, n3.x);
+		n3 = n3 * (1.0f / sqrtf(n3.x * n3.x + n3.y * n3.y));
+	}
+
+	// Vertices.
+	Vector2 v1, v2, v3;
+	Vector2 n1, n2, n3; // unit inward normals
+	Vector3 t1, t2, t3;
+
+	// Deltas.
+	Vector3 dx, dy;
+
+	float sign;
+	bool valid;
+};
+
+enum Mode {
+	Mode_Nearest,
+	Mode_Antialiased
+};
+
+// Process the given triangle. Returns false if rasterization was interrupted by the callback.
+static bool drawTriangle(Mode mode, Vector2::Arg extents, bool enableScissors, const Vector2 v[3], SamplingCallback cb, void *param) {
+	Triangle tri(v[0], v[1], v[2], Vector3(1, 0, 0), Vector3(0, 1, 0), Vector3(0, 0, 1));
+	// @@ It would be nice to have a conservative drawing mode that enlarges the triangle extents by one texel and is able to handle degenerate triangles.
+	// @@ Maybe the simplest thing to do would be raster triangle edges.
+	if (tri.valid) {
+		if (mode == Mode_Antialiased) {
+			return tri.drawAA(extents, enableScissors, cb, param);
+		}
+		if (mode == Mode_Nearest) {
+			return tri.draw(extents, enableScissors, cb, param);
+		}
+	}
+	return true;
+}
+
+// Process the given quad. Returns false if rasterization was interrupted by the callback.
+static bool drawQuad(Mode mode, Vector2::Arg extents, bool enableScissors, const Vector2 v[4], SamplingCallback cb, void *param) {
+	bool sign0 = triangleArea2(v[0], v[1], v[2]) > 0.0f;
+	bool sign1 = triangleArea2(v[0], v[2], v[3]) > 0.0f;
+	// Divide the quad into two non overlapping triangles.
+	if (sign0 == sign1) {
+		Triangle tri0(v[0], v[1], v[2], Vector3(0, 0, 0), Vector3(1, 0, 0), Vector3(1, 1, 0));
+		Triangle tri1(v[0], v[2], v[3], Vector3(0, 0, 0), Vector3(1, 1, 0), Vector3(0, 1, 0));
+		if (tri0.valid && tri1.valid) {
+			if (mode == Mode_Antialiased) {
+				return tri0.drawAA(extents, enableScissors, cb, param) && tri1.drawAA(extents, enableScissors, cb, param);
+			} else {
+				return tri0.draw(extents, enableScissors, cb, param) && tri1.draw(extents, enableScissors, cb, param);
+			}
+		}
+	} else {
+		Triangle tri0(v[0], v[1], v[3], Vector3(0, 0, 0), Vector3(1, 0, 0), Vector3(0, 1, 0));
+		Triangle tri1(v[1], v[2], v[3], Vector3(1, 0, 0), Vector3(1, 1, 0), Vector3(0, 1, 0));
+		if (tri0.valid && tri1.valid) {
+			if (mode == Mode_Antialiased) {
+				return tri0.drawAA(extents, enableScissors, cb, param) && tri1.drawAA(extents, enableScissors, cb, param);
+			} else {
+				return tri0.draw(extents, enableScissors, cb, param) && tri1.draw(extents, enableScissors, cb, param);
+			}
+		}
+	}
+	return true;
+}
+} // namespace raster
+
+// Full and sparse vector and matrix classes. BLAS subset.
+// Pseudo-BLAS interface.
+namespace sparse {
+enum Transpose {
+	NoTransposed = 0,
+	Transposed = 1
+};
+
+/**
+* Sparse matrix class. The matrix is assumed to be sparse and to have
+* very few non-zero elements, for this reason it's stored in indexed
+* format. To multiply column vectors efficiently, the matrix stores
+* the elements in indexed-column order, there is a list of indexed
+* elements for each row of the matrix. As with the FullVector the
+* dimension of the matrix is constant.
+**/
+class Matrix {
+public:
+	// An element of the sparse array.
+	struct Coefficient {
+		uint32_t x; // column
+		float v; // value
+	};
+
+	Matrix(uint32_t d) :
+			m_width(d) { m_array.resize(d); }
+	Matrix(uint32_t w, uint32_t h) :
+			m_width(w) { m_array.resize(h); }
+	Matrix(const Matrix &m) :
+			m_width(m.m_width) { m_array = m.m_array; }
+
+	const Matrix &operator=(const Matrix &m) {
+		xaAssert(width() == m.width());
+		xaAssert(height() == m.height());
+		m_array = m.m_array;
+		return *this;
+	}
+
+	uint32_t width() const { return m_width; }
+	uint32_t height() const { return m_array.size(); }
+	bool isSquare() const { return width() == height(); }
+
+	// x is column, y is row
+	float getCoefficient(uint32_t x, uint32_t y) const {
+		xaDebugAssert(x < width());
+		xaDebugAssert(y < height());
+		const uint32_t count = m_array[y].size();
+		for (uint32_t i = 0; i < count; i++) {
+			if (m_array[y][i].x == x) return m_array[y][i].v;
+		}
+		return 0.0f;
+	}
+
+	void setCoefficient(uint32_t x, uint32_t y, float f) {
+		xaDebugAssert(x < width());
+		xaDebugAssert(y < height());
+		const uint32_t count = m_array[y].size();
+		for (uint32_t i = 0; i < count; i++) {
+			if (m_array[y][i].x == x) {
+				m_array[y][i].v = f;
+				return;
+			}
+		}
+		if (f != 0.0f) {
+			Coefficient c = { x, f };
+			m_array[y].push_back(c);
+		}
+	}
+
+	float dotRow(uint32_t y, const FullVector &v) const {
+		xaDebugAssert(y < height());
+		const uint32_t count = m_array[y].size();
+		float sum = 0;
+		for (uint32_t i = 0; i < count; i++) {
+			sum += m_array[y][i].v * v[m_array[y][i].x];
+		}
+		return sum;
+	}
+
+	void madRow(uint32_t y, float alpha, FullVector &v) const {
+		xaDebugAssert(y < height());
+		const uint32_t count = m_array[y].size();
+		for (uint32_t i = 0; i < count; i++) {
+			v[m_array[y][i].x] += alpha * m_array[y][i].v;
+		}
+	}
+
+	void clearRow(uint32_t y) {
+		xaDebugAssert(y < height());
+		m_array[y].clear();
+	}
+
+	void scaleRow(uint32_t y, float f) {
+		xaDebugAssert(y < height());
+		const uint32_t count = m_array[y].size();
+		for (uint32_t i = 0; i < count; i++) {
+			m_array[y][i].v *= f;
+		}
+	}
+
+	const std::vector<Coefficient> &getRow(uint32_t y) const { return m_array[y]; }
+
+private:
+	/// Number of columns.
+	const uint32_t m_width;
+
+	/// Array of matrix elements.
+	std::vector<std::vector<Coefficient> > m_array;
+};
+
+// y = a * x + y
+static void saxpy(float a, const FullVector &x, FullVector &y) {
+	xaDebugAssert(x.dimension() == y.dimension());
+	const uint32_t dim = x.dimension();
+	for (uint32_t i = 0; i < dim; i++) {
+		y[i] += a * x[i];
+	}
+}
+
+static void copy(const FullVector &x, FullVector &y) {
+	xaDebugAssert(x.dimension() == y.dimension());
+	const uint32_t dim = x.dimension();
+	for (uint32_t i = 0; i < dim; i++) {
+		y[i] = x[i];
+	}
+}
+
+static void scal(float a, FullVector &x) {
+	const uint32_t dim = x.dimension();
+	for (uint32_t i = 0; i < dim; i++) {
+		x[i] *= a;
+	}
+}
+
+static float dot(const FullVector &x, const FullVector &y) {
+	xaDebugAssert(x.dimension() == y.dimension());
+	const uint32_t dim = x.dimension();
+	float sum = 0;
+	for (uint32_t i = 0; i < dim; i++) {
+		sum += x[i] * y[i];
+	}
+	return sum;
+}
+
+static void mult(Transpose TM, const Matrix &M, const FullVector &x, FullVector &y) {
+	const uint32_t w = M.width();
+	const uint32_t h = M.height();
+	if (TM == Transposed) {
+		xaDebugAssert(h == x.dimension());
+		xaDebugAssert(w == y.dimension());
+		y.fill(0.0f);
+		for (uint32_t i = 0; i < h; i++) {
+			M.madRow(i, x[i], y);
+		}
+	} else {
+		xaDebugAssert(w == x.dimension());
+		xaDebugAssert(h == y.dimension());
+		for (uint32_t i = 0; i < h; i++) {
+			y[i] = M.dotRow(i, x);
+		}
+	}
+}
+
+// y = M * x
+static void mult(const Matrix &M, const FullVector &x, FullVector &y) {
+	mult(NoTransposed, M, x, y);
+}
+
+static void sgemv(float alpha, Transpose TA, const Matrix &A, const FullVector &x, float beta, FullVector &y) {
+	const uint32_t w = A.width();
+	const uint32_t h = A.height();
+	if (TA == Transposed) {
+		xaDebugAssert(h == x.dimension());
+		xaDebugAssert(w == y.dimension());
+		for (uint32_t i = 0; i < h; i++) {
+			A.madRow(i, alpha * x[i], y);
+		}
+	} else {
+		xaDebugAssert(w == x.dimension());
+		xaDebugAssert(h == y.dimension());
+		for (uint32_t i = 0; i < h; i++) {
+			y[i] = alpha * A.dotRow(i, x) + beta * y[i];
+		}
+	}
+}
+
+// y = alpha*A*x + beta*y
+static void sgemv(float alpha, const Matrix &A, const FullVector &x, float beta, FullVector &y) {
+	sgemv(alpha, NoTransposed, A, x, beta, y);
+}
+
+// dot y-row of A by x-column of B
+static float dotRowColumn(int y, const Matrix &A, int x, const Matrix &B) {
+	const std::vector<Matrix::Coefficient> &row = A.getRow(y);
+	const uint32_t count = row.size();
+	float sum = 0.0f;
+	for (uint32_t i = 0; i < count; i++) {
+		const Matrix::Coefficient &c = row[i];
+		sum += c.v * B.getCoefficient(x, c.x);
+	}
+	return sum;
+}
+
+// dot y-row of A by x-row of B
+static float dotRowRow(int y, const Matrix &A, int x, const Matrix &B) {
+	const std::vector<Matrix::Coefficient> &row = A.getRow(y);
+	const uint32_t count = row.size();
+	float sum = 0.0f;
+	for (uint32_t i = 0; i < count; i++) {
+		const Matrix::Coefficient &c = row[i];
+		sum += c.v * B.getCoefficient(c.x, x);
+	}
+	return sum;
+}
+
+// dot y-column of A by x-column of B
+static float dotColumnColumn(int y, const Matrix &A, int x, const Matrix &B) {
+	xaDebugAssert(A.height() == B.height());
+	const uint32_t h = A.height();
+	float sum = 0.0f;
+	for (uint32_t i = 0; i < h; i++) {
+		sum += A.getCoefficient(y, i) * B.getCoefficient(x, i);
+	}
+	return sum;
+}
+
+static void transpose(const Matrix &A, Matrix &B) {
+	xaDebugAssert(A.width() == B.height());
+	xaDebugAssert(B.width() == A.height());
+	const uint32_t w = A.width();
+	for (uint32_t x = 0; x < w; x++) {
+		B.clearRow(x);
+	}
+	const uint32_t h = A.height();
+	for (uint32_t y = 0; y < h; y++) {
+		const std::vector<Matrix::Coefficient> &row = A.getRow(y);
+		const uint32_t count = row.size();
+		for (uint32_t i = 0; i < count; i++) {
+			const Matrix::Coefficient &c = row[i];
+			xaDebugAssert(c.x < w);
+			B.setCoefficient(y, c.x, c.v);
+		}
+	}
+}
+
+static void sgemm(float alpha, Transpose TA, const Matrix &A, Transpose TB, const Matrix &B, float beta, Matrix &C) {
+	const uint32_t w = C.width();
+	const uint32_t h = C.height();
+	uint32_t aw = (TA == NoTransposed) ? A.width() : A.height();
+	uint32_t ah = (TA == NoTransposed) ? A.height() : A.width();
+	uint32_t bw = (TB == NoTransposed) ? B.width() : B.height();
+	uint32_t bh = (TB == NoTransposed) ? B.height() : B.width();
+	xaDebugAssert(aw == bh);
+	xaDebugAssert(bw == ah);
+	xaDebugAssert(w == bw);
+	xaDebugAssert(h == ah);
+#ifdef NDEBUG
+	aw = ah = bw = bh = 0; // silence unused parameter warning
+#endif
+	for (uint32_t y = 0; y < h; y++) {
+		for (uint32_t x = 0; x < w; x++) {
+			float c = beta * C.getCoefficient(x, y);
+			if (TA == NoTransposed && TB == NoTransposed) {
+				// dot y-row of A by x-column of B.
+				c += alpha * dotRowColumn(y, A, x, B);
+			} else if (TA == Transposed && TB == Transposed) {
+				// dot y-column of A by x-row of B.
+				c += alpha * dotRowColumn(x, B, y, A);
+			} else if (TA == Transposed && TB == NoTransposed) {
+				// dot y-column of A by x-column of B.
+				c += alpha * dotColumnColumn(y, A, x, B);
+			} else if (TA == NoTransposed && TB == Transposed) {
+				// dot y-row of A by x-row of B.
+				c += alpha * dotRowRow(y, A, x, B);
+			}
+			C.setCoefficient(x, y, c);
+		}
+	}
+}
+
+static void mult(Transpose TA, const Matrix &A, Transpose TB, const Matrix &B, Matrix &C) {
+	sgemm(1.0f, TA, A, TB, B, 0.0f, C);
+}
+
+// C = A * B
+static void mult(const Matrix &A, const Matrix &B, Matrix &C) {
+	mult(NoTransposed, A, NoTransposed, B, C);
+}
+
+} // namespace sparse
+
+class JacobiPreconditioner {
+public:
+	JacobiPreconditioner(const sparse::Matrix &M, bool symmetric) :
+			m_inverseDiagonal(M.width()) {
+		xaAssert(M.isSquare());
+		for (uint32_t x = 0; x < M.width(); x++) {
+			float elem = M.getCoefficient(x, x);
+			//xaDebugAssert( elem != 0.0f ); // This can be zero in the presence of zero area triangles.
+			if (symmetric) {
+				m_inverseDiagonal[x] = (elem != 0) ? 1.0f / sqrtf(fabsf(elem)) : 1.0f;
+			} else {
+				m_inverseDiagonal[x] = (elem != 0) ? 1.0f / elem : 1.0f;
+			}
+		}
+	}
+
+	void apply(const FullVector &x, FullVector &y) const {
+		xaDebugAssert(x.dimension() == m_inverseDiagonal.dimension());
+		xaDebugAssert(y.dimension() == m_inverseDiagonal.dimension());
+		// @@ Wrap vector component-wise product into a separate function.
+		const uint32_t D = x.dimension();
+		for (uint32_t i = 0; i < D; i++) {
+			y[i] = m_inverseDiagonal[i] * x[i];
+		}
+	}
+
+private:
+	FullVector m_inverseDiagonal;
+};
+
+// Linear solvers.
+class Solver {
+public:
+	// Solve the symmetric system: At·A·x = At·b
+	static bool LeastSquaresSolver(const sparse::Matrix &A, const FullVector &b, FullVector &x, float epsilon = 1e-5f) {
+		xaDebugAssert(A.width() == x.dimension());
+		xaDebugAssert(A.height() == b.dimension());
+		xaDebugAssert(A.height() >= A.width()); // @@ If height == width we could solve it directly...
+		const uint32_t D = A.width();
+		sparse::Matrix At(A.height(), A.width());
+		sparse::transpose(A, At);
+		FullVector Atb(D);
+		sparse::mult(At, b, Atb);
+		sparse::Matrix AtA(D);
+		sparse::mult(At, A, AtA);
+		return SymmetricSolver(AtA, Atb, x, epsilon);
+	}
+
+	// See section 10.4.3 in: Mesh Parameterization: Theory and Practice, Siggraph Course Notes, August 2007
+	static bool LeastSquaresSolver(const sparse::Matrix &A, const FullVector &b, FullVector &x, const uint32_t *lockedParameters, uint32_t lockedCount, float epsilon = 1e-5f) {
+		xaDebugAssert(A.width() == x.dimension());
+		xaDebugAssert(A.height() == b.dimension());
+		xaDebugAssert(A.height() >= A.width() - lockedCount);
+		// @@ This is not the most efficient way of building a system with reduced degrees of freedom. It would be faster to do it on the fly.
+		const uint32_t D = A.width() - lockedCount;
+		xaDebugAssert(D > 0);
+		// Compute: b - Al * xl
+		FullVector b_Alxl(b);
+		for (uint32_t y = 0; y < A.height(); y++) {
+			const uint32_t count = A.getRow(y).size();
+			for (uint32_t e = 0; e < count; e++) {
+				uint32_t column = A.getRow(y)[e].x;
+				bool isFree = true;
+				for (uint32_t i = 0; i < lockedCount; i++) {
+					isFree &= (lockedParameters[i] != column);
+				}
+				if (!isFree) {
+					b_Alxl[y] -= x[column] * A.getRow(y)[e].v;
+				}
+			}
+		}
+		// Remove locked columns from A.
+		sparse::Matrix Af(D, A.height());
+		for (uint32_t y = 0; y < A.height(); y++) {
+			const uint32_t count = A.getRow(y).size();
+			for (uint32_t e = 0; e < count; e++) {
+				uint32_t column = A.getRow(y)[e].x;
+				uint32_t ix = column;
+				bool isFree = true;
+				for (uint32_t i = 0; i < lockedCount; i++) {
+					isFree &= (lockedParameters[i] != column);
+					if (column > lockedParameters[i]) ix--; // shift columns
+				}
+				if (isFree) {
+					Af.setCoefficient(ix, y, A.getRow(y)[e].v);
+				}
+			}
+		}
+		// Remove elements from x
+		FullVector xf(D);
+		for (uint32_t i = 0, j = 0; i < A.width(); i++) {
+			bool isFree = true;
+			for (uint32_t l = 0; l < lockedCount; l++) {
+				isFree &= (lockedParameters[l] != i);
+			}
+			if (isFree) {
+				xf[j++] = x[i];
+			}
+		}
+		// Solve reduced system.
+		bool result = LeastSquaresSolver(Af, b_Alxl, xf, epsilon);
+		// Copy results back to x.
+		for (uint32_t i = 0, j = 0; i < A.width(); i++) {
+			bool isFree = true;
+			for (uint32_t l = 0; l < lockedCount; l++) {
+				isFree &= (lockedParameters[l] != i);
+			}
+			if (isFree) {
+				x[i] = xf[j++];
+			}
+		}
+		return result;
+	}
+
+private:
+	/**
+	* Compute the solution of the sparse linear system Ab=x using the Conjugate
+	* Gradient method.
+	*
+	* Solving sparse linear systems:
+	* (1)		A·x = b
+	*
+	* The conjugate gradient algorithm solves (1) only in the case that A is
+	* symmetric and positive definite. It is based on the idea of minimizing the
+	* function
+	*
+	* (2)		f(x) = 1/2·x·A·x - b·x
+	*
+	* This function is minimized when its gradient
+	*
+	* (3)		df = A·x - b
+	*
+	* is zero, which is equivalent to (1). The minimization is carried out by
+	* generating a succession of search directions p.k and improved minimizers x.k.
+	* At each stage a quantity alfa.k is found that minimizes f(x.k + alfa.k·p.k),
+	* and x.k+1 is set equal to the new point x.k + alfa.k·p.k. The p.k and x.k are
+	* built up in such a way that x.k+1 is also the minimizer of f over the whole
+	* vector space of directions already taken, {p.1, p.2, . . . , p.k}. After N
+	* iterations you arrive at the minimizer over the entire vector space, i.e., the
+	* solution to (1).
+	*
+	* For a really good explanation of the method see:
+	*
+	* "An Introduction to the Conjugate Gradient Method Without the Agonizing Pain",
+	* Jonhathan Richard Shewchuk.
+	*
+	**/
+	static bool ConjugateGradientSolver(const sparse::Matrix &A, const FullVector &b, FullVector &x, float epsilon) {
+		xaDebugAssert(A.isSquare());
+		xaDebugAssert(A.width() == b.dimension());
+		xaDebugAssert(A.width() == x.dimension());
+		int i = 0;
+		const int D = A.width();
+		const int i_max = 4 * D; // Convergence should be linear, but in some cases, it's not.
+		FullVector r(D); // residual
+		FullVector p(D); // search direction
+		FullVector q(D); //
+		float delta_0;
+		float delta_old;
+		float delta_new;
+		float alpha;
+		float beta;
+		// r = b - A·x;
+		sparse::copy(b, r);
+		sparse::sgemv(-1, A, x, 1, r);
+		// p = r;
+		sparse::copy(r, p);
+		delta_new = sparse::dot(r, r);
+		delta_0 = delta_new;
+		while (i < i_max && delta_new > epsilon * epsilon * delta_0) {
+			i++;
+			// q = A·p
+			mult(A, p, q);
+			// alpha = delta_new / p·q
+			alpha = delta_new / sparse::dot(p, q);
+			// x = alfa·p + x
+			sparse::saxpy(alpha, p, x);
+			if ((i & 31) == 0) { // recompute r after 32 steps
+				// r = b - A·x
+				sparse::copy(b, r);
+				sparse::sgemv(-1, A, x, 1, r);
+			} else {
+				// r = r - alpha·q
+				sparse::saxpy(-alpha, q, r);
+			}
+			delta_old = delta_new;
+			delta_new = sparse::dot(r, r);
+			beta = delta_new / delta_old;
+			// p = beta·p + r
+			sparse::scal(beta, p);
+			sparse::saxpy(1, r, p);
+		}
+		return delta_new <= epsilon * epsilon * delta_0;
+	}
+
+	// Conjugate gradient with preconditioner.
+	static bool ConjugateGradientSolver(const JacobiPreconditioner &preconditioner, const sparse::Matrix &A, const FullVector &b, FullVector &x, float epsilon) {
+		xaDebugAssert(A.isSquare());
+		xaDebugAssert(A.width() == b.dimension());
+		xaDebugAssert(A.width() == x.dimension());
+		int i = 0;
+		const int D = A.width();
+		const int i_max = 4 * D; // Convergence should be linear, but in some cases, it's not.
+		FullVector r(D); // residual
+		FullVector p(D); // search direction
+		FullVector q(D); //
+		FullVector s(D); // preconditioned
+		float delta_0;
+		float delta_old;
+		float delta_new;
+		float alpha;
+		float beta;
+		// r = b - A·x
+		sparse::copy(b, r);
+		sparse::sgemv(-1, A, x, 1, r);
+		// p = M^-1 · r
+		preconditioner.apply(r, p);
+		delta_new = sparse::dot(r, p);
+		delta_0 = delta_new;
+		while (i < i_max && delta_new > epsilon * epsilon * delta_0) {
+			i++;
+			// q = A·p
+			mult(A, p, q);
+			// alpha = delta_new / p·q
+			alpha = delta_new / sparse::dot(p, q);
+			// x = alfa·p + x
+			sparse::saxpy(alpha, p, x);
+			if ((i & 31) == 0) { // recompute r after 32 steps
+				// r = b - A·x
+				sparse::copy(b, r);
+				sparse::sgemv(-1, A, x, 1, r);
+			} else {
+				// r = r - alfa·q
+				sparse::saxpy(-alpha, q, r);
+			}
+			// s = M^-1 · r
+			preconditioner.apply(r, s);
+			delta_old = delta_new;
+			delta_new = sparse::dot(r, s);
+			beta = delta_new / delta_old;
+			// p = s + beta·p
+			sparse::scal(beta, p);
+			sparse::saxpy(1, s, p);
+		}
+		return delta_new <= epsilon * epsilon * delta_0;
+	}
+
+	static bool SymmetricSolver(const sparse::Matrix &A, const FullVector &b, FullVector &x, float epsilon = 1e-5f) {
+		xaDebugAssert(A.height() == A.width());
+		xaDebugAssert(A.height() == b.dimension());
+		xaDebugAssert(b.dimension() == x.dimension());
+		JacobiPreconditioner jacobi(A, true);
+		return ConjugateGradientSolver(jacobi, A, b, x, epsilon);
+	}
+};
+
+namespace param {
+class Atlas;
+class Chart;
+
+// Fast sweep in 3 directions
+static bool findApproximateDiameterVertices(halfedge::Mesh *mesh, halfedge::Vertex **a, halfedge::Vertex **b) {
+	xaDebugAssert(mesh != NULL);
+	xaDebugAssert(a != NULL);
+	xaDebugAssert(b != NULL);
+	const uint32_t vertexCount = mesh->vertexCount();
+	halfedge::Vertex *minVertex[3];
+	halfedge::Vertex *maxVertex[3];
+	minVertex[0] = minVertex[1] = minVertex[2] = NULL;
+	maxVertex[0] = maxVertex[1] = maxVertex[2] = NULL;
+	for (uint32_t v = 1; v < vertexCount; v++) {
+		halfedge::Vertex *vertex = mesh->vertexAt(v);
+		xaDebugAssert(vertex != NULL);
+		if (vertex->isBoundary()) {
+			minVertex[0] = minVertex[1] = minVertex[2] = vertex;
+			maxVertex[0] = maxVertex[1] = maxVertex[2] = vertex;
+			break;
+		}
+	}
+	if (minVertex[0] == NULL) {
+		// Input mesh has not boundaries.
+		return false;
+	}
+	for (uint32_t v = 1; v < vertexCount; v++) {
+		halfedge::Vertex *vertex = mesh->vertexAt(v);
+		xaDebugAssert(vertex != NULL);
+		if (!vertex->isBoundary()) {
+			// Skip interior vertices.
+			continue;
+		}
+		if (vertex->pos.x < minVertex[0]->pos.x)
+			minVertex[0] = vertex;
+		else if (vertex->pos.x > maxVertex[0]->pos.x)
+			maxVertex[0] = vertex;
+		if (vertex->pos.y < minVertex[1]->pos.y)
+			minVertex[1] = vertex;
+		else if (vertex->pos.y > maxVertex[1]->pos.y)
+			maxVertex[1] = vertex;
+		if (vertex->pos.z < minVertex[2]->pos.z)
+			minVertex[2] = vertex;
+		else if (vertex->pos.z > maxVertex[2]->pos.z)
+			maxVertex[2] = vertex;
+	}
+	float lengths[3];
+	for (int i = 0; i < 3; i++) {
+		lengths[i] = length(minVertex[i]->pos - maxVertex[i]->pos);
+	}
+	if (lengths[0] > lengths[1] && lengths[0] > lengths[2]) {
+		*a = minVertex[0];
+		*b = maxVertex[0];
+	} else if (lengths[1] > lengths[2]) {
+		*a = minVertex[1];
+		*b = maxVertex[1];
+	} else {
+		*a = minVertex[2];
+		*b = maxVertex[2];
+	}
+	return true;
+}
+
+// Conformal relations from Brecht Van Lommel (based on ABF):
+
+static float vec_angle_cos(Vector3::Arg v1, Vector3::Arg v2, Vector3::Arg v3) {
+	Vector3 d1 = v1 - v2;
+	Vector3 d2 = v3 - v2;
+	return clamp(dot(d1, d2) / (length(d1) * length(d2)), -1.0f, 1.0f);
+}
+
+static float vec_angle(Vector3::Arg v1, Vector3::Arg v2, Vector3::Arg v3) {
+	float dot = vec_angle_cos(v1, v2, v3);
+	return acosf(dot);
+}
+
+static void triangle_angles(Vector3::Arg v1, Vector3::Arg v2, Vector3::Arg v3, float *a1, float *a2, float *a3) {
+	*a1 = vec_angle(v3, v1, v2);
+	*a2 = vec_angle(v1, v2, v3);
+	*a3 = PI - *a2 - *a1;
+}
+
+static void setup_abf_relations(sparse::Matrix &A, int row, const halfedge::Vertex *v0, const halfedge::Vertex *v1, const halfedge::Vertex *v2) {
+	int id0 = v0->id;
+	int id1 = v1->id;
+	int id2 = v2->id;
+	Vector3 p0 = v0->pos;
+	Vector3 p1 = v1->pos;
+	Vector3 p2 = v2->pos;
+	// @@ IC: Wouldn't it be more accurate to return cos and compute 1-cos^2?
+	// It does indeed seem to be a little bit more robust.
+	// @@ Need to revisit this more carefully!
+	float a0, a1, a2;
+	triangle_angles(p0, p1, p2, &a0, &a1, &a2);
+	float s0 = sinf(a0);
+	float s1 = sinf(a1);
+	float s2 = sinf(a2);
+	if (s1 > s0 && s1 > s2) {
+		std::swap(s1, s2);
+		std::swap(s0, s1);
+		std::swap(a1, a2);
+		std::swap(a0, a1);
+		std::swap(id1, id2);
+		std::swap(id0, id1);
+	} else if (s0 > s1 && s0 > s2) {
+		std::swap(s0, s2);
+		std::swap(s0, s1);
+		std::swap(a0, a2);
+		std::swap(a0, a1);
+		std::swap(id0, id2);
+		std::swap(id0, id1);
+	}
+	float c0 = cosf(a0);
+	float ratio = (s2 == 0.0f) ? 1.0f : s1 / s2;
+	float cosine = c0 * ratio;
+	float sine = s0 * ratio;
+	// Note  : 2*id + 0 --> u
+	//         2*id + 1 --> v
+	int u0_id = 2 * id0 + 0;
+	int v0_id = 2 * id0 + 1;
+	int u1_id = 2 * id1 + 0;
+	int v1_id = 2 * id1 + 1;
+	int u2_id = 2 * id2 + 0;
+	int v2_id = 2 * id2 + 1;
+	// Real part
+	A.setCoefficient(u0_id, 2 * row + 0, cosine - 1.0f);
+	A.setCoefficient(v0_id, 2 * row + 0, -sine);
+	A.setCoefficient(u1_id, 2 * row + 0, -cosine);
+	A.setCoefficient(v1_id, 2 * row + 0, sine);
+	A.setCoefficient(u2_id, 2 * row + 0, 1);
+	// Imaginary part
+	A.setCoefficient(u0_id, 2 * row + 1, sine);
+	A.setCoefficient(v0_id, 2 * row + 1, cosine - 1.0f);
+	A.setCoefficient(u1_id, 2 * row + 1, -sine);
+	A.setCoefficient(v1_id, 2 * row + 1, -cosine);
+	A.setCoefficient(v2_id, 2 * row + 1, 1);
+}
+
+bool computeLeastSquaresConformalMap(halfedge::Mesh *mesh) {
+	xaDebugAssert(mesh != NULL);
+	// For this to work properly, mesh should not have colocals that have the same
+	// attributes, unless you want the vertices to actually have different texcoords.
+	const uint32_t vertexCount = mesh->vertexCount();
+	const uint32_t D = 2 * vertexCount;
+	const uint32_t N = 2 * halfedge::countMeshTriangles(mesh);
+	// N is the number of equations (one per triangle)
+	// D is the number of variables (one per vertex; there are 2 pinned vertices).
+	if (N < D - 4) {
+		return false;
+	}
+	sparse::Matrix A(D, N);
+	FullVector b(N);
+	FullVector x(D);
+	// Fill b:
+	b.fill(0.0f);
+	// Fill x:
+	halfedge::Vertex *v0;
+	halfedge::Vertex *v1;
+	if (!findApproximateDiameterVertices(mesh, &v0, &v1)) {
+		// Mesh has no boundaries.
+		return false;
+	}
+	if (v0->tex == v1->tex) {
+		// LSCM expects an existing parameterization.
+		return false;
+	}
+	for (uint32_t v = 0; v < vertexCount; v++) {
+		halfedge::Vertex *vertex = mesh->vertexAt(v);
+		xaDebugAssert(vertex != NULL);
+		// Initial solution.
+		x[2 * v + 0] = vertex->tex.x;
+		x[2 * v + 1] = vertex->tex.y;
+	}
+	// Fill A:
+	const uint32_t faceCount = mesh->faceCount();
+	for (uint32_t f = 0, t = 0; f < faceCount; f++) {
+		const halfedge::Face *face = mesh->faceAt(f);
+		xaDebugAssert(face != NULL);
+		xaDebugAssert(face->edgeCount() == 3);
+		const halfedge::Vertex *vertex0 = NULL;
+		for (halfedge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) {
+			const halfedge::Edge *edge = it.current();
+			xaAssert(edge != NULL);
+			if (vertex0 == NULL) {
+				vertex0 = edge->vertex;
+			} else if (edge->next->vertex != vertex0) {
+				const halfedge::Vertex *vertex1 = edge->from();
+				const halfedge::Vertex *vertex2 = edge->to();
+				setup_abf_relations(A, t, vertex0, vertex1, vertex2);
+				//setup_conformal_map_relations(A, t, vertex0, vertex1, vertex2);
+				t++;
+			}
+		}
+	}
+	const uint32_t lockedParameters[] = {
+		2 * v0->id + 0,
+		2 * v0->id + 1,
+		2 * v1->id + 0,
+		2 * v1->id + 1
+	};
+	// Solve
+	Solver::LeastSquaresSolver(A, b, x, lockedParameters, 4, 0.000001f);
+	// Map x back to texcoords:
+	for (uint32_t v = 0; v < vertexCount; v++) {
+		halfedge::Vertex *vertex = mesh->vertexAt(v);
+		xaDebugAssert(vertex != NULL);
+		vertex->tex = Vector2(x[2 * v + 0], x[2 * v + 1]);
+	}
+	return true;
+}
+
+bool computeOrthogonalProjectionMap(halfedge::Mesh *mesh) {
+	Vector3 axis[2];
+	uint32_t vertexCount = mesh->vertexCount();
+	std::vector<Vector3> points(vertexCount);
+	points.resize(vertexCount);
+	for (uint32_t i = 0; i < vertexCount; i++) {
+		points[i] = mesh->vertexAt(i)->pos;
+	}
+	// Avoid redundant computations.
+	float matrix[6];
+	Fit::computeCovariance(vertexCount, points.data(), matrix);
+	if (matrix[0] == 0 && matrix[3] == 0 && matrix[5] == 0) {
+		return false;
+	}
+	float eigenValues[3];
+	Vector3 eigenVectors[3];
+	if (!Fit::eigenSolveSymmetric3(matrix, eigenValues, eigenVectors)) {
+		return false;
+	}
+	axis[0] = normalize(eigenVectors[0]);
+	axis[1] = normalize(eigenVectors[1]);
+	// Project vertices to plane.
+	for (halfedge::Mesh::VertexIterator it(mesh->vertices()); !it.isDone(); it.advance()) {
+		halfedge::Vertex *vertex = it.current();
+		vertex->tex.x = dot(axis[0], vertex->pos);
+		vertex->tex.y = dot(axis[1], vertex->pos);
+	}
+	return true;
+}
+
+void computeSingleFaceMap(halfedge::Mesh *mesh) {
+	xaDebugAssert(mesh != NULL);
+	xaDebugAssert(mesh->faceCount() == 1);
+	halfedge::Face *face = mesh->faceAt(0);
+	xaAssert(face != NULL);
+	Vector3 p0 = face->edge->from()->pos;
+	Vector3 p1 = face->edge->to()->pos;
+	Vector3 X = normalizeSafe(p1 - p0, Vector3(0.0f), 0.0f);
+	Vector3 Z = face->normal();
+	Vector3 Y = normalizeSafe(cross(Z, X), Vector3(0.0f), 0.0f);
+	uint32_t i = 0;
+	for (halfedge::Face::EdgeIterator it(face->edges()); !it.isDone(); it.advance(), i++) {
+		halfedge::Vertex *vertex = it.vertex();
+		xaAssert(vertex != NULL);
+		if (i == 0) {
+			vertex->tex = Vector2(0);
+		} else {
+			Vector3 pn = vertex->pos;
+			float xn = dot((pn - p0), X);
+			float yn = dot((pn - p0), Y);
+			vertex->tex = Vector2(xn, yn);
+		}
+	}
+}
+
+// Dummy implementation of a priority queue using sort at insertion.
+// - Insertion is o(n)
+// - Smallest element goes at the end, so that popping it is o(1).
+// - Resorting is n*log(n)
+// @@ Number of elements in the queue is usually small, and we'd have to rebalance often. I'm not sure it's worth implementing a heap.
+// @@ Searcing at removal would remove the need for sorting when priorities change.
+struct PriorityQueue {
+	PriorityQueue(uint32_t size = UINT_MAX) :
+			maxSize(size) {}
+
+	void push(float priority, uint32_t face) {
+		uint32_t i = 0;
+		const uint32_t count = pairs.size();
+		for (; i < count; i++) {
+			if (pairs[i].priority > priority) break;
+		}
+		Pair p = { priority, face };
+		pairs.insert(pairs.begin() + i, p);
+		if (pairs.size() > maxSize) {
+			pairs.erase(pairs.begin());
+		}
+	}
+
+	// push face out of order, to be sorted later.
+	void push(uint32_t face) {
+		Pair p = { 0.0f, face };
+		pairs.push_back(p);
+	}
+
+	uint32_t pop() {
+		uint32_t f = pairs.back().face;
+		pairs.pop_back();
+		return f;
+	}
+
+	void sort() {
+		//sort(pairs); // @@ My intro sort appears to be much slower than it should!
+		std::sort(pairs.begin(), pairs.end());
+	}
+
+	void clear() {
+		pairs.clear();
+	}
+
+	uint32_t count() const {
+		return pairs.size();
+	}
+
+	float firstPriority() const {
+		return pairs.back().priority;
+	}
+
+	const uint32_t maxSize;
+
+	struct Pair {
+		bool operator<(const Pair &p) const {
+			return priority > p.priority; // !! Sort in inverse priority order!
+		}
+
+		float priority;
+		uint32_t face;
+	};
+
+	std::vector<Pair> pairs;
+};
+
+struct ChartBuildData {
+	ChartBuildData(int p_id) :
+			id(p_id) {
+		planeNormal = Vector3(0);
+		centroid = Vector3(0);
+		coneAxis = Vector3(0);
+		coneAngle = 0;
+		area = 0;
+		boundaryLength = 0;
+		normalSum = Vector3(0);
+		centroidSum = Vector3(0);
+	}
+
+	int id;
+
+	// Proxy info:
+	Vector3 planeNormal;
+	Vector3 centroid;
+	Vector3 coneAxis;
+	float coneAngle;
+
+	float area;
+	float boundaryLength;
+	Vector3 normalSum;
+	Vector3 centroidSum;
+
+	std::vector<uint32_t> seeds; // @@ These could be a pointers to the halfedge faces directly.
+	std::vector<uint32_t> faces;
+	PriorityQueue candidates;
+};
+
+struct AtlasBuilder {
+	AtlasBuilder(const halfedge::Mesh *m) :
+			mesh(m),
+			facesLeft(m->faceCount()) {
+		const uint32_t faceCount = m->faceCount();
+		faceChartArray.resize(faceCount, -1);
+		faceCandidateArray.resize(faceCount, (uint32_t)-1);
+		// @@ Floyd for the whole mesh is too slow. We could compute floyd progressively per patch as the patch grows. We need a better solution to compute most central faces.
+		//computeShortestPaths();
+		// Precompute edge lengths and face areas.
+		uint32_t edgeCount = m->edgeCount();
+		edgeLengths.resize(edgeCount);
+		for (uint32_t i = 0; i < edgeCount; i++) {
+			uint32_t id = m->edgeAt(i)->id;
+			xaDebugAssert(id / 2 == i);
+#ifdef NDEBUG
+			id = 0; // silence unused parameter warning
+#endif
+			edgeLengths[i] = m->edgeAt(i)->length();
+		}
+		faceAreas.resize(faceCount);
+		for (uint32_t i = 0; i < faceCount; i++) {
+			faceAreas[i] = m->faceAt(i)->area();
+		}
+	}
+
+	~AtlasBuilder() {
+		const uint32_t chartCount = chartArray.size();
+		for (uint32_t i = 0; i < chartCount; i++) {
+			delete chartArray[i];
+		}
+	}
+
+	void markUnchartedFaces(const std::vector<uint32_t> &unchartedFaces) {
+		const uint32_t unchartedFaceCount = unchartedFaces.size();
+		for (uint32_t i = 0; i < unchartedFaceCount; i++) {
+			uint32_t f = unchartedFaces[i];
+			faceChartArray[f] = -2;
+			//faceCandidateArray[f] = -2; // @@ ?
+			removeCandidate(f);
+		}
+		xaDebugAssert(facesLeft >= unchartedFaceCount);
+		facesLeft -= unchartedFaceCount;
+	}
+
+	void computeShortestPaths() {
+		const uint32_t faceCount = mesh->faceCount();
+		shortestPaths.resize(faceCount * faceCount, FLT_MAX);
+		// Fill edges:
+		for (uint32_t i = 0; i < faceCount; i++) {
+			shortestPaths[i * faceCount + i] = 0.0f;
+			const halfedge::Face *face_i = mesh->faceAt(i);
+			Vector3 centroid_i = face_i->centroid();
+			for (halfedge::Face::ConstEdgeIterator it(face_i->edges()); !it.isDone(); it.advance()) {
+				const halfedge::Edge *edge = it.current();
+				if (!edge->isBoundary()) {
+					const halfedge::Face *face_j = edge->pair->face;
+					uint32_t j = face_j->id;
+					Vector3 centroid_j = face_j->centroid();
+					shortestPaths[i * faceCount + j] = shortestPaths[j * faceCount + i] = length(centroid_i - centroid_j);
+				}
+			}
+		}
+		// Use Floyd-Warshall algorithm to compute all paths:
+		for (uint32_t k = 0; k < faceCount; k++) {
+			for (uint32_t i = 0; i < faceCount; i++) {
+				for (uint32_t j = 0; j < faceCount; j++) {
+					shortestPaths[i * faceCount + j] = std::min(shortestPaths[i * faceCount + j], shortestPaths[i * faceCount + k] + shortestPaths[k * faceCount + j]);
+				}
+			}
+		}
+	}
+
+	void placeSeeds(float threshold, uint32_t maxSeedCount) {
+		// Instead of using a predefiened number of seeds:
+		// - Add seeds one by one, growing chart until a certain treshold.
+		// - Undo charts and restart growing process.
+		// @@ How can we give preference to faces far from sharp features as in the LSCM paper?
+		//   - those points can be found using a simple flood filling algorithm.
+		//   - how do we weight the probabilities?
+		for (uint32_t i = 0; i < maxSeedCount; i++) {
+			if (facesLeft == 0) {
+				// No faces left, stop creating seeds.
+				break;
+			}
+			createRandomChart(threshold);
+		}
+	}
+
+	void createRandomChart(float threshold) {
+		ChartBuildData *chart = new ChartBuildData(chartArray.size());
+		chartArray.push_back(chart);
+		// Pick random face that is not used by any chart yet.
+		uint32_t randomFaceIdx = rand.getRange(facesLeft - 1);
+		uint32_t i = 0;
+		for (uint32_t f = 0; f != randomFaceIdx; f++, i++) {
+			while (faceChartArray[i] != -1)
+				i++;
+		}
+		while (faceChartArray[i] != -1)
+			i++;
+		chart->seeds.push_back(i);
+		addFaceToChart(chart, i, true);
+		// Grow the chart as much as possible within the given threshold.
+		growChart(chart, threshold * 0.5f, facesLeft);
+		//growCharts(threshold - threshold * 0.75f / chartCount(), facesLeft);
+	}
+
+	void addFaceToChart(ChartBuildData *chart, uint32_t f, bool recomputeProxy = false) {
+		// Add face to chart.
+		chart->faces.push_back(f);
+		xaDebugAssert(faceChartArray[f] == -1);
+		faceChartArray[f] = chart->id;
+		facesLeft--;
+		// Update area and boundary length.
+		chart->area = evaluateChartArea(chart, f);
+		chart->boundaryLength = evaluateBoundaryLength(chart, f);
+		chart->normalSum = evaluateChartNormalSum(chart, f);
+		chart->centroidSum = evaluateChartCentroidSum(chart, f);
+		if (recomputeProxy) {
+			// Update proxy and candidate's priorities.
+			updateProxy(chart);
+		}
+		// Update candidates.
+		removeCandidate(f);
+		updateCandidates(chart, f);
+		updatePriorities(chart);
+	}
+
+	// Returns true if any of the charts can grow more.
+	bool growCharts(float threshold, uint32_t faceCount) {
+		// Using one global list.
+		faceCount = std::min(faceCount, facesLeft);
+		for (uint32_t i = 0; i < faceCount; i++) {
+			const Candidate &candidate = getBestCandidate();
+			if (candidate.metric > threshold) {
+				return false; // Can't grow more.
+			}
+			addFaceToChart(candidate.chart, candidate.face);
+		}
+		return facesLeft != 0; // Can continue growing.
+	}
+
+	bool growChart(ChartBuildData *chart, float threshold, uint32_t faceCount) {
+		// Try to add faceCount faces within threshold to chart.
+		for (uint32_t i = 0; i < faceCount;) {
+			if (chart->candidates.count() == 0 || chart->candidates.firstPriority() > threshold) {
+				return false;
+			}
+			uint32_t f = chart->candidates.pop();
+			if (faceChartArray[f] == -1) {
+				addFaceToChart(chart, f);
+				i++;
+			}
+		}
+		if (chart->candidates.count() == 0 || chart->candidates.firstPriority() > threshold) {
+			return false;
+		}
+		return true;
+	}
+
+	void resetCharts() {
+		const uint32_t faceCount = mesh->faceCount();
+		for (uint32_t i = 0; i < faceCount; i++) {
+			faceChartArray[i] = -1;
+			faceCandidateArray[i] = (uint32_t)-1;
+		}
+		facesLeft = faceCount;
+		candidateArray.clear();
+		const uint32_t chartCount = chartArray.size();
+		for (uint32_t i = 0; i < chartCount; i++) {
+			ChartBuildData *chart = chartArray[i];
+			const uint32_t seed = chart->seeds.back();
+			chart->area = 0.0f;
+			chart->boundaryLength = 0.0f;
+			chart->normalSum = Vector3(0);
+			chart->centroidSum = Vector3(0);
+			chart->faces.clear();
+			chart->candidates.clear();
+			addFaceToChart(chart, seed);
+		}
+	}
+
+	void updateCandidates(ChartBuildData *chart, uint32_t f) {
+		const halfedge::Face *face = mesh->faceAt(f);
+		// Traverse neighboring faces, add the ones that do not belong to any chart yet.
+		for (halfedge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) {
+			const halfedge::Edge *edge = it.current()->pair;
+			if (!edge->isBoundary()) {
+				uint32_t faceId = edge->face->id;
+				if (faceChartArray[faceId] == -1) {
+					chart->candidates.push(faceId);
+				}
+			}
+		}
+	}
+
+	void updateProxies() {
+		const uint32_t chartCount = chartArray.size();
+		for (uint32_t i = 0; i < chartCount; i++) {
+			updateProxy(chartArray[i]);
+		}
+	}
+
+	void updateProxy(ChartBuildData *chart) {
+		//#pragma message(NV_FILE_LINE "TODO: Use best fit plane instead of average normal.")
+		chart->planeNormal = normalizeSafe(chart->normalSum, Vector3(0), 0.0f);
+		chart->centroid = chart->centroidSum / float(chart->faces.size());
+	}
+
+	bool relocateSeeds() {
+		bool anySeedChanged = false;
+		const uint32_t chartCount = chartArray.size();
+		for (uint32_t i = 0; i < chartCount; i++) {
+			if (relocateSeed(chartArray[i])) {
+				anySeedChanged = true;
+			}
+		}
+		return anySeedChanged;
+	}
+
+	bool relocateSeed(ChartBuildData *chart) {
+		Vector3 centroid = computeChartCentroid(chart);
+		const uint32_t N = 10; // @@ Hardcoded to 10?
+		PriorityQueue bestTriangles(N);
+		// Find the first N triangles that fit the proxy best.
+		const uint32_t faceCount = chart->faces.size();
+		for (uint32_t i = 0; i < faceCount; i++) {
+			float priority = evaluateProxyFitMetric(chart, chart->faces[i]);
+			bestTriangles.push(priority, chart->faces[i]);
+		}
+		// Of those, choose the most central triangle.
+		uint32_t mostCentral;
+		float maxDistance = -1;
+		const uint32_t bestCount = bestTriangles.count();
+		for (uint32_t i = 0; i < bestCount; i++) {
+			const halfedge::Face *face = mesh->faceAt(bestTriangles.pairs[i].face);
+			Vector3 faceCentroid = face->triangleCenter();
+			float distance = length(centroid - faceCentroid);
+			if (distance > maxDistance) {
+				maxDistance = distance;
+				mostCentral = bestTriangles.pairs[i].face;
+			}
+		}
+		xaDebugAssert(maxDistance >= 0);
+		// In order to prevent k-means cyles we record all the previously chosen seeds.
+		uint32_t index = std::find(chart->seeds.begin(), chart->seeds.end(), mostCentral) - chart->seeds.begin();
+		if (index < chart->seeds.size()) {
+			// Move new seed to the end of the seed array.
+			uint32_t last = chart->seeds.size() - 1;
+			std::swap(chart->seeds[index], chart->seeds[last]);
+			return false;
+		} else {
+			// Append new seed.
+			chart->seeds.push_back(mostCentral);
+			return true;
+		}
+	}
+
+	void updatePriorities(ChartBuildData *chart) {
+		// Re-evaluate candidate priorities.
+		uint32_t candidateCount = chart->candidates.count();
+		for (uint32_t i = 0; i < candidateCount; i++) {
+			chart->candidates.pairs[i].priority = evaluatePriority(chart, chart->candidates.pairs[i].face);
+			if (faceChartArray[chart->candidates.pairs[i].face] == -1) {
+				updateCandidate(chart, chart->candidates.pairs[i].face, chart->candidates.pairs[i].priority);
+			}
+		}
+		// Sort candidates.
+		chart->candidates.sort();
+	}
+
+	// Evaluate combined metric.
+	float evaluatePriority(ChartBuildData *chart, uint32_t face) {
+		// Estimate boundary length and area:
+		float newBoundaryLength = evaluateBoundaryLength(chart, face);
+		float newChartArea = evaluateChartArea(chart, face);
+		float F = evaluateProxyFitMetric(chart, face);
+		float C = evaluateRoundnessMetric(chart, face, newBoundaryLength, newChartArea);
+		float P = evaluateStraightnessMetric(chart, face);
+		// Penalize faces that cross seams, reward faces that close seams or reach boundaries.
+		float N = evaluateNormalSeamMetric(chart, face);
+		float T = evaluateTextureSeamMetric(chart, face);
+		//float R = evaluateCompletenessMetric(chart, face);
+		//float D = evaluateDihedralAngleMetric(chart, face);
+		// @@ Add a metric based on local dihedral angle.
+		// @@ Tweaking the normal and texture seam metrics.
+		// - Cause more impedance. Never cross 90 degree edges.
+		// -
+		float cost = float(
+				options.proxyFitMetricWeight * F +
+				options.roundnessMetricWeight * C +
+				options.straightnessMetricWeight * P +
+				options.normalSeamMetricWeight * N +
+				options.textureSeamMetricWeight * T);
+		// Enforce limits strictly:
+		if (newChartArea > options.maxChartArea) cost = FLT_MAX;
+		if (newBoundaryLength > options.maxBoundaryLength) cost = FLT_MAX;
+		// Make sure normal seams are fully respected:
+		if (options.normalSeamMetricWeight >= 1000 && N != 0) cost = FLT_MAX;
+		xaAssert(std::isfinite(cost));
+		return cost;
+	}
+
+	// Returns a value in [0-1].
+	float evaluateProxyFitMetric(ChartBuildData *chart, uint32_t f) {
+		const halfedge::Face *face = mesh->faceAt(f);
+		Vector3 faceNormal = face->triangleNormal();
+		// Use plane fitting metric for now:
+		return 1 - dot(faceNormal, chart->planeNormal); // @@ normal deviations should be weighted by face area
+	}
+
+	float evaluateRoundnessMetric(ChartBuildData *chart, uint32_t /*face*/, float newBoundaryLength, float newChartArea) {
+		float roundness = square(chart->boundaryLength) / chart->area;
+		float newRoundness = square(newBoundaryLength) / newChartArea;
+		if (newRoundness > roundness) {
+			return square(newBoundaryLength) / (newChartArea * 4 * PI);
+		} else {
+			// Offer no impedance to faces that improve roundness.
+			return 0;
+		}
+	}
+
+	float evaluateStraightnessMetric(ChartBuildData *chart, uint32_t f) {
+		float l_out = 0.0f;
+		float l_in = 0.0f;
+		const halfedge::Face *face = mesh->faceAt(f);
+		for (halfedge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) {
+			const halfedge::Edge *edge = it.current();
+			float l = edgeLengths[edge->id / 2];
+			if (edge->isBoundary()) {
+				l_out += l;
+			} else {
+				uint32_t neighborFaceId = edge->pair->face->id;
+				if (faceChartArray[neighborFaceId] != chart->id) {
+					l_out += l;
+				} else {
+					l_in += l;
+				}
+			}
+		}
+		xaDebugAssert(l_in != 0.0f); // Candidate face must be adjacent to chart. @@ This is not true if the input mesh has zero-length edges.
+		float ratio = (l_out - l_in) / (l_out + l_in);
+		return std::min(ratio, 0.0f); // Only use the straightness metric to close gaps.
+	}
+
+	float evaluateNormalSeamMetric(ChartBuildData *chart, uint32_t f) {
+		float seamFactor = 0.0f;
+		float totalLength = 0.0f;
+		const halfedge::Face *face = mesh->faceAt(f);
+		for (halfedge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) {
+			const halfedge::Edge *edge = it.current();
+			if (edge->isBoundary()) {
+				continue;
+			}
+			const uint32_t neighborFaceId = edge->pair->face->id;
+			if (faceChartArray[neighborFaceId] != chart->id) {
+				continue;
+			}
+			//float l = edge->length();
+			float l = edgeLengths[edge->id / 2];
+			totalLength += l;
+			if (!edge->isSeam()) {
+				continue;
+			}
+			// Make sure it's a normal seam.
+			if (edge->isNormalSeam()) {
+				float d0 = clamp(dot(edge->vertex->nor, edge->pair->next->vertex->nor), 0.0f, 1.0f);
+				float d1 = clamp(dot(edge->next->vertex->nor, edge->pair->vertex->nor), 0.0f, 1.0f);
+				l *= 1 - (d0 + d1) * 0.5f;
+				seamFactor += l;
+			}
+		}
+		if (seamFactor == 0) return 0.0f;
+		return seamFactor / totalLength;
+	}
+
+	float evaluateTextureSeamMetric(ChartBuildData *chart, uint32_t f) {
+		float seamLength = 0.0f;
+		float totalLength = 0.0f;
+		const halfedge::Face *face = mesh->faceAt(f);
+		for (halfedge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) {
+			const halfedge::Edge *edge = it.current();
+			if (edge->isBoundary()) {
+				continue;
+			}
+			const uint32_t neighborFaceId = edge->pair->face->id;
+			if (faceChartArray[neighborFaceId] != chart->id) {
+				continue;
+			}
+			//float l = edge->length();
+			float l = edgeLengths[edge->id / 2];
+			totalLength += l;
+			if (!edge->isSeam()) {
+				continue;
+			}
+			// Make sure it's a texture seam.
+			if (edge->isTextureSeam()) {
+				seamLength += l;
+			}
+		}
+		if (seamLength == 0.0f) {
+			return 0.0f; // Avoid division by zero.
+		}
+		return seamLength / totalLength;
+	}
+
+	float evaluateChartArea(ChartBuildData *chart, uint32_t f) {
+		const halfedge::Face *face = mesh->faceAt(f);
+		return chart->area + faceAreas[face->id];
+	}
+
+	float evaluateBoundaryLength(ChartBuildData *chart, uint32_t f) {
+		float boundaryLength = chart->boundaryLength;
+		// Add new edges, subtract edges shared with the chart.
+		const halfedge::Face *face = mesh->faceAt(f);
+		for (halfedge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) {
+			const halfedge::Edge *edge = it.current();
+			//float edgeLength = edge->length();
+			float edgeLength = edgeLengths[edge->id / 2];
+			if (edge->isBoundary()) {
+				boundaryLength += edgeLength;
+			} else {
+				uint32_t neighborFaceId = edge->pair->face->id;
+				if (faceChartArray[neighborFaceId] != chart->id) {
+					boundaryLength += edgeLength;
+				} else {
+					boundaryLength -= edgeLength;
+				}
+			}
+		}
+		return std::max(0.0f, boundaryLength); // @@ Hack!
+	}
+
+	Vector3 evaluateChartNormalSum(ChartBuildData *chart, uint32_t f) {
+		const halfedge::Face *face = mesh->faceAt(f);
+		return chart->normalSum + face->triangleNormalAreaScaled();
+	}
+
+	Vector3 evaluateChartCentroidSum(ChartBuildData *chart, uint32_t f) {
+		const halfedge::Face *face = mesh->faceAt(f);
+		return chart->centroidSum + face->centroid();
+	}
+
+	Vector3 computeChartCentroid(const ChartBuildData *chart) {
+		Vector3 centroid(0);
+		const uint32_t faceCount = chart->faces.size();
+		for (uint32_t i = 0; i < faceCount; i++) {
+			const halfedge::Face *face = mesh->faceAt(chart->faces[i]);
+			centroid += face->triangleCenter();
+		}
+		return centroid / float(faceCount);
+	}
+
+	void fillHoles(float threshold) {
+		while (facesLeft > 0)
+			createRandomChart(threshold);
+	}
+
+	void mergeCharts() {
+		std::vector<float> sharedBoundaryLengths;
+		const uint32_t chartCount = chartArray.size();
+		for (int c = chartCount - 1; c >= 0; c--) {
+			sharedBoundaryLengths.clear();
+			sharedBoundaryLengths.resize(chartCount, 0.0f);
+			ChartBuildData *chart = chartArray[c];
+			float externalBoundary = 0.0f;
+			const uint32_t faceCount = chart->faces.size();
+			for (uint32_t i = 0; i < faceCount; i++) {
+				uint32_t f = chart->faces[i];
+				const halfedge::Face *face = mesh->faceAt(f);
+				for (halfedge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) {
+					const halfedge::Edge *edge = it.current();
+					//float l = edge->length();
+					float l = edgeLengths[edge->id / 2];
+					if (edge->isBoundary()) {
+						externalBoundary += l;
+					} else {
+						uint32_t neighborFace = edge->pair->face->id;
+						uint32_t neighborChart = faceChartArray[neighborFace];
+						if (neighborChart != (uint32_t)c) {
+							if ((edge->isSeam() && (edge->isNormalSeam() || edge->isTextureSeam())) || neighborChart == -2) {
+								externalBoundary += l;
+							} else {
+								sharedBoundaryLengths[neighborChart] += l;
+							}
+						}
+					}
+				}
+			}
+			for (int cc = chartCount - 1; cc >= 0; cc--) {
+				if (cc == c)
+					continue;
+				ChartBuildData *chart2 = chartArray[cc];
+				if (chart2 == NULL)
+					continue;
+				if (sharedBoundaryLengths[cc] > 0.8 * std::max(0.0f, chart->boundaryLength - externalBoundary)) {
+					// Try to avoid degenerate configurations.
+					if (chart2->boundaryLength > sharedBoundaryLengths[cc]) {
+						if (dot(chart2->planeNormal, chart->planeNormal) > -0.25) {
+							mergeChart(chart2, chart, sharedBoundaryLengths[cc]);
+							delete chart;
+							chartArray[c] = NULL;
+							break;
+						}
+					}
+				}
+				if (sharedBoundaryLengths[cc] > 0.20 * std::max(0.0f, chart->boundaryLength - externalBoundary)) {
+					// Compare proxies.
+					if (dot(chart2->planeNormal, chart->planeNormal) > 0) {
+						mergeChart(chart2, chart, sharedBoundaryLengths[cc]);
+						delete chart;
+						chartArray[c] = NULL;
+						break;
+					}
+				}
+			}
+		}
+		// Remove deleted charts.
+		for (int c = 0; c < int32_t(chartArray.size()); /*do not increment if removed*/) {
+			if (chartArray[c] == NULL) {
+				chartArray.erase(chartArray.begin() + c);
+				// Update faceChartArray.
+				const uint32_t faceCount = faceChartArray.size();
+				for (uint32_t i = 0; i < faceCount; i++) {
+					xaDebugAssert(faceChartArray[i] != -1);
+					xaDebugAssert(faceChartArray[i] != c);
+					xaDebugAssert(faceChartArray[i] <= int32_t(chartArray.size()));
+					if (faceChartArray[i] > c) {
+						faceChartArray[i]--;
+					}
+				}
+			} else {
+				chartArray[c]->id = c;
+				c++;
+			}
+		}
+	}
+
+	// @@ Cleanup.
+	struct Candidate {
+		uint32_t face;
+		ChartBuildData *chart;
+		float metric;
+	};
+
+	// @@ Get N best candidates in one pass.
+	const Candidate &getBestCandidate() const {
+		uint32_t best = 0;
+		float bestCandidateMetric = FLT_MAX;
+		const uint32_t candidateCount = candidateArray.size();
+		xaAssert(candidateCount > 0);
+		for (uint32_t i = 0; i < candidateCount; i++) {
+			const Candidate &candidate = candidateArray[i];
+			if (candidate.metric < bestCandidateMetric) {
+				bestCandidateMetric = candidate.metric;
+				best = i;
+			}
+		}
+		return candidateArray[best];
+	}
+
+	void removeCandidate(uint32_t f) {
+		int c = faceCandidateArray[f];
+		if (c != -1) {
+			faceCandidateArray[f] = (uint32_t)-1;
+			if (c == int(candidateArray.size() - 1)) {
+				candidateArray.pop_back();
+			} else {
+				// Replace with last.
+				candidateArray[c] = candidateArray[candidateArray.size() - 1];
+				candidateArray.pop_back();
+				faceCandidateArray[candidateArray[c].face] = c;
+			}
+		}
+	}
+
+	void updateCandidate(ChartBuildData *chart, uint32_t f, float metric) {
+		if (faceCandidateArray[f] == -1) {
+			const uint32_t index = candidateArray.size();
+			faceCandidateArray[f] = index;
+			candidateArray.resize(index + 1);
+			candidateArray[index].face = f;
+			candidateArray[index].chart = chart;
+			candidateArray[index].metric = metric;
+		} else {
+			int c = faceCandidateArray[f];
+			xaDebugAssert(c != -1);
+			Candidate &candidate = candidateArray[c];
+			xaDebugAssert(candidate.face == f);
+			if (metric < candidate.metric || chart == candidate.chart) {
+				candidate.metric = metric;
+				candidate.chart = chart;
+			}
+		}
+	}
+
+	void mergeChart(ChartBuildData *owner, ChartBuildData *chart, float sharedBoundaryLength) {
+		const uint32_t faceCount = chart->faces.size();
+		for (uint32_t i = 0; i < faceCount; i++) {
+			uint32_t f = chart->faces[i];
+			xaDebugAssert(faceChartArray[f] == chart->id);
+			faceChartArray[f] = owner->id;
+			owner->faces.push_back(f);
+		}
+		// Update adjacencies?
+		owner->area += chart->area;
+		owner->boundaryLength += chart->boundaryLength - sharedBoundaryLength;
+		owner->normalSum += chart->normalSum;
+		owner->centroidSum += chart->centroidSum;
+		updateProxy(owner);
+	}
+
+	uint32_t chartCount() const { return chartArray.size(); }
+	const std::vector<uint32_t> &chartFaces(uint32_t i) const { return chartArray[i]->faces; }
+
+	const halfedge::Mesh *mesh;
+	uint32_t facesLeft;
+	std::vector<int> faceChartArray;
+	std::vector<ChartBuildData *> chartArray;
+	std::vector<float> shortestPaths;
+	std::vector<float> edgeLengths;
+	std::vector<float> faceAreas;
+	std::vector<Candidate> candidateArray; //
+	std::vector<uint32_t> faceCandidateArray; // Map face index to candidate index.
+	MTRand rand;
+	CharterOptions options;
+};
+
+/// A chart is a connected set of faces with a certain topology (usually a disk).
+class Chart {
+public:
+	Chart() :
+			m_isDisk(false),
+			m_isVertexMapped(false) {}
+
+	void build(const halfedge::Mesh *originalMesh, const std::vector<uint32_t> &faceArray) {
+		// Copy face indices.
+		m_faceArray = faceArray;
+		const uint32_t meshVertexCount = originalMesh->vertexCount();
+		m_chartMesh.reset(new halfedge::Mesh());
+		m_unifiedMesh.reset(new halfedge::Mesh());
+		std::vector<uint32_t> chartMeshIndices(meshVertexCount, (uint32_t)~0);
+		std::vector<uint32_t> unifiedMeshIndices(meshVertexCount, (uint32_t)~0);
+		// Add vertices.
+		const uint32_t faceCount = faceArray.size();
+		for (uint32_t f = 0; f < faceCount; f++) {
+			const halfedge::Face *face = originalMesh->faceAt(faceArray[f]);
+			xaDebugAssert(face != NULL);
+			for (halfedge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) {
+				const halfedge::Vertex *vertex = it.current()->vertex;
+				const halfedge::Vertex *unifiedVertex = vertex->firstColocal();
+				if (unifiedMeshIndices[unifiedVertex->id] == ~0) {
+					unifiedMeshIndices[unifiedVertex->id] = m_unifiedMesh->vertexCount();
+					xaDebugAssert(vertex->pos == unifiedVertex->pos);
+					m_unifiedMesh->addVertex(vertex->pos);
+				}
+				if (chartMeshIndices[vertex->id] == ~0) {
+					chartMeshIndices[vertex->id] = m_chartMesh->vertexCount();
+					m_chartToOriginalMap.push_back(vertex->original_id);
+					m_chartToUnifiedMap.push_back(unifiedMeshIndices[unifiedVertex->id]);
+					halfedge::Vertex *v = m_chartMesh->addVertex(vertex->pos);
+					v->nor = vertex->nor;
+					v->tex = vertex->tex;
+				}
+			}
+		}
+		// This is ignoring the canonical map:
+		// - Is it really necessary to link colocals?
+		m_chartMesh->linkColocals();
+		//m_unifiedMesh->linkColocals();  // Not strictly necessary, no colocals in the unified mesh. # Wrong.
+		// This check is not valid anymore, if the original mesh vertices were linked with a canonical map, then it might have
+		// some colocal vertices that were unlinked. So, the unified mesh might have some duplicate vertices, because firstColocal()
+		// is not guaranteed to return the same vertex for two colocal vertices.
+		//xaAssert(m_chartMesh->colocalVertexCount() == m_unifiedMesh->vertexCount());
+		// Is that OK? What happens in meshes were that happens? Does anything break? Apparently not...
+		std::vector<uint32_t> faceIndices;
+		faceIndices.reserve(7);
+		// Add faces.
+		for (uint32_t f = 0; f < faceCount; f++) {
+			const halfedge::Face *face = originalMesh->faceAt(faceArray[f]);
+			xaDebugAssert(face != NULL);
+			faceIndices.clear();
+			for (halfedge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) {
+				const halfedge::Vertex *vertex = it.current()->vertex;
+				xaDebugAssert(vertex != NULL);
+				faceIndices.push_back(chartMeshIndices[vertex->id]);
+			}
+			m_chartMesh->addFace(faceIndices);
+			faceIndices.clear();
+			for (halfedge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) {
+				const halfedge::Vertex *vertex = it.current()->vertex;
+				xaDebugAssert(vertex != NULL);
+				vertex = vertex->firstColocal();
+				faceIndices.push_back(unifiedMeshIndices[vertex->id]);
+			}
+			m_unifiedMesh->addFace(faceIndices);
+		}
+		m_chartMesh->linkBoundary();
+		m_unifiedMesh->linkBoundary();
+		//exportMesh(m_unifiedMesh.ptr(), "debug_input.obj");
+		if (m_unifiedMesh->splitBoundaryEdges()) {
+			m_unifiedMesh.reset(halfedge::unifyVertices(m_unifiedMesh.get()));
+		}
+		//exportMesh(m_unifiedMesh.ptr(), "debug_split.obj");
+		// Closing the holes is not always the best solution and does not fix all the problems.
+		// We need to do some analysis of the holes and the genus to:
+		// - Find cuts that reduce genus.
+		// - Find cuts to connect holes.
+		// - Use minimal spanning trees or seamster.
+		if (!closeHoles()) {
+			/*static int pieceCount = 0;
+			StringBuilder fileName;
+			fileName.format("debug_hole_%d.obj", pieceCount++);
+			exportMesh(m_unifiedMesh.ptr(), fileName.str());*/
+		}
+		m_unifiedMesh.reset(halfedge::triangulate(m_unifiedMesh.get()));
+		//exportMesh(m_unifiedMesh.ptr(), "debug_triangulated.obj");
+		// Analyze chart topology.
+		halfedge::MeshTopology topology(m_unifiedMesh.get());
+		m_isDisk = topology.isDisk();
+	}
+
+	void buildVertexMap(const halfedge::Mesh *originalMesh, const std::vector<uint32_t> &unchartedMaterialArray) {
+		xaAssert(m_chartMesh.get() == NULL && m_unifiedMesh.get() == NULL);
+		m_isVertexMapped = true;
+		// Build face indices.
+		m_faceArray.clear();
+		const uint32_t meshFaceCount = originalMesh->faceCount();
+		for (uint32_t f = 0; f < meshFaceCount; f++) {
+			const halfedge::Face *face = originalMesh->faceAt(f);
+			if (std::find(unchartedMaterialArray.begin(), unchartedMaterialArray.end(), face->material) != unchartedMaterialArray.end()) {
+				m_faceArray.push_back(f);
+			}
+		}
+		const uint32_t faceCount = m_faceArray.size();
+		if (faceCount == 0) {
+			return;
+		}
+		// @@ The chartMesh construction is basically the same as with regular charts, don't duplicate!
+		const uint32_t meshVertexCount = originalMesh->vertexCount();
+		m_chartMesh.reset(new halfedge::Mesh());
+		std::vector<uint32_t> chartMeshIndices(meshVertexCount, (uint32_t)~0);
+		// Vertex map mesh only has disconnected vertices.
+		for (uint32_t f = 0; f < faceCount; f++) {
+			const halfedge::Face *face = originalMesh->faceAt(m_faceArray[f]);
+			xaDebugAssert(face != NULL);
+			for (halfedge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) {
+				const halfedge::Vertex *vertex = it.current()->vertex;
+				if (chartMeshIndices[vertex->id] == ~0) {
+					chartMeshIndices[vertex->id] = m_chartMesh->vertexCount();
+					m_chartToOriginalMap.push_back(vertex->original_id);
+					halfedge::Vertex *v = m_chartMesh->addVertex(vertex->pos);
+					v->nor = vertex->nor;
+					v->tex = vertex->tex; // @@ Not necessary.
+				}
+			}
+		}
+		// @@ Link colocals using the original mesh canonical map? Build canonical map on the fly? Do we need to link colocals at all for this?
+		//m_chartMesh->linkColocals();
+		std::vector<uint32_t> faceIndices;
+		faceIndices.reserve(7);
+		// Add faces.
+		for (uint32_t f = 0; f < faceCount; f++) {
+			const halfedge::Face *face = originalMesh->faceAt(m_faceArray[f]);
+			xaDebugAssert(face != NULL);
+			faceIndices.clear();
+			for (halfedge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) {
+				const halfedge::Vertex *vertex = it.current()->vertex;
+				xaDebugAssert(vertex != NULL);
+				xaDebugAssert(chartMeshIndices[vertex->id] != ~0);
+				faceIndices.push_back(chartMeshIndices[vertex->id]);
+			}
+			halfedge::Face *new_face = m_chartMesh->addFace(faceIndices);
+			xaDebugAssert(new_face != NULL);
+#ifdef NDEBUG
+			new_face = NULL; // silence unused parameter warning
+#endif
+		}
+		m_chartMesh->linkBoundary();
+		const uint32_t chartVertexCount = m_chartMesh->vertexCount();
+		Box bounds;
+		bounds.clearBounds();
+		for (uint32_t i = 0; i < chartVertexCount; i++) {
+			halfedge::Vertex *vertex = m_chartMesh->vertexAt(i);
+			bounds.addPointToBounds(vertex->pos);
+		}
+		ProximityGrid grid;
+		grid.init(bounds, chartVertexCount);
+		for (uint32_t i = 0; i < chartVertexCount; i++) {
+			halfedge::Vertex *vertex = m_chartMesh->vertexAt(i);
+			grid.add(vertex->pos, i);
+		}
+		uint32_t texelCount = 0;
+		const float positionThreshold = 0.01f;
+		const float normalThreshold = 0.01f;
+		uint32_t verticesVisited = 0;
+		uint32_t cellsVisited = 0;
+		std::vector<int> vertexIndexArray(chartVertexCount, -1); // Init all indices to -1.
+		// Traverse vertices in morton order. @@ It may be more interesting to sort them based on orientation.
+		const uint32_t cellCodeCount = grid.mortonCount();
+		for (uint32_t cellCode = 0; cellCode < cellCodeCount; cellCode++) {
+			int cell = grid.mortonIndex(cellCode);
+			if (cell < 0) continue;
+			cellsVisited++;
+			const std::vector<uint32_t> &indexArray = grid.cellArray[cell].indexArray;
+			for (uint32_t i = 0; i < indexArray.size(); i++) {
+				uint32_t idx = indexArray[i];
+				halfedge::Vertex *vertex = m_chartMesh->vertexAt(idx);
+				xaDebugAssert(vertexIndexArray[idx] == -1);
+				std::vector<uint32_t> neighbors;
+				grid.gather(vertex->pos, positionThreshold, /*ref*/ neighbors);
+				// Compare against all nearby vertices, cluster greedily.
+				for (uint32_t j = 0; j < neighbors.size(); j++) {
+					uint32_t otherIdx = neighbors[j];
+					if (vertexIndexArray[otherIdx] != -1) {
+						halfedge::Vertex *otherVertex = m_chartMesh->vertexAt(otherIdx);
+						if (distance(vertex->pos, otherVertex->pos) < positionThreshold &&
+								distance(vertex->nor, otherVertex->nor) < normalThreshold) {
+							vertexIndexArray[idx] = vertexIndexArray[otherIdx];
+							break;
+						}
+					}
+				}
+				// If index not assigned, assign new one.
+				if (vertexIndexArray[idx] == -1) {
+					vertexIndexArray[idx] = texelCount++;
+				}
+				verticesVisited++;
+			}
+		}
+		xaDebugAssert(cellsVisited == grid.cellArray.size());
+		xaDebugAssert(verticesVisited == chartVertexCount);
+		vertexMapWidth = ftoi_ceil(sqrtf(float(texelCount)));
+		vertexMapWidth = (vertexMapWidth + 3) & ~3; // Width aligned to 4.
+		vertexMapHeight = vertexMapWidth == 0 ? 0 : (texelCount + vertexMapWidth - 1) / vertexMapWidth;
+		//vertexMapHeight = (vertexMapHeight + 3) & ~3;                           // Height aligned to 4.
+		xaDebugAssert(vertexMapWidth >= vertexMapHeight);
+		xaPrint("Reduced vertex count from %d to %d.\n", chartVertexCount, texelCount);
+		// Lay down the clustered vertices in morton order.
+		std::vector<uint32_t> texelCodes(texelCount);
+		// For each texel, assign one morton code.
+		uint32_t texelCode = 0;
+		for (uint32_t i = 0; i < texelCount; i++) {
+			uint32_t x, y;
+			do {
+				x = morton::decodeMorton2X(texelCode);
+				y = morton::decodeMorton2Y(texelCode);
+				texelCode++;
+			} while (x >= uint32_t(vertexMapWidth) || y >= uint32_t(vertexMapHeight));
+			texelCodes[i] = texelCode - 1;
+		}
+		for (uint32_t i = 0; i < chartVertexCount; i++) {
+			halfedge::Vertex *vertex = m_chartMesh->vertexAt(i);
+			int idx = vertexIndexArray[i];
+			if (idx != -1) {
+				uint32_t tc = texelCodes[idx];
+				uint32_t x = morton::decodeMorton2X(tc);
+				uint32_t y = morton::decodeMorton2Y(tc);
+				vertex->tex.x = float(x);
+				vertex->tex.y = float(y);
+			}
+		}
+	}
+
+	bool closeHoles() {
+		xaDebugAssert(!m_isVertexMapped);
+		std::vector<halfedge::Edge *> boundaryEdges;
+		getBoundaryEdges(m_unifiedMesh.get(), boundaryEdges);
+		uint32_t boundaryCount = boundaryEdges.size();
+		if (boundaryCount <= 1) {
+			// Nothing to close.
+			return true;
+		}
+		// Compute lengths and areas.
+		std::vector<float> boundaryLengths;
+		for (uint32_t i = 0; i < boundaryCount; i++) {
+			const halfedge::Edge *startEdge = boundaryEdges[i];
+			xaAssert(startEdge->face == NULL);
+			//float boundaryEdgeCount = 0;
+			float boundaryLength = 0.0f;
+			//Vector3 boundaryCentroid(zero);
+			const halfedge::Edge *edge = startEdge;
+			do {
+				Vector3 t0 = edge->from()->pos;
+				Vector3 t1 = edge->to()->pos;
+				//boundaryEdgeCount++;
+				boundaryLength += length(t1 - t0);
+				//boundaryCentroid += edge->vertex()->pos;
+				edge = edge->next;
+			} while (edge != startEdge);
+			boundaryLengths.push_back(boundaryLength);
+			//boundaryCentroids.append(boundaryCentroid / boundaryEdgeCount);
+		}
+		// Find disk boundary.
+		uint32_t diskBoundary = 0;
+		float maxLength = boundaryLengths[0];
+		for (uint32_t i = 1; i < boundaryCount; i++) {
+			if (boundaryLengths[i] > maxLength) {
+				maxLength = boundaryLengths[i];
+				diskBoundary = i;
+			}
+		}
+		// Close holes.
+		for (uint32_t i = 0; i < boundaryCount; i++) {
+			if (diskBoundary == i) {
+				// Skip disk boundary.
+				continue;
+			}
+			halfedge::Edge *startEdge = boundaryEdges[i];
+			xaDebugAssert(startEdge != NULL);
+			xaDebugAssert(startEdge->face == NULL);
+			std::vector<halfedge::Vertex *> vertexLoop;
+			std::vector<halfedge::Edge *> edgeLoop;
+			halfedge::Edge *edge = startEdge;
+			do {
+				halfedge::Vertex *vertex = edge->next->vertex; // edge->to()
+				uint32_t j;
+				for (j = 0; j < vertexLoop.size(); j++) {
+					if (vertex->isColocal(vertexLoop[j])) {
+						break;
+					}
+				}
+				bool isCrossing = (j != vertexLoop.size());
+				if (isCrossing) {
+					halfedge::Edge *prev = edgeLoop[j]; // Previous edge before the loop.
+					halfedge::Edge *next = edge->next; // Next edge after the loop.
+					xaDebugAssert(prev->to()->isColocal(next->from()));
+					// Close loop.
+					edgeLoop.push_back(edge);
+					closeLoop(j + 1, edgeLoop);
+					// Link boundary loop.
+					prev->setNext(next);
+					vertex->setEdge(next);
+					// Start over again.
+					vertexLoop.clear();
+					edgeLoop.clear();
+					edge = startEdge;
+					vertex = edge->to();
+				}
+				vertexLoop.push_back(vertex);
+				edgeLoop.push_back(edge);
+				edge = edge->next;
+			} while (edge != startEdge);
+			closeLoop(0, edgeLoop);
+		}
+		getBoundaryEdges(m_unifiedMesh.get(), boundaryEdges);
+		boundaryCount = boundaryEdges.size();
+		xaDebugAssert(boundaryCount == 1);
+		return boundaryCount == 1;
+	}
+
+	bool isDisk() const {
+		return m_isDisk;
+	}
+	bool isVertexMapped() const {
+		return m_isVertexMapped;
+	}
+
+	uint32_t vertexCount() const {
+		return m_chartMesh->vertexCount();
+	}
+	uint32_t colocalVertexCount() const {
+		return m_unifiedMesh->vertexCount();
+	}
+
+	uint32_t faceCount() const {
+		return m_faceArray.size();
+	}
+	uint32_t faceAt(uint32_t i) const {
+		return m_faceArray[i];
+	}
+
+	const halfedge::Mesh *chartMesh() const {
+		return m_chartMesh.get();
+	}
+	halfedge::Mesh *chartMesh() {
+		return m_chartMesh.get();
+	}
+	const halfedge::Mesh *unifiedMesh() const {
+		return m_unifiedMesh.get();
+	}
+	halfedge::Mesh *unifiedMesh() {
+		return m_unifiedMesh.get();
+	}
+
+	//uint32_t vertexIndex(uint32_t i) const { return m_vertexIndexArray[i]; }
+
+	uint32_t mapChartVertexToOriginalVertex(uint32_t i) const {
+		return m_chartToOriginalMap[i];
+	}
+	uint32_t mapChartVertexToUnifiedVertex(uint32_t i) const {
+		return m_chartToUnifiedMap[i];
+	}
+
+	const std::vector<uint32_t> &faceArray() const {
+		return m_faceArray;
+	}
+
+	// Transfer parameterization from unified mesh to chart mesh.
+	void transferParameterization() {
+		xaDebugAssert(!m_isVertexMapped);
+		uint32_t vertexCount = m_chartMesh->vertexCount();
+		for (uint32_t v = 0; v < vertexCount; v++) {
+			halfedge::Vertex *vertex = m_chartMesh->vertexAt(v);
+			halfedge::Vertex *unifiedVertex = m_unifiedMesh->vertexAt(mapChartVertexToUnifiedVertex(v));
+			vertex->tex = unifiedVertex->tex;
+		}
+	}
+
+	float computeSurfaceArea() const {
+		return halfedge::computeSurfaceArea(m_chartMesh.get()) * scale;
+	}
+
+	float computeParametricArea() const {
+		// This only makes sense in parameterized meshes.
+		xaDebugAssert(m_isDisk);
+		xaDebugAssert(!m_isVertexMapped);
+		return halfedge::computeParametricArea(m_chartMesh.get());
+	}
+
+	Vector2 computeParametricBounds() const {
+		// This only makes sense in parameterized meshes.
+		xaDebugAssert(m_isDisk);
+		xaDebugAssert(!m_isVertexMapped);
+		Box bounds;
+		bounds.clearBounds();
+		uint32_t vertexCount = m_chartMesh->vertexCount();
+		for (uint32_t v = 0; v < vertexCount; v++) {
+			halfedge::Vertex *vertex = m_chartMesh->vertexAt(v);
+			bounds.addPointToBounds(Vector3(vertex->tex, 0));
+		}
+		return bounds.extents().xy();
+	}
+
+	float scale = 1.0f;
+	uint32_t vertexMapWidth;
+	uint32_t vertexMapHeight;
+	bool blockAligned = true;
+
+private:
+	bool closeLoop(uint32_t start, const std::vector<halfedge::Edge *> &loop) {
+		const uint32_t vertexCount = loop.size() - start;
+		xaDebugAssert(vertexCount >= 3);
+		if (vertexCount < 3) return false;
+		xaDebugAssert(loop[start]->vertex->isColocal(loop[start + vertexCount - 1]->to()));
+		// If the hole is planar, then we add a single face that will be properly triangulated later.
+		// If the hole is not planar, we add a triangle fan with a vertex at the hole centroid.
+		// This is still a bit of a hack. There surely are better hole filling algorithms out there.
+		std::vector<Vector3> points(vertexCount);
+		for (uint32_t i = 0; i < vertexCount; i++) {
+			points[i] = loop[start + i]->vertex->pos;
+		}
+		bool isPlanar = Fit::isPlanar(vertexCount, points.data());
+		if (isPlanar) {
+			// Add face and connect edges.
+			halfedge::Face *face = m_unifiedMesh->addFace();
+			for (uint32_t i = 0; i < vertexCount; i++) {
+				halfedge::Edge *edge = loop[start + i];
+				edge->face = face;
+				edge->setNext(loop[start + (i + 1) % vertexCount]);
+			}
+			face->edge = loop[start];
+			xaDebugAssert(face->isValid());
+		} else {
+			// If the polygon is not planar, we just cross our fingers, and hope this will work:
+			// Compute boundary centroid:
+			Vector3 centroidPos(0);
+			for (uint32_t i = 0; i < vertexCount; i++) {
+				centroidPos += points[i];
+			}
+			centroidPos *= (1.0f / vertexCount);
+			halfedge::Vertex *centroid = m_unifiedMesh->addVertex(centroidPos);
+			// Add one pair of edges for each boundary vertex.
+			for (uint32_t j = vertexCount - 1, i = 0; i < vertexCount; j = i++) {
+				halfedge::Face *face = m_unifiedMesh->addFace(centroid->id, loop[start + j]->vertex->id, loop[start + i]->vertex->id);
+				xaDebugAssert(face != NULL);
+#ifdef NDEBUG
+				face = NULL; // silence unused parameter warning
+#endif
+			}
+		}
+		return true;
+	}
+
+	static void getBoundaryEdges(halfedge::Mesh *mesh, std::vector<halfedge::Edge *> &boundaryEdges) {
+		xaDebugAssert(mesh != NULL);
+		const uint32_t edgeCount = mesh->edgeCount();
+		BitArray bitFlags(edgeCount);
+		bitFlags.clearAll();
+		boundaryEdges.clear();
+		// Search for boundary edges. Mark all the edges that belong to the same boundary.
+		for (uint32_t e = 0; e < edgeCount; e++) {
+			halfedge::Edge *startEdge = mesh->edgeAt(e);
+			if (startEdge != NULL && startEdge->isBoundary() && bitFlags.bitAt(e) == false) {
+				xaDebugAssert(startEdge->face != NULL);
+				xaDebugAssert(startEdge->pair->face == NULL);
+				startEdge = startEdge->pair;
+				const halfedge::Edge *edge = startEdge;
+				do {
+					xaDebugAssert(edge->face == NULL);
+					xaDebugAssert(bitFlags.bitAt(edge->id / 2) == false);
+					bitFlags.setBitAt(edge->id / 2);
+					edge = edge->next;
+				} while (startEdge != edge);
+				boundaryEdges.push_back(startEdge);
+			}
+		}
+	}
+
+	// Chart mesh.
+	std::auto_ptr<halfedge::Mesh> m_chartMesh;
+
+	std::auto_ptr<halfedge::Mesh> m_unifiedMesh;
+	bool m_isDisk;
+	bool m_isVertexMapped;
+
+	// List of faces of the original mesh that belong to this chart.
+	std::vector<uint32_t> m_faceArray;
+
+	// Map vertices of the chart mesh to vertices of the original mesh.
+	std::vector<uint32_t> m_chartToOriginalMap;
+
+	std::vector<uint32_t> m_chartToUnifiedMap;
+};
+
+// Estimate quality of existing parameterization.
+class ParameterizationQuality {
+public:
+	ParameterizationQuality() {
+		m_totalTriangleCount = 0;
+		m_flippedTriangleCount = 0;
+		m_zeroAreaTriangleCount = 0;
+		m_parametricArea = 0.0f;
+		m_geometricArea = 0.0f;
+		m_stretchMetric = 0.0f;
+		m_maxStretchMetric = 0.0f;
+		m_conformalMetric = 0.0f;
+		m_authalicMetric = 0.0f;
+	}
+
+	ParameterizationQuality(const halfedge::Mesh *mesh) {
+		xaDebugAssert(mesh != NULL);
+		m_totalTriangleCount = 0;
+		m_flippedTriangleCount = 0;
+		m_zeroAreaTriangleCount = 0;
+		m_parametricArea = 0.0f;
+		m_geometricArea = 0.0f;
+		m_stretchMetric = 0.0f;
+		m_maxStretchMetric = 0.0f;
+		m_conformalMetric = 0.0f;
+		m_authalicMetric = 0.0f;
+		const uint32_t faceCount = mesh->faceCount();
+		for (uint32_t f = 0; f < faceCount; f++) {
+			const halfedge::Face *face = mesh->faceAt(f);
+			const halfedge::Vertex *vertex0 = NULL;
+			Vector3 p[3];
+			Vector2 t[3];
+			for (halfedge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) {
+				const halfedge::Edge *edge = it.current();
+				if (vertex0 == NULL) {
+					vertex0 = edge->vertex;
+					p[0] = vertex0->pos;
+					t[0] = vertex0->tex;
+				} else if (edge->to() != vertex0) {
+					p[1] = edge->from()->pos;
+					p[2] = edge->to()->pos;
+					t[1] = edge->from()->tex;
+					t[2] = edge->to()->tex;
+					processTriangle(p, t);
+				}
+			}
+		}
+		if (m_flippedTriangleCount + m_zeroAreaTriangleCount == faceCount) {
+			// If all triangles are flipped, then none is.
+			m_flippedTriangleCount = 0;
+		}
+		xaDebugAssert(std::isfinite(m_parametricArea) && m_parametricArea >= 0);
+		xaDebugAssert(std::isfinite(m_geometricArea) && m_geometricArea >= 0);
+		xaDebugAssert(std::isfinite(m_stretchMetric));
+		xaDebugAssert(std::isfinite(m_maxStretchMetric));
+		xaDebugAssert(std::isfinite(m_conformalMetric));
+		xaDebugAssert(std::isfinite(m_authalicMetric));
+	}
+
+	bool isValid() const {
+		return m_flippedTriangleCount == 0; // @@ Does not test for self-overlaps.
+	}
+
+	float rmsStretchMetric() const {
+		if (m_geometricArea == 0) return 0.0f;
+		float normFactor = sqrtf(m_parametricArea / m_geometricArea);
+		return sqrtf(m_stretchMetric / m_geometricArea) * normFactor;
+	}
+
+	float maxStretchMetric() const {
+		if (m_geometricArea == 0) return 0.0f;
+		float normFactor = sqrtf(m_parametricArea / m_geometricArea);
+		return m_maxStretchMetric * normFactor;
+	}
+
+	float rmsConformalMetric() const {
+		if (m_geometricArea == 0) return 0.0f;
+		return sqrtf(m_conformalMetric / m_geometricArea);
+	}
+
+	float maxAuthalicMetric() const {
+		if (m_geometricArea == 0) return 0.0f;
+		return sqrtf(m_authalicMetric / m_geometricArea);
+	}
+
+	void operator+=(const ParameterizationQuality &pq) {
+		m_totalTriangleCount += pq.m_totalTriangleCount;
+		m_flippedTriangleCount += pq.m_flippedTriangleCount;
+		m_zeroAreaTriangleCount += pq.m_zeroAreaTriangleCount;
+		m_parametricArea += pq.m_parametricArea;
+		m_geometricArea += pq.m_geometricArea;
+		m_stretchMetric += pq.m_stretchMetric;
+		m_maxStretchMetric = std::max(m_maxStretchMetric, pq.m_maxStretchMetric);
+		m_conformalMetric += pq.m_conformalMetric;
+		m_authalicMetric += pq.m_authalicMetric;
+	}
+
+private:
+	void processTriangle(Vector3 q[3], Vector2 p[3]) {
+		m_totalTriangleCount++;
+		// Evaluate texture stretch metric. See:
+		// - "Texture Mapping Progressive Meshes", Sander, Snyder, Gortler & Hoppe
+		// - "Mesh Parameterization: Theory and Practice", Siggraph'07 Course Notes, Hormann, Levy & Sheffer.
+		float t1 = p[0].x;
+		float s1 = p[0].y;
+		float t2 = p[1].x;
+		float s2 = p[1].y;
+		float t3 = p[2].x;
+		float s3 = p[2].y;
+		float geometricArea = length(cross(q[1] - q[0], q[2] - q[0])) / 2;
+		float parametricArea = ((s2 - s1) * (t3 - t1) - (s3 - s1) * (t2 - t1)) / 2;
+		if (isZero(parametricArea)) {
+			m_zeroAreaTriangleCount++;
+			return;
+		}
+		Vector3 Ss = (q[0] * (t2 - t3) + q[1] * (t3 - t1) + q[2] * (t1 - t2)) / (2 * parametricArea);
+		Vector3 St = (q[0] * (s3 - s2) + q[1] * (s1 - s3) + q[2] * (s2 - s1)) / (2 * parametricArea);
+		float a = dot(Ss, Ss); // E
+		float b = dot(Ss, St); // F
+		float c = dot(St, St); // G
+		// Compute eigen-values of the first fundamental form:
+		float sigma1 = sqrtf(0.5f * std::max(0.0f, a + c - sqrtf(square(a - c) + 4 * square(b)))); // gamma uppercase, min eigenvalue.
+		float sigma2 = sqrtf(0.5f * std::max(0.0f, a + c + sqrtf(square(a - c) + 4 * square(b)))); // gamma lowercase, max eigenvalue.
+		xaAssert(sigma2 >= sigma1);
+		// isometric: sigma1 = sigma2 = 1
+		// conformal: sigma1 / sigma2 = 1
+		// authalic: sigma1 * sigma2 = 1
+		float rmsStretch = sqrtf((a + c) * 0.5f);
+		float rmsStretch2 = sqrtf((square(sigma1) + square(sigma2)) * 0.5f);
+		xaDebugAssert(equal(rmsStretch, rmsStretch2, 0.01f));
+#ifdef NDEBUG
+		rmsStretch2 = 0; // silence unused parameter warning
+#endif
+		if (parametricArea < 0.0f) {
+			// Count flipped triangles.
+			m_flippedTriangleCount++;
+			parametricArea = fabsf(parametricArea);
+		}
+		m_stretchMetric += square(rmsStretch) * geometricArea;
+		m_maxStretchMetric = std::max(m_maxStretchMetric, sigma2);
+		if (!isZero(sigma1, 0.000001f)) {
+			// sigma1 is zero when geometricArea is zero.
+			m_conformalMetric += (sigma2 / sigma1) * geometricArea;
+		}
+		m_authalicMetric += (sigma1 * sigma2) * geometricArea;
+		// Accumulate total areas.
+		m_geometricArea += geometricArea;
+		m_parametricArea += parametricArea;
+		//triangleConformalEnergy(q, p);
+	}
+
+	uint32_t m_totalTriangleCount;
+	uint32_t m_flippedTriangleCount;
+	uint32_t m_zeroAreaTriangleCount;
+	float m_parametricArea;
+	float m_geometricArea;
+	float m_stretchMetric;
+	float m_maxStretchMetric;
+	float m_conformalMetric;
+	float m_authalicMetric;
+};
+
+// Set of charts corresponding to a single mesh.
+class MeshCharts {
+public:
+	MeshCharts(const halfedge::Mesh *mesh) :
+			m_mesh(mesh) {}
+
+	~MeshCharts() {
+		for (size_t i = 0; i < m_chartArray.size(); i++)
+			delete m_chartArray[i];
+	}
+
+	uint32_t chartCount() const {
+		return m_chartArray.size();
+	}
+	uint32_t vertexCount() const {
+		return m_totalVertexCount;
+	}
+
+	const Chart *chartAt(uint32_t i) const {
+		return m_chartArray[i];
+	}
+	Chart *chartAt(uint32_t i) {
+		return m_chartArray[i];
+	}
+
+	// Extract the charts of the input mesh.
+	void extractCharts() {
+		const uint32_t faceCount = m_mesh->faceCount();
+		int first = 0;
+		std::vector<uint32_t> queue;
+		queue.reserve(faceCount);
+		BitArray bitFlags(faceCount);
+		bitFlags.clearAll();
+		for (uint32_t f = 0; f < faceCount; f++) {
+			if (bitFlags.bitAt(f) == false) {
+				// Start new patch. Reset queue.
+				first = 0;
+				queue.clear();
+				queue.push_back(f);
+				bitFlags.setBitAt(f);
+				while (first != (int)queue.size()) {
+					const halfedge::Face *face = m_mesh->faceAt(queue[first]);
+					// Visit face neighbors of queue[first]
+					for (halfedge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) {
+						const halfedge::Edge *edge = it.current();
+						xaDebugAssert(edge->pair != NULL);
+						if (!edge->isBoundary() && /*!edge->isSeam()*/
+								//!(edge->from()->tex() != edge->pair()->to()->tex() || edge->to()->tex() != edge->pair()->from()->tex()))
+								!(edge->from() != edge->pair->to() || edge->to() != edge->pair->from())) { // Preserve existing seams (not just texture seams).
+							const halfedge::Face *neighborFace = edge->pair->face;
+							xaDebugAssert(neighborFace != NULL);
+							if (bitFlags.bitAt(neighborFace->id) == false) {
+								queue.push_back(neighborFace->id);
+								bitFlags.setBitAt(neighborFace->id);
+							}
+						}
+					}
+					first++;
+				}
+				Chart *chart = new Chart();
+				chart->build(m_mesh, queue);
+				m_chartArray.push_back(chart);
+			}
+		}
+	}
+
+	/*
+	Compute charts using a simple segmentation algorithm.
+
+	LSCM:
+	- identify sharp features using local dihedral angles.
+	- identify seed faces farthest from sharp features.
+	- grow charts from these seeds.
+
+	MCGIM:
+	- phase 1: chart growth
+	  - grow all charts simultaneously using dijkstra search on the dual graph of the mesh.
+	  - graph edges are weighted based on planarity metric.
+	  - metric uses distance to global chart normal.
+	  - terminate when all faces have been assigned.
+	- phase 2: seed computation:
+	  - place new seed of the chart at the most interior face.
+	  - most interior is evaluated using distance metric only.
+
+	- method repeates the two phases, until the location of the seeds does not change.
+	  - cycles are detected by recording all the previous seeds and chartification terminates.
+
+	D-Charts:
+
+	- Uniaxial conic metric:
+	  - N_c = axis of the generalized cone that best fits the chart. (cone can a be cylinder or a plane).
+	  - omega_c = angle between the face normals and the axis.
+	  - Fitting error between chart C and tringle t: F(c,t) = (N_c*n_t - cos(omega_c))^2
+
+	- Compactness metrics:
+	  - Roundness:
+		- C(c,t) = pi * D(S_c,t)^2 / A_c
+		- S_c = chart seed.
+		- D(S_c,t) = length of the shortest path inside the chart betwen S_c and t.
+		- A_c = chart area.
+	  - Straightness:
+		- P(c,t) = l_out(c,t) / l_in(c,t)
+		- l_out(c,t) = lenght of the edges not shared between C and t.
+		- l_in(c,t) = lenght of the edges shared between C and t.
+
+	- Combined metric:
+	  - Cost(c,t) = F(c,t)^alpha + C(c,t)^beta + P(c,t)^gamma
+	  - alpha = 1, beta = 0.7, gamma = 0.5
+
+	Our basic approach:
+	- Just one iteration of k-means?
+	- Avoid dijkstra by greedily growing charts until a threshold is met. Increase threshold and repeat until no faces left.
+	- If distortion metric is too high, split chart, add two seeds.
+	- If chart size is low, try removing chart.
+
+	Postprocess:
+	- If topology is not disk:
+	  - Fill holes, if new faces fit proxy.
+	  - Find best cut, otherwise.
+	- After parameterization:
+	  - If boundary self-intersects:
+		- cut chart along the closest two diametral boundary vertices, repeat parametrization.
+		- what if the overlap is on an appendix? How do we find that out and cut appropiately?
+		  - emphasize roundness metrics to prevent those cases.
+	  - If interior self-overlaps: preserve boundary parameterization and use mean-value map.
+	*/
+	void computeCharts(const CharterOptions &options, const std::vector<uint32_t> &unchartedMaterialArray) {
+		Chart *vertexMap = NULL;
+		if (unchartedMaterialArray.size() != 0) {
+			vertexMap = new Chart();
+			vertexMap->buildVertexMap(m_mesh, unchartedMaterialArray);
+			if (vertexMap->faceCount() == 0) {
+				delete vertexMap;
+				vertexMap = NULL;
+			}
+		}
+		AtlasBuilder builder(m_mesh);
+		if (vertexMap != NULL) {
+			// Mark faces that do not need to be charted.
+			builder.markUnchartedFaces(vertexMap->faceArray());
+			m_chartArray.push_back(vertexMap);
+		}
+		if (builder.facesLeft != 0) {
+			// Tweak these values:
+			const float maxThreshold = 2;
+			const uint32_t growFaceCount = 32;
+			const uint32_t maxIterations = 4;
+			builder.options = options;
+			//builder.options.proxyFitMetricWeight *= 0.75; // relax proxy fit weight during initial seed placement.
+			//builder.options.roundnessMetricWeight = 0;
+			//builder.options.straightnessMetricWeight = 0;
+			// This seems a reasonable estimate.
+			uint32_t maxSeedCount = std::max(6U, builder.facesLeft);
+			// Create initial charts greedely.
+			xaPrint("### Placing seeds\n");
+			builder.placeSeeds(maxThreshold, maxSeedCount);
+			xaPrint("###   Placed %d seeds (max = %d)\n", builder.chartCount(), maxSeedCount);
+			builder.updateProxies();
+			builder.mergeCharts();
+#if 1
+			xaPrint("### Relocating seeds\n");
+			builder.relocateSeeds();
+			xaPrint("### Reset charts\n");
+			builder.resetCharts();
+			if (vertexMap != NULL) {
+				builder.markUnchartedFaces(vertexMap->faceArray());
+			}
+			builder.options = options;
+			xaPrint("### Growing charts\n");
+			// Restart process growing charts in parallel.
+			uint32_t iteration = 0;
+			while (true) {
+				if (!builder.growCharts(maxThreshold, growFaceCount)) {
+					xaPrint("### Can't grow anymore\n");
+					// If charts cannot grow more: fill holes, merge charts, relocate seeds and start new iteration.
+					xaPrint("### Filling holes\n");
+					builder.fillHoles(maxThreshold);
+					xaPrint("###   Using %d charts now\n", builder.chartCount());
+					builder.updateProxies();
+					xaPrint("### Merging charts\n");
+					builder.mergeCharts();
+					xaPrint("###   Using %d charts now\n", builder.chartCount());
+					xaPrint("### Reseeding\n");
+					if (!builder.relocateSeeds()) {
+						xaPrint("### Cannot relocate seeds anymore\n");
+						// Done!
+						break;
+					}
+					if (iteration == maxIterations) {
+						xaPrint("### Reached iteration limit\n");
+						break;
+					}
+					iteration++;
+					xaPrint("### Reset charts\n");
+					builder.resetCharts();
+					if (vertexMap != NULL) {
+						builder.markUnchartedFaces(vertexMap->faceArray());
+					}
+					xaPrint("### Growing charts\n");
+				}
+			};
+#endif
+			// Make sure no holes are left!
+			xaDebugAssert(builder.facesLeft == 0);
+			const uint32_t chartCount = builder.chartArray.size();
+			for (uint32_t i = 0; i < chartCount; i++) {
+				Chart *chart = new Chart();
+				m_chartArray.push_back(chart);
+				chart->build(m_mesh, builder.chartFaces(i));
+			}
+		}
+		const uint32_t chartCount = m_chartArray.size();
+		// Build face indices.
+		m_faceChart.resize(m_mesh->faceCount());
+		m_faceIndex.resize(m_mesh->faceCount());
+		for (uint32_t i = 0; i < chartCount; i++) {
+			const Chart *chart = m_chartArray[i];
+			const uint32_t faceCount = chart->faceCount();
+			for (uint32_t f = 0; f < faceCount; f++) {
+				uint32_t idx = chart->faceAt(f);
+				m_faceChart[idx] = i;
+				m_faceIndex[idx] = f;
+			}
+		}
+		// Build an exclusive prefix sum of the chart vertex counts.
+		m_chartVertexCountPrefixSum.resize(chartCount);
+		if (chartCount > 0) {
+			m_chartVertexCountPrefixSum[0] = 0;
+			for (uint32_t i = 1; i < chartCount; i++) {
+				const Chart *chart = m_chartArray[i - 1];
+				m_chartVertexCountPrefixSum[i] = m_chartVertexCountPrefixSum[i - 1] + chart->vertexCount();
+			}
+			m_totalVertexCount = m_chartVertexCountPrefixSum[chartCount - 1] + m_chartArray[chartCount - 1]->vertexCount();
+		} else {
+			m_totalVertexCount = 0;
+		}
+	}
+
+	void parameterizeCharts() {
+		ParameterizationQuality globalParameterizationQuality;
+		// Parameterize the charts.
+		uint32_t diskCount = 0;
+		const uint32_t chartCount = m_chartArray.size();
+		for (uint32_t i = 0; i < chartCount; i++) {
+			Chart *chart = m_chartArray[i];
+
+			bool isValid = false;
+
+			if (chart->isVertexMapped()) {
+				continue;
+			}
+
+			if (chart->isDisk()) {
+				diskCount++;
+				ParameterizationQuality chartParameterizationQuality;
+				if (chart->faceCount() == 1) {
+					computeSingleFaceMap(chart->unifiedMesh());
+					chartParameterizationQuality = ParameterizationQuality(chart->unifiedMesh());
+				} else {
+					computeOrthogonalProjectionMap(chart->unifiedMesh());
+					ParameterizationQuality orthogonalQuality(chart->unifiedMesh());
+					computeLeastSquaresConformalMap(chart->unifiedMesh());
+					ParameterizationQuality lscmQuality(chart->unifiedMesh());
+					chartParameterizationQuality = lscmQuality;
+				}
+				isValid = chartParameterizationQuality.isValid();
+				if (!isValid) {
+					xaPrint("*** Invalid parameterization.\n");
+				}
+				// @@ Check that parameterization quality is above a certain threshold.
+				// @@ Detect boundary self-intersections.
+				globalParameterizationQuality += chartParameterizationQuality;
+			}
+
+			// Transfer parameterization from unified mesh to chart mesh.
+			chart->transferParameterization();
+		}
+		xaPrint("  Parameterized %d/%d charts.\n", diskCount, chartCount);
+		xaPrint("  RMS stretch metric: %f\n", globalParameterizationQuality.rmsStretchMetric());
+		xaPrint("  MAX stretch metric: %f\n", globalParameterizationQuality.maxStretchMetric());
+		xaPrint("  RMS conformal metric: %f\n", globalParameterizationQuality.rmsConformalMetric());
+		xaPrint("  RMS authalic metric: %f\n", globalParameterizationQuality.maxAuthalicMetric());
+	}
+
+	uint32_t faceChartAt(uint32_t i) const {
+		return m_faceChart[i];
+	}
+	uint32_t faceIndexWithinChartAt(uint32_t i) const {
+		return m_faceIndex[i];
+	}
+
+	uint32_t vertexCountBeforeChartAt(uint32_t i) const {
+		return m_chartVertexCountPrefixSum[i];
+	}
+
+private:
+	const halfedge::Mesh *m_mesh;
+
+	std::vector<Chart *> m_chartArray;
+
+	std::vector<uint32_t> m_chartVertexCountPrefixSum;
+	uint32_t m_totalVertexCount;
+
+	std::vector<uint32_t> m_faceChart; // the chart of every face of the input mesh.
+	std::vector<uint32_t> m_faceIndex; // the index within the chart for every face of the input mesh.
+};
+
+/// An atlas is a set of charts.
+class Atlas {
+public:
+	~Atlas() {
+		for (size_t i = 0; i < m_meshChartsArray.size(); i++)
+			delete m_meshChartsArray[i];
+	}
+
+	uint32_t meshCount() const {
+		return m_meshChartsArray.size();
+	}
+
+	const MeshCharts *meshAt(uint32_t i) const {
+		return m_meshChartsArray[i];
+	}
+
+	MeshCharts *meshAt(uint32_t i) {
+		return m_meshChartsArray[i];
+	}
+
+	uint32_t chartCount() const {
+		uint32_t count = 0;
+		for (uint32_t c = 0; c < m_meshChartsArray.size(); c++) {
+			count += m_meshChartsArray[c]->chartCount();
+		}
+		return count;
+	}
+
+	const Chart *chartAt(uint32_t i) const {
+		for (uint32_t c = 0; c < m_meshChartsArray.size(); c++) {
+			uint32_t count = m_meshChartsArray[c]->chartCount();
+			if (i < count) {
+				return m_meshChartsArray[c]->chartAt(i);
+			}
+			i -= count;
+		}
+		return NULL;
+	}
+
+	Chart *chartAt(uint32_t i) {
+		for (uint32_t c = 0; c < m_meshChartsArray.size(); c++) {
+			uint32_t count = m_meshChartsArray[c]->chartCount();
+			if (i < count) {
+				return m_meshChartsArray[c]->chartAt(i);
+			}
+			i -= count;
+		}
+		return NULL;
+	}
+
+	// Add mesh charts and takes ownership.
+	// Extract the charts and add to this atlas.
+	void addMeshCharts(MeshCharts *meshCharts) {
+		m_meshChartsArray.push_back(meshCharts);
+	}
+
+	void extractCharts(const halfedge::Mesh *mesh) {
+		MeshCharts *meshCharts = new MeshCharts(mesh);
+		meshCharts->extractCharts();
+		addMeshCharts(meshCharts);
+	}
+
+	void computeCharts(const halfedge::Mesh *mesh, const CharterOptions &options, const std::vector<uint32_t> &unchartedMaterialArray) {
+		MeshCharts *meshCharts = new MeshCharts(mesh);
+		meshCharts->computeCharts(options, unchartedMaterialArray);
+		addMeshCharts(meshCharts);
+	}
+
+	void parameterizeCharts() {
+		for (uint32_t i = 0; i < m_meshChartsArray.size(); i++) {
+			m_meshChartsArray[i]->parameterizeCharts();
+		}
+	}
+
+private:
+	std::vector<MeshCharts *> m_meshChartsArray;
+};
+
+struct AtlasPacker {
+	AtlasPacker(Atlas *atlas) :
+			m_atlas(atlas),
+			m_width(0),
+			m_height(0) {
+		// Save the original uvs.
+		m_originalChartUvs.resize(m_atlas->chartCount());
+		for (uint32_t i = 0; i < m_atlas->chartCount(); i++) {
+			const halfedge::Mesh *mesh = atlas->chartAt(i)->chartMesh();
+			m_originalChartUvs[i].resize(mesh->vertexCount());
+			for (uint32_t j = 0; j < mesh->vertexCount(); j++)
+				m_originalChartUvs[i][j] = mesh->vertexAt(j)->tex;
+		}
+	}
+
+	uint32_t getWidth() const { return m_width; }
+	uint32_t getHeight() const { return m_height; }
+
+	// Pack charts in the smallest possible rectangle.
+	void packCharts(const PackerOptions &options) {
+		const uint32_t chartCount = m_atlas->chartCount();
+		if (chartCount == 0) return;
+		float texelsPerUnit = 1;
+		if (options.method == PackMethod::TexelArea)
+			texelsPerUnit = options.texelArea;
+		for (int iteration = 0;; iteration++) {
+			m_rand = MTRand();
+			std::vector<float> chartOrderArray(chartCount);
+			std::vector<Vector2> chartExtents(chartCount);
+			float meshArea = 0;
+			for (uint32_t c = 0; c < chartCount; c++) {
+				Chart *chart = m_atlas->chartAt(c);
+				if (!chart->isVertexMapped() && !chart->isDisk()) {
+					chartOrderArray[c] = 0;
+					// Skip non-disks.
+					continue;
+				}
+				Vector2 extents(0.0f);
+				if (chart->isVertexMapped()) {
+					// Arrange vertices in a rectangle.
+					extents.x = float(chart->vertexMapWidth);
+					extents.y = float(chart->vertexMapHeight);
+				} else {
+					// Compute surface area to sort charts.
+					float chartArea = chart->computeSurfaceArea();
+					meshArea += chartArea;
+					//chartOrderArray[c] = chartArea;
+					// Compute chart scale
+					float parametricArea = fabsf(chart->computeParametricArea()); // @@ There doesn't seem to be anything preventing parametric area to be negative.
+					if (parametricArea < NV_EPSILON) {
+						// When the parametric area is too small we use a rough approximation to prevent divisions by very small numbers.
+						Vector2 bounds = chart->computeParametricBounds();
+						parametricArea = bounds.x * bounds.y;
+					}
+					float scale = (chartArea / parametricArea) * texelsPerUnit;
+					if (parametricArea == 0) { // < NV_EPSILON)
+						scale = 0;
+					}
+					xaAssert(std::isfinite(scale));
+					// Compute bounding box of chart.
+					Vector2 majorAxis, minorAxis, origin, end;
+					computeBoundingBox(chart, &majorAxis, &minorAxis, &origin, &end);
+					xaAssert(isFinite(majorAxis) && isFinite(minorAxis) && isFinite(origin));
+					// Sort charts by perimeter. @@ This is sometimes producing somewhat unexpected results. Is this right?
+					//chartOrderArray[c] = ((end.x - origin.x) + (end.y - origin.y)) * scale;
+					// Translate, rotate and scale vertices. Compute extents.
+					halfedge::Mesh *mesh = chart->chartMesh();
+					const uint32_t vertexCount = mesh->vertexCount();
+					for (uint32_t i = 0; i < vertexCount; i++) {
+						halfedge::Vertex *vertex = mesh->vertexAt(i);
+						//Vector2 t = vertex->tex - origin;
+						Vector2 tmp;
+						tmp.x = dot(vertex->tex, majorAxis);
+						tmp.y = dot(vertex->tex, minorAxis);
+						tmp -= origin;
+						tmp *= scale;
+						if (tmp.x < 0 || tmp.y < 0) {
+							xaPrint("tmp: %f %f\n", tmp.x, tmp.y);
+							xaPrint("scale: %f\n", scale);
+							xaPrint("origin: %f %f\n", origin.x, origin.y);
+							xaPrint("majorAxis: %f %f\n", majorAxis.x, majorAxis.y);
+							xaPrint("minorAxis: %f %f\n", minorAxis.x, minorAxis.y);
+							xaDebugAssert(false);
+						}
+						//xaAssert(tmp.x >= 0 && tmp.y >= 0);
+						vertex->tex = tmp;
+						xaAssert(std::isfinite(vertex->tex.x) && std::isfinite(vertex->tex.y));
+						extents = max(extents, tmp);
+					}
+					xaDebugAssert(extents.x >= 0 && extents.y >= 0);
+					// Limit chart size.
+					if (extents.x > 1024 || extents.y > 1024) {
+						float limit = std::max(extents.x, extents.y);
+						scale = 1024 / (limit + 1);
+						for (uint32_t i = 0; i < vertexCount; i++) {
+							halfedge::Vertex *vertex = mesh->vertexAt(i);
+							vertex->tex *= scale;
+						}
+						extents *= scale;
+						xaDebugAssert(extents.x <= 1024 && extents.y <= 1024);
+					}
+					// Scale the charts to use the entire texel area available. So, if the width is 0.1 we could scale it to 1 without increasing the lightmap usage and making a better
+					// use of it. In many cases this also improves the look of the seams, since vertices on the chart boundaries have more chances of being aligned with the texel centers.
+					float scale_x = 1.0f;
+					float scale_y = 1.0f;
+					float divide_x = 1.0f;
+					float divide_y = 1.0f;
+					if (extents.x > 0) {
+						int cw = ftoi_ceil(extents.x);
+						if (options.blockAlign && chart->blockAligned) {
+							// Align all chart extents to 4x4 blocks, but taking padding into account.
+							if (options.conservative) {
+								cw = align(cw + 2, 4) - 2;
+							} else {
+								cw = align(cw + 1, 4) - 1;
+							}
+						}
+						scale_x = (float(cw) - NV_EPSILON);
+						divide_x = extents.x;
+						extents.x = float(cw);
+					}
+					if (extents.y > 0) {
+						int ch = ftoi_ceil(extents.y);
+						if (options.blockAlign && chart->blockAligned) {
+							// Align all chart extents to 4x4 blocks, but taking padding into account.
+							if (options.conservative) {
+								ch = align(ch + 2, 4) - 2;
+							} else {
+								ch = align(ch + 1, 4) - 1;
+							}
+						}
+						scale_y = (float(ch) - NV_EPSILON);
+						divide_y = extents.y;
+						extents.y = float(ch);
+					}
+					for (uint32_t v = 0; v < vertexCount; v++) {
+						halfedge::Vertex *vertex = mesh->vertexAt(v);
+						vertex->tex.x /= divide_x;
+						vertex->tex.y /= divide_y;
+						vertex->tex.x *= scale_x;
+						vertex->tex.y *= scale_y;
+						xaAssert(std::isfinite(vertex->tex.x) && std::isfinite(vertex->tex.y));
+					}
+				}
+				chartExtents[c] = extents;
+				// Sort charts by perimeter.
+				chartOrderArray[c] = extents.x + extents.y;
+			}
+			// @@ We can try to improve compression of small charts by sorting them by proximity like we do with vertex samples.
+			// @@ How to do that? One idea: compute chart centroid, insert into grid, compute morton index of the cell, sort based on morton index.
+			// @@ We would sort by morton index, first, then quantize the chart sizes, so that all small charts have the same size, and sort by size preserving the morton order.
+			//xaPrint("Sorting charts.\n");
+			// Sort charts by area.
+			m_radix = RadixSort();
+			m_radix.sort(chartOrderArray);
+			const uint32_t *ranks = m_radix.ranks();
+			// First iteration - guess texelsPerUnit.
+			if (options.method != PackMethod::TexelArea && iteration == 0) {
+				// Estimate size of the map based on the mesh surface area and given texel scale.
+				const float texelCount = std::max(1.0f, meshArea * square(texelsPerUnit) / 0.75f); // Assume 75% utilization.
+				texelsPerUnit = sqrt((options.resolution * options.resolution) / texelCount);
+				resetUvs();
+				continue;
+			}
+			// Init bit map.
+			m_bitmap.clearAll();
+			m_bitmap.resize(options.resolution, options.resolution, false);
+			int w = 0;
+			int h = 0;
+			// Add sorted charts to bitmap.
+			for (uint32_t i = 0; i < chartCount; i++) {
+				uint32_t c = ranks[chartCount - i - 1]; // largest chart first
+				Chart *chart = m_atlas->chartAt(c);
+				if (!chart->isVertexMapped() && !chart->isDisk()) continue;
+				//float scale_x = 1;
+				//float scale_y = 1;
+				BitMap chart_bitmap;
+				if (chart->isVertexMapped()) {
+					chart->blockAligned = false;
+					// Init all bits to 1.
+					chart_bitmap.resize(ftoi_ceil(chartExtents[c].x), ftoi_ceil(chartExtents[c].y), /*initValue=*/true);
+					// @@ Another alternative would be to try to map each vertex to a different texel trying to fill all the available unused texels.
+				} else {
+					// @@ Add special cases for dot and line charts. @@ Lightmap rasterizer also needs to handle these special cases.
+					// @@ We could also have a special case for chart quads. If the quad surface <= 4 texels, align vertices with texel centers and do not add padding. May be very useful for foliage.
+					// @@ In general we could reduce the padding of all charts by one texel by using a rasterizer that takes into account the 2-texel footprint of the tent bilinear filter. For example,
+					// if we have a chart that is less than 1 texel wide currently we add one texel to the left and one texel to the right creating a 3-texel-wide bitmap. However, if we know that the
+					// chart is only 1 texel wide we could align it so that it only touches the footprint of two texels:
+					//      |   |      <- Touches texels 0, 1 and 2.
+					//    |   |        <- Only touches texels 0 and 1.
+					// \   \ / \ /   /
+					//  \   X   X   /
+					//   \ / \ / \ /
+					//    V   V   V
+					//    0   1   2
+					if (options.conservative) {
+						// Init all bits to 0.
+						chart_bitmap.resize(ftoi_ceil(chartExtents[c].x) + 1 + options.padding, ftoi_ceil(chartExtents[c].y) + 1 + options.padding, /*initValue=*/false); // + 2 to add padding on both sides.
+						// Rasterize chart and dilate.
+						drawChartBitmapDilate(chart, &chart_bitmap, options.padding);
+					} else {
+						// Init all bits to 0.
+						chart_bitmap.resize(ftoi_ceil(chartExtents[c].x) + 1, ftoi_ceil(chartExtents[c].y) + 1, /*initValue=*/false); // Add half a texels on each side.
+						// Rasterize chart and dilate.
+						drawChartBitmap(chart, &chart_bitmap, Vector2(1), Vector2(0.5));
+					}
+				}
+				int best_x, best_y;
+				int best_cw, best_ch; // Includes padding now.
+				int best_r;
+				findChartLocation(options.quality, &chart_bitmap, chartExtents[c], w, h, &best_x, &best_y, &best_cw, &best_ch, &best_r, chart->blockAligned);
+				/*if (w < best_x + best_cw || h < best_y + best_ch)
+				{
+					xaPrint("Resize extents to (%d, %d).\n", best_x + best_cw, best_y + best_ch);
+				}*/
+				// Update parametric extents.
+				w = std::max(w, best_x + best_cw);
+				h = std::max(h, best_y + best_ch);
+				w = align(w, 4);
+				h = align(h, 4);
+				// Resize bitmap if necessary.
+				if (uint32_t(w) > m_bitmap.width() || uint32_t(h) > m_bitmap.height()) {
+					//xaPrint("Resize bitmap (%d, %d).\n", nextPowerOfTwo(w), nextPowerOfTwo(h));
+					m_bitmap.resize(nextPowerOfTwo(uint32_t(w)), nextPowerOfTwo(uint32_t(h)), false);
+				}
+				//xaPrint("Add chart at (%d, %d).\n", best_x, best_y);
+				addChart(&chart_bitmap, w, h, best_x, best_y, best_r);
+				//float best_angle = 2 * PI * best_r;
+				// Translate and rotate chart texture coordinates.
+				halfedge::Mesh *mesh = chart->chartMesh();
+				const uint32_t vertexCount = mesh->vertexCount();
+				for (uint32_t v = 0; v < vertexCount; v++) {
+					halfedge::Vertex *vertex = mesh->vertexAt(v);
+					Vector2 t = vertex->tex;
+					if (best_r) std::swap(t.x, t.y);
+					//vertex->tex.x = best_x + t.x * cosf(best_angle) - t.y * sinf(best_angle);
+					//vertex->tex.y = best_y + t.x * sinf(best_angle) + t.y * cosf(best_angle);
+					vertex->tex.x = best_x + t.x + 0.5f;
+					vertex->tex.y = best_y + t.y + 0.5f;
+					xaAssert(vertex->tex.x >= 0 && vertex->tex.y >= 0);
+					xaAssert(std::isfinite(vertex->tex.x) && std::isfinite(vertex->tex.y));
+				}
+			}
+			//w -= padding - 1; // Leave one pixel border!
+			//h -= padding - 1;
+			m_width = std::max(0, w);
+			m_height = std::max(0, h);
+			xaAssert(isAligned(m_width, 4));
+			xaAssert(isAligned(m_height, 4));
+			if (options.method == PackMethod::ExactResolution) {
+				texelsPerUnit *= sqrt((options.resolution * options.resolution) / (float)(m_width * m_height));
+				if (iteration > 1 && m_width <= options.resolution && m_height <= options.resolution) {
+					m_width = m_height = options.resolution;
+					return;
+				}
+				resetUvs();
+			} else {
+				return;
+			}
+		}
+	}
+
+	float computeAtlasUtilization() const {
+		const uint32_t w = m_width;
+		const uint32_t h = m_height;
+		xaDebugAssert(w <= m_bitmap.width());
+		xaDebugAssert(h <= m_bitmap.height());
+		uint32_t count = 0;
+		for (uint32_t y = 0; y < h; y++) {
+			for (uint32_t x = 0; x < w; x++) {
+				count += m_bitmap.bitAt(x, y);
+			}
+		}
+		return float(count) / (w * h);
+	}
+
+private:
+	void resetUvs() {
+		for (uint32_t i = 0; i < m_atlas->chartCount(); i++) {
+			halfedge::Mesh *mesh = m_atlas->chartAt(i)->chartMesh();
+			for (uint32_t j = 0; j < mesh->vertexCount(); j++)
+				mesh->vertexAt(j)->tex = m_originalChartUvs[i][j];
+		}
+	}
+
+	// IC: Brute force is slow, and random may take too much time to converge. We start inserting large charts in a small atlas. Using brute force is lame, because most of the space
+	// is occupied at this point. At the end we have many small charts and a large atlas with sparse holes. Finding those holes randomly is slow. A better approach would be to
+	// start stacking large charts as if they were tetris pieces. Once charts get small try to place them randomly. It may be interesting to try a intermediate strategy, first try
+	// along one axis and then try exhaustively along that axis.
+	void findChartLocation(int quality, const BitMap *bitmap, Vector2::Arg extents, int w, int h, int *best_x, int *best_y, int *best_w, int *best_h, int *best_r, bool blockAligned) {
+		int attempts = 256;
+		if (quality == 1) attempts = 4096;
+		if (quality == 2) attempts = 2048;
+		if (quality == 3) attempts = 1024;
+		if (quality == 4) attempts = 512;
+		if (quality == 0 || w * h < attempts) {
+			findChartLocation_bruteForce(bitmap, extents, w, h, best_x, best_y, best_w, best_h, best_r, blockAligned);
+		} else {
+			findChartLocation_random(bitmap, extents, w, h, best_x, best_y, best_w, best_h, best_r, attempts, blockAligned);
+		}
+	}
+
+	void findChartLocation_bruteForce(const BitMap *bitmap, Vector2::Arg /*extents*/, int w, int h, int *best_x, int *best_y, int *best_w, int *best_h, int *best_r, bool blockAligned) {
+		const int BLOCK_SIZE = 4;
+		int best_metric = INT_MAX;
+		int step_size = blockAligned ? BLOCK_SIZE : 1;
+		// Try two different orientations.
+		for (int r = 0; r < 2; r++) {
+			int cw = bitmap->width();
+			int ch = bitmap->height();
+			if (r & 1) std::swap(cw, ch);
+			for (int y = 0; y <= h + 1; y += step_size) { // + 1 to extend atlas in case atlas full.
+				for (int x = 0; x <= w + 1; x += step_size) { // + 1 not really necessary here.
+					// Early out.
+					int area = std::max(w, x + cw) * std::max(h, y + ch);
+					//int perimeter = max(w, x+cw) + max(h, y+ch);
+					int extents = std::max(std::max(w, x + cw), std::max(h, y + ch));
+					int metric = extents * extents + area;
+					if (metric > best_metric) {
+						continue;
+					}
+					if (metric == best_metric && std::max(x, y) >= std::max(*best_x, *best_y)) {
+						// If metric is the same, pick the one closest to the origin.
+						continue;
+					}
+					if (canAddChart(bitmap, w, h, x, y, r)) {
+						best_metric = metric;
+						*best_x = x;
+						*best_y = y;
+						*best_w = cw;
+						*best_h = ch;
+						*best_r = r;
+						if (area == w * h) {
+							// Chart is completely inside, do not look at any other location.
+							goto done;
+						}
+					}
+				}
+			}
+		}
+	done:
+		xaDebugAssert(best_metric != INT_MAX);
+	}
+
+	void findChartLocation_random(const BitMap *bitmap, Vector2::Arg /*extents*/, int w, int h, int *best_x, int *best_y, int *best_w, int *best_h, int *best_r, int minTrialCount, bool blockAligned) {
+		const int BLOCK_SIZE = 4;
+		int best_metric = INT_MAX;
+		for (int i = 0; i < minTrialCount || best_metric == INT_MAX; i++) {
+			int r = m_rand.getRange(1);
+			int x = m_rand.getRange(w + 1); // + 1 to extend atlas in case atlas full. We may want to use a higher number to increase probability of extending atlas.
+			int y = m_rand.getRange(h + 1); // + 1 to extend atlas in case atlas full.
+			if (blockAligned) {
+				x = align(x, BLOCK_SIZE);
+				y = align(y, BLOCK_SIZE);
+			}
+			int cw = bitmap->width();
+			int ch = bitmap->height();
+			if (r & 1) std::swap(cw, ch);
+			// Early out.
+			int area = std::max(w, x + cw) * std::max(h, y + ch);
+			//int perimeter = max(w, x+cw) + max(h, y+ch);
+			int extents = std::max(std::max(w, x + cw), std::max(h, y + ch));
+			int metric = extents * extents + area;
+			if (metric > best_metric) {
+				continue;
+			}
+			if (metric == best_metric && std::min(x, y) > std::min(*best_x, *best_y)) {
+				// If metric is the same, pick the one closest to the origin.
+				continue;
+			}
+			if (canAddChart(bitmap, w, h, x, y, r)) {
+				best_metric = metric;
+				*best_x = x;
+				*best_y = y;
+				*best_w = cw;
+				*best_h = ch;
+				*best_r = r;
+				if (area == w * h) {
+					// Chart is completely inside, do not look at any other location.
+					break;
+				}
+			}
+		}
+	}
+
+	void drawChartBitmapDilate(const Chart *chart, BitMap *bitmap, int padding) {
+		const int w = bitmap->width();
+		const int h = bitmap->height();
+		const Vector2 extents = Vector2(float(w), float(h));
+		// Rasterize chart faces, check that all bits are not set.
+		const uint32_t faceCount = chart->faceCount();
+		for (uint32_t f = 0; f < faceCount; f++) {
+			const halfedge::Face *face = chart->chartMesh()->faceAt(f);
+			Vector2 vertices[4];
+			uint32_t edgeCount = 0;
+			for (halfedge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) {
+				if (edgeCount < 4) {
+					vertices[edgeCount] = it.vertex()->tex + Vector2(0.5) + Vector2(float(padding), float(padding));
+				}
+				edgeCount++;
+			}
+			if (edgeCount == 3) {
+				raster::drawTriangle(raster::Mode_Antialiased, extents, true, vertices, AtlasPacker::setBitsCallback, bitmap);
+			} else {
+				raster::drawQuad(raster::Mode_Antialiased, extents, true, vertices, AtlasPacker::setBitsCallback, bitmap);
+			}
+		}
+		// Expand chart by padding pixels. (dilation)
+		BitMap tmp(w, h);
+		for (int i = 0; i < padding; i++) {
+			tmp.clearAll();
+			for (int y = 0; y < h; y++) {
+				for (int x = 0; x < w; x++) {
+					bool b = bitmap->bitAt(x, y);
+					if (!b) {
+						if (x > 0) {
+							b |= bitmap->bitAt(x - 1, y);
+							if (y > 0) b |= bitmap->bitAt(x - 1, y - 1);
+							if (y < h - 1) b |= bitmap->bitAt(x - 1, y + 1);
+						}
+						if (y > 0) b |= bitmap->bitAt(x, y - 1);
+						if (y < h - 1) b |= bitmap->bitAt(x, y + 1);
+						if (x < w - 1) {
+							b |= bitmap->bitAt(x + 1, y);
+							if (y > 0) b |= bitmap->bitAt(x + 1, y - 1);
+							if (y < h - 1) b |= bitmap->bitAt(x + 1, y + 1);
+						}
+					}
+					if (b) tmp.setBitAt(x, y);
+				}
+			}
+			std::swap(tmp, *bitmap);
+		}
+	}
+
+	void drawChartBitmap(const Chart *chart, BitMap *bitmap, const Vector2 &scale, const Vector2 &offset) {
+		const int w = bitmap->width();
+		const int h = bitmap->height();
+		const Vector2 extents = Vector2(float(w), float(h));
+		static const Vector2 pad[4] = {
+			Vector2(-0.5, -0.5),
+			Vector2(0.5, -0.5),
+			Vector2(-0.5, 0.5),
+			Vector2(0.5, 0.5)
+		};
+		// Rasterize 4 times to add proper padding.
+		for (int i = 0; i < 4; i++) {
+			// Rasterize chart faces, check that all bits are not set.
+			const uint32_t faceCount = chart->chartMesh()->faceCount();
+			for (uint32_t f = 0; f < faceCount; f++) {
+				const halfedge::Face *face = chart->chartMesh()->faceAt(f);
+				Vector2 vertices[4];
+				uint32_t edgeCount = 0;
+				for (halfedge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) {
+					if (edgeCount < 4) {
+						vertices[edgeCount] = it.vertex()->tex * scale + offset + pad[i];
+						xaAssert(ftoi_ceil(vertices[edgeCount].x) >= 0);
+						xaAssert(ftoi_ceil(vertices[edgeCount].y) >= 0);
+						xaAssert(ftoi_ceil(vertices[edgeCount].x) <= w);
+						xaAssert(ftoi_ceil(vertices[edgeCount].y) <= h);
+					}
+					edgeCount++;
+				}
+				if (edgeCount == 3) {
+					raster::drawTriangle(raster::Mode_Antialiased, extents, /*enableScissors=*/true, vertices, AtlasPacker::setBitsCallback, bitmap);
+				} else {
+					raster::drawQuad(raster::Mode_Antialiased, extents, /*enableScissors=*/true, vertices, AtlasPacker::setBitsCallback, bitmap);
+				}
+			}
+		}
+		// Expand chart by padding pixels. (dilation)
+		BitMap tmp(w, h);
+		tmp.clearAll();
+		for (int y = 0; y < h; y++) {
+			for (int x = 0; x < w; x++) {
+				bool b = bitmap->bitAt(x, y);
+				if (!b) {
+					if (x > 0) {
+						b |= bitmap->bitAt(x - 1, y);
+						if (y > 0) b |= bitmap->bitAt(x - 1, y - 1);
+						if (y < h - 1) b |= bitmap->bitAt(x - 1, y + 1);
+					}
+					if (y > 0) b |= bitmap->bitAt(x, y - 1);
+					if (y < h - 1) b |= bitmap->bitAt(x, y + 1);
+					if (x < w - 1) {
+						b |= bitmap->bitAt(x + 1, y);
+						if (y > 0) b |= bitmap->bitAt(x + 1, y - 1);
+						if (y < h - 1) b |= bitmap->bitAt(x + 1, y + 1);
+					}
+				}
+				if (b) tmp.setBitAt(x, y);
+			}
+		}
+		std::swap(tmp, *bitmap);
+	}
+
+	bool canAddChart(const BitMap *bitmap, int atlas_w, int atlas_h, int offset_x, int offset_y, int r) {
+		xaDebugAssert(r == 0 || r == 1);
+		// Check whether the two bitmaps overlap.
+		const int w = bitmap->width();
+		const int h = bitmap->height();
+		if (r == 0) {
+			for (int y = 0; y < h; y++) {
+				int yy = y + offset_y;
+				if (yy >= 0) {
+					for (int x = 0; x < w; x++) {
+						int xx = x + offset_x;
+						if (xx >= 0) {
+							if (bitmap->bitAt(x, y)) {
+								if (xx < atlas_w && yy < atlas_h) {
+									if (m_bitmap.bitAt(xx, yy)) return false;
+								}
+							}
+						}
+					}
+				}
+			}
+		} else if (r == 1) {
+			for (int y = 0; y < h; y++) {
+				int xx = y + offset_x;
+				if (xx >= 0) {
+					for (int x = 0; x < w; x++) {
+						int yy = x + offset_y;
+						if (yy >= 0) {
+							if (bitmap->bitAt(x, y)) {
+								if (xx < atlas_w && yy < atlas_h) {
+									if (m_bitmap.bitAt(xx, yy)) return false;
+								}
+							}
+						}
+					}
+				}
+			}
+		}
+		return true;
+	}
+
+	void addChart(const BitMap *bitmap, int atlas_w, int atlas_h, int offset_x, int offset_y, int r) {
+		xaDebugAssert(r == 0 || r == 1);
+		// Check whether the two bitmaps overlap.
+		const int w = bitmap->width();
+		const int h = bitmap->height();
+		if (r == 0) {
+			for (int y = 0; y < h; y++) {
+				int yy = y + offset_y;
+				if (yy >= 0) {
+					for (int x = 0; x < w; x++) {
+						int xx = x + offset_x;
+						if (xx >= 0) {
+							if (bitmap->bitAt(x, y)) {
+								if (xx < atlas_w && yy < atlas_h) {
+									xaDebugAssert(m_bitmap.bitAt(xx, yy) == false);
+									m_bitmap.setBitAt(xx, yy);
+								}
+							}
+						}
+					}
+				}
+			}
+		} else if (r == 1) {
+			for (int y = 0; y < h; y++) {
+				int xx = y + offset_x;
+				if (xx >= 0) {
+					for (int x = 0; x < w; x++) {
+						int yy = x + offset_y;
+						if (yy >= 0) {
+							if (bitmap->bitAt(x, y)) {
+								if (xx < atlas_w && yy < atlas_h) {
+									xaDebugAssert(m_bitmap.bitAt(xx, yy) == false);
+									m_bitmap.setBitAt(xx, yy);
+								}
+							}
+						}
+					}
+				}
+			}
+		}
+	}
+
+	static bool setBitsCallback(void *param, int x, int y, Vector3::Arg, Vector3::Arg, Vector3::Arg, float area) {
+		BitMap *bitmap = (BitMap *)param;
+		if (area > 0.0) {
+			bitmap->setBitAt(x, y);
+		}
+		return true;
+	}
+
+	// Compute the convex hull using Graham Scan.
+	static void convexHull(const std::vector<Vector2> &input, std::vector<Vector2> &output, float epsilon) {
+		const uint32_t inputCount = input.size();
+		std::vector<float> coords(inputCount);
+		for (uint32_t i = 0; i < inputCount; i++) {
+			coords[i] = input[i].x;
+		}
+		RadixSort radix;
+		radix.sort(coords);
+		const uint32_t *ranks = radix.ranks();
+		std::vector<Vector2> top;
+		top.reserve(inputCount);
+		std::vector<Vector2> bottom;
+		bottom.reserve(inputCount);
+		Vector2 P = input[ranks[0]];
+		Vector2 Q = input[ranks[inputCount - 1]];
+		float topy = std::max(P.y, Q.y);
+		float boty = std::min(P.y, Q.y);
+		for (uint32_t i = 0; i < inputCount; i++) {
+			Vector2 p = input[ranks[i]];
+			if (p.y >= boty) top.push_back(p);
+		}
+		for (uint32_t i = 0; i < inputCount; i++) {
+			Vector2 p = input[ranks[inputCount - 1 - i]];
+			if (p.y <= topy) bottom.push_back(p);
+		}
+		// Filter top list.
+		output.clear();
+		output.push_back(top[0]);
+		output.push_back(top[1]);
+		for (uint32_t i = 2; i < top.size();) {
+			Vector2 a = output[output.size() - 2];
+			Vector2 b = output[output.size() - 1];
+			Vector2 c = top[i];
+			float area = triangleArea(a, b, c);
+			if (area >= -epsilon) {
+				output.pop_back();
+			}
+			if (area < -epsilon || output.size() == 1) {
+				output.push_back(c);
+				i++;
+			}
+		}
+		uint32_t top_count = output.size();
+		output.push_back(bottom[1]);
+		// Filter bottom list.
+		for (uint32_t i = 2; i < bottom.size();) {
+			Vector2 a = output[output.size() - 2];
+			Vector2 b = output[output.size() - 1];
+			Vector2 c = bottom[i];
+			float area = triangleArea(a, b, c);
+			if (area >= -epsilon) {
+				output.pop_back();
+			}
+			if (area < -epsilon || output.size() == top_count) {
+				output.push_back(c);
+				i++;
+			}
+		}
+		// Remove duplicate element.
+		xaDebugAssert(output.front() == output.back());
+		output.pop_back();
+	}
+
+	// This should compute convex hull and use rotating calipers to find the best box. Currently it uses a brute force method.
+	static void computeBoundingBox(Chart *chart, Vector2 *majorAxis, Vector2 *minorAxis, Vector2 *minCorner, Vector2 *maxCorner) {
+		// Compute list of boundary points.
+		std::vector<Vector2> points;
+		points.reserve(16);
+		halfedge::Mesh *mesh = chart->chartMesh();
+		const uint32_t vertexCount = mesh->vertexCount();
+		for (uint32_t i = 0; i < vertexCount; i++) {
+			halfedge::Vertex *vertex = mesh->vertexAt(i);
+			if (vertex->isBoundary()) {
+				points.push_back(vertex->tex);
+			}
+		}
+		xaDebugAssert(points.size() > 0);
+		std::vector<Vector2> hull;
+		convexHull(points, hull, 0.00001f);
+		// @@ Ideally I should use rotating calipers to find the best box. Using brute force for now.
+		float best_area = FLT_MAX;
+		Vector2 best_min;
+		Vector2 best_max;
+		Vector2 best_axis;
+		const uint32_t hullCount = hull.size();
+		for (uint32_t i = 0, j = hullCount - 1; i < hullCount; j = i, i++) {
+			if (equal(hull[i], hull[j])) {
+				continue;
+			}
+			Vector2 axis = normalize(hull[i] - hull[j], 0.0f);
+			xaDebugAssert(isFinite(axis));
+			// Compute bounding box.
+			Vector2 box_min(FLT_MAX, FLT_MAX);
+			Vector2 box_max(-FLT_MAX, -FLT_MAX);
+			for (uint32_t v = 0; v < hullCount; v++) {
+				Vector2 point = hull[v];
+				float x = dot(axis, point);
+				if (x < box_min.x) box_min.x = x;
+				if (x > box_max.x) box_max.x = x;
+				float y = dot(Vector2(-axis.y, axis.x), point);
+				if (y < box_min.y) box_min.y = y;
+				if (y > box_max.y) box_max.y = y;
+			}
+			// Compute box area.
+			float area = (box_max.x - box_min.x) * (box_max.y - box_min.y);
+			if (area < best_area) {
+				best_area = area;
+				best_min = box_min;
+				best_max = box_max;
+				best_axis = axis;
+			}
+		}
+		// Consider all points, not only boundary points, in case the input chart is malformed.
+		for (uint32_t i = 0; i < vertexCount; i++) {
+			halfedge::Vertex *vertex = mesh->vertexAt(i);
+			Vector2 point = vertex->tex;
+			float x = dot(best_axis, point);
+			if (x < best_min.x) best_min.x = x;
+			if (x > best_max.x) best_max.x = x;
+			float y = dot(Vector2(-best_axis.y, best_axis.x), point);
+			if (y < best_min.y) best_min.y = y;
+			if (y > best_max.y) best_max.y = y;
+		}
+		*majorAxis = best_axis;
+		*minorAxis = Vector2(-best_axis.y, best_axis.x);
+		*minCorner = best_min;
+		*maxCorner = best_max;
+	}
+
+	Atlas *m_atlas;
+	BitMap m_bitmap;
+	RadixSort m_radix;
+	uint32_t m_width;
+	uint32_t m_height;
+	MTRand m_rand;
+	std::vector<std::vector<Vector2> > m_originalChartUvs;
+};
+
+} // namespace param
+} // namespace internal
+
+struct Atlas {
+	internal::param::Atlas atlas;
+	std::vector<internal::halfedge::Mesh *> heMeshes;
+	uint32_t width = 0;
+	uint32_t height = 0;
+	OutputMesh **outputMeshes = NULL;
+};
+
+void SetPrint(PrintFunc print) {
+	internal::s_print = print;
+}
+
+Atlas *Create() {
+	Atlas *atlas = new Atlas();
+	return atlas;
+}
+
+void Destroy(Atlas *atlas) {
+	xaAssert(atlas);
+	for (int i = 0; i < (int)atlas->heMeshes.size(); i++) {
+		delete atlas->heMeshes[i];
+		if (atlas->outputMeshes) {
+			OutputMesh *outputMesh = atlas->outputMeshes[i];
+			for (uint32_t j = 0; j < outputMesh->chartCount; j++)
+				delete[] outputMesh->chartArray[j].indexArray;
+			delete[] outputMesh->chartArray;
+			delete[] outputMesh->vertexArray;
+			delete[] outputMesh->indexArray;
+			delete outputMesh;
+		}
+	}
+	delete[] atlas->outputMeshes;
+	delete atlas;
+}
+
+static internal::Vector3 DecodePosition(const InputMesh &mesh, uint32_t index) {
+	xaAssert(mesh.vertexPositionData);
+	return *((const internal::Vector3 *)&((const uint8_t *)mesh.vertexPositionData)[mesh.vertexPositionStride * index]);
+}
+
+static internal::Vector3 DecodeNormal(const InputMesh &mesh, uint32_t index) {
+	xaAssert(mesh.vertexNormalData);
+	return *((const internal::Vector3 *)&((const uint8_t *)mesh.vertexNormalData)[mesh.vertexNormalStride * index]);
+}
+
+static internal::Vector2 DecodeUv(const InputMesh &mesh, uint32_t index) {
+	xaAssert(mesh.vertexUvData);
+	return *((const internal::Vector2 *)&((const uint8_t *)mesh.vertexUvData)[mesh.vertexUvStride * index]);
+}
+
+static uint32_t DecodeIndex(IndexFormat::Enum format, const void *indexData, uint32_t i) {
+	if (format == IndexFormat::HalfFloat)
+		return (uint32_t)((const uint16_t *)indexData)[i];
+	return ((const uint32_t *)indexData)[i];
+}
+
+static float EdgeLength(internal::Vector3 pos1, internal::Vector3 pos2) {
+	return internal::length(pos2 - pos1);
+}
+
+AddMeshError AddMesh(Atlas *atlas, const InputMesh &mesh, bool useColocalVertices) {
+	xaAssert(atlas);
+	AddMeshError error;
+	error.code = AddMeshErrorCode::Success;
+	error.face = error.index0 = error.index1 = UINT32_MAX;
+	// Expecting triangle faces.
+	if ((mesh.indexCount % 3) != 0) {
+		error.code = AddMeshErrorCode::InvalidIndexCount;
+		return error;
+	}
+	// Check if any index is out of range.
+	for (uint32_t j = 0; j < mesh.indexCount; j++) {
+		const uint32_t index = DecodeIndex(mesh.indexFormat, mesh.indexData, j);
+		if (index < 0 || index >= mesh.vertexCount) {
+			error.code = AddMeshErrorCode::IndexOutOfRange;
+			error.index0 = index;
+			return error;
+		}
+	}
+	// Build half edge mesh.
+	internal::halfedge::Mesh *heMesh = new internal::halfedge::Mesh;
+	std::vector<uint32_t> canonicalMap;
+	canonicalMap.reserve(mesh.vertexCount);
+	for (uint32_t i = 0; i < mesh.vertexCount; i++) {
+		internal::halfedge::Vertex *vertex = heMesh->addVertex(DecodePosition(mesh, i));
+		if (mesh.vertexNormalData)
+			vertex->nor = DecodeNormal(mesh, i);
+		if (mesh.vertexUvData)
+			vertex->tex = DecodeUv(mesh, i);
+		// Link colocals. You probably want to do this more efficiently! Sort by one axis or use a hash or grid.
+		uint32_t firstColocal = i;
+		if (useColocalVertices) {
+			for (uint32_t j = 0; j < i; j++) {
+				if (vertex->pos != DecodePosition(mesh, j))
+					continue;
+#if 0
+				if (mesh.vertexNormalData && vertex->nor != DecodeNormal(mesh, j))
+					continue;
+#endif
+				if (mesh.vertexUvData && vertex->tex != DecodeUv(mesh, j))
+					continue;
+				firstColocal = j;
+				break;
+			}
+		}
+		canonicalMap.push_back(firstColocal);
+	}
+	heMesh->linkColocalsWithCanonicalMap(canonicalMap);
+	for (uint32_t i = 0; i < mesh.indexCount / 3; i++) {
+		uint32_t tri[3];
+		for (int j = 0; j < 3; j++)
+			tri[j] = DecodeIndex(mesh.indexFormat, mesh.indexData, i * 3 + j);
+		// Check for zero length edges.
+		for (int j = 0; j < 3; j++) {
+			const uint32_t edges[6] = { 0, 1, 1, 2, 2, 0 };
+			const uint32_t index1 = tri[edges[j * 2 + 0]];
+			const uint32_t index2 = tri[edges[j * 2 + 1]];
+			const internal::Vector3 pos1 = DecodePosition(mesh, index1);
+			const internal::Vector3 pos2 = DecodePosition(mesh, index2);
+			if (EdgeLength(pos1, pos2) <= 0.0f) {
+				delete heMesh;
+				error.code = AddMeshErrorCode::ZeroLengthEdge;
+				error.face = i;
+				error.index0 = index1;
+				error.index1 = index2;
+				return error;
+			}
+		}
+		// Check for zero area faces.
+		{
+			const internal::Vector3 a = DecodePosition(mesh, tri[0]);
+			const internal::Vector3 b = DecodePosition(mesh, tri[1]);
+			const internal::Vector3 c = DecodePosition(mesh, tri[2]);
+			const float area = internal::length(internal::cross(b - a, c - a)) * 0.5f;
+			if (area <= 0.0f) {
+				delete heMesh;
+				error.code = AddMeshErrorCode::ZeroAreaFace;
+				error.face = i;
+				return error;
+			}
+		}
+		internal::halfedge::Face *face = heMesh->addFace(tri[0], tri[1], tri[2]);
+
+		if (!face && heMesh->errorCode == internal::halfedge::Mesh::ErrorCode::AlreadyAddedEdge) {
+			//there is still hope for this, no reason to not add, at least add as separate
+			face = heMesh->addUniqueFace(tri[0], tri[1], tri[2]);
+		}
+
+		if (!face) {
+			//continue;
+
+			if (heMesh->errorCode == internal::halfedge::Mesh::ErrorCode::AlreadyAddedEdge) {
+				error.code = AddMeshErrorCode::AlreadyAddedEdge;
+			} else if (heMesh->errorCode == internal::halfedge::Mesh::ErrorCode::DegenerateColocalEdge) {
+				error.code = AddMeshErrorCode::DegenerateColocalEdge;
+			} else if (heMesh->errorCode == internal::halfedge::Mesh::ErrorCode::DegenerateEdge) {
+				error.code = AddMeshErrorCode::DegenerateEdge;
+			} else if (heMesh->errorCode == internal::halfedge::Mesh::ErrorCode::DuplicateEdge) {
+				error.code = AddMeshErrorCode::DuplicateEdge;
+			}
+			error.face = i;
+			error.index0 = heMesh->errorIndex0;
+			error.index1 = heMesh->errorIndex1;
+			delete heMesh;
+			return error;
+		}
+		if (mesh.faceMaterialData)
+			face->material = mesh.faceMaterialData[i];
+	}
+	heMesh->linkBoundary();
+	atlas->heMeshes.push_back(heMesh);
+	return error;
+}
+
+void Generate(Atlas *atlas, CharterOptions charterOptions, PackerOptions packerOptions) {
+	xaAssert(atlas);
+	xaAssert(packerOptions.texelArea > 0);
+	// Chart meshes.
+	for (int i = 0; i < (int)atlas->heMeshes.size(); i++) {
+		std::vector<uint32_t> uncharted_materials;
+		atlas->atlas.computeCharts(atlas->heMeshes[i], charterOptions, uncharted_materials);
+	}
+	atlas->atlas.parameterizeCharts();
+	internal::param::AtlasPacker packer(&atlas->atlas);
+	packer.packCharts(packerOptions);
+	//float utilization = return packer.computeAtlasUtilization();
+	atlas->width = packer.getWidth();
+	atlas->height = packer.getHeight();
+	// Build output meshes.
+	atlas->outputMeshes = new OutputMesh *[atlas->heMeshes.size()];
+	for (int i = 0; i < (int)atlas->heMeshes.size(); i++) {
+		const internal::halfedge::Mesh *heMesh = atlas->heMeshes[i];
+		OutputMesh *outputMesh = atlas->outputMeshes[i] = new OutputMesh;
+		const internal::param::MeshCharts *charts = atlas->atlas.meshAt(i);
+		// Vertices.
+		outputMesh->vertexCount = charts->vertexCount();
+		outputMesh->vertexArray = new OutputVertex[outputMesh->vertexCount];
+		for (uint32_t i = 0; i < charts->chartCount(); i++) {
+			const internal::param::Chart *chart = charts->chartAt(i);
+			const uint32_t vertexOffset = charts->vertexCountBeforeChartAt(i);
+			for (uint32_t v = 0; v < chart->vertexCount(); v++) {
+				OutputVertex &output_vertex = outputMesh->vertexArray[vertexOffset + v];
+				output_vertex.xref = chart->mapChartVertexToOriginalVertex(v);
+				internal::Vector2 uv = chart->chartMesh()->vertexAt(v)->tex;
+				output_vertex.uv[0] = uv.x;
+				output_vertex.uv[1] = uv.y;
+			}
+		}
+		// Indices.
+		outputMesh->indexCount = heMesh->faceCount() * 3;
+		outputMesh->indexArray = new uint32_t[outputMesh->indexCount];
+		for (uint32_t f = 0; f < heMesh->faceCount(); f++) {
+			const uint32_t c = charts->faceChartAt(f);
+			const uint32_t i = charts->faceIndexWithinChartAt(f);
+			const uint32_t vertexOffset = charts->vertexCountBeforeChartAt(c);
+			const internal::param::Chart *chart = charts->chartAt(c);
+			xaDebugAssert(i < chart->chartMesh()->faceCount());
+			xaDebugAssert(chart->faceAt(i) == f);
+			const internal::halfedge::Face *face = chart->chartMesh()->faceAt(i);
+			const internal::halfedge::Edge *edge = face->edge;
+			outputMesh->indexArray[3 * f + 0] = vertexOffset + edge->vertex->id;
+			outputMesh->indexArray[3 * f + 1] = vertexOffset + edge->next->vertex->id;
+			outputMesh->indexArray[3 * f + 2] = vertexOffset + edge->next->next->vertex->id;
+		}
+		// Charts.
+		outputMesh->chartCount = charts->chartCount();
+		outputMesh->chartArray = new OutputChart[outputMesh->chartCount];
+		for (uint32_t i = 0; i < charts->chartCount(); i++) {
+			OutputChart *outputChart = &outputMesh->chartArray[i];
+			const internal::param::Chart *chart = charts->chartAt(i);
+			const uint32_t vertexOffset = charts->vertexCountBeforeChartAt(i);
+			const internal::halfedge::Mesh *mesh = chart->chartMesh();
+			outputChart->indexCount = mesh->faceCount() * 3;
+			outputChart->indexArray = new uint32_t[outputChart->indexCount];
+			for (uint32_t j = 0; j < mesh->faceCount(); j++) {
+				const internal::halfedge::Face *face = mesh->faceAt(j);
+				const internal::halfedge::Edge *edge = face->edge;
+				outputChart->indexArray[3 * j + 0] = vertexOffset + edge->vertex->id;
+				outputChart->indexArray[3 * j + 1] = vertexOffset + edge->next->vertex->id;
+				outputChart->indexArray[3 * j + 2] = vertexOffset + edge->next->next->vertex->id;
+			}
+		}
+	}
+}
+
+uint32_t GetWidth(const Atlas *atlas) {
+	xaAssert(atlas);
+	return atlas->width;
+}
+
+uint32_t GetHeight(const Atlas *atlas) {
+	xaAssert(atlas);
+	return atlas->height;
+}
+
+uint32_t GetNumCharts(const Atlas *atlas) {
+	xaAssert(atlas);
+	return atlas->atlas.chartCount();
+}
+
+const OutputMesh *const *GetOutputMeshes(const Atlas *atlas) {
+	xaAssert(atlas);
+	return atlas->outputMeshes;
+}
+
+const char *StringForEnum(AddMeshErrorCode::Enum error) {
+	if (error == AddMeshErrorCode::AlreadyAddedEdge)
+		return "already added edge";
+	if (error == AddMeshErrorCode::DegenerateColocalEdge)
+		return "degenerate colocal edge";
+	if (error == AddMeshErrorCode::DegenerateEdge)
+		return "degenerate edge";
+	if (error == AddMeshErrorCode::DuplicateEdge)
+		return "duplicate edge";
+	if (error == AddMeshErrorCode::IndexOutOfRange)
+		return "index out of range";
+	if (error == AddMeshErrorCode::InvalidIndexCount)
+		return "invalid index count";
+	if (error == AddMeshErrorCode::ZeroAreaFace)
+		return "zero area face";
+	if (error == AddMeshErrorCode::ZeroLengthEdge)
+		return "zero length edge";
+	return "success";
+}
+
+} // namespace xatlas