PBRT阅读：第十五章 蒙特卡罗积分 2: 提高效率 （第15.6-15.7节）

http://www.opengpu.org/forum.php?mod=viewthread&tid=7364&fromuid=10107

15.6 对光源采样

因为来自光源的直接照射是给定点上的反射光的主要的贡献部分， 所以非常有必要能够对该点附近的直接光照值不为零的地方进行方向采样。考虑一下被一个小球面光源照射下的漫反射表面：如果用BSDF的采样分布进行方向采样，效率就有可能非常低下，因为对这个点而言，光源只在一个小圆锥内的方向上可见。更好的方法则是基于光源的采样分布进行采样。例如，采样例程可以只在小球可能可见的那些方向上采样。这个球可能在这些方向上都被遮住或部分被遮住，但通常不太可能用增加采样密度来解决这个问题。

这一节用我们已经熟知的蒙特卡罗方法来解决这一问题。更进一步地，它显示了对光源上的出射光线的采样技术的推导和实现。这项技术对双向光传输算法（如光子映射和双向路径追踪）是至关重要的。

15.6.1 基本接口

所有Light类必须实现三个主要的跟采样有关的函数。其中两个，Light::Sample_L()和Light::Pdf()，可以跟BxDF::Sample_f()和BxDF::Pdf()相对应；它们也从一个点出发的球面方向分布上生成采样。在这里，这个分布并不试着匹配该点的BSDF，而是去匹配来自光源的直接照射的入射辐射亮度分布。

<Light Interface> +=

virtual Spectrum Sample_L(const Point &p, float u1,

float u2, Vector *wi, float *pdf,

VisibilityTester *vis) const = 0;

<Light Interface> +=

virtual float Pdf(const Point &p, const Vector &wi) const = 0;

这些例程有一些变型使用的是点p上的法向量。那些可以根据投影立体角来采样的光源可以使用这些例程，由于考虑了余弦衰减，就进一步地减小了方差。那些不使用表面法向量的例程对于体积散射是必需的，因为没有表面法向量可用，所以下面的例程不适用于体积散射。它们的缺省实现只是调用那些不使用表面法向量的版本，所以如果光源没有特殊的采样算法就不必再实现它们了。

<Light Interface> +=

virtual Spectrum Sample_L(const Point &p, const Normal &n,

float u1, float u2, Vector * wi , float *pdf,

VisibilityTester * visibility) const {

return Sample L(p, u1, u2, wi , pdf, visibility);

}

<Light Interface> +=

virtual float Pdf(const Point &p, const Normal &n,

const Vector &wi) const {

return Pdf(p, wi);

}

第三个光源采样函数很特别。它是在“离开”光源的光线分布上进行光线采样。这个函数用来支持那些考虑光线路径的算法，例如第16.5节的光子映射。它所返回的PDF值应该用基于光源表面积上的密度和关于立体角的密度的乘积来表达。该函数在本节中的实现将在两个分布上分别采样并计算它们的乘积，而不是采用构造一个联合的4维分布并采样的方法。

<Light Interface> +=

virtual Spectrum Sample_L(const Scene *scene, float u1,

float u2, float u3, float u4,

Ray *ray, float *pdf) const = 0;

15.6.2 带奇点的光源

跟理想镜面反射和透射相似，用delta函数定义的光源也可以很自然地被放入这个采样框架之中，当然调用它们的采样函数的例程要相当地小心，因为它们所返回的辐射亮度和PDF值带有delta分布。在大多数情况下，在做估计值计算时，这些delta分布会很自然地被抵消掉，而在多重重要性采样时就要注意处理这种情况，就象前面的BSDF那样。

点光源

点光源是这样用delta分布来描述的：在一个被照点上，点光源只在一个方向上对其照明。所以，采样问题就变得很简单。这里的在入射照明方向上调用的Sample_L()函数只是回调第13.2节所实现的基本Sample_L()函数。实际上，对于所有用delta分布来描述的光源都可以用这种方法来实现其Sample_L()函数。

<PointLight Method Definitions> +=

