diff --git a/library/PolyVoxCore/CMakeLists.txt b/library/PolyVoxCore/CMakeLists.txt
index 18097abe..f483964c 100644
--- a/library/PolyVoxCore/CMakeLists.txt
+++ b/library/PolyVoxCore/CMakeLists.txt
@@ -36,6 +36,8 @@ SET(CORE_INC_FILES
 	include/Material.h
 	include/MaterialDensityPair.h
 	include/Mesh.h
+	include/MeshDecimator.h
+	include/MeshDecimator.inl
 	include/MeshEdge.h
 	include/MeshFace.h
 	include/MeshVertex.h
diff --git a/library/PolyVoxCore/include/MeshDecimator.h b/library/PolyVoxCore/include/MeshDecimator.h
index 2de5428e..7c422b61 100644
--- a/library/PolyVoxCore/include/MeshDecimator.h
+++ b/library/PolyVoxCore/include/MeshDecimator.h
@@ -24,8 +24,38 @@ freely, subject to the following restrictions:
 #ifndef __PolyVox_MeshDecimator_H__
 #define __PolyVox_MeshDecimator_H__
 
+#include <bitset>
+
 namespace PolyVox
 {
+	enum POLYVOXCORE_API NormalFlags
+	{
+		NF_NORMAL_NEG_X,
+		NF_NORMAL_POS_X,
+		NF_NORMAL_NEG_Y,
+		NF_NORMAL_POS_Y,
+		NF_NORMAL_NEG_Z,
+		NF_NORMAL_POS_Z,
+		NF_NO_OF_FLAGS
+	};
+
+	enum Stages
+	{
+		STAGE_FACE,
+		STAGE_EDGE,
+		NO_OF_STAGES
+	};
+
+	struct VertexMetadata
+	{
+		bool hasDuplicate;
+		uint64_t materialKey;
+		list<uint32_t> trianglesUsingVertex;
+		int noOfDifferentNormals;
+		Vector3DFloat normal;
+		std::bitset<NF_NO_OF_FLAGS> m_bNormalFlags;
+	};
+
 	template <typename VertexType>
 	class MeshDecimator
 	{
@@ -35,12 +65,18 @@ namespace PolyVox
 		void execute();
 
 	private:
+
+		void fillVertexMetadata(std::vector<VertexMetadata>& vecVertexMetadata);
+
+		int countZeros(void);
+
 		SurfaceMesh<VertexType>* m_pInputMesh;
 		//SurfaceMesh<PositionMaterial>* pMeshOutput;
 
 		void countNoOfNeighboursUsingMaterial(void);
 		uint32_t performDecimationPass(float fMinDotProductForCollapse);
 		bool isSubset(std::bitset<VF_NO_OF_FLAGS> a, std::bitset<VF_NO_OF_FLAGS> b);
+		bool isSubsetCubic(std::bitset<NF_NO_OF_FLAGS> a, std::bitset<NF_NO_OF_FLAGS> b);
 
 		bool canCollapseEdge(uint32_t uSrc, uint32_t uDest);
 
@@ -50,14 +86,21 @@ namespace PolyVox
 		vector<bool> vertexLocked;
 		vector<uint32_t> vertexMapper;
 
-		vector< list<uint32_t> > trianglesUsingVertex;
-		vector<int> noOfDifferentNormals;
+		//vector< list<uint32_t> > trianglesUsingVertexCurrently;
 
-		vector<uint64_t> materialKey;
+		vector<int> vecOfTriCts;
+		vector<Vector3DFloat> vecOfTriNormals;
 
-		vector<bool> hasDuplicate;
+		int m_NoOfZeros;
+
+		//vector<int> noOfDifferentNormals;
+
+		vector<VertexMetadata> m_vecInitialVertexMetadata;
+		vector<VertexMetadata> m_vecCurrentVertexMetadata;
 
 		float fMinDotProductForCollapse;
+
+		uint8_t m_uStage;
 	};
 }
 
diff --git a/library/PolyVoxCore/include/MeshDecimator.inl b/library/PolyVoxCore/include/MeshDecimator.inl
index 367b437b..d87e6a50 100644
--- a/library/PolyVoxCore/include/MeshDecimator.inl
+++ b/library/PolyVoxCore/include/MeshDecimator.inl
@@ -38,12 +38,15 @@ namespace PolyVox
 		// determine when material boundary edges can collapse.
 		countNoOfNeighboursUsingMaterial();
 
+		fillVertexMetadata(m_vecInitialVertexMetadata);
+
+		m_uStage = STAGE_FACE;
 		uint32_t noOfEdgesCollapsed;
 		do
 		{
 			noOfEdgesCollapsed = performDecimationPass(fMinDotProductForCollapse);
 			m_pInputMesh->removeDegenerateTris();	
-			m_pInputMesh->removeUnusedVertices();
+			//m_pInputMesh->removeUnusedVertices();
 		}while(noOfEdgesCollapsed > 0);
 
 		//Decimation will have invalidated LOD levels.
@@ -55,58 +58,144 @@ namespace PolyVox
 	}
 
 	template <typename VertexType>
