PBRT笔记(9)——贴图

采样与抗锯齿

当高分辨率贴图被缩小时，贴图会出现严重的混淆现象。虽然第7章中的非均匀采样技术可以减少这种混叠的视觉影响，但是更好的解决方案是实现基于分辨率进行采样的纹理函数。

可以在使用贴图时先对贴图的分辨率进行判断，避免采样高分辨率贴图。
为了解决贴图采样函数造成的混淆问题，我们必须解决以下两个问题：

必须计算贴图空间的采样率，以及获得贴图分辨率，之后就可以计算出屏幕空间的采样率，最后为了获得图元表面的采样率就必须对贴图进行采样，以获得贴图采样率。
对于给定的贴图采样率，必须使用采样理论去引导计算贴图值，不能有高频的贴图。
寻找合适的图片采样率
考虑一个定义在场景图元表面上任意位置的贴图函数，其形参为图元表面的位置信息，T(p)。在忽略可见性与遮挡的情况下，可以表达为这个函数在贴图坐标(x,y)上，即T(f(x,y))。其中将贴图坐标映射到图元表面坐标。

有一个简单的想法，将这个2D贴图函数T(s,t),贴到一个正方体Z轴方向上，并与4个点(0,0,0)、(1,0,0)、(1,1,0)、(0,1,0)对齐。如此一来， $s=P_x \quad t=P_y$ ,它与屏幕坐标的相对关系可以表达为(xr,yr为屏幕图像的分辨率)： $s=\frac{x}{x_r} \quad t=\frac{y}{y_r}$ ，这里就是说明图元贴图UV坐标、屏幕UV坐标以以及映射关系。对于复杂场景、相机位置以及贴图坐标映射来说，获取准确的相对关系比较困难（对应坐标映射到贴图坐标）。对于抗锯齿来说只需要找出像素采样位置的变化与图像上某一点纹理采样位置的变化之间的关系。
这个关系可以使用一个近似的偏导数方程来表示：

$f(x',y') \approx f(x,y)+(x'-x)\frac{\partial_f}{\partial_x}+(y'-y)\frac{\partial_f}{\partial_y}$

如果在x’-x与y’-y方向上的变化率较小，以上的近似公式就是合理的。更重要的是用以上计算出的结果就可以直接算出贴图采样率。例如：在四边形中, $\partial_s/\partial_x=1/x_r \quad \partial_s / \partial_y =0 \quad \partial_t/\partial_x=0 \quad \partial_t /\partial_y=1/y_r$

计算偏导数值的逻辑在RayDifferential类中。使用GenerateRayDifferential()进行初始化，它除了包含了当前位置的光线外，还有x，y方向分别偏移1个像素的两条辅助光线，这两条光线不参与几何图元的相交计算。

之后将计算各种偏导∂p/∂x、∂p/∂y、∂u/∂x、∂v/∂x、∂u/∂y、∂v/∂y。以上值会在ComputeDifferentials()中计算。图10.3很好的展示了计算过程。

void SurfaceInteraction::ComputeDifferentials(
    const RayDifferential &ray) const {
    if (ray.hasDifferentials) {
        //计算微分偏移光线与平面的交点px与py，之后计算微分
        Float d = Dot(n, Vector3f(p.x, p.y, p.z));
        Float tx =
            -(Dot(n, Vector3f(ray.rxOrigin)) - d) / Dot(n, ray.rxDirection);
        if (std::isinf(tx) || std::isnan(tx)) goto fail;
        Point3f px = ray.rxOrigin + tx * ray.rxDirection;
        Float ty =
            -(Dot(n, Vector3f(ray.ryOrigin)) - d) / Dot(n, ray.ryDirection);
        if (std::isinf(ty) || std::isnan(ty)) goto fail;
        Point3f py = ray.ryOrigin + ty * ray.ryDirection;
        dpdx = px - p;
        dpdy = py - p;

        int dim[2];
        if (std::abs(n.x) > std::abs(n.y) && std::abs(n.x) > std::abs(n.z)) {
            dim[0] = 1;
            dim[1] = 2;
        } else if (std::abs(n.y) > std::abs(n.z)) {
            dim[0] = 0;
            dim[1] = 2;
        } else {
            dim[0] = 0;
            dim[1] = 1;
        }

        Float A[2][2] = {{dpdu[dim[0]], dpdv[dim[0]]},
                         {dpdu[dim[1]], dpdv[dim[1]]}};
        Float Bx[2] = {px[dim[0]] - p[dim[0]], px[dim[1]] - p[dim[1]]};
        Float By[2] = {py[dim[0]] - p[dim[0]], py[dim[1]] - p[dim[1]]};
        if (!SolveLinearSystem2x2(A, Bx, &dudx, &dvdx)) dudx = dvdx = 0;
        if (!SolveLinearSystem2x2(A, By, &dudy, &dvdy)) dudy = dvdy = 0;
    } else {
    fail:
        dudx = dvdx = 0;
        dudy = dvdy = 0;
        dpdx = dpdy = Vector3f(0, 0, 0);
    }
}

