Debugging spectral tracking

spectral_tracking.cpp

float3 direct_light(const Ray& ray, const float3& p) {
	// return lighting at this point... 
}

float3 spectral_tracking(const Ray& ray, const Material& material, float3& transmittance) {
	const float d = ray.max_t - ray.min_t;

	float3 w(1.f);
	float t = 0.f;

	const float mt = math::max_component(material.max_extinction());
	while (true) {
		const float r = rng_.random_float();
		t = t -std::log(1.f - r) / mt;
		if (t > d) {
			transmittance = w;
			return float3(0.f);
		}		
		
		const float3 p = ray.point(ray.min_t + t);

		const float3 sigma_a = material.absorption(p);

		const float3 sigma_s = material.scattering(p);

		const float3 sigma_t = sigma_a + sigma_s;

		const float3 sigma_n = float3(mt) - sigma_t;

		const float msa = math::max_component(sigma_a);
		const float mss = math::max_component(sigma_s);
		const float msn = math::max_component(sigma_n);
		const float c = 1.f / (msa + mss + msn);

		const float pa = msa * c;
		const float ps = mss * c;
		const float pn = msn * c;

		const float3 wa = (sigma_a / (mt * pa));
		const float3 ws = (sigma_s / (mt * ps));
		const float3 wn = (sigma_n / (mt * pn));

		const float r2 = rng_.random_float();
		if (r2 < pa) {
			transmittance = float3(0.f);
			return float3(0.f);
		} else if (r2 < 1.f - pn) {
			transmittance = float3(0.f);
			return w * ws * direct_light(ray, p, worker);
		} else {
			w *= wn;
		}
	}
}