-	uint32_t MeshDecimator<VertexType>::performDecimationPass(float fMinDotProductForCollapse)
+	void MeshDecimator<VertexType>::fillVertexMetadata(std::vector<VertexMetadata>& vecVertexMetadata)
 	{
-		hasDuplicate.clear();
+		vecVertexMetadata.clear();
+		vecVertexMetadata.resize(m_pInputMesh->m_vecVertices.size());
+		//Initialise the metadata
+		for(int ct = 0; ct < vecVertexMetadata.size(); ct++)
+		{
+			vecVertexMetadata[ct].hasDuplicate = false;
+			vecVertexMetadata[ct].materialKey = 0;
+			vecVertexMetadata[ct].trianglesUsingVertex.clear();
+			vecVertexMetadata[ct].noOfDifferentNormals = 0;
+		}
+
+		//Determine triangles using each vertex
+		/*trianglesUsingVertex.clear();
+		trianglesUsingVertex.resize(m_pInputMesh->m_vecVertices.size());*/
+		for(int ct = 0; ct < m_pInputMesh->m_vecTriangleIndices.size(); ct++)
+		{
+			int triangle = ct / 3;
+
+			vecVertexMetadata[m_pInputMesh->m_vecTriangleIndices[ct]].trianglesUsingVertex.push_back(triangle);
+		}
+
+		/*hasDuplicate.clear();
 		hasDuplicate.resize(m_pInputMesh->m_vecVertices.size());
-		std::fill(hasDuplicate.begin(), hasDuplicate.end(), false);
+		std::fill(hasDuplicate.begin(), hasDuplicate.end(), false);*/
 		for(int outerCt = 0; outerCt < m_pInputMesh->m_vecVertices.size()-1; outerCt++)
 		{
 			for(int innerCt = outerCt+1; innerCt < m_pInputMesh->m_vecVertices.size(); innerCt++)
 			{
 				if((m_pInputMesh->m_vecVertices[innerCt].position - m_pInputMesh->m_vecVertices[outerCt].position).lengthSquared() < 0.001f)
 				{
-					hasDuplicate[innerCt] = true;
-					hasDuplicate[outerCt] = true;
+					vecVertexMetadata[innerCt].hasDuplicate = true;
+					vecVertexMetadata[outerCt].hasDuplicate = true;
 				}
 			}
 		}
-
-		int duplicates = 0;
-		for(int ct = 0; ct < hasDuplicate.size(); ct++)
-		{
-			if(hasDuplicate[ct])
-			{
-				duplicates++;
-			}
-		}
-		std::cout << duplicates << std::endl;
 		
-		materialKey.clear();
+		/*materialKey.clear();
 		materialKey.resize(m_pInputMesh->m_vecVertices.size());
-		std::fill(materialKey.begin(), materialKey.end(), 0);
+		std::fill(materialKey.begin(), materialKey.end(), 0);*/
 		for(int ct = 0; ct < m_pInputMesh->m_vecTriangleIndices.size(); ct++)
 		{
 			uint32_t vertex = m_pInputMesh->m_vecTriangleIndices[ct];
 
 			//NOTE: uint8_t may not always be large engouh?
 			uint8_t uMaterial = m_pInputMesh->m_vecVertices[vertex].material;
-			materialKey[vertex] <<= 8;
-			materialKey[vertex] |= uMaterial;
+			vecVertexMetadata[vertex].materialKey <<= 8;
+			vecVertexMetadata[vertex].materialKey |= uMaterial;
 		}
 
+		// Each triangle exists in this vector once.
+		vecOfTriCts.clear();
+		vecOfTriCts.resize(m_pInputMesh->m_vecTriangleIndices.size() / 3);
+		for(int triCt = 0; triCt < vecOfTriCts.size(); triCt++)
+		{
+			vecOfTriCts[triCt] = triCt;
+		}
 
+		vecOfTriNormals.clear();
+		vecOfTriNormals.resize(vecOfTriCts.size());
+		for(int ct = 0; ct < vecOfTriCts.size(); ct++)
+		{
+			int triCt = vecOfTriCts[ct];
+			int v0 = m_pInputMesh->m_vecTriangleIndices[triCt * 3 + 0];
+			int v1 = m_pInputMesh->m_vecTriangleIndices[triCt * 3 + 1];
+			int v2 = m_pInputMesh->m_vecTriangleIndices[triCt * 3 + 2];
+
+			//Handle degenerates?
+
+			Vector3DFloat v0v1 = m_pInputMesh->m_vecVertices[v1].position - m_pInputMesh->m_vecVertices[v0].position;
+			Vector3DFloat v0v2 = m_pInputMesh->m_vecVertices[v2].position - m_pInputMesh->m_vecVertices[v0].position;
+			Vector3DFloat normal = v0v1.cross(v0v2);
+
+			normal.normalise();
+
+			vecOfTriNormals[ct] = normal;
+		}
+
+		//noOfDifferentNormals.clear();
+		//noOfDifferentNormals.resize(m_pInputMesh->m_vecVertices.size());
+		//std::fill(vecVertexMetadata.noOfDifferentNormals.begin(), vecVertexMetadata.noOfDifferentNormals.end(), 0);
+		for(int ct = 0; ct < m_pInputMesh->m_vecVertices.size(); ct++)
+		{
+			Vector3DFloat sumOfNormals(0.0f,0.0f,0.0f);
+			for(list<uint32_t>::const_iterator iter = vecVertexMetadata[ct].trianglesUsingVertex.cbegin(); iter != vecVertexMetadata[ct].trianglesUsingVertex.cend(); iter++)
+			{
+				sumOfNormals += vecOfTriNormals[*iter];
+			}
+
+			vecVertexMetadata[ct].noOfDifferentNormals = 0;
+			if(abs(sumOfNormals.getX()) > 0.001)
+				vecVertexMetadata[ct].noOfDifferentNormals++;
+			if(abs(sumOfNormals.getY()) > 0.001)
+				vecVertexMetadata[ct].noOfDifferentNormals++;
+			if(abs(sumOfNormals.getZ()) > 0.001)
+				vecVertexMetadata[ct].noOfDifferentNormals++;
+
+			if(sumOfNormals.getX() < -0.001)
+				vecVertexMetadata[ct].m_bNormalFlags.set(NF_NORMAL_NEG_X);
+			if(sumOfNormals.getX() > 0.001)
+				vecVertexMetadata[ct].m_bNormalFlags.set(NF_NORMAL_POS_X);
+			if(sumOfNormals.getY() < -0.001)
+				vecVertexMetadata[ct].m_bNormalFlags.set(NF_NORMAL_NEG_Y);
+			if(sumOfNormals.getY() > 0.001)
+				vecVertexMetadata[ct].m_bNormalFlags.set(NF_NORMAL_POS_Y);
+			if(sumOfNormals.getZ() < -0.001)
+				vecVertexMetadata[ct].m_bNormalFlags.set(NF_NORMAL_NEG_Z);
+			if(sumOfNormals.getZ() > 0.001)
+				vecVertexMetadata[ct].m_bNormalFlags.set(NF_NORMAL_POS_Z);
+
+			vecVertexMetadata[ct].normal = sumOfNormals;
+			vecVertexMetadata[ct].normal.normalise();
+
+
+			/*if(vecVertexMetadata[ct].noOfDifferentNormals == 3)
+				std::cout << "CORNER at " << ct <<std::endl;*/
+
+			//assert(noOfDifferentNormals[ct] > 0);
+		}
+
+		//std::cout << "----------" <<std::endl;
+	}
+
+	template <typename VertexType>
+	uint32_t MeshDecimator<VertexType>::performDecimationPass(float fMinDotProductForCollapse)
+	{
 		// I'm using a vector of lists here, rather than a vector of sets,
 		// because I don't believe that duplicates should occur. But this
 		// might be worth checking if we have problems in the future.
-		trianglesUsingVertex.clear();
-		trianglesUsingVertex.resize(m_pInputMesh->m_vecVertices.size());
+		/*trianglesUsingVertexCurrently.clear();
+		trianglesUsingVertexCurrently.resize(m_pInputMesh->m_vecVertices.size());
 		for(int ct = 0; ct < m_pInputMesh->m_vecTriangleIndices.size(); ct++)
 		{
 			int triangle = ct / 3;
 
-			trianglesUsingVertex[m_pInputMesh->m_vecTriangleIndices[ct]].push_back(triangle);
-		}
+			trianglesUsingVertexCurrently[m_pInputMesh->m_vecTriangleIndices[ct]].push_back(triangle);
+		}*/
 
 		// Count how many edges we have collapsed
 		uint32_t noOfEdgesCollapsed = 0;
@@ -131,59 +220,7 @@ namespace PolyVox
 			vertexLocked[ct] = false;
 		}
 
