/*
	This file contains the code listings for the book
	Foundations of Game Engine Development, Volume 2: Rendering
	by Eric Lengyel.

	This code is copyrighted. If you own a copy of FGED2, then you
	have the author's express permission to use this code in whatever
	way you see fit as long as you don't republish it or pass it off
	as your own work.

	If you use any of this code in a commercial product, such as a
	video game, then please give attribution to the author in the credits.
*/


// Listing 5.1

struct ColorRGBA
{
	float		red, green, blue, alpha;

	ColorRGBA() = default;

	ColorRGBA(float r, float g, float b, float a = 1.0F)
	{
		red = r;
		green = g;
		blue = b;
		alpha = a;
	}

	ColorRGBA& operator *=(float s)
	{
		red *= s;
		green *= s;
		blue *= s;
		alpha *= s;
		return (*this);
	}

	ColorRGBA& operator /=(float s)
	{
		s = 1.0F / s;
		red *= s;
		green *= s;
		blue *= s;
		alpha *= s;
		return (*this);
	}

	ColorRGBA& operator +=(const ColorRGBA& c)
	{
		red += c.red;
		green += c.green;
		blue += c.blue;
		alpha += c.alpha;
	}

	ColorRGBA& operator -=(const ColorRGBA& c)
	{
		red -= c.red;
		green -= c.green;
		blue -= c.blue;
		alpha -= c.alpha;
	}

	ColorRGBA& operator *=(const ColorRGBA& c)
	{
		red *= c.red;
		green *= c.green;
		blue *= c.blue;
		alpha *= c.alpha;
	}
};

inline ColorRGBA operator *(const ColorRGBA& c, float s)
{
	return (ColorRGBA(c.red * s, c.green * s, c.blue * s, c.alpha * s));
}

inline ColorRGBA operator /(const ColorRGBA& c, float s)
{
	s = 1.0F / s;
	return (ColorRGBA(c.red * s, c.green * s, c.blue * s, c.alpha * s));
}

inline ColorRGBA operator +(const ColorRGBA& a, const ColorRGBA& b)
{
	return (ColorRGBA(a.red + b.red, a.green + b.green, a.blue + b.blue, a.alpha + b.alpha));
}

inline ColorRGBA operator -(const ColorRGBA& a, const ColorRGBA& b)
{
	return (ColorRGBA(a.red - b.red, a.green - b.green, a.blue - b.blue, a.alpha - b.alpha));
}

inline ColorRGBA operator *(const ColorRGBA& a, const ColorRGBA& b)
{
	return (ColorRGBA(a.red * b.red, a.green * b.green, a.blue * b.blue, a.alpha * b.alpha));
}

// Listing 5.2

struct Triangle
{
	uint16		vertexIndex[3];
};

// Listing 5.3

uniform float4    mvp[4];

float4 TransformVertex(float3 position)
{
	return (float4(dot(position, mvp[0].xyz) + mvp[0].w,
	               dot(position, mvp[1].xyz) + mvp[1].w,
	               dot(position, mvp[2].xyz) + mvp[2].w,
	               dot(position, mvp[3].xyz) + mvp[3].w));
}

// Listing 6.1

Matrix4D MakeFrustumProjection(float fovy, float s, float n, float f)
{
	float g = 1.0F / tan(fovy * 0.5F);
	float k = f / (f - n);

	return (Matrix4D(g / s, 0.0F, 0.0F, 0.0F,
	                  0.0F,  g,   0.0F, 0.0F,
	                  0.0F, 0.0F,  k,  -n * k,
	                  0.0F, 0.0F, 1.0F, 0.0F));
}

// Listing 6.2

Matrix4D MakeInfiniteProjection(float fovy, float s, float n, float e)
{
	float g = 1.0F / tan(fovy * 0.5F);
	e = 1.0F - e;

	return (Matrix4D(g / s, 0.0F, 0.0F, 0.0F,
	                  0.0F,  g,   0.0F, 0.0F,
	                  0.0F, 0.0F,  e,  -n * e,
	                  0.0F, 0.0F, 1.0F, 0.0F));
}

// Listing 6.3

Matrix4D MakeRevFrustumProjection(float fovy, float s, float n, float f)
{
	float g = 1.0F / tan(fovy * 0.5F);
	float k = n / (n - f);

	return (Matrix4D(g / s, 0.0F, 0.0F, 0.0F,
	                  0.0F,  g,   0.0F, 0.0F,
	                  0.0F, 0.0F,  k,  -f * k,
	                  0.0F, 0.0F, 1.0F, 0.0F));
}

// Listing 6.4

Matrix4D MakeRevInfiniteProjection(float fovy, float s, float n, float e)
{
	float g = 1.0F / tan(fovy * 0.5F);

	return (Matrix4D(g / s, 0.0F, 0.0F,    0.0F,
	                  0.0F,  g,   0.0F,    0.0F,
	                  0.0F, 0.0F,  e,   n * (1.0F - e),
	                  0.0F, 0.0F, 1.0F,    0.0F));
}

// Listing 6.5

Matrix4D MakeOrthoProjection(float l, float r, float t, float b, float n, float f)
{
	float w_inv = 1.0F / (r - l);
	float h_inv = 1.0F / (b - t);
	float d_inv = 1.0F / (f - n);

	return (Matrix4D(2.0F * w_inv,   0.0F,      0.0F, -(r + l) * w_inv,
	                    0.0F,     2.0F * h_inv, 0.0F, -(b + t) * h_inv,
	                    0.0F,        0.0F,     d_inv,   -n * d_inv,
	                    0.0F,        0.0F,      0.0F,      1.0F));
}

// Listing 6.6

void ModifyProjectionNearPlane(Matrix4D& P, const Plane& k)
{
	Vector4D vcamera((sgn(k.x) - P(0,2)) / P(0,0),
	                 (sgn(k.y) - P(1,2)) / P(1,1),
	                 1.0F, (1.0F - P(2,2)) / P(2,3));

	float m = 1.0F / Dot(k, vcamera);
	P(2,0) = m * k.x;
	P(2,1) = m * k.y;
	P(2,2) = m * k.z;
	P(2,3) = m * k.w;
}

void ModifyRevProjectionNearPlane(Matrix4D& R, const Plane& k)
{
	Vector4D vcamera((sgn(k.x) - R(0,2)) / R(0,0),
	                 (sgn(k.y) - R(1,2)) / R(1,1),
	                 1.0F, R(2,2) / R(2,3));

	float m = -1.0F / Dot(k, vcamera);
	R(2,0) = m * k.x;
	R(2,1) = m * k.y;
	R(2,2) = m * k.z + 1.0F;
	R(2,3) = m * k.w;
}

// Listing 7.1

uniform float3	diffuseColor;		// (rho / pi) * C_diffuse
uniform float3	ambientColor;		// pi * C_ambient
uniform float3	lightColor;			// C_illum

float3 CalculateDiffuseReflection(float3 n, float3 l)
{
	float3 directColor = lightColor * saturate(dot(n, l));
	return ((ambientColor + directColor) * diffuseColor);
}

// Listing 7.2

uniform float3	specularColor;		// C_specular
uniform float3	lightColor;			// C_illum

float3 CalculateSpecularReflection(float3 n, float3 h, float alpha, float nl)
{
	float highlight = pow(saturate(dot(n, h)), alpha) * float(nl > 0.0);
	return (lightColor * specularColor * highlight);
}

// Listing 7.3

uniform TextureCube		environmentMap;
uniform float3			Menv[3];

float4 SampleEnvironmentMap(float3 n, float3 v)
{
	float3 r = n * (2.0 * dot(n, v)) - v;
	float3 texcoord = float3(dot(Menv[0], r), dot(Menv[1], r), dot(Menv[2], r));
	return (texture(environmentMap, texcoord));
}

// Listing 7.4

void CalculateTangents(int32 triangleCount, const Triangle *triangleArray, int32 vertexCount, const Point3D *vertexArray, const Vector3D *normalArray, const Point2D *texcoordArray, Vector4D *tangentArray)
{
	// Allocate temporary storage for tangents and bitangents and initialize to zeros.
	Vector3D *tangent = new Vector3D[vertexCount * 2];
	Vector3D *bitangent = tangent + vertexCount;
	for (int32 i = 0; i < vertexCount; i++)
	{
		tangent[i].Set(0.0F, 0.0F, 0.0F);
		bitangent[i].Set(0.0F, 0.0F, 0.0F);
	}

	// Calculate tangent and bitangent for each triangle and add to all three vertices.
	for (int32 k = 0; k < triangleCount; k++)
	{
		int32 i0 = triangleArray[k].index[0];
		int32 i1 = triangleArray[k].index[1];
		int32 i2 = triangleArray[k].index[2];
		const Point3D& p0 = vertexArray[i0];
		const Point3D& p1 = vertexArray[i1];
		const Point3D& p2 = vertexArray[i2];
		const Point2D& w0 = texcoordArray[i0];
		const Point2D& w1 = texcoordArray[i1];
		const Point2D& w2 = texcoordArray[i2];

		Vector3D e1 = p1 - p0, e2 = p2 - p0;
		float x1 = w1.x - w0.x, x2 = w2.x - w0.x;
		float y1 = w1.y - w0.y, y2 = w2.y - w0.y;

		float r = 1.0F / (x1 * y2 - x2 * y1);
		Vector3D t = (e1 * y2 - e2 * y1) * r;
		Vector3D b = (e2 * x1 - e1 * x2) * r;

		tangent[i0] += t;
		tangent[i1] += t;
		tangent[i2] += t;
		bitangent[i0] += b;
		bitangent[i1] += b;
		bitangent[i2] += b;
	}

	// Orthonormalize each tangent and calculate the handedness.
	for (int32 i = 0; i < vertexCount; i++)
	{
		const Vector3D& t = tangent[i];
		const Vector3D& b = bitangent[i];
		const Vector3D& n = normalArray[i];
		tangentArray[i].xyz() = Normalize(Reject(t, n));
		tangentArray[i].w = (Dot(Cross(t, b), n) > 0.0F) ? 1.0F : -1.0F;
	}

	delete[] tangent;
}

// Listing 7.5

void ConstructNormalMap(const float *heightMap, Vector3D *normalMap, int32 width, int32 height)
{
	for (int32 y = 0; y < height; y++)
	{
		int32 ym1 = (y - 1) & (height - 1);
		int32 yp1 = (y + 1) & (height - 1);

		const float *centerRow = heightMap + y * width;
		const float *upperRow = heightMap + ym1 * width;
		const float *lowerRow = heightMap + yp1 * width;

		for (int32 x = 0; x < width; x++)
		{
			int32 xm1 = (x - 1) & (width - 1);
			int32 xp1 = (x + 1) & (width - 1);

			// Calculate slopes.
			float dx = (centerRow[xp1] - centerRow[xm1]) * 0.5F;
			float dy = (lowerRow[x] - upperRow[x]) * 0.5F;

			// Normalize and clamp.
			float nz = 1.0F / sqrt(dx * dx + dy * dy + 1.0F);
			float nx = fmin(fmax(-dx * nz, -1.0F), 1.0F);
			float ny = fmin(fmax(-dy * nz, -1.0F), 1.0F);
			normalMap[x].Set(nx, ny, nz);
		}

		normalMap += width;
	}
}

// Listing 7.6

uniform Texture2D		normalMap;

float3 FetchNormalVector(float2 texcoord)
{
	float2 m = texture(normalMap, texcoord).xy;
	return (float3(m, sqrt(1.0 - m.x * m.x - m.y * m.y)));
}

// Listing 7.7

uniform float3 cameraPosition;		// Object-space camera position.
uniform float3 lightPosition;		// Object-space light position.

void CalculateTangentSpaceVL(float3 position, float3 normal, float4 tangent, out float3 vtan, out float3 ltan)
{
	float3 bitangent = cross(normal, tangent.xyz) * tangent.w;
	float3 v = cameraPosition - position;
	float3 l = lightPosition - position;
	vtan = float3(dot(tangent, v), dot(bitangent, v), dot(normal, v));
	ltan = float3(dot(tangent, l), dot(bitangent, l), dot(normal, l));
}

// Listing 7.8

uniform Texture2D		normalMap;

float3 FetchObjectNormalVector(float2 texcoord, float3 normal, float3 tangent, float sigma)
{
	float3 m = FetchNormalVector(texcoord);
	float3 n = normalize(normal);
	float3 t = normalize(tangent - n * dot(tangent, n));
	float3 b = cross(n, t) * sigma;
	return (t * m.x + b * m.y + n * m.z);
}

// Listing 7.9

uniform Texture2D		parallaxMap;

float2 ApplyParallaxOffset(float2 texcoord, float3 vdir, float2 scale)
{
	float2 pdir = vdir.xy * scale;
	for (int i = 0; i < 4; i++)
	{
		// Fetch n.z * h from the parallax map.
		float parallax = texture(parallaxMap, texcoord).x;
		texcoord += pdir * parallax;
	}

	return (texcoord);
}

// Listing 7.10

void ConstructHorizonMap(const float *heightMap, ColorRGBA *horizonMap, float *ambientMap, float ambientPower, int32 width, int32 height)
{
	constexpr int kAngleCount = 32;		// Must be at least 16 and a power of 2.
	constexpr float kAngleIndex = float(kAngleCount) / two_pi;
	constexpr int kHorizonRadius = 16;

	for (int32 y = 0; y < height; y++)
	{
		const float *centerRow = heightMap + y * width;
		for (int32 x = 0; x < width; x++)
		{
			// Get central height. Initialize max squared tangent array to all zeros.
			float h0 = centerRow[x];
			float maxTan2[kAngleCount] = {};

			// Search neighborhood for larger heights.
			for (int32 j = -kHorizonRadius + 1; j < kHorizonRadius; j++)
			{
				const float *row = heightMap + ((y + j) & (height - 1)) * width;
				for (int32 i = -kHorizonRadius + 1; i < kHorizonRadius; i++)
				{
					int32 r2 = i * i + j * j;
					if ((r2 < kHorizonRadius * kHorizonRadius) && (r2 != 0))
					{
						float dh = row[(x + i) & (width - 1)] - h0;
						if (dh > 0.0F)
						{
							// Larger height found. Apply to array entries.
							float direction = atan2(float(j), float(i));
							float delta = atan(0.7071F / sqrt(float(r2)));
							int32 minIndex = int32(floor((direction - delta) * kAngleIndex));
							int32 maxIndex = int32(ceil((direction + delta) * kAngleIndex));

							// Calculate squared tangent with Equation (7.53).
							float t = dh * dh / float(r2);
							for (int32 n = minIndex; n <= maxIndex; n++)
							{
								int32 m = n & (kAngleCount - 1);
								maxTan2[m] = fmax(maxTan2[m], t);
							}
						}
					}
				}
			}

			// Generate eight channels of horizon map.
			ColorRGBA *layerData = horizonMap;
			for (int32 layer = 0; layer < 2; layer++)
			{
				ColorRGBA color(0.0F, 0.0F, 0.0F, 0.0F);
				int32 firstIndex = kAngleCount / 16 + layer * (kAngleCount / 2);
				int32 lastIndex = firstIndex + kAngleCount / 8;

				for (int32 index = firstIndex; index <= lastIndex; index++)
				{
					float tr = maxTan2[(index - kAngleCount / 8) & (kAngleCount - 1)];
					float tg = maxTan2[index];
					float tb = maxTan2[index + kAngleCount / 8];
					float ta = maxTan2[(index + kAngleCount / 4) & (kAngleCount - 1)];

					color.red += sqrt(tr / (tr + 1.0F));
					color.green += sqrt(tg / (tg + 1.0F));
					color.blue += sqrt(tb / (tb + 1.0F));
					color.alpha += sqrt(ta / (ta + 1.0F));
				}

				layerData[x] = color / float(kAngleCount / 8 + 1);
				layerData += width * height;
			}

			// Generate ambient light factor.
			float sum = 0.0F;
			for (int32 k = 0; k < kAngleCount; k++) sum += 1.0F / sqrt(maxTan2[k] + 1.0F);
			ambientMap[x] = pow(sum * (1.0F / float(kAngleCount)), ambientPower);
		}

		horizonMap += width;
		ambientMap += width;
	}
}

// Listing 7.11

void GenerateHorizonCube(ColorRGBA *texel)
{
	for (int face = 0; face < 6; face++)
	{
		for (float y = -0.9375F; y < 1.0F; y += 0.125F)
		{
			for (float x = -0.9375F; x < 1.0F; x += 0.125F)
			{
				Vector2D	v;

				float r = 1.0F / sqrt(1.0F + x * x + y * y);
				switch (face)
				{
					case 0:
						v.Set(r, -y * r);
						break;
					case 1:
						v.Set(-r, -y * r);
						break;
					case 2:
						v.Set(x * r, r);
						break;
					case 3:
						v.Set(x * r, -r);
						break;
					case 4:
						v.Set(x * r, -y * r);
						break;
					case 5:
						v.Set(-x * r, -y * r);
						break;
				}

				float t = atan2(v.y, v.x) / pi_over_4;

				float red = 0.0F;
				float green = 0.0F;
				float blue = 0.0F;
				float alpha = 0.0F;

				if (t < -3.0F)
				{
					red = t + 3.0F;
					green = -4.0F - t;
				}
				else if (t < -2.0F)
				{
					green = t + 2.0F;
					blue = -3.0F - t;
				}
				else if (t < -1.0F)
				{
					blue = t + 1.0F;
					alpha = -2.0F - t;
				}
				else if (t <  0.0F)
				{
					alpha = t;
					red = t + 1.0F;
				}
				else if (t <  1.0F)
				{
					red = 1.0F - t;
					green = t;
				}
				else if (t <  2.0F)
				{
					green = 2.0F - t;
					blue = t - 1.0F;
				}
				else if (t <  3.0F)
				{
					blue = 3.0F - t;
					alpha = t - 2.0F;
				}
				else
				{
					alpha = 4.0F - t;
					red = 3.0F - t;
				}

				texel->Set(red, green, blue, alpha);
				texel++;
			}
		}
	}
}

