计算反射/折射的屏幕空间是否需要二阶偏导数？

Question

我编写了一个基本的光线追踪器来跟踪屏幕空间。每个片段都有一个相关的像素半径。当射线dir 被一些distance 从eye 挤压时撞击几何体，我计算该撞击的法线向量N，并将其与另外四个射线组合。在伪代码中：

def distance := shortestDistanceToSurface(sdf, eye, dir, pixelRadius)
def p := eye + dir * distance
def N := estimateNormal(sdf, p)
def glance := distance * glsl.dot(dir, N)
def dx := (dirPX / glsl.dot(dirPX, N) - dirNX / glsl.dot(dirNX, N)) * glance
def dy := (dirPY / glsl.dot(dirPY, N) - dirNY / glsl.dot(dirNY, N)) * glance

这里，dirPX、dirNX、dirPY 和 dirNY 是在屏幕空间的四个方向中的每一个方向上偏移了dir 像素半径的光线，但仍然瞄准相同参照点。这给出了dx 和dy，它们是像素上的偏导数，表示当光线穿过屏幕空间时命中沿几何体表面移动的速率。

因为我跟踪屏幕空间，我可以使用pre-filtered samplers, as discussed by Inigo Quilez。他们看起来很棒。但是，现在我想添加反射（和折射），这意味着我需要递归，我不确定如何计算这些光线并跟踪屏幕空间。

根本的问题是，为了弄清楚几何体上某个位置反射了什么颜色的光，我不仅需要采集点样本，还需要检查反射的整个屏幕空间。我可以使用偏导数在几何上给我四个新点，它们近似于一个椭圆，这是原始像素从屏幕上的投影：

def px := dx * pixelRadius
def py := dy * pixelRadius
def pPX := p + px
def pNX := p - px
def pPY := p + py
def pNY := p - py 

我可以通过将椭圆弄成圆形来计算近似的像素半径。我知道这会破坏某些理想的各向异性模糊，但是没有作弊的光线追踪器是什么？

def nextRadius := (glsl.length(dx) * glsl.length(dy)).squareRoot() * pixelRadius

但是，我不知道在哪里将这些点反映到几何中；我不知道将它们的光线聚焦在哪里。如果我必须选择焦点，那么它将是任意的，并且取决于几何体反映其自身图像的位置，这可能会任意模糊或摩尔反射图像。

我需要取二阶偏导数吗？我可以像一阶导数一样逼近它们，然后我可以用它们来调整正常的N，只需稍加改变，就像点击p。然后法线引导椭圆的焦点，并将其映射到近似圆锥截面。我担心三件事：

1. 我担心做一些额外的向量加法和乘法的成本，这可能可以忽略不计；
2. 还有关于精度损失（在进行这些廉价衍生品时已经非常糟糕）是否会在多次反射中损失过大；
3. 最后，我应该如何处理屏幕空间爆炸的情况；当我有一个镜像球体时，我应该如何在反射空间的大楔形上进行采样，例如将棋盘格图案平均成灰色？

虽然不用担心，但我根本不知道如何获取四个向量并快速为它们拟合一个令人信服的锥体，但这可能只是花一些时间在白板上做代数的问题。

编辑：在 John Amanatides 1984 年的论文 Ray Tracing with Cones 中，曲率信息确实被计算出来，并用于将估计的锥体拟合到反射光线上。在 Homan Igehy 1999 年的论文Tracing Ray Differentials 中，仅使用一阶导数，而明确忽略了二阶导数。

可能还有其他选择吗？我已经尝试在一次反射后丢弃像素半径并仅采集点样本，它们看起来很糟糕，有很多锯齿和噪音。也许存在可以基于每种材料计算的视场或景深近似值。像往常一样，多重采样会有所帮助，但我想要一个分析解决方案，这样我就不会不必要地浪费这么多 CPU。

（sdf 是 signed distance function，我正在做球体追踪；同样的程序既计算距离又计算法线。glsl 是 GLSL 标准库。）

Answer 1

我不会接受我自己的答案，但我会解释一下我做了什么，以便我可以暂时放下这个问题。我最终采用了类似于 Amanatides 的方法，计算每条射线周围的锥体。

每次计算法线向量时，我也会计算平均曲率。使用众所周知的技巧计算法线向量。让p 成为我的问题中的 SDF 的零，让epsilon 成为用于numerical differentiation 的合理偏移量，然后让vp 和vn 成为向量，其分量是附近 SDF 的评估p 但每个组件都受到epsilon 的干扰。在伪代码中：

def justX := V(1.0, 0.0, 0.0)
def justY := V(0.0, 1.0, 0.0)
def justZ := V(0.0, 0.0, 1.0)
def vp := V(sdf(p + justX * epsilon), sdf(p + justY * epsilon), sdf(p + justZ * epsilon))
def vn := V(sdf(p - justX * epsilon), sdf(p - justY * epsilon), sdf(p - justZ * epsilon))

现在，通过巧妙地滥用finite difference coefficients，我们可以同时计算一阶导数和二阶导数。第三个系数sdf(p)，假设已经为零。这给出了我们的法线向量N，即梯度，以及我们的平均曲率H，即Laplacian。

def N := glsl.normalize(vp - vn)
def H := sum(vp + vn) / (epsilon * epsilon)

我们可以从平均曲率估计高斯曲率。虽然平均曲率告诉我们的锥体扩大或收缩多少，但高斯曲率始终是非负的，并且测量了在锥体的相交区域中添加了多少额外面积 (spherical excess)。用K代替H给出高斯曲率，代入后：

def K := H * H

现在我们准备调整片段宽度的计算。让我们假设，除了屏幕空间中的pixelRadius 和从相机到几何体的distance，我们还有dradius，像素半径随距离的变化率。我们可以像以前一样用点积来计算扫视因子，三角函数也是类似的：

def glance := glsl.dot(dir, N).abs().reciprocal()
def fradius := (pixelRadius + dradius * distance) * glance * (1.0 + K)
def fwidth := fradius * 2.0

现在我们有了fwidth，就像在 GLSL 中一样。最后，当需要反思时，我们需要通过将二阶导数曲率积分到一阶导数中来调整半径的变化：

def dradiusNew := dradius + H * fradius

结果令人满意。碎片可能有点太大了；很难判断某些东西是否过于模糊或只是适当地抗锯齿。我想知道 Amanatides 是否使用一组不同的方程来处理曲率；我在白板上花费了太长时间，最终得到了令人愉快的简单操作。

此图像未经过超采样；每个像素都有一个片段，一个射线和一个锥体。