-		
-		// Each triangle exists in this vector once.
-		vector<int> vecOfTriCts(m_pInputMesh->m_vecTriangleIndices.size() / 3);
-		for(int triCt = 0; triCt < vecOfTriCts.size(); triCt++)
-		{
-			vecOfTriCts[triCt] = triCt;
-		}
-
-		// It *may* be beneficial to randomise the order in which triangles
-		// are processed to get a more uniform distribution off collapses and 
-		// more equally sized triangles at the end. This need more testing really.
-		//random_shuffle(vecOfTriCts.begin(), vecOfTriCts.end());
-
-		vector<Vector3DFloat> vecOfTriNormals(vecOfTriCts.size());
-		for(int ct = 0; ct < vecOfTriCts.size(); ct++)
-		{
-			int triCt = vecOfTriCts[ct];
-			int v0 = m_pInputMesh->m_vecTriangleIndices[triCt * 3 + 0];
-			int v1 = m_pInputMesh->m_vecTriangleIndices[triCt * 3 + 1];
-			int v2 = m_pInputMesh->m_vecTriangleIndices[triCt * 3 + 2];
-
-			//Handle degenerates?
-
-			Vector3DFloat v0v1 = m_pInputMesh->m_vecVertices[v1].position - m_pInputMesh->m_vecVertices[v0].position;
-			Vector3DFloat v0v2 = m_pInputMesh->m_vecVertices[v2].position - m_pInputMesh->m_vecVertices[v0].position;
-			Vector3DFloat normal = v0v1.cross(v0v2);
-
-			vecOfTriNormals[ct] = normal;
-		}
-
-		noOfDifferentNormals.clear();
-		noOfDifferentNormals.resize(m_pInputMesh->m_vecVertices.size());
-		std::fill(noOfDifferentNormals.begin(), noOfDifferentNormals.end(), 0);
-		for(int ct = 0; ct < m_pInputMesh->m_vecVertices.size(); ct++)
-		{
-			//const list<uint32_t>& trianglesUsingVertex = trianglesUsingVertex[ct];
-
-			Vector3DFloat sumOfNormals(0.0f,0.0f,0.0f);
-			for(list<uint32_t>::const_iterator iter = trianglesUsingVertex[ct].cbegin(); iter != trianglesUsingVertex[ct].cend(); iter++)
-			{
-				sumOfNormals += vecOfTriNormals[*iter];
-			}
-
-			noOfDifferentNormals[ct] = 0;
-			if(abs(sumOfNormals.getX()) > 0.001)
-				noOfDifferentNormals[ct]++;
-			if(abs(sumOfNormals.getY()) > 0.001)
-				noOfDifferentNormals[ct]++;
-			if(abs(sumOfNormals.getZ()) > 0.001)
-				noOfDifferentNormals[ct]++;
-
-			assert(noOfDifferentNormals[ct] > 0);
-		}
+		fillVertexMetadata(m_vecCurrentVertexMetadata);
 
 		//For each triange...
 		for(int ctIter = 0; ctIter < vecOfTriCts.size(); ctIter++)
@@ -213,6 +250,8 @@ namespace PolyVox
 			}
 		}
 
+		m_NoOfZeros = countZeros();
+
 		if(noOfEdgesCollapsed > 0)
 		{
 			//Fix up the indices
@@ -227,6 +266,8 @@ namespace PolyVox
 			}
 		}
 
+		m_NoOfZeros = countZeros();
+
 		return noOfEdgesCollapsed;
 	}
 
@@ -289,210 +330,231 @@ namespace PolyVox
 		return result;
 	}
 