// Listing 7.12

uniform Texture2DArray   horizonMap;
uniform TextureCube      weightCube;

float ApplyHorizonMap(float2 texcoord, float3 ldir)
{
	const float kShadowHardness = 8.0;

	// Read horizon channel factors from cube map.
	float4 weights = texture(weightCube, ldir);

	// Extract positive and negative weights for horizon map layers 0 and 1.
	float4 w0 = saturate(weights);
	float4 w1 = saturate(-weights);

	// Sample the horizon map and multiply by the weights for each layer.
	float s0 = dot(texture(horizonMap, float3(texcoord, 0.0)), w0);
	float s1 = dot(texture(horizonMap, float3(texcoord, 1.0)), w1);

	// Return lighting factor calculated with Equation (7.58).
	return (saturate((ldir.z - (s0 + s1)) * kShadowHardness + 1.0));
}

// Listing 8.1

uniform float2	attenConst;		// (r0^2 / (rmax2 - r0^2), rmax2)

float CalculateInverseSquareAttenuation(float3 p, float3 l)
{
	float3 ldir = l - p;
	float r2 = dot(ldir, ldir);
	return (saturate(attenConst.x * (attenConst.y / r2 - 1.0)));
}

// Listing 8.2

uniform float3	attenConst;		// (-k2/rmax2, 1/(1 - exp(-k2)), exp(-k2)/(1 - exp(-k2)))

float CalculateExponentialAttenuation(float3 p, float3 l)
{
	float3 ldir = l - p;
	float r2 = dot(ldir, ldir);
	return (saturate(exp(r2 * attenConst.x) * attenConst.y - attenConst.z));
}

// Listing 8.3

uniform float2	attenConst;		// (1 / rmax2, 2 / rmax)

float CalculateSmoothAttenuation(float3 p, float3 l)
{
	float3 ldir = l - p;
	float r2 = dot(ldir, ldir);
	return (saturate(r2 * attenConst.x * (sqrt(r2) * attenConst.y - 3.0) + 1.0));
}

// Listing 8.4

bool CalculateScissorRect(const Point3D& center, float radius, float projectionDistance, float aspectRatio, Point2D *rectMin, Point2D *rectMax)
{
	// Initialize the rectangle to the full viewport.
	float xmin = -1.0F;
	float xmax = 1.0F;
	float ymin = -1.0F;
	float ymax = 1.0F;

	// Make sure the sphere isn't behind the camera.
	if (center.z > -radius)
	{
		float lx2 = center.x * center.x;
		float ly2 = center.y * center.y;
		float lz2 = center.z * center.z;
		float r2 = radius * radius;

		// The factors mx and my project points into the viewport.
		float my = projectionDistance / center.z;
		float mx = my / aspectRatio;

		float a = lz2 - r2;
		if (fabs(a) > FLT_MIN)		// Quadratic case when a != 0.
		{
			float f = lz2 * r2;
			float a_inv = 1.0F / a;
			float d = f * (a + lx2);
			if (d > 0.0F)				// Discriminant positive for x extents.
			{
				d = sqrt(d);
				float b = r2 * center.x;
				float x1 = (center.x + (b - d) * a_inv) * mx;
				float x2 = (center.x + (b + d) * a_inv) * mx;

				if (a > 0.0F)			// Both tangencies valid.
				{
					xmin = fmax(xmin, x1); xmax = fmin(xmax, x2);
				}
				else 					// One tangency valid.
				{
					if (center.x > 0.0F) xmin = fmax(xmin, fmax(x1, x2));
					else xmax = fmin(xmax, fmin(x1, x2));
				}
			}

			d = f * (a + ly2);
			if (d > 0.0F)				// Discriminant positive for y extents.
			{
				d = sqrt(d);
				float b = r2 * center.y;
				float y1 = (center.y + (b - d) * a_inv) * my;
				float y2 = (center.y + (b + d) * a_inv) * my;

				if (a > 0.0F)			// Both tangencies valid.
				{
					ymin = fmax(ymin, y1);
					ymax = fmin(ymax, y2);
				}
				else					// One tangency valid.
				{
					if (center.y > 0.0F) ymin = fmax(ymin, fmax(y1, y2));
					else ymax = fmin(ymax, fmin(y1, y2));
				}
			}
		}
		else		// Linear case when a == 0.
		{
			if (fabs(center.x) > FLT_MIN)		// Tangency valid for x.
			{
				float x = (center.x - (lx2 + r2) * 0.5F / center.x) * mx;
				if (center.x > 0.0F) xmin = fmax(xmin, x);
				else xmax = fmin(xmax, x);	
			}

			if (fabs(center.y) > FLT_MIN)		// Tangency valid for y.
			{
				float y = (center.y - (ly2 + r2) * 0.5F / center.y) * my;
				if (center.y > 0.0F) ymin = fmax(ymin, y);
				else ymax = fmin(ymax, y);	
			}
		}
	}

	rectMin->Set(xmin, ymin);
	rectMax->Set(xmax, ymax);
	return ((xmin < xmax) && (ymin < ymax));
}

// Listing 8.5

bool CalculateDepthBounds(const Matrix4D& P, float n, const Point3D& center, float radius, float *minDepth, float *maxDepth)
{
	float zmin = fmin(P(2,2) + P(2,3) / fmax(center.z - radius, n), 1.0F);
	float zmax = fmin(P(2,2) + P(2,3) / fmax(center.z + radius, n), 1.0F);
	*minDepth = zmin;
	*maxDepth = zmax;
	return (zmin < zmax);
}

bool CalculateRevDepthBounds(const Matrix4D& R, float n, const Point3D& center, float radius, float *minDepth, float *maxDepth)
{
	float zmin = fmax(R(2,2) + R(2,3) / fmax(center.z + radius, n), 0.0F);
	float zmax = fmax(R(2,2) + R(2,3) / fmax(center.z - radius, n), 0.0F);
	*minDepth = zmin;
	*maxDepth = zmax;
	return (zmin < zmax);
}

// Listing 8.6

void CalculateObliqueDepthBounds(const Matrix4D& P, const Point3D& center, float radius, float *minDepth, float *maxDepth)
{
	float zmin = 0.0F;
	float zmax = 1.0F;

	float kz = P(2,2);
	float kxy2 = P(2,0) * P(2,0) + P(2,1) * P(2,1);
	float kxyz2 = kxy2 + kz * kz;
	float klxyw = P(2,0) * center.x + P(2,1) * center.y + P(2,3);
	float kl = klxyw + kz * center.z;    // k dot l
	float r2 = radius * radius;

	float a = center.z * center.z - r2;
	float b = kl * center.z - r2 * kz;
	float c = kl * kl - r2 * kxyz2;
	if (fabs(a) > FLT_MIN)
	{
		float d = (klxyw * klxyw + a * kxy2) * r2;
		if (d > FLT_MIN)
		{
			// Bounding sphere does not intersect red cylinder in figure.
			float a_inv = 1.0F / a;
			float f = sqrt(d) * a_inv;
			if (a > 0.0F)
			{
				// Bounding sphere is completely in front of x-y plane.
				// Both roots correspond to valid depth bounds.
				zmin = fmin(fmax(b * a_inv - f, 0.0F), 1.0F);
				zmax = fmin(fmax(b * a_inv + f, 0.0F), 1.0F);
			}
			else if (c > 0.0F)
			{
				// Bounding sphere is completely in front of near plane.
				// Roots must have different signs. Positive root is zmin.
				zmin = fmax(b * a_inv - f, 0.0F);
			}
			else
			{
				float hl = kl * kz - kxyz2 * center.z;    // h dot l
				if (b < -FLT_MIN)
				{
					// Light position is inside striped parallelogram.
					if (hl < 0.0F) zmax = b * a_inv + f;
					else zmin = b * a_inv - f;
				}
				else if (b > FLT_MIN)
				{
					// If light position is between points C and D, it's not visible.
					if (hl > 0.0F) zmax = 0.0F;
				}
				else
				{
					// If light position is at point C, it's not visible.
					if (kl < 0.0F) zmax = 0.0F;
				}
			}
		}
	}
	else if (center.z > 0.0F)  // The light position is on upper boundary of green band.
	{
		if (fabs(b) > FLT_MIN)
		{
			float zdev = c / b * 0.5F;
			if (b < 0.0F) zmax = fmin(zdev, 1.0F);
			else zmin = fmax(zdev, 0.0F);
		}
	}
	else
	{
		zmax = 0.0F;	// Light position is on lower boundary of green band.
	}

	*minDepth = zmin;
	*maxDepth = zmax;
	return (zmin < zmax);
}

// Listing 8.7

uniform Texture2DShadow		shadowTexture;
uniform float				shadowOffset;	

float CalculateSpotShadow(float4 shadowCoord)
{
	// Project the texture coordinates and fetch the center sample.
	float3 p = shadowCoord.xyz / shadowCoord.w;
	float light = texture(shadowTexture, p);

	// Fetch four more samples at diagonal offsets.
	p.xy -= shadowOffset;
	light += texture(shadowTexture, p.xy, p.z);
	p.x += shadowOffset * 2.0;
	light += texture(shadowTexture, p.xy, p.z);
	p.y += shadowOffset * 2.0;
	light += texture(shadowTexture, p.xy, p.z);
	p.x -= shadowOffset * 2.0;
	light += texture(shadowTexture, p.xy, p.z);

	return (light * 0.2);	// Return average value.
}

// Listing 8.8

uniform TextureCubeShadow		shadowTexture;
uniform float					shadowOffset;		// 2.0 * delta
uniform float2					depthTransform;		// (m22, m23)

float CalculatePointShadow(float3 lightCoord)
{
	float3 absq = abs(lightCoord);
	float mxy = max(absq.x, absq.y);
	float m = max(mxy, absq.z);			// Value of largest component.

	// Calculate offset vectors.
	float offset = shadowOffset * m;
	float dxy = (mxy > absq.z) ? offset : 0.0;
	float dx = (absq.x > absq.y) ? dxy : 0.0;
	float2 oxy = float2(offset - dx, dx);
	float2 oyz = float2(offset - dxy, dxy);

	float3 limit = float3(m, m, m);
	limit.xy -= oxy * (1.0 / 1024.0);	// Epsilon = 1/1024.
	limit.yz -= oyz * (1.0 / 1024.0);

	// Calculate projected depth and fetch the center sample.
	float depth = depthTransform.x + depthTransform.y / m;
	float light = texture(shadowTexture, lightCoord, depth);

	// Fetch four more samples at diagonal offsets.
	lightCoord.xy -= oxy;
	lightCoord.yz -= oyz;
	light += texture(shadowTexture, clamp(lightCoord, -limit, limit), depth);
	lightCoord.xy += oxy * 2.0;
	light += texture(shadowTexture, clamp(lightCoord, -limit, limit), depth);
	lightCoord.yz += oyz * 2.0;
	light += texture(shadowTexture, clamp(lightCoord, -limit, limit), depth);
	lightCoord.xy -= oxy * 2.0;
	light += texture(shadowTexture, clamp(lightCoord, -limit, limit), depth);

	return (light * 0.2);	// Return average value.
}

// Listing 8.9

uniform Texture2DArrayShadow	shadowTexture;
uniform float4					shadowOffset[2];
uniform float3					cascadeScale[3];
uniform float3					cascadeOffset[3];

float CalculateInfiniteShadow(float3 cascadeCoord0, float3 cascadeBlend)
{
	float3		p1, p2;

	// Apply scales and offsets to get texcoords in all four cascades.
	float3 cascadeCoord1 = cascadeCoord0 * cascadeScale[0] + cascadeOffset[0];
	float3 cascadeCoord2 = cascadeCoord0 * cascadeScale[1] + cascadeOffset[1];
	float3 cascadeCoord3 = cascadeCoord0 * cascadeScale[2] + cascadeOffset[2];
	
	// Calculate layer indices i and j.
	bool beyondCascade2 = (cascadeBlend.y >= 0.0);
	bool beyondCascade3 = (cascadeBlend.z >= 0.0);
	p1.z = float(beyondCascade2) * 2.0;
	p2.z = float(beyondCascade3) * 2.0 + 1.0;

	// Select texture coordinates.
	float2 shadowCoord1 = (beyondCascade2) ? cascadeCoord2.xy : cascadeCoord0.xy;
	float2 shadowCoord2 = (beyondCascade3) ? cascadeCoord3.xy : cascadeCoord1.xy;
	float depth1 = (beyondCascade2) ? cascadeCoord2.z : cascadeCoord0.z;
	float depth2 = (beyondCascade3) ? saturate(cascadeCoord3.z) : cascadeCoord1.z;

	// Calculate blend weight w.
	float3 blend = saturate(cascadeBlend);
	float weight = (beyondCascade2) ? blend.y - blend.z : 1.0 - blend.x;

	// Fetch four samples from the first cascade.
	p1.xy = shadowCoord1 + shadowOffset[0].xy;
	float light1 = texture(shadowTexture, p1, depth1);
	p1.xy = shadowCoord1 + shadowOffset[0].zw;
	light1 += texture(shadowTexture, p1, depth1);
	p1.xy = shadowCoord1 + shadowOffset[1].xy;
	light1 += texture(shadowTexture, p1, depth1);
	p1.xy = shadowCoord1 + shadowOffset[1].zw;
	light1 += texture(shadowTexture, p1, depth1);

	// Fetch four samples from the second cascade.
	p2.xy = shadowCoord2 + shadowOffset[0].xy;
	float light2 = texture(shadowTexture, p2, depth2);
	p2.xy = shadowCoord2 + shadowOffset[0].zw;
	light2 += texture(shadowTexture, p2, depth2);
	p2.xy = shadowCoord2 + shadowOffset[1].xy;
	light2 += texture(shadowTexture, p2, depth2);
	p2.xy = shadowCoord2 + shadowOffset[1].zw;
	light2 += texture(shadowTexture, p2, depth2);

	// Return blended average value.
	return (lerp(light2, light1, weight) * 0.25);
}

// Listing 8.10

struct Edge
{
	uint16		vertexIndex[2];
	uint16		faceIndex[2];
};

// Listing 8.11

int32 BuildEdgeArray(int32 vertexCount, int32 triangleCount, const Triangle *triangleArray, Edge *edgeArray)
{
	// Initialize all edge lists to empty.
	uint16 *firstEdge = new uint16[vertexCount + triangleCount * 2];
	uint16 *nextEdge = firstEdge + vertexCount;
	for (int32 k = 0; k < vertexCount; k++) firstEdge[k] = 0xFFFF;

	int32 edgeCount = 0;
	const Triangle *triangle = triangleArray;

	// Identify all edges that have increasing vertex indices in CCW direction.
	for (int32 k = 0; k < triangleCount; k++)
	{
		uint16 i1 = triangle[k].index[2];
		for (int32 v = 0; v < 3; v++)
		{
			uint16 i2 = triangle[k].index[v];
			if (i1 < i2)
			{
				Edge *edge = &edgeArray[edgeCount];
				edge->vertexIndex[0] = i1;
				edge->vertexIndex[1] = i2;
				edge->faceIndex[0] = uint16(k);
				edge->faceIndex[1] = uint16(k);

				// Add the edge to the front of the list for the first vertex.
				nextEdge[edgeCount] = firstEdge[i1];
				firstEdge[i1] = edgeCount++;
			}

			i1 = i2;
		}
	}

	// Match all edges to the triangles for which they are wound clockwise.
	for (int32 k = 0; k < triangleCount; k++)
	{
		uint16 i1 = triangle[k].index[2];
		for (int32 v = 0; v < 3; v++)
		{
			uint16 i2 = triangle[k].index[v];
			if (i1 > i2)
			{
				for (uint16 e = firstEdge[i2]; e != 0xFFFF; e = nextEdge[e])
				{
					Edge *edge = &edgeArray[e];
					if ((edge->vertexIndex[1] == i1) && (edge->faceIndex[0] == edge->faceIndex[1]))
					{
						edge->faceIndex[1] = uint16(k);
						break;
					}
				}
			}

			i1 = i2;
		}
	}

	delete[] firstEdge;
	return (edgeCount);
}

// Listing 8.12

uniform float4		mvp[4];
uniform float4		lightPosition;

float4 ExtrudeSilhouette(float4 position)
{
	float3 v = lerp(position.xyz * lightPosition.w - lightPosition.xyz, position.xyz, position.w);

	return (float4(dot(v, mvp[0].xyz) + position.w * mvp[0].w,
	               dot(v, mvp[1].xyz) + position.w * mvp[1].w,
	               dot(v, mvp[2].xyz) + position.w * mvp[2].w,
	               dot(v, mvp[3].xyz) + position.w * mvp[3].w));
}

// Listing 8.13

uniform float4		mvp[4];
uniform float4		lightPosition;

