Whitted光线追踪实现

Whitted光线追踪是的光线追踪中最经典示例之一。该技术由特纳·惠特（Turner Whitted）于1980年在Siggraph论文“用于阴影显示的改进的照明模型”中首次详细介绍。在Whitted撰写论文之前，大多数程序已经能够模拟漫反射和镜面反射的外观。著名的Phong照明模型已经众所周知。但是，模拟复杂的反射和折射尚未完成。 Whitted建议使用光线跟踪来解决此问题。如果对象A是类似镜子的表面，并且对象B位于其顶部，则我们希望看到B反射到A的反射。如果A不是平面，则没有简单的解决方案可以计算出反射。如果B也是反射性的，事情就变得更加艰难。镜面不断反射自身的图像，从而引起电影中经常看到的“无穷大空间”效果。如下图效果：

我们举一个例子，下图显示了一个示例：

其中球体在图像最终到达眼睛之前被两个反射镜反射了两次。那么透明对象呢？我们应该看到透明的表面，但是由于折射现象，诸如水或玻璃之类的材料会使光线弯曲。反射方向仅取决于表面方向和入射光方向。可以使用斯涅尔定律计算折射方向，该方向取决于表面方向（表面法线），入射光方向和材料折射率（水大约为1.3，玻璃大约为1.5）。

Whitted建议使用这些定律来计算光线与反射或透明表面相交时的反射和折射方向，并遵循这些光线的路径来找出它们将相交的对象的颜色。这些反射和折射光线的交点可能会发生三种情况：

1：如果相交点处的表面不透明且具有漫反射特性，我们要做的就是使用光照模型（例如Phong模型）来计算相交点处对象的颜色。此过程还涉及朝场景中的每个光源方向投射光线，以查找该点是否在阴影中。这些射线称为阴影射线。
2：如果表面是类似镜子的表面，我们只需在相交点处跟踪另一条反射射线即可。
3：如果曲面是透明曲面，则在交点处投射一条反射射线和一条折射射线。

下图就是包含反射和折射的例子(其中S为阴影射线，R为反射射线，T为透射射线)：

如何定义表面的颜色呢？我们观察下图：

如果P1产生与P2相交的二次射线，并且在P2处产生另一个与P3的另一个对象相交的二次射线，则P3的颜色变为P2的颜色，P2的颜色又变为P1的颜色。透明物体的折射率越高，这些镜面反射越强。另外，入射角越大，反射的光越多（这是由于菲涅耳效应引起的）。对于折射射线，使用相同的方法。如果光线被折射，则P1折射，然后在P2折射，然后最终在P3入射到散射对象，P2呈现P3的颜色，P1呈现P2的颜色，如图下图所示：

我们现在来看一下Whitted光线跟踪的算法流程：

Whitted的算法使用后向跟踪来模拟反射和折射。我们从眼睛开始，投射主光线。
找到第一个可见表面，并判断表面的种类：
如果物体表面不透明并属于漫反射特性，则在相交点计算该表面的颜色。
如果物体表面不透明但属于镜面反射特性，则仅投射反射射线，并执行2.
如果物体透明，则投射一条透射射线和一条反射射线，并分别对其执行2.
如果投射光线不与任何物体相交，则返回背景颜色。

Whitted光线跟踪的主要消耗在2.中的射线求交，以及判断表面种类后的射线投射递归。避免无限递归，比图全是镜子的场景，我们通常会设置一个递归深度。

Whitted Ray-Tracing的实现

Whitted光线追踪算法的代码主要分为以下4块：

从摄像机的每个像素发射一条主射线。(创建主射线)
根据发射的主射线和场景对象进行交互象。(投射射线交互)
根据相交的对象，以及交点获得表面属性(法线，纹理坐标，以及表面特性(如类似镜面只进行完美反射)等)。(获取表面属性)
根据获得的表面特性，判断是否需要进一步发射射线。如玻璃，表面特性即有完美反射特性，又有折射特性，则需要递归发射两条新射线，进行跟踪计算。(递归光线跟踪)

创建主射线

该部分实现如下图所示功能：

