【WebGPU实时光追美少女】BVH与射线场景求交插值

少女dtysky

世界Skill

时刻2021.09.26

路标

本系列文章设计的所有代码均已开源，Github仓库在这里：dtysky/webgpu-renderer。

在上一篇文章管线组织与GBuffer中，我提到了路径追踪的第一步是射线和场景的求交，进一步引出了BVH的概念。在这一章，我将会论述如何利用BVH对场景中的所有三角形进行高效得分割和求交，并且给出得出交点后属性插值的方法。

不知道为何最近有点抑郁，所以更新效率低了一些...

求交部分的代码都在buildin/shaders/ray-tracing/hitTest.chunk.wgsl内。

射线和场景求交

首先要讨论的是射线和场景的求交，或者用在物理引擎中的叫法碰撞测试hitTest。我们并不需要先论述BVH这样的加速结构，因为无论哪种加速结构，到最后总是要进行射线和三角形的求交和插值。

射线和三角形描述

为了求交计算，需要先对射线和三角形进行描述。首先是射线，描述一条射线有很多种方式，这里使用最便于GPU计算的结构：

struct Ray {
  origin: vec3<f32>;
  dir: vec3<f32>;
  invDir: vec3<f32>;
};

通过一个原点origin和方向dir来使用方程p = origin + t * dir描述射线，invDir是方向向量的倒数，便于后续使用。有了射线的结构后，考虑如何生成射线，但这一部分我将放在下一章来讲。

接下来是三角形，在前面的文章我论述过如何将整个场景的数据进行合并，其中就包括了顶点数据indexData的合并。一个三角形由三个顶点描述，所以只需要知道三角形三个顶点的索引即可。这里先不讨论BVH，所以就认为我们是在顺序、每三个顶点遍历整个indexData即可。得到顶点索引后，便可以直接从storage buffer中获得对应的顶点属性，比如u_positions.value[index0]。

求交

有了射线和三角形的信息，便可以进行求交计算。射线和三角形的求交思路不止一种，一个最简单和直观的思路为：

先求射线和三角形所在平面的交点。
再判断交点是否在三角形内部。

但由于需要首先确定三角形所在平面，这么做实际运算效率其实不高，实际上我们可以直接通过方程来描述三角形，和射线方程联立整理，来简化计算：

空间射线方程：P = O + t * D
空间三角形方程：P = u * E1 + v * E2，其中E1 = V1 - V0和E2 = V2 - V0是两条边的向量，u、v是权重，V1、V2和V3是三个顶点。

联立后得到方程：O + t * D = u * E1 + v * E2，这是一个三元方程组，自然可以化简写成矩阵形式：

[-D, E1, E2] * [t, u, v] = T，其中T = O - V0，[-D, E1, E2]是3x3列主序矩阵，[t, u, v]是转置（赖得打公式了...

而经过克莱姆法则和混合积公式推导，记det = D x E2 · E1最后可写为：

t = (T x E1 * E2) / det
u = (D x E2 * T) / det
v = (T x E1 * E2) / det

如此一来，便可以求出这三个系数，并在求取的过程中顺便将不满足条件的情况判断出来：