float4 ExtrudeDarkCap(float3 position)
{
	float3 v = position - lightPosition.xyz;
	return (float4(dot(v, mvp[0].xyz),
	               dot(v, mvp[1].xyz),
	               dot(v, mvp[2].xyz),
				   dot(v, mvp[3].xyz)));
}

// Listing 8.14

uniform float3		fogColor;
uniform float		fogDensity;

float3 ApplyFog(float3 shadedColor, float3 v)
{
	float f = exp(-fogDensity * length(v));
	return (lerp(fogColor, shadedColor, f));
}

// Listing 8.15

uniform float3		fogColor;
uniform float		fogDensity;

float3 ApplyHalfspaceFog(float3 shadedColor, float3 v, float fv, float u1, float u2)
{
	const float kFogEpsilon = 0.0001;

	float x = min(u2, 0.0);
	float tau = 0.5 * fogDensity * length(v) * (u1 - x * x / (abs(fv) + kFogEpsilon));
	return (lerp(fogColor, shadedColor, exp(tau)));
}

// Listing 9.1

int32 ClipPolygon(int32 vertexCount, const Point3D *vertex, const Plane& plane, float *location, Point3D *result)
{
	const float kPolygonEpsilon = 0.001F;
	int32 positiveCount = 0, negativeCount = 0;

	// Calculate the signed distance to plane for all vertices.
	for (int32 a = 0; a < vertexCount; a++)
	{
		float d = Dot(plane, vertex[a]);
		location[a] = d;

		if (d > kPolygonEpsilon) positiveCount++;
		else if (d < -kPolygonEpsilon) negativeCount++;
	}

	if (negativeCount == 0)
	{
		// No vertices on negative side of plane. Copy original polygon to result.
		for (int32 a = 0; a < vertexCount; a++) result[a] = vertex[a];
		return (vertexCount);
	}
	else if (positiveCount == 0)
	{
		// No vertices on positive side of plane.
		return (0);
	}

	// Loop through all edges, starting with edge from last vertex to first vertex.
	int32 resultCount = 0;
	const Point3D *p1 = &vertex[vertexCount - 1];
	float d1 = location[vertexCount - 1];
	for (int32 index = 0; index < vertexCount; index++)
	{
		const Point3D *p2 = &vertex[index];
		float d2 = location[index];
		if (d2 < -kPolygonEpsilon)
		{
			// Current vertex is on negative side of plane.
			if (d1 > kPolygonEpsilon)
			{
				// Preceding vertex is on positive side of plane.
				float t = d1 / (d1 - d2);
				result[resultCount++] = *p1 * (1.0F - t) + *p2 * t;
			}
		}
		else
		{
			// Current vertex is on positive side of plane or in plane.
			if ((d2 > kPolygonEpsilon) && (d1 < -kPolygonEpsilon))
			{
				// Current vertex on positive side, and preceding vertex on negative side.
				float t = d2 / (d2 - d1);
				result[resultCount++] = *p2 * (1.0F - t) + *p1 * t;
			}

			result[resultCount++] = *p2;
		}

		p1 = p2;
		d1 = d2;
	}

	return (resultCount);
}

// Listing 9.2

constexpr int32 kMaxPolyhedronVertexCount   = 28;
constexpr int32 kMaxPolyhedronFaceCount     = 16;
constexpr int32 kMaxPolyhedronEdgeCount     = (kMaxPolyhedronFaceCount - 2) * 3;
constexpr int32 kMaxPolyhedronFaceEdgeCount = kMaxPolyhedronFaceCount - 1;

struct Edge
{
	uint8			vertexIndex[2];
	uint8			faceIndex[2];
};

struct Face
{
	uint8			edgeCount;
	uint8			edgeIndex[kMaxPolyhedronFaceEdgeCount];
};

struct Polyhedron
{
	uint8			vertexCount;
	uint8			edgeCount;
	uint8			faceCount;
	Point3D			vertex[kMaxPolyhedronVertexCount];
	Edge			edge[kMaxPolyhedronEdgeCount];
	Face			face[kMaxPolyhedronFaceCount];
	Plane			plane[kMaxPolyhedronFaceCount];
};

// Listings 9.3, 9.4, 9.5, and 9.6

bool ClipPolyhedron(const Polyhedron *polyhedron, const Plane& plane, Polyhedron *result)
{
	float		vertexLocation[kMaxPolyhedronVertexCount];
	int8		vertexCode[kMaxPolyhedronVertexCount];
	int8		edgeCode[kMaxPolyhedronEdgeCount];
	uint8		vertexRemap[kMaxPolyhedronVertexCount];
	uint8		edgeRemap[kMaxPolyhedronEdgeCount];
	uint8		faceRemap[kMaxPolyhedronFaceCount];
	uint8		planeEdgeTable[kMaxPolyhedronFaceEdgeCount];

	const float kPolyhedronEpsilon = 0.001F;
	int32 minCode = 6;
	int32 maxCode = 0;

	// Classify vertices.
	uint32 vertexCount = polyhedron->vertexCount;
	for (uint32 a = 0; a < vertexCount; a++)
	{
		vertexRemap[a] = 0xFF;
		float d = Dot(plane, polyhedron->vertex[a]);
		vertexLocation[a] = d;

		int8 code = (d > -kPolyhedronEpsilon) + (d > kPolyhedronEpsilon) * 2;
		minCode = min(minCode, code);
		maxCode = max(maxCode, code);
		vertexCode[a] = code;
	}

	if (minCode != 0)
	{
		*result = *polyhedron;			// No vertices on negative side of clip plane.
		return (true);
	}

	if (maxCode <= 1) return (false);	// No vertices on positive side of clip plane.

	// Classify edges.
	uint32 edgeCount = polyhedron->edgeCount;
	for (uint32 a = 0; a < edgeCount; a++)
	{
		edgeRemap[a] = 0xFF;
		const Edge *edge = &polyhedron->edge[a];
		edgeCode[a] = int8(vertexCode[edge->vertexIndex[0]] + vertexCode[edge->vertexIndex[1]]);
	}

	// Determine which faces will be in result.
	uint32 resultFaceCount = 0;
	uint32 faceCount = polyhedron->faceCount;
	for (uint32 a = 0; a < faceCount; a++)
	{
		faceRemap[a] = 0xFF;
		const Face *face = &polyhedron->face[a];
		uint32 faceEdgeCount = face->edgeCount;
		for (uint32 b = 0; b < faceEdgeCount; b++)
		{
			if (edgeCode[face->edgeIndex[b]] >= 3)
			{
				// Face has a vertex on the positive side of the plane.
				result->plane[resultFaceCount] = polyhedron->plane[a];
				faceRemap[a] = uint8(resultFaceCount++);
				break;
			}
		}
	}

	uint32 resultVertexCount = 0, resultEdgeCount = 0;
	for (uint32 a = 0; a < edgeCount; a++)
	{
		if (edgeCode[a] >= 2)
		{
			// The edge is not completely clipped away.
			const Edge *edge = &polyhedron->edge[a];
			Edge *resultEdge = &result->edge[resultEdgeCount];
			edgeRemap[a] = uint8(resultEdgeCount++);

			resultEdge->faceIndex[0] = faceRemap[edge->faceIndex[0]];
			resultEdge->faceIndex[1] = faceRemap[edge->faceIndex[1]];

			// Loop over both vertices of edge.
			for (int32 i = 0; i < 2; i++)
			{
				uint8 vertexIndex = edge->vertexIndex[i];
				if (vertexCode[vertexIndex] != 0)
				{
					// This vertex on positive side of plane or in plane.
					uint8 remappedVertexIndex = vertexRemap[vertexIndex];
					if (remappedVertexIndex == 0xFF)
					{
						remappedVertexIndex = resultVertexCount++;
						vertexRemap[vertexIndex] = remappedVertexIndex;
						result->vertex[remappedVertexIndex] = polyhedron->vertex[vertexIndex];
					}

					resultEdge->vertexIndex[i] = remappedVertexIndex;
				}
				else
				{
					// This vertex on negative side, and other vertex on positive side.
					uint8 otherVertexIndex = edge->vertexIndex[1 - i];
					const Point3D& p1 = polyhedron->vertex[vertexIndex];
					const Point3D& p2 = polyhedron->vertex[otherVertexIndex];
					float d1 = vertexLocation[vertexIndex];
					float d2 = vertexLocation[otherVertexIndex];
					float t = d2 / (d2 - d1);
					result->vertex[resultVertexCount] = p2 * (1.0F - t) + p1 * t;
					resultEdge->vertexIndex[i] = uint8(resultVertexCount++);
				}
			}
		}
	}

	uint32 planeEdgeCount = 0;
	for (uint32 a = 0; a < faceCount; a++)
	{
		uint8 remappedFaceIndex = faceRemap[a];
		if (remappedFaceIndex != 0xFF)
		{
			// The face is not completely clipped away.
			Edge *newEdge = nullptr;
			uint8 newEdgeIndex = 0xFF;

			const Face *face = &polyhedron->face[a];
			uint32 faceEdgeCount = face->edgeCount;
			Face *resultFace = &result->face[remappedFaceIndex];
			uint32 resultFaceEdgeCount = 0;

			for (uint32 b = 0; b < faceEdgeCount; b++)		// Loop over face's original edges.
			{
				uint8 edgeIndex = face->edgeIndex[b];
				int32 code = edgeCode[edgeIndex];
				if (code & 1)	
				{
					// One endpoint on negative side of plane, and other either
					// on positive side (code == 3) or in plane (code == 1).
					if (!newEdge)
					{
						// At this point, we know we need a new edge.
						newEdgeIndex = resultEdgeCount;
						newEdge = &result->edge[resultEdgeCount];
						planeEdgeTable[planeEdgeCount++] = uint8(resultEdgeCount++);
						*newEdge = Edge{{0xFF, 0xFF}, {remappedFaceIndex, 0xFF}};
					}

					const Edge *edge = &polyhedron->edge[edgeIndex];
					bool ccw = (edge->faceIndex[0] == a);
					bool insertEdge = ccw ^ (vertexCode[edge->vertexIndex[0]] == 0);

					if (code == 3)	// Original edge has been clipped.
					{
						uint8 remappedEdgeIndex = edgeRemap[edgeIndex];
						resultFace->edgeIndex[resultFaceEdgeCount++] = remappedEdgeIndex;
						const Edge *resultEdge = &result->edge[remappedEdgeIndex];
						if (insertEdge)
						{
							newEdge->vertexIndex[0] = resultEdge->vertexIndex[ccw];
							resultFace->edgeIndex[resultFaceEdgeCount++] = newEdgeIndex;
						}
						else
						{
							newEdge->vertexIndex[1] = resultEdge->vertexIndex[!ccw];
						}
					}
					else		// Original edge has been deleted, code == 1.
					{
						if (insertEdge)
						{
							newEdge->vertexIndex[0] = vertexRemap[edge->vertexIndex[!ccw]];
							resultFace->edgeIndex[resultFaceEdgeCount++] = newEdgeIndex;
						}
						else
						{
							newEdge->vertexIndex[1] = vertexRemap[edge->vertexIndex[ccw]];
						}
					}
				}
				else if (code != 0)
				{
					// Neither endpoint is on the negative side of the clipping plane.
					uint8 remappedEdgeIndex = edgeRemap[edgeIndex];
					resultFace->edgeIndex[resultFaceEdgeCount++] = remappedEdgeIndex;
					if (code == 2) planeEdgeTable[planeEdgeCount++] = remappedEdgeIndex;
				}
			}

			if ((newEdge) && (max(newEdge->vertexIndex[0], newEdge->vertexIndex[1]) == 0xFF))
			{
				// The input polyhedron was invalid.
				*result = *polyhedron;
				return (true);
			}

			resultFace->edgeCount = uint8(resultFaceEdgeCount);
		}
	}

	if (planeEdgeCount > 2)
	{
		result->plane[resultFaceCount] = plane;
		Face *resultFace = &result->face[resultFaceCount];
		resultFace->edgeCount = uint8(planeEdgeCount);

		for (uint32 a = 0; a < planeEdgeCount; a++)
		{
			uint8 edgeIndex = planeEdgeTable[a];
			resultFace->edgeIndex[a] = edgeIndex;

			Edge *resultEdge = &result->edge[edgeIndex];
			uint8 k = (resultEdge->faceIndex[1] == 0xFF);
			resultEdge->faceIndex[k] = uint8(resultFaceCount);
		}

		resultFaceCount++;
	}

	result->vertexCount = uint8(resultVertexCount);
	result->edgeCount = uint8(resultEdgeCount);
	result->faceCount = uint8(resultFaceCount);
	return (true);
}

// Listing 9.7

float CalculateDiameter(int32 vertexCount, const Point3D *vertex, int32 *a, int32 *b)
{
	constexpr int32 kDirectionCount = 13;

	static const float direction[kDirectionCount][3] =
	{
		{1.0F, 0.0F, 0.0F}, {0.0F, 1.0F, 0.0F}, {0.0F, 0.0F, 1.0F},
		{1.0F, 1.0F, 0.0F}, {1.0F, 0.0F, 1.0F}, {0.0F, 1.0F, 1.0F},
		{1.0F, -1.0F, 0.0F}, {1.0F, 0.0F, -1.0F}, {0.0F, 1.0F, -1.0F},
		{1.0F, 1.0F, 1.0F}, {1.0F, -1.0F, 1.0F}, {1.0F, 1.0F, -1.0F}, {1.0F, -1.0F, -1.0F}
	};

	float		dmin[kDirectionCount], dmax[kDirectionCount];
	int32		imin[kDirectionCount], imax[kDirectionCount];

	// Find min and max dot products for each direction and record vertex indices.
	for (int32 j = 0; j < kDirectionCount; j++)
	{
		const float *u = direction[j];
		dmin[j] = dmax[j] = u[0] * vertex[0].x + u[1] * vertex[0].y + u[2] * vertex[0].z;
		imin[j] = imax[j] = 0;

		for (int32 i = 1; i < vertexCount; i++)
		{
			float d = u[0] * vertex[i].x + u[1] * vertex[i].y + u[2] * vertex[i].z;
			if (d < dmin[j])
			{
				dmin[j] = d;
				imin[j] = i;
			}
			else if (d > dmax[j])
			{
				dmax[j] = d;
				imax[j] = i;
			}
		}
	}

	// Find direction for which vertices at min and max extents are furthest apart.
	float d2 = SquaredMagnitude(vertex[imax[0]] - vertex[imin[0]]);

	int32 k = 0;
	for (int32 j = 1; j < kDirectionCount; j++)
	{
		float m2 = SquaredMagnitude(vertex[imax[j]] - vertex[imin[j]]);
		if (m2 > d2)
		{
			d2 = m2;
			k = j;
		}
	}

	*a = imin[k];
	*b = imax[k];
	return (d2);
}

// Listing 9.8

float CalculateBoundingSphere(int32 vertexCount, const Point3D *vertex, Point3D *center)
{
	int32		a, b;

	// Determine initial center and radius.
	float d2 = CalculateDiameter(vertexCount, vertex, &a, &b);
	*center = (vertex[a] + vertex[b]) * 0.5F;
	float radius = sqrt(d2) * 0.5F;

	// Make pass through vertices and adjust sphere as necessary.
	for (int32 i = 0; i < vertexCount; i++)
	{
		Vector3D pv = vertex[i] - *center;
		float m2 = SquaredMagnitude(pv);
		if (m2 > radius * radius)
		{
			Point3D q = *center - (pv * (radius / sqrt(m2)));
			*center = (q + vertex[i]) * 0.5F;
			radius = Magnitude(q - *center);
		}
	}

	return (radius);
}

// Listing 9.9

void CalculateAxisAlignedBoundingBox(int32 vertexCount, const Point3D *vertex, Point3D *center, Vector3D *size)
{
	Point3D vmin = vertex[0];
	Point32 vmax = vertex[0];
	for (int32 i = 1; i < vertexCount; i++)
	{
		vmin = Min(vmin, vertex[i]);
		vmax = Max(vmax, vertex[i]);
	}

	*center = (vmin + vmax) * 0.5F;
	*size = (vmax - vmin) * 0.5F;
}

// Listing 9.10

Vector3D MakePerpendicularVector(const Vector3D& v)
{
	float x = Fabs(v.x);
	float y = Fabs(v.y);
	float z = Fabs(v.z);
	if (z < Fmin(x, y)) return (Vector3D(v.y, -v.x, 0.0F));
	if (y < x) return (Vector3D(-v.z, 0.0F, v.x));
	return (Vector3D(0.0F, v.z, -v.y));
}

// Listing 9.11

void CalculateSecondaryDiameter(int32 vertexCount, const Point3D *vertex, const Vector3D& axis, int32 *a, int32 *b)
{
	constexpr int32 kDirectionCount = 4;

	static const float direction[kDirectionCount][2] =
	{
		{1.0F, 0.0F}, {0.0F, 1.0F}, {1.0F, 1.0F}, {1.0F, -1.0F}
	};

	float		dmin[kDirectionCount], dmax[kDirectionCount];
	int32		imin[kDirectionCount], imax[kDirectionCount];

	// Create vectors x and y perpendicular to the primary axis.
	Vector3D x = MakePerpendicularVector(axis);
	Vector3D y = Cross(axis, x);

	// Find min and max dot products for each direction and record vertex indices.
	for (int32 j = 0; j < kDirectionCount; j++)
	{
		Vector3D t = x * direction[j][0] + y * direction[j][1];
		dmin[j] = dmax[j] = Dot(t, vertex[0]);
		imin[j] = imax[j] = 0;

		for (int32 i = 1; i < vertexCount; i++)
		{
			float d = Dot(t, vertex[i]);
			if (d < dmin[j])
			{
				dmin[j] = d;
				imin[j] = i;
			}
			else if (d > dmax[j])
			{
				dmax[j] = d;
				imax[j] = i;
			}
		}
	}

	// Find diameter in plane perpendicular to primary axis.
	Vector3D dv = vertex[imax[0]] - vertex[imin[0]];
	float d2 = SquaredMagnitude(dv - axis * Dot(dv, axis));

	int32 k = 0;
	for (int32 j = 1; j < kDirectionCount; j++)
	{
		dv = vertex[imax[j]] - vertex[imin[j]];
		float m2 = SquaredMagnitude(dv - axis * Dot(dv, axis));
		if (m2 > d2)
		{
			d2 = m2;
			k = j;
		}
	}

	*a = imin[k];
	*b = imax[k];
}

