↓ 跳过正文
TSSSNI’S BLOG

pbrt-v4 Ep. VI: 几何形状

·18718 字·38 分钟
Graphics Rendering Pbrt
目录

pbrt通过Shape抽象出光线相交, 包围盒等接口, 其余形状无关的功能由Primitive封装. 本章主要介绍Shape.

基础接口
#

包围结构
#

Shape中Bounds()接口返回包围盒, NormalBounds()返回法线方向的包围锥.

光线-包围结构相交
#

Bounds3<T>::IntersectP接口提供该功能. 光线与各个轴上的近平面与远平面相交, 三轴近平面上最大的tt即为tneart_{near}, 同理可得tfart_{far}, 若tfar<tneart_{far} < t_{near}则与包围盒不相交. 为提高效率, pbrt支持传入光线方向倒数与表示各个轴方向是否为负的向量, 避免重复计算. 对于NaN, 由于任何NaN参与的逻辑运算都为否, 这里不需要特殊处理.

相交测试
#

Shape的派生类需要实现Intersect接口, 传入的光线位于渲染空间, 返回ShapeIntersection, 代表最近的相交点. 同样的派生类实现IntersectP用于判断是否相交而非具体相交细节.

// ShapeIntersection Definition
struct ShapeIntersection {
    SurfaceInteraction intr;
    Float tHit;
    std::string ToString() const;
};

相交坐标空间
#

大部分形状的相交在渲染空间中计算, 部分形状需要在本地空间表示, 例如球. 变换不会影响相交的tt.

面判断
#

光栅化中通常只考虑物体的正面, pbrt中为避免光追失效正反面都需要判断相交.

面积
#

派生类需要实现Area接口来提供面积信息, 用于将Shape作为面积光源.

采样
#

采样返回当前采样点的几何信息.

// ShapeSample Definition
struct ShapeSample {
    Interaction intr;
    Float pdf;
    std::string ToString() const;
};

Shape的采样可以采用类似于二维概率密度函数的形式, 调用者需要保证传入的坐标在形状范围内.

PBRT_CPU_GPU inline pstd::optional<ShapeSample> Sample(Point2f u) const;

PBRT_CPU_GPU inline Float PDF(const Interaction &) const;

同样也可以提供参考点, 例如将Shape作为面积光照, 某点需要在Shape上采样来计算辐射亮度, 此时提供参考点或立体角的采样会更加合适. 采样的实现可能需要判断当前立体角对应的光线是否相交.

PBRT_CPU_GPU inline pstd::optional<ShapeSample> Sample(const ShapeSampleContext &ctx,
                                                           Point2f u) const;
PBRT_CPU_GPU inline Float PDF(const ShapeSampleContext &ctx, Vector3f wi) const;

球体
#

球面坐标可以转化为uv方程以用于纹理映射, 也可以用于表示不完整的球.

ϕ=uϕmax⁡θ=θmin⁡+v(θmax⁡−θmin⁡) \begin{equation} \begin{aligned} \phi &= u \phi_{\max}\\ \theta &= \theta_{\min} + v(\theta_{\max} - \theta_{\min}) \end{aligned} \end{equation}

由于球体使用本地空间, 构造函数需要传入变换. pbrt没有使用添加动画的变换, 动画由Primitive处理.

包围结构
#

由于球体可能是不完整的, pbrt会计算出对应的zmin⁡z_{\min}和zmax⁡z_{\max}, x轴与y轴不做额外计算.

PBRT_CPU_GPU Bounds3f Sphere::Bounds() const {
    return (*renderFromObject)(
        Bounds3f(Point3f(-radius, -radius, zMin), Point3f(radius, radius, zMax)));
}

表面面积
#

定义在[a,b][a,b]上的曲线f(x)f(x)绕x轴形成的形状的表面积可以表示为下式.