// 判断是否相交，如果相交，返回所有的顶点属性
fn triangleHitTest(ray: Ray, leaf: BVHLeaf) -> FragmentInfo {
  var info: FragmentInfo;
  let indexes: vec3<u32> = leaf.indexes;
  // 索引出三个顶点
  info.p0 = u_positions.value[indexes[0]];
  info.p1 = u_positions.value[indexes[1]];
  info.p2 = u_positions.value[indexes[2]];

  let e1: vec3<f32> = info.p1 - info.p0;
  let e2: vec3<f32> = info.p2 - info.p0;
  let p: vec3<f32> = cross(ray.dir, e2);
  var det: f32 = dot(e1, p);
  var t: vec3<f32> = ray.origin - info.p0;

  // 保证`det`为正
  if (det < 0.) {
    t = -t;
    det = -det;
  }

  // `det`为0，三向量共面，算作不相交
  if (det < 0.0001) {
    return info;
  }

  let u: f32 = dot(t, p);

  // `u / det`必须在(0,1)
  if (u < 0. || u > det) {
    return info;
  }

  let q: vec3<f32> = cross(t, e1);
  let v: f32 = dot(ray.dir, q);

  // `v / det`和`(1 - v - u) / det`必须在(0,1)
  if (v < 0. || v + u - det > 0.) {
    return info;
  }

  let lt: f32 = dot(e2, q);

  // lt即为射线方程的`t`，由于是射线，其必须大于0
  if (lt < 0.) {
    return info;
  }

  let invDet: f32 = 1. / det;
  info.weights = vec3<f32>(0., u, v) * invDet;
  info.weights.x = 1. - info.weights.y - info.weights.z;

  info.hit = true;
  info.t = lt * invDet;
  info.hitPoint = ray.origin + ray.dir * info.t;
  info.uv0 = u_uvs.value[indexes.x];
  info.uv1 = u_uvs.value[indexes.y];
  info.uv2 = u_uvs.value[indexes.z];
  info.n0 = u_normals.value[indexes.x];
  info.n1 = u_normals.value[indexes.y];
  info.n2 = u_normals.value[indexes.z];
  info.meshIndex = u_meshMatIndexes.value[indexes.x].x;
  info.matIndex = u_meshMatIndexes.value[indexes.x].y;
  let metallicRoughnessFactorNormalScaleMaterialType: vec4<f32> = material.u_metallicRoughnessFactorNormalScaleMaterialTypes[info.matIndex];
  info.matType = u32(metallicRoughnessFactorNormalScaleMaterialType[3]);
}

如果这些测试都通过，则射线和三角形相交，而我们也顺便得到了u、v、det，通过它们，也就顺便获得了重心坐标。于此同时，我还将属于这个三角形的三个顶点属性都存了下来，以供下一个阶段做插值。

重心坐标插值

重心坐标插值这个概念在传统光栅化渲染器的栅格化一步也会提到，即当某个像素被一个平面三角形覆盖时，需要通过三角形的三个顶点属性算出当前像素实际的顶点属性。而对于路径追踪渲染器，由于同样是用三角形描述场景，所以也需要这么一个插值的过程，来计算出交点实际的顶点属性值。

所谓重心坐标，就是上面的info.weights。对于一个三角形内的点，它们都在(0,1)之间，记录了三角形的每个顶点到交点的距离的比例，从这个角度而言，它们实际上也是这个交点相对于三角形三个顶点的权重。

// 通过相交点的参数，求取最终插值后的结果并返回
fn fillHitPoint(frag: FragmentInfo, ray: Ray) -> HitPoint {
  var info: HitPoint;

  info.hit = true;
  info.hited = frag.t;
  info.meshIndex = frag.meshIndex;
  info.matIndex = frag.matIndex;
  let metallicRoughnessFactorNormalScaleMaterialType: vec4<f32> = material.u_metallicRoughnessFactorNormalScaleMaterialTypes[frag.matIndex];
  info.position = frag.p0 * frag.weights[0] + frag.p1 * frag.weights[1] + frag.p2 * frag.weights[2];
  let uv: vec2<f32> = frag.uv0 * frag.weights[0] + frag.uv1 * frag.weights[1] + frag.uv2 * frag.weights[2];
  let textureIds: vec4<i32> = material.u_matId2TexturesId[frag.matIndex];  
  let faceNormal: vec3<f32> = getFaceNormal(frag);
  info.sign = sign(dot(faceNormal, -ray.dir));
  info.normal = info.sign * getNormal(frag, uv, textureIds[1], metallicRoughnessFactorNormalScaleMaterialType[2]);
  let baseColor: vec4<f32> = getBaseColor(material.u_baseColorFactors[frag.matIndex], textureIds[0], uv);
  info.baseColor = baseColor.rgb;
  info.glass = baseColor.a;
  info.matType = frag.matType;
  info.isSpecGloss = isMatSpecGloss(frag.matType);
  info.isGlass = isMatGlass(frag.matType);

  if (info.isSpecGloss) {
    let specularGlossinessFactors: vec4<f32> = material.u_specularGlossinessFactors[frag.matIndex];
    info.spec = getSpecular(specularGlossinessFactors.xyz, textureIds[2], uv).rgb;
    info.gloss = getGlossiness(specularGlossinessFactors[3], textureIds[2], uv);
  } else {
    info.metal = getMetallic(metallicRoughnessFactorNormalScaleMaterialType[0], textureIds[2], uv);
    info.rough = getRoughness(metallicRoughnessFactorNormalScaleMaterialType[1], textureIds[2], uv);
  }

  info.pbrData = pbrPrepareData(info.isSpecGloss, info.baseColor, info.metal, info.rough, info.spec, info.gloss);

  return info;
}