// Listing 9.12

void FindExtremalVertices(int32 vertexCount, const Point3D *vertex, const Plane& plane, int32 *e, int32 *f)
{
	*e = 0;
	*f = 0;

	float dmin = Dot(plane, vertex[0]);
	float dmax = dmin;

	for (int32 i = 1; i < vertexCount; i++)
	{
		float m = Dot(plane, vertex[i]);
		if (m < dmin)
		{
			dmin = m;
			*e = i;
		}
		else if (m > dmax)
		{
			dmax = m;
			*f = i;
		}
	}
}

void GetPrimaryBoxDirections(int32 vertexCount, const Point3D *vertex, int32 a, int32 b, Vector3D *direction)
{
	int32 c = 0;
	direction[0] = vertex[b] - vertex[a];
	float dmax = DistPointLine(vertex[0], vertex[a], direction[0]);
	for (int32 i = 1; i < vertexCount; i++)
	{
		float m = DistPointLine(vertex[i], vertex[a], direction[0]);
		if (m > dmax)
		{
			dmax = m;
			c = i;
		}
	}

	direction[1] = vertex[c] - vertex[a];
	direction[2] = vertex[c] - vertex[b];
	Vector3D normal = Cross(direction[0], direction[1]);
	Plane plane(normal, -Dot(normal, vertex[a]));

	int32 e, f;

	FindExtremalVertices(vertexCount, vertex, plane, &e, &f);
	direction[3] = vertex[e] - vertex[a];
	direction[4] = vertex[e] - vertex[b];
	direction[5] = vertex[e] - vertex[c];
	direction[6] = vertex[f] - vertex[a];
	direction[7] = vertex[f] - vertex[b];
	direction[8] = vertex[f] - vertex[c];
}

void GetSecondaryBoxDirections(int32 vertexCount, const Point3D *vertex, const Vector3D& axis, int32 a, int32 b, Vector3D *direction)
{
	direction[0] = vertex[b] - vertex[a];
	Vector3D normal = Cross(axis, direction[0]);
	Plane plane(normal, -Dot(normal, vertex[a]));

	int32 e, f;

	FindExtremalVertices(vertexCount, vertex, plane, &e, &f);
	direction[1] = vertex[e] - vertex[a];
	direction[2] = vertex[e] - vertex[b];
	direction[3] = vertex[f] - vertex[a];
	direction[4] = vertex[f] - vertex[b];

	for (int32 j = 0; j < 5; j++) direction[j] -= axis * Dot(direction[j], axis);
}

// Listing 9.13

void CalculateOrientedBoundingBox(int32 vertexCount, const Point3D *vertex, Point3D *center, Vector3D *size, Vector3D *axis)
{
	int32			a, b;
	Vector3D		primaryDirection[9], secondaryDirection[5];

	CalculateDiameter(vertexCount, vertex, &a, &b);
	GetPrimaryBoxDirections(vertexCount, vertex, a, b, primaryDirection);

	float area = FLT_MAX;
	for (int32 k = 0; k < 9; k++)		// Loop over all candidates for primary axis.
	{
		Vector3D s = Normalize(primaryDirection[k]);
		CalculateSecondaryDiameter(vertexCount, vertex, s, &a, &b);
		GetSecondaryBoxDirections(vertexCount, vertex, s, a, b, secondaryDirection);

		for (int32 j = 0; j < 5; j++)		// Loop over all candidates for secondary axis.
		{
			Vector3D t = Normalize(secondaryDirection[j]), u = Cross(s, t);
			float smin = Dot(s, vertex[0]);
			float smax = smin;
			float tmin = Dot(t, vertex[0]);
			float tmax = tmin;
			float umin = Dot(u, vertex[0]);
			float umax = umin;

			for (int32 i = 1; i < vertexCount; i++)
			{
				float ds = Dot(s, vertex[i]);
				float dt = Dot(t, vertex[i]);
				float du = Dot(u, vertex[i]);
				smin = Fmin(smin, ds);
				smax = Fmax(smax, ds);
				tmin = Fmin(tmin, dt);
				tmax = Fmax(tmax, dt);
				umin = Fmin(umin, du);
				umax = Fmax(umax, du);
			}

			float hx = (smax - smin) * 0.5F;
			float hy = (tmax - tmin) * 0.5F,
			float hz = (umax - umin) * 0.5F;

			// Calculate one-eighth surface area and see if it's better.
			float m = hx * hy + hy * hz + hz * hx;
			if (m < area)
			{
				*center = (s * (smin + smax) + t * (tmin + tmax) + u * (umin + umax)) * 0.5F;
				size->Set(hx, hy, hz);
				axis[0] = s;
				axis[1] = t;
				axis[2] = u;
				area = m;
			}
		}
	}
}

// Listing 9.14

void BuildFrustumPolyhedron(const Transform4D& Mcam, float g, float s, float n, float f, Polyhedron *polyhedron)
{
	polyhedron->vertexCount = 8;
	polyhedron->edgeCount = 12;
	polyhedron->faceCount = 6;

	// Generate vertices for the near side.
	float y = n / g;
	float x = y * s;
	polyhedron->vertex[0] = Mcam * Point3D(x, y, n);
	polyhedron->vertex[1] = Mcam * Point3D(x, -y, n);
	polyhedron->vertex[2] = Mcam * Point3D(-x, -y, n);
	polyhedron->vertex[3] = Mcam * Point3D(-x, y, n);

	// Generate vertices for the far side.
	y = f / g;
	x = y * s;
	polyhedron->vertex[4] = Mcam * Point3D(x, y, f);
	polyhedron->vertex[5] = Mcam * Point3D(x, -y, f);
	polyhedron->vertex[6] = Mcam * Point3D(-x, -y, f);
	polyhedron->vertex[7] = Mcam * Point3D(-x, y, f);

	// Generate lateral planes.
	Transform4D inverse = Inverse(Mcam);
	float mx = 1.0F / sqrt(g * g + s * s);
	float my = 1.0F / sqrt(g * g + 1.0F);
	polyhedron->plane[0] = Plane(-g * mx, 0.0F, s * mx, 0.0F) * inverse;
	polyhedron->plane[1] = Plane(0.0F, g * my, my, 0.0F) * inverse;
	polyhedron->plane[2] = Plane(g * mx, 0.0F, s * mx, 0.0F) * inverse;
	polyhedron->plane[3] = Plane(0.0F, -g * my, my, 0.0F) * inverse;

	// Generate near and far planes.
	float d = Dot(Mcam[2], Mcam[3]);
	polyhedron->plane[4].Set(Mcam[2], -(d + n));
	polyhedron->plane[5].Set(-Mcam[2], d + f);

	// Generate all edges and lateral faces.
	Edge *edge = polyhedron->edge;
	Face *face = polyhedron->face;
	for (int32 i = 0; i < 4; i++, edge++, face++)
	{
		edge[0].vertexIndex[0] = uint8(i);
		edge[0].vertexIndex[1] = uint8(i + 4);
		edge[0].faceIndex[0] = uint8(i);
		edge[0].faceIndex[1] = uint8((i - 1) & 3);

		edge[4].vertexIndex[0] = uint8(i);
		edge[4].vertexIndex[1] = uint8((i + 1) & 3);
		edge[4].faceIndex[0] = 4;
		edge[4].faceIndex[1] = uint8(i);

		edge[8].vertexIndex[0] = uint8(((i + 1) & 3) + 4);
		edge[8].vertexIndex[1] = uint8(i + 4);
		edge[8].faceIndex[0] = 5;
		edge[8].faceIndex[1] = uint8(i);

		face->edgeCount = 4;
		face->edgeIndex[0] = uint8(i);
		face->edgeIndex[1] = uint8((i + 1) & 3);
		face->edgeIndex[2] = uint8(i + 4);
		face->edgeIndex[3] = uint8(i + 8);
	}

	// Generate near and far faces.
	face[0].edgeCount = face[1].edgeCount = 4;
	face[0].edgeIndex[0] = 4;
	face[0].edgeIndex[1] = 5;
	face[0].edgeIndex[2] = 6;
	face[0].edgeIndex[3] = 7;
	face[1].edgeIndex[0] = 8;
	face[1].edgeIndex[1] = 9;
	face[1].edgeIndex[2] = 10;
	face[1].edgeIndex[3] = 11;
}

// Listing 9.15

bool SphereVisible(int32 planeCount, const Plane *planeArray, const Point3D& center, float radius)
{
	float negativeRadius = -radius;
	for (int32 i = 0; i < planeCount; i++)
	{
		if (Dot(planeArray[i], center) <= negativeRadius) return (false);
	}

	return (true);
}

// Listing 9.16

bool OrientedBoxVisible(int32 planeCount, const Plane *planeArray, const Point3D& center, const Vector3D& size, const Vector3D *axis)
{
	for (int32 i = 0; i < planeCount; i++)
	{
		const Plane& g = planeArray[i];
		float rg = fabs(Dot(g, axis[0]) * size.x) + fabs(Dot(g, axis[1]) * size.y) + fabs(Dot(g, axis[2]) * size.z);
		if (Dot(g, center) <= -rg) return (false);
	}

	return (true);
}

bool AxisAlignedBoxVisible(int32 planeCount, const Plane *planeArray, const Point3D& center, const Vector3D& size)
{
	for (int32 i = 0; i < planeCount; i++)
	{
		const Plane& g = planeArray[i];
		float rg = fabs(g.x * size.x) + fabs(g.y * size.y) + fabs(g.z * size.z);
		if (Dot(g, center) <= -rg) return (false);
	}

	return (true);
}

// Listing 9.17

bool OrientedBoxIlluminated(const Point3D& lightPosition, float rmax, const Point3D& center, const Vector3D& size, const Vector3D *axis)
{
	Vector3D v = center - lightPosition;
	float vs = fabs(Dot(v, axis[0]));
	float vt = fabs(Dot(v, axis[1]));
	float vu = fabs(Dot(v, axis[2]));

	float v2 = Dot(v, v);
	float m = (size.x * vs + size.y * vt + size.z * vu) * rsqrt(v2) + rmax;
	if (v2 >= m * m) return (false);

	return (fmax(fmax(vs - size.x, vt - size.y), vu - size.z) < rmax);
}

// Listing 9.18

int32 CalculateShadowRegion(const Polyhedron *polyhedron, const Vector4D& lightPosition, Plane *shadowPlane)
{
	const float kShadowRegionEpsilon = 1.0e-6F;
	bool frontArray[kMaxPolyhedronFaceCount];

	// Classify faces of polyhedron and record back planes.
	int32 shadowPlaneCount = 0;
	int32 cameraPlaneCount = polyhedron->faceCount;
	for (int32 i = 0; i < cameraPlaneCount; i++)
	{
		const Plane& plane = polyhedron->plane[i];
		frontArray[i] = (Dot(plane, lightPosition) < 0.0F);
		if (!frontArray[i]) shadowPlane[shadowPlaneCount++] = plane;
	}

	// Construct planes containing silhouette edges and light position.
	const Edge *edge = polyhedron->edge;
	int32 edgeCount = polyhedron->edgeCount;
	for (int32 i = 0; i < edgeCount; i++, edge++)
	{
		bool front = frontArray[edge->faceIndex[0]];
		if (front ^ frontArray[edge->faceIndex[1]])
		{
			// This edge is on the silhouette.
			const Point3D& v0 = polyhedron->vertex[edge->vertexIndex[0]];
			const Point3D& v1 = polyhedron->vertex[edge->vertexIndex[1]];
			Vector3D n = Cross(lightPosition.xyz() - v0 * lightPosition.w, v1 - v0);

			// Make sure plane is not degenerate.
			float m = SquaredMagnitude(n);
			if (m > kShadowRegionEpsilon)
			{
				// Normalize and point inward.
				n *= ((front) ? 1.0F : -1.0F) / sqrt(m);
				shadowPlane[shadowPlaneCount].Set(n, -Dot(n, v0));
				if (++shadowPlaneCount == kMaxPolyhedronFaceCount) break;
			}
		}
	}

	return (shadowPlaneCount);
}

// Listing 9.19

void BuildPortalRegion(int32 vertexCount, const Point3D *portalVertex, const Plane& portalPlane, const Plane& backPlane, const Point3D& cameraPosition, Polyhedron *polyhedron)
{
	polyhedron->vertexCount = vertexCount * 2;
	polyhedron->edgeCount = vertexCount * 3;
	polyhedron->faceCount = vertexCount + 2;

	Point3D *v = polyhedron->vertex;
	Point3D *w = v + vertexCount;
	Plane *plane = polyhedron->plane;

	// Calculate lateral planes and vertex positions.
	float bc = Dot(backPlane, cameraPosition);
	Vector3D u0 = portalVertex[vertexCount - 1] - cameraPosition;
	for (int32 k = 0; k < vertexCount; k++)
	{
		Vector3D u1 = portalVertex[k] - cameraPosition;
		Vector3D normal = Normalize(Cross(u1, u0));
		plane[k].Set(normal, -Dot(normal, cameraPosition));

		v[k] = portalVertex[k];
		w[k] = cameraPosition - u1 * (bc / Dot(backPlane, u1));
		u0 = u1;
	}

	// Generate front and back planes.
	plane[vertexCount] = -portalPlane;
	plane[vertexCount + 1] = backPlane;

	// Generate all edges and lateral faces.
	int32 i = vertexCount - 1;
	Edge *edge = polyhedron->edge;
	Face *face = polyhedron->face;
	for (int32 j = 0; j < vertexCount; i = j++, edge++, face++)
	{
		int32 k = j + 1;
		k &= (k - vertexCount) >> 8;   // if j + 1 >= n, then k = 0.

		edge[0].vertexIndex[0] = uint8(j);
		edge[0].vertexIndex[1] = uint8(j + vertexCount);
		edge[0].faceIndex[1] = uint8(j);
		edge[0].faceIndex[0] = uint8(k);

		edge[vertexCount].vertexIndex[0] = uint8(j);
		edge[vertexCount].vertexIndex[1] = uint8(k);
		edge[vertexCount].faceIndex[1] = uint8(k);
		edge[vertexCount].faceIndex[0] = uint8(vertexCount);

		edge[vertexCount * 2].vertexIndex[0] = uint8(k + vertexCount);
		edge[vertexCount * 2].vertexIndex[1] = uint8(j + vertexCount);
		edge[vertexCount * 2].faceIndex[1] = uint8(k);
		edge[vertexCount * 2].faceIndex[0] = uint8(vertexCount + 1);

		face->edgeCount = 4;
		face->edgeIndex[0] = uint8(j);
		face->edgeIndex[1] = uint8(i + vertexCount * 2);
		face->edgeIndex[2] = uint8(i);
		face->edgeIndex[3] = uint8(i + vertexCount);
	}

	// Generate front and back faces.
	face[0].edgeCount = vertexCount;
	face[1].edgeCount = vertexCount;
	for (int32 k = 0; k < vertexCount; k++)
	{
		face[0].edgeIndex[k] = uint8(k + vertexCount);
		face[1].edgeIndex[k] = uint8(k + vertexCount * 2);
	}
}

// Listing 9.20

const uint8 occluderPolygonIndex[43] =		// 3-bit vertex count, 5-bit polygon index.
{
	0x00, 0x80, 0x81, 0, 0x82, 0xC9, 0xC8, 0, 0x83, 0xC7, 0xC6, 0, 0, 0, 0, 0,
	0x84, 0xCF, 0xCE, 0, 0xD1, 0xD9, 0xD8, 0, 0xD0, 0xD7, 0xD6, 0, 0, 0, 0, 0,
	0x85, 0xCB, 0xCA, 0, 0xCD, 0xD5, 0xD4, 0, 0xCC, 0xD3, 0xD2
};

const uint8 occluderVertexIndex[26][6] =	// Vertex indices for all 26 polygons.
{
	{1, 3, 7, 5}, {2, 0, 4, 6}, {3, 2, 6, 7}, {0, 1, 5, 4}, {4, 5, 7, 6}, {1, 0, 2, 3},
	{2, 0, 1, 5, 4, 6}, {0, 1, 3, 7, 5, 4}, {3, 2, 0, 4, 6, 7}, {1, 3, 2, 6, 7, 5},
	{1, 0, 4, 6, 2, 3}, {5, 1, 0, 2, 3, 7}, {4, 0, 2, 3, 1, 5}, {0, 2, 6, 7, 3, 1},
	{0, 4, 5, 7, 6, 2}, {4, 5, 1, 3, 7, 6}, {1, 5, 7, 6, 4, 0}, {5, 7, 3, 2, 6, 4},
	{3, 1, 5, 4, 6, 2}, {2, 3, 7, 5, 4, 0}, {1, 0, 4, 6, 7, 3}, {0, 2, 6, 7, 5, 1},
	{7, 6, 2, 0, 1, 5}, {6, 4, 0, 1, 3, 7}, {5, 7, 3, 2, 0, 4}, {4, 5, 1, 3, 2, 6}
};

