PBRT_V2 总结记录 <45> Texture Sampling Rate 和 DifferentialGeometry->ComputeDifferentials

理论基础：

Fortunately, for texture antialiasing, we don’t need to be able to evaluate f (x, y) for arbitrary (x, y), but just need to find the relationship between changes in pixel sample position and the resulting change in texture
sample position at a particular point on the image. This relationship is given by the partial
derivatives of this function, ∂f/∂x and ∂f/∂y. For example, these can be used to find
a first-order approximation to the value of f ,

If these partial derivatives are changing slowly with respect to the distances x’ − x and
y‘ − y, this is a reasonable approximation. More importantly, the values of these partial
derivatives give an approximation to the change in texture sample position for a shift of
one pixel in the x and y directions, respectively, and thus directly yield the texture sampling
rate.

（实现Texture的抗锯齿，我门需要去发现，假如改变了像素采样点的位置， 那么会导致 纹理采样点位置 发生怎么样的变化 ，个人理解就是，当像素采样点的间距变化了，那么纹理采样点的间距会怎么变化呢？需要找出这两者之间的关系，

那么 偏导数就可以 说明到这样的关系，例如之后会 计算 ∂u/∂x,  ∂u/∂y, ∂v/∂x,  ∂v/∂y, 就 ∂u/∂x 偏导数的 意义就是，当前像素的采样点 偏移了 Δx, 那么 对应 偏移了多少 Δu，这样就可以把 像素采样点的位置变化 与 纹理采样点的位置变化关联起来了）

 

 

实现：

DifferentialGeometry 类：

// DifferentialGeometry Declarations
struct DifferentialGeometry {
    DifferentialGeometry() { 
        u = v = dudx = dvdx = dudy = dvdy = 0.; 
        shape = NULL; 
    }
    // DifferentialGeometry Public Methods
    DifferentialGeometry(const Point &P, const Vector &DPDU,
            const Vector &DPDV, const Normal &DNDU,
            const Normal &DNDV, float uu, float vv,
            const Shape *sh);
    void ComputeDifferentials(const RayDifferential &r) const;

    // DifferentialGeometry Public Data
    Point p;
    Normal nn;
    float u, v;
    const Shape *shape;
    Vector dpdu, dpdv;
    Normal dndu, dndv;
    mutable Vector dpdx, dpdy;
    mutable float dudx, dvdx, dudy, dvdy;
};

 

1. 
void DifferentialGeometry::ComputeDifferentials(
        const RayDifferential &ray) const 

作用：

（ComputeDifferentials 方法里面，主要是计算 dpdx, dpdy， dudx, dvdx, dudy, dvdy， 其他的 dpdu, dpdv, dndu, dndv ,都是在  virtual bool Shape::Intersect(const Ray &ray, float *tHit, float *rayEpsilon, DifferentialGeometry *dg) const 中计算，需要注意一下这个方法的调用时机，在 Intersection::GetBSDF() 里面直接调用，但是 在 Material GetBSDF()之前，主要是在使用Texture之前，先把数据计算好。

BSDF *Intersection::GetBSDF(const RayDifferential &ray,
                            MemoryArena &arena) const {
    PBRT_STARTED_BSDF_SHADING(const_cast<RayDifferential *>(&ray));
    dg.ComputeDifferentials(ray);
    BSDF *bsdf = primitive->GetBSDF(dg, ObjectToWorld, arena);
    PBRT_FINISHED_BSDF_SHADING(const_cast<RayDifferential *>(&ray), bsdf);
    return bsdf;
}

）

The DifferentialGeometry::ComputeDifferentials() method computes these values.  dpdx, dpdy，dudx, dvdx, dudy, dvdy;
It is called by Intersection::GetBSDF() and Intersection::GetBSSRDF() before the
Material GetBSDF() or GetBSSRDF() method is called so that these values will be available
for any texture evaluation routines that are called by the material. Because ray differentials
aren’t available for all rays traced by the system

 

2.  dpdx, dpdy， dudx, dvdx, dudy, dvdy 表示：

（dx，dy 表示的是 屏幕x轴，y轴偏移一个 像素采样点，

dpdx，dpdy 表示 屏幕中x，y轴偏移一个像素采样点，对应世界坐标点p 会怎么偏移，