我们从摄像机点到当前像素发射一条射线，该射线获得的最终颜色会保存在当前像素。值得注意的是，若我们物体表面特性既不是完美镜面反射特性，也不是折射特性时，我们定义其为默认特性(这里使用phong模型着色)。当表面为默认特性时，我们判断该点到光源位置是否有障碍。如果有，则表明该点处于阴影中，如下所示：

从物体表面的点到光源的射线，我们称为阴影射线，该部分会在代码块4.中详细介绍。当我们从每个像素发射射线，并求得最后像素的颜色，我们就能得到一张渲染完成的图片，下图演示了这一过程。

我们在投射每一条射线时，需要清楚其起点是相机的起点，方向是从相机的起点到像素中心的向量。要计算点在像素中心的位置，我们需要转换像素坐标从栅格空间到世界空间。该过程可以理解为我们光栅化的逆过程：

我们首先需要吧像素归一化，即将像素空间转换为NDC空间(归一化设备坐标):

因为我们希望最终的摄影机射线穿过像素的中间。所以我们要向像素位置添加一个小的偏移(0.5)，这会使NDC空间约束在[0,1]之间，虽然这和我们在OpenGL中常见的[-1,1]有所不同，但是我们后续获得的屏幕空间是[-1,1]即可。我们需要位于图像左侧的像素应具有负的x坐标，而位于右侧的像素应具有正的x坐标：

但是y值要注意，对于位于x轴上方的像素，y为负，而对于位于下方的像素，y为正（反之亦然）。下面的公式将纠正此问题：

但是，我们的宽高比往往都不是1:1，现在，让我们看一下图像的尺寸为7 x 5像素的情况：

用宽度除以图像的高度得到1.4。在屏幕空间中定义像素坐标时，它们的范围为[-1,1]。但是，这些像素新会沿x轴将像素压缩。为了使它们（像素）再次变为正方形，我们需要将像素的x坐标乘以图像纵横比，在这种情况下为1.4。请注意，此操作将y像素坐标（在屏幕空间中）保持不变。它们仍然在[-1,1]范围内，但x像素坐标现在在[-1.4,1.4]范围内，我们变换后的屏幕空间坐标公式：

最后，我们需要考虑视口范围。到目前为止，屏幕空间中定义的任何点的y坐标都在[-1，1]范围内。我们通常设定图像平面与相机原点相距1个单位。如果从侧面看摄像机，则可以通过将摄像机的原点连接到胶片平面的顶部和底部边缘来绘制三角形：

我们发现改变 $\alpha$ 角度会影响BC的长度。就可以获得：

那么我们最后的屏幕空间x，y坐标就可以通过以下公式计算：

我们就可以写出第一部分的代码：

//options：设定，objects：场景内所有对象，lights场景内所以灯光
void render(const Options &options,const std::vector<std::unique_ptr<Object>> &objects,const std::vector<std::unique_ptr<Light>> &lights)
{//创建图像Vec3f *framebuffer = new Vec3f[options.width * options.height];Vec3f *pix = framebuffer;//为每个像素投射射线部分float scale = tan(deg2rad(options.fov * 0.5));float imageAspectRatio = options.width / (float)options.height;Vec3f orig(0);for (uint32_t j = 0; j < options.height; ++j) {for (uint32_t i = 0; i < options.width; ++i) {// 生成主要射线float x = (2 * (i + 0.5) / (float)options.width - 1) * imageAspectRatio * scale;float y = (1 - 2 * (j + 0.5) / (float)options.height) * scale;Vec3f dir = normalize(Vec3f(x, y, -1));//投射射线*(pix++) = castRay(orig, dir, objects, lights, options, 0);}}// 保存图像为ppmstd::ofstream ofs;ofs.open("./out.ppm");ofs << "P6\n" << options.width << " " << options.height << "\n255\n";for (uint32_t i = 0; i < options.height * options.width; ++i) {char r = (char)(255 * clamp(0, 1, framebuffer[i].x));char g = (char)(255 * clamp(0, 1, framebuffer[i].y));char b = (char)(255 * clamp(0, 1, framebuffer[i].z));ofs << r << g << b;}ofs.close();delete [] framebuffer;
}

我们这里直接使摄像机在世界坐标(0,0)朝向vec3(0,0,-1)构建屏幕空间，这样就省去了屏幕空间向世界空间的转换。

投射射线交互

