PBRT中的误差舍入管理

由于计算机二进制的特性，计算机不可能精确表示所有实数。计算机用浮点数来代替实数。浮点数有固定的存储要求，但是使用时始终会引入误差。在进行多次计算后，误差便会一直积累。在光线求交应用中，误差积累是很致命的。可能会导致完全遗漏有效的相交点，例如，即使精确的相交值为正，但计算出的相交点的值为负。此外，计算出的射线形状相交点可能在形状的实际表面之上或之下。如下图所示：

可以发现当从阴影射线和反射射线计算出的交点开始跟踪新射线时，如果射线原点在实际表面以下，则可能会发现与该表面的不正确相交。相反，如果原点离表面太远，则阴影和反射可能会分离。

解决光线跟踪中此问题的典型做法是将生成的光线偏移固定的偏移值，而忽略沿光线p + td的交点比某个tmin值更近的情况。但这会有问题，如下图所示：

若反射角非常大，则此方法需要相当高的tmin值才能有效地工作。那么如果因为误差已经在表面上的点将偏移到更远的距离，这就使阴影大大的偏离了其正确位置。

我们利用浮点算法的基础思想，去了解产生误差的原因，以及如何保守的定位射线原点。以便永远不会发现不正确的自相交，同时将射线原点保持在非常接近实际相交点的位置，从而最大程度地减少错误的遗漏。我们不需要其他固定的偏移值。

浮点表示

IEEE标准指定以符号位表示32位浮点数，以指数表示8位，以有效数表示23位。对于8位，指数e的范围是0到255。实际使用的指数eb为：

实际浮点表达形式如下：

s是1位的符号位。m为23位有效位，且其拥有隐式前导位1(1.m表示)。可以用eb来替换e-127。例如，对于浮点数6.5被写为 $1.101_{2}*2^{2}$ ，其中下标2为2进制表示。小数点右边第i个数为 $1/2^{i}$ 因此6.5的计算过程如下：

eb=2，即e=129=10000001。m=10100000000000000000000。

因此，对于值6.5，该值的二进制内存表示为：

同样，浮点值1.0的二进制表示为：

根据这个表达形式，在整个范围内，两个相邻幂之间可表示的浮点之间的间距是均匀的(它可以表示有效位的一个增量)。在范围 $[2^{e},2^{e+1}]$ 中，间距为：

例如，对于介于1和2之间的浮点数，e = 0，并且浮点值之间的间隔为 $2^{-23}\approx 1.19209...*10^{-7}$ 。此间距也称为最后一个单位的大小(“ ulp”)。

到目前为止，我们已经描述了表示形式。那么0应该如何表示呢？因此将最小指数e=0或eb=-127留给特殊处理。有了这个指数，浮点值将被解释为在有效位数中没有隐式的前导位(即有效位不再是1.m，而是0.m)，这意味着有效位数所有位都是0会产生：

可以看到无论s位是正还是负，结果都为0。这就表示了IEEE标准的-0=0。若有前导位1(1.m)的情况下，最小的数为：

而消除前导位后，可表示的最小值为：

最大指数e = 255也保留用于特殊处理。因此，可以表示的最大常规浮点值具有e = 254（或eb = 127），大约为:

在e = 255时，如果有效位全为零，则根据符号位，该值对应于正无穷大或负无穷大。例如，在执行像浮点运算那样的运算时会产生无穷大的值。无穷大的算术运算导致无穷大。为了进行比较，正无穷大大于任何非无穷大，负无穷大也类似。

MaxFloat和Infinity常量分别初始化为可表示的最大值和“无穷大”浮点值。我们使它们可以在单独的常量中使用，因此使用这些值的代码无需使用冗长的C ++标准库调用来获取其值。

使用eb = 255时，非零有效位对应于特殊的NaN值，这些值是由诸如取负数的平方根或尝试计算0/0之类的运算得出的。NaN本身的任何算术运算始终返回NaN。所以程序运行结果出现NaN值，则表示中间过程某处出现了问题。与NaN值进行的任何比较均返回false，所以我们可以根据!(x==x)来判断该值是否为NaN。

程序实现

如果我们想利用位去操作浮点数，比如利用uint32_t来控制浮点的32位数。由于指针不能转换为另一种类型，我们可以使用memcpy()函数。该函数将类型转换为void*之后进行处理。我们可以利用以下代码完成浮点到位或者位到浮点的转换：


