diff --git a/library/PolyVoxCore/include/PolyVoxCore/AmbientOcclusionCalculator.h b/library/PolyVoxCore/include/PolyVoxCore/AmbientOcclusionCalculator.h
index 0cfdef5b..f5b75bf4 100644
--- a/library/PolyVoxCore/include/PolyVoxCore/AmbientOcclusionCalculator.h
+++ b/library/PolyVoxCore/include/PolyVoxCore/AmbientOcclusionCalculator.h
@@ -67,6 +67,131 @@ namespace PolyVox
 
 		polyvox_function<bool(const typename VolumeType::VoxelType& voxel)> m_funcIsTransparent;
 	};
+
+	class AmbientOcclusionCalculatorRaycastCallback
+	{
+	public:
+		AmbientOcclusionCalculatorRaycastCallback()
+		{
+		}
+
+		bool operator()(const SimpleVolume<uint8_t>::Sampler& sampler)
+		{
+			return sampler.getVoxel() == 0;
+		}
+	};
+
+	template<typename VolumeType>
+	void calculateAmbientOcclusion(VolumeType* volInput, Array<3, uint8_t>* arrayResult, Region region, float fRayLength, uint8_t uNoOfSamplesPerOutputElement, polyvox_function<bool(const typename VolumeType::VoxelType& voxel)> funcIsTransparent)
+	{
+		Region m_region = region;
+		typename VolumeType::Sampler m_sampVolume(volInput);
+		VolumeType* m_volInput = volInput;
+		Array<3, uint8_t>* m_arrayResult = arrayResult;
+		float m_fRayLength = fRayLength;
+
+		uint8_t m_uNoOfSamplesPerOutputElement = uNoOfSamplesPerOutputElement;
+
+		uint16_t mRandomUnitVectorIndex = 0;
+		uint16_t mRandomVectorIndex = 0;
+		uint16_t mIndexIncreament;
+
+		polyvox_function<bool(const typename VolumeType::VoxelType& voxel)> m_funcIsTransparent = funcIsTransparent;
+
+		//Make sure that the size of the volume is an exact multiple of the size of the array.
+		assert(m_volInput->getWidth() % arrayResult->getDimension(0) == 0);
+		assert(m_volInput->getHeight() % arrayResult->getDimension(1) == 0);
+		assert(m_volInput->getDepth() % arrayResult->getDimension(2) == 0);
+
+		//Our initial indices. It doesn't matter exactly what we set here, but the code below makes 
+		//sure they are different for different regions which helps reduce tiling patterns in the results.
+		mRandomUnitVectorIndex += m_region.getLowerCorner().getX() + m_region.getLowerCorner().getY() + m_region.getLowerCorner().getZ();
+		mRandomVectorIndex += m_region.getLowerCorner().getX() + m_region.getLowerCorner().getY() + m_region.getLowerCorner().getZ();
+
+		//This value helps us jump around in the array a bit more, so the
+		//nth 'random' value isn't always followed by the n+1th 'random' value.
+		mIndexIncreament = 1;
+
+		const int iRatioX = m_volInput->getWidth()  / m_arrayResult->getDimension(0);
+		const int iRatioY = m_volInput->getHeight() / m_arrayResult->getDimension(1);
+		const int iRatioZ = m_volInput->getDepth()  / m_arrayResult->getDimension(2);
+
+		const float fRatioX = iRatioX;
+		const float fRatioY = iRatioY;
+		const float fRatioZ = iRatioZ;
+		const Vector3DFloat v3dRatio(fRatioX, fRatioY, fRatioZ);
+
+		const float fHalfRatioX = fRatioX * 0.5f;
+		const float fHalfRatioY = fRatioY * 0.5f;
+		const float fHalfRatioZ = fRatioZ * 0.5f;
+		const Vector3DFloat v3dHalfRatio(fHalfRatioX, fHalfRatioY, fHalfRatioZ);
+
+		const Vector3DFloat v3dOffset(0.5f,0.5f,0.5f);
+
+		//RaycastResult raycastResult;
+		//Raycast<VolumeType> raycast(m_volInput, Vector3DFloat(0.0f,0.0f,0.0f), Vector3DFloat(1.0f,1.0f,1.0f), raycastResult, polyvox_bind(&PolyVox::AmbientOcclusionCalculator<VolumeType>::raycastCallback, this, std::placeholders::_1));
+
+		//This loop iterates over the bottom-lower-left voxel in each of the cells in the output array
+		for(uint16_t z = m_region.getLowerCorner().getZ(); z <= m_region.getUpperCorner().getZ(); z += iRatioZ)
+		{
+			for(uint16_t y = m_region.getLowerCorner().getY(); y <= m_region.getUpperCorner().getY(); y += iRatioY)
+			{
+				for(uint16_t x = m_region.getLowerCorner().getX(); x <= m_region.getUpperCorner().getX(); x += iRatioX)
+				{
+					//Compute a start position corresponding to 
+					//the centre of the cell in the output array.
+					Vector3DFloat v3dStart(x, y, z);
+					v3dStart -= v3dOffset;
+					v3dStart += v3dHalfRatio;
+
+					//Keep track of how many rays did not hit anything
+					uint8_t uVisibleDirections = 0;
+
+					for(int ct = 0; ct < m_uNoOfSamplesPerOutputElement; ct++)
+					{						
+						//We take a random vector with components going from -1 to 1 and scale it to go from -halfRatio to +halfRatio.
+						//This jitter value moves our sample point from the center of the array cell to somewhere else in the array cell
+						Vector3DFloat v3dJitter = randomVectors[(mRandomVectorIndex += (++mIndexIncreament)) % 1019]; //Prime number helps avoid repetition on sucessive loops.
+						v3dJitter *= v3dHalfRatio;
+						const Vector3DFloat v3dRayStart = v3dStart + v3dJitter;
+
+						Vector3DFloat v3dRayDirection = randomUnitVectors[(mRandomUnitVectorIndex += (++mIndexIncreament)) % 1021]; //Differenct prime number.
+						v3dRayDirection *= m_fRayLength;
+						
+						/*raycast.setStart(v3dRayStart);
+						raycast.setDirection(v3dRayDirection);
+						raycast.execute();
+
+						if(raycastResult.foundIntersection == false)
+						{
+							++uVisibleDirections;
+						}*/
+
+						MyRaycastResult result = raycastWithDirection(m_volInput, v3dRayStart, v3dRayDirection, AmbientOcclusionCalculatorRaycastCallback());
+						if(result == MyRaycastResults::Completed)
+						{
+							++uVisibleDirections;
+						}
+					}
+
+					float fVisibility;
+					if(m_uNoOfSamplesPerOutputElement == 0)
+					{
+						//The user might request zero samples (I've done this in the past while debugging - I don't want to
+						//wait for ambient occlusion but I do want as valid result for rendering). Avoid the divide by zero.
+						fVisibility = 1.0f;
+					}
+					else
+					{
+						fVisibility = static_cast<float>(uVisibleDirections) / static_cast<float>(m_uNoOfSamplesPerOutputElement);
+						assert((fVisibility >= 0.0f) && (fVisibility <= 1.0f));
+					}
+
+					(*m_arrayResult)[z / iRatioZ][y / iRatioY][x / iRatioX] = static_cast<uint8_t>(255.0f * fVisibility);
+				}
+			}
+		}
+	}
 }
 
 #include "PolyVoxCore/AmbientOcclusionCalculator.inl"