castRay函数是Whitted光追的主体函数，根据上述对Whitted光追算法对讲解，这里我们的代码为：

//orig：投射光线起点，dir：投射光线方向，objects：场景所有对象，lights：所有灯光，options：设定，depth：最大递归深度
Vec3f castRay(const Vec3f &orig, const Vec3f &dir,const std::vector<std::unique_ptr<Object>> &objects,const std::vector<std::unique_ptr<Light>> &lights,const Options &options,uint32_t depth,bool test = false)
{//超过最大递归深度，直接返回背景颜色if (depth > options.maxDepth) {return options.backgroundColor;}//光线击中点的颜色Vec3f hitColor = options.backgroundColor;//击中对象的表面属性float tnear = kInfinity;Vec2f uv;//重心坐标uint32_t index = 0;Object *hitObject = nullptr;if (trace(orig, dir, objects, tnear, index, uv, &hitObject)) {//光线几何求交Vec3f hitPoint = orig + dir * tnear;Vec3f N; // 法线Vec2f st; // 纹理坐标hitObject->getSurfaceProperties(hitPoint, dir, index, uv, N, st);//获取击中表面属性Vec3f tmp = hitPoint;switch (hitObject->materialType) {case REFLECTION_AND_REFRACTION:{//...}case REFLECTION:{//...}default:{//...}}}return hitColor;
}

这里的代码表示了：

castRay在递归深度到达阈值时直接返回背景颜色。若未到达则利用函数trace()和场景里的所有对象进行几何求交计算，并使用求交计算获得的对象hitObject获得表面属性(法线N，纹理坐标st)。然后根据击中对象材质的不同，进行不同的操作。

获取表面属性

获取表面属性主要由函数trace()和getSurfaceProperties()完成。trace进行射线求交计算。getSurfaceProperties为求交计算后返回的对象成员函数，用来获得表面属性。

我们看trace()函数：

//orig：射线起点，dir：射线方向，objects：场景所有对象，tNear：最近交点的距离，index：三角形于三角网格对象的索引，uv：三角形的重心坐标，hitObject：击中对象
bool trace(const Vec3f &orig, const Vec3f &dir,const std::vector<std::unique_ptr<Object>> &objects,float &tNear, uint32_t &index, Vec2f &uv, Object **hitObject)
{*hitObject = nullptr;for (uint32_t k = 0; k < objects.size(); ++k) {float tNearK = kInfinity;uint32_t indexK;Vec2f uvK;if (objects[k]->intersect(orig, dir, tNearK, indexK, uvK) && tNearK < tNear) {*hitObject = objects[k].get();tNear = tNearK;index = indexK;uv = uvK;}}return (*hitObject != nullptr);
}

该函数就对场景所有对象遍历，判断是否相交，如果相交且距离更近，则更新信息。对象求交在我https://blog.csdn.net/qq_39300235/article/details/105232626这篇文章有详细介绍，这里就直接贴上代码：