+	template <typename VertexType>
+	bool MeshDecimator<VertexType>::isSubsetCubic(std::bitset<NF_NO_OF_FLAGS> a, std::bitset<NF_NO_OF_FLAGS> b)
+	{
+		bool result = true;
+
+		for(int ct = 0; ct < NF_NO_OF_FLAGS; ct++) //Start at '1' to skip material flag
+		{
+			if(a.test(ct))
+			{
+				if(b.test(ct) == false)
+				{
+					result = false;
+					break;
+				}
+			}
+		}
+
+		return result;
+	}
+
 	//template <typename VertexType>
 	bool MeshDecimator<PositionMaterialNormal>::canCollapseEdge(uint32_t v0, uint32_t v1)
 	{
 		//A vertex will be locked if it has already been involved in a collapse this pass.
-				if(vertexLocked[v0] || vertexLocked[v1])
+		if(vertexLocked[v0] || vertexLocked[v1])
+		{
+			return false;
+		}
+
+		if(m_pInputMesh->m_vecVertices[v0].getMaterial() != m_pInputMesh->m_vecVertices[v1].getMaterial())
+		{
+			return false;
+		}
+
+		//For now, don't collapse vertices on material edges...
+		if(m_pInputMesh->m_vecVertices[v0].isOnMaterialEdge() || m_pInputMesh->m_vecVertices[v1].isOnMaterialEdge())
+		{
+			if(true)
+			{
+				bool pass = false;						
+
+				bool allMatch = false;
+
+				// On the original undecimated mesh a material boundary vertex on a straight edge will
+				// have four neighbours with the same material. If it's on a corner it will have a
+				// different number. We only collapse straight edges to avoid changingthe shape of the
+				// material boundary.
+				if(m_vecNoOfNeighboursUsingMaterial[v0] == m_vecNoOfNeighboursUsingMaterial[v1])
+				{
+					if(m_vecNoOfNeighboursUsingMaterial[v0] == 4)
+					{
+						allMatch = true;
+					}
+				}
+
+				bool movementValid = false;
+				Vector3DFloat movement = m_pInputMesh->m_vecVertices[v1].getPosition() - m_pInputMesh->m_vecVertices[v0].getPosition();
+				movement.normalise();
+				if(movement.dot(Vector3DFloat(0,0,1)) > 0.999)
+				{
+					movementValid = true;
+				}
+
+				if(movement.dot(Vector3DFloat(0,1,0)) > 0.999)
+				{
+					movementValid = true;
+				}
+
+				if(movement.dot(Vector3DFloat(1,0,0)) > 0.999)
+				{
+					movementValid = true;
+				}
+
+				if(movement.dot(Vector3DFloat(0,0,-1)) > 0.999)
+				{
+					movementValid = true;
+				}
+
+				if(movement.dot(Vector3DFloat(0,-1,0)) > 0.999)
+				{
+					movementValid = true;
+				}
+
+				if(movement.dot(Vector3DFloat(-1,0,0)) > 0.999)
+				{
+					movementValid = true;
+				}
+
+				if(movementValid && allMatch)
+				{
+					pass = true;
+				}
+
+				if(!pass)
 				{
 					return false;
 				}
+			}
+			else //Material collapses not allowed
+			{
+				return false;
+			}
+		}
 
-				if(m_pInputMesh->m_vecVertices[v0].getMaterial() != m_pInputMesh->m_vecVertices[v1].getMaterial())
+		// Vertices on the geometrical edge of surface meshes need special handling. 
+		// We check for this by whether any of the edge flags are set.
+		if(m_pInputMesh->m_vecVertices[v0].m_bFlags.any() || m_pInputMesh->m_vecVertices[v1].m_bFlags.any())
+		{
+			// Assume we can't collapse until we prove otherwise...
+			bool bCollapseGeometryEdgePair = false;
+
+			// We can collapse normal vertices onto edge vertices, and edge vertices
+			// onto corner vertices, but not vice-versa. Hence we check whether all
+			// the edge flags in the source vertex are also set in the destination vertex.
+			if(isSubset(m_pInputMesh->m_vecVertices[v0].m_bFlags, m_pInputMesh->m_vecVertices[v1].m_bFlags))
+			{
+				// In general adjacent regions surface meshes may collapse differently
+				// and this can cause cracks. We solve this by only allowing the collapse
+				// is the normals are exactly the same. We do not use the user provided
+				// tolerence here (but do allow for floating point error).
+				if(m_pInputMesh->m_vecVertices[v0].getNormal().dot(m_pInputMesh->m_vecVertices[v1].getNormal()) > 0.999)
 				{
-					return false;
+					// Ok, this pair can collapse.
+					bCollapseGeometryEdgePair = true;
 				}
+			}
 
-				//For now, don't collapse vertices on material edges...
-				if(m_pInputMesh->m_vecVertices[v0].isOnMaterialEdge() || m_pInputMesh->m_vecVertices[v1].isOnMaterialEdge())
-				{
-					if(true)
-					{
-						bool pass = false;						
+			// Use the result.
+			if(!bCollapseGeometryEdgePair)
+			{
+				return false;
+			}
+		}
 
-						bool allMatch = false;
+		//Check the normals are within the threashold.
+		if(m_pInputMesh->m_vecVertices[v0].getNormal().dot(m_pInputMesh->m_vecVertices[v1].getNormal()) < fMinDotProductForCollapse)
+		{
+			return false;
+		}
 
-						// On the original undecimated mesh a material boundary vertex on a straight edge will
-						// have four neighbours with the same material. If it's on a corner it will have a
-						// different number. We only collapse straight edges to avoid changingthe shape of the
-						// material boundary.
-						if(m_vecNoOfNeighboursUsingMaterial[v0] == m_vecNoOfNeighboursUsingMaterial[v1])
-						{
-							if(m_vecNoOfNeighboursUsingMaterial[v0] == 4)
-							{
-								allMatch = true;
-							}
-						}
+		////////////////////////////////////////////////////////////////////////////////
+		//The last test is whether we will flip any of the faces
 
-						bool movementValid = false;
-						Vector3DFloat movement = m_pInputMesh->m_vecVertices[v1].getPosition() - m_pInputMesh->m_vecVertices[v0].getPosition();
-						movement.normalise();
-						if(movement.dot(Vector3DFloat(0,0,1)) > 0.999)
-						{
-							movementValid = true;
-						}
+		bool faceFlipped = false;
+		//list<uint32_t> triangles = trianglesUsingVertexCurrently[v0];
+		list<uint32_t> triangles = m_vecCurrentVertexMetadata[v0].trianglesUsingVertex;
+		/*set<uint32_t> triangles;
+		std::set_union(trianglesUsingVertex[v0].begin(), trianglesUsingVertex[v0].end(),
+			trianglesUsingVertex[v1].begin(), trianglesUsingVertex[v1].end(),
+			std::inserter(triangles, triangles.begin()));*/
 
-						if(movement.dot(Vector3DFloat(0,1,0)) > 0.999)
-						{
-							movementValid = true;
-						}
-
-						if(movement.dot(Vector3DFloat(1,0,0)) > 0.999)
-						{
-							movementValid = true;
-						}
-
-						if(movement.dot(Vector3DFloat(0,0,-1)) > 0.999)
-						{
-							movementValid = true;
-						}
-
-						if(movement.dot(Vector3DFloat(0,-1,0)) > 0.999)
-						{
-							movementValid = true;
-						}
-
-						if(movement.dot(Vector3DFloat(-1,0,0)) > 0.999)
-						{
-							movementValid = true;
-						}
-
-						if(movementValid && allMatch)
-						{
-							pass = true;
-						}
-
-						if(!pass)
-						{
-							return false;
-						}
-					}
-					else //Material collapses not allowed
-					{
-						return false;
-					}
-				}
-
-				// Vertices on the geometrical edge of surface meshes need special handling. 
-				// We check for this by whether any of the edge flags are set.
-				if(m_pInputMesh->m_vecVertices[v0].m_bFlags.any() || m_pInputMesh->m_vecVertices[v1].m_bFlags.any())
-				{
-					// Assume we can't collapse until we prove otherwise...
-					bool bCollapseGeometryEdgePair = false;
-
-					// We can collapse normal vertices onto edge vertices, and edge vertices
-					// onto corner vertices, but not vice-versa. Hence we check whether all
-					// the edge flags in the source vertex are also set in the destination vertex.
-					if(isSubset(m_pInputMesh->m_vecVertices[v0].m_bFlags, m_pInputMesh->m_vecVertices[v1].m_bFlags))
-					{
-						// In general adjacent regions surface meshes may collapse differently
-						// and this can cause cracks. We solve this by only allowing the collapse
-						// is the normals are exactly the same. We do not use the user provided
-						// tolerence here (but do allow for floating point error).
-						if(m_pInputMesh->m_vecVertices[v0].getNormal().dot(m_pInputMesh->m_vecVertices[v1].getNormal()) > 0.999)
-						{
-							// Ok, this pair can collapse.
-							bCollapseGeometryEdgePair = true;
-						}
-					}
-
-					// Use the result.
-					if(!bCollapseGeometryEdgePair)
-					{
-						return false;
-					}
-				}
-
-				//Check the normals are within the threashold.
-				if(m_pInputMesh->m_vecVertices[v0].getNormal().dot(m_pInputMesh->m_vecVertices[v1].getNormal()) < fMinDotProductForCollapse)
-				{
-					return false;
-				}
-
-				////////////////////////////////////////////////////////////////////////////////
-				//The last test is whether we will flip any of the faces
-
-				bool faceFlipped = false;
-				list<uint32_t> triangles = trianglesUsingVertex[v0];
-				/*set<uint32_t> triangles;
-				std::set_union(trianglesUsingVertex[v0].begin(), trianglesUsingVertex[v0].end(),
-					trianglesUsingVertex[v1].begin(), trianglesUsingVertex[v1].end(),
-					std::inserter(triangles, triangles.begin()));*/
-
-				for(list<uint32_t>::iterator triIter = triangles.begin(); triIter != triangles.end(); triIter++)
-				{
-					uint32_t tri = *triIter;
+		for(list<uint32_t>::iterator triIter = triangles.begin(); triIter != triangles.end(); triIter++)
+		{
+			uint32_t tri = *triIter;
 					
-					uint32_t v0Old = m_pInputMesh->m_vecTriangleIndices[tri * 3];
-					uint32_t v1Old = m_pInputMesh->m_vecTriangleIndices[tri * 3 + 1];
-					uint32_t v2Old = m_pInputMesh->m_vecTriangleIndices[tri * 3 + 2];
+			uint32_t v0Old = m_pInputMesh->m_vecTriangleIndices[tri * 3];
+			uint32_t v1Old = m_pInputMesh->m_vecTriangleIndices[tri * 3 + 1];
+			uint32_t v2Old = m_pInputMesh->m_vecTriangleIndices[tri * 3 + 2];
 
-					//Check if degenerate
-					if((v0Old == v1Old) || (v1Old == v2Old) || (v2Old == v0Old))
-					{
-						continue;
-					}
+			//Check if degenerate
+			if((v0Old == v1Old) || (v1Old == v2Old) || (v2Old == v0Old))
+			{
+				continue;
+			}
 
-					uint32_t v0New = v0Old;
-					uint32_t v1New = v1Old;
-					uint32_t v2New = v2Old;
+			uint32_t v0New = v0Old;
+			uint32_t v1New = v1Old;
+			uint32_t v2New = v2Old;
 
-					if(v0New == v0)
-						v0New = v1;
-					if(v1New == v0)
-						v1New = v1;
-					if(v2New == v0)
-						v2New = v1;
+			if(v0New == v0)
+				v0New = v1;
+			if(v1New == v0)
+				v1New = v1;
+			if(v2New == v0)
+				v2New = v1;
 
-					//Check if degenerate
-					if((v0New == v1New) || (v1New == v2New) || (v2New == v0New))
-					{
-						continue;
-					}
+			//Check if degenerate
+			if((v0New == v1New) || (v1New == v2New) || (v2New == v0New))
+			{
+				continue;
+			}
 
-					Vector3DFloat v0OldPos = m_pInputMesh->m_vecVertices[vertexMapper[v0Old]].getPosition();
-					Vector3DFloat v1OldPos = m_pInputMesh->m_vecVertices[vertexMapper[v1Old]].getPosition();
-					Vector3DFloat v2OldPos = m_pInputMesh->m_vecVertices[vertexMapper[v2Old]].getPosition();
+			Vector3DFloat v0OldPos = m_pInputMesh->m_vecVertices[vertexMapper[v0Old]].getPosition();
+			Vector3DFloat v1OldPos = m_pInputMesh->m_vecVertices[vertexMapper[v1Old]].getPosition();
+			Vector3DFloat v2OldPos = m_pInputMesh->m_vecVertices[vertexMapper[v2Old]].getPosition();
 
-					Vector3DFloat v0NewPos = m_pInputMesh->m_vecVertices[vertexMapper[v0New]].getPosition();
-					Vector3DFloat v1NewPos = m_pInputMesh->m_vecVertices[vertexMapper[v1New]].getPosition();
-					Vector3DFloat v2NewPos = m_pInputMesh->m_vecVertices[vertexMapper[v2New]].getPosition();
+			Vector3DFloat v0NewPos = m_pInputMesh->m_vecVertices[vertexMapper[v0New]].getPosition();
+			Vector3DFloat v1NewPos = m_pInputMesh->m_vecVertices[vertexMapper[v1New]].getPosition();
+			Vector3DFloat v2NewPos = m_pInputMesh->m_vecVertices[vertexMapper[v2New]].getPosition();
 
-					/*Vector3DFloat v0OldPos = m_vecVertices[v0Old].getPosition();
-					Vector3DFloat v1OldPos = m_vecVertices[v1Old].getPosition();
-					Vector3DFloat v2OldPos = m_vecVertices[v2Old].getPosition();
+			/*Vector3DFloat v0OldPos = m_vecVertices[v0Old].getPosition();
+			Vector3DFloat v1OldPos = m_vecVertices[v1Old].getPosition();
+			Vector3DFloat v2OldPos = m_vecVertices[v2Old].getPosition();
 
-					Vector3DFloat v0NewPos = m_vecVertices[v0New].getPosition();
-					Vector3DFloat v1NewPos = m_vecVertices[v1New].getPosition();
-					Vector3DFloat v2NewPos = m_vecVertices[v2New].getPosition();*/
+			Vector3DFloat v0NewPos = m_vecVertices[v0New].getPosition();
+			Vector3DFloat v1NewPos = m_vecVertices[v1New].getPosition();
+			Vector3DFloat v2NewPos = m_vecVertices[v2New].getPosition();*/
 
-					Vector3DFloat OldNormal = (v1OldPos - v0OldPos).cross(v2OldPos - v1OldPos);
-					Vector3DFloat NewNormal = (v1NewPos - v0NewPos).cross(v2NewPos - v1NewPos);
+			Vector3DFloat OldNormal = (v1OldPos - v0OldPos).cross(v2OldPos - v1OldPos);
+			Vector3DFloat NewNormal = (v1NewPos - v0NewPos).cross(v2NewPos - v1NewPos);
 
-					OldNormal.normalise();
-					NewNormal.normalise();
+			OldNormal.normalise();
+			NewNormal.normalise();
 
-					// Note for after holiday - We are still getting faces flipping despite the following test. I tried changing
-					// the 0.0 to 0.9 (which should still let coplanar faces merge) but oddly nothing then merged. Investigate this.
-					float dotProduct = OldNormal.dot(NewNormal);
-					//cout << dotProduct << endl;
-					if(dotProduct < 0.9f)
-					{
-						//cout << "   Face flipped!!" << endl;
+			// Note for after holiday - We are still getting faces flipping despite the following test. I tried changing
+			// the 0.0 to 0.9 (which should still let coplanar faces merge) but oddly nothing then merged. Investigate this.
+			float dotProduct = OldNormal.dot(NewNormal);
+			//cout << dotProduct << endl;
+			if(dotProduct < 0.9f)
+			{
+				//cout << "   Face flipped!!" << endl;
 
-						faceFlipped = true;
+				faceFlipped = true;
 
-						/*vertexLocked[v0] = true;
-						vertexLocked[v1] = true;*/
+				/*vertexLocked[v0] = true;
+				vertexLocked[v1] = true;*/
 
-						break;
-					}
-				}
+				break;
+			}
+		}
 
-				if(faceFlipped == true)
-				{
-					return false;
-				}
+		if(faceFlipped == true)
+		{
+			return false;
+		}
 		return true;
 	}
 