const float occluderVertexPosition[8][3] =	// Normalized vertex coords for unit cube.
{
	{0.0F, 0.0F, 0.0F}, {1.0F, 0.0F, 0.0F}, {0.0F, 1.0F, 0.0F}, {1.0F, 1.0F, 0.0F},
	{0.0F, 0.0F, 1.0F}, {1.0F, 0.0F, 1.0F}, {0.0F, 1.0F, 1.0F}, {1.0F, 1.0F, 1.0F}
};

int32 MakeOcclusionRegion(const Vector3D& size, const Plane *frustumPlane, const Transform4D& Mocc, const Transform4D& Mcam, Plane *occluderPlane)
{
	Point3D		polygonVertex[2][10];
	float		vertexDistance[2][4];
	float		vertexLocation[10];

	const float kOccluderEpsilon = 0.002F;
	uint32 occlusionCode = 0;
	int32 planeCount = 0;

	// Transform the camera position into the occluder's object space.
	Transform4D m = Inverse(Mocc);
	Point3D cameraPosition = m * Mcam.GetTranslation();

	// Calculate 6-bit occlusion code and generate front planes.
	uint32 axisCode = 0x01;
	for (int32 i = 0; i < 3; i++, axisCode <<= 2)
	{
		if (cameraPosition[i] > size[i])
		{
			occlusionCode |= axisCode;
			occluderPlane[planeCount++].Set(-m(i,0), -m(i,1), -m(i,2), size[i] - m(i,3));
		}
		else if (cameraPosition[i] < 0.0F)
		{
			occlusionCode |= axisCode << 1;
			occluderPlane[planeCount++].Set(m(i,0), m(i,1), m(i,2), m(i,3));
		}
	}

	// Look up silhouette polygon with occlusion code.
	uint32 polygonIndex = occluderPolygonIndex[occlusionCode];
	const uint8 *vertexIndex = occluderVertexIndex[polygonIndex & 0x1F];
	int32 vertexCount = polygonIndex >> 5;

	// Generate silhouette vertices in camera space.
	Transform4D McamInverse = Inverse(Mcam);
	Transform4D t = McamInverse * Mocc;
	for (int32 i = 0; i < vertexCount; i++)
	{
		const float *p = occluderVertexPosition[vertexIndex[i]];
		polygonVertex[1][i] = t * Point3D(p[0] * size.x, p[1] * size.y, p[2] * size.z);
	}

	// Clip silhouette to lateral planes of view frustum.
	const Point3D *vertex = polygonVertex[1];
	for (int32 k = 0; k < 4; k++)
	{
		const Plane& plane = frustumPlane[k];
		Point3D *result = polygonVertex[k & 1];
		vertexCount = ClipPolygon(vertexCount, vertex, plane, vertexLocation, result);
		vertex = result;
	}

	if (vertexCount < 3) return (0);

	// Generate occlusion region planes in world space.
	const Vector3D *v1 = &vertex[vertexCount - 1];
	for (int32 k = 0; k < 4; k++) vertexDistance[0][k] = Dot(frustumPlane[k].GetNormal(), *v1);

	for (int32 i = 0; i < vertexCount; i++)
	{
		bool cull = false; int32 j = i & 1;
		const Vector3D *v2 = &vertex[i];
		Vector3D planeNormal = Normalize(Cross(*v2, *v1));
		for (int32 k = 0; k < 4; k++)
		{
			const Vector3D& frustumNormal = frustumPlane[k].GetNormal();
			float d = Dot(frustumNormal, *v2);
			vertexDistance[j ^ 1][k] = d;

			// Cull edge lying in frustum plane, but only if its extrusion points inward.
			if ((fmax(d, vertexDistance[j][k]) < kOccluderEpsilon) && (Dot(planeNormal, frustumNormal) > 0.0F)) cull = true;
		}

		if (!cull) occluderPlane[planeCount++] = Plane(planeNormal, 0.0F) * McamInverse;
		v1 = v2;
	}

	return (planeCount);
}

// Listing 9.21

int32 CalculateParallelFogOcclusionPlanes(const Plane& fogPlane, const Point3D& cameraPosition, float fogDensity, float maxOpticalDepth, Plane *occlusionPlane)
{
	float z0 = Dot(fogPlane, cameraPosition);
	float z02 = z0 * z0;
	float sigma = 2.0F * maxOpticalDepth / fogDensity;

	// Calculate the plane below the camera, which always exists.
	float zmin = -sqrt(z02 + sigma);
	occlusionPlane[0].Set(-fogPlane.GetNormal(), zmin - fogPlane.w);

	// Return if the plane above camera does not exist.
	if (z02 < sigma) return (1);

	// Calculate the plane above the camera.
	float zmax = -sqrt(z02 - sigma);
	occlusionPlane[1].Set(fogPlane.GetNormal(), fogPlane.w - zmax);
	return (2);
}

// Listing 9.22

bool CalculatePerpendicularFogOcclusionPlane(const Plane& fogPlane, const Point3D& cameraPosition, const Vector3D& viewDirection, float fogDensity, float maxOpticalDepth, Plane *occlusionPlane)
{
	// Project view direction onto fog plane.
	Vector3D normal = Reject(viewDirection, fogPlane.GetNormal());
	float n2 = Dot(normal, normal);
	if (n2 > FLT_MIN)
	{
		// View direction has nonzero horizontal component.
		float z0 = Dot(fogPlane, cameraPosition);
		float z02 = z0 * z0;
		float z0_inv2 = 1.0F / z02;
		float sigma = 2.0F * maxOpticalDepth / fogDensity;
		float sigma2 = sigma * sigma;
		float m = sigma2 * z0_inv2 * z0_inv2;
		float r = 0.0F;

		if (m < 1.6875F)
		{
			// When m < 27/16, a root of h(u) exists.
			float u = 1.0F - m * 0.125F;
			float u2 = u * u;

			// Apply second iteration of Newton's method.
			u -= (((u + 2.0F) * u2 - 2.0F) * u + m - 1.0F) / ((u * 4.0F + 6.0F) * u2 - 2.0F);

			// Plug root into Equation (9.44).
			float up1 = u + 1.0F;
			float um1 = u - 1.0F;
			r = sqrt(sigma2 * z0_inv2 / (up1 * up1) - um1 * um1 * z02);
		}

		if (m > 1.0F)
		{
			// Calculate occlusion distance with Equation (9.50).
			// Take larger value of r when both solutions are valid.
			r = fmax(-z0 * sqrt(m - 1.0F), r);
		}

		// Construct occlusion plane perpendicular to fog plane.
		normal *= 1.0F / sqrt(n2);
		occlusionPlane->Set(normal, -Dot(normal, cameraPosition) - r);
		return (true);
	}

	return (false);
}

// Listing 10.1

uniform float3 cameraPosition;
uniform float3 cameraDown;

float3 CalculateSphericalBillboardVertexPosition(float3 center, float2 billboard)
{
	float3 n = normalize(cameraPosition - center);
	float3 a = normalize(cross(n, cameraDown));
	float3 b = cross(n, a);
	return (center + a * billboard.x + b * billboard.y);
}

// Listing 10.2

uniform float3 cameraPosition;

float3 CalculateCylindricalBillboardVertexPosition(float3 center, float billboard)
{
	float2 a = float2(center.y - cameraPosition.y, cameraPosition.x - center.x);
	a *= rsqrt(max(dot(a, a), 0.0001));
	return (float3(center.xy + a * billboard, center.z));
}

// Listing 10.3

uniform float3 cameraPosition;

float3 CalculatePolyboardVertexPosition(float3 center, float3 tangent, float billboard)
{
	float3 a = cross(tangent, cameraPosition - center);
	a *= rsqrt(max(dot(a, a), 0.0001));
	return (center + a * billboard);
}

// Listing 10.4

float4 CalculateStructureOutput(float z)
{
	float h = asfloat(asuint(z) & 0xFFFFE000U);
	return (float4(ddx(z), ddy(z), h, z - h));
}

// Listing 10.5

uniform TextureRect		structureBuffer;

float CalculateDepthFade(float2 pixelCoord, float z, float scale)
{
	float2 depth = texture(structureBuffer, pixelCoord).zw;
	float delta = depth.x + depth.y - z;
	return (saturate(scale * delta));
}

// Listing 10.6

uniform TextureRect		structureBuffer;
uniform float3			cameraPosition, cameraView;
uniform float			R2, recipR2, recip3R2, normalizer;

float CalculateHaloBrightness(float3 pobject, float2 pixelCoord)
{
	float3 vdir = cameraPosition - pobject;
	float v2 = dot(vdir, vdir);
	float p2 = dot(pobject, pobject);
	float pv = -dot(pobject, vdir);
	float m = sqrt(max(pv * pv - v2 * (p2 - R2), 0.0));

	// Read z0 from the structure buffer.
	float2 depth = texture(structureBuffer, pixelCoord).zw;
	float t0 = 1.0 + (depth.x + depth.y) / dot(cameraView, vdir);

	// Calculate clamped limits of integration.
	float t1 = clamp((pv - m) / v2, t0, 1.0);
	float t2 = clamp((pv + m) / v2, t0, 1.0);
	float u1 = t1 * t1;
	float u2 = t2 * t2;

	// Evaluate density integral, normalize, and square.
	float B = ((1.0 - p2 * recipR2) * (t2 - t1) + pv * recipR2 * (u2 - u1) - v2 * recip3R2 * (t2 * u2 - t1 * u1)) * normalizer;
	return (B * B * v2);
}

// Listing 10.7

uniform TextureRect		structureBuffer;
uniform float3			cameraView;
uniform float			shaftSigma, shaftRho0, shaftTau, normalizer;

float CalculateShaftBrightness(float pz, float3 vdir, float2 pixelCoord, float t1, float t2)
{
	// Read z0 from the structure buffer, calculate t0, and clamp to [t0,1].
	float2 depth = texture(structureBuffer, pixelCoord).zw;
	float t0 = 1.0 + (depth.x + depth.y) / dot(cameraView, vdir);
	t1 = clamp(t1, t0, 1.0);
	t2 = clamp(t2, t0, 1.0);

	// Limit to range where density is not negative.
	float tlim = (shaftTau - pz) / vdir.z;
	if (vdir.z * shaftSigma < 0.0)
	{
		t1 = min(t1, tlim);
		t2 = min(t2, tlim);
	}
	else
	{
		t1 = max(t1, tlim);
		t2 = max(t2, tlim);
	}

	// Evaluate density integral, normalize, and square.
	float B = (shaftSigma * (pz + vdir.z * ((t1 + t2) * 0.5)) + shaftRho0) * (t2 - t1) * normalizer;
	return (B * B * dot(vdir, vdir));
}

// Listing 10.8

uniform float3			cameraPosition;
uniform float			rx2, ry2, rx2ry2;

float CalculateCylinderShaftBrightness(float3 pobject, float2 pixelCoord)
{
	float3 vdir = cameraPosition - pobject;
	float2 v2 = vdir.xy * vdir.xy;
	float2 p2 = pobject.xy * pobject.xy;

	// Calculate quadratic coefficients.
	float a = ry2 * v2.x + rx2 * v2.y;
	float b = -ry2 * pobject.x * vdir.x - rx2 * pobject.y * vdir.y;
	float c = ry2 * p2.x + rx2 * p2.y - rx2ry2;
	float m = sqrt(max(b * b - a * c, 0.0));

	// Calculate limits and integrate.
	float t1 = max((b - m) / a, 0.0);
	float t2 = max((b + m) / a, 0.0);
	return (CalculateShaftBrightness(pobject.z, vdir, pixelCoord, t1, t2));
}

// Listing 10.9

uniform float3			cameraPosition;
uniform float			sx, sy;

float CalculateBoxShaftBrightness(float3 pobject, float2 pixelCoord)
{
	float3 vdir = cameraPosition - pobject;
	float t1 = 0.0;
	float t2 = 1.0;

	// Find intersections with planes perpendicular to x axis.
	float a = -pobject.x / vdir.x;
	float b = (sx - pobject.x) / vdir.x;
	if (vdir.x > 0.0)
	{
		t1 = max(t1, a);
		t2 = min(t2, b);
	}
	else
	{
		t1 = max(t1, b);
		t2 = min(t2, a);
	}

	// Find intersections with planes perpendicular to y axis.
	a = -pobject.y / vdir.y;
	b = (sy - pobject.y) / vdir.y;
	if (vdir.y > 0.0)
	{
		t1 = max(t1, a);
		t2 = min(t2, b);
	}
	else
	{
		t1 = max(t1, b);
		t2 = min(t2, a);
	}

	return (CalculateShaftBrightness(pobject.z, vdir, pixelCoord, t1, t2));
}

// Listing 10.10

uniform TextureRect		structureBuffer;
uniform Texture2D		rotationTexture;
uniform float			vectorScale, intensity;

float CalculateAmbientOcclusion(float2 pixelCoord)
{
	const float kTangentTau = 0.03125;

	// These are the offset vectors used for the four samples.
	const float dx[4] = {0.1, 0.0, -0.3, 0.0};
	const float dy[4] = {0.0, 0.2, 0.0, -0.4};

	// Sample the structure buffer at the central pixel.
	float4 structure = texture(structureBuffer, pixelCoord);
	float z0 = structure.z + structure.w;

	// Calculate the normal vector.
	float scale = vectorScale * z0;
	float3 normal = normalize(float3(structure.xy, -scale));
	scale = 1.0 / scale;

	// Fetch a cos/sin pair from the 4x4 rotation texture.
	float2 rot = texture(rotationTexture, pixelCoord * 0.25).xy;
	float occlusion = 0.0;
	float weight = 0.0;

	for (int i = 0; i < 4; i++)
	{
		float3		v;

		// Calculate the rotated offset vector for this sample.
		v.x = rot.x * dx[i] - rot.y * dy[i];
		v.y = rot.y * dx[i] + rot.x * dy[i];

		// Fetch the depth from the structure buffer at the offset location.
		float2 depth = texture(structureBuffer, (pixelCoord + v.xy * scale)).zw;
		v.z = depth.x + depth.y - z0;

		// Calculate w(v) and f(v), and accumulate H(v) = w(v)f(v).
		float d = dot(normal, v);
		float w = saturate(1.0 - d * 0.5);
		float c = saturate(d * rsqrt(dot(v, v)) - kTangentTau);

		occlusion += w - w * sqrt(1.0 - c * c);
		weight += w;
	}

	// Return the ambient light factor.
	return (1.0 - occlusion * intensity / max(weight, 0.0001));
}

// Listing 10.11

uniform TextureRect		structureBuffer;
uniform TextureRect		occlusionBuffer;

float BlurAmbientOcclusion(float2 pixelCoord)
{
	const float kDepthDelta = 0.0078125;

	// Use depth and gradient to calculate a valid range for the blur samples.
	float4 structure = texture(structureBuffer, pixelCoord);
	float range = (max(abs(structure.x), abs(structure.y)) + kDepthDelta) * 1.5;
	float z0 = structure.z + structure.w;

	float2 sample = float2(0.0, 1.0);
	float3 occlusion = float3(0.0, 0.0, 0.0);
	for (int j = 0; j < 2; j++)
	{
		float y = float(j * 2) - 0.5;
		for (int i = 0; i < 2; i++)
		{
			float x = float(i * 2) - 0.5;

			// Fetch a filtered sample and accumulate.
			float2 sampleCoord = pixelCoord + float2(x, y);
			sample.x = texture(occlusionBuffer, sampleCoord).x;
			occlusion.z += sample.x;

			// If depth at sample is in range, acculumate the occlusion value.
			float2 depth = texture(structureBuffer, sampleCoord).zw;
			if (abs(depth.x + depth.y - z0) < range) occlusion.xy += sample;
		}
	}

	// Divide the accumulated occlusion value by the number of samples that passed.
	return ((occlusion.y > 0.0) ? occlusion.x / occlusion.y : occlusion.z * 0.25);
}

// Listing 10.12

uniform float3		frustumParams;      // (dmin * s/g, dmin * 1/g, dmin)
uniform float3		shadowMatrix[3];	

void CalculateAtmosphericShadowingRays(float2 position, out float2 cameraRay, out float3 shadowRay)
{
	// Calculate point on camera-space plane z = dmin.
	float3 q = float3(position.xy * frustumParams.xy, frustumParams.z);
	cameraRay = q.xy;

	// Transform camera-space ray direction into shadow space for cascade 0.
	shadowRay.x = dot(q, shadowMatrix[0]);
	shadowRay.y = dot(q, shadowMatrix[1]);
	shadowRay.z = dot(q, shadowMatrix[2]);
}

// Listing 10.13

uniform TextureRect				structureBuffer;
uniform Texture2DArrayShadow	shadowTexture;
uniform Texture2D				noiseTexture;
uniform float					minAtmosphereDepth;			// dmin
uniform float					atmosphereDepthRatio;		// dmax / dmin
uniform float					atmosphereBrightness;		// lambda * (dmax - dmin) / (n + 1)
uniform float					maxCascadeDepth[4];
uniform float3					cascadeScale[3];
uniform float3					cascadeOffset[3];
uniform float3					shadowCameraPosition;
uniform float3					cameraLightDirection;
uniform float3					anisotropyConst;			// (1 - g, 1 + g * g, 2 * g)
uniform float2					noiseShift;