以上是插值、求得最后结果的逻辑，可以见到像是position、uv之类的顶点属性都使用了重心坐标进行插值，而后也是用插值后的结果来采样贴图、计算材质属性的。

加速结构

如果不考虑性能，理论上上面的算法已经足够满足需求，我们只需要在每次发射射线的时候用该射线和场景内所有三角形求交，得到所有相交的点，然后找到其中离射线原点最近的即可。然而这是不现实的，在当前的这个时代，一个场景有上万上十万个三角形都挺正常，在实时的情况下即便是1SPP开销都是巨大的，所以这就需要我们有办法来加速整个求交的过程，也就引出了接下来要论述的——加速结构。

常见加速结构

加速结构的本质在于场景划分，通过合适的算法划分场景，来达到索引优化，最终减少真正的计算。场景划分不仅仅用于路径追踪，实际上在传统的实时渲染中也有很大的作用，比如ForwardAdd中的多光源实现等等。

如图，常见的加速结构有八叉树、KD树、BSP树等等。但他们都属于空间划分，而空间划分对于我们的场景有很大的劣势：

物体可能存在与多个划分区域，划分不彻底。
为了建立正确的划分，需要求取区域和三角形的相交。

所以在实时光追的引擎中，大家一般都使用BVH，即层次包围盒：

如图，这种结构不针对空间，而是根据物体的层次来对场景做划分。这样可以保证每个三角形都存在于唯一的划分区域内，同时整棵树比较平衡、最大深度低，而且构建十分简单。但当然世上没有免费的优化，其代价为：

由于不是空间划分，所以包围盒之间可能重叠，必须完全遍历整棵树才能得到最终结果，某些情况下相对求交效率可能低一些。
虽然构建效率高，但结构变了需要完全重建。不过这一部分已经被现代算法优化得很好了，有GPU算法（但本文不会提及）。

不过相对来讲，BVH仍然是当下最适合实时路径追踪的加速结构。

构造BVH

构造BVH可以在CPU进行，也可以在GPU进行（使用莫顿码），但由于个人精力有限，所以这里只实现了传统CPU的版本。

CPU部分构造BVH的代码全部在extension/BVH.ts内。

这里我不会讲所有的构建的代码都事无巨细列出来，只会说一些思路，具体的可以参考实际代码。

简单来讲，为了构建一个BVH，我们至少需要三步：

通过当前合并好的三角形，求取每个三角形的包围盒（世界空间）。
通过这些包围盒，构建二叉树，树的叶子节点为实际的三角形信息，其他为包围盒信息。
将树状结构打平为紧密的数组来存储。

所以我也讲这三步实现为了三个方法。

求取包围盒

首先需要一个数据结构来描述包围盒，其声明为：

export class Bounds {
  // 最大边界
  public max: Float32Array;
  // 最小边界
  public min: Float32Array;
  // 中点
  get center(): Float32Array;
  // 尺寸
  get size(): number;
  // 长度最大的轴
  get maxExtends(): EAxis;
  // 包围盒的表面积
  get surfaceArea(): number ;

  // 生成一个空的包围盒
  public initEmpty(): Bounds;
  // 从顶点生成一个包围盒
  public fromVertexes(v1: Float32Array, v2: Float32Array, v3: Float32Array): Bounds;
  // 扩张一个包围盒
  public update(v: Float32Array): Bounds;
  // 合并另一个包围盒，求并
  public mergeBounds(bounds: Bounds): Bounds;
  // 一个点是否在包围盒内
  public pointIn(p: Float32Array): boolean;
  // 获取一个点与包围盒某个轴最小值的偏移比例，一般0~1
  public getOffset(axis: EAxis, v: Float32Array): number;
}

这个类也是接下来所有操作的基础。有了它，便可以很简单地求得所有三角形的包围盒了：