A=ϕmax⁡∫abf(x)ds=ϕmax⁡∫abf(x)dx2+dy2=ϕmax⁡∫abf(x)1+f′2(x)dx \begin{equation} \begin{aligned} A &= \phi_{\max} \int_a^b f(x) \mathrm{d}s\\ &= \phi_{\max} \int_a^b f(x) \sqrt{\mathrm{d}x^2 + \mathrm{d}y^2}\\ &= \phi_{\max} \int_a^b f(x) \sqrt{1 + f'^2(x)} \mathrm{d}x \end{aligned} \end{equation}

将球面视为绕z轴旋转, 可以得到下式.

A=ϕmax⁡∫zmin⁡zmax⁡r2−z21+z2r2−z2dz=ϕmax⁡r(zmax⁡−zmin⁡) \begin{equation} \begin{aligned} A &= \phi_{\max} \int_{z_{\min}}^{z_{\max}} \sqrt{r^2 - z^2} \sqrt{1 + \frac{z^2}{r^2 - z^2}} \mathrm{d}z\\ &= \phi_{\max} r (z_{\max} - z{\min}) \end{aligned} \end{equation}

相交测试
#

球体首先执行BasicIntersect, 获取本地空间下的相交结果, 返回QuadricIntersection, 别的二次曲面形状也会使用这个结构体. 之后再执行InteractionFromIntersection, 这样分离使得IntersectP只需要执行BasicIntersect, 且在GPU上可以先对所有可相交的形状执行BasicIntersect, 找到最近点后在执行InteractionFromIntersection.

struct QuadricIntersection {
    Float tHit;
    Point3f pObj;
    Float phi;
};

通过相交点计算出uv后pbrt会计算偏导数, 空间位置的偏导数根据球面坐标的定义计算, 法线的偏导数通过微分几何中常用的Weingarten方程计算.

上述过程执行完成后即可构造SurfaceInteraction.

bool flipNormal = reverseOrientation ^ transformSwapsHandedness;
Vector3f woObject = (*objectFromRender)(wo);
return (*renderFromObject)(SurfaceInteraction(Point3fi(pHit, pError),
                                              Point2f(u, v), woObject, dpdu, dpdv,
                                              dndu, dndv, time, flipNormal));

表面采样
#

本地空间下的球面均匀采样较为简单, pdf为面积的倒数. 对于有参考点的采样, pbrt会将采样点朝向球心偏移, 以避免位于球面附近的点因为浮点误差而位于错误的一侧, 同时pdf会从面积微分转为立体角微分.

Point3f pCenter = (*renderFromObject)(Point3f(0, 0, 0));
Point3f pOrigin = ctx.OffsetRayOrigin(pCenter);

球外的参考点会与球面形成一个可见锥体, 与参考点到球心形成的中心向量可以形成的最大角度见下式. 此时可以在可见锥体内均匀采样立体角方向, 获取与中心向量的夹角θ\theta和绕中心向量的旋转角ω\omega即可. pbrt在[cos⁡θmax,1][\cos\theta_{max}, 1]这一范围内采样.

θmax⁡=arcsin(r∣p−pc∣) \begin{equation} \theta_{\max} = \text{arcsin}(\frac{r}{|p-p_c|}) \end{equation}

对于θmax\theta_{max}较小的情况, 由cos⁡θ\cos\theta计算sin⁡θ\sin\theta会导致浮点误差, 因为cos⁡θ\cos\theta接近1, 1附近的浮点精度与0附近相比是不够的. 在角度较小的情况下根据sin⁡θ\sin\theta的一阶泰勒展开可以得到sin⁡θ≈θ\sin\theta \approx \theta, 此时可以得到足够的精度. 采样完成后可以计算得到交点.

if (sin2ThetaMax < 0.00068523f /* sin^2(1.5 deg) */) {
    // Compute cone sample via Taylor series expansion for small angles
    sin2Theta = sin2ThetaMax * u[0];
    cosTheta = std::sqrt(1 - sin2Theta);
    oneMinusCosThetaMax = sin2ThetaMax / 2;
}

对于pdf, 首先需要满足如下的关系.

1=∫Φp(ω)dω=∫02π∫0θmax⁡p(ϕ)p(θ)sin⁡θdθdω \begin{equation} \begin{aligned} 1 &= \int_{\Phi} p(\omega) \mathrm{d}\omega\\ &= \int_0^{2\pi} \int_0^{\theta_{\max}} p(\phi)p(\theta) \sin\theta \mathrm{d}\theta \mathrm{d}\omega \end{aligned} \end{equation}

已知p(ϕ)=12πp(\phi)=\frac{1}{2\pi}, p(θ)p(\theta)为均匀分布, 则有如下关系.

1=c∫0θmax⁡sin⁡θdθ=c(1−cos⁡θmax⁡) \begin{equation} \begin{aligned} 1 &= c \int_0^{\theta_{\max}} \sin\theta \mathrm{d}\theta\\ &= c (1 - \cos\theta_{\max}) \end{aligned} \end{equation}

此时可以得到最终的pdf.

p(ω)=12π(1−cos⁡θmax⁡) \begin{equation} p(\omega) = \frac{1}{2\pi(1 - \cos\theta_{\max})} \end{equation}

圆柱
#

圆柱的参数方程如下, z是圆柱的高度.

ϕ=uϕmax⁡x=rcos⁡ϕy=rsin⁡ϕz=zmin⁡+v(zmax⁡−zmin⁡) \begin{equation} \begin{aligned} \phi &= u \phi_{\max}\\ x &= r \cos\phi\\ y &= r \sin\phi\\ z &= z_{\min} + v(z_{\max} - z_{\min}) \end{aligned} \end{equation}

包围结构
# 

PBRT_CPU_GPU Bounds3f Cylinder::Bounds() const {
    return (*renderFromObject)(
        Bounds3f({-radius, -radius, zMin}, {radius, radius, zMax}));
}

表面面积
# 

PBRT_CPU_GPU
Float Area() const { return (zMax - zMin) * radius * phiMax; }

相交测试
#

通过x2+y2=r2x^2 + y^2 = r^2计算相交, u=ϕϕmax⁡u = \frac{\phi}{\phi_{\max}}, v=z−zmin⁡zmax⁡−zmin⁡v = \frac{z - z_{\min}}{z_{\max} - z_{\min}}.

表面采样
#

圆柱体没有特定的根据参考点的立体角采样的方法, 而是先进行面积采样再转为立体角, 根据立体角与面积的转换公式可以获取pdf.

p(ω)=1Ar2cos⁡θ \begin{equation} p(\omega) = \frac{1}{A} \frac{r^2}{\cos\theta} \end{equation}

圆盘
#

圆盘参数方程如下, 本地空间中圆盘位于xy平面, rir_i代表inner radius.

ϕ=uϕmax⁡x=((1−v)r+vri)cos⁡ϕy=((1−v)r+vri)sin⁡ϕz=h \begin{equation} \begin{aligned} \phi &= u \phi_{\max}\\ x &= ((1-v)r + vr_i)\cos\phi\\ y &= ((1-v)r + vr_i)\sin\phi\\ z &= h \end{aligned} \end{equation}

包围结构
# 

PBRT_CPU_GPU Bounds3f Disk::Bounds() const {
    return (*renderFromObject)(
        Bounds3f(Point3f(-radius, -radius, height), Point3f(radius, radius, height)));
}

PBRT_CPU_GPU DirectionCone Disk::NormalBounds() const {
    Normal3f n = (*renderFromObject)(Normal3f(0, 0, 1));
    if (reverseOrientation)
        n = -n;
    return DirectionCone(Vector3f(n));
}

表面面积
#

A=ϕmax⁡2(r2−ri2) \begin{equation} A = \frac{\phi_{\max}}{2}(r^2 - r^2_i) \end{equation}

相交测试
#

计算出z位于h时光线对应的t值以及判断该点是否位于圆盘内部即可, u=ϕϕmax⁡,v=r−rhitr−riu=\frac{\phi}{\phi_{\max}}, v=\frac{r - r_{\text{hit}}}{r - r_i}.

表面采样
#

对于圆盘采样, 如果在rr和θ\theta上分别均匀采样, 由于半径的增长导致ϕ\phi采样点对应的位置越来越分散, 采样点会集中在圆心附近.

我们需要满足p(x,y)=1πp(x,y)=\frac{1}{\pi}, 即p(r,θ)=rπp(r,\theta)=\frac{r}{\pi}, 此时可以得到如下关系以及根据逆变换法得到的采样方程.

p(r)=∫02πp(r,θ)dθ=2rp(θ∣r)=p(r,θ)p(r)=12π \begin{equation} \begin{aligned} p(r) &= \int_0^{2\pi} p(r,\theta) \mathrm{d}\theta = 2r\\ p(\theta|r) &= \frac{p(r, \theta)}{p(r)} = \frac{1}{2\pi} \end{aligned} \end{equation} r=ϵ1θ=2πϵ2 \begin{equation} \begin{aligned} r &= \sqrt{\epsilon_1}\\ \theta &= 2\pi\epsilon_2 \end{aligned} \end{equation}

此时由于rr的非线性, 导致u,v对应的面积的不均匀, 这影响了分层抽样的效果. pbrt采用同心映射将[−1,1][-1,1]对应的xy平面映射到圆上. 下式分别为∣x∣>∣y∣|x|>|y|与∣y∣≥∣x∣|y| \ge |x|的情况.

r=xθ=π4yx \begin{equation} \begin{aligned} r &= x\\ \theta &= \frac{\pi}{4} \frac{y}{x} \end{aligned} \end{equation} r=yθ=π2−π4xy \begin{equation} \begin{aligned} r &= y\\ \theta &= \frac{\pi}{2} - \frac{\pi}{4} \frac{x} {y} \end{aligned} \end{equation}

三角网格
#

根据Euler–Poincaré公式, 闭合几何体的顶点数VV, 边数EE与面数FF满足如下关系, 其中gg代表几何体上洞的数量.

V−E+F=2(1−g) \begin{equation} V - E + F = 2 (1 - g) \end{equation}

对于三角形网格还有如下关系, 面数较多时gg可以忽略不计.

E=32FF≈2V \begin{equation} \begin{aligned} E &= \frac{3}{2} F\\ F &\approx 2 V \end{aligned} \end{equation}

网格存储
#

TriangleMesh类的声明如下. 尽管构造函数传入的是vector, 内部仍然使用指针来存储数据, 以此来达到多个共享几何数据的效果. 尽管pbrt提供了instance功能, 像高精度地形的存储仍然需要这种方式来节省内存.

class TriangleMesh {
  public:
    // TriangleMesh Public Methods
    TriangleMesh(const Transform &renderFromObject, bool reverseOrientation,
                 std::vector<int> vertexIndices, std::vector<Point3f> p,
                 std::vector<Vector3f> S, std::vector<Normal3f> N,
                 std::vector<Point2f> uv, std::vector<int> faceIndices, Allocator alloc);

    std::string ToString() const;

    bool WritePLY(std::string filename) const;

    static void Init(Allocator alloc);

    // TriangleMesh Public Members
    int nTriangles, nVertices;
    const int *vertexIndices = nullptr;
    const Point3f *p = nullptr;
    const Normal3f *n = nullptr;
    const Vector3f *s = nullptr;
    const Point2f *uv = nullptr;
    const int *faceIndices = nullptr;
    bool reverseOrientation, transformSwapsHandedness;
};

pbrt通过BufferCache类来将vector转为指针, 它会查找cache中是否已经有一份相同的数据并返回它的指针. 为了防止mutex导致的低并发效率, pbrt根据数据的hash将它存储在不同的cache中. pbrt中通过全局对象来访问BufferCache.

const T *LookupOrAdd(pstd::span<const T> buf, Allocator alloc) {
    ++nBufferCacheLookups;
    // Return pointer to data if _buf_ contents are already in the cache
    Buffer lookupBuffer(buf.data(), buf.size());
    int shardIndex = uint32_t(lookupBuffer.hash) >> (32 - logShards);
    DCHECK(shardIndex >= 0 && shardIndex < nShards);
    mutex[shardIndex].lock_shared();
    if (auto iter = cache[shardIndex].find(lookupBuffer);
        iter != cache[shardIndex].end()) {
        const T *ptr = iter->ptr;
        mutex[shardIndex].unlock_shared();
        DCHECK(std::memcmp(buf.data(), iter->ptr, buf.size() * sizeof(T)) == 0);
        ++nBufferCacheHits;
        redundantBufferBytes += buf.size() * sizeof(T);
        return ptr;
    }

    // Add _buf_ contents to cache and return pointer to cached copy
    mutex[shardIndex].unlock_shared();
    T *ptr = alloc.allocate_object<T>(buf.size());
    std::copy(buf.begin(), buf.end(), ptr);
    bytesUsed += buf.size() * sizeof(T);
    mutex[shardIndex].lock();
    // Handle the case of another thread adding the buffer first
    if (auto iter = cache[shardIndex].find(lookupBuffer);
        iter != cache[shardIndex].end()) {
        const T *cachePtr = iter->ptr;
        mutex[shardIndex].unlock();
        alloc.deallocate_object(ptr, buf.size());
        ++nBufferCacheHits;
        redundantBufferBytes += buf.size() * sizeof(T);
        return cachePtr;
    }

    cache[shardIndex].insert(Buffer(ptr, buf.size()));
    mutex[shardIndex].unlock();
    return ptr;
}

Buffer类用于存储指针并在构造时计算hash.

struct Buffer {
    // BufferCache::Buffer Public Methods
    Buffer() = default;
    Buffer(const T *ptr, size_t size) : ptr(ptr), size(size) {
        hash = HashBuffer(ptr, size);
    }

    bool operator==(const Buffer &b) const {
        return size == b.size && hash == b.hash &&
                std::memcmp(ptr, b.ptr, size * sizeof(T)) == 0;
    }

    const T *ptr = nullptr;
    size_t size = 0, hash;
}

对于三角网格pbrt会在构造时将顶点位置, 法线等属性转到渲染空间中, 然后再存储到cache, 以此节省后续访问时的计算开销. 显然cache miss会增加, 但是计算开销减小了.

三角形类
#

为节省空间开销, Triangle类只存储对应的TriangleMesh的序号以及自身在网格中的序号.

Triangle() = default;
Triangle(int meshIndex, int triIndex) : meshIndex(meshIndex), triIndex(triIndex) {}

NormalBounds计算三角面的几何法线, 如果几何法线与顶点法线的插值结果不指向同一侧, 几何法线会被调整.

PBRT_CPU_GPU DirectionCone Triangle::NormalBounds() const {
    // Get triangle vertices in _p0_, _p1_, and _p2_
    const TriangleMesh *mesh = GetMesh();
    const int *v = &mesh->vertexIndices[3 * triIndex];
    Point3f p0 = mesh->p[v[0]], p1 = mesh->p[v[1]], p2 = mesh->p[v[2]];

    Normal3f n = Normalize(Normal3f(Cross(p1 - p0, p2 - p0)));
    // Ensure correct orientation of geometric normal for normal bounds
    if (mesh->n) {
        Normal3f ns(mesh->n[v[0]] + mesh->n[v[1]] + mesh->n[v[2]]);
        n = FaceForward(n, ns);
    } else if (mesh->reverseOrientation ^ mesh->transformSwapsHandedness)
        n *= -1;

    return DirectionCone(Vector3f(n));
}

相交测试
#

三角形的光线测试是一个独立函数, 而非实现了Shape的接口. pbrt首先变换到以光线起始点为原点, 光线方向为zz轴的空间, 这样可以避免恰好打到边上的光线被认为不相交, 后续章节会解释. 由于浮点精度问题, 需要判断转换后的三角形是否退化为线段.

变换分为三步, 首先平移变换将光线起始点移到原点, 其次将绝对值最大的轴置换为zz轴以保证zz非0, 最后通过切变变换将光线方向与z轴对齐. 切变变换定义如下, 用于将光线方向中xx,yy置0, zz置1, 与旋转相比切变的开销较小. 经过变换后, 只需要判断(0,0)(0, 0)是否位于三角形在xyxy平面的投影上.

S=(10−dxdz001−dydz0001dz00001) \begin{equation} \begin{aligned} S = \begin{pmatrix} 1 & 0 & -\frac{\mathbf{d}_x}{\mathbf{d}_z} & 0\\ 0 & 1 & -\frac{\mathbf{d}_y}{\mathbf{d}_z} & 0\\ 0 & 0 & \frac{1}{\mathbf{d}_z} & 0\\ 0 & 0 & 0 & 1 \end{pmatrix} \end{aligned} \end{equation}

叉乘结果具有方向, 利用行列式可以构造边方程, 用于判断判断点pp在边p0p1p_0p_1的哪一侧. 若某点位于三角形三条边的同侧, 可以认为它位于三角形内部. 由于行列式也可以用于计算两条边组成的平行四边形的面积, 边方程同时计算了重心坐标.

e(p)=∣px−p0xpy−p0yp1x−p0xp1y−p0y∣ \begin{equation} \begin{aligned} e(p) = \begin{vmatrix} p_x - p_{0x} & p_y - p_{0y}\\ p_{1x} - p_{0x} & p_{1y} - p_{0y} \end{vmatrix} \end{aligned} \end{equation}

为避免除以分母带来的误差, pbrt计算交点时不使用面积归一化边方程结果. 经过变换后zz与tt相同, 相交测试完成后可计算重心坐标与相交的tt值.

// Compute scaled hit distance to triangle and test against ray $t$ range
p0t.z *= Sz;
p1t.z *= Sz;
p2t.z *= Sz;
Float tScaled = e0 * p0t.z + e1 * p1t.z + e2 * p2t.z;
if (det < 0 && (tScaled >= 0 || tScaled < tMax * det))
    return {};
else if (det > 0 && (tScaled <= 0 || tScaled > tMax * det))
    return {};

// Compute barycentric coordinates and $t$ value for triangle intersection
Float invDet = 1 / det;
Float b0 = e0 * invDet, b1 = e1 * invDet, b2 = e2 * invDet;
Float t = tScaled * invDet;

三角形内部各个值根据重心坐标线性插值, 包括uvuv也是线性变化, 因此各个值在uvuv上的偏导数在三角形内部相同. 因此求解线性方程即可获取位置, 法线等关于uvuv的偏导数. 位置偏导数主要用于法线纹理, 法线偏导数则用于视差纹理.

(u0−u2v0−v2u1−u2v1−v2)(∂p∂u∂p∂v)=(p0−p2p1−p2) \begin{equation} \begin{aligned} \begin{pmatrix} u_0 - u_2 & v_0 - v_2\\ u_1 - u_2 & v_1 - v_2 \end{pmatrix} \begin{pmatrix} \frac{\partial p}{\partial u}\\ \frac{\partial p}{\partial v} \end{pmatrix} = \begin{pmatrix} p_0 - p_2\\ p_1 - p_2 \end{pmatrix} \end{aligned} \end{equation}

pbrt每次相交都会重新计算偏导数, 以时间换空间. 如果没有uvuv值pbrt会使用默认值. SurfaceInteraction将几何法线初始化为p在uvuv偏导的叉乘, 三角网格由于贴图uvuv的缘故无法采用这种方式, pbrt采用边的叉乘来得到几何边线. 若顶点拥有法线与切线, pbrt会采用插值结果, 且采用类似的方式计算法线在uvuv的偏导.

表面采样
#

采样通过将均匀采样结果转化为重心坐标来实现, 因此在任意形状的三角形上采样均不影响结果. 通过将正方形中关于对角线对称的采样结果视为同一个点可以得到均匀的采样分布, 但这会使得原本相距较远的采样点变为同一个采样点, 这会影响部分采样器的效果.

pbrt采用如下转换方式. 更直观地说, 这将y=1−xy=1-x,x=0x=0和y=0y=0围成的三角形根据y=xy=x分割, 每部分映射到y=xy=x,y=1−xy=1-x和x=0x=0或y=0y=0围成的三角形中. 此时Jacobi行列式为常数, 所以具有面积保持的特征.

f(x,y)=(x−δ,y−δ),δ={x2x<yy2y≤x \begin{equation} \begin{aligned} f(x, y) = (x - \delta, y - \delta), \delta = \begin{cases} \frac{x}{2} & x < y\\ \frac{y}{2} & y \le x \end{cases} \end{aligned} \end{equation}

假设三角形光源各个点发射相同的光, BSDF为常数值, 在三角形面积上均匀采样, 利用面积与立体角微分的转换可以得到如下的采样公式, 其中VV为可见性方程, pp为观察点位置, p′p'为光源上的采样点的位置. 由于分母上的平方项, 距离光源较近的物体会有较大的采样误差, 因此直接将立体角采样转为面积采样并不合适. 对于立体角范围过小和过大的三角形, 由于浮点误差pbrt仍然采用面积采样.

ρLe1/A(V(p,p′)∣cos⁡θ′∣∣cos⁡θ1∣∥p′−p∥2) \begin{equation} \frac{\rho L_e}{1/A}\left(V(p, p')|\cos\theta'|\frac{|\cos\theta_1|}{\Vert p' - p \Vert^2} \right) \end{equation}

对于其他情况, pbrt使用Stratified Sampling of Spherical Triangles提出的球面三角形采样. 首先根据参考点构建单位球, 将三角形投影到单位球上形成球面三角形abcabc, 在abcabc上执行均匀面积采样. 各个边均可与球心构成平面, 计算平面法线, 各个面法线的夹角的补角为平面夹角, 等价于球面三角形的内角, 使用α\alpha, β\beta, γ\gamma分别表示.

实际采样算法如下: 令abcabc面积为A, 首先根据ϵ0\epsilon_0获取在ac⌢\overset{\frown}{ac}上的点c′c', 使形成的三角形abc′abc'面积为ϵ0A\epsilon_0 A. 然后令在c′c'处与bc′⌢\overset{\frown}{bc'}垂直的微分底边与bb形成的微分三角形面积为Abc′⌢A_{\overset{\frown}{bc'}}, 在bc′⌢\overset{\frown}{bc'}上根据ϵ1\epsilon_1采样, 获取c′′c''使得构成的微分三角形Abc′′⌢A_{\overset{\frown}{bc''}}面积为ϵ1Abc′⌢\epsilon_1A_{\overset{\frown}{bc'}}.

由于Aπ=α+β+γ=A+πA_\pi=\alpha+\beta+\gamma=A+\pi, pbrt直接使用AπA_\pi执行采样, 以减少小三角形计算Aπ−πA_\pi-\pi导致的浮点误差, 这通过均匀采样并在π\pi与AπA_\pi间插值实现.

abc′abc'的内角为α\alpha,β′\beta',γ′\gamma', α\alpha与Aπ′A'_\pi已知, β′\beta'与γ′\gamma'未知. 求解c′c'需求解ac′⌢\overset{\frown}{ac'}的长度, 此时根据球面三角公式可得如下关系.

cos⁡β′=−cos⁡γ′cos⁡α+sin⁡γ′sin⁡αcos⁡ac′⌢cos⁡ac′⌢=cos⁡β′+cos⁡(Aπ′−α−β′)cos⁡αsin⁡(Aπ′−α−β′)sin⁡α=cos⁡β′+cos⁡(ϕ−β′)cos⁡αsin⁡(ϕ−β′)sin⁡α=cos⁡β′+(cos⁡ϕcos⁡β′+sin⁡ϕsin⁡β′)cos⁡α(sin⁡ϕcos⁡β′−cos⁡ϕsin⁡β′)sin⁡α \begin{equation} \begin{aligned} \cos\beta'&=-\cos\gamma'\cos\alpha+\sin\gamma'\sin\alpha\cos\overset{\frown}{ac'}\\ \cos\overset{\frown}{ac'} &=\frac{\cos\beta'+\cos(A'_\pi-\alpha-\beta')\cos\alpha}{\sin(A'_\pi-\alpha-\beta')\sin\alpha}\\ &=\frac{\cos\beta'+\cos(\phi-\beta')\cos\alpha}{\sin(\phi-\beta')\sin\alpha}\\ &=\frac{\cos\beta'+(\cos\phi\cos\beta' + \sin\phi\sin\beta')\cos\alpha}{(\sin\phi\cos\beta'-\cos\phi\sin\beta')\sin\alpha}\\ \end{aligned} \end{equation}

此时β′\beta'仍然未知, 引入已知量ab⌢\overset{\frown}{ab}, 可得如下关系, 其中k1k_1,k2k_2皆为已知量.

cos⁡(ϕ−β′)=−cos⁡β′cos⁡α+sin⁡β′sin⁡αcos⁡ab⌢0=(cos⁡ϕ+cos⁡α)cos⁡β′+(sin⁡ϕ−sin⁡αcos⁡ab⌢)sin⁡β′=k1cos⁡β′+k2sin⁡β′cos⁡β′=±k2k12+k22sin⁡β′=∓k1k12+k22 \begin{equation} \begin{aligned} \cos(\phi-\beta')&=-\cos\beta'\cos\alpha+\sin\beta'\sin\alpha\cos\overset{\frown}{ab}\\ 0 &=(\cos\phi+\cos\alpha)\cos\beta'+(\sin\phi-\sin\alpha\cos\overset{\frown}{ab})\sin\beta'\\ &=k_1\cos\beta'+k_2\sin\beta'\\ \cos\beta'&=\frac{\pm k_2}{\sqrt{k_1^2+k_2^2}}\\ \sin\beta'&=\frac{\mp k_1}{\sqrt{k_1^2+k_2^2}} \end{aligned} \end{equation}

将上式代入, cos⁡ac′⌢\cos\overset{\frown}{ac'}求解结果如下, 获取ac′⌢\overset{\frown}{ac'}后即可计算c′c'的实际位置.

cos⁡ac′⌢=k2+(k2cos⁡ϕ−k1sin⁡ϕ)cos⁡α(k2sin⁡ϕ+k1cos⁡ϕ)sin⁡αc′=cos⁡ac′⌢a+sin⁡ac′⌢c⊥ \begin{equation} \begin{aligned} \cos\overset{\frown}{ac'}&=\frac{k_2+(k_2\cos\phi-k_1\sin\phi)\cos\alpha}{(k_2\sin\phi+k_1\cos\phi)\sin\alpha}\\ \mathbf{c'}&=\cos\overset{\frown}{ac'}\mathbf{a}+\sin\overset{\frown}{ac'}\mathbf{c_\perp} \end{aligned} \end{equation}

将bb变换到球的北极, 根据微分立体角转微分面积的公式, 以及经过bb的两平面的夹角β\beta与所占方位角ϕ\phi相等这一关系, c′′c''处与bc′⌢\overset{\frown}{bc'}垂直的微分圆弧长度为sin⁡bc′⌢dbc′⌢dβ\sin\overset{\frown}{bc'}d\overset{\frown}{bc'}\mathrm{d}\beta. 此时可根据bc′⌢\overset{\frown}{bc'}上的积分求出微分三角形的面积.

Abc′′⌢=∫0bc′′⌢sin⁡θdθdβ=(1−cos⁡bc′′⌢)dβ \begin{equation} \begin{aligned} A_{\overset{\frown}{bc''}}=\int_0^{\overset{\frown}{bc''}} \sin\theta \mathrm{d}\theta \mathrm{d}\beta=(1-\cos\overset{\frown}{bc''})\mathrm{d}\beta \end{aligned} \end{equation}

根据Abc′′⌢A_{\overset{\frown}{bc''}}与Abc′⌢A_{\overset{\frown}{bc'}}的比值, 可得根据ϵ1\epsilon_1采样c′′c''的方式, 向量c′′\mathbf{c''}即为采样方向ω\omega. 由于不再需要判断光线是否相交, 此时根据o+td=(1−b1−b2)v0+b1v1+b2v2o+t\mathbf{d}=(1-b_1-b_2)v_0+b_1v_1+b_2v_2求解重心坐标即可.

ϵ1=1−cos⁡bc′′⌢1−cos⁡bc′⌢cos⁡bc′′⌢=1−ϵ1(1−cos⁡bc′⌢) \begin{equation} \begin{aligned} \epsilon_1&=\frac{1-\cos\overset{\frown}{bc''}}{1-\cos\overset{\frown}{bc'}}\\ \cos\overset{\frown}{bc''}&=1-\epsilon_1(1-\cos\overset{\frown}{bc'}) \end{aligned} \end{equation}

完整采样过程实现在SampleSphericalTriangle.

// Sampling Function Definitions
PBRT_CPU_GPU pstd::array<Float, 3> SampleSphericalTriangle(const pstd::array<Point3f, 3> &v, Point3f p,
                                              Point2f u, Float *pdf) {
    if (pdf)
        *pdf = 0;
    // Compute vectors _a_, _b_, and _c_ to spherical triangle vertices
    Vector3f a(v[0] - p), b(v[1] - p), c(v[2] - p);
    CHECK_GT(LengthSquared(a), 0);
    CHECK_GT(LengthSquared(b), 0);
    CHECK_GT(LengthSquared(c), 0);
    a = Normalize(a);
    b = Normalize(b);
    c = Normalize(c);

    // Compute normalized cross products of all direction pairs
    Vector3f n_ab = Cross(a, b), n_bc = Cross(b, c), n_ca = Cross(c, a);
    if (LengthSquared(n_ab) == 0 || LengthSquared(n_bc) == 0 || LengthSquared(n_ca) == 0)
        return {};
    n_ab = Normalize(n_ab);
    n_bc = Normalize(n_bc);
    n_ca = Normalize(n_ca);

    // Find angles $\alpha$, $\beta$, and $\gamma$ at spherical triangle vertices
    Float alpha = AngleBetween(n_ab, -n_ca);
    Float beta = AngleBetween(n_bc, -n_ab);
    Float gamma = AngleBetween(n_ca, -n_bc);

    // Uniformly sample triangle area $A$ to compute $A'$
    Float A_pi = alpha + beta + gamma;
    Float Ap_pi = Lerp(u[0], Pi, A_pi);
    if (pdf) {
        Float A = A_pi - Pi;
        *pdf = (A <= 0) ? 0 : 1 / A;
    }

    // Find $\cos\beta'$ for point along _b_ for sampled area
    Float cosAlpha = std::cos(alpha), sinAlpha = std::sin(alpha);
    Float sinPhi = std::sin(Ap_pi) * cosAlpha - std::cos(Ap_pi) * sinAlpha;
    Float cosPhi = std::cos(Ap_pi) * cosAlpha + std::sin(Ap_pi) * sinAlpha;
    Float k1 = cosPhi + cosAlpha;
    Float k2 = sinPhi - sinAlpha * Dot(a, b) /* cos c */;
    Float cosBp = (k2 + (DifferenceOfProducts(k2, cosPhi, k1, sinPhi)) * cosAlpha) /
                  ((SumOfProducts(k2, sinPhi, k1, cosPhi)) * sinAlpha);
    // Happens if the triangle basically covers the entire hemisphere.
    // We currently depend on calling code to detect this case, which
    // is sort of ugly/unfortunate.
    CHECK(!IsNaN(cosBp));
    cosBp = Clamp(cosBp, -1, 1);

    // Sample $c'$ along the arc between $a$ and $c$
    Float sinBp = SafeSqrt(1 - Sqr(cosBp));
    Vector3f cp = cosBp * a + sinBp * Normalize(GramSchmidt(c, a));

    // Compute sampled spherical triangle direction and return barycentrics
    Float cosTheta = 1 - u[1] * (1 - Dot(cp, b));
    Float sinTheta = SafeSqrt(1 - Sqr(cosTheta));
    Vector3f w = cosTheta * b + sinTheta * Normalize(GramSchmidt(cp, b));
    // Find barycentric coordinates for sampled direction _w_
    Vector3f e1 = v[1] - v[0], e2 = v[2] - v[0];
    Vector3f s1 = Cross(w, e2);
    Float divisor = Dot(s1, e1);
    CHECK_RARE(1e-6, divisor == 0);
    if (divisor == 0) {
        // This happens with triangles that cover (nearly) the whole
        // hemisphere.
        return {1.f / 3.f, 1.f / 3.f, 1.f / 3.f};
    }
    Float invDivisor = 1 / divisor;
    Vector3f s = p - v[0];
    Float b1 = Dot(s, s1) * invDivisor;
    Float b2 = Dot(w, Cross(s, e1)) * invDivisor;

    // Return clamped barycentrics for sampled direction
    b1 = Clamp(b1, 0, 1);
    b2 = Clamp(b2, 0, 1);
    if (b1 + b2 > 1) {
        b1 /= b1 + b2;
        b2 /= b1 + b2;
    }
    return {Float(1 - b1 - b2), Float(b1), Float(b2)};
}

对于需要根据采样方向倒推PDF的情况, 由于c′′c''已知, 根据bc′′⌢\overset{\frown}{bc''}与ac⌢\overset{\frown}{ac}的交点可求得c′c'. 交点为圆弧构成的平面的相交直线与球的交点, 直线方向通过两平面法线叉乘求得, 选择与a+c\mathbf{a}+\mathbf{c}同向的交点, 同时已知交线会经过球心, 此时可以确定c′c'以及abc′abc'的面积, 以此求得ϵ0\epsilon_0. ϵ1\epsilon_1根据概率公式直接求出即可.

pbrt实现了InvertSphericalTriangleSample执行倒推过程, 对于ac′ac'长度较小的情况, 为避免浮点误差β′\beta'较小时认为ϵ0=0\epsilon_0=0.

// Via Jim Arvo's SphTri.C
PBRT_CPU_GPU Point2f InvertSphericalTriangleSample(const pstd::array<Point3f, 3> &v, Point3f p,
                                      Vector3f w) {
    // Compute vectors _a_, _b_, and _c_ to spherical triangle vertices
    Vector3f a(v[0] - p), b(v[1] - p), c(v[2] - p);
    CHECK_GT(LengthSquared(a), 0);
    CHECK_GT(LengthSquared(b), 0);
    CHECK_GT(LengthSquared(c), 0);
    a = Normalize(a);
    b = Normalize(b);
    c = Normalize(c);

    // Compute normalized cross products of all direction pairs
    Vector3f n_ab = Cross(a, b), n_bc = Cross(b, c), n_ca = Cross(c, a);
    if (LengthSquared(n_ab) == 0 || LengthSquared(n_bc) == 0 || LengthSquared(n_ca) == 0)
        return {};
    n_ab = Normalize(n_ab);
    n_bc = Normalize(n_bc);
    n_ca = Normalize(n_ca);

    // Find angles $\alpha$, $\beta$, and $\gamma$ at spherical triangle vertices
    Float alpha = AngleBetween(n_ab, -n_ca);
    Float beta = AngleBetween(n_bc, -n_ab);
    Float gamma = AngleBetween(n_ca, -n_bc);

    // Find vertex $\VEC{c'}$ along $\VEC{a}\VEC{c}$ arc for $\w{}$
    Vector3f cp = Normalize(Cross(Cross(b, w), Cross(c, a)));
    if (Dot(cp, a + c) < 0)
        cp = -cp;

    // Invert uniform area sampling to find _u0_
    Float u0;
    if (Dot(a, cp) > 0.99999847691f /* 0.1 degrees */)
        u0 = 0;
    else {
        // Compute area $A'$ of subtriangle
        Vector3f n_cpb = Cross(cp, b), n_acp = Cross(a, cp);
        CHECK_RARE(1e-5, LengthSquared(n_cpb) == 0 || LengthSquared(n_acp) == 0);
        if (LengthSquared(n_cpb) == 0 || LengthSquared(n_acp) == 0)
            return Point2f(0.5, 0.5);
        n_cpb = Normalize(n_cpb);
        n_acp = Normalize(n_acp);
        Float Ap = alpha + AngleBetween(n_ab, n_cpb) + AngleBetween(n_acp, -n_cpb) - Pi;

        // Compute sample _u0_ that gives the area $A'$
        Float A = alpha + beta + gamma - Pi;
        u0 = Ap / A;
    }

    // Invert arc sampling to find _u1_ and return result
    Float u1 = (1 - Dot(w, b)) / (1 - Dot(cp, b));
    return Point2f(Clamp(u0, 0, 1), Clamp(u1, 0, 1));
}

Kajiya方程中的cos⁡θ\cos\theta项同样也会影响方差, pbrt将其包括在pdf中. 由于在ϵ0ϵ1\epsilon_0\epsilon_1上它是平滑变化的, 已知aa对应(0,1)(0, 1), bb对应(0,0)(0, 0)与(1,0)(1, 0), cc对应(1,1)(1, 1), 通过在单位正方形将顶点对应的cos⁡θ\cos\theta作为采样权重, 内部采样权重为顶点双线性插值结果, 可执行重要性抽样, 将双线性采样拆解为两次线性采样即可.

// Sample spherical triangle from reference point
// Apply warp product sampling for cosine factor at reference point
Float pdf = 1;
if (ctx.ns != Normal3f(0, 0, 0)) {
    // Compute $\cos\theta$-based weights _w_ at sample domain corners
    Point3f rp = ctx.p();
    Vector3f wi[3] = {Normalize(p0 - rp), Normalize(p1 - rp), Normalize(p2 - rp)};
    pstd::array<Float, 4> w =
        pstd::array<Float, 4>{std::max<Float>(0.01, AbsDot(ctx.ns, wi[1])),
                              std::max<Float>(0.01, AbsDot(ctx.ns, wi[1])),
                              std::max<Float>(0.01, AbsDot(ctx.ns, wi[0])),
                              std::max<Float>(0.01, AbsDot(ctx.ns, wi[2]))};

    u = SampleBilinear(u, w);
    DCHECK(u[0] >= 0 && u[0] < 1 && u[1] >= 0 && u[1] < 1);
    pdf = BilinearPDF(u, w);
}

PBRT_CPU_GPU inline Point2f SampleBilinear(Point2f u, pstd::span<const Float> w) {
    DCHECK_EQ(4, w.size());
    Point2f p;
    // Sample $y$ for bilinear marginal distribution
    p.y = SampleLinear(u[1], w[0] + w[1], w[2] + w[3]);

    // Sample $x$ for bilinear conditional distribution
    p.x = SampleLinear(u[0], Lerp(p.y, w[0], w[2]), Lerp(p.y, w[1], w[3]));

    return p;
}

双线性片
#

由四个点组成的面可以覆盖[0,1]2[0,1]^2的uv空间的形状即位双线性片, 表面上的点对应的值可以通过插值获取. 与三角形网格类似, BilinearPatchMesh内部数据存储在cache中, 构造时变换到渲染空间, BilinearPatch中只存储序号.

f(u,v)=(1−u)(1−v)p0,0+u(1−v)p1,0+(1−u)vp0,1+uvp1,1∂p∂u=(1−v)(p1,0−p0,0)+v(p1,1−p0,1)∂p∂v=(1−u)(p1,0−p0,0)+u(p1,1−p0,1) \begin{equation} \begin{aligned} f(u,v) &= (1 - u)(1 - v)p_{0,0} + u (1 - v) p_{1,0} + (1 - u) v p_{0,1} + u v p_{1,1}\\ \frac{\partial p}{\partial u} &= (1 - v)(p_{1,0} - p_{0,0}) + v(p_{1,1} - p_{0,1})\\ \frac{\partial p}{\partial v} &= (1 - u)(p_{1,0} - p_{0,0}) + u(p_{1,1} - p_{0,1}) \end{aligned} \end{equation}

双线性片很多时候以矩形的形式出现, 这种情况可以简化很多计算, pbrt通过判断某个点到第四个点的向量是否与前三个点组成的平面的法线垂直来检查四个点是否共面, 通过判断四个点到他们的中心距离是否相等来判断是否是矩形.

PBRT_CPU_GPU
bool IsRectangle(const BilinearPatchMesh *mesh) const {
    // Get bilinear patch vertices in _p00_, _p01_, _p10_, and _p11_
    const int *v = &mesh->vertexIndices[4 * blpIndex];
    Point3f p00 = mesh->p[v[0]], p10 = mesh->p[v[1]];
    Point3f p01 = mesh->p[v[2]], p11 = mesh->p[v[3]];

    if (p00 == p01 || p01 == p11 || p11 == p10 || p10 == p00)
        return false;
    // Check if bilinear patch vertices are coplanar
    Normal3f n(Normalize(Cross(p10 - p00, p01 - p00)));
    if (AbsDot(Normalize(p11 - p00), n) > 1e-5f)
        return false;

    // Check if planar vertices form a rectangle
    Point3f pCenter = (p00 + p01 + p10 + p11) / 4;
    Float d2[4] = {DistanceSquared(p00, pCenter), DistanceSquared(p01, pCenter),
                    DistanceSquared(p10, pCenter), DistanceSquared(p11, pCenter)};
    for (int i = 1; i < 4; ++i)
        if (std::abs(d2[i] - d2[0]) / d2[0] > 1e-4f)
            return false;
    return true;
}

双线性片的法线包围锥体由四个点处法线的均值决定.

// Compute average normal and return normal bounds for patch
Vector3f n = Normalize(n00 + n10 + n01 + n11);
Float cosTheta = std::min({Dot(n, n00), Dot(n, n01), Dot(n, n10), Dot(n, n11)});
return DirectionCone(n, Clamp(cosTheta, -1, 1));

相交测试
#

首先找到一条对uu进行插值形成的直线使得它到光线的距离为0, 计算两条不平行的直线的距离首先求二者叉乘获取法线, 再与两条直线分别组成两个平行的平面, 求二者的距离即可. 最终得到的uu的二次方程系数如下.

a=((p1,0−p0,0)×(p0,1−p1,1))⋅dc=(p0,0−o)×d)⋅(p0,1−p0,0)b=(p1,0−o)×d)⋅(p1,1−p1,0)−(a+c) \begin{equation} \begin{aligned} a &= ((p_{1,0} - p_{0,0}) \times (p_{0,1} - p_{1, 1})) \cdot \mathbf{d}\\ c &= (p_{0,0} - \mathbf{o}) \times \mathbf{d}) \cdot (p_{0,1} - p_{0,0})\\ b &= (p_{1,0} - \mathbf{o}) \times \mathbf{d}) \cdot (p_{1,1} - p_{1,0}) - (a + c) \end{aligned} \end{equation}

若uu有解则可以按如下方式得到vv与光线的tt值.

fu(v)=(1−v)pu,0+vpu,1 \begin{equation} \begin{aligned} f_u(v) &= (1 - v)p_{u,0} + vp_{u,1}\\ \end{aligned} \end{equation} t=det(fu(0)−o,fu(1)−fu(0),d×(fu(1)−fu(0)))∥d×(fu(1)−fu(0))∥2v=det(fu(0)−o,d,d×(fu(1)−fu(0)))∥d×(fu(1)−fu(0))∥2 \begin{equation} \begin{aligned} t &= \frac{\text{det}(f_u(0) - \mathbf{o}, f_u(1) - f_u(0), \mathbf{d} \times (f_u(1) - f_u(0)))}{\Vert \mathbf{d} \times (f_u(1) - f_u(0)) \Vert^2}\\ v &= \frac{\text{det}(f_u(0) - \mathbf{o}, \mathbf{d}, \mathbf{d} \times (f_u(1) - f_u(0)))}{\Vert \mathbf{d} \times (f_u(1) - f_u(0)) \Vert^2} \end{aligned} \end{equation}

为了不与BilinearPatch本身的uvuv混淆, 贴图的uvuv用stst来表示. 由于SurfaceInteraction中的位置在uvuv上的偏导实际上是在sv上的偏导, pbrt通过链式法则求得该值.

∂p∂u=p1,0+v(p1,1−p1,0)−p0,0−v(p0,1−p0,0)∂s∂u=s1,0+v(s1,1−s1,0)−s0,0−v(s0,1−s0,0) \begin{equation} \begin{aligned} \frac{\partial p}{\partial u} &= p_{1,0} + v (p_{1,1} - p{1,0}) - p_{0, 0} - v (p_{0, 1} - p{0, 0})\\ \frac{\partial s}{\partial u} &= s_{1,0} + v (s_{1,1} - s{1,0}) - s_{0, 0} - v (s_{0, 1} - s{0, 0})\\ \end{aligned} \end{equation}

表面采样
#

pbrt中双线性片采样的pdf与当前的微分面积相关, 这使得uvuv变化较为剧烈处不会聚集过多的采样点.

// Compute PDF for sampling the $(u,v)$ coordinates given by _intr.uv_
Float pdf;
if (mesh->imageDistribution)
    pdf = mesh->imageDistribution->PDF(uv);
else if (!IsRectangle(mesh)) {
    // Initialize _w_ array with differential area at bilinear patch corners
    pstd::array<Float, 4> w = {
        Length(Cross(p10 - p00, p01 - p00)), Length(Cross(p10 - p00, p11 - p10)),
        Length(Cross(p01 - p00, p11 - p01)), Length(Cross(p11 - p10, p11 - p01))};

    pdf = BilinearPDF(uv, w);
} else
    pdf = 1;

// Find $\dpdu$ and $\dpdv$ at bilinear patch $(u,v)$
Point3f pu0 = Lerp(uv[1], p00, p01), pu1 = Lerp(uv[1], p10, p11);
Vector3f dpdu = pu1 - pu0;
Vector3f dpdv = Lerp(uv[0], p01, p11) - Lerp(uv[0], p00, p10);

// Return final bilinear patch area sampling PDF
return pdf / Length(Cross(dpdu, dpdv));

对于非矩形, 有采样分布贴图或面积较小的双线性片采用面积采样转化为立体角采样的方法, 否则将矩形投影到球上形成球面矩形再进行采样. 对于矩形的uvuv, 可以通过将点投影到边上并根据边长归一化获得.

曲线
#

pbrt的Curve类采用具有宽度的一维三次Bézier曲线, 可以渲染为平面, 圆柱体或带状曲面. 对于带状曲面, pbrt会提供顶点的法线用于后续的插值, 支持只表示曲线的某一部分.

p(u)=(1−u)3p0+3(1−u)2up1+3(1−u)u2p2+u3p3 \begin{equation} p(u) = (1 - u)^3 p_0 + 3(1 - u)^2 u p_1 + 3(1 - u)u^2 p_2 + u^3 p_3 \end{equation}

包围结构
#

pbrt首先计算控制点形成的包围盒, 再根据曲线宽度扩展包围盒. 法线包围结构返回的是完整的球这一保守结果.

PBRT_CPU_GPU Bounds3f Curve::Bounds() const {
    pstd::span<const Point3f> cpSpan(common->cpObj);
    Bounds3f objBounds = BoundCubicBezier(cpSpan, uMin, uMax);
    // Expand _objBounds_ by maximum curve width over $u$ range
    Float width[2] = {Lerp(uMin, common->width[0], common->width[1]),
                      Lerp(uMax, common->width[0], common->width[1])};
    objBounds = Expand(objBounds, std::max(width[0], width[1]) * 0.5f);

    return (*common->renderFromObject)(objBounds);
}

相交测试
#

Curve类的相交测试接口通过调用IntersectRay成员函数来实现, 通过不断二分判断当前曲线包围盒是否与光线相交, 最终使得曲线区域接近线性线段方便计算相交.

将曲线转换到光线空间使用lookAt矩阵实现, 相机正上方对应的向量会被设置为与曲线两端连成的向量垂直, 这使得曲线位于xyxy平面上且与xx轴接近平行, 可以获取范围更小的包围盒以加速二分.

pbrt根据variation diminishing判断细分曲线所需要的深度, 深度为0时控制点形成的曲线接近直线. 由于可能出现与曲线拥有多个交点的情况, 非阴影光线递归求交时若深度不为0二分的两个分支都会搜索.

// Compute refinement depth for curve, _maxDepth_
Float L0 = 0;
for (int i = 0; i < 2; ++i)
    L0 = std::max(
        L0, std::max(std::max(std::abs(cp[i].x - 2 * cp[i + 1].x + cp[i + 2].x),
                                std::abs(cp[i].y - 2 * cp[i + 1].y + cp[i + 2].y)),
                        std::abs(cp[i].z - 2 * cp[i + 1].z + cp[i + 2].z)));
int maxDepth = 0;
if (L0 > 0) {
    Float eps = std::max(common->width[0], common->width[1]) * .05f;  // width / 20
    // Compute log base 4 by dividing log2 in half.
    int r0 = Log2Int(1.41421356237f * 6.f * L0 / (8.f * eps)) / 2;
    maxDepth = Clamp(r0, 0, 10);
}

深度为0判断相交时, 首先根据控制点获得首尾两侧与曲线垂直的直线, 计算在xyxy平面上对于当前原点的边函数, 以判断光线交点是否在曲线范围内.

// Intersect ray with curve segment
// Test ray against segment endpoint boundaries
// Test sample point against tangent perpendicular at curve start
Float edge = (cp[1].y - cp[0].y) * -cp[0].y + cp[0].x * (cp[0].x - cp[1].x);
if (edge < 0)
    return false;

// Test sample point against tangent perpendicular at curve end
edge = (cp[2].y - cp[3].y) * -cp[3].y + cp[3].x * (cp[3].x - cp[2].x);
if (edge < 0)
    return false;

对于带状曲面, pbrt会根据球面插值得到的法线调整条带宽度.

// Compute $u$ coordinate of curve intersection point and _hitWidth_
Float u = Clamp(Lerp(w, u0, u1), u0, u1);
Float hitWidth = Lerp(u, common->width[0], common->width[1]);
Normal3f nHit;
if (common->type == CurveType::Ribbon) {
    // Scale _hitWidth_ based on ribbon orientation
    if (common->normalAngle == 0)
        nHit = common->n[0];
    else {
        Float sin0 =
            std::sin((1 - u) * common->normalAngle) * common->invSinNormalAngle;
        Float sin1 =
            std::sin(u * common->normalAngle) * common->invSinNormalAngle;
        nHit = sin0 * common->n[0] + sin1 * common->n[1];
    }
    hitWidth *= AbsDot(nHit, ray.d) / rayLength;
}

曲线上的uvuv空间的vv与点在与曲线uu处的垂线上的位置有关.

// Initialize _SurfaceInteraction_ _intr_ for curve intersection
// Compute $v$ coordinate of curve intersection point
Float ptCurveDist = std::sqrt(ptCurveDist2);
Float edgeFunc = dpcdw.x * -pc.y + pc.x * dpcdw.y;
Float v = (edgeFunc > 0) ? 0.5f + ptCurveDist / hitWidth
                         : 0.5f - ptCurveDist / hitWidth;

u上的偏微分可以直接由Bézier曲线的定义获取. 对于条状曲线, 由于法线已知, 可以直接计算v上的偏微分. 对于平面, v上的偏微分方向沿着垂线, 且长度与当前宽度相等. 对于圆柱体, v上的偏微分会绕u上的偏微分旋转.

// Compute $\dpdu$ and $\dpdv$ for curve intersection
Vector3f dpdu, dpdv;
EvaluateCubicBezier(pstd::MakeConstSpan(common->cpObj), u, &dpdu);
CHECK_NE(Vector3f(0, 0, 0), dpdu);
if (common->type == CurveType::Ribbon)
    dpdv = Normalize(Cross(nHit, dpdu)) * hitWidth;
else {
    // Compute curve $\dpdv$ for flat and cylinder curves
    Vector3f dpduPlane = objectFromRay.ApplyInverse(dpdu);
    Vector3f dpdvPlane =
        Normalize(Vector3f(-dpduPlane.y, dpduPlane.x, 0)) * hitWidth;
    if (common->type == CurveType::Cylinder) {
        // Rotate _dpdvPlane_ to give cylindrical appearance
        Float theta = Lerp(v, -90, 90);
        Transform rot = Rotate(-theta, dpduPlane);
        dpdvPlane = rot(dpdvPlane);
    }
    dpdv = objectFromRay(dpdvPlane);
}

浮点精度
#

浮点精度对渲染结果具有很大的影响, 典型的就是实时渲染中的depth bias, 光追中的浮点误差通常发生在求交中. 比较简单的做法是通过添加ϵ\epsilon值来做offset, pbrt的浮点误差处理不需要使用这种方式.

浮点算数
#

算数运算
#

默认当前浮点运算采用的是IEEE754的默认舍入模式, 即4以下向下取, 6以上想上去, 5舍入到最近的偶数, 这种方式可以减少四舍五入时5总是向上舍入带来的系统误差.

对于指数为ee的浮点数, 各个浮点数之间的最小距离为2e−232^{e-23}, 这个距离被称为最低有效位的量级(ulp). 对于浮点加法的运算结果, 令其误差范围为[(a+b)(1−ϵ),(a+b)(1+ϵ)][(a+b)(1-\epsilon),(a+b)(1+\epsilon)], 由于此时ulp为(a+b)2−23(a+b)2^{-23}, 可知∣ϵ∣≤2−24≈5.960464…×10−8|\epsilon| \le 2^{-24} \approx 5.960464\ldots \times 10^{-8}.

工具代码
#

pbrt支持获取浮点数中的下一个值, 即当前指数下距离最近的两个浮点数, 除NaN, 无穷大等特殊情况, 这通过直接修改最后一位bit来实现.

PBRT_CPU_GPU
inline float NextFloatUp(float v) {
    // Handle infinity and negative zero for _NextFloatUp()_
    if (IsInf(v) && v > 0.f)
        return v;
    if (v == -0.f)
        v = 0.f;

    // Advance _v_ to next higher float
    uint32_t ui = FloatToBits(v);
    if (v >= 0)
        ++ui;
    else
        --ui;
    return BitsToFloat(ui);
}

误差传播
#

绝对误差与相对误差的定义如下.

δa=∣a~−a∣δr=∣a~−aa∣=δaa \begin{equation} \begin{aligned} \delta_a &= |\tilde{a} - a|\\ \delta_r &= |\frac{\tilde{a} - a}{a}| = \frac{\delta_a}{a} \end{aligned} \end{equation}

舍入后的浮点数与对应的实数具有如下关系.

a~∈a±σa=a(1±σr) \begin{equation} \tilde{a} \in a \pm \sigma_a = a(1 \pm \sigma_r) \end{equation}

进行多次浮点加法的舍入结果如下.

(((a⊕b)⊕c)⊕d)∈(a+b)(1±ϵm)3+c(1±ϵm)2+d(1±ϵm) \begin{equation} (((a \oplus b) \oplus c) \oplus d) \in (a+b)(1 \pm \epsilon_m)^3 + c(1 \pm \epsilon_m)^2 + d(1 \pm \epsilon_m) \end{equation}

对于高阶误差的幂具有如下的不等式.

(1±ϵm)n≤(1±(n+1)ϵm) \begin{equation} (1 \pm \epsilon_m)^n \le (1 \pm (n + 1) \epsilon_m) \end{equation}

此时简化舍入后的结果, 获得的最大浮点误差如下.

4ϵm∣a+b∣+3ϵm∣c∣+2ϵm∣d∣ \begin{equation} 4 \epsilon_m |a + b| + 3 \epsilon_m |c| + 2 \epsilon_m |d| \end{equation}

上述方法被称为前向误差分析, 也可以通过判断输入得到同一个输出的输入范围来进行后向误差分析, 这不是很适用于几何计算.

下述不等式具有更精确的上界.

(1±ϵm)n≤1+γnγn=nϵm1−nϵm \begin{equation} \begin{aligned} (1 \pm \epsilon_m)^n \le 1 + \gamma_n\\ \gamma_n = \frac{n \epsilon_m}{1 - n \epsilon_m} \end{aligned} \end{equation}

可以用上式来表示计算时输入本身所带有的误差, 对于乘法和加法可得如下结果, 可以看出乘法的误差比加法小很多, 且当加法输入绝对值相近但符号相反时误差会很大, 这种现象被称为灾难性抵消.

a(1±γi)⊗b(1±γj)∈ab(1±γi+j+1)δr=γi+j+1 \begin{equation} \begin{aligned} a(1 \pm \gamma_i) \otimes b(1 \pm \gamma_j) &\in ab(1 \pm \gamma_{i + j + 1})\\ \delta_r &= \gamma_{i+j+1} \end{aligned} \end{equation} a(1±γi)⊕b(1±γj)∈a(1±γi+1)+b(1±γj+1)δr=∣a∣γi+1+∣b∣γj+1a+b \begin{equation} \begin{aligned} a(1 \pm \gamma_i) \oplus b(1 \pm \gamma_j) &\in a(1 \pm \gamma_{i+1}) + b(1 \pm \gamma_{j+1})\\ \delta_r &= \frac{|a|\gamma_{i+1} + |b|\gamma_{j+1}}{a + b} \end{aligned} \end{equation}

误差分析
#

pbrt通过Interval类提供误差分析的功能, 每次执行运算时它都会执行误差区间累积的计算.

保守光线-包围结构相交
#

光线相交计算的误差如下. 若tmin⁡t_{\min}与tmax⁡t_{\max}的误差区间重合, pbrt会给tmax⁡t_{\max}增加2γ3tmax⁡2 \gamma_3 t_{\max}来保守的确定光线与物体相交.

t=(x⊖ox)⊗(1⊘dx)∈x−oxdx(1±ϵ)3x−oxdx∈t(1±γ3)δr=γ3∣t∣ \begin{equation} \begin{aligned} t = (x \ominus o_x) &\otimes (1 \oslash \mathbf{d}_x) \in \frac{x - o_x}{\mathbf{d}_x}(1 \pm \epsilon)^3\\ \frac{x - o_x}{\mathbf{d}_x} &\in t (1 \pm \gamma_3)\\ \delta_r &= \gamma_3 |t| \end{aligned} \end{equation}

精确二次方程判别式
#

球和圆柱等物体的相交计算需要用到二次方程判别式b2−4acb^2 - 4ac, 当物体较远时若b2≈4acb^2 \approx 4ac会导致灾难性抵消现象的出现.

球体和圆柱体的二次方程判别式可以转为如下形式, 此时不再需要计算cc, 因此误差更小.

b2−4ac=4a(b24a−c)=4a((o⋅d^)2−(o⋅o)+r2)=4a(−∥o−b2ad∥2+r2) \begin{equation} \begin{aligned} b^2 - 4ac &= 4a\left(\frac{b^2}{4a} - c\right)\\ &= 4a\left((\mathbf{o} \cdot \hat{\mathbf{d}})^2 - (\mathbf{o} \cdot \mathbf{o}) + r^2\right)\\ &= 4a\left(-\Vert \mathbf{o} - \frac{b}{2a}\mathbf{d} \Vert^2 + r^2\right) \end{aligned} \end{equation} b2−4ac=4a(r2−∥oxy−b2adxy∥2) \begin{equation} b^2 - 4ac = 4a \left(r^2 - \Vert \mathbf{o_{xy}} - \frac{b}{2a} \mathbf{d_{xy}} \Vert^2\right) \end{equation}

稳定三角形相交
#

由于边方程只需要判断符号, 浮点误差的影响不大. 对于边方程为0的情况暂时没有很好的精度误差检查方法, 需要用double重新计算, 实际使用中几乎不可能有这种情况.

边界交点误差
#

若tt本身带有误差, 光线相交的绝对误差如下.

δa=γ1∣ox∣+δt(1+γ2)∣dx∣+γ2∣thitdx∣ \begin{equation} \delta_a = \gamma_1|o_x|+\delta_t(1 + \gamma_2)|\mathbf{d_x}|+\gamma_2|t_{\text{hit}}\mathbf{d_x}| \end{equation}

二次曲面重投影
#

计算出与二次曲面的交点后可以进行重投影来减小误差, 例如下式中对于球面可以通过半径来重投影, 误差为γ5\gamma_5. 圆柱只需要在xyxy上投影, 误差为γ3\gamma_3. 圆盘由于只需要设置zz值, 它相当于是没有误差的.

x′=xrx2+y2+z2 \begin{equation} x' = x \frac{r}{x^2 + y^2 + z^2} \end{equation}

三角形参数评估
#

采用边方程计算重心坐标后再进行插值的误差如下.

δa=γ7(∣b0x0∣+∣b1x1∣+∣b2x2∣) \begin{equation} \delta_a = \gamma_7(|b_0x_0| + |b_1x_1| + |b_2x_2|) \end{equation}

双线性片参数评估
#

双线性插值误差如下.

γ6(∣x0,0∣+∣x0,1∣+∣x1,0∣+∣x1,1∣) \begin{equation} \gamma_6(|x_{0,0}|+|x_{0,1}|+|x_{1,0}|+|x_{1,1}|) \end{equation}

曲线参数评估
#

为避免离开曲线的光线与其重新相交, 误差会被设置为曲线的宽度, 若曲线宽度过宽用双线性片来代替是更好的选择.

变换误差
#

变换后的x′x'误差如下, 下式分别是原值不包括和包括误差的情况.

γ3(∣m0,0x∣+∣m0,1y∣+∣m0,2∣z+∣m0,3∣) \begin{equation} \gamma_3(|m_{0,0}x|+|m_{0,1}y|+|m_{0,2}|z+|m_{0,3}|) \end{equation} (γ3+1)(∣m0,0∣δx+∣m0,1∣δy+∣m0,2∣δz)+γ3(∣m0,0x∣+∣m0,1y∣+∣m0,2z∣+∣m0,3∣) \begin{equation} (\gamma_3+1)(|m_{0,0}|\delta_x + |m_{0,1}|\delta_y + |m_{0,2}|\delta_z) + \gamma_3(|m_{0,0}x| + |m_{0,1}y| + |m_{0,2}z| + |m_{0,3}|) \end{equation}

稳定生成光线起点
#

我们需要将某个面附近的光线起点沿着法线移动, 使得它位于交点的误差范围之外, 以此使用最小的偏移量避免错误的相交, 保证了阴影, 反射等效果的质量. 偏移向量定义如下, 为避免向下舍入pbrt会将值移动一个ulp.

n⋅(δx,δy,δz) \begin{equation} \mathbf{n} \cdot (\delta_x, \delta_y, \delta_z) \end{equation} 
for (int i = 0; i < 3; ++i) {
    if (offset[i] > 0)      po[i] = NextFloatUp(po[i]);
    else if (offset[i] < 0) po[i] = NextFloatDown(po[i]);
}

在阴影与光源求交的计算中有可能因为距离光源过近的物体导致错误相交, pbrt通过提前停止来求交来简单的处理.

constexpr Float ShadowEpsilon = 0.0001f;

为避免变换带来的误差, 每次变换中光线起点都会移动到误差边界.

if (Float lengthSquared = LengthSquared(d); lengthSquared > 0) {
    Float dt = Dot(Abs(d), o.Error()) / lengthSquared;
    o += d * dt;
    if (tMax)
        *tMax -= dt;
}

避免光线起点后方的相交
#

由于误差实际为负的tt值可能计算结果为正, 对于部分形状pbrt实现了高效的保守误差分析.

三角形
#

采用边方程重心坐标计算出的tt误差如下, 若不超过这个值则认为交点在后方.

δe=2(γ2max⁡i∣xi∣max⁡i∣yi∣+δymax⁡i∣xi∣+δxmax⁡i∣yi∣)δt=3(γ3max⁡i∣ei∣∣zi∣+δemax⁡i∣zi∣+δzmax⁡i∣ei∣) \begin{equation} \begin{aligned} \delta_e &= 2(\gamma_2 \max_i|x_i|\max_i|y_i| + \delta_y\max_i|x_i|+\delta_x\max_i|y_i|)\\ \delta_t &= 3(\gamma_3 \max_i|e_i||z_i| + \delta_e\max_i|z_i|+\delta_z\max_i|e_i|) \end{aligned} \end{equation}

双线性片
#

pbrt通过计算一个简单的ϵ\epsilon来解决后方相交误差.

cmax⁡(v)=max⁡(vx,vy,vz)ϵ=γ10(cmax⁡(∣o∣)+cmax⁡(∣d∣)+cmax⁡(∣p0,0∣)+cmax⁡(∣p0,1∣)+cmax⁡(∣p1,0∣)+cmax⁡(∣p1,1∣)) \begin{equation} \begin{aligned} c_{\max}(\mathbf{v}) &= \max(\mathbf{v_x}, \mathbf{v_y}, \mathbf{v_z})\\ \epsilon &= \gamma_10 (c_{\max}(|\mathbf{o}|) + c_{\max}(|\mathbf{d}|)\\ &+ c_{\max}(|\mathbf{p_{0,0}}|)+ c_{\max}(|\mathbf{p_{0,1}}|)\\ &+ c_{\max}(|\mathbf{p_{1,0}}|)+ c_{\max}(|\mathbf{p_{1,1}}|)) \end{aligned} \end{equation}
← pbrt-v4 Ep. V: 相机模型 pbrt-v4 Ep. VII: 图元与加速结构 →