路径追踪的理论与实现：代码实现

上面两篇博客主要讲述了路径追踪渲染的理论，这次来展示一下它的代码实现。这部分代码主要是我参考PBRT写的，没有实现光的折射，但也足以帮助理解路径追踪算法的原理。
首先，是遍历图片中的所有像素，对每个像素进行采样,注意需要将得到的RGB值进行gamma矫正：

for (int i = 0; i < width; i++)
{
    for (int j = 0; j < height; j++)
    {	
	Pixel p = sample_pixel(camera, i, j, width, height);
	XMFLOAT3 rgb = Spectrum::XYZToRGB(p.xyz);
	if (p.filterWeightSum > 0.0)
	{
		rgb.x /= p.filterWeightSum;
		rgb.y /= p.filterWeightSum;
		rgb.z /= p.filterWeightSum;
	}
	int r = int(MathHelper::Clamp<float>(GammaCorrect(rgb.x) *255.0f + 0.5f, 0.0, 255.0));
        int g = int(MathHelper::Clamp<float>(GammaCorrect(rgb.y) *255.0f + 0.5f, 0.0, 255.0));
        int b = int(MathHelper::Clamp<float>(GammaCorrect(rgb.z) *255.0f + 0.5f, 0.0, 255.0));
	image.setPixelColor(QPoint(i, height - 1 - j), QColor(r, g, b));
	}
}

具体采样的代码如下，就是在像素内部随机选一个点，打出一条光线计算光强。注意还需要计算一下采样点的权重，在这里我选用了三角滤波进行加权：

Pixel PathTracingRenderer::sample_pixel(Camera* camera, int x, int y, int width, int height)
{
	float sample_x = 0.0f;
	float sample_y = 0.0f;
	Spectrum r;
	Pixel p;
	for (int i = 0; i < sample_count; i++)
	{
		sample_x = x + generateRandomFloat();
		sample_y = y + generateRandomFloat();
		float w = TriangleFilterEval(sample_x - x - 0.5f, sample_y - y - 0.5f, 1.0f);
		Ray ray = camera->getRay(sample_x / width, sample_y / height);
		Spectrum radiance = Li(ray);
		p.xyz = MathHelper::AddFloat3(p.xyz, RGBToXYZ((radiance * w).getFloat3()));
		p.filterWeightSum += w;
	}
	return p;
}

接下来我们深入函数Li，看看究竟是怎么计算光强的。由于代码有点多，所以我将要点都予以注释：

Spectrum PathTracingRenderer::Li(const Ray& r)
{
	bool specularBounce = false;
	Spectrum L, beta(1.0f, 1.0f, 1.0f);
	Ray ray(r);
        // 逐步增加路径条数，计算每次路径的P值
	for (int bounce = 0; bounce < max_bounce; bounce++)
	{
	    IntersectInfo it;
            //对场景中的物体做碰撞检测，可以利用BVH和KD-Tree加快速度
	    g_pGlobalSys->cast_ray_to_get_intersection(ray, it);
	    if (bounce == 0 || specularBounce)
	    {
		if (!it.isSurfaceInteraction())
		{
		    //如果没有碰撞到物体，计算所有光的环境光强
                    auto& lights = g_pGlobalSys->objectManager.getAllLights();
		    for each (Light* light in lights)
		        L += (beta* light->Le(ray));
		}
		else
		{
		    //如果碰撞到的是AreaLight，则计算AreaLight的光强（Le函数内部会剔除掉非AreaLight）
                    L += beta * it.Le(MathHelper::NegativeFloat3(ray.direction));
		}
	}

	if (!it.isSurfaceInteraction() || bounce >= max_bounce)
	    break;
        //根据碰撞到物体的材质赋予其BSDF
	it.ComputeScatteringFunctions();
        //计算直接光照的强度，相当于物体本身的发光
	L += beta * UniformSampleOneLight(it);
	XMFLOAT3 wo(-ray.direction.x, -ray.direction.y, -ray.direction.z), wi;
	float pdf;
	BxDFType flags;
	float f1 = generateRandomFloat(), f2 = generateRandomFloat();
	//根据BRDF采样出射方向wi和pdf
        Spectrum f = it.bsdf->Sample_f(wo, &wi, XMFLOAT2(f1, f2), &pdf, BSDF_ALL, &flags);
	if (f.isBlack() || pdf == 0.0f)
	    break;
	beta *= (f * abs(MathHelper::DotFloat3(wi, it.normal))/pdf);
	specularBounce = (flags&BxDFType::BSDF_SPECULAR) != 0;
	ray = it.spawnRay(wi);

	// 俄罗斯轮盘计算权重和决定是否要将路径终止
	if (bounce > 3)
	{
	    float q = std::max<float>(0.05f, 1.0 - beta.y());
	    if (generateRandomFloat() < q)
	        break;
	    beta /= (1.0 - q);
	}
    }
    return L;
}