float CalculateAtmosShadowing(float2 pixelCoord, float2 cameraRay, float3 shadowRay)
{
	float4		sampleCoord;

	const float m = 0.25;						// Slope at first sample.
	const float dt = 1.0 / 64.0;				// Delta t = 1 / n for n samples.
	const float dw = 2.0 * (1.0 - m) * dt;		// Change in weight per step.

	// Fetch depth of solid surface from structure buffer at pixel location.
	float2 strc = texture(structureBuffer, pixelCoord * 2.0).zw;
	float depth = strc.x + strc.y;

	// Calculate scale making camera ray have length of dmin.
	float dmin = minAtmosphereDepth;
	float invlength = rsqrt(dot(cameraRay, cameraRay) + dmin * dmin);
	float scale = dmin * invlength;

	// Calculate begin and end points in shadow space.
	float3 p1 = shadowRay * scale;
	float3 p2 = p1 * atmosphereDepthRatio + shadowCameraPosition;
	p1 += shadowCameraPosition;

	// Calculate begin and end depths in camera space.
	float z1 = dmin * scale;
	float z2 = z1 * atmosphereDepthRatio;

	// Calculate Henyey-Greenstein function and multiply by density.
	float h = anisotropyConst.x * rsqrt(anisotropyConst.y - anisotropyConst.z * dot(float3(cameraRay, dmin), cameraLightDirection) * invlength);
	float intensity = h * h * h * atmosphereBrightness;

	// Fetch random offset from 32x32 noise texture.
	float t = texture(noiseTexture, pixelCoord * 0.03125 + noiseShift).x * dt;
	float tmax = t + 1.0;
	float atm = 0.0;			// Atmsophere accumulator.
	float weight = m;			// First weight is always m.

	// Start with cascade 0.
	sampleCoord.z = 0.0;
	float zmax = min(depth, maxCascadeDepth[0]);

	for (; t <= tmax; t += dt)
	{
		float u = (t * (1.0 - m) + m) * t;
		float z = lerp(z1, z2, u);
		if (z > zmax) break;

		sampleCoord.xyw = lerp(p1, p2, u);							// Calculate position.
		atm += texture(shadowTexture, sampleCoord) * weight;		// Accumulate sample.
		weight += dw;												// Increase weight.
	}

	for (int cascade = 1; cascade < 4; cascade++)
	{
		// Handle remaining cascades.
		sampleCoord.z = float(cascade);
		zmax = min(depth, maxCascadeDepth[cascade]);

		for (; t <= tmax; t += dt)
		{
			float u = (t * (1.0 - m) + m) * t;
			float z = lerp(z1, z2, u);
			if (z > zmax) break;

			// Calculate position and transform into current cascade space.
			sampleCoord.xyw = lerp(p1, p2, u) * cascadeScale[cascade - 1].xyz + cascadeOffset[cascade - 1].xyz;
			atm += texture(shadowTexture, sampleCoord) * weight;
			weight += dw;
		}
	}

	return (atm * intensity);
}

// Listing 10.14

uniform TextureRect		atmosphereBuffer;
uniform float3			lightPosition;

float GetAtmosphereIntensity(float2 pixelCoord)
{
	float2 center = pixelCoord * 0.5;
	float2 direction = normalize(lightPosition.xy - center * lightPosition.z);
	float a = texture(atmosphereBuffer, center).x;
	float b = texture(atmosphereBuffer, center + direction).x;
	float c = texture(atmosphereBuffer, center - direction).x;
	return (min(max(min(a, b), c), max(a, b)));
}

// Listing 10.15

uniform float4 newMotionMatrix[3];
uniform float4 oldMotionMatrix[3];

void TransformMotionBlurPositions(float3 pobject, out float3 pnew, out float3 pold)
{
	pnew.x = dot(newMotionMatrix[0].xyz, pobject) + newMotionMatrix[0].w;
	pnew.y = dot(newMotionMatrix[1].xyz, pobject) + newMotionMatrix[1].w;
	pnew.z = dot(newMotionMatrix[2].xyz, pobject) + newMotionMatrix[2].w;

	pold.x = dot(oldMotionMatrix[0].xyz, pobject) + oldMotionMatrix[0].w;
	pold.y = dot(oldMotionMatrix[1].xyz, pobject) + oldMotionMatrix[1].w;
	pold.z = dot(oldMotionMatrix[2].xyz, pobject) + oldMotionMatrix[2].w;
}

// Listing 10.16

uniform float velocityScale;

float2 CalculateMotionBlurVelocity(float3 pnew, float3 pold)
{
	float2 velocity = (pnew.xy / pnew.z - pold.xy / pold.z) * velocityScale;
	return (velocity / max(length(velocity), 1.0) * 0.5 + 0.5);
}

// Listing 10.17

uniform TextureRect		colorBuffer;
uniform TextureRect		velocityBuffer;
uniform float			vstep;

float3 ApplySimpleMotionBlur(float2 pixelCoord)
{
	// Read color buffer and velocity buffer at center pixel.
	float3 color = texture(colorBuffer, pixelCoord).xyz;
	float2 velocity = texture(velocityBuffer, pixelCoord).xy * 2.0 - 1.0;

	// Add 8 more samples along velocity direction.
	for (int i = 1; i <= 4; i++)
	{
		float dp = float(i) * vstep;
		color += texture(colorBuffer, pixelCoord + velocity * dp).xyz;
		color += texture(colorBuffer, pixelCoord - velocity * dp).xyz;
	}

	// Return average of all samples.
	return (color * 0.1111111);
}

// Listing 10.18

uniform TextureRect		colorBuffer;
uniform TextureRect		structureBuffer;
uniform TextureRect		velocityBuffer;
uniform float			rmax, vstep;

float3 ApplyComplexMotionBlur(float2 pixelCoord)
{
	const float kDepthDelta = 0.0078125;
	const float kVeloSigma = 0.0625;

	// Read color buffer, structure buffer, and velocity buffer at center pixel.
	float4 color = float4(texture(colorBuffer, pixelCoord).xyz, 1.0);
	float4 structure = texture(structureBuffer, pixelCoord);
	float2 velocity = texture(velocityBuffer, pixelCoord).xy * 2.0 - 1.0;

	// Use gradient to calculate minimum depth for other samples.
	float zmin = structure.z + structure.w - (abs(dot(velocity, structure.xy)) + kDepthDelta) * rmax;

	float4 sample = float4(0.0, 0.0, 0.0, 1.0);
	for (int i = 1; i <= 4; i++)
	{
		float dp = float(i) * vstep;

		// Read all three buffers at sample location, and subtract center velocity.
		float2 sampleCoord = pixelCoord + velocity * dp;
		sample.xyz = texture(colorBuffer, sampleCoord).xyz;
		float2 depth = texture(structureBuffer, sampleCoord).zw;
		float2 dv = texture(velocityBuffer, sampleCoord).xy * 2.0 - 1.0 - velocity;

		// Add sample only if depth greater than minimum or velocities similar enough.
		if ((depth.x + depth.y > zmin) || (dot(dv, dv) < kVeloSigma)) color += sample;

		// Repeat in opposite direction.
		sampleCoord = pixelCoord - velocity * dp;
		sample.xyz = texture(colorBuffer, sampleCoord).xyz;
		depth = texture(structureBuffer, sampleCoord).zw;

		dv = texture(velocityBuffer, sampleCoord).xy * 2.0 - 1.0 - velocity;
		if ((depth.x + depth.y > zmin) || (dot(dv, dv) < kVeloSigma)) color += sample;
	}

	// Divide color sum by number of passing samples.
	return (color.xyz / color.w);
}

// Listing 10.19

bool CalculateDeviceSpaceExtents(const Point3D& bmin, const Point3D& bmax, const Matrix4D& mvp, Point2D *vmin, Point2D *vmax)
{
	static const int8 edgeVertexIndex[24] =
	{
		0, 1, 1, 2, 2, 3, 3, 0, 4, 5, 5, 6, 6, 7, 7, 4, 0, 4, 1, 5, 2, 6, 3, 7
	};

	// Transform bounding box vertices into clip space.
	Vector4D vertex[8];
	vertex[0] = mvp * Point3D(bmin.x, bmin.y, bmin.z);
	vertex[1] = mvp * Point3D(bmax.x, bmin.y, bmin.z);
	vertex[2] = mvp * Point3D(bmax.x, bmax.y, bmin.z);
	vertex[3] = mvp * Point3D(bmin.x, bmax.y, bmin.z);
	vertex[4] = mvp * Point3D(bmin.x, bmin.y, bmax.z);
	vertex[5] = mvp * Point3D(bmax.x, bmin.y, bmax.z);
	vertex[6] = mvp * Point3D(bmax.x, bmax.y, bmax.z);
	vertex[7] = mvp * Point3D(bmin.x, bmax.y, bmax.z);

	// Initialize device-space rectangle.
	float xmin = FLT_MAX, xmax = -FLT_MAX;
	float ymin = FLT_MAX, ymax = -FLT_MAX;

	const int8 *edge = edgeVertexIndex;
	for (int32 i = 0; i < 12; i++, edge += 2)
	{
		Vector4D p1 = vertex[edge[0]], p2 = vertex[edge[1]];

		// Clip edge against near plane.
		if (p1.z < 0.0F)
		{
			if (p2.z < 0.0F) continue;		// Edge completely clipped away.
			Vector4D dp = p1 - p2;
			p1 -= dp * (p1.z / dp.z);
		}
		else if (p2.z < 0.0F)
		{
			Vector4D dp = p2 - p1;
			p2 -= dp * (p2.z / dp.z);
		}

		// Perform perspective divide.
		float f1 = 1.0F / p1.w;
		float f2 = 1.0F / p2.w;
		p1.x *= f1; p1.y *= f1;
		p2.x *= f2; p2.y *= f2;

		// Update device-space rectangle.
		xmin = fmin(xmin, fmin(p1.x, p2.x));
		xmax = fmax(xmax, fmax(p1.x, p2.x));
		ymin = fmin(ymin, fmin(p1.y, p2.y));
		ymax = fmax(ymax, fmax(p1.y, p2.y));
	}

	if (xmin < xmax)		// At least one edge was visible.
	{
		vmin->Set(xmin, ymin);
		vmax->Set(xmax, ymax);
		return (true);
	}

	return (false);
}

// Listing 10.20

typedef int8 Voxel;

// The GetVoxel() function loads a single value from the scalar field at coords (i,j,k).
inline Voxel GetVoxel(const Voxel *field, int32 n, int32 m, int32 i, int32 j, int32 k)
{
	return (field[(k * m + j) * n + i]);
}

uint32 LoadCell(const Voxel *field, int32 n, int32 m, int32 i, int32 j, int32 k, Voxel *distance)
{
	distance[0] = GetVoxel(field, n, m, i, j, k);
	distance[1] = GetVoxel(field, n, m, i + 1, j, k);
	distance[2] = GetVoxel(field, n, m, i, j + 1, k);
	distance[3] = GetVoxel(field, n, m, i + 1, j + 1, k);
	distance[4] = GetVoxel(field, n, m, i, j, k + 1);
	distance[5] = GetVoxel(field, n, m, i + 1, j, k + 1);
	distance[6] = GetVoxel(field, n, m, i, j + 1, k + 1);
	distance[7] = GetVoxel(field, n, m, i + 1, j + 1, k + 1);

	// Concatenate sign bits of the voxel values to form the case index for the cell.
	return (((distance[0] >> 7) & 0x01) | ((distance[1] >> 6) & 0x02)
	      | ((distance[2] >> 5) & 0x04) | ((distance[3] >> 4) & 0x08)
	      | ((distance[4] >> 3) & 0x10) | ((distance[5] >> 2) & 0x20)
	      | ((distance[6] >> 1) & 0x40) |  (distance[7] & 0x80));
}

// Listing 10.21

const uint8 equivClassTable[256] =
{
	// Equivalence class index for each of the 256 possible cases.

	0x00, 0x01, 0x01, 0x03, 0x01, 0x03, 0x02, 0x04, 0x01, 0x02, 0x03, 0x04, 0x03, 0x04, 0x04, 0x03,
	0x01, 0x03, 0x02, 0x04, 0x02, 0x04, 0x06, 0x0C, 0x02, 0x05, 0x05, 0x0B, 0x05, 0x0A, 0x07, 0x04,
	0x01, 0x02, 0x03, 0x04, 0x02, 0x05, 0x05, 0x0A, 0x02, 0x06, 0x04, 0x0C, 0x05, 0x07, 0x0B, 0x04,
	0x03, 0x04, 0x04, 0x03, 0x05, 0x0B, 0x07, 0x04, 0x05, 0x07, 0x0A, 0x04, 0x08, 0x0E, 0x0E, 0x03,
	0x01, 0x02, 0x02, 0x05, 0x03, 0x04, 0x05, 0x0B, 0x02, 0x06, 0x05, 0x07, 0x04, 0x0C, 0x0A, 0x04,
	0x03, 0x04, 0x05, 0x0A, 0x04, 0x03, 0x07, 0x04, 0x05, 0x07, 0x08, 0x0E, 0x0B, 0x04, 0x0E, 0x03,
	0x02, 0x06, 0x05, 0x07, 0x05, 0x07, 0x08, 0x0E, 0x06, 0x09, 0x07, 0x0F, 0x07, 0x0F, 0x0E, 0x0D,
	0x04, 0x0C, 0x0B, 0x04, 0x0A, 0x04, 0x0E, 0x03, 0x07, 0x0F, 0x0E, 0x0D, 0x0E, 0x0D, 0x02, 0x01,
	0x01, 0x02, 0x02, 0x05, 0x02, 0x05, 0x06, 0x07, 0x03, 0x05, 0x04, 0x0A, 0x04, 0x0B, 0x0C, 0x04,
	0x02, 0x05, 0x06, 0x07, 0x06, 0x07, 0x09, 0x0F, 0x05, 0x08, 0x07, 0x0E, 0x07, 0x0E, 0x0F, 0x0D,
	0x03, 0x05, 0x04, 0x0B, 0x05, 0x08, 0x07, 0x0E, 0x04, 0x07, 0x03, 0x04, 0x0A, 0x0E, 0x04, 0x03,
	0x04, 0x0A, 0x0C, 0x04, 0x07, 0x0E, 0x0F, 0x0D, 0x0B, 0x0E, 0x04, 0x03, 0x0E, 0x02, 0x0D, 0x01,
	0x03, 0x05, 0x05, 0x08, 0x04, 0x0A, 0x07, 0x0E, 0x04, 0x07, 0x0B, 0x0E, 0x03, 0x04, 0x04, 0x03,
	0x04, 0x0B, 0x07, 0x0E, 0x0C, 0x04, 0x0F, 0x0D, 0x0A, 0x0E, 0x0E, 0x02, 0x04, 0x03, 0x0D, 0x01,
	0x04, 0x07, 0x0A, 0x0E, 0x0B, 0x0E, 0x0E, 0x02, 0x0C, 0x0F, 0x04, 0x0D, 0x04, 0x0D, 0x03, 0x01,
	0x03, 0x04, 0x04, 0x03, 0x04, 0x03, 0x0D, 0x01, 0x04, 0x0D, 0x03, 0x01, 0x03, 0x01, 0x01, 0x00
};

struct ClassData
{
	uint8		geometryCounts;		// Vertex and triangle counts.
	uint8		vertexIndex[15];	// List of 3-15 vertex indices.
};

const ClassData classGeometryTable[16] =
{
	// Triangulation data for each equivalence class.

	{0x00, {}},
	{0x31, {0, 1, 2}},
	{0x62, {0, 1, 2, 3, 4, 5}},
	{0x42, {0, 1, 2, 0, 2, 3}},
	{0x53, {0, 1, 4, 1, 3, 4, 1, 2, 3}},
	{0x73, {0, 1, 2, 0, 2, 3, 4, 5, 6}},
	{0x93, {0, 1, 2, 3, 4, 5, 6, 7, 8}},
	{0x84, {0, 1, 4, 1, 3, 4, 1, 2, 3, 5, 6, 7}},
	{0x84, {0, 1, 2, 0, 2, 3, 4, 5, 6, 4, 6, 7}},
	{0xC4, {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11}},
	{0x64, {0, 4, 5, 0, 1, 4, 1, 3, 4, 1, 2, 3}},
	{0x64, {0, 5, 4, 0, 4, 1, 1, 4, 3, 1, 3, 2}},
	{0x64, {0, 4, 5, 0, 3, 4, 0, 1, 3, 1, 2, 3}},
	{0x64, {0, 1, 2, 0, 2, 3, 0, 3, 4, 0, 4, 5}},
	{0x75, {0, 1, 2, 0, 2, 3, 0, 3, 4, 0, 4, 5, 0, 5, 6}},
	{0x95, {0, 4, 5, 0, 3, 4, 0, 1, 3, 1, 2, 3, 6, 7, 8}}
};