Spectrum PointLight::Sample_L(const Point &p,float u1 ,

float u2, Vector *wi, float *pdf,

VisibilityTester * visibility ) const {

*pdf = l.f;

return Sample_L(p, wi , visibility ) ;

}

由于它有delta分布，所以它的Pdf()函数返回0。这是跟这个事实相一致的：即其它的采样例程无法在无限小的光源上随机地采样一个点。

<PointLight Method Definitions> +=

float PointLight::Pdf(const Point &, const Vector &) const {

return 0.;

}

对离开光源的光线进行采样的函数也很简单。光源原点必须在光源位置上，这附近的密度可以用一个delta分布来描述。我们在球面上对方向进行均匀采样，采样的总密度是这两个密度的乘积。跟往常一样，我们将忽略掉不应该包含在PDF的delta分布，因为它会跟采样例程所返回的辐射亮度值中相应的delta项相抵消。

<PointLight Method Definitions> +=

Spectrum PointLight::Sample_L(const Scene *scene, float u1,

float u2, float u3, float u4,

Ray *ray, float *pdf) const {

ray->o = lightPos;

ray->d = UniformSampleSphere(ul, u2);

*pdf = UniformSpherePdf();

return Intensity;

}

聚光灯

SpotLight类也是将请求传递给第13.2.1节的Sample_L()，并且它的Pdf()函数返回0。

在聚光灯上某个合理的分布上对出射光线进行采样会更有趣一些。虽然我们可以象对付点光源那样在球面上做均匀方向采样，但这个分布很可能跟聚光灯的实际光线分布相违背。例如，如果光源的光束角很窄，那么在许多方向采样上光源并没有照明作用。我们应该在一个光源投射光线的方向锥上做均匀采样。虽然采样分布没有考虑光束边界附近的衰减现象，但这在实际应用中并不是个问题。

PDF p(θ,Φ)是可分类的，其中p(Φ) = 1/(2π)。我们只需要找到θ的分布，并在该分布上采样。我们希望在从中心方向到最大光束角度这个范围之内做均匀采样：

所以有 p(θ) = sinθ( 1 - cosθmax)

<MC Function Definitions> +=

COREDLL float UniformConePdf(float cosThetaMax) {

return 1.f / (2.f * M_PI * (1.f - cosThetaMax));

}

然后我们对PDF积分求得CDF，在用下面的采样技术：

cosθ = (1 - ξ) + ξ cosθmax

这里有两个UniformSampleCone()实现了这个采样技术。第一个在围绕(0,0,1)轴进行采样，另一个没有列在这里，它使用

的是采样坐标系的三个坐标轴向量，并按照z轴进行采样。

<MC Function Definitions> +=

Vector UniformSampleCone(float u1 , float u2, float costhetamax) {

float costheta = Lerp(u1, costhetamax, 1.f );

float sintheta = sqrt(l.f - costheta*costheta);

float phi = u2 * 2.f * M_PI;

return Vector(cosf(phi) * sintheta,

sinf(phi) * sintheta,

costheta);

}

光源的光线采样函数将利用这些工具函数在聚光灯的光束圆锥内进行采样。

<SpotLight Method Definitions> +=

Spectrum SpotLight::Sample_L(const Scene *scene, float u1,

float u2, float u3, float u4,

Ray *ray, float *pdf) const {

ray->o = lightPos;

Vector v = UniformSampleCone(u1, u2, cosTotalWidth);

ray->d = LightToWorld(v);

*pdf = UniformConePdf(cosTotalWidth);

return Intensity * Falloff(ray->d);

}

投射光源和配光图光源

投射光源和配光图光源的采样例程在本质上是跟SpotLight和PointLight是一样的。当对出射光进行采样时，ProjectionLight对包住图像贴图的光锥进行均匀采样，而对于配光图光源而言，只需对单位球面进行均匀采样。

方向光源

从一个给定点上的方向光源进行入射光采样是很容易的。就像点光源那样，第12章实现了它的Sample_L()函数，其PDF函数返回0。

从光源的出射光线分布中进行光线采样更有趣一些。光线的方向是由一个delta分布来确定的。它必须是光源的反方向。而它的原点却有无限多个的3D点可供选择。我们如何进行选择并计算它的密度函数呢？