贴图过滤函数

对于贴图采样率，有必要去除超过纹理采样率Nyquist极限的频率。这时候就需要使用过滤器了，当然本人没《信号与系统》基础，这里直接跳了。

反射与透射的光线微分

为了计算图元表面交点处的反射或透射。我们需要两条经过偏移的微分光线。寻找这些光线需要一些技巧。Igehy(1999)观察到,比起通过x、y方向上的偏移来计算微分光线，如果已知反射方向ωi相对于贴图像素采样上（x、y方向）的变化率，则可以进行以下近似计算：

$\omega=\omega_i+\frac{\partial_{\omega_i}}{\partial_x}$

将反射公式偏导化可得：

$\frac{\partial_{\omega_i}}{\partial_x}=\frac{\partial}{\partial_x}(-\omega_o+2(\omega_o \cdot n)n)=-\frac{\partial_{\omega_o}}{\partial_x}+2((\omega_o \cdot n)\frac{\partial_n}{\partial_x}+\frac{\partial_{(\omega_o \cdot n)}}{\partial_x}n)$

利用点乘性质可得：

$\frac{\partial_{\omega_o \cdot n}}{\partial_x}=\frac{\partial_{\omega_o}}{\partial_x}\cdot n+\omega_o \cdot \frac{\partial_n}{\partial_x}$

以下代码在：SamplerIntegrator::SpecularReflect中

//计算反射微分光线的方向
RayDifferential rd = isect.SpawnRay(wi);
    if (ray.hasDifferentials) {
        rd.hasDifferentials = true;
        rd.rxOrigin = isect.p + isect.dpdx;
        rd.ryOrigin = isect.p + isect.dpdy;
        //计算x、y分量上的方向
        Normal3f dndx = isect.shading.dndu * isect.dudx +
                        isect.shading.dndv * isect.dvdx;
        Normal3f dndy = isect.shading.dndu * isect.dudy +
                        isect.shading.dndv * isect.dvdy;
        Vector3f dwodx = -ray.rxDirection - wo,
                 dwody = -ray.ryDirection - wo;
        Float dDNdx = Dot(dwodx, ns) + Dot(wo, dndx);
        Float dDNdy = Dot(dwody, ns) + Dot(wo, dndy);
        rd.rxDirection =wi - dwodx + 2.f * Vector3f(Dot(wo, ns) * dndx + dDNdx * ns);
        rd.ryDirection =wi - dwody + 2.f * Vector3f(Dot(wo, ns) * dndy + dDNdy * ns);
    }
    return f * Li(rd, scene, sampler, arena, depth + 1) * AbsDot(wi, ns) /pdf;

贴图类与接口

基类为Texture，唯一接口Evaluate(),根据图元表面位置获取对应贴图坐标的值。
子类ConstantTexture，位于constant.h中。返回相同值，因为这个特性，所以他不需要抗锯齿。
子类ScaleTexture，位于scale.h中。返回两个贴图值的乘机。
子类MixTexture，位于mix.h中。使用一个alpha贴图值，返回两个贴图值的线性插值。
子类BilerpTexture，位于bilerp.h中，设定四个贴图坐标，返回以此计算的双线性插值。

ImageTexture

子类ImageTexture，位于imagemap.h中，以二位数组的方式存储贴图函数的采样点（将像素图的每一个像素值当成一个采样点），之后再利用这些采样点还原这个贴图函数。

与上述子类不同，它存储与返回的数据都是参数化的。除了贴图文件名外，还有伽马矫正、mipmap、过滤器（抗锯齿用）等参数。这些都被存储在TexInfo结构体中。

贴图内存管理

PBRT维护了一个map用来管理贴图资源，避免资源被多次加载的情况。

