Update from book source. No functional change.

src/pbrt/util/colorspace.cpp

@@ -20,26 +20,17 @@ extern const SquareMatrix<3> LMSFromXYZ, XYZFromLMS;
 
 // RGBColorSpace Method Definitions
 RGBColorSpace::RGBColorSpace(Point2f r, Point2f g, Point2f b, SpectrumHandle illuminant,
-                             const RGBToSpectrumTable *rgbToSpectrumTable,
-                             Allocator alloc)
-    : r(r),
-      g(g),
-      b(b),
-      illuminant(illuminant, alloc),
-      rgbToSpectrumTable(rgbToSpectrumTable) {
-    // Compute whitepoint primaries
+                             const RGBToSpectrumTable *rgbToSpec, Allocator alloc)
+    : r(r), g(g), b(b), illuminant(illuminant, alloc), rgbToSpectrumTable(rgbToSpec) {
+    // Compute whitepoint primaries and XYZ coordinates
     XYZ W = SpectrumToXYZ(illuminant);
     w = W.xy();
-
-    // Compute XYZ coordinates of primaries
     XYZ R = XYZ::FromxyY(r), G = XYZ::FromxyY(g), B = XYZ::FromxyY(b);
 
-    // Compute _XYZFromRGB_ matrix for color space
+    // Initialize XYZ color space conversion matrices
     SquareMatrix<3> rgb(R.X, G.X, B.X, R.Y, G.Y, B.Y, R.Z, G.Z, B.Z);
     XYZ C = *Inverse(rgb) * W;
     XYZFromRGB = rgb * SquareMatrix<3>::Diag(C[0], C[1], C[2]);
-
-    // Compute _RGBFromXYZ_ matrix for color space
     RGBFromXYZ = *Inverse(XYZFromRGB);
 }
 

src/pbrt/util/transform.cpp