protected _setupBoundsInfo = (worldPositions: Float32Array, indexes: Uint32Array) => {
  this._boundsInfos = [];

  for (let i = 0; i < indexes.length; i += 3) {
    const idxes = indexes.slice(i, i + 3);

    // 从`worldPositions`中解出每个顶点需要的数据
    copyTypedArray(3, tmpV1, 0, worldPositions, idxes[0] * 4);
    copyTypedArray(3, tmpV2, 0, worldPositions, idxes[1] * 4);
    copyTypedArray(3, tmpV3, 0, worldPositions, idxes[2] * 4);

    // 构建包围盒
    const bounds = new Bounds().fromVertexes(tmpV1, tmpV2, tmpV3);
    // 用于构建树的信息，包括包围盒本身，以及对应的三角形顶点索引
    this._boundsInfos.push({indexes: idxes, bounds});
  }
}

构建二叉树

有了包围盒后，便可以利用它们构建二叉树，构建BVH树的基本原理非常简单，首先我们要定义两个数据结构：结点BVHNode和叶子BVHLeaf：

export interface IBVHNode {
  // 沿哪个轴划分
  axis: EAxis;
  // 所有子节点合并的包围盒
  bounds: Bounds;
  // 左右子节点
  child0: IBVHNode | IBVHLeaf;
  child1: IBVHNode | IBVHLeaf;
  // 结点在树中的深度（层数）
  depth: number;
}

// 用于存储数个按策略不可划分的三角形
export interface IBVHLeaf {
  // 包含三角形的起始索引
  infoStart: number;
  // 包含三角形的结束索引
  infoEnd: number;
  // 所有子节点的包围盒
  bounds: Bounds;
}

最基础的BVH的构建本质上是一个递归，自顶而下不断二分，直到无法继续分割。同时在这个过程中将所有子节点的包围盒合并到父节点，最终无法分割的节点记为叶子节点BVHLeaf:

protected _buildRecursive(start: number, end: number, depth: number): IBVHNode | IBVHLeaf {
  const {_boundsInfos} = this;
  const bounds = new Bounds().initEmpty();
  // 合并所有子节点的包围盒
  for (let i = start; i < end; i += 1) {
    bounds.mergeBounds(_boundsInfos[i].bounds);
  }

  const nPrimitives = end - start;
  if (nPrimitives <= this._maxPrimitivesPerLeaf) {
    // 如果当前三角形数量不可再分割，直接返回叶子节点
    return { infoStart: start, infoEnd: end, bounds };
  } else {
    // 否则，通过一定规则算出新的分割点`mid`，进行下一轮递归
    // 这里还可能执行一些图元顺序的调整，来达到最有分割，分割有不同的策略

    // 递归计算左右子节点
    const child0 = this._buildRecursive(start, mid, depth + 1);
    const child1 = this._buildRecursive(mid, end, depth + 1);

    // 返回新的`BVHNode`
    return {
      bounds: new Bounds().initEmpty().mergeBounds(child0.bounds).mergeBounds(child1.bounds),
      axis: dim, child0, child1, depth
    };
  }
}

可见原理确实很简单，实际上对于这种传统的BVH构建算法，最精髓的在于注释中写道的分割策略。考察一个BVH的构建是否优秀，一般需要评估两点：

在相同的图元数量下，BVH的层数是否最低。
每一层两个子节点中的包围盒重叠程度是否足够低。

围绕这两点有各种不同的算法来权衡，在本项目的实现中，使用的是SAH - Surface Area Heuristic，即表面积启发式算法。这是一种从概率角度的优化，从理论上来讲，一条射线和场景求交的代价，主要由和图元即三角形的代价来决定，一次求交中计算的图元数量越多，也就意味着效率越低。而和图元求交的概率，又可以用包含图元的包围盒的重叠程度来表征，进一步可以用包围盒的表面积来表征，最大程度减少包围盒的表面积，就是SAH的目的，所以它实际上是解决了上面的第二个问题。

如图，实现上来讲，我们可以先求得所有包围盒中心所构成的最大包围盒，以此包围盒最长的轴作为用于切分的轴。如果包围盒已然不可划分（min=max），则直接返回，接着如果目前图元数量小于一定级别，则使用简单的包围盒中心位置做快速选择nth_element划分；再否则就进行正常的SAH，这里使用了一种基于均分的buckets（桶）的算法）——将图元按照中心位置划分到不同的桶内，再从低到高不断合并这些桶来求得划分代价，接着得出最低的划分代价，最终按照这个最低代价的桶为中点进行划分。实现详见代码。