const uint16 vertexCodeTable[256][12] =
{
	// List of 3-12 vertex codes for each of the 256 possible cases.
	// The meanings of the bit fields are shown in Figure 10.33.

	{},
	{0xA188, 0x9050, 0x72E0},
	{0xA188, 0x65E9, 0x8359},
	{0x9050, 0x72E0, 0x65E9, 0x8359},
	{0x9050, 0x849A, 0x18F2},
	{0x72E0, 0xA188, 0x849A, 0x18F2},
	{0xA188, 0x65E9, 0x8359, 0x9050, 0x849A, 0x18F2},
	{0x849A, 0x18F2, 0x72E0, 0x65E9, 0x8359},
	{0x8359, 0x0BFB, 0x849A},
	{0xA188, 0x9050, 0x72E0, 0x849A, 0x8359, 0x0BFB},
	{0xA188, 0x65E9, 0x0BFB, 0x849A},
	{0x9050, 0x72E0, 0x65E9, 0x0BFB, 0x849A},
	{0x9050, 0x8359, 0x0BFB, 0x18F2},
	{0x8359, 0x0BFB, 0x18F2, 0x72E0, 0xA188},
	{0xA188, 0x65E9, 0x0BFB, 0x18F2, 0x9050},
	{0x72E0, 0x65E9, 0x0BFB, 0x18F2},
	{0x72E0, 0x1674, 0x27AC},
	{0xA188, 0x9050, 0x1674, 0x27AC},
	{0xA188, 0x65E9, 0x8359, 0x72E0, 0x1674, 0x27AC},
	{0x65E9, 0x8359, 0x9050, 0x1674, 0x27AC},
	{0x9050, 0x849A, 0x18F2, 0x72E0, 0x1674, 0x27AC},
	{0x1674, 0x27AC, 0xA188, 0x849A, 0x18F2},
	{0x72E0, 0x1674, 0x27AC, 0xA188, 0x65E9, 0x8359, 0x9050, 0x849A, 0x18F2},
	{0x849A, 0x18F2, 0x1674, 0x27AC, 0x65E9, 0x8359},
	{0x849A, 0x8359, 0x0BFB, 0x72E0, 0x1674, 0x27AC},
	{0xA188, 0x9050, 0x1674, 0x27AC, 0x849A, 0x8359, 0x0BFB},
	{0x849A, 0xA188, 0x65E9, 0x0BFB, 0x72E0, 0x1674, 0x27AC},
	{0x849A, 0x0BFB, 0x65E9, 0x27AC, 0x1674, 0x9050},
	{0x9050, 0x8359, 0x0BFB, 0x18F2, 0x72E0, 0x1674, 0x27AC},
	{0x8359, 0x0BFB, 0x18F2, 0x1674, 0x27AC, 0xA188},
	{0xA188, 0x65E9, 0x0BFB, 0x18F2, 0x9050, 0x72E0, 0x1674, 0x27AC},
	{0x27AC, 0x65E9, 0x0BFB, 0x18F2, 0x1674},
	{0x65E9, 0x27AC, 0x097D},
	{0xA188, 0x9050, 0x72E0, 0x65E9, 0x27AC, 0x097D},
	{0x8359, 0xA188, 0x27AC, 0x097D},
	{0x27AC, 0x097D, 0x8359, 0x9050, 0x72E0},
	{0x9050, 0x849A, 0x18F2, 0x65E9, 0x27AC, 0x097D},
	{0xA188, 0x849A, 0x18F2, 0x72E0, 0x65E9, 0x27AC, 0x097D},
	{0xA188, 0x27AC, 0x097D, 0x8359, 0x9050, 0x849A, 0x18F2},
	{0x849A, 0x18F2, 0x72E0, 0x27AC, 0x097D, 0x8359},
	{0x849A, 0x8359, 0x0BFB, 0x65E9, 0x27AC, 0x097D},
	{0xA188, 0x9050, 0x72E0, 0x849A, 0x8359, 0x0BFB, 0x65E9, 0x27AC, 0x097D},
	{0x0BFB, 0x849A, 0xA188, 0x27AC, 0x097D},
	{0x9050, 0x72E0, 0x27AC, 0x097D, 0x0BFB, 0x849A},
	{0x9050, 0x8359, 0x0BFB, 0x18F2, 0x65E9, 0x27AC, 0x097D},
	{0x8359, 0x0BFB, 0x18F2, 0x72E0, 0xA188, 0x65E9, 0x27AC, 0x097D},
	{0x9050, 0x18F2, 0x0BFB, 0x097D, 0x27AC, 0xA188},
	{0x097D, 0x0BFB, 0x18F2, 0x72E0, 0x27AC},
	{0x65E9, 0x72E0, 0x1674, 0x097D},
	{0xA188, 0x9050, 0x1674, 0x097D, 0x65E9},
	{0x72E0, 0x1674, 0x097D, 0x8359, 0xA188},
	{0x8359, 0x9050, 0x1674, 0x097D},
	{0x65E9, 0x72E0, 0x1674, 0x097D, 0x9050, 0x849A, 0x18F2},
	{0x18F2, 0x849A, 0xA188, 0x65E9, 0x097D, 0x1674},
	{0x72E0, 0x1674, 0x097D, 0x8359, 0xA188, 0x9050, 0x849A, 0x18F2},
	{0x18F2, 0x1674, 0x097D, 0x8359, 0x849A},
	{0x65E9, 0x72E0, 0x1674, 0x097D, 0x849A, 0x8359, 0x0BFB},
	{0xA188, 0x9050, 0x1674, 0x097D, 0x65E9, 0x849A, 0x8359, 0x0BFB},
	{0x72E0, 0x1674, 0x097D, 0x0BFB, 0x849A, 0xA188},
	{0x849A, 0x9050, 0x1674, 0x097D, 0x0BFB},
	{0x65E9, 0x72E0, 0x1674, 0x097D, 0x9050, 0x8359, 0x0BFB, 0x18F2},
	{0xA188, 0x8359, 0x0BFB, 0x18F2, 0x1674, 0x097D, 0x65E9},
	{0xA188, 0x72E0, 0x1674, 0x097D, 0x0BFB, 0x18F2, 0x9050},
	{0x18F2, 0x1674, 0x097D, 0x0BFB},
	{0x18F2, 0x0ABE, 0x1674},
	{0xA188, 0x9050, 0x72E0, 0x18F2, 0x0ABE, 0x1674},
	{0xA188, 0x65E9, 0x8359, 0x18F2, 0x0ABE, 0x1674},
	{0x9050, 0x72E0, 0x65E9, 0x8359, 0x18F2, 0x0ABE, 0x1674},
	{0x9050, 0x849A, 0x0ABE, 0x1674},
	{0x72E0, 0xA188, 0x849A, 0x0ABE, 0x1674},
	{0x9050, 0x849A, 0x0ABE, 0x1674, 0xA188, 0x65E9, 0x8359},
	{0x1674, 0x0ABE, 0x849A, 0x8359, 0x65E9, 0x72E0},
	{0x8359, 0x0BFB, 0x849A, 0x18F2, 0x0ABE, 0x1674},
	{0xA188, 0x9050, 0x72E0, 0x849A, 0x8359, 0x0BFB, 0x18F2, 0x0ABE, 0x1674},
	{0xA188, 0x65E9, 0x0BFB, 0x849A, 0x18F2, 0x0ABE, 0x1674},
	{0x9050, 0x72E0, 0x65E9, 0x0BFB, 0x849A, 0x18F2, 0x0ABE, 0x1674},
	{0x0ABE, 0x1674, 0x9050, 0x8359, 0x0BFB},
	{0xA188, 0x8359, 0x0BFB, 0x0ABE, 0x1674, 0x72E0},
	{0xA188, 0x65E9, 0x0BFB, 0x0ABE, 0x1674, 0x9050},
	{0x1674, 0x72E0, 0x65E9, 0x0BFB, 0x0ABE},
	{0x72E0, 0x18F2, 0x0ABE, 0x27AC},
	{0x18F2, 0x0ABE, 0x27AC, 0xA188, 0x9050},
	{0x72E0, 0x18F2, 0x0ABE, 0x27AC, 0xA188, 0x65E9, 0x8359},
	{0x18F2, 0x0ABE, 0x27AC, 0x65E9, 0x8359, 0x9050},
	{0x9050, 0x849A, 0x0ABE, 0x27AC, 0x72E0},
	{0xA188, 0x849A, 0x0ABE, 0x27AC},
	{0x9050, 0x849A, 0x0ABE, 0x27AC, 0x72E0, 0xA188, 0x65E9, 0x8359},
	{0x8359, 0x849A, 0x0ABE, 0x27AC, 0x65E9},
	{0x72E0, 0x18F2, 0x0ABE, 0x27AC, 0x849A, 0x8359, 0x0BFB},
	{0x18F2, 0x0ABE, 0x27AC, 0xA188, 0x9050, 0x849A, 0x8359, 0x0BFB},
	{0x72E0, 0x18F2, 0x0ABE, 0x27AC, 0x849A, 0xA188, 0x65E9, 0x0BFB},
	{0x9050, 0x18F2, 0x0ABE, 0x27AC, 0x65E9, 0x0BFB, 0x849A},
	{0x72E0, 0x27AC, 0x0ABE, 0x0BFB, 0x8359, 0x9050},
	{0x0BFB, 0x0ABE, 0x27AC, 0xA188, 0x8359},
	{0x9050, 0xA188, 0x65E9, 0x0BFB, 0x0ABE, 0x27AC, 0x72E0},
	{0x65E9, 0x0BFB, 0x0ABE, 0x27AC},
	{0x65E9, 0x27AC, 0x097D, 0x18F2, 0x0ABE, 0x1674},
	{0xA188, 0x9050, 0x72E0, 0x65E9, 0x27AC, 0x097D, 0x18F2, 0x0ABE, 0x1674},
	{0xA188, 0x27AC, 0x097D, 0x8359, 0x18F2, 0x0ABE, 0x1674},
	{0x27AC, 0x097D, 0x8359, 0x9050, 0x72E0, 0x18F2, 0x0ABE, 0x1674},
	{0x849A, 0x0ABE, 0x1674, 0x9050, 0x65E9, 0x27AC, 0x097D},
	{0x72E0, 0xA188, 0x849A, 0x0ABE, 0x1674, 0x65E9, 0x27AC, 0x097D},
	{0x849A, 0x0ABE, 0x1674, 0x9050, 0xA188, 0x27AC, 0x097D, 0x8359},
	{0x72E0, 0x27AC, 0x097D, 0x8359, 0x849A, 0x0ABE, 0x1674},
	{0x849A, 0x8359, 0x0BFB, 0x65E9, 0x27AC, 0x097D, 0x18F2, 0x0ABE, 0x1674},
	{0xA188, 0x9050, 0x72E0, 0x849A, 0x8359, 0x0BFB, 0x65E9, 0x27AC, 0x097D, 0x18F2, 0x0ABE, 0x1674},
	{0x0BFB, 0x849A, 0xA188, 0x27AC, 0x097D, 0x18F2, 0x0ABE, 0x1674},
	{0x849A, 0x9050, 0x72E0, 0x27AC, 0x097D, 0x0BFB, 0x18F2, 0x0ABE, 0x1674},
	{0x0ABE, 0x1674, 0x9050, 0x8359, 0x0BFB, 0x65E9, 0x27AC, 0x097D},
	{0xA188, 0x8359, 0x0BFB, 0x0ABE, 0x1674, 0x72E0, 0x65E9, 0x27AC, 0x097D},
	{0x0BFB, 0x0ABE, 0x1674, 0x9050, 0xA188, 0x27AC, 0x097D},
	{0x72E0, 0x27AC, 0x097D, 0x0BFB, 0x0ABE, 0x1674},
	{0x097D, 0x65E9, 0x72E0, 0x18F2, 0x0ABE},
	{0x0ABE, 0x097D, 0x65E9, 0xA188, 0x9050, 0x18F2},
	{0x0ABE, 0x18F2, 0x72E0, 0xA188, 0x8359, 0x097D},
	{0x0ABE, 0x097D, 0x8359, 0x9050, 0x18F2},
	{0x9050, 0x849A, 0x0ABE, 0x097D, 0x65E9, 0x72E0},
	{0x65E9, 0xA188, 0x849A, 0x0ABE, 0x097D},
	{0x72E0, 0x9050, 0x849A, 0x0ABE, 0x097D, 0x8359, 0xA188},
	{0x8359, 0x849A, 0x0ABE, 0x097D},
	{0x097D, 0x65E9, 0x72E0, 0x18F2, 0x0ABE, 0x849A, 0x8359, 0x0BFB},
	{0x097D, 0x65E9, 0xA188, 0x9050, 0x18F2, 0x0ABE, 0x849A, 0x8359, 0x0BFB},
	{0x097D, 0x0BFB, 0x849A, 0xA188, 0x72E0, 0x18F2, 0x0ABE},
	{0x9050, 0x18F2, 0x0ABE, 0x097D, 0x0BFB, 0x849A},
	{0x0ABE, 0x097D, 0x65E9, 0x72E0, 0x9050, 0x8359, 0x0BFB},
	{0xA188, 0x8359, 0x0BFB, 0x0ABE, 0x097D, 0x65E9},
	{0xA188, 0x72E0, 0x9050, 0x0BFB, 0x0ABE, 0x097D},
	{0x0BFB, 0x0ABE, 0x097D},
	{0x0BFB, 0x097D, 0x0ABE},
	{0xA188, 0x9050, 0x72E0, 0x0BFB, 0x097D, 0x0ABE},
	{0xA188, 0x65E9, 0x8359, 0x0BFB, 0x097D, 0x0ABE},
	{0x9050, 0x72E0, 0x65E9, 0x8359, 0x0BFB, 0x097D, 0x0ABE},
	{0x9050, 0x849A, 0x18F2, 0x0BFB, 0x097D, 0x0ABE},
	{0xA188, 0x849A, 0x18F2, 0x72E0, 0x0BFB, 0x097D, 0x0ABE},
	{0xA188, 0x65E9, 0x8359, 0x9050, 0x849A, 0x18F2, 0x0BFB, 0x097D, 0x0ABE},
	{0x849A, 0x18F2, 0x72E0, 0x65E9, 0x8359, 0x0BFB, 0x097D, 0x0ABE},
	{0x8359, 0x097D, 0x0ABE, 0x849A},
	{0x849A, 0x8359, 0x097D, 0x0ABE, 0xA188, 0x9050, 0x72E0},
	{0x097D, 0x0ABE, 0x849A, 0xA188, 0x65E9},
	{0x72E0, 0x65E9, 0x097D, 0x0ABE, 0x849A, 0x9050},
	{0x18F2, 0x9050, 0x8359, 0x097D, 0x0ABE},
	{0x097D, 0x8359, 0xA188, 0x72E0, 0x18F2, 0x0ABE},
	{0x18F2, 0x9050, 0xA188, 0x65E9, 0x097D, 0x0ABE},
	{0x0ABE, 0x18F2, 0x72E0, 0x65E9, 0x097D},
	{0x72E0, 0x1674, 0x27AC, 0x0BFB, 0x097D, 0x0ABE},
	{0xA188, 0x9050, 0x1674, 0x27AC, 0x0BFB, 0x097D, 0x0ABE},
	{0xA188, 0x65E9, 0x8359, 0x72E0, 0x1674, 0x27AC, 0x0BFB, 0x097D, 0x0ABE},
	{0x65E9, 0x8359, 0x9050, 0x1674, 0x27AC, 0x0BFB, 0x097D, 0x0ABE},
	{0x9050, 0x849A, 0x18F2, 0x72E0, 0x1674, 0x27AC, 0x0BFB, 0x097D, 0x0ABE},
	{0x1674, 0x27AC, 0xA188, 0x849A, 0x18F2, 0x0BFB, 0x097D, 0x0ABE},
	{0xA188, 0x65E9, 0x8359, 0x9050, 0x849A, 0x18F2, 0x72E0, 0x1674, 0x27AC, 0x0BFB, 0x097D, 0x0ABE},
	{0x8359, 0x849A, 0x18F2, 0x1674, 0x27AC, 0x65E9, 0x0BFB, 0x097D, 0x0ABE},
	{0x849A, 0x8359, 0x097D, 0x0ABE, 0x72E0, 0x1674, 0x27AC},
	{0xA188, 0x9050, 0x1674, 0x27AC, 0x849A, 0x8359, 0x097D, 0x0ABE},
	{0x097D, 0x0ABE, 0x849A, 0xA188, 0x65E9, 0x72E0, 0x1674, 0x27AC},
	{0x65E9, 0x097D, 0x0ABE, 0x849A, 0x9050, 0x1674, 0x27AC},
	{0x18F2, 0x9050, 0x8359, 0x097D, 0x0ABE, 0x72E0, 0x1674, 0x27AC},
	{0x18F2, 0x1674, 0x27AC, 0xA188, 0x8359, 0x097D, 0x0ABE},
	{0x9050, 0xA188, 0x65E9, 0x097D, 0x0ABE, 0x18F2, 0x72E0, 0x1674, 0x27AC},
	{0x18F2, 0x1674, 0x27AC, 0x65E9, 0x097D, 0x0ABE},
	{0x65E9, 0x27AC, 0x0ABE, 0x0BFB},
	{0x65E9, 0x27AC, 0x0ABE, 0x0BFB, 0xA188, 0x9050, 0x72E0},
	{0x8359, 0xA188, 0x27AC, 0x0ABE, 0x0BFB},
	{0x9050, 0x8359, 0x0BFB, 0x0ABE, 0x27AC, 0x72E0},
	{0x65E9, 0x27AC, 0x0ABE, 0x0BFB, 0x9050, 0x849A, 0x18F2},
	{0xA188, 0x849A, 0x18F2, 0x72E0, 0x0BFB, 0x65E9, 0x27AC, 0x0ABE},
	{0x8359, 0xA188, 0x27AC, 0x0ABE, 0x0BFB, 0x9050, 0x849A, 0x18F2},
	{0x8359, 0x849A, 0x18F2, 0x72E0, 0x27AC, 0x0ABE, 0x0BFB},
	{0x65E9, 0x27AC, 0x0ABE, 0x849A, 0x8359},
	{0x65E9, 0x27AC, 0x0ABE, 0x849A, 0x8359, 0xA188, 0x9050, 0x72E0},
	{0xA188, 0x27AC, 0x0ABE, 0x849A},
	{0x72E0, 0x27AC, 0x0ABE, 0x849A, 0x9050},
	{0x9050, 0x8359, 0x65E9, 0x27AC, 0x0ABE, 0x18F2},
	{0x8359, 0x65E9, 0x27AC, 0x0ABE, 0x18F2, 0x72E0, 0xA188},
	{0x9050, 0xA188, 0x27AC, 0x0ABE, 0x18F2},
	{0x72E0, 0x27AC, 0x0ABE, 0x18F2},
	{0x0ABE, 0x0BFB, 0x65E9, 0x72E0, 0x1674},
	{0x9050, 0x1674, 0x0ABE, 0x0BFB, 0x65E9, 0xA188},
	{0x72E0, 0x1674, 0x0ABE, 0x0BFB, 0x8359, 0xA188},
	{0x0BFB, 0x8359, 0x9050, 0x1674, 0x0ABE},
	{0x0ABE, 0x0BFB, 0x65E9, 0x72E0, 0x1674, 0x9050, 0x849A, 0x18F2},
	{0x1674, 0x0ABE, 0x0BFB, 0x65E9, 0xA188, 0x849A, 0x18F2},
	{0x0ABE, 0x0BFB, 0x8359, 0xA188, 0x72E0, 0x1674, 0x9050, 0x849A, 0x18F2},
	{0x8359, 0x849A, 0x18F2, 0x1674, 0x0ABE, 0x0BFB},
	{0x72E0, 0x65E9, 0x8359, 0x849A, 0x0ABE, 0x1674},
	{0x65E9, 0xA188, 0x9050, 0x1674, 0x0ABE, 0x849A, 0x8359},
	{0x1674, 0x0ABE, 0x849A, 0xA188, 0x72E0},
	{0x9050, 0x1674, 0x0ABE, 0x849A},
	{0x0ABE, 0x18F2, 0x9050, 0x8359, 0x65E9, 0x72E0, 0x1674},
	{0xA188, 0x8359, 0x65E9, 0x18F2, 0x1674, 0x0ABE},
	{0xA188, 0x72E0, 0x1674, 0x0ABE, 0x18F2, 0x9050},
	{0x18F2, 0x1674, 0x0ABE},
	{0x18F2, 0x0BFB, 0x097D, 0x1674},
	{0x0BFB, 0x097D, 0x1674, 0x18F2, 0xA188, 0x9050, 0x72E0},
	{0x0BFB, 0x097D, 0x1674, 0x18F2, 0xA188, 0x65E9, 0x8359},
	{0x8359, 0x9050, 0x72E0, 0x65E9, 0x18F2, 0x0BFB, 0x097D, 0x1674},
	{0x0BFB, 0x097D, 0x1674, 0x9050, 0x849A},
	{0xA188, 0x849A, 0x0BFB, 0x097D, 0x1674, 0x72E0},
	{0x0BFB, 0x097D, 0x1674, 0x9050, 0x849A, 0xA188, 0x65E9, 0x8359},
	{0x849A, 0x0BFB, 0x097D, 0x1674, 0x72E0, 0x65E9, 0x8359},
	{0x849A, 0x8359, 0x097D, 0x1674, 0x18F2},
	{0x849A, 0x8359, 0x097D, 0x1674, 0x18F2, 0xA188, 0x9050, 0x72E0},
	{0x1674, 0x097D, 0x65E9, 0xA188, 0x849A, 0x18F2},
	{0x849A, 0x9050, 0x72E0, 0x65E9, 0x097D, 0x1674, 0x18F2},
	{0x8359, 0x097D, 0x1674, 0x9050},
	{0xA188, 0x8359, 0x097D, 0x1674, 0x72E0},
	{0x65E9, 0x097D, 0x1674, 0x9050, 0xA188},
	{0x65E9, 0x097D, 0x1674, 0x72E0},
	{0x27AC, 0x72E0, 0x18F2, 0x0BFB, 0x097D},
	{0xA188, 0x27AC, 0x097D, 0x0BFB, 0x18F2, 0x9050},
	{0x27AC, 0x72E0, 0x18F2, 0x0BFB, 0x097D, 0xA188, 0x65E9, 0x8359},
	{0x27AC, 0x65E9, 0x8359, 0x9050, 0x18F2, 0x0BFB, 0x097D},
	{0x849A, 0x0BFB, 0x097D, 0x27AC, 0x72E0, 0x9050},
	{0x097D, 0x27AC, 0xA188, 0x849A, 0x0BFB},
	{0x27AC, 0x72E0, 0x9050, 0x849A, 0x0BFB, 0x097D, 0x8359, 0xA188, 0x65E9},
	{0x849A, 0x0BFB, 0x097D, 0x27AC, 0x65E9, 0x8359},
	{0x8359, 0x097D, 0x27AC, 0x72E0, 0x18F2, 0x849A},
	{0x18F2, 0x849A, 0x8359, 0x097D, 0x27AC, 0xA188, 0x9050},
	{0x097D, 0x27AC, 0x72E0, 0x18F2, 0x849A, 0xA188, 0x65E9},
	{0x9050, 0x18F2, 0x849A, 0x65E9, 0x097D, 0x27AC},
	{0x72E0, 0x9050, 0x8359, 0x097D, 0x27AC},
	{0x8359, 0x097D, 0x27AC, 0xA188},
	{0x9050, 0xA188, 0x65E9, 0x097D, 0x27AC, 0x72E0},
	{0x65E9, 0x097D, 0x27AC},
	{0x1674, 0x18F2, 0x0BFB, 0x65E9, 0x27AC},
	{0x1674, 0x18F2, 0x0BFB, 0x65E9, 0x27AC, 0xA188, 0x9050, 0x72E0},
	{0xA188, 0x27AC, 0x1674, 0x18F2, 0x0BFB, 0x8359},
	{0x27AC, 0x1674, 0x18F2, 0x0BFB, 0x8359, 0x9050, 0x72E0},
	{0x9050, 0x1674, 0x27AC, 0x65E9, 0x0BFB, 0x849A},
	{0x1674, 0x72E0, 0xA188, 0x849A, 0x0BFB, 0x65E9, 0x27AC},
	{0x0BFB, 0x8359, 0xA188, 0x27AC, 0x1674, 0x9050, 0x849A},
	{0x849A, 0x0BFB, 0x8359, 0x72E0, 0x27AC, 0x1674},
	{0x8359, 0x65E9, 0x27AC, 0x1674, 0x18F2, 0x849A},
	{0x1674, 0x18F2, 0x849A, 0x8359, 0x65E9, 0x27AC, 0xA188, 0x9050, 0x72E0},
	{0x18F2, 0x849A, 0xA188, 0x27AC, 0x1674},
	{0x849A, 0x9050, 0x72E0, 0x27AC, 0x1674, 0x18F2},
	{0x27AC, 0x1674, 0x9050, 0x8359, 0x65E9},
	{0x8359, 0x65E9, 0x27AC, 0x1674, 0x72E0, 0xA188},
	{0xA188, 0x27AC, 0x1674, 0x9050},
	{0x72E0, 0x27AC, 0x1674},
	{0x72E0, 0x18F2, 0x0BFB, 0x65E9},
	{0x9050, 0x18F2, 0x0BFB, 0x65E9, 0xA188},
	{0xA188, 0x72E0, 0x18F2, 0x0BFB, 0x8359},
	{0x9050, 0x18F2, 0x0BFB, 0x8359},
	{0x849A, 0x0BFB, 0x65E9, 0x72E0, 0x9050},
	{0xA188, 0x849A, 0x0BFB, 0x65E9},
	{0x72E0, 0x9050, 0x849A, 0x0BFB, 0x8359, 0xA188},
	{0x8359, 0x849A, 0x0BFB},
	{0x8359, 0x65E9, 0x72E0, 0x18F2, 0x849A},
	{0x18F2, 0x849A, 0x8359, 0x65E9, 0xA188, 0x9050},
	{0x72E0, 0x18F2, 0x849A, 0xA188},
	{0x9050, 0x18F2, 0x849A},
	{0x9050, 0x8359, 0x65E9, 0x72E0},
	{0xA188, 0x8359, 0x65E9},
	{0xA188, 0x72E0, 0x9050},
	{}
};