@@ -500,146 +562,226 @@ namespace PolyVox
 	bool MeshDecimator<PositionMaterial>::canCollapseEdge(uint32_t v0, uint32_t v1)
 	{
 		//A vertex will be locked if it has already been involved in a collapse this pass.
-				if(vertexLocked[v0] || vertexLocked[v1])
+		if(vertexLocked[v0] || vertexLocked[v1])
+		{
+			return false;
+		}
+
+		if(m_pInputMesh->m_vecVertices[v0].getMaterial() != m_pInputMesh->m_vecVertices[v1].getMaterial())
+		{
+			return false;
+		}
+
+		if(m_vecInitialVertexMetadata[v0].m_bNormalFlags.count() == 3) //Corner
+		{
+			return false;
+		}
+
+		if(m_vecInitialVertexMetadata[v0].m_bNormalFlags.count() == 2) //Edge
+		{
+			if(isSubsetCubic(m_vecInitialVertexMetadata[v0].m_bNormalFlags, m_vecInitialVertexMetadata[v1].m_bNormalFlags) == false)
+			{
+				return false;
+			}
+		}
+
+		if(m_vecInitialVertexMetadata[v0].m_bNormalFlags.count() == 1) //Face
+		{
+			if(isSubsetCubic(m_vecInitialVertexMetadata[v0].m_bNormalFlags, m_vecInitialVertexMetadata[v1].m_bNormalFlags) == false)
+			{
+				return false;
+			}
+		}
+
+		/*if(m_vecInitialVertexMetadata[v0].normal.dot(m_vecInitialVertexMetadata[v1].normal) < 0.999f)
+		{
+			return false;
+		}*/
+
+		if(m_vecInitialVertexMetadata[v0].hasDuplicate)
+		{
+			bool canMerge = false;
+
+			/*if(m_vecInitialVertexMetadata[v1].hasDuplicate)
+			{
+				if(isSubsetCubic(m_vecInitialVertexMetadata[v0].m_bNormalFlags, m_vecInitialVertexMetadata[v1].m_bNormalFlags))
 				{
-					return false;
+					canMerge = true;
 				}
+			}*/
 
-				if(m_pInputMesh->m_vecVertices[v0].getMaterial() != m_pInputMesh->m_vecVertices[v1].getMaterial())
-				{
-					return false;
-				}
+			return canMerge;
+		}
 
-				if(noOfDifferentNormals[v0] == 3)
-				{
-					return false;
-				}
+		Vector3DFloat offset = static_cast<Vector3DFloat>(m_pInputMesh->m_Region.getLowerCorner());
+		bool v0Inside = m_pInputMesh->m_Region.containsPoint(m_pInputMesh->m_vecVertices[v0].getPosition() + offset);
+		bool v1Inside = m_pInputMesh->m_Region.containsPoint(m_pInputMesh->m_vecVertices[v1].getPosition() + offset);
+		if(v0Inside == false)
+		{
+			return false;
+		}
 
-				/*if(materialKey[v0] != materialKey[v1])
-				{
-					return false;
-				}*/
 
-				/*if((hasDuplicate[v0]) || (hasDuplicate[v1]))
-				{
-					return false;
-				}*/
+		////////////////////////////////////////////////////////////////////////////////
+		//The last test is whether we will flip any of the faces
 
-				if(hasDuplicate[v0])
-				{
-					//if(!hasDuplicate[v1])
-					{
-						return false;
-					}
-				}
+		bool faceFlipped = false;
+		//list<uint32_t> triangles = trianglesUsingVertexCurrently[v0];
+		list<uint32_t> triangles = m_vecCurrentVertexMetadata[v0].trianglesUsingVertex;
+		/*set<uint32_t> triangles;
+		std::set_union(trianglesUsingVertex[v0].begin(), trianglesUsingVertex[v0].end(),
+			trianglesUsingVertex[v1].begin(), trianglesUsingVertex[v1].end(),
+			std::inserter(triangles, triangles.begin()));*/
 
-				if(noOfDifferentNormals[v0] > noOfDifferentNormals[v1])
-				{
-					return false;
-				}
-
-				Vector3DFloat offset = static_cast<Vector3DFloat>(m_pInputMesh->m_Region.getLowerCorner());
-				bool v0Inside = m_pInputMesh->m_Region.containsPoint(m_pInputMesh->m_vecVertices[v0].getPosition() + offset);
-				bool v1Inside = m_pInputMesh->m_Region.containsPoint(m_pInputMesh->m_vecVertices[v1].getPosition() + offset);
-				if(v0Inside == false)
-				{
-					/*if(v1Inside == false)
-					{
-						if(noOfDifferentNormals[v0] > 1)
-						{
-							return false;
-						}
-					}
-					else*/
-					{
-						return false;
-					}
-				}
-
-				////////////////////////////////////////////////////////////////////////////////
-				//The last test is whether we will flip any of the faces
-
-				bool faceFlipped = false;
-				list<uint32_t> triangles = trianglesUsingVertex[v0];
-				/*set<uint32_t> triangles;
-				std::set_union(trianglesUsingVertex[v0].begin(), trianglesUsingVertex[v0].end(),
-					trianglesUsingVertex[v1].begin(), trianglesUsingVertex[v1].end(),
-					std::inserter(triangles, triangles.begin()));*/
-
-				for(list<uint32_t>::iterator triIter = triangles.begin(); triIter != triangles.end(); triIter++)
-				{
-					uint32_t tri = *triIter;
+		for(list<uint32_t>::iterator triIter = triangles.begin(); triIter != triangles.end(); triIter++)
+		{
+			uint32_t tri = *triIter;
 					
-					uint32_t v0Old = m_pInputMesh->m_vecTriangleIndices[tri * 3];
-					uint32_t v1Old = m_pInputMesh->m_vecTriangleIndices[tri * 3 + 1];
-					uint32_t v2Old = m_pInputMesh->m_vecTriangleIndices[tri * 3 + 2];
+			uint32_t v0Old = m_pInputMesh->m_vecTriangleIndices[tri * 3];
+			uint32_t v1Old = m_pInputMesh->m_vecTriangleIndices[tri * 3 + 1];
+			uint32_t v2Old = m_pInputMesh->m_vecTriangleIndices[tri * 3 + 2];
 
-					//Check if degenerate
-					if((v0Old == v1Old) || (v1Old == v2Old) || (v2Old == v0Old))
-					{
-						continue;
-					}
+			//Check if degenerate
+			if((v0Old == v1Old) || (v1Old == v2Old) || (v2Old == v0Old))
+			{
+				continue;
+			}
 
-					uint32_t v0New = v0Old;
-					uint32_t v1New = v1Old;
-					uint32_t v2New = v2Old;
+			uint32_t v0New = v0Old;
+			uint32_t v1New = v1Old;
+			uint32_t v2New = v2Old;
 
-					if(v0New == v0)
-						v0New = v1;
-					if(v1New == v0)
-						v1New = v1;
-					if(v2New == v0)
-						v2New = v1;
+			if(v0New == v0)
+				v0New = v1;
+			if(v1New == v0)
+				v1New = v1;
+			if(v2New == v0)
+				v2New = v1;
 
-					//Check if degenerate
-					if((v0New == v1New) || (v1New == v2New) || (v2New == v0New))
-					{
-						continue;
-					}
+			//Check if degenerate
+			if((v0New == v1New) || (v1New == v2New) || (v2New == v0New))
+			{
+				continue;
+			}
 
-					Vector3DFloat v0OldPos = m_pInputMesh->m_vecVertices[vertexMapper[v0Old]].getPosition();
-					Vector3DFloat v1OldPos = m_pInputMesh->m_vecVertices[vertexMapper[v1Old]].getPosition();
-					Vector3DFloat v2OldPos = m_pInputMesh->m_vecVertices[vertexMapper[v2Old]].getPosition();
+			Vector3DFloat v0OldPos = m_pInputMesh->m_vecVertices[vertexMapper[v0Old]].getPosition();
+			Vector3DFloat v1OldPos = m_pInputMesh->m_vecVertices[vertexMapper[v1Old]].getPosition();
+			Vector3DFloat v2OldPos = m_pInputMesh->m_vecVertices[vertexMapper[v2Old]].getPosition();
 
-					Vector3DFloat v0NewPos = m_pInputMesh->m_vecVertices[vertexMapper[v0New]].getPosition();
-					Vector3DFloat v1NewPos = m_pInputMesh->m_vecVertices[vertexMapper[v1New]].getPosition();
-					Vector3DFloat v2NewPos = m_pInputMesh->m_vecVertices[vertexMapper[v2New]].getPosition();
+			Vector3DFloat v0NewPos = m_pInputMesh->m_vecVertices[vertexMapper[v0New]].getPosition();
+			Vector3DFloat v1NewPos = m_pInputMesh->m_vecVertices[vertexMapper[v1New]].getPosition();
+			Vector3DFloat v2NewPos = m_pInputMesh->m_vecVertices[vertexMapper[v2New]].getPosition();
 
-					/*Vector3DFloat v0OldPos = m_vecVertices[v0Old].getPosition();
-					Vector3DFloat v1OldPos = m_vecVertices[v1Old].getPosition();
-					Vector3DFloat v2OldPos = m_vecVertices[v2Old].getPosition();
+			/*Vector3DFloat v0OldPos = m_vecVertices[v0Old].getPosition();
+			Vector3DFloat v1OldPos = m_vecVertices[v1Old].getPosition();
+			Vector3DFloat v2OldPos = m_vecVertices[v2Old].getPosition();
 
-					Vector3DFloat v0NewPos = m_vecVertices[v0New].getPosition();
-					Vector3DFloat v1NewPos = m_vecVertices[v1New].getPosition();
-					Vector3DFloat v2NewPos = m_vecVertices[v2New].getPosition();*/
+			Vector3DFloat v0NewPos = m_vecVertices[v0New].getPosition();
+			Vector3DFloat v1NewPos = m_vecVertices[v1New].getPosition();
+			Vector3DFloat v2NewPos = m_vecVertices[v2New].getPosition();*/
 
-					Vector3DFloat OldNormal = (v1OldPos - v0OldPos).cross(v2OldPos - v1OldPos);
-					Vector3DFloat NewNormal = (v1NewPos - v0NewPos).cross(v2NewPos - v1NewPos);
+			Vector3DFloat OldNormal = (v1OldPos - v0OldPos).cross(v2OldPos - v1OldPos);
+			Vector3DFloat NewNormal = (v1NewPos - v0NewPos).cross(v2NewPos - v1NewPos);
 
-					OldNormal.normalise();
-					NewNormal.normalise();
+			OldNormal.normalise();
+			NewNormal.normalise();
 
-					// Note for after holiday - We are still getting faces flipping despite the following test. I tried changing
-					// the 0.0 to 0.9 (which should still let coplanar faces merge) but oddly nothing then merged. Investigate this.
-					float dotProduct = OldNormal.dot(NewNormal);
-					//cout << dotProduct << endl;
-					if(dotProduct < 0.9f)
-					{
-						//cout << "   Face flipped!!" << endl;
+			// Note for after holiday - We are still getting faces flipping despite the following test. I tried changing
+			// the 0.0 to 0.9 (which should still let coplanar faces merge) but oddly nothing then merged. Investigate this.
+			float dotProduct = OldNormal.dot(NewNormal);
+			//cout << dotProduct << endl;
+			if(dotProduct < 0.9f)
+			{
+				//cout << "   Face flipped!!" << endl;
 
-						faceFlipped = true;
+				faceFlipped = true;
 
-						/*vertexLocked[v0] = true;
-						vertexLocked[v1] = true;*/
+				/*vertexLocked[v0] = true;
+				vertexLocked[v1] = true;*/
 
-						break;
-					}
-				}
+				break;
+			}
+		}
 
-				if(faceFlipped == true)
-				{
-					return false;
-				}
+		if(faceFlipped == true)
+		{
+			return false;
+		}
+		
+
+		/*if(m_vecInitialVertexMetadata[v0].trianglesUsingVertex.size() != m_vecInitialVertexMetadata[v1].trianglesUsingVertex.size())
+		{
+			return false;
+		}*/
+
+		/*if(m_vecInitialVertexMetadata[v0].noOfDifferentNormals == 3)
+		{
+			return false;
+		}*/
+
+		/*if(m_vecInitialVertexMetadata[v1].noOfDifferentNormals == 3)
+		{
+			return false;
+		}*/
+
+		/*if(m_vecInitialVertexMetadata[v0].noOfDifferentNormals > m_vecInitialVertexMetadata[v1].noOfDifferentNormals)
+		{
+			return false;
+		}*/
+
+		/*if(m_vecInitialVertexMetadata[v0].materialKey != m_vecInitialVertexMetadata[v1].materialKey)
+		{
+			return false;
+		}*/
+
+		/*if((hasDuplicate[v0]) || (hasDuplicate[v1]))
+		{
+			return false;
+		}*/
+
+		/*if(m_vecInitialVertexMetadata[v0].hasDuplicate)
+		{
+			//if(!hasDuplicate[v1])
+			{
+				return false;
+			}
+		}*/
 		return true;
 	}