打平为数组

有了BVH树还不够，因为最后一步要送到CS内进行最后的求交，所以必须要将其组装成一个storage buffer，这就要求我们设计一种紧凑的数据结构，来将整棵树存储到一个ArrayBuffer中。而我设计的结构如下:

对于Node，使用uint32 x 1来存储子节点的类型（最低位，0为Node，1为Leaf）和偏移，float x 3存储子节点的包围盒，共两个子节点，消耗8个float的长度。
对于Leaf，使用uint32 x 1来存储接下来有几个三角形属于同一个节点，uint32 x 3存储当前三角形的三个顶点索引。

有了结构，整个Buffer的构建也就比较简单，同样是递归构建：

protected _traverseNode = (
  node: IBVHNode | IBVHLeaf,
  info: {maxDepth: number, nodes: number[], leaves: number[]},
  depth: number = 1,
  parentOffset: number = -1,
  childIndex: number = 0
) => {
  info.maxDepth = Math.max(depth, info.maxDepth);
  const {nodes, leaves} = info;
  const {_boundsInfos, _bvhNodes, _bvhLeaves} = this;

  if (isBVHLeaf(node)) {
    _bvhLeaves.push(node);
    if (parentOffset >= 0) {
      nodes[parentOffset * 8 + childIndex] = (1 << 31) | (leaves.length / 4);
    }

    const count = node.infoEnd - node.infoStart;
    for (let i = node.infoStart; i < node.infoEnd; i += 1) {
      const idxes = _boundsInfos[i].indexes;

      leaves.push(count, idxes[0], idxes[1], idxes[2]);
    }
  } else {
    _bvhNodes.push(node);
    const bounds = node.bounds;
    const nodeOffset = nodes.length / 8;

    if (parentOffset >= 0) {
      nodes[parentOffset * 8 + childIndex] = nodeOffset * 2;
    }

    nodes.push(
      0, bounds.min[0], bounds.min[1], bounds.min[2],
      0, bounds.max[0], bounds.max[1], bounds.max[2]
    );

    this._traverseNode(node.child0, info, depth + 1, nodeOffset, 0);
    this._traverseNode(node.child1, info, depth + 1, nodeOffset, 4);
  }
}

注意这里我们构建好的Buffer正好是16字节对齐的。

BVH求交

构建了BVH后，使用u_bvh作为storage buffer传给CS来进行求交。BVH的求交远离比较简单巧妙，由于BVH本质上是一种AABB，同时和射线也并不需要得到真正的交点，只需要得到是否相交的结果。所以可以将其视为三对互相平行的平板，如图：

可以将射线投影到三个平面，最终需要满足最大的xmin - origin小于最小的xmax - origin，这个原理也好理解，如果将前者理解为光线进入盒子的时间，后者理解为光线离开盒子的时间，那么我们最终要保证进入时间小于离开时间：

fn boxHitTest(ray: Ray, max: vec3<f32>, min: vec3<f32>) -> f32 {
  let t1: vec3<f32> = (min - ray.origin) * ray.invDir;
  let t2: vec3<f32> = (max - ray.origin) * ray.invDir;
  let tvmin: vec3<f32> = min(t1, t2);
  let tvmax: vec3<f32> = max(t1, t2);
  let tmin: f32 = max(tvmin.x, max(tvmin.y, tvmin.z));
  let tmax: f32 = min(tvmax.x, min(tvmax.y, tvmax.z));

  if (tmax - tmin < 0.) {
    return -1.;
  }

  if (tmin > 0.) {
    return tmin;
  }

  return tmax;
}

以上方法返回的结果如果大于0，则表示射线和BVH相交，可见比起三角形求交，是十分高效的。而其中的max和min则是通过对数据的decode得来的：

struct Child {
  isLeaf: bool;
  offset: u32;
};

fn decodeChild(index: u32) -> Child {
  return Child((index >> 31u) != 0u, (index << 1u) >> 1u);
}