我们想得到的特性是：在场景中光线的交点被远处的光源以均匀的概率所照明。其中一个方法是构造一个跟场景包围球半径相同的圆盘，它的法向量就是光源的方向，然后在圆盘上用ConcertricSampleDisck()随机地选择一个点。一旦选定这个点后，就沿光源方向后移场景包围球半径的距离，并作为光线的原点，这样原点就落在了包围球之外，但却可以跟包围球相交。

这是一个可行的采样方案，因为根据构造过程可知所有进入包围球的光线都有非零概率。采样密度的面积成分是均匀的，因而等于圆盘的面积的倒数。这个密度是由基于光源方向的delta分布而给定的，因而并不包括在返回的PDF值中。

<DistantLight Method Definitions> +=

Spectrum DistantLight::Sample_L(const Scene *scene,

const LightSample &ls, float u1, float u2, float time,

Ray *ray, Normal *Ns, float *pdf) const {

<Choose point on disk oriented toward infinite light direction>

<Set ray origin and direction for infinite light ray>

*Ns = (Normal)ray->d;

*pdf = 1.f / (M_PI * worldRadius * worldRadius);

return L;

}

选择有向圆盘上的点需要使用一点向量代数的知识。我们先用两个跟圆盘法向量相垂直的向量构造一个坐标系(v1,v2,n)，如图。给定一个单位圆盘上的随机点(d1,d2)，就可以就可以计算出它的位置距离圆心的偏移 d1 v1 + d2 v2，从而得到它的位置。

<Choose point on disk oriented toward infinite light direction> =

Point worldCenter;

float worldRadius;

scene->WorldBound().BoundingSphere(&worldCenter, &worldRadius);

Vector v1, v2;

CoordinateSystem(lightDir, &v1, &v2);

float d1, d2;

ConcentricSampleDisk(ls.uPos[0], ls.uPos[1], &d1, &d2);

Point Pdisk = worldCenter + worldRadius * (d1 * v1 + d2 * v2);

最后，将点沿光源方向进行偏移，完成对光线的初始化。

<Set ray origin and direction for infinite light ray> =

ray->o = Pdisk + worldRadius * lightDir;

ray->d = -lightDir;

15.6.3 面积光源

我们知道面光源是根据一个发光的形体来定义的。因此，为了恰当地在这样的光源上进行入射光采样，就需要在形体的表面上进行采样。为此，我们在Shape类中加入采样函数来对表面上的随机点进行采样。这样，AreaLight的采样函数就可以调用这些函数了。

对形体采样

有两个形体采样函数，名称都是Shape::Sample()。第一个函数用某些关于表面面积的采样分布在表面上选择点，并返回所选点的位置和表面法向量。例如，在光子贴图算法中，就可以用这个函数在光源的表面上进行采样。

<Shape Interface> +=

virtual Point Sample(float u1, float u2, Normal *Ns) const {

Severe("Unimplemented Shape::Sample() method called");

return Point();

}

实现这个函数的Shape类几乎总是根据表面面积来做均匀采样。这样一来，我们就可以提供一个缺省的Shape::Pdf()实现，它返回面积的倒数。

<Shape Interface> +=

virtual float Pdf(const Point &Pshape) const {

return 1.f / Area();

}

第二种采样函数使用被积分表面上的点作为参数。这个函数对于光照特别有用，因为调用者可以传入被照明的点，这样一来形体类可以保证只对（对于该点而言）可见的那部分进行采样。它的缺省实现忽略了这个点而只是调用上一个采样函数。

<Shape Interface> +=

virtual Point Sample(const Point &P, float u1, float u2,

Normal *Ns) const {

return Sample(ul, u2, Ns);

}

这两个函数有一个重要区别：第一个使用基于形体表面面积的概率密度函数来生成采样点，而第二个使用的是基于点p的立体角的密度函数。这个区别还是很必要的，因为它可以用于生成入射方向的面光源采样例程，所采样的方向是从光源上的一点p'到给定点p，即p' - p。由于pbrt的积分器在对直接光照积分进行求值时是对起始于p的方向进行积分的，所以用关于p的立体角的采样密度就很方便。因此，第二个Pdf()函数的缺省实现是将面积概率密度转换成立体角概率密度。

