最近一个月重写了我第一篇论文《Sheared Polygonal Texture filtering》的OpenGL实现，在这个过程中总结了一些关于float精度的缺陷。为了便于了解这些问题，首先简单介绍一下我的这篇论文，然后在针对float精度引起了问题做一个总结

Texture Mapping

如图，当我们看屏幕时，屏幕上的像素投影到纹理空间，其中的透视投影Perspective Projection称为纹理映射，这样，屏幕空间的正方形像素在纹理空间对应了一个四边形，我们称该四边形为该像素的footprint，其中可能覆盖了部分的texels。当然，这里有一个有意思的讨论：像素是什么？这里假设一个像素是正方形，但也可能是圆形，或者液晶屏下是RGB的阵列等不同的信号，What is a pixel, this is a question.

Texture filter

texture mapping之后，我们需要确定该像素对应的颜色或者值。纹理空间下的任一纹理，我们获取该footprint下的所有texels的weighted value，然后赋给该像素，这个过程称为texture filter。

通常，我们会通过super sample遍历所有的texel来求解这个过程，这对应的是一个点查询。而目前常用的Mipmap或Anisotropic Mipmap则是找到该footprint对应的一个规律形状，比如正方形或矩形，预计算中我们已经保存了该规则形状的值，所以只需要一次面查询就可以得到该值。所以，点查询和基于预计算的面查询在性能和质量之间拉扯，顾此失彼，此事古难全。而EWA则先找到对应的椭圆形状，然后进行点查询，通常是一种高质量且性能可以接受的方式，通常应用在离线渲染中，还不能满足实时渲染的要求。

这也是论文的动机，既然点查询和面查询各有利弊，那是否存在一种线查询，能够兼顾质量和性能？

SPTF & SSAT

如此，我们的问题还是求解四边形的权重面积，一个面积积分的问题。而格林公式告诉我们，面积分可以转为线积分求解，只需要我们在四边形的edge走一圈，对应的directed area求和。这就是该论文的方法 SPTF，这样，一个面积分变成了四个边的线积分。问题变成了求解任意线（段）。

常见的SAT告诉我们如果求解一条水平的直线，在成都生活多年的我，自然想到了，只需要一个shear变换，我也能把斜的掰直。同时因为边界外都是空白区域，黄色梯形区域的面积和等同于其对应的矩形，所以，只需要一个shear，我们就可以通过SAT求解任意线段，而我们只需要预计算这个Sheared SAT，方便我们在SPTF中查询对应的SSAT。

Terrible float

这个算法本身是简单易懂的，包括两大部分：预计算时的SSAT和实时渲染中的SPTF两部分。在应用中需要考虑连续空间的算法如何在离散下的实现问题 ”Continuous Math on Discrete Hardware“。下面，我会给出应用中遇到的Float precision相关的缺陷。之前公众号的文章《Precision？Precision！》中介绍过Float的规范，有兴趣的可以查阅。

HDR

首先，我在 GPU 中实现了并行的 Horizontal SAT 列求和，思路来自 Fast SAT 那篇论文。为了验证结果正确性，我使用了一张全红纹理（值为 1）进行测试，并将 SAT 结果导出为 HDR 格式进行检查。小尺寸纹理下结果完全正确，但当换成 1024 的大纹理后，底部数值较大区域开始出现错误。

一开始我怀疑是 texture coordinate 的 UV 精度问题。因为小纹理下精度足够，而大纹理可能因 UV 的细微误差导致采样偏移，于是我改用 texture fetch，结果问题依旧存在。随后又排查了纹理越界读取、OpenGL 浮点纹理格式设置等各种可能性，依然毫无结果。这个问题耗费了我大半天。

后来我发现，错误会在数值超过 $256$ 后开始出现。进一步调试时，发现 GPU 内存中的变量其实是正确的，这让我一度怀疑是 HDR 看图软件显示错误。接着我想到，会不会是 HDR 格式本身的精度限制导致了问题。