//对象虚基类
class Object
{public:Object() :materialType(DIFFUSE_AND_GLOSSY),ior(1.3), Kd(0.8), Ks(0.2), diffuseColor(0.2), specularExponent(25) {}virtual ~Object() {}virtual bool intersect(const Vec3f &, const Vec3f &, float &, uint32_t &, Vec2f &) const = 0;virtual void getSurfaceProperties(const Vec3f &, const Vec3f &, const uint32_t &, const Vec2f &, Vec3f &, Vec2f &) const = 0;virtual Vec3f evalDiffuseColor(const Vec2f &) const { return diffuseColor; }// material propertiesMaterialType materialType;float ior;float Kd, Ks;Vec3f diffuseColor;float specularExponent;
};//圆的解析几何求接
bool solveQuadratic(const float &a, const float &b, const float &c, float &x0, float &x1)
{float discr = b * b - 4 * a * c;if (discr < 0) return false;else if (discr == 0) x0 = x1 = - 0.5 * b / a;else {float q = (b > 0) ?-0.5 * (b + sqrt(discr)) :-0.5 * (b - sqrt(discr));x0 = q / a;x1 = c / q;}if (x0 > x1) std::swap(x0, x1);return true;
}//圆
class Sphere : public Object{
public:Sphere(const Vec3f &c, const float &r) : center(c), radius(r), radius2(r * r) {}bool intersect(const Vec3f &orig, const Vec3f &dir, float &tnear, uint32_t &index, Vec2f &uv) const{// analytic solutionVec3f L = orig - center;float a = dotProduct(dir, dir);float b = 2 * dotProduct(dir, L);float c = dotProduct(L, L) - radius2;float t0, t1;if (!solveQuadratic(a, b, c, t0, t1)) return false;if (t0 < 0) t0 = t1;if (t0 < 0) return false;tnear = t0;return true;}void getSurfaceProperties(const Vec3f &P, const Vec3f &I, const uint32_t &index, const Vec2f &uv, Vec3f &N, Vec2f &st) const{ N = normalize(P - center); }Vec3f center;float radius, radius2;
};bool rayTriangleIntersect(const Vec3f &v0, const Vec3f &v1, const Vec3f &v2,const Vec3f &orig, const Vec3f &dir,float &tnear, float &u, float &v)
{Vec3f edge1 = v1 - v0;Vec3f edge2 = v2 - v0;Vec3f pvec = crossProduct(dir, edge2);//P=D叉乘E2float det = dotProduct(edge1, pvec);//P点乘E1if (det == 0 || det < 0) return false;//不存在的情况(平行或背面)Vec3f tvec = orig - v0;u = dotProduct(tvec, pvec);if (u < 0 || u > det) return false;//u最后的值只能在[0,1]Vec3f qvec = crossProduct(tvec, edge1);//Q=T叉乘E1v = dotProduct(dir, qvec);if (v < 0 || u + v > det) return false;//u+v最后的值只能在[0,1]float invDet = 1 / det;tnear = dotProduct(edge2, qvec) * invDet;u *= invDet;v *= invDet;return true;
}//三角网格
class MeshTriangle : public Object
{
public:MeshTriangle(const Vec3f *verts,const uint32_t *vertsIndex,const uint32_t &numTris,const Vec2f *st){uint32_t maxIndex = 0;for (uint32_t i = 0; i < numTris * 3; ++i)if (vertsIndex[i] > maxIndex) maxIndex = vertsIndex[i];maxIndex += 1;vertices = std::unique_ptr<Vec3f[]>(new Vec3f[maxIndex]);memcpy(vertices.get(), verts, sizeof(Vec3f) * maxIndex);vertexIndex = std::unique_ptr<uint32_t[]>(new uint32_t[numTris * 3]);memcpy(vertexIndex.get(), vertsIndex, sizeof(uint32_t) * numTris * 3);numTriangles = numTris;stCoordinates = std::unique_ptr<Vec2f[]>(new Vec2f[maxIndex]);memcpy(stCoordinates.get(), st, sizeof(Vec2f) * maxIndex);}bool intersect(const Vec3f &orig, const Vec3f &dir, float &tnear, uint32_t &index, Vec2f &uv) const{bool intersect = false;for (uint32_t k = 0; k < numTriangles; ++k) {const Vec3f & v0 = vertices[vertexIndex[k * 3]];const Vec3f & v1 = vertices[vertexIndex[k * 3 + 1]];const Vec3f & v2 = vertices[vertexIndex[k * 3 + 2]];float t, u, v;if (rayTriangleIntersect(v0, v1, v2, orig, dir, t, u, v) && t < tnear) {tnear = t;uv.x = u;uv.y = v;index = k;intersect |= true;}}return intersect;}void getSurfaceProperties(const Vec3f &P, const Vec3f &I, const uint32_t &index, const Vec2f &uv, Vec3f &N, Vec2f &st) const{const Vec3f &v0 = vertices[vertexIndex[index * 3]];const Vec3f &v1 = vertices[vertexIndex[index * 3 + 1]];const Vec3f &v2 = vertices[vertexIndex[index * 3 + 2]];Vec3f e0 = normalize(v1 - v0);Vec3f e1 = normalize(v2 - v1);N = normalize(crossProduct(e0, e1));const Vec2f &st0 = stCoordinates[vertexIndex[index * 3]];const Vec2f &st1 = stCoordinates[vertexIndex[index * 3 + 1]];const Vec2f &st2 = stCoordinates[vertexIndex[index * 3 + 2]];st = st0 * (1 - uv.x - uv.y) + st1 * uv.x + st2 * uv.y;}Vec3f evalDiffuseColor(const Vec2f &st) const{float scale = 5;float pattern = (fmodf(st.x * scale, 1) > 0.5) ^ (fmodf(st.y * scale, 1) > 0.5);return mix(Vec3f(0.815, 0.235, 0.031), Vec3f(0.937, 0.937, 0.231), pattern);}std::unique_ptr<Vec3f[]> vertices;uint32_t numTriangles;std::unique_ptr<uint32_t[]> vertexIndex;std::unique_ptr<Vec2f[]> stCoordinates;
};

上面我们也看到getSurfaceProperties函数获取法线和纹理坐标。最后就是我们的光线递归投射了。

递归光线跟踪

在讲递归光线跟踪之前，我们先提一下反射，折射和菲涅耳现象。

反射

我们知道反射现象是入射角等于反射角，我们如何通过入射角计算出反射角呢，根据下图：

可以得到：

向量B就是I在N上的投影(这里我们的向量都是单位向量)，方向与N相反：

这里的 $cos(\theta )$ 可以计算为I点乘N(小于0)。要注意的是这里的 $\theta$ 是I和N的夹角(大于90度)，把B带入公式可以计算A和R：

反射代码就出来了：

Vec3f reflect(const Vec3f &I, const Vec3f &N)
{ return I - 2 * dotProduct(I, N) * N;
}

折射

当光线经过两种不同的介质时会发生折射现象，其光线方向会改变。折射图示：

向量T为入射光线I经过折射生成的新向量。根据经过的介质不同，我们折射角度也不同，影响折射的介质属性称为折射率，它们的关系由Snell定律给出：

图示中 $\eta _{1}$ 为空气的折射率， $\eta _{2}$ 为水的折射率。如何去计算T向量呢？我们看下图：

我们设定：

其中A和B为：

我们知道I和T是属于同一平面的，我们的M就是在这一平面内且与N垂直的单位向量：

其中C是I在N上的投影，方向与N一致：

这里的 $\theta _{1}$ 是-I和N的夹角( $\theta _{1}$ 小于90度)。将上述公式合并可得：

根据Snell折射定律：

带入T，最后得出：

我们定义：

化简得：

但是我们可能射线会从内部射出：

因此我们的法线在处理这种情况需要反转，且交换折射率。这里我们要注意的是，这里要确保 $cos(\theta _{1})$ 始终为正，所以我们代码为：

//ior：材料的折射率
Vec3f refract(const Vec3f &I, const Vec3f &N, const float &ior)
{float cosi = clamp(-1, 1, dotProduct(I, N));float etai = 1, etat = ior;Vec3f n = N;//如果cosi<0,表示从外部射入，否则表示从内部射出if (cosi < 0) { cosi = -cosi; } else { std::swap(etai, etat); n= -N; }float eta = etai / etat;float k = 1 - eta * eta * (1 - cosi * cosi);//当k小于零则表示全反射return k < 0 ? 0 : eta * I + (eta * cosi - sqrtf(k)) * n;
}

最后我们还需要考虑全反射的情况，如果k<0，则表示全反射，无折射射线。

菲涅耳

透明物体(例如玻璃或水)既折射又反射。它们反射的光量与透射的光量实际上取决于入射角。当入射角减小时，透射的光量增加。并且由于按照能量守恒定律，反射光量加上折射光量必须等于入射光总量(不考虑光损失)，因此可以推断出，入射角增大时，反射光量会增加，角度接近90度时最高可达100％。从技术上讲，玻璃球的边缘是100％反射的。但是在球体的中心，该球体仅反射约6％的入射光。我们可以看一幅图：

可以看到，离我们越近的水面，投射现象越强烈，离我们越远的水面反射现象越强烈。反射光与折射光的数量可以使用菲涅耳方程来计算。光由两个垂直波组成，我们称其为平行偏振光和垂直偏振光。我们需要使用两个不同的方程式计算这两个波的反射光比率，并对结果求平均值。两个菲涅耳方程为：

通过求平均值我们可以得到反射的比率：

项η1，η2是两种介质的折射率。项cosθ1和cosθ2分别是入射角和折射角。如前所述，由于能量守恒，折射光的比率可以简单地计算为：

我们可以写出代码：

void fresnel(const Vec3f &I, const Vec3f &N, const float &ior, float &kr)
{float cosi = clamp(-1, 1, dotProduct(I, N));float etai = 1, etat = ior;if (cosi > 0) {  std::swap(etai, etat); }// Snell's定律计算sintfloat sint = etai / etat * sqrtf(std::max(0.f, 1 - cosi * cosi));// 判断是否全反射if (sint >= 1) {kr = 1;//全反射}else {float cost = sqrtf(std::max(0.f, 1 - sint * sint));cosi = fabsf(cosi);float Rs = ((etat * cosi) - (etai * cost)) / ((etat * cosi) + (etai * cost));float Rp = ((etai * cosi) - (etat * cost)) / ((etai * cosi) + (etat * cost));kr = (Rs * Rs + Rp * Rp) / 2;}}

菲涅耳也需要考虑全反射的情况，如果sinθ2大于1，那么我们就有全反射的情况。在此特定情况下，无需计算菲涅耳公式。我们可以将FR设置为1。而且如果发现入射光线在折射率偏大的对象内部，则需要交换折射率。可以通过测试表面法线和入射射线方向之间的角度的余弦符号(cosθ1的符号)来再次完成此操作。

现在可以来实行我们的Whitted光线跟踪的递归部分了，我们先看处理表面是既有反射也有折射(玻璃)的情况：

case REFLECTION_AND_REFRACTION:{Vec3f reflectionDirection = normalize(reflect(dir, N));Vec3f refractionDirection = normalize(refract(dir, N, hitObject->ior));Vec3f reflectionRayOrig = (dotProduct(reflectionDirection, N) < 0) ?hitPoint - N * options.bias :hitPoint + N * options.bias;Vec3f refractionRayOrig = (dotProduct(refractionDirection, N) < 0) ?hitPoint - N * options.bias :hitPoint + N * options.bias;Vec3f reflectionColor = castRay(reflectionRayOrig, reflectionDirection, objects, lights, options, depth + 1, 1);Vec3f refractionColor = castRay(refractionRayOrig, refractionDirection, objects, lights, options, depth + 1, 1);float kr;fresnel(dir, N, hitObject->ior, kr);hitColor = reflectionColor * kr + refractionColor * (1 - kr);break;}

我们计算折射和反射向量，和它们的起点。然后递归使用castRay函数从新去跟踪折射射线和反射射线。获取的颜色值，利用菲涅耳计算出kr为反射颜色和折射颜色分配权重。这里要注意的是我们获取向量起点的方式做了一些操作。因为计算机精度的问题，直接使用起点产生噪声值，如下图：

这是因为表面的有些点由于损失精度，到了表面内部，我们可以用一个小小的偏移解决这个问题：

通过沿法线方向偏移之后可以避免这种情况，但是我们的射线在考虑到反射或者折射时，偏移方向会有所不同：

我们发现在内侧需要向内偏移，在外侧需要向外偏移，我们便可以用cos(θ)的符号来判断内外侧。即出射光线和法线点乘，小于零为内侧，大于零为外侧。

同样我们的表面只有反射属性的情况代码为：

case REFLECTION:{Vec3f reflectionDirection = normalize(reflect(dir, N));Vec3f reflectionRayOrig = (dotProduct(reflectionDirection, N) < 0) ?hitPoint - N * options.bias :hitPoint + N * options.bias;hitColor = castRay(reflectionRayOrig, reflectionDirection, objects, lights, options, depth + 1, 1);break;}

接下来就是处理默认材质情况。默认材质具有漫反射和镜面反射属性，我们使用phong光照模型来处理。直接上代码：

default:{Vec3f diffuse = 0, specularColor = 0;Vec3f shadowPointOrig = (dotProduct(dir, N) < 0) ?hitPoint + N * options.bias :hitPoint - N * options.bias;for (uint32_t i = 0; i < lights.size(); ++i) {Vec3f lightDir = lights[i]->position - hitPoint;float lightDistance2 = dotProduct(lightDir, lightDir);lightDir = normalize(lightDir);float LdotN = std::max(0.f, dotProduct(lightDir, N));Object *shadowHitObject = nullptr;float tNearShadow = kInfinity;//是否处于阴影bool inShadow = trace(shadowPointOrig, lightDir, objects, tNearShadow, index, uv, &shadowHitObject) &&tNearShadow * tNearShadow < lightDistance2;//漫反射部分diffuse += (1 - inShadow) * lights[i]->intensity * LdotN* hitObject->evalDiffuseColor(st);Vec3f reflectionDirection = reflect(-lightDir, N);//镜面反射部分specularColor += (1 - inShadow) * powf(std::max(0.f, -dotProduct(reflectionDirection, dir)), hitObject->specularExponent) * lights[i]->intensity;}hitColor = diffuse * hitObject->Kd + specularColor * hitObject->Ks;break;}

现在我们的完整castRay代码为：

Vec3f castRay(const Vec3f &orig, const Vec3f &dir,const std::vector<std::unique_ptr<Object>> &objects,const std::vector<std::unique_ptr<Light>> &lights,const Options &options,uint32_t depth,bool test = false)
{if (depth > options.maxDepth) {return options.backgroundColor;}Vec3f hitColor = options.backgroundColor;float tnear = kInfinity;Vec2f uv;uint32_t index = 0;Object *hitObject = nullptr;if (trace(orig, dir, objects, tnear, index, uv, &hitObject)) {Vec3f hitPoint = orig + dir * tnear;Vec3f N; // 法线Vec2f st; // 纹理坐标hitObject->getSurfaceProperties(hitPoint, dir, index, uv, N, st);Vec3f tmp = hitPoint;switch (hitObject->materialType) {case REFLECTION_AND_REFRACTION:{Vec3f reflectionDirection = normalize(reflect(dir, N));Vec3f refractionDirection = normalize(refract(dir, N, hitObject->ior));Vec3f reflectionRayOrig = (dotProduct(reflectionDirection, N) < 0) ?hitPoint - N * options.bias :hitPoint + N * options.bias;Vec3f refractionRayOrig = (dotProduct(refractionDirection, N) < 0) ?hitPoint - N * options.bias :hitPoint + N * options.bias;Vec3f reflectionColor = castRay(reflectionRayOrig, reflectionDirection, objects, lights, options, depth + 1, 1);Vec3f refractionColor = castRay(refractionRayOrig, refractionDirection, objects, lights, options, depth + 1, 1);float kr;fresnel(dir, N, hitObject->ior, kr);hitColor = reflectionColor * kr + refractionColor * (1 - kr);break;}case REFLECTION:{Vec3f reflectionDirection = normalize(reflect(dir, N));Vec3f reflectionRayOrig = (dotProduct(reflectionDirection, N) < 0) ?hitPoint - N * options.bias :hitPoint + N * options.bias;hitColor = castRay(reflectionRayOrig, reflectionDirection, objects, lights, options, depth + 1, 1);break;}default:{Vec3f diffuse = 0, specularColor = 0;Vec3f shadowPointOrig = (dotProduct(dir, N) < 0) ?hitPoint + N * options.bias :hitPoint - N * options.bias;for (uint32_t i = 0; i < lights.size(); ++i) {Vec3f lightDir = lights[i]->position - hitPoint;float lightDistance2 = dotProduct(lightDir, lightDir);lightDir = normalize(lightDir);float LdotN = std::max(0.f, dotProduct(lightDir, N));Object *shadowHitObject = nullptr;float tNearShadow = kInfinity;//是否处于阴影bool inShadow = trace(shadowPointOrig, lightDir, objects, tNearShadow, index, uv, &shadowHitObject) &&tNearShadow * tNearShadow < lightDistance2;//漫反射部分diffuse += (1 - inShadow) * lights[i]->intensity * LdotN* hitObject->evalDiffuseColor(st);Vec3f reflectionDirection = reflect(-lightDir, N);//镜面反射部分specularColor += (1 - inShadow) * powf(std::max(0.f, -dotProduct(reflectionDirection, dir)), hitObject->specularExponent) * lights[i]->intensity;}hitColor = diffuse * hitObject->Kd + specularColor * hitObject->Ks;break;}}}return hitColor;
}

我们whitted光追的实现效果图：

其中有两个光源，三个默认材质球和一个玻璃球。