inline uint32_t FloatToBits(float f) {uint32_t ui;memcpy(&ui, &f, sizeof(float));return ui;
}inline float BitsToFloat(uint32_t ui) {float f;memcpy(&f, &ui, sizeof(uint32_t));return f;
}

这些转换可用于浮点值向上增大到下一个更大的浮点值，或向下减小到下一个较小的可表示浮点值。我们看一下递增的代码：

inline float NextFloatUp(float v) {// 处理∞和-0.0的情况if (std::isinf(v) && v > 0.) return v;if (v == -0.f) v = 0.f;// 获取间隔一个ulp的更大的浮点值uint32_t ui = FloatToBits(v);if (v >= 0)++ui;else--ui;return BitsToFloat(ui);
}inline float NextFloatDown(float v) {if (std::isinf(v) && v < 0.) return v;if (v == 0.f) v = -0.f;uint32_t ui = FloatToBits(v);if (v > 0)--ui;else++ui;return BitsToFloat(ui);
}

有两种重要的特殊情况：如果v是正无穷大，则此函数将不变地返回v。由于-0.0和0.0的位模式不相邻，所以负零会向前跳到正零，然后继续执行提高有效位数的代码。

给定浮点值，我们希望将有效位数增加1，如果结果溢出，则有效位数将重置为零，而指数则增加1。正好因为指数位于有效数字上方的高位，所以如果有效位数均全为1，然后将1加到有效数字的低位将导致向指数位进位。另请注意，当可表示的最大有限浮点值的位表示递增时，将得到正浮点无穷大的位模式。递减函数是相同的仅仅数值上相反。

浮点算术运算

使用带圆圈的运算符表示浮点算术运算，并使用sqrt表示浮点平方根，这些精度保证可以表示为：

round()位四舍五入操作。舍入误差的界限也可以用实数间隔表示：例如加法，我们可以说舍入结果在一个间隔内：

$\epsilon$ 为我们的一个ulp间隔。但我们引入的 $\epsilon$ 在精度上还是有所问题，如下图：

在加减一个 $\epsilon$ 的区间经常包含了两个浮点值，容易舍入到错误的浮点值。如果我们使用半个ulp间隔，便可以精确到一个浮点值。一个间隔步长为 $(a+b)2^{-23}$ 。半个间隔步长为 $(a+b)2^{-24}$ ，我们称其为machine epsilo，用 $\epsilon_{m}$ 表示， $\epsilon_{m}=2^{24}\approx 5.960464...*10^{-8}$ 。因此我们有：

PBRT中的Efloat类管理误差

