From 4a42535f13fda057a77fdc0766b90e9de00b942f Mon Sep 17 00:00:00 2001
From: David Williams <david@volumesoffun.com>
Date: Wed, 23 Jul 2014 23:35:46 +0200
Subject: [PATCH 1/7] Added functions to perform 'octahedral' encoding and
 decoding of normals.

See http://jcgt.org/published/0003/02/01/paper-lowres.pdf
---
 examples/Basic/main.cpp                       | 42 +++++++++----
 examples/common/example.frag                  | 38 ++++++------
 examples/common/example.vert                  | 45 ++++++++------
 .../MarchingCubesSurfaceExtractor.h           | 62 ++++++++++++++++++-
 4 files changed, 135 insertions(+), 52 deletions(-)

diff --git a/examples/Basic/main.cpp b/examples/Basic/main.cpp
index a13da717..1f4b210e 100644
--- a/examples/Basic/main.cpp
+++ b/examples/Basic/main.cpp
@@ -50,14 +50,34 @@ void createSphereInVolume(SimpleVolume<uint8_t>& volData, float fRadius)
 				//And compute how far the current position is from the center of the volume
 				float fDistToCenter = (v3dCurrentPos - v3dVolCenter).length();
 
-				uint8_t uVoxelValue = 0;
+				// We actually want our volume to have high values in the center and low values as we move out, because our
+				// eath should be a solid sphere surrounded by empty space. If we invert the distance then this is a step in
+				// the right direction. We still have zero in the center, but lower (negative) values as we move out.
+				float density = -fDistToCenter;
 
-				//If the current voxel is less than 'radius' units from the center then we make it solid.
-				if(fDistToCenter <= fRadius)
-				{
-					//Our new voxel value
-					uVoxelValue = 255;
-				}
+				// By adding the 'planetRadius' we now have a function which starts at 'planetRadius' and still decreases as it
+				// moves out. The function passes through zero at a distance of 'planetRadius' and then continues do decrease
+				// as it gets even further out.
+				density += fRadius;
+
+				// Ideally we would like our final density value to be '255' for voxels inside the planet and '0' for voxels
+				// outside the planet. At the surface there should be a transition but this should occur not too quickly and
+				// not too slowly, as both of these will result in a jagged appearance to the mesh.
+				//
+				// We probably want the transition to occur over a few voxels, whereas it currently occurs over 255 voxels
+				// because it was derived from the distance. By scaling the density field we effectivly compress the rate
+				// at which it changes at the surface. We also make the center much too high and the outside very low, but
+				// we will clamp these to the corect range later.
+				//
+				// Note: You can try commenting out or changing the value on this line to see the effect it has.
+				density *= 50;
+
+				// Until now we've been defining our density field as if the threshold was at zero, with positive densities
+				// being solid and negative densities being empty. But actually Cubiquity operates on the range 0 to 255, and
+				// uses a threashold of 127 to decide where to place the generated surface.  Therefore we shift and clamp our
+				// density value and store it in a byte.
+				density += 127;
+				uint8_t uVoxelValue = (uint8_t)(clamp(density, 0.0f, 255.0f));
 
 				//Wrte the voxel value into the volume	
 				volData.setVoxelAt(x, y, z, uVoxelValue);
@@ -74,12 +94,12 @@ int main(int argc, char *argv[])
 	openGLWidget.show();
 
 	//Create an empty volume and then place a sphere in it
-	SimpleVolume<uint8_t> volData(PolyVox::Region(Vector3DInt32(0,0,0), Vector3DInt32(63, 63, 63)));
-	createSphereInVolume(volData, 30);
+	SimpleVolume<uint8_t> volData(PolyVox::Region(Vector3DInt32(0,0,0), Vector3DInt32(31, 31, 31)));
+	createSphereInVolume(volData, 15);
 
 	// Extract the surface for the specified region of the volume. Uncomment the line for the kind of surface extraction you want to see.
-	auto mesh = extractCubicMesh(&volData, volData.getEnclosingRegion());
-	//auto mesh = extractMarchingCubesMesh(&volData, volData.getEnclosingRegion());
+	//auto mesh = extractCubicMesh(&volData, volData.getEnclosingRegion());
+	auto mesh = extractMarchingCubesMesh(&volData, volData.getEnclosingRegion());
 
 	// The surface extractor outputs the mesh in an efficient compressed format which is not directly suitable for rendering. The easiest approach is to 
 	// decode this on the CPU as shown below, though more advanced applications can upload the compressed mesh to the GPU and decompress in shader code.
