【光线追踪系列十七】直接光源采样

本文主要参照 Ray Tracing: The Rest of Your Life，其中只是主要精炼光追相关理论，具体实现可参照原文。

前面文章中我们也提到要减少噪点，就是图片上的黑点点。简单说一下为什么会有那么多小点点，就是因为光线路径中没有触碰到光源，路径计算之后就会是黑色的点，我们可以通过发射大量的光线，比如计算每个像素点的时候发射8k~1w条采样光线进行路径计算；也可以路径计算方面做文章：比如加深路径计算递归深度；等等诸如此类。但是上述方法都是暴力解决法，相当耗时，我们可以运用数学对其进行优化，从而实现画质和效率的双面提升，即直接对光源进行采样！

一、理论

我们朝光源方向发送光线或者生成朝向光源的随机方向都是很容易实现的，但是我们需要知道的是，pdf（direction）是什么，如下图对光源采样。

对于一个光源区域A，如果我们均匀采样该区域，那么这个pdf就等于1/A(意思就是每个点的概率均等)。

因此，我们希望尽可能地往重要的方向采样，但在各方向均匀采样，并不会使得重要方向的采样数量变多。接下来，我们人为地让更多的射线往光照的方向采样。

我们可以向着光源的位置生成一个随机的方向，只需在光源材质上随机选取一个点，然后向着这个点的位置生成射线。但我们还需要知道这个pdf(direction)。假设光源的面积为A，如果对这个光源表面随机采样，则该光源表面的pdf为1/A，但如果要从一个定义了方向的单位球的球心为起点，发射射线到这个光源表面，则需要求解对应的pdf公式。

如上图所示，dA 表示光源表面面积的一小部分，设采样到这一小部分面积的概率为 p_q(q) * dA（采样比例乘以微分区域）,也就是 dA/A。

对于单位球而言，dw 为对应 dA 的单位球表面面积的一小部分，设采样到这一小部分的概率 p(direction) ⋅dw 。

则 dw 和 dA 有如下对应关系：

即：方位角微分区域：光源微分区域分成（球心到A中心距离平方）份，取其中的 cosα 代表的份额数。

因为采样到 dw 和 dA 的概率必须相等，则有：

故，

二、代码实现

上述理论整合为代码，如下：

float lightStartX = 213;
float lightEndX = 343;
float lightStartZ = 227;
float lightEndZ = 332;
float lightY = 554;

hittable *random_scene() {
    ...
    list[i++] = new xz_rect(lightStartX, lightEndX, lightStartZ, lightEndZ, lightY, light);
    ...
}

vec3 color(const ray& r, hittable *world, int depth) {
    hit_record rec;
    if (world->hit(r, 0.001, MAXFLOAT, rec)) {
        ray scattered;
        vec3 attenuation;
        vec3 emitted = rec.mat_ptr->emitted(rec.u, rec.v, rec.p);
        float pdf;
        vec3 albedo;
        if (depth < 50 && rec.mat_ptr->scatter(r, rec, albedo, scattered, pdf)) {
        
            //new start
            float sizeX = lightEndX - lightStartX;
            float sizeZ = lightEndZ - lightStartZ;
            float light_area = (lightEndX-lightStartX)*(lightEndZ-lightStartZ);
            vec3 on_light = vec3(lightStartX + random_double()*sizeX, 
                                lightY, 
                                lightStartZ + random_double()*sizeZ);
            vec3 to_light = on_light - rec.p;
            float distance_squared = to_light.squared_length();
            to_light.make_unit_vector();
            if (dot(to_light, rec.normal) < 0) return emitted;
            float light_cosine = fabs(to_light.y());
            if (light_cosine < 0.000001) return emitted;
            pdf = distance_squared / (light_cosine * light_area);
            scattered = ray(rec.p, to_light, r.time());
            //new end
            
            return emitted + albedo*rec.mat_ptr->scattering_pdf(r, rec, scattered)
                *color(scattered, world, depth+1) / pdf;
        }
        else
            return emitted;
    }
    else
        return vec3(0,0,0);
}

渲染结果如下：

上图为每像素采样5次，运行时间与噪点大幅减少，这个效果相当惊艳。这应该是渲染Cornell Box以来最大的突破了，利用重要性采样，将所有Lambertian材质的散射光线指向了光源的方向，减少了大量噪点，且当散射光线直射光源的时候，不再发生散射，大幅减少了射线转折的次数。

不过，天花板下面的的光源周围出现了一些噪点，因为光源是一个两面矩形，距离天花板还有一定的空间。我们可以将其修改为只向下发光：

修改material的虚函数：

class material {
public:
	virtual bool scatter(const ray& r_in,
		const hit_record& rec, vec3& albedo, ray& scattered, float& pdf) const {
		return false;
	}
	virtual float scattering_pdf(const ray& r_in, const hit_record& rec,
		const ray& scattered) const {
		return 0;
	}
	/*virtual vec3 emitted(float u, float v, const vec3& p) const {
		return vec3(0, 0, 0);
	}*/
	virtual vec3 emitted(const ray& r_in, const hit_record& rec, float u, float v, const vec3& p) const {
		return vec3(0, 0, 0);
	}
};

修改color函数：

vec3 color(const ray& r, hittable *world, int depth) {
	hit_record rec;
	if (world->hit(r, 0.001, FLT_MAX, rec)) {
		ray scattered;
		vec3 attenuation;
		//vec3 emitted = rec.mat_ptr->emitted(rec.u, rec.v, rec.p);
		vec3 emitted = rec.mat_ptr->emitted(r, rec, rec.u, rec.v, rec.p);
		float pdf;
		vec3 albedo;
		if (depth < 50 && rec.mat_ptr->scatter(r, rec, albedo, scattered, pdf)) {

			//new start
			float sizeX = lightEndX - lightStartX;
			float sizeZ = lightEndZ - lightStartZ;
			float light_area = (lightEndX - lightStartX)*(lightEndZ - lightStartZ);
			vec3 on_light = vec3(lightStartX + random_double()*sizeX,
				lightY,
				lightStartZ + random_double()*sizeZ);
			vec3 to_light = on_light - rec.p;
			float distance_squared = to_light.squared_length();
			to_light.make_unit_vector();
			if (dot(to_light, rec.normal) < 0) return emitted;
			float light_cosine = fabs(to_light.y());
			if (light_cosine < 0.000001) return emitted;
			pdf = distance_squared / (light_cosine * light_area);
			scattered = ray(rec.p, to_light, r.time());
			//new end

			return emitted + albedo * rec.mat_ptr->scattering_pdf(r, rec, scattered)
				*color(scattered, world, depth + 1) / pdf;
		}
		else
			return emitted;
	}
	else
		return vec3(0, 0, 0);
}

修改diffuse_light：

class diffuse_light : public material {
public:
	texture *emit;
	diffuse_light(texture *a) : emit(a) {}
	virtual bool scatter(const ray& r_in, const hit_record& rec, vec3& attenuation, ray& scattered) const {
		return false;
	}
	virtual vec3 emitted(const ray& r_in, const hit_record& rec, float u, float v, const vec3& p) const {
		
		if (dot(rec.normal, r_in.direction()) >= 0.0)
			return emit->value(u, v, p);
		else
			return vec3(0, 0, 0);
	}
	/*virtual vec3 emitted(float u, float v, const vec3& p) const {
		return emit->value(u, v, p);
	}*/
};

效果如下（采样数50，光源周围的噪点明显减少了）：