<Shape Interface> +=

virtual float Pdf(const Point &p, const Vector &wi) const {

<Intersect sample ray with area light geometry>

<Convert light sample weight to solid angle measure>

return pdf;

}

给定了点p和方向ωi，Pdf()函数确定光线(p, ωi)跟形体是否相交。如果光线跟形体根本没有交点，则相应的概率就是0，因为我们假定在方向采样时先在形体上取点。如果存在交点，就很方便地也得到了其微分几何信息。注意这个求交测试只是在光线跟单个的形体之间进行的，忽略了场景中其它的几何体，所以求交测试还是很快的。

<Intersect sample ray with area light geometry> =

Differential Geometry dgLight;

Ray ray(p, wi);

float thit;

if (!Intersect(ray, &thit, &dgLight)) return 0.;

为了计算关于立体角的PDF值，这个函数先技术关于表面积的PDF。为了将关于面积的密度函数转换成关于立体角的密度函数，需要乘以下列因子：

dωi / dA = r2/cosθo

其中θo是光源上的点到点p的方向跟光源上表面法向量的夹角，r2是两点距离的平方。

<Convert light sample weight to solid angle measure> =

float pdf = DistanceSquared(p, ray(thit)) /

(AbsDot(dgLight.nn, -wi) * Area());

面光源的采样函数

有了前面的采样函数，AreaLight::Sample_L()就很简单了。大部分的工作都交给了Shape类，AreaLight只需计算出射辐射亮度值。

<AreaLight Method Definitions> +=

Spectrum AreaLight::Sample_L(const Point &p, const Normal &n,

float u1, float u2, Vector *wi, float *pdf,

VisibilityTester Visibility) const {

Normal ns;

Point ps = shape->Sample(p, n, u1 , u2, &ns);

*wi = Normalize(ps - p);

*pdf = shape->Pdf(p, *wi);

visibility->SetSegment(p, ps);

return L(p, ns, -*wi);

}

因为我们要调用的Shape::Pdf()返回关于立体角的密度，故而AreaLight::Pdf()可以将其直接返回。

<AreaLight Method Definitions> +=

float AreaLight::Pdf(const Point &p, const Normal &k,

const Vector &wi) const {

return shape->Pdf(p, n, wi);

}

Light::Sample_L()和Light::Pdf()也提供了只传入点而没有传入法向量的接口，跟前面类似，略。

有了形体的采样函数之后，对面光源上出离光线进行采样的函数也是很容易实现的。第一个Shape::Sample()函数用来根据给定的密度函数来找光线原点。由于AreaLight在所有方向上均匀地辐射，所以这里用均匀的方向分布。由于光源只在其表面法向量所在的那一面发光，如果所生成的采样方向在法向量相反的那个半球，就要将它翻转，这样在没有照明的方向上的采样不至于被浪费。

用于光线采样的PDF是进行原点采样的PDF和1/2π（在半球上均匀方向采样的PDF)的乘积。

<AreaLight Method Definitions> +=

Spectrum AreaLight::Sample_L(const Scene *scene, float u1,

float u2, float u3, float u4,

Ray *ray, float *pdf) const {

Normal ns;

ray->o = shape->Sample(u1, u2, &ns);

ray->d = UniformSampleSphere(u3, u4);

if (Dot(ray->d, ns) < 0.) ray->d *= -1;

*pdf = shape->Pdf(ray->o) * INVJTWOPI;

return L(ray->o, ns, ray->d);

}

圆盘采样

Disk采样函数用到了同心圆采样函数来得到单位圆上的一个点，并将其按照给定的圆盘半径和高度进行比例和偏移变换。注意这里没有考虑部分圆盘的情况。

<Disk Public Methods> =

Point Disk::Sample(float u1, float u2, Normal *Ns) const {

Point p;

ConcentricSampleDisk(ul, u2, &p.x, &p.y);

p.x *= radius;

p.y *= radius;

p.z = height;

*Ns = Normalize(ObjectToWorld(Normal(0,0,1)));

if (reverseOrientation) *Ns *= - 1.f ;

return ObjectToWorld(p);

}