原因则是在HDR中采用RGBE的格式:

\(V_{max} = \max(R, G, B)\) \(M_{byte} = \text{floor}\left( \frac{V \times 256}{2^{E_{real}}} \right)\)

当$V=257$, pack:

\[E_{real} = \lceil \log_2(257) \rceil = \lceil 8.005 \rceil = 9\] \[E_{byte} = E_{real} + 128 = 137\] \[M_{byte} = \text{floor}(\frac{257 \times 256}{512}) = \text{floor}(128.5) = 128\]

unpack:

\[V = \frac{M_{byte}}{256} \times 2^{(E_{byte} - 128)} = 256\]

于是我改用 EXR 格式导出，结果一切正常。

在此之前，我一直以为 HDR 和 EXR 没有本质区别。严格来说，这并不是真正的 float 运算精度问题，更像一次人为事故，但其原理与 float 精度损失非常相似，也更容易帮助理解浮点数为何会出现精度误差。最终原因看起来简单得不可思议，却恰恰属于那种视而不见的问题。人类调试时最大的敌人，通常不是复杂性，而是万万想不到的自信。

floor()

错误现象如上图左侧高亮区域所示。后来调试发现，只有斜率为 对应的 footprint 会出现问题，其他斜率对应的形状都能得到正确面积。这一度让我开始认真思考“宇宙终极答案是不是 42”这种本不该出现在 GPU 调试过程里的哲学问题。毕竟，当一个 bug 只在某个特定浮点数上触发时，人类和玄学之间的距离就已经非常危险了。

更离谱的是，如果强制所有像素都使用斜率 $0.4$ 对应的 footprint，就会得到右图那种灾难性的结果。

当然，如果这是产品代码，其实完全可以加一个判断，直接屏蔽 $0.4$ 的情况。毕竟 GPU 调试本身就是一种精神消耗：要么靠猜测一个可能性，要么穷举所有可能性。而一个只在特定条件触发的缺陷，究竟值不值得继续深挖，本身就是工程问题的一部分。

万幸，我听过那首歌里唱的：“没有确定的以后，没有谁祝福我，反而想要勇敢接受。”很多事情都这样，魔鬼永远藏在细节里。就像周五同事说，他分手了，很难过，不知道该干什么。我笑着说，不管怎样，还是应该战术性地“舔”一下，至少表达你对这段感情的认真和执着。其实内心还是有点羡慕的。年轻真好，会把爱情看得那么伟大。

原因如上。float并不能精确表示 $0.05$，而$0.4$ 恰好又是 $0.05$ 的八倍。于是问题出现了：

在生成 SSAT 时，参与计算的 $0.4$ 实际上比真实值略小一点，因此算出的 offset 是 $384$；而在最终渲染阶段，使用的 $0.4$ 又比真实值略大一点，于是算出的 offset 变成了 $385$。读写位置因此发生错位，结果自然完全错误。

本质上，这并不是什么复杂算法错误，而是典型的浮点数精度误差导致的边界问题。只是它隐藏得极深，因为两边看起来都“像是” $0.4$。可惜 GPU 不相信“差不多”，它只认二进制位。所谓“失之毫厘，谬以千里”，大概就是如此。很多时候，真正让人难受的并不是彻底失败，而是只差一点点。就像 offset 从 $384$ 变成 $385$，就像人生里那些“如果当时再多一点”的瞬间。明明无限接近正确答案，却终究不是那个答案。于是，我在生成SSAT时把float改为double，从而保证offset的误差小于1，不会妨碍在offset整型的尺度。哎，人生如果能够调试，那该多好，如果你知道如何穿越，请你一定告诉我。

round()

错误就出现在这张纹理中的黑色区域，那里反而会产生类似 firefly 的亮点噪声。到了这个阶段，我已经搭建了一整套完善的调试系统：输出 N 张中间纹理，记录每个像素的位置、构成四边形的顶点、对应斜率，以及 GPU 算法的 CPU 对照版本。这样一来，我终于可以精确定位问题像素的 POS，并逐步回溯整个计算过程。