+
+	template <typename VertexType>
+	int MeshDecimator<VertexType>::countZeros(void)
+	{
+		int total = 0;
+		for(int ct = 0; ct < m_pInputMesh->m_vecTriangleIndices.size(); ct++)
+		{
+			if(m_pInputMesh->m_vecTriangleIndices[ct] == 0)
+			{
+				total++;
+			}
+		}
+		return total;
+	}
 }
+
+#ifdef BLAH
+		Vector3DFloat offset = static_cast<Vector3DFloat>(m_pInputMesh->m_Region.getLowerCorner());
+		bool v0Inside = m_pInputMesh->m_Region.containsPoint(m_pInputMesh->m_vecVertices[v0].getPosition() + offset);
+		bool v1Inside = m_pInputMesh->m_Region.containsPoint(m_pInputMesh->m_vecVertices[v1].getPosition() + offset);
+		if(v0Inside == false)
+		{
+			/*if(v1Inside == false)
+			{
+				//if(noOfDifferentNormals[v0] > 1)
+				if(trianglesUsingVertex[v0].size() != trianglesUsingVertex[v1].size())
+				//if(materialKey[v0] != materialKey[v1])
+				{
+					return false;
+				}
+			}
+			else*/
+			{
+				return false;
+			}
+		}
+#endif
\ No newline at end of file