float3 direct_light(const float3& position) {
	// direct light inside volume including phase function and beam transmittance
}

float3 homogeneous2(const Ray& ray, float3& transmittance) {
	float3 radiance(0.f);

	if (Algorithm::Tracking == algorithm) {
		const float3 sigma_a = material.absorption(float3::identity());

		const float3 sigma_s = material.scattering(float3::identity());

		const float3 extinction = sigma_a + sigma_s;

		const float3 scattering_albedo = sigma_s / extinction;

		transmittance = math::exp(-d * extinction);

		const float r = rng_.random_float();
		const float scatter_distance = -std::log(1.f - r * (1.f - average(transmittance))) / average(extinction);

		const float3 p = ray.point(scatter_distance);

		float3 l = direct_light(p);

		l *= (1.f - transmittance) * scattering_albedo;

		radiance = l;
	} else if (Algorithm::Delta_tracking == algorithm) {
		const float max_extinction = average(material.max_extinction());
		bool terminated = false;
		float t = 0.f;

		float3 p;
		float3 extinction;
		float3 scattering_albedo;

		do {
			const float r = rng_.random_float();
			t = t -std::log(1.f - r) / max_extinction;
			if (t > d) {
				break;
			}

			p = ray.point(t);

			const float3 sigma_a = material.absorption(p);

			const float3 sigma_s = material.scattering(p);

			extinction = sigma_a + sigma_s;
			const float r2 = rng_.random_float();
			if (r2 < average(extinction) / max_extinction) {
				terminated = true;

				scattering_albedo = sigma_s / extinction;
			}
		} while (!terminated);

		if (terminated) {
			float3 l = direct_light(p);

			l *= scattering_albedo * extinction;

			radiance = l;

			transmittance = float3(0.f);
		} else {
			radiance = float3(0.f);
			transmittance = float3(1.f);
		}
	} else if (Algorithm::Experiment == algorithm) {
		// Exactly the same as Delta-tracking,
		// except that L is multiplied by scattering_albedo, instead of sigma_s
		const float max_extinction = average(material.max_extinction());
		bool terminated = false;
		float t = 0.f;

		float3 p;
		float3 scattering_albedo;

		do {
			const float r = rng_.random_float();
			t = t -std::log(1.f - r) / max_extinction;
			if (t > d) {
				break;
			}

			p = ray.point(t);

			const float3 sigma_a = material.absorption(p);

			const float3 sigma_s = material.scattering(p);

			const float3 extinction = sigma_a + sigma_s;
			const float r2 = rng_.random_float();
			if (r2 < average(extinction) / max_extinction) {
				terminated = true;

				scattering_albedo = sigma_s / extinction;
			}
		} while (!terminated);

		if (terminated) {
			float3 l = direct_light(p);

			l *= scattering_albedo;

			radiance = l;

			transmittance = float3(0.f);
		} else {
			radiance = float3(0.f);
			transmittance = float3(1.f);
		}
	}

	return radiance;
}