这个系列是Physically Based Rendering: From Theory To Implementation的读书笔记，本节主要介绍PBRT中的纹理。
现实生活中物体的材质往往不是均匀分布在物体的表面而是在物体表面的不同位置发生变化，因此纹理(texture)对于渲染出具有真实感的图像有至关重要的作用。在PBRT中”纹理”是一个相对广泛的概念，实际上从某个域到另一个域的映射都可以认为是一种纹理。渲染中最常见的纹理是二维图像到三维模型的映射，我们可以在一张二维图像上记录模型不同位置所具有的不同材质，然后渲染时对模型的不同区域加载相应的材质，这样就能够获得更加真实的图像。
Sampling and Antialiasing
在Sampling and Reconstruction一节中我们介绍过纹理是产生走样的重要因素之一。幸运的是纹理所导致的走样是相对容易处理的，实际上只要使用合适的方法我们无需在每个像素进行多次采样就能够渲染出纹理没有走样的图像。
要去除纹理所导致的走样需要考虑2个问题：
- 确定纹理平面中的采样频率。
- 利用采样定律和采样频率来去除纹理中的高频信息。
Finding the Texture Sampling Rate
对于正交相机而言纹理的采样率是相对容易计算的。假设纹理图像位于NDC上由点(0, 0, 0)、(1, 0, 0)、(1, 1, 0)和(0, 1, 0)所定义的平面上，此时纹理坐标与图像坐标的转换关系为：
\[s = \frac{x}{x_r}, \ \ t = \frac{y}{y_r}\]其中$(x_r, y_r)$为图像的分辨率。显然此时纹理图像上的采样间隔为$(\frac{1}{x_r}, \frac{1}{y_r})$，因此任何变化频率大于它的纹理特征都需要去除掉才能保证没有走样。
不难发现纹理采样的频率与纹理平面到NDC上映射的息息相关。对于更加复杂的情况可能很难解析地计算它们之间的映射关系，但我们仍然可以使用泰勒展开来进行近似。对于定义在图像平面上的任意二元函数$f(x, y)$，它的一阶泰勒展开为：
\[f(x', y') \approx f(x, y) + \frac{\partial f}{\partial x} (x' - x) + \frac{\partial f}{\partial y} (y' - y)\]对于前面提到的正交相机的情况，偏导数项分别为$\frac{\partial s}{\partial x} = \frac{1}{x_r}$，$\frac{\partial s}{\partial y} = 0$，$\frac{\partial t}{\partial x} = 0$，$\frac{\partial t}{\partial y} = \frac{1}{y_r}$。显然这些偏导数对于估计纹理采样率起着重要的作用，而在PBRT中则使用光线微分RayDifferential类来存储这些偏导项。RayDifferential类继承自Ray，它包含原始光线和两条辅助光线。这两条辅助光线原点都位于相机中心，方向分别为原始光线在图像坐标上交点的x轴和y轴偏移1个像素。需要注意的是计算光线与场景求交时只考虑原始光线而不考虑两条辅助光线。
class RayDifferential : public Ray {
public:
// RayDifferential Public Data
bool hasDifferentials;
Point3f rxOrigin, ryOrigin;
Vector3f rxDirection, ryDirection;
};
计算这些偏导数的核心是假设光线与模型的交点在局部是一个平面。假设交点为$p$，则$p$点局部对应的平面方程为：
\[a x + b y + c z + d = 0\]其中$a = \mathbf{n}_x$，$b = \mathbf{n}_y$，$c = \mathbf{n}_z$，$d = - (\mathbf{n} \cdot p)$。这样我们就可以利用交点$p$和辅助光线来计算各个偏导数：假设辅助光线和平面的交点分别为$p_x$和$p_y$，则偏导数可以近似为：
\[\frac{\partial p}{\partial x} \approx p_x - p, \ \ \frac{\partial p}{\partial y} \approx p_y - p\]
更进一步，由于x轴和y轴方向的偏移量是1个像素，我们可以忽略偏移量而直接使用偏导数来进行一阶近似。因此这里唯一需要计算的项是辅助光线和平面的交点$p_x$和$p_y$，我们可以利用光线与平面相交的方程来获得。对于平面$a x + b y + c z + d = 0$，光线与平面相交时的时间参数$t$为：
\[t = \frac{-(a, b, c) \cdot o - d}{(a, b, c) \cdot \mathbf{d}}\]因此我们只需要带入$p$点的法向以及微分光线的方向向量解出时间参数$t$，再把$t$带入到微分光线中即可获得交点坐标。
最后来考虑纹理平面上的偏导数。对偏移后的点$p’$同样使用泰勒展开进行一阶近似可以得到：
\[p' = p + \frac{\partial p}{\partial u} \Delta_u + \frac{\partial p}{\partial v} \Delta_v\]上式也可以等价地表示为一个线性方程组：
\[\begin{bmatrix} p'_x - p_x \\ p'_y - p_y \\ p'_z - p_z \\ \end{bmatrix} = \begin{bmatrix} \frac{\partial p_x}{\partial u} & \frac{\partial p_x}{\partial v} \\ \frac{\partial p_y}{\partial u} & \frac{\partial p_y}{\partial v} \\ \frac{\partial p_z}{\partial u} & \frac{\partial p_z}{\partial v} \\ \end{bmatrix} \begin{bmatrix} \Delta_u \\ \Delta_v \end{bmatrix}\]我们使用辅助光线的交点$p_x$和$p_y$来代替$p’$，求解两个方程组即可得到纹理平面关于图像平面的偏导数$(\frac{\partial u}{\partial x}, \frac{\partial v}{\partial x})$和$(\frac{\partial u}{\partial y}, \frac{\partial v}{\partial y})$。
整个计算偏导数的过程封装在SurfaceInteraction::ComputeDifferentials()函数中：
void SurfaceInteraction::ComputeDifferentials(
const RayDifferential &ray) const {
if (ray.hasDifferentials) {
// Estimate screen space change in and
// Compute auxiliary intersection points with plane
Float d = -Dot(n, Vector3f(p.x, p.y, p.z));
Float tx = (-Dot(n, Vector3f(ray.rxOrigin)) - d) /
Dot(n, ray.rxDirection);
Point3f px = ray.rxOrigin + tx * ray.rxDirection;
Float ty = (-Dot(n, Vector3f(ray.ryOrigin)) - d) /
Dot(n, ray.ryDirection);
Point3f py = ray.ryOrigin + ty * ray.ryDirection;
dpdx = px - p;
dpdy = py - p;
// Compute offsets at auxiliary points
// Choose two dimensions to use for ray offset computation
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;
}
// Initialize A, Bx, and By matrices for offset computation
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 {
dudx = dvdx = 0;
dudy = dvdy = 0;
dpdx = dpdy = Vector3f(0, 0, 0);
}
}
Filtering Texture Functions
我们的目标是去除掉纹理图像中的高频成分，因此首先需要利用sinc()函数进行滤波以限制频率：
\[T_b' (x, y) = \int_{-\infty}^\infty \int_{-\infty}^\infty \text{sinc}(x') \text{sinc}(x') \ T'(f(x-x', y-y')) \ dx' \ dy'\]接下来在屏幕空间上利用一个滤波器$g(x, y)$进行卷积来获得最终的纹理：
\[T_f' (x, y) = \int_{-\text{yWidth/2}}^\text{yWidth/2} \int_{-\text{xWidth/2}}^\text{xWidth/2} g(x', y') \ T_b' (x-x', y-y') \ dx' \ dy'\]在实践中有时会完全忽略第二步只利用sinc()函数滤波后的结果来简化渲染流程，甚至是使用box filter来代替sinc()函数。这样做尽管会损失一些图像的质量，但可以极大地简化代码的复杂程度。
除了对纹理图像进行滤波外，如果已知纹理的高频成分还可以对这些频率取平均来减少高频成分的影响。最后一种常用的反走样方法是在采样点附近进行多次采样，并利用这些采样值的平均作为纹理值，这种方法称为超采样(supersampling)。超采样相当于在纹理空间提高了采样率，因此它的计算代价要比滤波的方法高一些。尽管如此，在纹理空间进行超采样仍然比直接提高图像分辨率要高效的多。
Ray Differentials for Specular Reflection and Transmission
光线微分技术在纹理滤波中非常高效，而且我们还可以把光线微分推广到镜面反射以及折射中。以镜面反射为例，微分后的光线起点分别为：
\[p_x = p + \frac{\partial p}{\partial x}, \ p_y = p + \frac{\partial p}{\partial y}\]
接下来考虑微分光线的出射方向。以x方向的微分光线为例，记反射光线的方向为$\omega_i$，对它进行一阶近似有：
\[\omega \approx \omega_i + \frac{\partial \omega_i}{\partial x}\]根据镜面反射入射立体角与出射立体角的关系：
\[\omega_i = -\omega_o + 2 (\omega_o \cdot \mathbf{n}) \mathbf{n}\]对等式两边进行微分得到：
\[\frac{\partial \omega_i}{\partial x} = - \frac{\partial \omega_o}{\partial x} + 2 \bigg( (\omega_o \cdot \mathbf{n}) \frac{\partial \mathbf{n}}{\partial x} + \frac{\partial (\omega_o \cdot \mathbf{n})}{\partial x} \mathbf{n} \bigg)\]其中利用内积的性质有：
\[\frac{\partial (\omega_o \cdot \mathbf{n})}{\partial x} = \frac{\partial \omega_o}{\partial x} \cdot \mathbf{n} + \omega_o \cdot \frac{\partial \mathbf{n}}{\partial x}\]整理上面的式子可以得到反射光线微分的计算代码：
RayDifferential rd = isect.SpawnRay(wi);
if (ray.hasDifferentials) {
rd.hasDifferentials = true;
rd.rxOrigin = isect.p + isect.dpdx;
rd.ryOrigin = isect.p + isect.dpdy;
// Compute differential reflected directions
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;
Vector3f 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);
}
Texture Coordinate Generation
接下来我们讨论纹理坐标的计算方法，即如何把曲面上的点映射到纹理平面上。需要说明的是对于非参数化的曲面使用解析的方法来计算纹理坐标都十分困难，实际上对于复杂曲面进行参数化仍然是图形学尤其是几何处理的难点之一。
PBRT中使用TextureMapping2D类作为各种纹理坐标计算的基类，它通过TextureMapping2D::Map()函数来计算光线与曲面交点处的纹理坐标以及相关的导数。为了便于区分，我们记纹理平面上的坐标为(s, t)而曲面参数化平面上的坐标为(u, v)。TextureMapping2D基类的定义如下：
class TextureMapping2D {
public:
virtual ~TextureMapping2D() { }
virtual Point2f Map(const SurfaceInteraction &si,
Vector2f *dstdx, Vector2f *dstdy) const = 0;
};
2D (u, v) Mapping
UVMapping2D类是最基本的纹理映射，它表示对参数化平面进行平移和缩放。
class UVMapping2D : public TextureMapping2D {
public:
UVMapping2D(Float su = 1, Float sv = 1, Float du = 0, Float dv = 0);
Point2f Map(const SurfaceInteraction &si, Vector2f *dstdx,
Vector2f *dstdy) const;
private:
const Float su, sv, du, dv;
};
UVMapping2D类的纹理坐标计算公式为：
利用链式法则容易得到纹理坐标关于图像坐标的导数：
\[\begin{aligned} \frac{\partial s}{\partial x} &= \frac{\partial u}{\partial x} \frac{\partial s}{\partial u} + \frac{\partial v}{\partial x} \frac{\partial s}{\partial v} \\ &= s_u \frac{\partial u}{\partial x} \end{aligned}\]其它的导数项与它类似。整理后得到UVMapping2D::Map()函数的定义如下：
Point2f UVMapping2D::Map(const SurfaceInteraction &si,
Vector2f *dstdx, Vector2f *dstdy) const {
// Compute texture differentials for 2D (u, v) mapping
*dstdx = Vector2f(su * si.dudx, sv * si.dvdx);
*dstdy = Vector2f(su * si.dudy, sv * si.dvdy);
return Point2f(su * si.uv[0] + du,
sv * si.uv[1] + dv);
}
Spherical Mapping
SphericalMapping2D类表示对球面进行纹理映射。为了适配不同位置不同朝向的曲面，SphericalMapping2D类额外保存了世界坐标系到纹理平面的映射WorldToTexture，同时还定义了成员函数SphericalMapping2D::sphere()来计算参数化平面上的坐标：
class SphericalMapping2D : public TextureMapping2D {
public:
SphericalMapping2D(const Transform &WorldToTexture)
: WorldToTexture(WorldToTexture) { }
Point2f Map(const SurfaceInteraction &si, Vector2f *dstdx,
Vector2f *dstdy) const;
private:
Point2f sphere(const Point3f &p) const;
const Transform WorldToTexture;
};
Point2f SphericalMapping2D::sphere(const Point3f &p) const {
Vector3f vec = Normalize(WorldToTexture(p) - Point3f(0,0,0));
Float theta = SphericalTheta(vec), phi = SphericalPhi(vec);
return Point2f(theta * InvPi, phi * Inv2Pi);
}
在计算导数时，SphericalMapping2D类没有使用解析形式的导数而是使用了数值导数来进行近似：
Point2f SphericalMapping2D::Map(const SurfaceInteraction &si,
Vector2f *dstdx, Vector2f *dstdy) const {
Point2f st = sphere(si.p);
// Compute texture coordinate differentials for sphere mapping
const Float delta = .1f;
Point2f stDeltaX = sphere(si.p + delta * si.dpdx);
*dstdx = (stDeltaX - st) / delta;
Point2f stDeltaY = sphere(si.p + delta * si.dpdy);
*dstdy = (stDeltaY - st) / delta;
// Handle sphere mapping discontinuity for coordinate differentials
if ((*dstdx)[1] > .5) (*dstdx)[1] = 1 - (*dstdx)[1];
else if ((*dstdx)[1] < -.5f) (*dstdx)[1] = -((*dstdx)[1] + 1);
if ((*dstdy)[1] > .5) (*dstdy)[1] = 1 - (*dstdy)[1];
else if ((*dstdy)[1] < -.5f) (*dstdy)[1] = -((*dstdy)[1] + 1);
return st;
}
Cylindrical Mapping
对于柱面的情况，我们使用CylindricalMapping2D类来进行描述。
class CylindricalMapping2D : public TextureMapping2D {
public:
CylindricalMapping2D(const Transform &WorldToTexture)
: WorldToTexture(WorldToTexture) { }
Point2f Map(const SurfaceInteraction &si, Vector2f *dstdx,
Vector2f *dstdy) const;
private:
Point2f cylinder(const Point3f &p) const {
Vector3f vec = Normalize(WorldToTexture(p) - Point3f(0,0,0));
return Point2f((Pi + std::atan2(vec.y, vec.x)) * Inv2Pi, vec.z);
}
const Transform WorldToTexture;
};
类似于SphericalMapping2D这里同样使用了数值导数：
Point2f CylindricalMapping2D::Map(const SurfaceInteraction &si,
Vector2f *dstdx, Vector2f *dstdy) const {
Point2f st = cylinder(si.p);
// Compute texture coordinate differentials for cylinder mapping
const Float delta = .01f;
Point2f stDeltaX = cylinder(si.p + delta * si.dpdx);
*dstdx = (stDeltaX - st) / delta;
if ((*dstdx)[1] > .5) (*dstdx)[1] = 1.f - (*dstdx)[1];
else if ((*dstdx)[1] < -.5f) (*dstdx)[1] = -((*dstdx)[1] + 1);
Point2f stDeltaY = cylinder(si.p + delta * si.dpdy);
*dstdy = (stDeltaY - st) / delta;
if ((*dstdy)[1] > .5) (*dstdy)[1] = 1.f - (*dstdy)[1];
else if ((*dstdy)[1] < -.5f) (*dstdy)[1] = -((*dstdy)[1] + 1);
return st;
}
Planar Mapping
最后一种常用的纹理映射是二维平面到二维平面的映射，我们使用PlanarMapping2D类来进行表示。二维平面到二维平面的映射可以表示为一个仿射变换，因此这里使用两个向量$\mathbf{v}_s$和$\mathbf{v}_t$以及偏移量$d_s$和$d_t$来描述整个仿射变换。
class PlanarMapping2D : public TextureMapping2D {
public:
Point2f Map(const SurfaceInteraction &si, Vector2f *dstdx,
Vector2f *dstdy) const;
PlanarMapping2D(const Vector3f &vs, const Vector3f &vt,
Float ds = 0, Float dt = 0)
: vs(vs), vt(vt), ds(ds), dt(dt) { }
private:
const Vector3f vs, vt;
const Float ds, dt;
};
由于仿射变换是一个线性变换，我们可以使用解析的方式来计算纹理坐标以及各个导数：
Point2f PlanarMapping2D::Map(const SurfaceInteraction &si,
Vector2f *dstdx, Vector2f *dstdy) const {
Vector3f vec(si.p);
*dstdx = Vector2f(Dot(si.dpdx, vs), Dot(si.dpdx, vt));
*dstdy = Vector2f(Dot(si.dpdy, vs), Dot(si.dpdy, vt));
return Point2f(ds + Dot(vec, vs), dt + Dot(vec, vt));
}
3D Mapping
本节最后引入了三维纹理映射的计算方法。从广义纹理的概念来讲，三维纹理仅仅是把二维纹理坐标推广到了三维空间上，因此对于基类只需要修改Map()函数即可。
class TextureMapping3D {
public:
virtual ~TextureMapping3D() { }
virtual Point3f Map(const SurfaceInteraction &si,
Vector3f *dpdx, Vector3f *dpdy) const = 0;
};
最简单的三维纹理映射是使用一个线性变换将三维空间中的点映射到纹理空间，我们使用TransformMapping3D来描述这种映射过程。
class TransformMapping3D : public TextureMapping3D {
public:
Point3f TransformMapping3D::Map(const SurfaceInteraction &si,
Vector3f *dpdx, Vector3f *dpdy) const {
*dpdx = WorldToTexture(si.dpdx);
*dpdy = WorldToTexture(si.dpdy);
return WorldToTexture(si.p);
}
private:
const Transform WorldToTexture;
};
Texture Interface and Basic Textures
PBRT中纹理类都继承自基类Texture，通过Texture::Evaluate()函数来获得纹理值。
template <typename T>
class Texture {
public:
virtual T Evaluate(const SurfaceInteraction &) const = 0;
virtual ~Texture() { }
};
Constant Texture
最简单的纹理类是常量纹理ConstantTexture，它直接返回常量value。
template <typename T>
class ConstantTexture : public Texture<T> {
public:
ConstantTexture(const T &value) : value(value) { }
T Evaluate(const SurfaceInteraction &) const {
return value;
}
private:
T value;
};
Scale Texture
ScaleTexture类表示将两个纹理相乘的值作为新的纹理。
template <typename T1, typename T2>
class ScaleTexture : public Texture<T2> {
public:
ScaleTexture(const std::shared_ptr<Texture<T1>> &tex1,
const std::shared_ptr<Texture<T2>> &tex2)
: tex1(tex1), tex2(tex2) { }
T2 Evaluate(const SurfaceInteraction &si) const {
return tex1->Evaluate(si) * tex2->Evaluate(si);
}
private:
std::shared_ptr<Texture<T1>> tex1;
std::shared_ptr<Texture<T2>> tex2;
};
Mix Textures
MixTexture是对ScaleTexture的推广，它通过一个额外的纹理amount实现对两个给定纹理tex1和tex2的线性插值。
template <typename T>
class MixTexture : public Texture<T> {
public:
MixTexture(const std::shared_ptr<Texture<T>> &tex1,
const std::shared_ptr<Texture<T>> &tex2,
const std::shared_ptr<Texture<Float>> &amount)
: tex1(tex1), tex2(tex2), amount(amount) { }
T Evaluate(const SurfaceInteraction &si) const {
T t1 = tex1->Evaluate(si), t2 = tex2->Evaluate(si);
Float amt = amount->Evaluate(si);
return (1 - amt) * t1 + amt * t2;
}
private:
std::shared_ptr<Texture<T>> tex1, tex2;
std::shared_ptr<Texture<Float>> amount;
};
Bilinear Interpolation
在计算纹理时需要使用双线性插值来提取非整数坐标上的纹理值，在这种情况下可以使用BilerpTexture类。BilerpTexture类接收一个TextureMapping2D对象来计算纹理坐标(s, t)，(s, t)处的纹理值由(0, 0)、(0, 1)、(1, 0)、(1, 1)四个点通过插值来进行计算。
template <typename T>
class BilerpTexture : public Texture<T> {
public:
BilerpTexture(std::unique_ptr<TextureMapping2D> mapping, const T &v00,
const T &v01, const T &v10, const T &v11)
: mapping(std::move(mapping)), v00(v00), v01(v01), v10(v10),
v11(v11) { }
T Evaluate(const SurfaceInteraction &si) const {
Vector2f dstdx, dstdy;
Point2f st = mapping->Map(si, &dstdx, &dstdy);
return (1-st[0]) * (1-st[1]) * v00 + (1-st[0]) * (st[1]) * v01 +
( st[0]) * (1-st[1]) * v10 + ( st[0]) * (st[1]) * v11;
}
private:
std::unique_ptr<TextureMapping2D> mapping;
const T v00, v01, v10, v11;
};
Image Texture
上面介绍的纹理计算方法只适用于比较简单的纹理，对于更加复杂的情况我们往往会使用一张图片来记录纹理，然后在渲染时加载它来获得曲面上的纹理值，这样的纹理表示方法称为纹理映射(texture map)。在PBRT中使用ImageTexture类来表示纹理映射，它同样继承自Texture基类但比前面介绍过的派生类要复杂一些。
template <typename Tmemory, typename Treturn>
class ImageTexture : public Texture<Treturn> {
public:
// ImageTexture Public Methods
ImageTexture(std::unique_ptr<TextureMapping2D> m, const std::string &filename,
bool doTri, Float maxAniso, ImageWrap wm, Float scale, bool gamma);
static void ClearCache();
Treturn Evaluate(const SurfaceInteraction &si) const;
private:
// ImageTexture Private Methods
static MIPMap<Tmemory> *GetTexture(const std::string &filename,
bool doTrilinear, Float maxAniso, ImageWrap wm, Float scale, bool gamma);
static void convertIn(const RGBSpectrum &from, RGBSpectrum *to, Float scale, bool gamma);
static void convertIn(const RGBSpectrum &from, Float *to, Float scale, bool gamma);
static void convertOut(const RGBSpectrum &from, Spectrum *to);
static void convertOut(Float from, Float *to);
// ImageTexture Private Data
std::unique_ptr<TextureMapping2D> mapping;
MIPMap<Tmemory> *mipmap;
static std::map<TexInfo, std::unique_ptr<MIPMap<Tmemory>>> textures;
};
在初始化ImageTexture时会将纹理图像加载到内存中，同时初始化一个MIPMap对象。
template <typename Tmemory, typename Treturn>
ImageTexture<Tmemory, Treturn>::ImageTexture(
std::unique_ptr<TextureMapping2D> mapping,
const std::string &filename, bool doTrilinear, Float maxAniso,
ImageWrap wrapMode, Float scale, bool gamma)
: mapping(std::move(mapping)) {
mipmap = GetTexture(filename, doTrilinear, maxAniso,
wrapMode, scale, gamma);
}
Texture Memory Management
在渲染时同一张纹理图可能会被多次访问到。为了提高内存的访问效率，ImageTexture通过静态成员textures来保存所有的纹理图像。在调用ImageTexture::GetTexture()函数访问纹理时，如果已经在内存中加载了纹理则返回纹理的MIPMap，否则新建一个MIPMap并放入内存中。
template <typename Tmemory, typename Treturn>
MIPMap<Tmemory> *ImageTexture<Tmemory, Treturn>::GetTexture(
const std::string &filename,
bool doTrilinear, Float maxAniso, ImageWrap wrap, Float scale,
bool gamma) {
// Return MIPMap from texture cache if present
TexInfo texInfo(filename, doTrilinear, maxAniso, wrap, scale, gamma);
if (textures.find(texInfo) != textures.end())
return textures[texInfo].get();
// Create MIPMap for filename
Point2i resolution;
std::unique_ptr<RGBSpectrum[]> texels = ReadImage(filename, &resolution);
MIPMap<Tmemory> *mipmap = nullptr;
if (texels) {
// Convert texels to type Tmemory and create MIPMap
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;
}
其中TexInfo是一个简单的结构体，它保存了纹理的文件名已经相关的滤波信息。当需要读取图像到内存中时会调用ReadImage()函数将图像数据加载到texels中并利用convertIn()函数对图像进行矫正。
static void convertIn(const RGBSpectrum &from, RGBSpectrum *to,
Float scale, bool gamma) {
for (int i = 0; i < RGBSpectrum::nSamples; ++i)
(*to)[i] = scale * (gamma ? InverseGammaCorrect(from[i])
: from[i]);
}
static void convertIn(const RGBSpectrum &from, Float *to,
Float scale, bool gamma) {
*to = scale * (gamma ? InverseGammaCorrect(from.y())
: from.y());
}
这里需要补充一些Gamma矫正的内容。PBRT使用了sRGB标准来描述显示颜色与像素值之间的对应关系，对于[0, 1]范围内的像素值$x$经过Gamma矫正后的值$\gamma(x)$才是显示设备实际显示的值：
\[\gamma(x) = \left\{ \begin{aligned} & 12.92 \ x, & & x \leq 0.0031308 \\ & 1.055 \ x^{1/2.4}, & & x \gt 0.0031308 \end{aligned} \right.\]这样我们可以定义Gamma矫正和反矫正函数：
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;
}
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);
}
在加载图像时需要使用InverseGammaCorrect()恢复像素的原始值，而在输出图像时则需要使用GammaCorrect()将像素值转换成sRGB标准输出。
最后，当渲染工作结束后需要调用ImageTexture::ClearCache()来清空缓存。
static void ClearCache() {
textures.erase(textures.begin(), textures.end());
}
ImageTexture Evaluation
ImageTexture类同样使用Evaluate()函数来计算纹理值。计算时会利用MIPMap进行反走样滤波，然后再通过convertOut()函数返回纹理值。
Treturn Evaluate(const SurfaceInteraction &si) const {
Vector2f dstdx, dstdy;
Point2f st = mapping->Map(si, &dstdx, &dstdy);
Tmemory mem = mipmap->Lookup(st, dstdx, dstdy);
Treturn ret;
convertOut(mem, &ret);
return ret;
}
static void convertOut(const RGBSpectrum &from, Spectrum *to) {
Float rgb[3];
from.ToRGB(rgb);
*to = Spectrum::FromRGB(rgb);
}
static void convertOut(Float from, Float *to) {
*to = from;
}
MIP Maps
对于纹理图像我们需要使用一些反走样技术才能获得高质量的渲染结果。和图像滤波的反走样类似，对于超过Nyquist极限的频率成分我们需要将它们提前滤掉。和在渲染中对图像滤波不同的是，对纹理进行滤波是相对容易的：我们无需发射光线对场景进行渲染，同时纹理图像上的频率是确定的，我们也无需考虑频率无穷大的情况。因此我们只需要对纹理图像进行一些预处理即可完成反走样滤波。
纹理图像滤波的难点在于渲染时每个像素对应纹理图像上的采样频率是不一样的。每个像素的采样频率由模型的几何、纹理坐标映射、相机投影、图像采样频率等很多因素共同决定，因此我们无法使用固定的采样频率而需要使用一些方法来动态计算图像不同像素的采样频率。
在PBRT中使用MIPMap类来对图像上的不同像素进行反走样，它的定义如下：
enum class ImageWrap { Repeat, Black, Clamp };
template <typename T>
class MIPMap {
public:
// MIPMap Public Methods
MIPMap(const Point2i &resolution, const T *data, bool doTri = false,
Float maxAniso = 8.f, ImageWrap wrapMode = ImageWrap::Repeat);
int Width() const { return resolution[0]; }
int Height() const { return resolution[1]; }
int Levels() const { return pyramid.size(); }
const T &Texel(int level, int s, int t) const;
T Lookup(const Point2f &st, Float width = 0.f) const;
T Lookup(const Point2f &st, Vector2f dstdx, Vector2f dstdy) const;
private:
// MIPMap Private Methods
std::unique_ptr<ResampleWeight[]> resampleWeights(int oldRes, int newRes);
Float clamp(Float v) { return Clamp(v, 0.f, Infinity); }
RGBSpectrum clamp(const RGBSpectrum &v) { return v.Clamp(0.f, Infinity); }
SampledSpectrum clamp(const SampledSpectrum &v) { return v.Clamp(0.f, Infinity); }
T triangle(int level, const Point2f &st) const;
T EWA(int level, Point2f st, Vector2f dst0, Vector2f dst1) const;
// MIPMap Private Data
const bool doTrilinear;
const Float maxAnisotropy;
const ImageWrap wrapMode;
Point2i resolution;
std::vector<std::unique_ptr<BlockedArray<T>>> pyramid;
static constexpr int WeightLUTSize = 128;
static Float weightLut[WeightLUTSize];
};
简单来说，MIPMap使用了图像金字塔来表示纹理图像。其中原始纹理图像位于最底层，它上面每一层图像的高和宽都缩减到下一层图像的$\frac{1}{2}$，而最上层只包含1个像素。这样在计算纹理值时只需要选择相应的层上对应纹理坐标出的值即可。在初始化时MIPMap会显式构造出每一层上的纹理图像，具体的构造函数如下：
struct ResampleWeight {
int firstTexel;
Float weight[4];
};
template <typename T>
MIPMap<T>::MIPMap(const Point2i &res, const T *img, bool doTrilinear,
Float maxAnisotropy, ImageWrap wrapMode)
: doTrilinear(doTrilinear), maxAnisotropy(maxAnisotropy),
wrapMode(wrapMode), resolution(res) {
std::unique_ptr<T[]> resampledImage = nullptr;
if (!IsPowerOf2(resolution[0]) || !IsPowerOf2(resolution[1])) {
// Resample image to power-of-two resolution
Point2i resPow2(RoundUpPow2(resolution[0]), RoundUpPow2(resolution[1]));
// Resample image in s direction
std::unique_ptr<ResampleWeight[]> sWeights = resampleWeights(resolution[0], resPow2[0]);
resampledImage.reset(new T[resPow2[0] * resPow2[1]]);
// Apply sWeights to zoom in s direction
ParallelFor(
[&](int t) {
for (int s = 0; s < resPow2[0]; ++s) {
// Compute texel (s, t) in s-zoomed image
resampledImage[t * resPow2[0] + s] = 0.f;
for (int j = 0; j < 4; ++j) {
int origS = sWeights[s].firstTexel + j;
if (wrapMode == ImageWrap::Repeat)
origS = Mod(origS, resolution[0]);
else if (wrapMode == ImageWrap::Clamp)
origS = Clamp(origS, 0, resolution[0] - 1);
if (origS >= 0 && origS < (int)resolution[0])
resampledImage[t * resPow2[0] + s] +=
sWeights[s].weight[j] * img[t * resolution[0] + origS];
}
}
}, resolution[1], 16);
// Resample image in t direction
std::unique_ptr<ResampleWeight[]> tWeights = resampleWeights(resolution[1], resPow2[1]);
std::vector<T *> resampleBufs;
int nThreads = MaxThreadIndex();
for (int i = 0; i < nThreads; ++i)
resampleBufs.push_back(new T[resPow2[1]]);
ParallelFor(
[&](int s) {
T *workData = resampleBufs[threadIndex];
for (int t = 0; t < resPow2[1]; ++t) {
workData[t] = 0.f;
for (int j = 0; j < 4; ++j) {
int offset = tWeights[t].firstTexel + j;
if (wrapMode == ImageWrap::Repeat) offset = Mod(offset, resolution[1]);
else if (wrapMode == ImageWrap::Clamp) offset = Clamp(offset, 0, (int)resolution[1]-1);
if (offset >= 0 && offset < (int)resolution[1])
workData[t] += tWeights[t].weight[j] *
resampledImage[offset * resPow2[0] + s];
}
}
for (int t = 0; t < resPow2[1]; ++t)
resampledImage[t * resPow2[0] + s] = clamp(workData[t]);
}, resPow2[0], 32);
for (auto ptr : resampleBufs)
delete[] ptr;
resolution = resPow2;
}
// Initialize levels of MIPMap from image
int nLevels = 1 + Log2Int(std::max(resolution[0], resolution[1]));
pyramid.resize(nLevels);
// Initialize most detailed level of MIPMap
pyramid[0].reset(new BlockedArray<T>(resolution[0], resolution[1],
resampledImage ? resampledImage.get() : img));
for (int i = 1; i < nLevels; ++i) {
// Initialize i MIPMap level from i-1 level
int sRes = std::max(1, pyramid[i - 1]->uSize() / 2);
int tRes = std::max(1, pyramid[i - 1]->vSize() / 2);
pyramid[i].reset(new BlockedArray<T>(sRes, tRes));
// Filter four texels from finer level of pyramid
ParallelFor(
[&](int t) {
for (int s = 0; s < sRes; ++s)
(*pyramid[i])(s, t) = .25f *
(Texel(i-1, 2*s, 2*t) + Texel(i-1, 2*s+1, 2*t) +
Texel(i-1, 2*s, 2*t+1) + Texel(i-1, 2*s+1, 2*t+1));
}, tRes, 16);
}
// Initialize EWA filter weights if needed
if (weightLut[0] == 0.) {
for (int i = 0; i < WeightLUTSize; ++i) {
Float alpha = 2;
Float r2 = Float(i) / Float(WeightLUTSize - 1);
weightLut[i] = std::exp(-alpha * r2) - std::exp(-alpha);
}
}
}
如果原始纹理图像的高宽为2的幂则直接通过上采样构造出上一层的纹理，其中每个像素为下一层对应$2 \times 2$区域像素的平均值；否则需要首先把原始纹理图像的尺寸拓展到2的幂再构造金字塔。在拓展图像尺寸时使用了重采样的方法，我们需要在s和t两个方向上分别进行重采样来保证采样后纹理图像的尺寸均为2的幂。
在具体进行重采样时每次选择采样点最近的4个点进行加权平均来作为采样后的值，这四个点的权重由MIPMap::resampleWeights()函数来计算：
std::unique_ptr<ResampleWeight[]> resampleWeights(int oldRes, int newRes) {
Assert(newRes >= oldRes);
std::unique_ptr<ResampleWeight[]> wt(new ResampleWeight[newRes]);
Float filterwidth = 2.f;
for (int i = 0; i < newRes; ++i) {
// Compute image resampling weights for ith texel
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);
}
// Normalize filter weights for texel resampling
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;
}
在访问MIPMap对象时可以通过MIPMap::Texel()函数来获取指定层数l上(s, t)位置处的纹理值：
template <typename T>
const T &MIPMap<T>::Texel(int level, int s, int t) const {
const BlockedArray<T> &l = *pyramid[level];
// Compute texel accounting for boundary conditions
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);
}
Isotropic Triangle Filter
在ImageTexture::Evaluate()函数中通过MIPMap::Lookup()来获得纹理值。PBRT提供了两个Lookup()函数，其中比较简单的是isotropic triangle filter，它的基本思想是选择某一层的纹理图像使得采样点(s, t)对应的窗口恰好包含4个纹理样本点，然后再通过插值来计算出(s, t)处的纹理值。triangle filter要求纹理图像必须是一个正方形图像，同时滤波器的窗口也是一个正方形，因此triangle filter是一个各向同性(isotropic)滤波器。
纹理图像的层数l由窗口大小w决定，它们之间的关系为：
\[\frac{1}{w} = 2^{\text{nLevels}-1-l}\]这样计算出层数l后即可在对应的纹理图像上进行插值。在大多数情况下这样解出的层数l都不是整数，因此还需要在两个相邻层上进行插值，如下图所示。
整理之后可以得到isotropic triangle filter对应的Lookup()函数：
template <typename T>
T MIPMap<T>::Lookup(const Point2f &st, Float width) const {
// Compute MIPMap level for trilinear filtering
Float level = Levels() - 1 + Log2(std::max(width, (Float)1e-8));
// Perform trilinear interpolation at appropriate MIPMap level
if (level < 0)
return triangle(0, st);
else if (level >= Levels() - 1)
return Texel(Levels() - 1, 0, 0);
else {
int iLevel = std::floor(level);
Float delta = level - iLevel;
return Lerp(delta, triangle(iLevel, st), triangle(iLevel + 1, st));
}
}
每一层上具体的插值计算由MIPMap::triangle()函数来实现，它表示对图像使用triangle filter，每个像素的权重为：
template <typename T>
T MIPMap<T>::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);
}
Elliptically Weighted Average
另一种比较复杂的方法是EWA算法(elliptically weighted average)。EWA算法会考虑纹理坐标的微分在纹理图像上的两个主方向，然后通过Gaussian核对纹理进行滤波。和triangle filter不同的是，EWA算法可以考虑任意方向上纹理坐标的变化，因此它是一个各向异性(anisotropic)的滤波器。因此EWA算法可以极大地提高纹理滤波的质量，它也是目前最好的滤波算法之一。
EWA算法通过MIPMap::Lookup()函数来进行调用，它的定义如下：
template <typename T>
T MIPMap<T>::Lookup(const Point2f &st, Vector2f dst0,
Vector2f dst1) const {
if (doTrilinear) {
Float width = std::max(std::max(std::abs(dst0[0]),
std::abs(dst0[1])),
std::max(std::abs(dst1[0]),
std::abs(dst1[1])));
return Lookup(st, 2 * width);
}
// Compute ellipse minor and major axes
if (dst0.LengthSquared() < dst1.LengthSquared())
std::swap(dst0, dst1);
Float majorLength = dst0.Length();
Float minorLength = dst1.Length();
// Clamp ellipse eccentricity if too large
if (minorLength * maxAnisotropy < majorLength && minorLength > 0) {
Float scale = majorLength / (minorLength * maxAnisotropy);
dst1 *= scale;
minorLength *= scale;
}
if (minorLength == 0)
return triangle(0, st);
// Choose level of detail for EWA lookup and perform EWA filtering
Float lod = std::max((Float)0, Levels() - (Float)1 + Log2(minorLength));
int ilod = std::floor(lod);
return Lerp(lod - ilod, EWA(ilod, st, dst0, dst1),
EWA(ilod + 1, st, dst0, dst1));
}
EWA算法同样需要在金字塔的不同层上进行插值。计算层号的方法与triangle filter相同，但在每一层上则需要调用MIPMap::EWA()函数来计算滤波后的纹理值。
template <typename T>
T MIPMap<T>::EWA(int level, Point2f st, Vector2f dst0,
Vector2f dst1) const {
if (level >= Levels()) return Texel(Levels() - 1, 0, 0);
// Convert EWA coordinates to appropriate scale for level
st[0] = st[0] * pyramid[level]->uSize() - 0.5f;
st[1] = st[1] * pyramid[level]->vSize() - 0.5f;
dst0[0] *= pyramid[level]->uSize();
dst0[1] *= pyramid[level]->vSize();
dst1[0] *= pyramid[level]->uSize();
dst1[1] *= pyramid[level]->vSize();
// Compute ellipse coefficients to bound EWA filter region
Float A = dst0[1] * dst0[1] + dst1[1] * dst1[1] + 1;
Float B = -2 * (dst0[0] * dst0[1] + dst1[0] * dst1[1]);
Float C = dst0[0] * dst0[0] + dst1[0] * dst1[0] + 1;
Float invF = 1 / (A * C - B * B * 0.25f);
A *= invF;
B *= invF;
C *= invF;
// Compute the ellipse’s bounding box in texture space
Float det = -B * B + 4 * A * C;
Float invDet = 1 / det;
Float uSqrt = std::sqrt(det * C), vSqrt = std::sqrt(A * det);
int s0 = std::ceil (st[0] - 2 * invDet * uSqrt);
int s1 = std::floor(st[0] + 2 * invDet * uSqrt);
int t0 = std::ceil (st[1] - 2 * invDet * vSqrt);
int t1 = std::floor(st[1] + 2 * invDet * vSqrt);
// Scan over ellipse bound and compute quadratic equation
T sum(0.f);
Float sumWts = 0;
for (int it = t0; it <= t1; ++it) {
Float tt = it - st[1];
for (int is = s0; is <= s1; ++is) {
Float ss = is - st[0];
// Compute squared radius and filter texel if inside ellipse
Float r2 = A * ss * ss + B * ss * tt + C * tt * tt;
if (r2 < 1) {
int index = std::min((int)(r2 * WeightLUTSize),
WeightLUTSize - 1);
Float weight = weightLut[index];
sum += Texel(level, is, it) * weight;
sumWts += weight;
}
}
}
return sum / sumWts;
}
MIPMap::EWA()的核心是对纹理图像上椭圆区域内的纹理值进行加权平均。因此首先需要把规范化坐标映射到纹理图像上确定椭圆窗口的范围，然后利用椭圆方程判断像素是否位于滤波器窗口中。对于窗口中的像素，每个像素的权重即为标准高斯核在(s’-s, t’-t)处的值。
Solid and Procedural Texturing
这一节中我们来介绍三维纹理。前面介绍的纹理映射都是二维到二维平面的映射，我们可以把二维纹理平面推广到三维这样就得到了三维空间中的纹理映射。实际上三维纹理映射并没有那么的抽象，如果把物体的空间坐标视为纹理空间，那么物体的空间位置就构成了天然的三维纹理映射。
三维纹理的一大难点在于如何进行存储。由于增加了一个维度，直接在内存中存储三维纹理会极大地提高内存需求；同时，这种显式的纹理存储也不利于设计师去设计实际的纹理。因此三维纹理一般是通过程序化(procedural)的形式来进行存储，即我们不会记录纹理的值而是根据需要来生成纹理。这种形式的纹理有很多好处，比如说它节省了大量的内存空间同时我们可以获得任意精度的纹理形状而无需考虑纹理图像分辨率。程序纹理当然也存在一些缺陷，比方说很难直接对这样的纹理进行反走样滤波，同时程序纹理也需要更多的计算资源。
UV Texture
UV纹理是一种简单的程序纹理，我们直接把曲面的(u, v)坐标作为纹理颜色的R和G通道值即可。这种纹理在调试曲面的参数化时非常好用。
class UVTexture : public Texture<Spectrum> {
public:
UVTexture(std::unique_ptr<TextureMapping2D> mapping) : mapping(std::move(mapping)) { }
Spectrum Evaluate(const SurfaceInteraction &si) const {
Vector2f dstdx, dstdy;
Point2f st = mapping->Map(si, &dstdx, &dstdy);
Float rgb[3] =
{ st[0] - std::floor(st[0]), st[1] - std::floor(st[1]), 0 };
return Spectrum::FromRGB(rgb);
}
private:
std::unique_ptr<TextureMapping2D> mapping;
};
Checkerboard
棋盘格纹理是一种常见的纹理形式，它的特点是曲面上会出现颜色交替的格子。在PBRT的实现中格子颜色不仅限于黑白两种，实际上可以使用任意两种纹理来构造棋盘格。
计算曲面上任意点的纹理值非常简单，只需要根据曲面(u, v)坐标对两个纹理进行轮换即可。
// Point sample Checkerboard2DTexture
if (((int)std::floor(st[0]) + (int)std::floor(st[1])) % 2 == 0)
return tex1->Evaluate(si);
return tex2->Evaluate(si);
对棋盘格纹理进行反走样滤波会复杂一些。我们首先需要计算纹理平面关于图像平面的导数，然后在纹理平面上使用一个矩形框把对应的区域框起来。如果这个框都在某个纹理中则直接返回该纹理上的纹理值，否则还需要进行滤波。
这里我们使用一个指示函数来指示从哪个纹理上计算纹理值：
\[c(x) = \left\{ \begin{aligned} & 0, & & \lfloor x \rfloor \ \text{is even} \\ & 1, & & \text{otherwise} \end{aligned} \right.\]$c(x)$积分的平均值表示使用1号纹理的比例：
\[\int_0^x c(x) \ dx = \lfloor x/2 \rfloor + 2 \max (x/2 - \lfloor x/2 \rfloor - 0.5, 0)\]我们再对两个纹理进行加权平均就完成了滤波：
auto bumpInt = [](Float x) {
return (int)std::floor(x / 2) +
2 * std::max(x / 2 - (int)std::floor(x / 2) - (Float)0.5, (Float)0); };
Float sint = (bumpInt(s1) - bumpInt(s0)) / (2 * ds);
Float tint = (bumpInt(t1) - bumpInt(t0)) / (2 * dt);
Float area2 = sint + tint - 2 * sint * tint;
if (ds > 1 || dt > 1) area2 = .5f;
return (1 - area2) * tex1->Evaluate(si) +
area2 * tex2->Evaluate(si);
完整的棋盘格纹理定义如下：
enum class AAMethod { None, ClosedForm };
template <typename T>
class Checkerboard2DTexture : public Texture<T> {
public:
Checkerboard2DTexture(std::unique_ptr<TextureMapping2D> mapping,
const std::shared_ptr<Texture<T>> &tex1,
const std::shared_ptr<Texture<T>> &tex2, AAMethod aaMethod)
: mapping(std::move(mapping)), tex1(tex1), tex2(tex2), aaMethod(aaMethod) { }
T Evaluate(const SurfaceInteraction &si) const {
Vector2f dstdx, dstdy;
Point2f st = mapping->Map(si, &dstdx, &dstdy);
if (aaMethod == AAMethod::None) {
// Point sample Checkerboard2DTexture
if (((int)std::floor(st[0]) + (int)std::floor(st[1])) % 2 == 0)
return tex1->Evaluate(si);
return tex2->Evaluate(si);
} else {
// Compute closed-form box-filtered Checkerboard2DTexture value
// Evaluate single check if filter is entirely inside one of them
Float ds = std::max(std::abs(dstdx[0]), std::abs(dstdy[0]));
Float dt = std::max(std::abs(dstdx[1]), std::abs(dstdy[1]));
Float s0 = st[0] - ds, s1 = st[0] + ds;
Float t0 = st[1] - dt, t1 = st[1] + dt;
if (std::floor(s0) == std::floor(s1) &&
std::floor(t0) == std::floor(t1)) {
// Point sample Checkerboard2DTexture>
if (((int)std::floor(st[0]) + (int)std::floor(st[1])) % 2 == 0)
return tex1->Evaluate(si);
return tex2->Evaluate(si);
}
// Apply box filter to checkerboard region
auto bumpInt = [](Float x) {
return (int)std::floor(x / 2) +
2 * std::max(x / 2 - (int)std::floor(x / 2) - (Float)0.5, (Float)0); };
Float sint = (bumpInt(s1) - bumpInt(s0)) / (2 * ds);
Float tint = (bumpInt(t1) - bumpInt(t0)) / (2 * dt);
Float area2 = sint + tint - 2 * sint * tint;
if (ds > 1 || dt > 1) area2 = .5f;
return (1 - area2) * tex1->Evaluate(si) + area2 * tex2->Evaluate(si);
}
}
private:
std::unique_ptr<TextureMapping2D> mapping;
const std::shared_ptr<Texture<T>> tex1, tex2;
const AAMethod aaMethod;
};
Solid Checkerboard
Checkerboard3DTexture类定义了三维棋盘格纹理。
template <typename T>
class Checkerboard3DTexture : public Texture<T> {
public:
Checkerboard3DTexture(std::unique_ptr<TextureMapping3D> mapping,
const std::shared_ptr<Texture<T>> &tex1,
const std::shared_ptr<Texture<T>> &tex2)
: mapping(std::move(mapping)), tex1(tex1), tex2(tex2) { }
T Evaluate(const SurfaceInteraction &si) const {
Vector3f dpdx, dpdy;
Point3f p = mapping->Map(si, &dpdx, &dpdy);
if (((int)std::floor(p.x) + (int)std::floor(p.y) +
(int)std::floor(p.z)) % 2 == 0)
return tex1->Evaluate(si);
else
return tex2->Evaluate(si);
}
private:
std::unique_ptr<TextureMapping3D> mapping;
std::shared_ptr<Texture<T>> tex1, tex2;
};
Noise
通过噪声我们可以为纹理图像添加更多的细节。显然我们不希望噪声是一个完全随机的内容，而是希望能够通过一些方法产生可控的噪声。PBRT定义了噪声函数Noise()来生成这样的可控随机噪声：对于给定的三维空间坐标(x, y, z)，Noise()会生成[-1, 1]区间上的伪随机数。生成的噪声具有有限带宽和最大频率，因此可以控制噪声的频率使得它不会超过采样频率的极限。
Perlin Noise
PRBT中的噪声函数源于Perlin噪声(Perlin noise)。它在整数坐标(x, y, z)上都取0，但每一点都有不同的梯度，然后通过插值就得到了空间中任意点的噪声。
生成Perlin噪声也非常简单。首先要计算(x, y, z)周围8个正数点处的梯度，然后计算梯度与指向采样点的向量的点乘表示各个梯度对采样点的影响，最后对这些点进行三线性插值即可得到噪声。
static constexpr int NoisePermSize = 256;
static int NoisePerm[2 * NoisePermSize] = {
151, 160, 137, 91, 90, 15, 131, 13, 201, 95, 96,
53, 194, 233, 7, 225, 140, 36, 103, 30, 69, 142,
// Remainder of the noise permutation table...
};
inline Float Grad(int x, int y, int z, Float dx, Float dy, Float dz) {
int h = NoisePerm[NoisePerm[NoisePerm[x] + y] + z];
h &= 15;
Float u = h < 8 || h == 12 || h == 13 ? dx : dy;
Float v = h < 4 || h == 12 || h == 13 ? dy : dz;
return ((h & 1) ? -u : u) + ((h & 2) ? -v : v);
}
inline Float NoiseWeight(Float t) {
Float t3 = t * t * t;
Float t4 = t3 * t;
return 6 * t4 * t - 15 * t4 + 10 * t3;
}
Float Noise(Float x, Float y, Float z) {
// Compute noise cell coordinates and offsets
int ix = std::floor(x), iy = std::floor(y), iz = std::floor(z);
Float dx = x - ix, dy = y - iy, dz = z - iz;
// Compute gradient weights
ix &= NoisePermSize - 1;
iy &= NoisePermSize - 1;
iz &= NoisePermSize - 1;
Float w000 = Grad(ix, iy, iz, dx, dy, dz);
Float w100 = Grad(ix+1, iy, iz, dx-1, dy, dz);
Float w010 = Grad(ix, iy+1, iz, dx, dy-1, dz);
Float w110 = Grad(ix+1, iy+1, iz, dx-1, dy-1, dz);
Float w001 = Grad(ix, iy, iz+1, dx, dy, dz-1);
Float w101 = Grad(ix+1, iy, iz+1, dx-1, dy, dz-1);
Float w011 = Grad(ix, iy+1, iz+1, dx, dy-1, dz-1);
Float w111 = Grad(ix+1, iy+1, iz+1, dx-1, dy-1, dz-1);
// Compute trilinear interpolation of weights
Float wx = NoiseWeight(dx), wy = NoiseWeight(dy), wz = NoiseWeight(dz);
Float x00 = Lerp(wx, w000, w100);
Float x10 = Lerp(wx, w010, w110);
Float x01 = Lerp(wx, w001, w101);
Float x11 = Lerp(wx, w011, w111);
Float y0 = Lerp(wy, x00, x10);
Float y1 = Lerp(wy, x01, x11);
return Lerp(wz, y0, y1);
}
Float Noise(const Point3f &p) { return Noise(p.x, p.y, p.z); }
Random Polka Dots
噪声函数的一个经典应用是生成波尔卡圆点(Polka dots)。我们可以把纹理平面划分成网格，然后假设每个网格中有50%的概率存在一个圆点，最后把纹理赋予个模型就产生了波尔卡圆点。
PBRT中使用DotsTexture类表示波尔卡圆点。
template <typename T>
class DotsTexture : public Texture<T> {
public:
DotsTexture(std::unique_ptr<TextureMapping2D> mapping,
const std::shared_ptr<Texture<T>> &outsideDot,
const std::shared_ptr<Texture<T>> &insideDot)
: mapping(std::move(mapping)), outsideDot(outsideDot),
insideDot(insideDot) { }
T Evaluate(const SurfaceInteraction &si) const {
// Compute cell indices for dots
Vector2f dstdx, dstdy;
Point2f st = mapping->Map(si, &dstdx, &dstdy);
int sCell = std::floor(st[0] + .5f), tCell = std::floor(st[1] + .5f);
// Return insideDot result if point is inside dot
if (Noise(sCell + .5f, tCell + .5f) > 0) {
Float radius = .35f;
Float maxShift = 0.5f - radius;
Float sCenter = sCell + maxShift * Noise(sCell + 1.5f, tCell + 2.8f);
Float tCenter = tCell + maxShift * Noise(sCell + 4.5f, tCell + 9.8f);
Vector2f dst = st - Point2f(sCenter, tCenter);
if (dst.LengthSquared() < radius*radius)
return insideDot->Evaluate(si);
}
return outsideDot->Evaluate(si);
}
private:
std::unique_ptr<TextureMapping2D> mapping;
std::shared_ptr<Texture<T>> outsideDot, insideDot;
};
Noise Idioms and Spectral Synthesis
在很多时候我们会通过缩放和加权的方式对已有的噪声进行组合并生成新的程序噪声，这样的过程称为谱生成(spectral synthesis)。生成噪声的谱等于原来噪声谱的加权和：
\[f_s (x) = \sum_i w_i f_i (s_i x)\]Fractional Brownian Motion
通常情况下权重$w_i$和缩放系数$s_i$构成等比数列$w_i = w_{i-1}/2$、$s_i = 2 s_{i-1}$。这样高频成分的权重会不断缩减，对生成噪声的影响也逐渐减少。当我们把这样的递推方式作用在Perlin噪声上就得到了分型布朗运动(fractional Brownian motion, fBm)。
inline Float SmoothStep(Float min, Float max, Float value) {
Float v = Clamp((value - min) / (max - min), 0, 1);
return v * v * (-2 * v + 3);
}
Float FBm(const Point3f &p, const Vector3f &dpdx, const Vector3f &dpdy,
Float omega, int maxOctaves) {
// Compute number of octaves for antialiased FBm
Float len2 = std::max(dpdx.LengthSquared(), dpdy.LengthSquared());
Float n = Clamp(-1 - .5f * Log2(len2), 0, maxOctaves);
int nInt = std::floor(n);
// Compute sum of octaves of noise for FBm
Float sum = 0, lambda = 1, o = 1;
for (int i = 0; i < nInt; ++i) {
sum += o * Noise(lambda * p);
lambda *= 1.99f;
o *= omega;
}
Float nPartial = n - nInt;
sum += o * SmoothStep(.3f, .7f, nPartial) * Noise(lambda * p);
return sum;
}
Turbulence
和fBm类似的另一种噪声是Turbulence噪声，它的组合公式为：
\[f_s (x) = \sum_i w_i \vert f_i (s_i x) \vert\]
Float Turbulence(const Point3f &p, const Vector3f &dpdx,
const Vector3f &dpdy, Float omega, int maxOctaves) {
// Compute number of octaves for antialiased FBm
Float len2 = std::max(dpdx.LengthSquared(), dpdy.LengthSquared());
Float n = Clamp(-1 - .5f * Log2(len2), 0, maxOctaves);
int nInt = std::floor(n);
// Compute sum of octaves of noise for turbulence
Float sum = 0, lambda = 1, o = 1;
for (int i = 0; i < nInt; ++i) {
sum += o * std::abs(Noise(lambda * p));
lambda *= 1.99f;
o *= omega;
}
// Account for contributions of clamped octaves in turbulence
Float nPartial = n - nInt;
sum += o * Lerp(SmoothStep(.3f, .7f, nPartial),
0.2, std::abs(Noise(lambda * p)));
for (int i = nInt; i < maxOctaves; ++i) {
sum += o * 0.2f;
o *= omega;
}
return sum;
}
Bumpy and Wrinkled Textures
利用FBm()和Turbulence()函数可以生成凹凸贴图来模拟物体表面起伏和褶皱的效果。
template <typename T>
class FBmTexture : public Texture<T> {
public:
FBmTexture(std::unique_ptr<TextureMapping3D> mapping, int octaves,
Float omega)
: mapping(std::move(mapping)), omega(omega), octaves(octaves) { }
T Evaluate(const SurfaceInteraction &si) const {
Vector3f dpdx, dpdy;
Point3f P = mapping->Map(si, &dpdx, &dpdy);
return FBm(P, dpdx, dpdy, omega, octaves);
}
private:
std::unique_ptr<TextureMapping3D> mapping;
const Float omega;
const int octaves;
};
template <typename T>
class WrinkledTexture : public Texture<T> {
public:
WrinkledTexture(std::unique_ptr<TextureMapping3D> mapping, int octaves,
Float omega)
: mapping(std::move(mapping)), omega(omega), octaves(octaves) {}
T Evaluate(const SurfaceInteraction &si) const {
Vector3f dpdx, dpdy;
Point3f p = mapping->Map(si, &dpdx, &dpdy);
return Turbulence(p, dpdx, dpdy, omega, octaves);
}
private:
std::unique_ptr<TextureMapping3D> mapping;
Float omega;
int octaves;
};
Windy Waves
FBm()函数还可以用来模拟水面涟漪的效果。
template <typename T>
class WindyTexture : public Texture<T> {
public:
WindyTexture(std::unique_ptr<TextureMapping3D> mapping)
: mapping(std::move(mapping)) { }
T Evaluate(const SurfaceInteraction &si) const {
Vector3f dpdx, dpdy;
Point3f P = mapping->Map(si, &dpdx, &dpdy);
Float windStrength = FBm(.1f * P, .1f * dpdx, .1f * dpdy, .5, 3);
Float waveHeight = FBm(P, dpdx, dpdy, .5, 6);
return std::abs(windStrength) * waveHeight;
}
private:
std::unique_ptr<TextureMapping3D> mapping;
};
Marble
噪声函数的一个经典应用是用来模拟大理石纹理。
class MarbleTexture : public Texture<Spectrum> {
public:
MarbleTexture(std::unique_ptr<TextureMapping3D> mapping, int octaves,
Float omega, Float scale, Float variation)
: mapping(std::move(mapping)), octaves(octaves), omega(omega),
scale(scale), variation(variation) { }
Spectrum Evaluate(const SurfaceInteraction &si) const {
Vector3f dpdx, dpdy;
Point3f p = mapping->Map(si, &dpdx, &dpdy);
p *= scale;
Float marble = p.y + variation *
FBm(p, scale * dpdx, scale * dpdy, omega, octaves);
Float t = .5f + .5f * std::sin(marble);
// Evaluate marble spline at t
static Float c[][3] = { { .58f, .58f, .6f }, { .58f, .58f, .6f }, { .58f, .58f, .6f },
{ .5f, .5f, .5f }, { .6f, .59f, .58f }, { .58f, .58f, .6f },
{ .58f, .58f, .6f }, { .2f, . 2f, .33f }, { .58f, .58f, .6f }, };
#define NC sizeof(c) / sizeof(c[0])
#define NSEG (NC-3)
int first = std::floor(t * NSEG);
t = (t * NSEG - first);
Spectrum c0 = Spectrum::FromRGB(c[first]);
Spectrum c1 = Spectrum::FromRGB(c[first+1]);
Spectrum c2 = Spectrum::FromRGB(c[first+2]);
Spectrum c3 = Spectrum::FromRGB(c[first+3]);
// Bezier spline evaluated with de Castilejau's algorithm
Spectrum s0 = (1.f - t) * c0 + t * c1;
Spectrum s1 = (1.f - t) * c1 + t * c2;
Spectrum s2 = (1.f - t) * c2 + t * c3;
s0 = (1.f - t) * s0 + t * s1;
s1 = (1.f - t) * s1 + t * s2;
// Extra scale of 1.5 to increase variation among colors
return 1.5f * ((1.f - t) * s0 + t * s1);
}
private:
std::unique_ptr<TextureMapping3D> mapping;
const int octaves;
const Float omega, scale, variation;
};
Reference
PREVIOUSMaterials