使用GetTexture()来获取贴图的Mipmap

    TexInfo texInfo(filename, doTrilinear, maxAniso, wrap, scale, gamma);
    if (textures.find(texInfo) != textures.end())
        return textures[texInfo].get();
    /*
    中间略
    */
    std::unique_ptr<RGBSpectrum[]> texels = ReadImage(filename, &resolution);
    /*
    中间略
    */
    MIPMap<Tmemory> *mipmap = nullptr;
    if (texels) {
        //readmage返回的是RGBSpectrum类型，所以就必须将其转化为mipmap的类型Tmemory
        std::unique_ptr<Tmemory[]> convertedTexels(
            new Tmemory[resolution.x * resolution.y]);
        for (int i = 0; i < resolution.x * resolution.y; ++i)
            convertIn(texels[i], &convertedTexels[i], scale, gamma);
        mipmap = new MIPMap<Tmemory>(resolution, convertedTexels.get(),
                                     doTrilinear, maxAniso, wrap);
    } else {
        // Create one-valued _MIPMap_
        Tmemory oneVal = scale;
        mipmap = new MIPMap<Tmemory>(Point2i(1, 1), &oneVal);
    }
    textures[texInfo].reset(mipmap);
    return mipmap;

convertIn函数中还涉及到了缩放以及gamma矫正，以便将像素值映射到指定范围。pbrt遵循sRGB标准，该标准规定了一条特定的曲线来匹配CRT显示器的显示。sRGB gamma曲线是一个分段函数，其低值为线性项，其大中型值为幂项。