圆柱面采样

在圆柱面上进行均匀采样也是很简单的。我们对高度和Φ值进行均匀采样。很直观地我们可以理解这个方法可行，因为圆柱面只是一个卷起来的长方形。

<Cylinder Public Methods> =

Point Cylinder::Sample(float u1, float u2, Normal *Ns) const

float z = Lerp(ul, zmin, zmax);

float t = u2 * phiMax;

Point p = Point(radius * cosf(t) , radius * sinf(t) , z);

*Ns == Normalize(ObjectToWorld(Normal (p.x, p.y, 0.)));

if (reverseOrientation) *Ns *= -1.f;

return ObjectToWorld(p);

}

三角形采样

前一章所定义的UniformSampleTriple()函数返回三角形内一个均匀采样点的重心坐标。利用重心坐标我们可以求出点的位置。

<TriangleMesh Method Definitions> +=

Point Triangle::Sample(float u l , float u2, Normal *Ns) const {

float b1, b2;

UniformSampleTriangle(u1, u2, &bl, &b2);

<Get triangle vertices in pi, p2, and p3>

Point p = b1 * pi + b2 * p2 + ( 1.f - b1 - b2) * p3;

Normal n = Normal(Cross(p2-pl, p3-pl));

*Ns = Normalize(n);

if (reverseOrientation) *Ns *= -1.f;

return p;

}

球面采样

跟Disk一样，这里的球面采样不考虑部分球面的情况。如果函数并不传入被照明的外部点，在球面上采样是非常简单的。它只需用UniformSampleSphere()函数生成单位球面上的点并按照球半径进行比例变换就可以了。

<Sphere Public Methods> =

Point Sphere::Sample(float u1, float u2, Normal *ns) const {

Point p = Point(0,0,0) + radius * UniformSampleSphere(u1, u2);

*ns = Normalize(ObjectToWorld(Normal(p.x, p.y, p.z)));

if (reverseOrientation) *ns *= -1.f;

return ObjectToWorld(p);

}

如果给定了被照明的点，我们就有更好的办法。虽然在球面上均匀采样可以得到正确的估算值，但更好的方法是不要在球面上那些对该点而言不可见的地方进行采样。采样例程要在从被着色点开始的包含球面的立体角区域内进行均匀的方向采样。它首先在这个方向锥内对从中心向量ωc开始的偏移角θ进行采样，然后再对绕该向量的旋转角Φ进行采样。

从被着色点p看球面，它所包含的角度是：

其中r是球半径，pc是球心。这里的实现均匀地在这个方向锥进行采样，然后计算方向采样跟球面的交点来得到光源上的采样位置。

<Sphere Method Definitions> +=

Point Sphere::Sample(const Point &p, float u1, float u2,

Normal *ns) const {

<Compute coordinate system for sphere sampling>

<Sample uniformly on sphere if p is inside it>

<Sample sphere uniformly inside subtended cone>

}

如果我们先计算出一个用于球面采样的坐标系，这个采样例程就会很简单：这个坐标系的z轴是球心到被照明点的向量。

<Compute coordinate system for sphere sampling> =

Point Pcenter = (*ObjectToWorld)(Point(0,0,0));

Vector wc = Normalize(Pcenter - p);

Vector wcX, wcY;

CoordinateSystem(wc, &wcX, &wcY);

我们应该当心那些位于球面之内的点。在这种情况下，就要对整个球面进行采样，因为在球面内整个球面都是可见的。注意我们用了一个很小的常数1e-4来做这项检查，这样可以避免当点很靠近球面时可能发生的错误。

<Sample uniformly on sphere if p is inside it> =

if (DistanceSquared(p, Pcenter) - radius*radius < 1e-4f)

return Sample(ul, u2, ns);

注意我们必须认真对待这里的精度误差问题。如果所生成的光线球面的边缘，那么Sphere::Intersect()就可能意外地返回false。在这种情况下，这里的实现就在被着色点到球心的连线上任取一点并返回。这样做使得采样例程有些偏差，但这所引起的误差还是极其微小的。