diff --git a/examples/common/example.frag b/examples/common/example.frag
index c9cbf215..56e7d9d6 100644
--- a/examples/common/example.frag
+++ b/examples/common/example.frag
@@ -1,19 +1,19 @@
-#version 130
-
-// Passed in from the vertex shader
-in vec4 worldPosition;
-in vec4 worldNormal;
-
-// the color that gets written to the display
-out vec4 outputColor;
-
-void main()
-{
-	// Again, for the purposes of these examples we cannot be sure that per-vertex normals are provided. A sensible fallback 
-	// is to use this little trick to compute per-fragment flat-shaded normals from the world positions using derivative operations.
-	vec3 normal = normalize(cross(dFdy(worldPosition.xyz), dFdx(worldPosition.xyz)));
-	
-	// We are just using the normal as the output color, and making it lighter so it looks a bit nicer. 
-	// Obviously a real shader would also do texuring, lighting, or whatever is required for the application.
-	outputColor = vec4(abs(normal) * 0.5 + vec3(0.5, 0.5, 0.5), 1.0);
-}
+#version 130
+
+// Passed in from the vertex shader
+in vec4 worldPosition;
+in vec3 worldNormal;
+
+// the color that gets written to the display
+out vec4 outputColor;
+
+void main()
+{
+	// Again, for the purposes of these examples we cannot be sure that per-vertex normals are provided. A sensible fallback 
+	// is to use this little trick to compute per-fragment flat-shaded normals from the world positions using derivative operations.
+	vec3 normal = normalize(cross(dFdy(worldPosition.xyz), dFdx(worldPosition.xyz)));
+	
+	// We are just using the normal as the output color, and making it lighter so it looks a bit nicer. 
+	// Obviously a real shader would also do texuring, lighting, or whatever is required for the application.
+	outputColor = vec4(abs(worldNormal.xyz), 1.0);
+}
diff --git a/examples/common/example.vert b/examples/common/example.vert
index 94d35f44..86f60296 100644
--- a/examples/common/example.vert
+++ b/examples/common/example.vert
@@ -1,20 +1,25 @@
-#version 140
-
-in vec4 position; // This will be the position of the vertex in model-space
-
-// The usual matrices are provided
-uniform mat4 cameraToClipMatrix;
-uniform mat4 worldToCameraMatrix;
-uniform mat4 modelToWorldMatrix;
-
-// This will be used by the fragment shader to calculate flat-shaded normals. This is an unconventional approach
-// but we use it in this example framework because not all surface extractor generate surface normals.
-out vec4 worldPosition;
-
-void main()
-{
-	// Standard sequence of OpenGL transformations.
-	worldPosition = modelToWorldMatrix * position;
-	vec4 cameraPosition = worldToCameraMatrix * worldPosition;
-	gl_Position = cameraToClipMatrix * cameraPosition;
-}
+#version 140
+
+in vec4 position; // This will be the position of the vertex in model-space
+in vec3 normal;
+
+// The usual matrices are provided
+uniform mat4 cameraToClipMatrix;
+uniform mat4 worldToCameraMatrix;
+uniform mat4 modelToWorldMatrix;
+
+// This will be used by the fragment shader to calculate flat-shaded normals. This is an unconventional approach
+// but we use it in this example framework because not all surface extractor generate surface normals.
+out vec4 worldPosition;
+out vec3 worldNormal;
+
+void main()
+{
+	// Standard sequence of OpenGL transformations.
+	worldPosition = modelToWorldMatrix * position;
+	vec4 cameraPosition = worldToCameraMatrix * worldPosition;
+	
+	worldNormal = normal;
+	
+	gl_Position = cameraToClipMatrix * cameraPosition;
+}
diff --git a/library/PolyVoxCore/include/PolyVoxCore/MarchingCubesSurfaceExtractor.h b/library/PolyVoxCore/include/PolyVoxCore/MarchingCubesSurfaceExtractor.h
index e927af7a..8317b4b4 100644
--- a/library/PolyVoxCore/include/PolyVoxCore/MarchingCubesSurfaceExtractor.h
+++ b/library/PolyVoxCore/include/PolyVoxCore/MarchingCubesSurfaceExtractor.h
@@ -50,7 +50,7 @@ namespace PolyVox
 		// Each component of the normal is encoded using 5 bits of this variable.
 		// The 16 bits are -xxxxxyyyyyzzzzz (note the left-most bit is currently 
 		// unused). Some extra shifting and scaling is required to make it signed.
-		uint16_t encodedNormal;
+		Vector2DFloat encodedNormal;
 
 		// User data
 		DataType data;
@@ -64,7 +64,7 @@ namespace PolyVox
 		return result;
 	}
 
-	inline uint16_t encodeNormal(const Vector3DFloat& normal)
+	/*inline uint16_t encodeNormal(const Vector3DFloat& normal)
 	{
 		Vector3DFloat v3dNormal = normal;
 		v3dNormal += Vector3DFloat(1.0f, 1.0f, 1.0f);
@@ -96,6 +96,64 @@ namespace PolyVox
 		result -= Vector3DFloat(1.0f, 1.0f, 1.0f);
 
 		return result;
+	}*/
+
+	// Returns ±1
+	float signNotZero(float v)
+	{
+		return v >= 0.0 ? +1.0 : -1.0;
+	}
+
+	Vector2DFloat signNotZero(Vector2DFloat v)
+	{
+		return Vector2DFloat((v.getX() >= 0.0) ? +1.0 : -1.0, (v.getY() >= 0.0) ? +1.0 : -1.0);
+	}
+
+	// Assume normalized input. Output is on [-1, 1] for each component.
+	Vector2DFloat float32x3_to_oct(Vector3DFloat v)
+	{
+		// Project the sphere onto the octahedron, and then onto the xy plane
+		Vector2DFloat p(v.getX(), v.getY());			
+		p = p * (1.0f / (abs(v.getX()) + abs(v.getY()) + abs(v.getZ())));
+
+		float refX = ((1.0f - abs(p.getY())) * signNotZero(p.getX()));
+		float refY = ((1.0f - abs(p.getX())) * signNotZero(p.getY()));
+
+		Vector2DFloat ref(refX, refY);
+
+		// Reflect the folds of the lower hemisphere over the diagonals
+		return (v.getZ() <= 0.0) ? ref : p;
+	}
+
+	Vector3DFloat oct_to_float32x3(Vector2DFloat e)
+	{
+		Vector3DFloat v = Vector3DFloat(e.getX(), e.getY(), 1.0 - abs(e.getX()) - abs(e.getY()));
+
+		float refX = ((1.0f - abs(v.getY())) * signNotZero(v.getX()));
+		float refY = ((1.0f - abs(v.getX())) * signNotZero(v.getY()));
+
+		Vector2DFloat ref(refX, refY);
+
+		if (v.getZ() < 0.0f)
+		{
+			//v.xy = (1.0 - abs(v.yx)) * signNotZero(v.xy);
+			v.setX(refX);
+			v.setY(refY);
+		}
+
+		v.normalise();
+
+		return v;
+	}
+
+	inline Vector2DFloat encodeNormal(const Vector3DFloat& normal)
+	{
+		return float32x3_to_oct(normal);
+	}
+
+	inline Vector3DFloat decode(const Vector2DFloat& encodedNormal)
+	{
+		return oct_to_float32x3(encodedNormal);
 	}
 
 	/// Decodes a MarchingCubesVertex by converting it into a regular Vertex which can then be directly used for rendering.

From 79c62be039d472a254683c70b1c110c53249925b Mon Sep 17 00:00:00 2001
From: David Williams <david@volumesoffun.com>
Date: Wed, 23 Jul 2014 23:50:04 +0200
Subject: [PATCH 2/7] Oct-encoded normal now stored as uint again.