值得一提的是，UniformSampleOneLight函数会随机挑选一个光源进行直接光照的计算，具体的计算过程用的是复合重要性采样的思想：

Spectrum PathTracingRenderer::UniformSampleOneLight(const IntersectInfo& it)
{
	std::vector<Light*> lights = g_pGlobalSys->objectManager.getAllLights();
	Light* sel_light = lights[rand() % (lights.size())];
	XMFLOAT2 uScattering(generateRandomFloat(), generateRandomFloat());
	XMFLOAT2 uLight(generateRandomFloat(), generateRandomFloat());
	float light_count = g_pGlobalSys->objectManager.getLightsCountParameter();
	return light_count * EstimateDirect(it, uScattering, sel_light, uLight);
}

Spectrum EstimateDirect(const IntersectInfo& it, XMFLOAT2 uScattering, Light* light, XMFLOAT2 uLight, bool specular)
{
	Spectrum ld;
	BxDFType bsdfFlags = specular ? BSDF_ALL : BxDFType(BSDF_ALL & ~BSDF_SPECULAR);
	// sample light source with multiple importance sampling
	XMFLOAT3 wi;
	float light_pdf = 0.0f, scattering_pdf = 0.0f;
	VisibilityTester vt;
	Spectrum li = light->sample_li(it, uLight, &wi, &light_pdf, vt);
	if (light_pdf > 0.0f && !li.isBlack())
	{
		Spectrum f;
		if (it.isSurfaceInteraction())
		{
			f = it.bsdf->f(it.wo, wi, bsdfFlags) * MathHelper::DotFloat3(wi, it.normal);
			scattering_pdf = it.bsdf->Pdf(it.wo, wi, bsdfFlags);
		}
		if (!f.isBlack())
		{
			if (!vt.unoccluded())
				li = Spectrum();
			if (!li.isBlack())
			{
				if (light->isDelta())
					ld += f * li / light_pdf;
				else
				{
					float weight = PowerHeuristic(1, light_pdf, 1, scattering_pdf);
					ld += f * li * weight / light_pdf;
				}
			}
		}
	}

	// sample BSDF with multiple importance sampling
	if (!light->isDelta())
	{
		Spectrum f;
		bool sampledSpecular = false;
		BxDFType sampled_type;
		if (it.isSurfaceInteraction())
		{ 
			f = it.bsdf->Sample_f(it.wo, &wi, uScattering, &scattering_pdf, bsdfFlags, &sampled_type);
			f = f * MathHelper::DotFloat3(wi, it.normal);
			sampledSpecular = sampled_type & BSDF_SPECULAR;
		}

		if (!f.isBlack() && scattering_pdf > 0.0f)
		{
			float weight = 1.0;
			if (!sampledSpecular)
			{
				light_pdf = light->Pdf_Li(it, wi);
				if (light_pdf == 0.0)
					return ld;
				weight = PowerHeuristic(1, scattering_pdf, 1, light_pdf);
			}
			IntersectInfo light_it;
			Ray ray = it.spawnRay(wi);
			g_pGlobalSys->cast_ray_to_get_intersection(ray, light_it);
			Spectrum li;
			if (light_it.isSurfaceInteraction())
			{
				if (light_it.obj->getType() == AREA_LIGHT)
					li = light_it.Le(MathHelper::NegativeFloat3(wi));
			}
			else
				li = light->Le(ray);
			if (!li.isBlack())
				ld = ld + f * li * weight / scattering_pdf;
		}
	}
	return ld;
}

主要的代码结构就是这样，接下来展示一下我利用上述代码实现的渲染结果：


经和PBRT比较，同样的场景和参数，渲染结果差距不大。