<Sample sphere uniformly inside subtended cone> =

float cosThetaMax =

sqrtf(max(0.f, 1.f - radius*radius / DistanceSquared(p, Pcenter)));

Differential Geometry dgSphere;

float thit;

Point ps;

Ray r( p , UniformSampleCone(u1, u2, cosThetaMax, wcX, wcY, wc));

if (!Intersect(r, &thit, &dgSphere))

ps = Pcenter - radius * wc;

else

ps = r(thit) ;

*ns = Normal(Normalize(ps - Pcenter));

if (reverseOrientation) *ns *= -1.f ;

return ps;

<Sphere Public Methods> +=

float Sphere::Pdf(const Point &p, const Vector &wi) const {

Point Pcenter = ObjectToWorld(Point(0,0,0));

<Return uniform weight if point inside sphere>

<Compute general sphere weight>

}

为了计算该采样例程所用的权值，我们必须先区分开球内和球外所使用的采样策略。如果着色点位于球内，就需要使用均匀采样策略。在这种情况下，这里的实现将Pdf()的调用传递给父类，而负类负责立体角的转换。

<Return uniform weight if point inside sphere> =

if (DistanceSquared(p, Pcenter) - radius*radius < le-4f)

return Shape::Pdf(p, wi ) ;

而在一般情况下，我们只要重新计算球面所包含的角度并调用UniformConePdf()。注意这里无需对采样测度进行转换，因为UniformConePdf()已经返回关于立体角的权值了。

<Compute general sphere weight> =

float cosThetaMax =

sqrtf(max(0.f, 1.f - radius*radius / DistanceSquared(p, Pcenter)));

return UniformConePdf(cosThetaMax);

15.6.4 SHAPESET采样

如果给定的形体需要加细成多个形体再进行采样，AreaLight构造器就会申请一个ShapeSet对象。ShapeSet使用跟Shape一样的接口，所以也必须需要实现采样例程。为了能够做到在整个形体上根据表面积进行均匀采样，我们要对形体集合中每个形体进行采样，而概率分布基于其面积跟全部形体面积之比。因此，ShapeSet构造器根据这个比率所确定的单个概率分布来要计算每个形体的离散的CDF。然后采样例程使用这个概率分布随机地选择形体进行采样。这里使用的是线性搜索，当ShapeSet所含的形体数目很小时效率还不错，但对于很大的形体集合而言却效率低下。

<ShapeSet Private Data> +=

float area;

vector<float> areaCDF;

<ShapeSet Public Methods> =