---
 .../MarchingCubesSurfaceExtractor.h           | 28 +++++++++++++++----
 1 file changed, 23 insertions(+), 5 deletions(-)

diff --git a/library/PolyVoxCore/include/PolyVoxCore/MarchingCubesSurfaceExtractor.h b/library/PolyVoxCore/include/PolyVoxCore/MarchingCubesSurfaceExtractor.h
index 8317b4b4..39ea68c2 100644
--- a/library/PolyVoxCore/include/PolyVoxCore/MarchingCubesSurfaceExtractor.h
+++ b/library/PolyVoxCore/include/PolyVoxCore/MarchingCubesSurfaceExtractor.h
@@ -50,7 +50,7 @@ namespace PolyVox
 		// Each component of the normal is encoded using 5 bits of this variable.
 		// The 16 bits are -xxxxxyyyyyzzzzz (note the left-most bit is currently 
 		// unused). Some extra shifting and scaling is required to make it signed.
-		Vector2DFloat encodedNormal;
+		uint16_t encodedNormal;
 
 		// User data
 		DataType data;
@@ -146,14 +146,32 @@ namespace PolyVox
 		return v;
 	}
 
-	inline Vector2DFloat encodeNormal(const Vector3DFloat& normal)
+	inline uint16_t encodeNormal(const Vector3DFloat& normal)
 	{
-		return float32x3_to_oct(normal);
+		Vector2DFloat floatResult = float32x3_to_oct(normal);
+
+		floatResult += Vector2DFloat(1.0f, 1.0f); // To range 0.0f to 2.0f
+		floatResult *= Vector2DFloat(127.5f, 127.5f); // To range 0.0f to 255.0f
+
+		uint16_t resultX = static_cast<uint16_t>(floatResult.getX() + 0.5f);
+		uint16_t resultY = static_cast<uint16_t>(floatResult.getY() + 0.5f);
+
+		resultX &= 0xFF;
+		resultY &= 0xFF;
+
+		return (resultX << 8) | resultY;
 	}
 
-	inline Vector3DFloat decode(const Vector2DFloat& encodedNormal)
+	inline Vector3DFloat decode(const uint16_t& encodedNormal)
 	{
-		return oct_to_float32x3(encodedNormal);
+		uint16_t x = (encodedNormal >> 8) & 0xFF;
+		uint16_t y = (encodedNormal     ) & 0xFF;
+		Vector2DFloat floatNormal(x, y);
+
+		floatNormal /= Vector2DFloat(127.5f, 127.5f);
+		floatNormal -= Vector2DFloat(1.0f, 1.0f);
+
+		return oct_to_float32x3(floatNormal);
 	}
 
 	/// Decodes a MarchingCubesVertex by converting it into a regular Vertex which can then be directly used for rendering.

From 972bc3a4564d08ed00e66ef758d5a1a0d729be78 Mon Sep 17 00:00:00 2001
From: David Williams <david@volumesoffun.com>
Date: Thu, 24 Jul 2014 15:13:08 +0200
Subject: [PATCH 3/7] Restructuring some code.

---
 .../MarchingCubesSurfaceExtractor.h           | 105 +++++++++---------
 1 file changed, 52 insertions(+), 53 deletions(-)

diff --git a/library/PolyVoxCore/include/PolyVoxCore/MarchingCubesSurfaceExtractor.h b/library/PolyVoxCore/include/PolyVoxCore/MarchingCubesSurfaceExtractor.h
index 39ea68c2..f5abaf0a 100644
--- a/library/PolyVoxCore/include/PolyVoxCore/MarchingCubesSurfaceExtractor.h
+++ b/library/PolyVoxCore/include/PolyVoxCore/MarchingCubesSurfaceExtractor.h
@@ -64,65 +64,37 @@ namespace PolyVox
 		return result;
 	}
 
-	/*inline uint16_t encodeNormal(const Vector3DFloat& normal)
-	{
-		Vector3DFloat v3dNormal = normal;
-		v3dNormal += Vector3DFloat(1.0f, 1.0f, 1.0f);
-		uint16_t encodedX = static_cast<uint16_t>(roundToNearestInteger(v3dNormal.getX() * 15.5f));
-		uint16_t encodedY = static_cast<uint16_t>(roundToNearestInteger(v3dNormal.getY() * 15.5f));
-		uint16_t encodedZ = static_cast<uint16_t>(roundToNearestInteger(v3dNormal.getZ() * 15.5f));
-		POLYVOX_ASSERT(encodedX < 32, "Encoded value out of range");
-		POLYVOX_ASSERT(encodedY < 32, "Encoded value out of range");
-		POLYVOX_ASSERT(encodedZ < 32, "Encoded value out of range");
-		uint16_t encodedNormal = (encodedX << 10) | (encodedY << 5) | encodedZ;
-		return encodedNormal;
-	}
-
-	/// Decodes a normal from a MarchingCubesVertex
-	inline Vector3DFloat decode(const uint16_t encodedNormal)
-	{
-		// Get normal components in the range 0 to 31
-		uint16_t x = (encodedNormal >> 10) & 0x1F;
-		uint16_t y = (encodedNormal >> 5) & 0x1F;
-		uint16_t z = (encodedNormal) & 0x1F;
-
-		// Build the resulting vector
-		Vector3DFloat result(x, y, z);
-
-		// Convert to range 0.0 to 2.0
-		result *= (1.0f / 15.5f); // Division is compile-time constant
-
-		// Convert to range -1.0 to 1.0
-		result -= Vector3DFloat(1.0f, 1.0f, 1.0f);
-
-		return result;
-	}*/
-
 	// Returns ±1
 	float signNotZero(float v)
 	{
-		return v >= 0.0 ? +1.0 : -1.0;
-	}
-
-	Vector2DFloat signNotZero(Vector2DFloat v)
-	{
-		return Vector2DFloat((v.getX() >= 0.0) ? +1.0 : -1.0, (v.getY() >= 0.0) ? +1.0 : -1.0);
+		return v >= 0.0f ? +1.0f : -1.0f;
 	}
 
 	// Assume normalized input. Output is on [-1, 1] for each component.
 	Vector2DFloat float32x3_to_oct(Vector3DFloat v)
 	{
-		// Project the sphere onto the octahedron, and then onto the xy plane
-		Vector2DFloat p(v.getX(), v.getY());			
-		p = p * (1.0f / (abs(v.getX()) + abs(v.getY()) + abs(v.getZ())));
+		// Get the input components
+		float vx = v.getX();
+		float vy = v.getY();
+		float vz = v.getZ();
 
-		float refX = ((1.0f - abs(p.getY())) * signNotZero(p.getX()));
-		float refY = ((1.0f - abs(p.getX())) * signNotZero(p.getY()));
-
-		Vector2DFloat ref(refX, refY);
+		// Project the sphere onto the octahedron, and then onto the xy plane					
+		float px = vx * (1.0f / (abs(vx) + abs(vy) + abs(vz)));
+		float py = vy * (1.0f / (abs(vx) + abs(vy) + abs(vz)));
 
 		// Reflect the folds of the lower hemisphere over the diagonals
-		return (v.getZ() <= 0.0) ? ref : p;
+		if (vz <= 0.0f)
+		{
+			float refx = ((1.0f - abs(py)) * signNotZero(px));
+			float refy = ((1.0f - abs(px)) * signNotZero(py));
+			px = refx;
+			py = refy;
+		}
+
+		Vector2DFloat p(px, py);
+
+		// Reflect the folds of the lower hemisphere over the diagonals
+		return p;
 	}
 
 	Vector3DFloat oct_to_float32x3(Vector2DFloat e)
