Comparing homogeneous single scattering algorithms from pixar course notes

tracking.cpp

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

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

	const float3 extinction = sigma_a + sigma_s;

	const float3 scattering_albedo = sigma_s / extinction;

	// length of the volume interval
	const float d = ray.max_t - ray.min_t;

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

	if (Algorithm::PBRT == algorithm) {
		// PBRT style raymarching
	} else if (Algorithm::Algorithm_1 == algorithm) {
		// Algorithm 1 on page 12
		const float r = rng_.random_float();
		const float scatter_distance = -std::log(1.f - r) / average(extinction);

		if (scatter_distance > d) {
			// boundary
			return float3(0.f);
		}

		const float r2 = rng_.random_float();
		if (r2 < spectrum::average(sigma_a) / average(extinction)) {
			// absorption/emission
			return float3(0.f);
		}

		const float3 p = ray.point(ray.min_t + scatter_distance);

		const float3 l = direct_light(p);

		radiance += l;
	} else if (Algorithm::Listing_10 == algorithm) {
		// Code listing 10 on page 35-36
		const float r = rng_.random_float();
		const float scatter_distance = -std::log(1.f - r * (1.f - average(transmittance))) / average(extinction);

		const float3 tr = math::exp(-scatter_distance * extinction);

		const float3 p = ray.point(ray.min_t + scatter_distance);

		float3 l = direct_light(p);

		l *= extinction * scattering_albedo * tr;

		radiance += l;
	} else if (Algorithm::Experiment == algorithm) {
		// basically PBRT raymarching with only one step
		const float r = rng_.random_float();
		const float scatter_distance = r * d;

		const float3 tr = math::exp(-scatter_distance * extinction);

		const float3 p = ray.point(ray.min_t + scatter_distance);

		float3 l = direct_light(p);

		l *= d * extinction * scattering_albedo * tr;

		radiance += l;
	}

	return radiance;
}