diff --git a/tests/TestAmbientOcclusionGenerator.cpp b/tests/TestAmbientOcclusionGenerator.cpp
index 7491f7c5..18402c96 100644
--- a/tests/TestAmbientOcclusionGenerator.cpp
+++ b/tests/TestAmbientOcclusionGenerator.cpp
@@ -61,11 +61,8 @@ void TestAmbientOcclusionGenerator::testExecute()
 	const int32_t g_uArraySideLength = g_uVolumeSideLength / 2;
 	Array<3, uint8_t> ambientOcclusionResult(ArraySizes(g_uArraySideLength)(g_uArraySideLength)(g_uArraySideLength));
 
-	//Create the ambient occlusion calculator
-	AmbientOcclusionCalculator< SimpleVolume<uint8_t> > calculator(&volData, &ambientOcclusionResult, volData.getEnclosingRegion(), 32.0f, 255, isVoxelTransparent);
-
-	//Execute the calculator
-	calculator.execute();
+	// Calculate the ambient occlusion values
+	calculateAmbientOcclusion(&volData, &ambientOcclusionResult, volData.getEnclosingRegion(), 32.0f, 255, isVoxelTransparent);
 	
 	//Check the results by sampling along a line though the centre of the volume. Because
 	//of the two walls we added, samples in the middle are darker than those at the edge.