@@ -148,17 +120,44 @@ namespace PolyVox
 
 	inline uint16_t encodeNormal(const Vector3DFloat& normal)
 	{
-		Vector2DFloat floatResult = float32x3_to_oct(normal);
+		// The first part of this function is based off the code in Listing 1 of http://jcgt.org/published/0003/02/01/
+		// It was rewritten in C++ and is restructued for the CPU rather than the GPU.
 
-		floatResult += Vector2DFloat(1.0f, 1.0f); // To range 0.0f to 2.0f
-		floatResult *= Vector2DFloat(127.5f, 127.5f); // To range 0.0f to 255.0f
+		// Get the input components
+		float vx = normal.getX();
+		float vy = normal.getY();
+		float vz = normal.getZ();
 
-		uint16_t resultX = static_cast<uint16_t>(floatResult.getX() + 0.5f);
-		uint16_t resultY = static_cast<uint16_t>(floatResult.getY() + 0.5f);
+		// Project the sphere onto the octahedron, and then onto the xy plane					
+		float px = vx * (1.0f / (abs(vx) + abs(vy) + abs(vz)));
+		float py = vy * (1.0f / (abs(vx) + abs(vy) + abs(vz)));
 
+		// Reflect the folds of the lower hemisphere over the diagonals.
+		if (vz <= 0.0f)
+		{
+			float refx = ((1.0f - abs(py)) * (px >= 0.0f ? +1.0f : -1.0f));
+			float refy = ((1.0f - abs(px)) * (py >= 0.0f ? +1.0f : -1.0f));
+			px = refx;
+			py = refy;
+		}
+
+		// The next part was not given in the paper. We map our two
+		// floats into two bytes and store them in a single uint16_t
+
+		// Move from range [-1.0f, 1.0f] to [0.0f, 255.0f]
+		px = (px + 1.0) * 127.5f;
+		py = (py + 1.0) * 127.5f;
+
+		// Convert to uints
+		uint16_t resultX = static_cast<uint16_t>(px + 0.5f);
+		uint16_t resultY = static_cast<uint16_t>(py + 0.5f);
+
+		// Make sure only the lower bits are set. Probably
+		// not necessary but we're just being careful really.
 		resultX &= 0xFF;
 		resultY &= 0xFF;
 
+		// Contatenate the bytes and return the result.
 		return (resultX << 8) | resultY;
 	}
 

From 98e722271ee9ca4d159b8ae67f1e21af339b5bc9 Mon Sep 17 00:00:00 2001
From: David Williams <david@volumesoffun.com>
Date: Thu, 24 Jul 2014 15:40:03 +0200
Subject: [PATCH 4/7] More restructuring code.

---
 .../MarchingCubesSurfaceExtractor.h           | 84 ++++++-------------
 1 file changed, 24 insertions(+), 60 deletions(-)

diff --git a/library/PolyVoxCore/include/PolyVoxCore/MarchingCubesSurfaceExtractor.h b/library/PolyVoxCore/include/PolyVoxCore/MarchingCubesSurfaceExtractor.h
index f5abaf0a..2ef1b103 100644
--- a/library/PolyVoxCore/include/PolyVoxCore/MarchingCubesSurfaceExtractor.h
+++ b/library/PolyVoxCore/include/PolyVoxCore/MarchingCubesSurfaceExtractor.h
@@ -64,60 +64,6 @@ namespace PolyVox
 		return result;
 	}
 