fn getBVHNodeInfo(offset: u32) -> BVHNode {
  var node: BVHNode;
  let data0: vec4<f32> = u_bvh.value[offset];
  let data1: vec4<f32> = u_bvh.value[offset + 1u];
  node.child0Index = bitcast<u32>(data0[0]);
  node.child1Index = bitcast<u32>(data1[0]);
  node.min = data0.yzw;
  node.max = data1.yzw;

  return node;
}

fn getBVHLeafInfo(offset: u32) -> BVHLeaf {
  var leaf: BVHLeaf;
  let data1: vec4<f32> = u_bvh.value[offset];
  leaf.primitives = bitcast<u32>(data1.x);
  leaf.indexes.x = bitcast<u32>(data1.y);
  leaf.indexes.y = bitcast<u32>(data1.z);
  leaf.indexes.z = bitcast<u32>(data1.w);

  return leaf;
}

综合实现

有了BVH、BVH求交算法、三角形求交算法、插值算法，便可以组装出最终的求交方案。但在这之前，我们需要先明确求交的分类：

普通射线求交：最普通的求交，用于间接光照，需要找到离射线原点最近的三角形。
阴影射线求交：阴影射线，用于直接光照，详细在后面文章会讲到，只需要判断射线和一定范围内的三角形有相交。
光源求交：和光源求交，是和三角形不同的另一种算法，在这里顺便也讲了。

一般来讲，普通求交要和光源求交联合起来，因为如果射线最先相交的是光源，还进行正常的表面计算会出现问题，阴影射线求交理论同理（若只有一个光源可以忽略）。

普通射线

由于要寻找的是离射线原点最近的交点，所以一定是遍历完整棵树才能停下，这可以用一次dfs实现，途中还可以直接剪枝。dfs最简单的是使用递归，但显然我们不能再wgsl里使用递归，所以只能自己维护一个递归栈，用迭代来实现了：

fn hitTest(ray: Ray) -> HitPoint {
  var hit: HitPoint;
  var fragInfo: FragmentInfo;
  // 初始化为最大长度
  fragInfo.t = MAX_RAY_LENGTH;
  var node: BVHNode;
  // 构建BVH深度大小的结点栈，栈中数据为结点索引
  var nodeStack: array<u32, ${BVH_DEPTH}>;
  nodeStack[0] = 0u;
  var stackDepth: i32 = 0;

  loop {
    // 当前深度小于10，栈中已无结点，结束
     if (stackDepth < 0) {
       break;
     }

    let child = decodeChild(nodeStack[stackDepth]);
    stackDepth = stackDepth - 1;

    if (child.isLeaf) {
      // 进行叶子的求交，内部本质上是进行了三角形求交，最终返回求交的信息
      let info: FragmentInfo = leafHitTest(ray, child.offset, false);

      if (info.hit && info.t < fragInfo.t) {
        // 如果交点比当前交点近，则使用替换
        fragInfo = info;
      }

      continue;
    }

    node = getBVHNodeInfo(child.offset);
    // 进行结点的求交
    let hited: f32 = boxHitTest(ray, node.max, node.min);

    if (hited < 0.) {
      continue;
    }

    // 相交，两个子节点入栈，进入下一次迭代
    stackDepth = stackDepth + 1;
    nodeStack[stackDepth] = node.child0Index;
    stackDepth = stackDepth + 1;
    nodeStack[stackDepth] = node.child1Index;
  }

  if (fragInfo.hit) {
    // 如果有交点，进行重心坐标插值，同时生成材质属性
    hit = fillHitPoint(fragInfo, ray);
  }

  return hit;
}

阴影射线

阴影射线求交和普通求交大差不差，唯一的区别在于需要给定一个初始的最大值maxT，即射线原点到光源之上采样点之间的距离，并且判断有相交后即可结束求交，立即返回：

