void raytrace(RGBA* pixelData, const Ray& rays, const int iteration, int w, int h)
{
	if(iteration > 3)
		return;

	// Initialisation state of no intersection, -1.
	int nearest[4] = { -1, -1, -1, -1 };

	int uniqueSpheres = 0;
	int spheresHit[4];

	// Set the nearest intersection to as large as possible.
	SSEInt nearestDistance = _mm_set_ps1(0xffffffff);

	// Used to store which ray hits which sphere in the vector.
	int sphereIndex = 0;

	// For every sphere in the scene see if the ray intersects it.
	for(vector<RTSphere>::const_iterator s = spheres.begin(); s != spheres.end(); s++)
	{
		// Intersect the packet of rays with the sphere and have the distance to 
		// the intersection point returned.
		SSEInt distance = s->IntersectTest(rays);

		for(int r = 0; r < 4; r++)
		{
			// If the intersection distance is greater than zero and it is nearer
			// than the previous nearest intersection, this is our nearest sphere.
			if(distance.m128_f32[r] > 0 && 
				distance.m128_f32[r] < nearestDistance.m128_f32[r])
			{
				nearest[r] = sphereIndex;	// Store index of the nearest sphere.
				nearestDistance.m128_f32[r] = distance.m128_f32[r];
			}
		}

		sphereIndex++;
	}

	for(int sh = 0; sh < uniqueSpheres; sh++)
	{
		int sphereIndex = spheresHit[sh];
		RTSphere& sphere = spheres[sphereIndex];	// The sphere to be tested.
		
		// Calculate the distance to intersection point minus the small amount
		// to avoid the intersection point actually intersecting the sphere.
		SSEInt mult = _mm_mul_ps(nearestDistance, _mm_set_ps1(1.0f - EPSILON));

		// The normalised ray vector multiplied by the distance to the
		// intersection point, less the small "nearly 1" multiplier - 'mult'.
		SSEInt extraX = _mm_mul_ps(rays.directionX, mult);
		SSEInt extraY = _mm_mul_ps(rays.directionY, mult);
		SSEInt extraZ = _mm_mul_ps(rays.directionZ, mult);

		// The intersection point is the ray's initial position plus the ray's
		// direction, after being multiplied by mult.
		SSEInt intPointX = _mm_add_ps(extraX, rays.positionX);
		SSEInt intPointY = _mm_add_ps(extraY, rays.positionY);
		SSEInt intPointZ = _mm_add_ps(extraZ, rays.positionZ);

		// For every light in the scene, we need to check if it is casting
		// light onto an object for both diffuse and specular lighting.
		for (std::vector<RTLight>::size_type i = 0; i != lights.size(); ++i)
		{
			RTLight light = lights[i];

			// The i'th light's position.
			SSEInt lightPosX = _mm_set_ps1(light.position.x);
			SSEInt lightPosY = _mm_set_ps1(light.position.y);
			SSEInt lightPosZ = _mm_set_ps1(light.position.z);

			// The i'th light's direction.
			SSEInt lightDirX = _mm_sub_ps(lightPosX, intPointX);
			SSEInt lightDirY = _mm_sub_ps(lightPosY, intPointY);
			SSEInt lightDirZ = _mm_sub_ps(lightPosZ, intPointZ);

			// Normalise the light direction.
			NormalizeSSE(lightDirX, lightDirY, lightDirZ);

			V3& spherePosition = sphere.GetPosition();

			// Calculate the normal at the intersection point, for a sphere this is
			// simply the intersection point - sphere centre, normalised.
			SSEInt sphereNormalX = _mm_sub_ps(intPointX, _mm_set_ps1(spherePosition.x));
			SSEInt sphereNormalY = _mm_sub_ps(intPointY, _mm_set_ps1(spherePosition.y));
			SSEInt sphereNormalZ = _mm_sub_ps(intPointZ, _mm_set_ps1(spherePosition.z));

			NormalizeSSE(sphereNormalX, sphereNormalY, sphereNormalZ);

			float reflectionFactor = sphere.GetReflection();

			if(reflectionFactor > 0)
			{
				SSEInt reflectionX, reflectionY, reflectionZ;

				ReflectSSE(rays.directionX, rays.directionY, rays.directionZ, 
							sphereNormalX, sphereNormalY, sphereNormalZ,
							reflectionX, reflectionY, reflectionZ);

				NormalizeSSE(reflectionX, reflectionY, reflectionZ);

				Ray reflectedPacket;
				
				reflectedPacket.directionX = reflectionX;
				reflectedPacket.directionY = reflectionY;
				reflectedPacket.directionZ = reflectionZ;
				reflectedPacket.positionX = intPointX;
				reflectedPacket.positionY = intPointY;
				reflectedPacket.positionZ = intPointZ;

				raytrace(pixelData, reflectedPacket, iteration + 1, w, h);
				
				for(int r = 0; r < 4; r++)
				{
					pixelData[r].blue *= reflectionFactor;
					pixelData[r].green *= reflectionFactor;
					pixelData[r].red *= reflectionFactor;
				}
			}

			// Calculate dot product for Diffuse Lighting.
			const SSEInt dotProduct = DotSSE(sphereNormalX, sphereNormalY, sphereNormalZ, 
										lightDirX, lightDirY, lightDirZ);

			// Since not every ray in the packet of four strikes this particular sphere
			// they need to be dealt with seperately right at the end.
			for(int r = 0; r < 4; r++)
			{
				// If this ray hit this sphere, and it's dot product indicates it is lit...
				if(nearest[r] == sphereIndex && dotProduct.m128_f32[r] > 0)
				{
					RGBA& sphereColour = sphere.GetColour();
					float sphereDiffuse = sphere.GetDiffuse();

					// Set the red, green and blue component in an int so as to avoid
					// overflowing a byte in case of "brighter than white" pixels.
					int red = pixelData[r].red + (sphereColour.red * dotProduct.m128_f32[r] * light.power * sphereDiffuse);
					int green = pixelData[r].green + (sphereColour.green * dotProduct.m128_f32[r] * light.power * sphereDiffuse);
					int blue = pixelData[r].blue + (sphereColour.blue * dotProduct.m128_f32[r] * light.power * sphereDiffuse);

					// Clamp the RGB components back within the 0 - 255 range.
					if(red > 255)
						red = 255;
					if(green > 255)
						green = 255;
					if(blue > 255)
						blue = 255;

					// Finally set the pixel data.
					pixelData[r].red = red; pixelData[r].green = green; pixelData[r].blue = blue;
				}
			}

			// We only need calculate the specular lighting component for this
			// sphere if it has a specular value greater than zero.
			if(sphere.GetSpecular() > 0)
			{
				// No shadowing yet...
				float specularShade = 1;

				// Calculate the vector from intersection point back to the light.
				SSEInt toLightVectorX = _mm_sub_ps(lightPosX, intPointX);
				SSEInt toLightVectorY = _mm_sub_ps(lightPosY, intPointY);
				SSEInt toLightVectorZ = _mm_sub_ps(lightPosZ, intPointZ);

				NormalizeSSE(toLightVectorX, toLightVectorY, toLightVectorZ);

				// Calculate the vector along which light will be reflected off the
				// surface of the sphere and store it in reflectionX/Y/Z.
				SSEInt reflectionX, reflectionY, reflectionZ;
				sphere.ReflectRayAtPoint(toLightVectorX, toLightVectorY, toLightVectorZ, 
											intPointX, intPointY, intPointZ, 
											reflectionX, reflectionY, reflectionZ);

				// Calculate the specular dot product 
				const SSEInt specDP = DotSSE(rays.directionX, rays.directionY, rays.directionZ, 
												reflectionX, reflectionY, reflectionZ);

				for(int r = 0; r < 4; r++)
				{
					// If this ray is hitting this sphere, and it's dot product > 0
					if(nearest[r] == sphereIndex && specDP.m128_f32[r] > 0)
					{
						const float specular = pow(specDP.m128_f32[r], 10) * sphere.GetSpecular();

						// Calculate the new RGB components, being sure not to overflow.
						int red = pixelData[r].red + (specular * 255);
						int green = pixelData[r].green + (specular * 255);
						int blue = pixelData[r].blue + (specular * 255);

						// Clamp back within 0 - 255 range.
						if(red > 255)
							red = 255;
						if(green > 255)
							green = 255;
						if(blue > 255)
							blue = 255;

						// Set the pixel data with the new values.
						pixelData[r].red = red; pixelData[r].green = green; pixelData[r].blue = blue;
					}
				}
			}
		}
	}
}