-	// Returns ±1
-	float signNotZero(float v)
-	{
-		return v >= 0.0f ? +1.0f : -1.0f;
-	}
-
-	// Assume normalized input. Output is on [-1, 1] for each component.
-	Vector2DFloat float32x3_to_oct(Vector3DFloat v)
-	{
-		// Get the input components
-		float vx = v.getX();
-		float vy = v.getY();
-		float vz = v.getZ();
-
-		// Project the sphere onto the octahedron, and then onto the xy plane					
-		float px = vx * (1.0f / (abs(vx) + abs(vy) + abs(vz)));
-		float py = vy * (1.0f / (abs(vx) + abs(vy) + abs(vz)));
-
-		// Reflect the folds of the lower hemisphere over the diagonals
-		if (vz <= 0.0f)
-		{
-			float refx = ((1.0f - abs(py)) * signNotZero(px));
-			float refy = ((1.0f - abs(px)) * signNotZero(py));
-			px = refx;
-			py = refy;
-		}
-
-		Vector2DFloat p(px, py);
-
-		// Reflect the folds of the lower hemisphere over the diagonals
-		return p;
-	}
-
-	Vector3DFloat oct_to_float32x3(Vector2DFloat e)
-	{
-		Vector3DFloat v = Vector3DFloat(e.getX(), e.getY(), 1.0 - abs(e.getX()) - abs(e.getY()));
-
-		float refX = ((1.0f - abs(v.getY())) * signNotZero(v.getX()));
-		float refY = ((1.0f - abs(v.getX())) * signNotZero(v.getY()));
-
-		Vector2DFloat ref(refX, refY);
-
-		if (v.getZ() < 0.0f)
-		{
-			//v.xy = (1.0 - abs(v.yx)) * signNotZero(v.xy);
-			v.setX(refX);
-			v.setY(refY);
-		}
-
-		v.normalise();
-
-		return v;
-	}
-
 	inline uint16_t encodeNormal(const Vector3DFloat& normal)
 	{
 		// The first part of this function is based off the code in Listing 1 of http://jcgt.org/published/0003/02/01/
@@ -163,14 +109,32 @@ namespace PolyVox
 
 	inline Vector3DFloat decode(const uint16_t& encodedNormal)
 	{
-		uint16_t x = (encodedNormal >> 8) & 0xFF;
-		uint16_t y = (encodedNormal     ) & 0xFF;
-		Vector2DFloat floatNormal(x, y);
+		// Extract the two bytes from the uint16_t.
+		uint16_t ux = (encodedNormal >> 8) & 0xFF;
+		uint16_t uy = (encodedNormal     ) & 0xFF;
 
-		floatNormal /= Vector2DFloat(127.5f, 127.5f);
-		floatNormal -= Vector2DFloat(1.0f, 1.0f);
+		// Convert to floats in the range [-1.0f, +1.0f].
+		float ex = ux / 127.5f - 1.0f;
+		float ey = uy / 127.5f - 1.0f;
 
-		return oct_to_float32x3(floatNormal);
+		// Reconstruct the origninal vector. This is a C++ implementation
+		// of Listing 2 of http://jcgt.org/published/0003/02/01/
+		float vx = ex;
+		float vy = ey;
+		float vz = 1.0f - abs(ex) - abs(ey);
+
+		if (vz < 0.0f)
+		{
+			float refX = ((1.0f - abs(vy)) * (vx >= 0.0f ? +1.0f : -1.0f));
+			float refY = ((1.0f - abs(vx)) * (vy >= 0.0f ? +1.0f : -1.0f));
+			vx = refX;
+			vy = refY;
+		}
+
+		// Normalise and return the result.
+		Vector3DFloat v(vx, vy, vz);
+		v.normalise();
+		return v;
 	}
 
 	/// Decodes a MarchingCubesVertex by converting it into a regular Vertex which can then be directly used for rendering.

From bfa14a32dfb4febe0a2d407e09a14712bb4b2bf5 Mon Sep 17 00:00:00 2001
From: David Williams <david@volumesoffun.com>
Date: Thu, 24 Jul 2014 16:55:11 +0200
Subject: [PATCH 5/7] Split 'decode()' function into several variants so it's
 not so heavily overloaded.

---
 examples/Basic/main.cpp                           |  2 +-
 examples/DecodeOnGPU/main.cpp                     |  2 +-
 examples/OpenGL/main.cpp                          |  2 +-
 examples/Paging/main.cpp                          |  2 +-
 examples/SmoothLOD/main.cpp                       |  4 ++--
 .../include/PolyVoxCore/CubicSurfaceExtractor.h   |  6 +++---
 .../PolyVoxCore/MarchingCubesSurfaceExtractor.h   | 15 +++++++--------
 library/PolyVoxCore/include/PolyVoxCore/Mesh.h    |  4 ++--
 8 files changed, 18 insertions(+), 19 deletions(-)

diff --git a/examples/Basic/main.cpp b/examples/Basic/main.cpp
index 1f4b210e..0770e269 100644
--- a/examples/Basic/main.cpp
+++ b/examples/Basic/main.cpp
@@ -103,7 +103,7 @@ int main(int argc, char *argv[])
 
 	// The surface extractor outputs the mesh in an efficient compressed format which is not directly suitable for rendering. The easiest approach is to 
 	// decode this on the CPU as shown below, though more advanced applications can upload the compressed mesh to the GPU and decompress in shader code.
-	auto decodedMesh = decode(mesh);
+	auto decodedMesh = decodeMesh(mesh);
 
 	//Pass the surface to the OpenGL window
 	openGLWidget.addMesh(decodedMesh);
diff --git a/examples/DecodeOnGPU/main.cpp b/examples/DecodeOnGPU/main.cpp
index 2c0fff6d..7f2d83c8 100644
--- a/examples/DecodeOnGPU/main.cpp
+++ b/examples/DecodeOnGPU/main.cpp
@@ -151,7 +151,7 @@ int main(int argc, char *argv[])
 
 	// The surface extractor outputs the mesh in an efficient compressed format which is not directly suitable for rendering. The easiest approach is to 
 	// decode this on the CPU as shown below, though more advanced applications can upload the compressed mesh to the GPU and decompress in shader code.
-	//auto decodedMesh = decode(mesh);
+	//auto decodedMesh = decodeMesh(mesh);
 
 	//Pass the surface to the OpenGL window
 	OpenGLMeshData meshData = buildOpenGLMeshData(mesh);
diff --git a/examples/OpenGL/main.cpp b/examples/OpenGL/main.cpp
index 12d1068f..bdf700e0 100644
--- a/examples/OpenGL/main.cpp
+++ b/examples/OpenGL/main.cpp
@@ -128,7 +128,7 @@ int main(int argc, char *argv[])
 				auto mesh = extractMarchingCubesMesh(&volData, regToExtract);
 
 				// The returned mesh needs to be decoded to be appropriate for GPU rendering.
-				auto decodedMesh = decode(mesh);
+				auto decodedMesh = decodeMesh(mesh);
 
 				// Pass the surface to the OpenGL window. Note that we are also passing an offset in this multi-mesh example. This is because
 				// the surface extractors return a mesh with 'local space' positions to reduce storage requirements and precision problems.
diff --git a/examples/Paging/main.cpp b/examples/Paging/main.cpp
index ee2d2d23..58e94b5d 100644
--- a/examples/Paging/main.cpp
+++ b/examples/Paging/main.cpp
@@ -189,7 +189,7 @@ int main(int argc, char *argv[])
 	auto mesh = extractCubicMesh(&volData, reg2);
 	std::cout << "#vertices: " << mesh.getNoOfVertices() << std::endl;
 
-	auto decodedMesh = decode(mesh);
+	auto decodedMesh = decodeMesh(mesh);
 
 	//Pass the surface to the OpenGL window
 	openGLWidget.addMesh(decodedMesh);
diff --git a/examples/SmoothLOD/main.cpp b/examples/SmoothLOD/main.cpp
index 01c2fdaf..0054c5ce 100644
--- a/examples/SmoothLOD/main.cpp
+++ b/examples/SmoothLOD/main.cpp
@@ -92,12 +92,12 @@ int main(int argc, char *argv[])
 	//Extract the surface
 	auto meshLowLOD = extractMarchingCubesMesh(&volDataLowLOD, volDataLowLOD.getEnclosingRegion());
 	// The returned mesh needs to be decoded to be appropriate for GPU rendering.
-	auto decodedMeshLowLOD = decode(meshLowLOD);
+	auto decodedMeshLowLOD = decodeMesh(meshLowLOD);
 
 	//Extract the surface
 	auto meshHighLOD = extractMarchingCubesMesh(&volData, PolyVox::Region(Vector3DInt32(30, 0, 0), Vector3DInt32(63, 63, 63)));
 	// The returned mesh needs to be decoded to be appropriate for GPU rendering.
-	auto decodedMeshHighLOD = decode(meshHighLOD);
+	auto decodedMeshHighLOD = decodeMesh(meshHighLOD);
 
 	//Pass the surface to the OpenGL window
 	openGLWidget.addMesh(decodedMeshHighLOD, Vector3DInt32(30, 0, 0));
diff --git a/library/PolyVoxCore/include/PolyVoxCore/CubicSurfaceExtractor.h b/library/PolyVoxCore/include/PolyVoxCore/CubicSurfaceExtractor.h
index 71571038..5e059e10 100644
--- a/library/PolyVoxCore/include/PolyVoxCore/CubicSurfaceExtractor.h
+++ b/library/PolyVoxCore/include/PolyVoxCore/CubicSurfaceExtractor.h
@@ -54,7 +54,7 @@ namespace PolyVox
 	};
 
 	/// Decodes a position from a CubicVertex