fn hitTestShadow(ray: Ray, maxT: f32) -> FragmentInfo {
  var fragInfo: FragmentInfo;
  var node: BVHNode;
  var nodeStack: array<u32, ${BVH_DEPTH}>;
  nodeStack[0] = 0u;
  var stackDepth: i32 = 0;

  loop {
     if (stackDepth < 0) {
       break;
     }

    let child = decodeChild(nodeStack[stackDepth]);
    stackDepth = stackDepth - 1;

    if (child.isLeaf) {
      let info: FragmentInfo = leafHitTest(ray, child.offset, true);

      if (info.hit && info.t < maxT && info.t > EPS) {
        // 判断在射线原点和光源采样点之间有遮挡，返回
        return info;
      }

      continue;
    }

    node = getBVHNodeInfo(child.offset);
    let hited: f32 = boxHitTest(ray, node.max, node.min);

    if (hited < 0.) {
      continue;
    }

    stackDepth = stackDepth + 1;
    nodeStack[stackDepth] = node.child0Index;
    stackDepth = stackDepth + 1;
    nodeStack[stackDepth] = node.child1Index;
  }

  return fragInfo;
}

光源

最后就是和光源的求交了，早在前面的文章WebGPU基础与简单渲染器中，我们就提到了本引擎是如何描述光源，尤其是面光源的：

struct LightInfo {
  lightType: u32;
  areaMode: u32;
  areaSize: vec2<f32>;
  color: vec4<f32>;
  worldTransform: mat4x4<f32>;
  worldTransformInverse: mat4x4<f32>;
};

而对于面光源（无论是什么形状），一个通用的做法是先求取射线和平面的交点，然后判断交点是否在光源内，这个也就是本文一开始提到的、低效的三维平面求交。但这里好在引擎已经提供了光源的worldTransformInverse，所以理论上我们可以认为面光源总是在自己的本地空间（XZ平面），此时只要反向将射线变换到光源的本地空间，问题就变得简单了很多：

求取射线与平面的交点被简化了，最终得到的是一个XZ平面的点。
判定交点和光源的关系被简化到了二维空间。

目前引擎支持disc和rect两种面光源，在平面上可以非常简单得判定包含关系：

fn hitTestXZPlane(ray: Ray, inverseMat: mat4x4<f32>) -> vec3<f32> {
  // 射线方向变换到光源本地空间
  let invDir: vec3<f32> = normalize((inverseMat * vec4<f32>(ray.dir, 0.)).xyz);
  // 光源本地空间法线的反向基向量
  let normal: vec3<f32> = vec3<f32>(0., 0., 1.);
  let dot: f32 = dot(invDir, normal);

  if (abs(dot) < EPS) {
    // 射线平行或背向光源
    return vec3<f32>(MAX_RAY_LENGTH, MAX_RAY_LENGTH, MAX_RAY_LENGTH);
  }

  // 射线原点变换到光源本地空间
  let invOrigin: vec3<f32> = (inverseMat * vec4<f32>(ray.origin, 1.)).xyz;
  // 求取相交段的射线步长
  let t: f32 = dot(-invOrigin, normal) / dot;

  if (t < EPS) {
    return vec3<f32>(MAX_RAY_LENGTH, MAX_RAY_LENGTH, MAX_RAY_LENGTH);
  }

  return vec3<f32>(invOrigin.xy + t * invDir.xy, t);
}

fn hitTestLights(ray: Ray) -> vec4<f32> {
  let areaLight: LightInfo = global.u_lightInfos[0];
  // xyz存光的颜色，w存射线相交段的`t`
  var res: vec4<f32> = vec4<f32>(areaLight.color.rgb, MAX_RAY_LENGTH);
  // 目前只支持一个光源
  if (areaLight.lightType != LIGHT_TYPE_AREA) {
    return res;
  }

  // 求得XZ平面上的交点，以及相交段射线的`t`
  let xyt: vec3<f32> = hitTestXZPlane(ray, areaLight.worldTransformInverse);

  var hit: bool = false;
  if (areaLight.areaMode == LIGHT_AREA_DISC) {
    // 判定是否在圆盘内
    hit = length(xyt.xy) < areaLight.areaSize.x;
  } else {
    // 判定是否在矩形内
    hit = all(abs(xyt.xy) < areaLight.areaSize / 2.);
  }

  if (hit) {
    res.a = xyt.z;
  }

  return res;
}

展示

以上就是本文的所有内容，但在最后，为了调试，我还做了一个功能来将将整个BVH的包围盒画出来，其实原理也很简单：

通过包围盒信息，生成line类型的图元数据（主要是索引数据是一个图元两个点，而非三角形的三个点）。
组装Mesh绘制。

由于太过简单这里就不详细讲了，直接看看最后的结果吧：

点击查看评论