同样的道理， dudx, dvdx, dudy, dvdy 表示 屏幕中x，y轴偏移一个像素采样点， 对应的 texture 的 (u,v) 坐标怎么偏移）

Because the differential rays are offset one pixel sample in each direction, there’s no need
to divide these differences by a Δ value, since Δ = 1.

 

3. void DifferentialGeometry::ComputeDifferentials(
        const RayDifferential &ray) const 

代码：

void DifferentialGeometry::ComputeDifferentials(
        const RayDifferential &ray) const {
    if (ray.hasDifferentials) {
        // Estimate screen space change in $\pt{}$ and $(u,v)$

        // Compute auxiliary intersection points with plane
        float d = -Dot(nn, Vector(p.x, p.y, p.z));
        Vector rxv(ray.rxOrigin.x, ray.rxOrigin.y, ray.rxOrigin.z);
        float tx = -(Dot(nn, rxv) + d) / Dot(nn, ray.rxDirection);
        if (isnan(tx)) goto fail;
        Point px = ray.rxOrigin + tx * ray.rxDirection;
        Vector ryv(ray.ryOrigin.x, ray.ryOrigin.y, ray.ryOrigin.z);
        float ty = -(Dot(nn, ryv) + d) / Dot(nn, ray.ryDirection);
        if (isnan(ty)) goto fail;
        Point py = ray.ryOrigin + ty * ray.ryDirection;
        dpdx = px - p;
        dpdy = py - p;

        // Compute $(u,v)$ offsets at auxiliary points

        // Initialize _A_, _Bx_, and _By_ matrices for offset computation
        float A[2][2], Bx[2], By[2];
        int axes[2];
        if (fabsf(nn.x) > fabsf(nn.y) && fabsf(nn.x) > fabsf(nn.z)) {
            axes[0] = 1; axes[1] = 2;
        }
        else if (fabsf(nn.y) > fabsf(nn.z)) {
            axes[0] = 0; axes[1] = 2;
        }
        else {
            axes[0] = 0; axes[1] = 1;
        }

        // Initialize matrices for chosen projection plane
        A[0][0] = dpdu[axes[0]];
        A[0][1] = dpdv[axes[0]];
        A[1][0] = dpdu[axes[1]];
        A[1][1] = dpdv[axes[1]];
        Bx[0] = px[axes[0]] - p[axes[0]];
        Bx[1] = px[axes[1]] - p[axes[1]];
        By[0] = py[axes[0]] - p[axes[0]];
        By[1] = py[axes[1]] - p[axes[1]];
        if (!SolveLinearSystem2x2(A, Bx, &dudx, &dvdx)) {
            dudx = 0.; dvdx = 0.;
        }
        if (!SolveLinearSystem2x2(A, By, &dudy, &dvdy)) {
            dudy = 0.; dvdy = 0.;
        }
    }
    else {
fail:
        dudx = dvdx = 0.;
        dudy = dvdy = 0.;
        dpdx = dpdy = Vector(0,0,0);
    }
}

 

细节

a. 计算 dpdx, dpdy

// Compute auxiliary intersection points with plane
        float d = -Dot(nn, Vector(p.x, p.y, p.z));
        Vector rxv(ray.rxOrigin.x, ray.rxOrigin.y, ray.rxOrigin.z);
        float tx = -(Dot(nn, rxv) + d) / Dot(nn, ray.rxDirection);
        if (isnan(tx)) goto fail;
        Point px = ray.rxOrigin + tx * ray.rxDirection;
        Vector ryv(ray.ryOrigin.x, ray.ryOrigin.y, ray.ryOrigin.z);
        float ty = -(Dot(nn, ryv) + d) / Dot(nn, ray.ryDirection);
        if (isnan(ty)) goto fail;
        Point py = ray.ryOrigin + ty * ray.ryDirection;
        dpdx = px - p;
        dpdy = py - p;

作用：

（计算 dpdx, dpdy，先确定一个平面，这个平面就是 这个相交点 的 位置 P 和法线 组成的，利用 ray-plane 相交的公式就可以计算出 两条 rays rx，ry 的交点 px, py，那么计算 dpdx，dpdy，用到了dx，dy，因为都是偏移一个像素采样点，所以dx，dy 都是1，那么 得到 ）