$\gamma= \left\{ \begin{array}{cc} 12.92x, & x\leq 0.0031308\\ 0, & x >0.0031308 \end{array} \right.$

inline Float GammaCorrect(Float value) {
    if (value <= 0.0031308f) return 12.92f * value;
    return 1.055f * std::pow(value, (Float)(1.f / 2.4f)) - 0.055f;
}

这个函数被用在WriteImage()函数中，用于写入8位sRGB图片数据。当然convertIn中用的是与之相反的操作：

inline Float InverseGammaCorrect(Float value) {
    if (value <= 0.04045f) return value * 1.f / 12.92f;
    return std::pow((value + 0.055f) * 1.f / 1.055f, (Float)2.4f);
}

如果图片渲染完成，会执行ClearCache(),将之前维护资源用的map清空。

mipmap

之前有说高分辨率图像显示在较小的区域中时会出现图像混叠现象，如果使用给定点，以及使用贴图空间采样率进行估算得到结果。这样做消耗比较高。而为了满足Nyquist标准，我们必须去除至少高于相邻采样距离两倍的频率。

贴图采样与重建不同，它对性能有着一定要求，同时其采样率也是不确定的，它会随着场景多边形、贴图坐标映射函数、摄像机投影等因素的变化而变化。

PBRT使用了两种方法实现mipmap,第一种是三线性插值法，快速而且容易实现，在早期显卡中广泛使用。第二种是椭圆加权平均算法，速度慢且复杂，但是效果好。

上述方法均采用了图像金字塔结构（生成分辨率从小到大的金字塔结构），与原始图片相比会多占用1/3空间。使用mipmap得确保图像分辨率为 $2^n$ 。

ImageWrap枚举的作用是：当提供的贴图坐标不在合法的[0,1]范围中时，传递给MIPMap构造函数的指定行为。

如果用于给予的贴图分辨率不是2的幂，PBRT则会讲图片分辨率调整为下一个2的幂。这里涉及到图片缩放涉及到采样与重建理论。我们想要从原始样本中重采样（新的采样位置），从而重建一个连续的图像函数。因为采样率提升了，所以我们不必担心因为采样不足而造成的图像混叠问题，我们只需要重新采样并且重建图像函数。

MipMap使用一个可分离的重构过滤器来完成这个任务。可分离过滤器可以写成以为1维过滤器的乘积f(x,y)=f(x)f(y)。实现重采样可以分为两个一维重采样：第一步：重采样s完成（s’,t）分辨率的图片,第二步重采样t完成（s’,t’）分辨率图片。(s,t)=>(s’,t’) 这样可以大大减少计算复杂度。

resampleWeights()方法确定所有原始像素对新的像素的贡献值权重值。它返回一个ResampleWeight结构体数组。这里的重构器会计算4个原始像素的贡献权重，因为4个像素紧挨在一起，所以只需要一个偏移值和一个权重数组。（以上内容都在构造函数中）

struct ResampleWeight {
    int firstTexel;
    Float weight[4];
};
std::unique_ptr<ResampleWeight[]> resampleWeights(int oldRes, int newRes) {
    CHECK_GE(newRes, oldRes);
    std::unique_ptr<ResampleWeight[]> wt(new ResampleWeight[newRes]);
    Float filterwidth = 2.f;
    for (int i = 0; i < newRes; ++i) {

        Float center = (i + .5f) * oldRes / newRes;
        wt[i].firstTexel = std::floor((center - filterwidth) + 0.5f);
        for (int j = 0; j < 4; ++j) {
            Float pos = wt[i].firstTexel + j + .5f;
            wt[i].weight[j] = Lanczos((pos - center) / filterwidth);
        }
        //规整化操作保证了图像亮度统一
        Float invSumWts = 1 / (wt[i].weight[0] + wt[i].weight[1] +
                               wt[i].weight[2] + wt[i].weight[3]);
        for (int j = 0; j < 4; ++j) wt[i].weight[j] *= invSumWts;
    }
    return wt;
}

获取了权重后，就会根据ImageWrap参数，进行计算，需要计算s与t方向。

由于储存图像使用了大量内存，且每次图像像素查找滤波值都需要读取8~20个像素，为了提升性能PBRT在这里使用了BlockedArray，具体详见附录A。

//存储了mipmap的图像金字塔结构
std::vector<std::unique_ptr<BlockedArray<T>>> pyramid;

第一层为原始图像（如果分辨率不为2的幂，则会存入重采样的图像）。
在展示如何初始化其余级别之前，我们先定义一个texel访问函数：MIPMap::Texel()返回给定离散整数值Texel位置的Texel值的引用。对于超出范围的则会根据wrapMode返回对应的值。

const T &MIPMap<T>::Texel(int level, int s, int t) const {
    CHECK_LT(level, pyramid.size());
    const BlockedArray<T> &l = *pyramid[level];
    switch (wrapMode) {
    case ImageWrap::Repeat:
        s = Mod(s, l.uSize());
        t = Mod(t, l.vSize());
        break;
    case ImageWrap::Clamp:
        s = Clamp(s, 0, l.uSize() - 1);
        t = Clamp(t, 0, l.vSize() - 1);
        break;
    case ImageWrap::Black: {
        static const T black = 0.f;
        if (s < 0 || s >= (int)l.uSize() || t < 0 || t >= (int)l.vSize())
            return black;
        break;
    }
    }
    return l(s, t);
}

最后使用盒式过滤器，生成所有级别的mipmap。这里使用Lanczos过滤器会得到更好的结果。

各项同性三角形过滤器

两个MIPMap::Lookup()方法中第一个方法是返回使用三角形滤波器移除高频信息后的图像，虽然无法生成高质量的结果，但速度比较快。该滤波器因为各项同性的关系不支持非正方形或非轴对称的范围。该滤波器的主要缺点是：在斜角度观察纹理时图像容易变模糊。因为不同的角度会导致采样率不一致。

像素间隔宽度为： $\frac{1}{w}=2^{nLevels-1-l}$

//那么可以求出l为：
Float level = Levels() - 1 + Log2(std::max(width, (Float)1e-8));

if (level < 0)
    return triangle(0, st);
else if (level >= Levels() - 1)
    return Texel(Levels() - 1, 0, 0);
else {
    //通过插值计算来实现不同mipmap级别过度效果
    int iLevel = std::floor(level);
    Float delta = level - iLevel;
    return Lerp(delta, triangle(iLevel, st), triangle(iLevel + 1, st));
}

为了计算图像纹理函数在任意(s,t)位置的值，MIPMAP::triangle寻找(s,t)周围四个最近的像素点，其权重为各自与(s,t)的距离。第一步：使用下面的两个像素点插值得到(s,t0),使用上面的两个像素插值得到(s,t1)。第二步：将上面得到的两个值进行插值计算得到(s,t)。（双线性插值计算）

MIPMap::triangle(int level, const Point2f &st) const {
    level = Clamp(level, 0, Levels() - 1);
    Float s = st[0] * pyramid[level]->uSize() - 0.5f;
    Float t = st[1] * pyramid[level]->vSize() - 0.5f;
    int s0 = std::floor(s), t0 = std::floor(t);
    Float ds = s - s0, dt = t - t0;
    return (1 - ds) * (1 - dt) * Texel(level, s0, t0) +
           (1 - ds) * dt * Texel(level, s0, t0 + 1) +
           ds * (1 - dt) * Texel(level, s0 + 1, t0) +
           ds * dt * Texel(level, s0 + 1, t0 + 1);
}