Point Sample(float u1, float u2, Normal *Ns) const {

float ls = RandomFloat();

u_int sn;

for (sn = 0; sn < shapes.size()-1; ++sn)

if (ls < areaCDF[sn]) break;

return shapes[sn]->Sample(u1, u2, Ns);

15.6.5 无限面积光源

InfiniteAreaLight可以被视为一个包含整个场景的无限大的球体，在所有的方向上提供照明。如果我们有被照点所在表面的法向量，就用一个在该点周围的余弦加权分布进行采样；否则，就使用一个均匀球面分布进行采样。这些采样密度函数都可以给出正确的结果，但并没有考虑光源的图像贴图所引起的辐射亮度分布的方向性变化。

<InfiniteAreaLight Method Definitions> +=

Spectrum InfiniteAreaLight::Sample_L(const Point &p,

const Normal &n, float u1, float u2, Vector *wi,

float *pdf, VisibilityTester Visibility) const {

<Sample cosine-weighted direction on unit sphere>

<Compute pdf for cosine-weighted infinite light direction>

<Transform direction to world space>

visibility->SetRay(p, *wi);

return Le(RayDifferential(p, *wi));

}

在单位球面上取一个余弦加权的采样点所用的方法几乎跟第14.5.3节中的Malley法相同，只是我们随机地选择上半球或下半球再进行采样。

<Sample cosine-weighted direction on unit sphere> =

float x, y, z;

ConcentricSampleDisk(u1, u2, &x, &y);

z = sqrtf(max(0.f, 1.f - x*x - y*y));

if (RandomFloatO < .5)

*wi = Vector(x, y, z );

采样的PDF的计算方法跟ConsineHemispherePdf()相似，只是要考虑到整个球体所含的立体角。

<Compute pdf for cosine-weighted infinite light direction> =

*pdf = fabsf(wi->z) * INV_TW0PI;

就像BSDF采样那样，在一个标准坐标系（法向量为(0,0,1)方向）下采样到一个方向以后，需要将被选择的点变换到世界坐标系中，我们可以使用以该法向量为z轴的任意坐标系来做到这一点：

<Transform direction to world space> =

Vector v1, v2;

CoordinateSystem(Normalize(Vector(n)), &v1, &v2);

*wi = Vector(v1.x * wi->x + v2.x * wi->y + n.x * wi->z,

v1.y * wi->x + v2.y * wi->y + n.y * wi->z,

v1.z * wi->x + v2.z * wi->y + n.z * wi->z);

Pdf()函数返回跟InfiniteAreaLight::Sample_L()相同的值：

<InfiniteAreaLight Method Definitions> +=

float InfiniteAreaLight::Pdf(const Point &, const Normal &n,

const Vector &wi) const {

return AbsDot(n, wi) * INV_TWOPI;

}

如果没有提供着色法向量（例如，如果在参与介质中一个点上的采样的情况），就在球面上做均匀的方向采样。方法很简单，从略。

生成从无限光源出发的随机光线有些麻烦，因为我们需要光线方向本身是在整个场景中是均匀分布的。Li（2003）证明了被一个球所包围的区域内的均匀分布的直线可以通过连接两个球面上均匀分布上的点而得到。这里的实现先找到包含整个场景的球面在使用这个方法来生成均匀分布的随机光线。

<InfiniteAreaLight Method Definitions> +=

Spectrum InfiniteAreaLight::Sample_L(const Scene *scene,

float u1, float u2, float u3, float u4,

Ray *ray, float *pdf) const {

<Choose two points p1 and p2 on scene bounding sphere>

<Construct ray between pi and p2>

<Compute InfiniteAreaLight ray weight>

return Le(RayDifferential(ray->o, -ray->d));

}

当然，这个光源的“球面”是包含整个场景的“隐含”的球面。为了使用球面采样例程，还需要显式地找到场景包围球的半径和球心再进行采样。

<Choose two points pi and p2 on scene bounding sphere> =

Point worldCenter;

float worldRadius;

scene->WorldBound().BoundingSphere(&worldCenter, &worldRadius);

worldRadius *= 1.01f;

Point p1 = worldCenter + worldRadius * UniformSampleSphere(u1, u2);

Point p2 = worldCenter + worldRadius * UniformSampleSphere(u3, u4);

一旦选定了场景包围球上的两个点p1,p2，就很容易用这两个点做一条光线：p1为原点，p2-p1为方向向量。

<Construct ray between p1 and p2> =

ray->o = p1;

ray->d = Normalize(p2-p1);

为了计算这些光线的PDF，第一个点的面积分布要乘以第二个点的方向分布。

<Compute InfiniteAreaLight ray weight> =

Vector to_center = Normalize(worldCenter - p1);

float costheta = AbsDot(to_center, ray->d);

*pdf = costheta / (4.f * M_PI * worldRadius * worldRadius);

15.7 体积散射

传输方程是描述在某个环境中参与介质对辐射亮度分布的影响的积分方程。在第17.1节中将介绍这个方程，而且第17章中的VolumeIntegrator类使用蒙特卡罗积分来计算这个积分的估算值。传输方程对于蒙特卡罗积分而言是个很有趣的例子，因为它的积分范围是无限的，其基本形式为：

由于它有指数项，应用于体积散射的这个积分具有有限值。

我们不可能从均匀分布中采样来估算这个积分，因为不存在范围是[0, ∞)的均匀PDF，故而必须使用其它形式的采样分布。在某些情况下，函数g(t)是常数，这时一个很自然的选择是使用基于指数项的采样分布。第13.3.1节中有一般采样过程的推导。如果从这个分布中取采样Ti，那么积分的估算公式为：

对于g(t)不是常量的更一般的情况，就可以任选一个常量c（例如可以基于g(t)的平均值），估算公式为：

因为在pbrt中VolumeRegion的范围是有限的，系统实际上并不需要在无限区域上估算这些积分，所以我们仍可以使用均匀采样分布。然而，当VolumeRegion的范围增大并且g(t)项增大时，均匀分布是很低效的，因为许多采样取自那些指数函数值很小的地方，没有使用基于指数项的重要性采样的方差就会增大。

15.7.1对相函数采样

在某些应用中，对用Henyey-Greenstein相函数所描述的分布进行采样是很有用处的，这里的应用可以是用来计算参与介质中直接光照的多重重要性采样，还有计算参与介质中多重散射效果的更一般性的算法。这个分布的PDF可以分离为θ分量和Φ分量，其中p(Φ) = 1/(2π)。θ的分布如下：

如果g不等于0；否则，cosθ = 1 -2ξ。

这里的例程实现用一个向量ω，一个非对称参数g，还有两个随机数作为参数，然后从这个分布中采样，并且利用所得的结果构造出相应的出射光线。

<Monte Carlo Function Definitions> +=

Vector SampleHG(const Vector &w, float g, float u1, float u2) {

float costheta;

if (fabsf(g) < 1e-3)

costheta = 1.f - 2.f * u1;

else {

float sqrTerm = (1.f - g * g) /

(1.f - g + 2.f * g * u1);

costheta = (1.f + g * g - sqrTerm * sqrTerm) / (2.f * g);

}

float sintheta = sqrtf(max(0.f, 1.f-costheta*costheta));

float phi = 2.f * M_PI * u2;

Vector v1, v2;

CoordinateSystem(w, &v1, &v2);

return SphericalDirection(sintheta, costheta, phi, v1, v2, w);

}

因为相函数已经是分布函数了，所以对给定方向进行采样的概率密度是相函数在给定方向对上的值。

<Monte Carlo Function Definitions> +=

float HGPdf(const Vector &w, const Vector &wp, float g) {

return PhaseHG(w, wp, g);

}

15.7.2 计算光学厚度

第12.1.3节介绍了光学厚度τ，它可以作为一条光线穿过参与介质的总密度的测量。对于有恒定散射性质的介质而言，光学厚度可以用Beer定律的封闭形式来计算。这里我们要实现一般性的DensityRegion::tau()函数，它计算光线穿过任意介质的光学厚度的估算值。

回忆一下τ的定义：

其估算公式为：

其中随机变量Ti取自于某个分布p。一个很自然的方法是使用分层采样，将从0到d的直线分成N层，每层中放入一个随机采样。回忆一下蒙特卡罗跟标准数值积分相比的强项在于它可以处理高维积分和不连续的被积函数。这里我们有一维积分并且有一个变化很光滑的被积函数。Pauly和他的合作者为此找到了更有效的积分技术，即生成一个分层模式，其中每层的采样有相同的偏移量：

Ti = (ξ+i)/dN

其中均匀随机变量ξ用于所有的采样。其结果就是相邻层的采样不会像普通分层采样那样挤在一起，从而减小了方差。

这里关于Tau()函数的实现就使用了这项技术，所以我们只用了一个随机变量u。而且这里间接地选定了采样个数N：调用者要传入stepSize，即每层的大小，该函数会求得光线贯穿介质的长度d，N就可以用下式求得：

N = d / stepSize

下面的实现并没有直接计算N，实际上，while循环为每个采样循环一次，直到点t0离开这个区域。

<Volume Scattering Definitions> +=

Spectrum DensityRegion::tau(const Ray &r, float stepSize,

float u) const {

float t0, t1;

float length = r.d.Length();

if (length == 0.f) return 0.f;

Ray rn(r.o, r.d / length, r.mint * length, r.maxt * length, r.time);

if (!IntersectP(rn, &t0, &t1)) return 0.;

Spectrum tau(0.);

t0 += u * stepSize;

while (t0 < t1) {

tau += sigma_t(rn(t0), -rn.d, r.time);

t0 += stepSize;

}

return tau * stepSize;

}