// Listing 10.22

struct CellStorage
{
	uint16		corner[7];
	uint16		edge[9];
};

inline uint16 *ReuseCornerVertex(int32 n, int32 i, int32 j, CellStorage *const (& deckStorage)[2], uint16 cornerIndex, uint16 deltaCode)
{
	// The corner index in the preceding cell is the sum of the original
	// corner index and the masked delta code.
	cornerIndex += deltaCode;

	// The three bits of the delta code indicate whether one should
	// be subtracted from the cell coords in the x, y, and z directions.
	int32 dx = deltaCode & 1;
	int32 dy = (deltaCode >> 1) & 1;
	int32 dz = deltaCode >> 2;

	// deckStorage[0] points to the current deck, and
	// deckStorage[1] points to the preceding deck
	CellStorage *deck = deckStorage[dz];

	// Return the address of the vertex index in the preceding cell.
	// The new corner index can never be zero.
	return (&deck[(j - dy) * n + (i - dx)].corner[cornerIndex - 1]);
}

inline uint16 ReuseEdgeVertex(int32 n, int32 i, int32 j, CellStorage *const (& deckStorage)[2], uint16 edgeIndex, uint16 deltaCode)
{
	// Edge index in preceding cell differs from original edge index by 3 for each of
	// the lowest three bits in the masked delta code, and by 6 for the highest bit.
	edgeIndex += ((deltaCode & 1) + ((deltaCode >> 1) & 1) + ((deltaCode >> 2) & 3)) * 3;

	// Bits 0, 1, and 3 of the delta code indicate whether one should
	// be subtracted from the cell coords in the x, y, and z directions.
	int32 dx = deltaCode & 1;
	int32 dy = (deltaCode >> 1) & 1;
	int32 dz = deltaCode >> 3;

	// deckStorage[0] points to the current deck, and
	// deckStorage[1] points to the preceding deck
	const CellStorage *deck = deckStorage[dz];

	// Return the vertex index stored in the preceding cell.
	// The new edge index can never be less than 3.
	return (deck[(j - dy) * n + (i - dx)].edge[edgeIndex - 3]);
}

// Listing 10.23

void ProcessCell(const Voxel *field, int32 n, int32 m, int32 i, int32 j, int32 k, CellStorage *const (& deckStorage)[2], uint32 deltaMask, int32& meshVertexCount, int32& meshTriangleCount, Integer3D *meshVertexArray, Triangle *meshTriangleArray)
{
	Voxel		distance[8];

	// Get storage for current cell and set vertex indices at corners
	// to invalid values so those not generated here won't get reused.
	CellStorage *cellStorage = &deckStorage[0][j * n + i];
	for (int32 a = 0; a < 7; a++) cellStorage->corner[a] = 0xFFFF;

	// Call LoadCell() to populate the distance array and get case index.
	uint32 caseIndex = LoadCell(field, n, m, i, j, k, distance);

	// Look up the equivalence class index and use it to look up
	// geometric data for this cell. No geometry if case is 0 or 255.
	int32 equivClass = equivClassTable[caseIndex];
	const ClassData *classData = &classGeometryTable[equivClass];
	uint32 geometryCounts = classData->geometryCounts;

	if (geometryCounts != 0)
	{
		uint16		cellVertexIndex[12];

		int32 vertexCount = geometryCounts >> 4;
		int32 triangleCount = geometryCounts & 0x0F;

		// Look up vertex codes using original case index.
		const uint16 *vertexCode = vertexCodeTable[caseIndex];

		// Duplicate middle bit of delta mask to construct 4-bit mask used for edges.
		uint16 edgeDeltaMask = ((deltaMask << 1) & 0x0C) | (deltaMask & 0x03);

		for (int32 a = 0; a < vertexCount; a++)
		{
			uint16			vertexIndex;
			uint8			corner[2];
			Integer3D		position[2];

			// Extract corner numbers from low 6 bits of vertex code.
			uint16 vcode = vertexCode[a];
			corner[0] = vcode & 0x07;
			corner[1] = (vcode >> 3) & 0x07;

			// Construct integer coordinates of edge's endpoints.
			position[0].x = i + (corner[0] & 1);
			position[0].y = j + ((corner[0] >> 1) & 1);
			position[0].z = k + ((corner[0] >> 2) & 1);
			position[1].x = i + (corner[1] & 1);
			position[1].y = j + ((corner[1] >> 1) & 1);
			position[1].z = k + ((corner[1] >> 2) & 1);

			// Calculate interpolation parameter with Equation (10.95).
			int32 d0 = distance[corner[0]];
			int32 d1 = distance[corner[1]];
			int32 t = (d1 << 8) / (d1 - d0);

			if ((t & 0x00FF) != 0)
			{
				// Vertex falls in the interior of an edge.
				// Extract edge index and delta code from vertex code.
				uint16 edgeIndex = (vcode >> 8) & 0x0F;
				uint16 deltaCode = (vcode >> 12) & edgeDeltaMask;

				if (deltaCode != 0)
				{
					// Reuse vertex from edge in preceding cell.
					vertexIndex = ReuseEdgeVertex(n, i, j, deckStorage, edgeIndex, deltaCode);
				}
				else
				{
					// Generate a new vertex with Equation (10.96).
					vertexIndex = meshVertexCount++;
					Integer3D *vertex = &meshVertexArray[vertexIndex];
					*vertex = position[0] * t + position[1] * (0x0100 - t);

					if (edgeIndex >= 3)
					{
						// Store vertex index for potential reuse later.
						cellStorage->edge[edgeIndex - 3] = vertexIndex;
					}
				}

				cellVertexIndex[a] = vertexIndex;
			}
			else
			{
				// Vertex falls exactly at the first corner of the cell if
				// t == 0, and at the second corner if t == 0x0100.
				uint8 c = (t == 0);
				uint8 cornerIndex = corner[c];

				// Corner vertex in preceding cell may not have been
				// generated, so we get address and look for valid index.
				uint16 *indexAddress = nullptr;
				uint16 deltaCode = (cornerIndex ^ 7) & deltaMask;
				if (deltaCode != 0)
				{
					// Reuse vertex from corner in preceding cell.
					indexAddress = ReuseCornerVertex(n, i, j, deckStorage, cornerIndex, deltaCode);
				}
				else if (cornerIndex != 0)
				{
					// Vertex will be stored for potential reuse later.
					indexAddress = &cellStorage->corner[cornerIndex - 1];
				}

				vertexIndex = (indexAddress) ? *indexAddress : 0xFFFF;
				if (vertexIndex == 0xFFFF)
				{
					// Vertex was not previously generated.
					vertexIndex = meshVertexCount++;
					if (indexAddress) *indexAddress = vertexIndex;

					// Shift corner position to add 8 bits of fraction.
					meshVertexArray[vertexIndex] = position[c] << 8;
				}

				cellVertexIndex[a] = vertexIndex;
			}
		}

		// Generate triangles for this cell using table data.
		const uint8 *classVertexIndex = classData->vertexIndex;
		Triangle *meshTriangle = &meshTriangleArray[meshTriangleCount];
		meshTriangleCount += triangleCount;

		for (int32 a = 0; a < triangleCount; a++)
		{
			meshTriangle[a].vertexIndex[0] = cellVertexIndex[classVertexIndex[0]];
			meshTriangle[a].vertexIndex[1] = cellVertexIndex[classVertexIndex[1]];
			meshTriangle[a].vertexIndex[2] = cellVertexIndex[classVertexIndex[2]];
			classVertexIndex += 3;
		}
	}
}

// Listing 10.24

void ExtractIsosurface(const Voxel *field, int32 n, int32 m, int32 h, int32 *meshVertexCount, int32 *meshTriangleCount, Integer3D *meshVertexArray, Triangle *meshTriangleArray)
{
	CellStorage		*deckStorage[2];

	// Allocate storage for two decks of history.
	CellStorage *precedingCellStorage = new CellStorage[n * m * 2];

	int32 vertexCount = 0;
	int32 triangleCount = 0;
	uint16 deltaMask = 0;

	for (int32 k = 0; k < h - 1; k++)
	{
		// Ping-pong between history decks.
		deckStorage[0] = &precedingCellStorage[n * m * (k & 1)];

		for (int32 j = 0; j < m - 1; j++)
		{
			for (int32 i = 0; i < n - 1; i++)
			{
				ProcessCell(field, n, m, i, j, k, deckStorage, deltaMask, vertexCount, triangleCount, meshVertexArray, meshTriangleArray);
				deltaMask |= 1;  					// Allow reuse in x direction.
			}

			deltaMask = (deltaMask | 2) & 6;		// Allow reuse in y direction, but not x.
		}

		deckStorage[1] = deckStorage[0];			// Current deck becomes preceding deck.
		deltaMask = 4;								// Allow reuse only in z direction.
	}

	delete[] precedingCellStorage;

	*meshVertexCount = vertexCount;
	*meshTriangleCount = triangleCount;
}