For this approximation,
we need the plane through the point intersected by the main ray that is tangent to the
surface. This plane is given by the implicit equation

where a = nx, b = ny, c = nz, and d =−(n · p). We can then compute the intersection
points px and py between the auxiliary rays rx and ry with this plane (Figure 10.3). These
new points give an approximation to the partial derivatives of position on the surface
∂p/∂x and ∂p/∂y, based on forward differences:

Because the differential rays are offset one pixel sample in each direction, there’s no need
to divide these differences by a Δ value, since Δ = 1.

The ray–plane intersection algorithm gives the t value where a ray described by origin o
and direction d intersects a plane described by ax + by + cz + d = 0:

To compute this value for the two auxiliary rays, the plane’s d coefficient is computed
first. It isn’t necessary to compute the a, b, and c coefficients, since they’re available in
dg.nn.We can then apply the formula directly.

Figure 10.3: By approximating the local surface geometry at the intersection point with the tangent
plane through p, approximations to the points at which the auxiliary rays rx and ry would intersect
the surface can be found by finding their intersection points with the tangent plane px and py.

 

2. 

float A[2][2], Bx[2], By[2];
        int axes[2];
        if (fabsf(nn.x) > fabsf(nn.y) && fabsf(nn.x) > fabsf(nn.z)) {
            axes[0] = 1; axes[1] = 2;
        }
        else if (fabsf(nn.y) > fabsf(nn.z)) {
            axes[0] = 0; axes[1] = 2;
        }
        else {
            axes[0] = 0; axes[1] = 1;
        }

        // Initialize matrices for chosen projection plane
        A[0][0] = dpdu[axes[0]];
        A[0][1] = dpdv[axes[0]];
        A[1][0] = dpdu[axes[1]];
        A[1][1] = dpdv[axes[1]];
        Bx[0] = px[axes[0]] - p[axes[0]];
        Bx[1] = px[axes[1]] - p[axes[1]];
        By[0] = py[axes[0]] - p[axes[0]];
        By[1] = py[axes[1]] - p[axes[1]];
        if (!SolveLinearSystem2x2(A, Bx, &dudx, &dvdx)) {
            dudx = 0.; dvdx = 0.;
        }
        if (!SolveLinearSystem2x2(A, By, &dudy, &dvdy)) {
            dudy = 0.; dvdy = 0.;
        }

作用：

（计算 dudx,dvdx, dvdy, dvdy,  值得注意的是，最后计算出来的 Δu，Δv 其实就是 dudx，dvdx, .... , 因为上面有说，Δ x 其实就是一个像素采样点的偏移，Δ 为1）

Figure 10.4: An estimate of the difference in (u, v) parametric coordinates from p to p' can be found
by finding the coordinates of p' with respect to the coordinate system defined by p, ∂p/∂u, and ∂p/∂v.

 

Given a position p' on the plane, we can compute its position with respect to the coordinate
system by

or, equivalently,

A solution to this linear system of equations for one of the auxiliary points p' = px
or p' = py gives the corresponding screen space partial derivatives (∂u/∂x, ∂v/∂x) or
(∂u/∂y, ∂v/∂y), respectively

 

This linear system has three equations with two unknowns(3个表达式，但是每一个表达式有2个未知数)—that is, it’s overconstrained（无解）.
We need to be careful since one of the equations may be degenerate—for example, if
∂p/∂u and ∂p/∂v are in the xy plane such that their z components are both zero, then
the third equation will be degenerate. To deal with this case, because we only need two
equations to solve the system, we’d like to choose the two that won’t lead to a degenerate
system. An easy way to do this is to take the cross product of ∂p/∂u and ∂p/∂v, see which
coordinate of the result has the largest magnitude, and use the other two. Their cross
product is already available in nn, so using this approach is straightforward. Even after
all this, it may happen that the linear system has no solution (usually due to the partial
derivatives not forming a coordinate system on the plane). In that case, all that can be
done is to return arbitrary values.

（3个表达式，但是每一个表达式有2个未知数，那么这个 linear system 是无解的，那么 解这个 linear system 的话，其实只需要2个表达式，一种方案就是， 利用 ∂p/∂u and ∂p/∂v 计算出 法线出来，取法线最小的2个轴，假如，法线的x,y 轴分量最小，那么表达式就取 x，y 的表达式）