-	inline Vector3DFloat decode(const Vector3DUint8& encodedPosition)
+	inline Vector3DFloat decodePosition(const Vector3DUint8& encodedPosition)
 	{
 		Vector3DFloat result(encodedPosition.getX(), encodedPosition.getY(), encodedPosition.getZ());
 		result -= 0.5f; // Apply the required offset
@@ -63,10 +63,10 @@ namespace PolyVox
 
 	/// Decodes a MarchingCubesVertex by converting it into a regular Vertex which can then be directly used for rendering.
 	template<typename DataType>
-	Vertex<DataType> decode(const CubicVertex<DataType>& cubicVertex)
+	Vertex<DataType> decodeVertex(const CubicVertex<DataType>& cubicVertex)
 	{
 		Vertex<DataType> result;
-		result.position = decode(cubicVertex.encodedPosition);
+		result.position = decodePosition(cubicVertex.encodedPosition);
 		result.normal.setElements(0.0f, 0.0f, 0.0f); // Currently not calculated
 		result.data = cubicVertex.data; // Data is not encoded
 		return result;
diff --git a/library/PolyVoxCore/include/PolyVoxCore/MarchingCubesSurfaceExtractor.h b/library/PolyVoxCore/include/PolyVoxCore/MarchingCubesSurfaceExtractor.h
index 2ef1b103..18d385f1 100644
--- a/library/PolyVoxCore/include/PolyVoxCore/MarchingCubesSurfaceExtractor.h
+++ b/library/PolyVoxCore/include/PolyVoxCore/MarchingCubesSurfaceExtractor.h
@@ -47,9 +47,8 @@ namespace PolyVox
 		// Each component of the position is stored using 8.8 fixed-point encoding.
 		Vector3DUint16 encodedPosition;
 
-		// Each component of the normal is encoded using 5 bits of this variable.
-		// The 16 bits are -xxxxxyyyyyzzzzz (note the left-most bit is currently 
-		// unused). Some extra shifting and scaling is required to make it signed.
+		// The normal is encoded as a 16-bit unsigned integer using the 'oct16'
+		// encoding described here: http://jcgt.org/published/0003/02/01/
 		uint16_t encodedNormal;
 
 		// User data
@@ -57,7 +56,7 @@ namespace PolyVox
 	};
 
 	/// Decodes a position from a MarchingCubesVertex
-	inline Vector3DFloat decode(const Vector3DUint16& encodedPosition)
+	inline Vector3DFloat decodePosition(const Vector3DUint16& encodedPosition)
 	{
 		Vector3DFloat result(encodedPosition.getX(), encodedPosition.getY(), encodedPosition.getZ());
 		result *= (1.0f / 256.0f); // Division is compile-time constant
@@ -107,7 +106,7 @@ namespace PolyVox
 		return (resultX << 8) | resultY;
 	}
 
-	inline Vector3DFloat decode(const uint16_t& encodedNormal)
+	inline Vector3DFloat decodeNormal(const uint16_t& encodedNormal)
 	{
 		// Extract the two bytes from the uint16_t.
 		uint16_t ux = (encodedNormal >> 8) & 0xFF;
@@ -139,11 +138,11 @@ namespace PolyVox
 
 	/// Decodes a MarchingCubesVertex by converting it into a regular Vertex which can then be directly used for rendering.
 	template<typename DataType>
-	Vertex<DataType> decode(const MarchingCubesVertex<DataType>& marchingCubesVertex)
+	Vertex<DataType> decodeVertex(const MarchingCubesVertex<DataType>& marchingCubesVertex)
 	{
 		Vertex<DataType> result;
-		result.position = decode(marchingCubesVertex.encodedPosition);
-		result.normal = decode(marchingCubesVertex.encodedNormal);
+		result.position = decodePosition(marchingCubesVertex.encodedPosition);
+		result.normal = decodeNormal(marchingCubesVertex.encodedNormal);
 		result.data = marchingCubesVertex.data; // Data is not encoded
 		return result;
 	}
diff --git a/library/PolyVoxCore/include/PolyVoxCore/Mesh.h b/library/PolyVoxCore/include/PolyVoxCore/Mesh.h
index a1c23b50..64356785 100644
--- a/library/PolyVoxCore/include/PolyVoxCore/Mesh.h
+++ b/library/PolyVoxCore/include/PolyVoxCore/Mesh.h
@@ -72,14 +72,14 @@ namespace PolyVox
 	};
 
 	template <typename MeshType>
-	Mesh< Vertex< typename MeshType::VertexType::DataType >, typename MeshType::IndexType > decode(const MeshType& mesh)
+	Mesh< Vertex< typename MeshType::VertexType::DataType >, typename MeshType::IndexType > decodeMesh(const MeshType& mesh)
 	{
 		Mesh< Vertex< typename MeshType::VertexType::DataType >, typename MeshType::IndexType > result;
 		result.m_vecVertices.resize(mesh.m_vecVertices.size());
 
 		for(typename MeshType::IndexType ct = 0; ct < mesh.m_vecVertices.size(); ct++)
 		{
-			result.m_vecVertices[ct] = decode(mesh.m_vecVertices[ct]);
+			result.m_vecVertices[ct] = decodeVertex(mesh.m_vecVertices[ct]);
 		}
 
 		result.m_vecTriangleIndices = mesh.m_vecTriangleIndices;

From 2b7ef5b9664fbeab221ba6c32e7043580d1a596d Mon Sep 17 00:00:00 2001
From: David Williams <david@volumesoffun.com>
Date: Thu, 24 Jul 2014 22:21:53 +0200
Subject: [PATCH 6/7] Reverted some changes which were just for testing the new
 normal encoding.

---
 examples/Basic/main.cpp      | 42 +++++++++------------------------
 examples/common/example.frag | 38 +++++++++++++++---------------
 examples/common/example.vert | 45 ++++++++++++++++--------------------
 3 files changed, 50 insertions(+), 75 deletions(-)

diff --git a/examples/Basic/main.cpp b/examples/Basic/main.cpp
index 0770e269..b697b189 100644
--- a/examples/Basic/main.cpp
+++ b/examples/Basic/main.cpp
@@ -50,34 +50,14 @@ void createSphereInVolume(SimpleVolume<uint8_t>& volData, float fRadius)
 				//And compute how far the current position is from the center of the volume
 				float fDistToCenter = (v3dCurrentPos - v3dVolCenter).length();
 
-				// We actually want our volume to have high values in the center and low values as we move out, because our
-				// eath should be a solid sphere surrounded by empty space. If we invert the distance then this is a step in
-				// the right direction. We still have zero in the center, but lower (negative) values as we move out.
-				float density = -fDistToCenter;
+				uint8_t uVoxelValue = 0;
 
-				// By adding the 'planetRadius' we now have a function which starts at 'planetRadius' and still decreases as it
-				// moves out. The function passes through zero at a distance of 'planetRadius' and then continues do decrease
-				// as it gets even further out.
-				density += fRadius;
-
-				// Ideally we would like our final density value to be '255' for voxels inside the planet and '0' for voxels
-				// outside the planet. At the surface there should be a transition but this should occur not too quickly and
-				// not too slowly, as both of these will result in a jagged appearance to the mesh.
-				//
-				// We probably want the transition to occur over a few voxels, whereas it currently occurs over 255 voxels
-				// because it was derived from the distance. By scaling the density field we effectivly compress the rate
-				// at which it changes at the surface. We also make the center much too high and the outside very low, but
-				// we will clamp these to the corect range later.
-				//
-				// Note: You can try commenting out or changing the value on this line to see the effect it has.
-				density *= 50;
-
-				// Until now we've been defining our density field as if the threshold was at zero, with positive densities
-				// being solid and negative densities being empty. But actually Cubiquity operates on the range 0 to 255, and
-				// uses a threashold of 127 to decide where to place the generated surface.  Therefore we shift and clamp our
-				// density value and store it in a byte.
-				density += 127;
-				uint8_t uVoxelValue = (uint8_t)(clamp(density, 0.0f, 255.0f));
+				//If the current voxel is less than 'radius' units from the center then we make it solid.
+				if(fDistToCenter <= fRadius)
+				{
+					//Our new voxel value
+					uVoxelValue = 255;
+				}
 
 				//Wrte the voxel value into the volume	
 				volData.setVoxelAt(x, y, z, uVoxelValue);
@@ -94,12 +74,12 @@ int main(int argc, char *argv[])
 	openGLWidget.show();
 
 	//Create an empty volume and then place a sphere in it
-	SimpleVolume<uint8_t> volData(PolyVox::Region(Vector3DInt32(0,0,0), Vector3DInt32(31, 31, 31)));
-	createSphereInVolume(volData, 15);
+	SimpleVolume<uint8_t> volData(PolyVox::Region(Vector3DInt32(0,0,0), Vector3DInt32(63, 63, 63)));
+	createSphereInVolume(volData, 30);
 
 	// Extract the surface for the specified region of the volume. Uncomment the line for the kind of surface extraction you want to see.
-	//auto mesh = extractCubicMesh(&volData, volData.getEnclosingRegion());
-	auto mesh = extractMarchingCubesMesh(&volData, volData.getEnclosingRegion());
+	auto mesh = extractCubicMesh(&volData, volData.getEnclosingRegion());
+	//auto mesh = extractMarchingCubesMesh(&volData, volData.getEnclosingRegion());
 
 	// The surface extractor outputs the mesh in an efficient compressed format which is not directly suitable for rendering. The easiest approach is to 
 	// decode this on the CPU as shown below, though more advanced applications can upload the compressed mesh to the GPU and decompress in shader code.
diff --git a/examples/common/example.frag b/examples/common/example.frag
index 56e7d9d6..c9cbf215 100644
--- a/examples/common/example.frag
+++ b/examples/common/example.frag
@@ -1,19 +1,19 @@
-#version 130
-
-// Passed in from the vertex shader
-in vec4 worldPosition;
-in vec3 worldNormal;
-
-// the color that gets written to the display
-out vec4 outputColor;
-
-void main()
-{
-	// Again, for the purposes of these examples we cannot be sure that per-vertex normals are provided. A sensible fallback 
-	// is to use this little trick to compute per-fragment flat-shaded normals from the world positions using derivative operations.
-	vec3 normal = normalize(cross(dFdy(worldPosition.xyz), dFdx(worldPosition.xyz)));
-	
-	// We are just using the normal as the output color, and making it lighter so it looks a bit nicer. 
-	// Obviously a real shader would also do texuring, lighting, or whatever is required for the application.
-	outputColor = vec4(abs(worldNormal.xyz), 1.0);
-}
+#version 130
+
+// Passed in from the vertex shader
+in vec4 worldPosition;
+in vec4 worldNormal;
+
+// the color that gets written to the display
+out vec4 outputColor;
+
+void main()
+{
+	// Again, for the purposes of these examples we cannot be sure that per-vertex normals are provided. A sensible fallback 
+	// is to use this little trick to compute per-fragment flat-shaded normals from the world positions using derivative operations.
+	vec3 normal = normalize(cross(dFdy(worldPosition.xyz), dFdx(worldPosition.xyz)));
+	
+	// We are just using the normal as the output color, and making it lighter so it looks a bit nicer. 
+	// Obviously a real shader would also do texuring, lighting, or whatever is required for the application.
+	outputColor = vec4(abs(normal) * 0.5 + vec3(0.5, 0.5, 0.5), 1.0);
+}
diff --git a/examples/common/example.vert b/examples/common/example.vert
index 86f60296..94d35f44 100644
--- a/examples/common/example.vert
+++ b/examples/common/example.vert
@@ -1,25 +1,20 @@
-#version 140
-
-in vec4 position; // This will be the position of the vertex in model-space
-in vec3 normal;
-
-// The usual matrices are provided
-uniform mat4 cameraToClipMatrix;
-uniform mat4 worldToCameraMatrix;
-uniform mat4 modelToWorldMatrix;
-
-// This will be used by the fragment shader to calculate flat-shaded normals. This is an unconventional approach
-// but we use it in this example framework because not all surface extractor generate surface normals.
-out vec4 worldPosition;
-out vec3 worldNormal;
-
-void main()
-{
-	// Standard sequence of OpenGL transformations.
-	worldPosition = modelToWorldMatrix * position;
-	vec4 cameraPosition = worldToCameraMatrix * worldPosition;
-	
-	worldNormal = normal;
-	
-	gl_Position = cameraToClipMatrix * cameraPosition;
-}
+#version 140
+
+in vec4 position; // This will be the position of the vertex in model-space
+
+// The usual matrices are provided
+uniform mat4 cameraToClipMatrix;
+uniform mat4 worldToCameraMatrix;
+uniform mat4 modelToWorldMatrix;
+
+// This will be used by the fragment shader to calculate flat-shaded normals. This is an unconventional approach
+// but we use it in this example framework because not all surface extractor generate surface normals.
+out vec4 worldPosition;
+
+void main()
+{
+	// Standard sequence of OpenGL transformations.
+	worldPosition = modelToWorldMatrix * position;
+	vec4 cameraPosition = worldToCameraMatrix * worldPosition;
+	gl_Position = cameraToClipMatrix * cameraPosition;
+}

From 112b277452f090afb649978e264f23204e74d2dc Mon Sep 17 00:00:00 2001
From: David Williams <david@volumesoffun.com>
Date: Thu, 24 Jul 2014 22:38:21 +0200
Subject: [PATCH 7/7] Updated GLSL decoder for new normal packing scheme.

---
 examples/DecodeOnGPU/decode.vert | 25 +++++++++++++++++++------
 1 file changed, 19 insertions(+), 6 deletions(-)

diff --git a/examples/DecodeOnGPU/decode.vert b/examples/DecodeOnGPU/decode.vert
index 5164468a..1ca80d09 100644
--- a/examples/DecodeOnGPU/decode.vert
+++ b/examples/DecodeOnGPU/decode.vert
@@ -13,17 +13,30 @@ uniform mat4 modelToWorldMatrix;
 out vec4 worldPosition;
 out vec4 worldNormal;
 
+// Returns +/- 1
+vec2 signNotZero(vec2 v)
+{
+	return vec2((v.x >= 0.0) ? +1.0 : -1.0, (v.y >= 0.0) ? +1.0 : -1.0);
+}
+
 void main()
 {
 	vec4 decodedPosition = position;
 	decodedPosition.xyz = decodedPosition.xyz * (1.0 / 256.0);
 	
-	uint encodedX = (normal >> 10u) & 0x1Fu;
-	uint encodedY = (normal >> 5u) & 0x1Fu;
-	uint encodedZ = (normal) & 0x1Fu;
-	worldNormal.xyz = vec3(encodedX, encodedY, encodedZ);
-	worldNormal.xyz = worldNormal.xyz / 15.5;
-	worldNormal.xyz = worldNormal.xyz - vec3(1.0, 1.0, 1.0);
+	//Get the encoded bytes of the normal
+	uint encodedX = (normal >> 8u) & 0xFFu;
+	uint encodedY = (normal) & 0xFFu;
+	
+	// Map to range [-1.0, +1.0]
+	vec2 e = vec2(encodedX, encodedY);
+	e = e * vec2(1.0 / 127.5, 1.0 / 127.5);
+	e = e - vec2(1.0, 1.0);
+	
+	// Now decode normal using listing 2 of http://jcgt.org/published/0003/02/01/
+	vec3 v = vec3(e.xy, 1.0 - abs(e.x) - abs(e.y));
+	if (v.z < 0) v.xy = (1.0 - abs(v.yx)) * signNotZero(v.xy);	
+	worldNormal.xyz = normalize(v);
 	worldNormal.w = 1.0;
 	
 	// Standard sequence of OpenGL transformations.