@@ -79,12 +79,12 @@ Transform RotateZ(Float theta) {
 // clang-format on
 
 Transform LookAt(const Point3f &pos, const Point3f &look, const Vector3f &up) {
-    SquareMatrix<4> cameraToWorld;
+    SquareMatrix<4> worldFromCamera;
     // Initialize fourth column of viewing matrix
-    cameraToWorld[0][3] = pos.x;
-    cameraToWorld[1][3] = pos.y;
-    cameraToWorld[2][3] = pos.z;
-    cameraToWorld[3][3] = 1;
+    worldFromCamera[0][3] = pos.x;
+    worldFromCamera[1][3] = pos.y;
+    worldFromCamera[2][3] = pos.z;
+    worldFromCamera[3][3] = 1;
 
     // Initialize first three columns of viewing matrix
     Vector3f dir = Normalize(look - pos);
@@ -95,24 +95,24 @@ Transform LookAt(const Point3f &pos, const Point3f &look, const Vector3f &up) {
                   up.x, up.y, up.z, dir.x, dir.y, dir.z);
     Vector3f right = Normalize(Cross(Normalize(up), dir));
     Vector3f newUp = Cross(dir, right);
-    cameraToWorld[0][0] = right.x;
-    cameraToWorld[1][0] = right.y;
-    cameraToWorld[2][0] = right.z;
-    cameraToWorld[3][0] = 0.;
-    cameraToWorld[0][1] = newUp.x;
-    cameraToWorld[1][1] = newUp.y;
-    cameraToWorld[2][1] = newUp.z;
-    cameraToWorld[3][1] = 0.;
-    cameraToWorld[0][2] = dir.x;
-    cameraToWorld[1][2] = dir.y;
-    cameraToWorld[2][2] = dir.z;
-    cameraToWorld[3][2] = 0.;
-
-    pstd::optional<SquareMatrix<4>> worldToCamera = Inverse(cameraToWorld);
+    worldFromCamera[0][0] = right.x;
+    worldFromCamera[1][0] = right.y;
+    worldFromCamera[2][0] = right.z;
+    worldFromCamera[3][0] = 0.;
+    worldFromCamera[0][1] = newUp.x;
+    worldFromCamera[1][1] = newUp.y;
+    worldFromCamera[2][1] = newUp.z;
+    worldFromCamera[3][1] = 0.;
+    worldFromCamera[0][2] = dir.x;
+    worldFromCamera[1][2] = dir.y;
+    worldFromCamera[2][2] = dir.z;
+    worldFromCamera[3][2] = 0.;
+
+    pstd::optional<SquareMatrix<4>> cameraFromWorld = Inverse(worldFromCamera);
 #ifdef PBRT_DEBUG_BUILD
-    DCHECK(worldToCamera);
+    DCHECK(cameraFromWorld);
 #endif
-    return Transform(*worldToCamera, cameraToWorld);
+    return Transform(*cameraFromWorld, worldFromCamera);
 }
 
 Transform Orthographic(Float zNear, Float zFar) {

src/pbrt/util/vecmath.h

@@ -877,8 +877,7 @@ PBRT_CPU_GPU inline auto Length(const Vector2<T> &v) -> typename TupleLength<T>:
 }
 
 template <typename T>
-PBRT_CPU_GPU inline auto Normalize(const Vector2<T> &v)
-    -> Vector2<typename TupleLength<T>::type> {
+PBRT_CPU_GPU inline auto Normalize(const Vector2<T> &v) {
     return v / Length(v);
 }
 
@@ -918,7 +917,7 @@ PBRT_CPU_GPU inline Vector3<T> Cross(const Normal3<T> &v1, const Vector3<T> &v2)
 
 template <typename T>
 PBRT_CPU_GPU inline T LengthSquared(const Vector3<T> &v) {
-    return v.x * v.x + v.y * v.y + v.z * v.z;
+    return Sqr(v.x) + Sqr(v.y) + Sqr(v.z);
 }
 
 template <typename T>
@@ -928,7 +927,7 @@ PBRT_CPU_GPU inline auto Length(const Vector3<T> &v) -> typename TupleLength<T>:
 }
 
 template <typename T>
-PBRT_CPU_GPU inline Vector3<T> Normalize(const Vector3<T> &v) {
+PBRT_CPU_GPU inline auto Normalize(const Vector3<T> &v) {
     return v / Length(v);
 }
 
@@ -963,8 +962,8 @@ PBRT_CPU_GPU inline Float AngleBetween(const Normal3<T> &a, const Normal3<T> &b)
 }
 
 template <typename T>
-PBRT_CPU_GPU inline Vector3<T> GramSchmidt(const Vector3<T> &a, const Vector3<T> &b) {
-    return a - Dot(a, b) * b;
+PBRT_CPU_GPU inline Vector3<T> GramSchmidt(const Vector3<T> &v, const Vector3<T> &w) {
+    return v - Dot(v, w) * w;
 }
 
 template <typename T>
@@ -1025,8 +1024,7 @@ PBRT_CPU_GPU inline auto Length(const Normal3<T> &n) -> typename TupleLength<T>:
 }
 
 template <typename T>
-PBRT_CPU_GPU inline auto Normalize(const Normal3<T> &n)
-    -> Normal3<typename TupleLength<T>::type> {
+PBRT_CPU_GPU inline auto Normalize(const Normal3<T> &n) {
     return n / Length(n);
 }
 
@@ -1509,7 +1507,7 @@ PBRT_CPU_GPU inline auto DistanceSquared(const Point3<T> &p, const Bounds3<U> &b
     TDist dx = std::max<TDist>({0, b.pMin.x - p.x, p.x - b.pMax.x});
     TDist dy = std::max<TDist>({0, b.pMin.y - p.y, p.y - b.pMax.y});
     TDist dz = std::max<TDist>({0, b.pMin.z - p.z, p.z - b.pMax.z});
-    return dx * dx + dy * dy + dz * dz;
+    return Sqr(dx) + Sqr(dy) + Sqr(dz);
 }
 
 template <typename T, typename U>
@@ -1620,21 +1618,59 @@ PBRT_CPU_GPU inline Bounds2<T> Union(const Bounds2<T> &b, const Point2<T> &p) {
 }
 
 // Spherical Geometry Inline Functions
-PBRT_CPU_GPU
-inline Vector3f SphericalDirection(Float sinTheta, Float cosTheta, Float phi) {
+PBRT_CPU_GPU inline Float SphericalTriangleArea(Vector3f a, Vector3f b, Vector3f c) {
+    // Compute normalized cross products of all direction pairs
+    Vector3f n_ab = Cross(a, b), n_bc = Cross(b, c), n_ca = Cross(c, a);
+    if (LengthSquared(n_ab) == 0 || LengthSquared(n_bc) == 0 || LengthSquared(n_ca) == 0)
+        return {};
+    n_ab = Normalize(n_ab);
+    n_bc = Normalize(n_bc);
+    n_ca = Normalize(n_ca);
+
+    // Compute angles $\alpha$, $\beta$, and $\gamma$ at spherical triangle vertices
+    Float alpha = AngleBetween(n_ab, -n_ca);
+    Float beta = AngleBetween(n_bc, -n_ab);
+    Float gamma = AngleBetween(n_ca, -n_bc);
+
+    return std::abs(alpha + beta + gamma - Pi);
+}
+
+PBRT_CPU_GPU inline Float SphericalQuadArea(Vector3f a, Vector3f b, Vector3f c,
+                                            Vector3f d);
+
+PBRT_CPU_GPU inline Float SphericalQuadArea(Vector3f a, Vector3f b, Vector3f c,
+                                            Vector3f d) {
+    Vector3f axb = Cross(a, b), bxc = Cross(b, c);
+    Vector3f cxd = Cross(c, d), dxa = Cross(d, a);
+    if (LengthSquared(axb) == 0 || LengthSquared(bxc) == 0 || LengthSquared(cxd) == 0 ||
+        LengthSquared(dxa) == 0)
+        return 0;
+    axb = Normalize(axb);
+    bxc = Normalize(bxc);
+    cxd = Normalize(cxd);
+    dxa = Normalize(dxa);
+
+    Float alpha = AngleBetween(dxa, -axb);
+    Float beta = AngleBetween(axb, -bxc);
+    Float gamma = AngleBetween(bxc, -cxd);
+    Float delta = AngleBetween(cxd, -dxa);
+
+    return std::abs(alpha + beta + gamma + delta - 2 * Pi);
+}
+
+PBRT_CPU_GPU inline Vector3f SphericalDirection(Float sinTheta, Float cosTheta,
+                                                Float phi) {
     DCHECK(sinTheta >= -1.0001 && sinTheta <= 1.0001);
     DCHECK(cosTheta >= -1.0001 && cosTheta <= 1.0001);
     return Vector3f(Clamp(sinTheta, -1, 1) * std::cos(phi),
                     Clamp(sinTheta, -1, 1) * std::sin(phi), Clamp(cosTheta, -1, 1));
 }
 
-PBRT_CPU_GPU
-inline Float SphericalTheta(const Vector3f &v) {
+PBRT_CPU_GPU inline Float SphericalTheta(const Vector3f &v) {
     return SafeACos(v.z);
 }
 
-PBRT_CPU_GPU
-inline Float SphericalPhi(const Vector3f &v) {
+PBRT_CPU_GPU inline Float SphericalPhi(const Vector3f &v) {
     Float p = std::atan2(v.y, v.x);
     return (p < 0) ? (p + 2 * Pi) : p;
 }
@@ -1643,7 +1679,7 @@ PBRT_CPU_GPU inline Float CosTheta(const Vector3f &w) {
     return w.z;
 }
 PBRT_CPU_GPU inline Float Cos2Theta(const Vector3f &w) {
-    return w.z * w.z;
+    return Sqr(w.z);
 }
 PBRT_CPU_GPU inline Float AbsCosTheta(const Vector3f &w) {
     return std::abs(w.z);
@@ -1663,13 +1699,11 @@ PBRT_CPU_GPU inline Float Tan2Theta(const Vector3f &w) {
     return Sin2Theta(w) / Cos2Theta(w);
 }
 
-PBRT_CPU_GPU
-inline Float CosPhi(const Vector3f &w) {
+PBRT_CPU_GPU inline Float CosPhi(const Vector3f &w) {
     Float sinTheta = SinTheta(w);
     return (sinTheta == 0) ? 1 : Clamp(w.x / sinTheta, -1, 1);
 }
-PBRT_CPU_GPU
-inline Float SinPhi(const Vector3f &w) {
+PBRT_CPU_GPU inline Float SinPhi(const Vector3f &w) {
     Float sinTheta = SinTheta(w);
     return (sinTheta == 0) ? 0 : Clamp(w.y / sinTheta, -1, 1);
 }
@@ -1682,48 +1716,6 @@ inline Float CosDPhi(const Vector3f &wa, const Vector3f &wb) {
     return Clamp((wa.x * wb.x + wa.y * wb.y) / std::sqrt(waxy * wbxy), -1, 1);
 }
 
-PBRT_CPU_GPU inline Float SphericalTriangleArea(const Vector3f &a, const Vector3f &b,
-                                                const Vector3f &c) {
-    // Compute normalized cross products of all direction pairs
-    Vector3f n_ab = Cross(a, b), n_bc = Cross(b, c), n_ca = Cross(c, a);
-    if (LengthSquared(n_ab) == 0 || LengthSquared(n_bc) == 0 || LengthSquared(n_ca) == 0)
-        return {};
-    n_ab = Normalize(n_ab);
-    n_bc = Normalize(n_bc);
-    n_ca = Normalize(n_ca);
-
-    // Compute angles $\alpha$, $\beta$, and $\gamma$ at spherical triangle vertices
-    Float alpha = AngleBetween(n_ab, -n_ca);
-    Float beta = AngleBetween(n_bc, -n_ab);
-    Float gamma = AngleBetween(n_ca, -n_bc);
-
-    return std::abs(alpha + beta + gamma - Pi);
-}
-
-PBRT_CPU_GPU inline Float SphericalQuadArea(const Vector3f &a, const Vector3f &b,
-                                            const Vector3f &c, const Vector3f &d);
-
-PBRT_CPU_GPU
-inline Float SphericalQuadArea(const Vector3f &a, const Vector3f &b, const Vector3f &c,
-                               const Vector3f &d) {
-    Vector3f axb = Cross(a, b), bxc = Cross(b, c);
-    Vector3f cxd = Cross(c, d), dxa = Cross(d, a);
-    if (LengthSquared(axb) == 0 || LengthSquared(bxc) == 0 || LengthSquared(cxd) == 0 ||
-        LengthSquared(dxa) == 0)
-        return 0;
-    axb = Normalize(axb);
-    bxc = Normalize(bxc);
-    cxd = Normalize(cxd);
-    dxa = Normalize(dxa);
-
-    Float alpha = AngleBetween(dxa, -axb);
-    Float beta = AngleBetween(axb, -bxc);
-    Float gamma = AngleBetween(bxc, -cxd);
-    Float delta = AngleBetween(cxd, -dxa);
-
-    return std::abs(alpha + beta + gamma + delta - 2 * Pi);
-}
-
 PBRT_CPU_GPU
 inline bool SameHemisphere(const Vector3f &w, const Vector3f &wp) {
     return w.z * wp.z > 0;
@@ -1740,14 +1732,15 @@ class OctahedralVector {
     // OctahedralVector Public Methods
     OctahedralVector() = default;
     PBRT_CPU_GPU
-    OctahedralVector(const Vector3f &v) {
-        Float invL1Norm = 1 / (std::abs(v.x) + std::abs(v.y) + std::abs(v.z));
-        if (v.z < 0.0f) {
-            x = Encode((1 - std::abs(v.y * invL1Norm)) * SignNotZero(v.x));
-            y = Encode((1 - std::abs(v.x * invL1Norm)) * SignNotZero(v.y));
+    OctahedralVector(Vector3f v) {
+        v /= std::abs(v.x) + std::abs(v.y) + std::abs(v.z);
+        if (v.z >= 0) {
+            x = Encode(v.x);
+            y = Encode(v.y);
         } else {
-            x = Encode(v.x * invL1Norm);
-            y = Encode(v.y * invL1Norm);
+            // Encode octahedral vector with $z < 0$
+            x = Encode((1 - std::abs(v.y)) * SignNotZero(v.x));
+            y = Encode((1 - std::abs(v.x)) * SignNotZero(v.y));
         }
     }
 
@@ -1757,11 +1750,13 @@ class OctahedralVector {
         v.x = -1 + 2 * (x / 65535.f);
         v.y = -1 + 2 * (y / 65535.f);
         v.z = 1 - (std::abs(v.x) + std::abs(v.y));
+        // Reparameterize directions in the $z<0$ portion of the octahedron
         if (v.z < 0) {
             Float xo = v.x;
             v.x = (1 - std::abs(v.y)) * SignNotZero(xo);
             v.y = (1 - std::abs(xo)) * SignNotZero(v.y);
         }
+
         return Normalize(v);
     }
 
@@ -1771,14 +1766,14 @@ class OctahedralVector {
 
   private:
     // OctahedralVector Private Methods
-    PBRT_CPU_GPU
-    static Float SignNotZero(Float v) { return (v < 0) ? -1 : 1; }
-
     PBRT_CPU_GPU
     static uint16_t Encode(Float f) {
         return std::round(Clamp((f + 1) / 2, 0, 1) * 65535.f);
     }
 
+    PBRT_CPU_GPU
+    static Float SignNotZero(Float v) { return (v < 0) ? -1 : 1; }
+
     // OctahedralVector Private Members
     uint16_t x, y;
 };