原因则是我用UINT替代float时的round问题。如上图，红线下的面积应该大于橙线下的面积，而在该纹理中有一段区域是黑色，所以应该是大于等于。但在 float 转 UINT 的过程中，四舍五入产生了 1 个单位的误差。橙线“五入”，红线“四舍”，导致下面的面积反而比上面大了 1。随后程序继续执行： 而由于使用的是 UINT，$-1$ 会直接下溢成一个极大的无符号整数，于是屏幕上便炸出了那些耀眼的 firefly。

这个问题，要解决就有点龌龊了： max(红，橙) - 橙。

Fixed point

遇到了各种各样的问题，人也会在挫折里慢慢长大，不再 “too simple, sometimes naive”。

后来，当我把 step size 调得非常小时，希望支持更多斜率来构造更精确的 footprint，又遇到了上图左侧那种错误。这一次，我第一反应就是：会不会是 step size 太小，导致在寻找最近斜率时出现了 float 精度问题？

例如，$0.1249999$本应属于：$0.12$,却因为浮点误差，被错误匹配成了$0.13$。问题可能并不在 SAT、本身算法或者积分面积，而是在最开始“寻找最近斜率”这一步，就已经选错了 footprint。

于是我把浮点数改成了定点数。虽然 step size 变得更小，但定点数能够精确表达这些离散值，因此成功绕开了这个问题。

终于，这次猜对了。

问题可能以分钟为单位很快解决，可我却并不开心。因为这种调试过程最可怕的地方就在于：你开始越来越熟悉 bug 的思维方式。你不再相信“应该没问题”，而是条件反射般地怀疑每一个 round、每一个 cast、每一个 0.000001。这种感觉就好像曾经爱过的那些歌曲，有一天再也不听，我想你能体会到这种失落吧。这个世界有时候确实糟糕透了，就好像每次看到儿子怎么这么笨时，都会小声的骂一句傻逼，然后默默的希望他永远的不要长大，我也怀念我的笨。

Line intersection

上图是问题的效果，小 quad footprint 会出现噪点，而 semi-parallelogram 完全稳定。这让问题看起来像是几何建模层面的理论误差，比如 shear 之后 UV 不再正交，线性插值去近似非线性变换，从而引入系统性偏差。而这个误差是很难量化的，无从分析也就无法判断问题是否可解。

后来意外发现问题出在线线求交上。

我们用 $(p,d)$ 表示一条线，其中 $p$ 是点，$d$ 是方向，也可写成 $ax+by+c=0$，再通过两线求交计算交点。算法本身是正确的，但在 $p$ 很大、两点又非常接近时（比如 4K 纹理中一个像素对应极小局部区域），浮点误差会被显著放大。

结果就是：线线求交中累积的数值误差，相对于两点间距离已经不可忽略，最终导致计算出的交点位置与 slope/step size 的分类不一致，从而引入误差。

于是，我们通过如上图的公式求解，通过两个p点的相对距离，仅需要较少的计算就可以获取该交点，误差也更小。

fma

而另一个方法则是fused multiply-add (fma)。这个方法是真好，如果你的数值计算中有ax+b的，大胆的用fma吧，又快有好的事情在这个物质世界已经不多了。如果哪天要纹身，我可能会考虑的一个公式就是$Ax+b$。如上，如果我用fma替代第一个版本的线线求交算法，也会解决该bug。

这一整个月踩的坑，本质都绕不开 float 精度问题。每个问题都不大，但不记录，下次遇到基本还是会再掉进去一次。

而且这些问题还挺“跨平台玄学”的：不同硬件、不同 GPU、甚至升级个 PyTorch 版本，float 行为都可能悄悄变一点。

所以这一个月最大的收获不是修 bug，而是把这些“看起来不值一提的小误差”系统性整理了一遍。 结论也很朴素：单元测试很重要。尤其是在你已经开始不信任人类时。

Whenever you use float, you create error.