PBRT中的Efloat类来简单管理误差舍入，其使用了区间来保守获得正确的射线原点。但其没有使用更精确的 $\epsilon_{m}$ 来计算区间，因为我们在计算表面交点时，保证交点不在表面下方，且不距离表面太远就能满足我们的需求了。过于精确的计算，可能带来太大的消耗，画面质量也不会有什么提升。Efloat最重要的三个私有成员float v, low, high;。其中v指计算后的值，low为区间下界，high为区间上界。每次进行新的计算，都需要跟新这三个值，我们来看对+号的运算符重载：

    EFloat operator+(EFloat ef) const {EFloat r;r.v = v + ef.v;
#ifndef NDEBUGr.vPrecise = vPrecise + ef.vPrecise;
#endif  // DEBUG// Interval arithemetic addition, with the result rounded away from// the value r.v in order to be conservative.r.low = NextFloatDown(LowerBound() + ef.LowerBound());r.high = NextFloatUp(UpperBound() + ef.UpperBound());r.Check();return r;}

其中vPrecise是有高精度要求才要考虑的操作，用double对其保存。我们可以看到在执行加法操作时，首先对两个Efloat数的v值求和，然后让其下界和上界分别相加。获得新的上界需要重新获取最近的(间隔一个ulp)更大的float值，新的下界也要获取最近的(间隔一个ulp)新的更小的float值。该方法能更保守的确定精确值就在该区间内。Check()函数来检查上下界是否出现NaN值或者无穷情况。

减法操作也是类似，但是要获取保守的上界和下界，要获得下界需要让被减数的下界减去减数的上界，获得上界则是用被减数上界减去被减数的下界。代码如下：

    EFloat operator-(EFloat ef) const {EFloat r;r.v = v - ef.v;
#ifndef NDEBUGr.vPrecise = vPrecise - ef.vPrecise;
#endifr.low = NextFloatDown(LowerBound() - ef.UpperBound());r.high = NextFloatUp(UpperBound() - ef.LowerBound());r.Check();return r;}

接下来看乘法和除法的操作：

EFloat operator*(EFloat ef) const {EFloat r;r.v = v * ef.v;
#ifndef NDEBUGr.vPrecise = vPrecise * ef.vPrecise;
#endifFloat prod[4] = {LowerBound() * ef.LowerBound(), UpperBound() * ef.LowerBound(),LowerBound() * ef.UpperBound(), UpperBound() * ef.UpperBound()};r.low = NextFloatDown(std::min(std::min(prod[0], prod[1]), std::min(prod[2], prod[3])));r.high = NextFloatUp(std::max(std::max(prod[0], prod[1]), std::max(prod[2], prod[3])));r.Check();return r;}EFloat operator/(EFloat ef) const {EFloat r;r.v = v / ef.v;
#ifndef NDEBUGr.vPrecise = vPrecise / ef.vPrecise;
#endifif (ef.low < 0 && ef.high > 0) {分母为0的情况r.low = -Infinity;r.high = Infinity;} else {Float div[4] = {LowerBound() / ef.LowerBound(), UpperBound() / ef.LowerBound(),LowerBound() / ef.UpperBound(), UpperBound() / ef.UpperBound()};r.low = NextFloatDown(std::min(std::min(div[0], div[1]), std::min(div[2], div[3])));r.high = NextFloatUp(std::max(std::max(div[0], div[1]), std::max(div[2], div[3])));}r.Check();return r;}

乘法和除法的保守方式是分别用上下界和被操作数的上下界进行计算，将最小值赋给low，最大值赋给high。值得注意的是除法要处理分母为0的情况。

误差传播

测量误差有两种方法是有用的：绝对误差和相对误差。我们计算机计算的四舍五入后的结果 $\widetilde{a}$ 和实际值a之间是不同的我们称绝对误差 $\delta _{a}$ 为:

相对误差 $\delta _{r}$ 是绝对误差和实际值的比：

使用相对误差的定义，我们可以将计算值 $\widetilde{a}$ 写成精确结果a的关系式：

如果计算r =（（（（a + b）+ c）+ d）：

有如下约束：

可以重新获得式子：

其中误差间隔为：

可以发现，如果我们在执行加法时，按照从小到达排序，可以使误差间隔降到最低。

如果计算r = (a + b) + (c + d）这样的类型呢？根据：

利用上述同样的方法计算，可以得出：

可以发现下面这种计算过程在d较小，a和b较大的情况下比上面的过程得出的误差间隔小很多。根据计算前值的不同，选择不同的方式可以获得更好的误差间隔。

Higham(2002)给出了一种更好的误差约束。他用 $1+\theta _{n}$ 来约束 $(1\pm \epsilon _{m})^{n}$ 。其中：

用 $\gamma_{n}$ 采用该边界：

获取该值的代码如下：

inline constexpr Float gamma(int n) {return (n * MachineEpsilon) / (1 - n * MachineEpsilon);
}

该方法PBRT在一些专门处理射线求交时会使用，但是在efloat本身的误差积累种并没有采用这种方法。虽然可能随着计算误差会逐渐增加，但并不会对输出图像造成不可忍受的影响。

值得注意的是，每次进行新的运算都会在原有误差基础上产生一个新误差： $(1\pm\ \epsilon_{m} )^{1}$ ，即 $(1\pm\gamma _{1} )$ ，我们来看乘法运算的例子：

加法：

加法可以约束在间隔：内，若a和b有相同的符号，则他们可以共享相同的 $\gamma _{n}$ 。

专用的射线求交误差管理方式

边界盒求交的浮点误差

在边界盒求交时，PBRT引入了 $\gamma_{n}$ 来避免因为误差导致无法判断相交的情况。计算光线和边界盒交点时，沿着光线进入边界盒的光线和离开光线的tmax找出沿光线的参数tmin。如果tmin < tmax，则光线穿过盒子；否则会错过它。使用浮点算法时，计算出的t值可能会有误差，如果纯粹由于舍入误差而导致计算出的tmin值大于tmax时，相交测试将错误地返回false。

用浮点运算计算光线交点t有：

同样可以写成：

如果在误差间隔内，tmin和tmax不重叠，那么肯定可以正确返回，但如果重叠了，则可能出现错误情况，如下图：

如果间隔确实重叠，则不可能知道t值的实际顺序。在这种情况下，在执行比较之前，将tmax增加两倍误差范围2g3tmax，可以确保在这种情况下我们保守地返回true。代码执行：

tFar *= 1 + 2 * gamma(3);

我们可以完整的看一下射线和边界盒(AABB)求交的代码：

template <typename T>
inline bool Bounds3<T>::IntersectP(const Ray &ray, const Vector3f &invDir,const int dirIsNeg[3]) const {const Bounds3f &bounds = *this;// 检查光线和x平面和y平面的相交情况Float tMin = (bounds[dirIsNeg[0]].x - ray.o.x) * invDir.x;Float tMax = (bounds[1 - dirIsNeg[0]].x - ray.o.x) * invDir.x;Float tyMin = (bounds[dirIsNeg[1]].y - ray.o.y) * invDir.y;Float tyMax = (bounds[1 - dirIsNeg[1]].y - ray.o.y) * invDir.y;// 利用gamma(3)确保在tmin和tmax在同一个间隔内返回truetMax *= 1 + 2 * gamma(3);tyMax *= 1 + 2 * gamma(3);if (tMin > tyMax || tyMin > tMax) return false;if (tyMin > tMin) tMin = tyMin;if (tyMax < tMax) tMax = tyMax;// 检查光线和x平面的相交情况Float tzMin = (bounds[dirIsNeg[2]].z - ray.o.z) * invDir.z;Float tzMax = (bounds[1 - dirIsNeg[2]].z - ray.o.z) * invDir.z;// 利用gamma(3)确保在tmin和tmax在同一个间隔内返回truetzMax *= 1 + 2 * gamma(3);if (tMin > tzMax || tzMin > tMax) return false;if (tzMin > tMin) tMin = tzMin;if (tzMax < tMax) tMax = tzMax;return (tMin < ray.tMax) && (tMax > 0);
}

二次曲面重投影

在计算射线和二次曲面的交点时，也会引入误差问题，重新计算相交会在计算上过于昂贵。这里可以引入重投影方法。虽然重投影也会有误差，但该方法产生的误差通常非常小。重投影的重点是检测可能导致重投影点落在表面以下的误差。

考虑射线-球相交：利用已计算的 $P=(x , y,z)$ 去重新投影计算 $P'=(x' , y',z')$ 。以计算x为例：

由于在对象空间，球体球心坐标为(0,0,0)，其半径为r。我们x'的浮点运算过程为：

由于 $x^{2},y^{2},z^{2}$ 都为正，所以可以合并 $\gamma$ ：

所以射线-球求交的误差计算代码为：

//pHit为vec3类型的击中点
pHit *= radius / Distance(pHit, Point3f(0, 0, 0));
Vector3f pError = gamma(5) * Abs((Vector3f)pHit);

三角形求交误差管理

三角形计算交点会使用他们三条边的边缘函数e0，e1，e2的值。重心坐标 $b_{i}$ 由这三个值计算得出：

三角形v'可以用重心坐标和三角形三个顶点计算得出：

我们v'的误差处理过程为：首先计算：

因为击中点在三角形内，所以e0，e1，e2有相同的符号，所以有：

如果现在考虑对三角形中与边缘函数值相对应的位置的x坐标进行插值，则有：

带入上述d求得的误差间隔，则有：

为了方便计算，我们做微小的误差扩张，则可计算出绝对误差在以下间隔内：

保守误差间隔为：

因此三角形求交的误差管理代码为：

Float xAbsSum = (std::abs(b0 * p0.x) + std::abs(b1 * p1.x) +std::abs(b2 * p2.x));
Float yAbsSum = (std::abs(b0 * p0.y) + std::abs(b1 * p1.y) +std::abs(b2 * p2.y));
Float zAbsSum = (std::abs(b0 * p0.z) + std::abs(b1 * p1.z) +std::abs(b2 * p2.z));
Vector3f pError = gamma(7) * Vector3f(xAbsSum, yAbsSum, zAbsSum);

形变矩阵的误差影响

变换矩阵transform在程序运行中是经常使用的。矩阵运算，即矩阵和vec3的点相乘会对x，y，z分别进行三次乘法运算并求和。这就自然的引入了误差。我们直接看浮点运算公式(以x为例,其次坐标w为1)：

计算出的绝对误差可由如下式子约束：

y'，z'可以同理得出。矩阵变换transform的代码为：

template <typename T> inline Point3<T>
Transform::operator()(const Point3<T> &p, Vector3<T> *pError) const {T x = p.x, y = p.y, z = p.z;//计算变换后的点pT xp = m.m[0][0] * x + m.m[0][1] * y + m.m[0][2] * z + m.m[0][3];T yp = m.m[1][0] * x + m.m[1][1] * y + m.m[1][2] * z + m.m[1][3];T zp = m.m[2][0] * x + m.m[2][1] * y + m.m[2][2] * z + m.m[2][3];T wp = m.m[3][0] * x + m.m[3][1] * y + m.m[3][2] * z + m.m[3][3];//计算误差T xAbsSum = (std::abs(m.m[0][0] * x) + std::abs(m.m[0][1] * y) +std::abs(m.m[0][2] * z) + std::abs(m.m[0][3]));T yAbsSum = (std::abs(m.m[1][0] * x) + std::abs(m.m[1][1] * y) +std::abs(m.m[1][2] * z) + std::abs(m.m[1][3]));T zAbsSum = (std::abs(m.m[2][0] * x) + std::abs(m.m[2][1] * y) +std::abs(m.m[2][2] * z) + std::abs(m.m[2][3]));*pError = gamma(3) * Vector3<T>(xAbsSum, yAbsSum, zAbsSum);if (wp == 1) return Point3<T>(xp, yp, zp);else         return Point3<T>(xp, yp, zp) / wp;}

安全的射线起点

计算出的交点及其误差范围为形成一个小3D框，该框限制了空间区域。我们知道精确的相交点必须在此盒子内部某处，因此曲面必须穿过盒子（至少足以显示相交的点）。如下图所示：

用这些框可以定位离开表面的射线的原点，使其始终位于表面的外侧，以免它们不正确地重新自相交。当追踪生成的光线离开交点p时，我们会偏移它们的原点，以确保它们超出错误框的边界，因此不会错误地与曲面重新自相交。

为了确保产生的射线源一定在表面的外侧，我们沿着法线移动了足够的距离，以使垂直于法线的平面在误差边界框之外。其中平面方程：

f(x,y,z)=0，表示了经过交点的平面。平面方程f的值能表示沿法线的偏移，该偏移给出了经过该点的平面。我们想在误差边界框的八个角处找到f的最大值；如果我们对平面加上和减去此偏移量进行偏移，则有两个不与误差框相交的平面。我们误差边界框由( $(\pm \delta _{x},\pm \delta _{y},\pm \delta _{z})$ )给出，最大的距离可以非常简单的计算出：

实现代码：

inline Point3f OffsetRayOrigin(const Point3f &p, const Vector3f &pError,const Normal3f &n, const Vector3f &w) {Float d = Dot(Abs(n), pError);Vector3f offset = d * Vector3f(n);if (Dot(w, n) < 0)offset = -offset;Point3f po = p + offset;//偏移后的保守误差处理 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]);}return po;
}

代码中，d由法线n的绝对值和误差的点乘计算得出。我们偏移后依旧会引入误差，因为执行了加号计算。所以为了获取保守边界，我们在将po分别位移一个ulp的上下界。

到此，阴影光线与区域光源仍存在一个相关的问题：我们希望找到形状与光源非常接近的任何相交并实际上将其遮挡，同时避免报告与光源表面的不正确相交。因此我们将阴影光线的tMax值设置为小于1，以便它们在光源表面之前停止。于是我们发射射线的代码为：

const Float ShadowEpsilon = 0.0001f;
Ray SpawnRayTo(const Point3f &p2) const {Point3f origin = OffsetRayOrigin(p, pError, n, p2 - p);Vector3f d = p2 - origin;return Ray(origin, d, 1 - ShadowEpsilon, time, GetMedium(d));
}

由于可能计算好的射线还会执行transform变换，所以我们还需要在transform变换射线时，处理误差情况，代码如下：

inline Ray Transform::operator()(const Ray &r) const {Vector3f oError;Point3f o = (*this)(r.o, &oError);Vector3f d = (*this)(r.d);// Offset ray origin to edge of error bounds and compute _tMax_Float lengthSquared = d.LengthSquared();Float tMax = r.tMax;if (lengthSquared > 0) {Float dt = Dot(Abs(d), oError) / lengthSquared;o += d * dt;tMax -= dt;}return Ray(o, d, tMax, r.time, r.medium);
}

我们给起点沿射线方向偏移一个保守距离，以及让tmax减去一个保守距离。保证我们的射线足够安全。