A.5 Mathematical Routines

This section describes a number of useful mathematical functions and classes that support basic operations in pbrt, such as solving small linear systems, manipulating matrices, and linear interpolation.

The Lerp() function linearly interpolates between the two provided values.

<<Global Inline Functions>>+=  
inline Float Lerp(Float t, Float v1, Float v2) { return (1 - t) * v1 + t * v2; }

A.5.1 Solving Quadratic Equations

The Quadratic() function finds solutions of the quadratic equation a t squared plus b t plus c equals 0 ; the Boolean return value indicates whether solutions were found.

<<Global Inline Functions>>+= 
inline bool Quadratic(Float a, Float b, Float c, Float *t0, Float *t1) { <<Find quadratic discriminant>> 
double discrim = (double)b * (double)b - 4 * (double)a * (double)c; if (discrim < 0) return false; double rootDiscrim = std::sqrt(discrim);
<<Compute quadratic t values>> 
double q; if (b < 0) q = -.5 * (b - rootDiscrim); else q = -.5 * (b + rootDiscrim); *t0 = q / a; *t1 = c / q; if (*t0 > *t1) std::swap(*t0, *t1); return true;
}

The implementation always uses double-precision floating-point values regardless of the type of Float in order to return a result with minimal floating-point error. If the discriminant ( b squared minus 4 a c ) is negative, then there are no real roots and the function returns false.

<<Find quadratic discriminant>>= 
double discrim = (double)b * (double)b - 4 * (double)a * (double)c; if (discrim < 0) return false; double rootDiscrim = std::sqrt(discrim);

The usual version of the quadratic equation can give poor numerical precision when b almost-equals plus-or-minus StartRoot b squared minus 4 a c EndRoot due to cancellation error. It can be rewritten algebraically to a more stable form:

StartLayout 1st Row 1st Column t 0 2nd Column equals StartFraction q Over a EndFraction 2nd Row 1st Column t 1 2nd Column equals StartFraction c Over q EndFraction comma EndLayout

where

q equals StartLayout Enlarged left-brace 1st Row 1st Column minus .5 left-parenthesis b minus StartRoot b squared minus 4 a c EndRoot right-parenthesis 2nd Column b less-than 0 comma 2nd Row 1st Column minus .5 left-parenthesis b plus StartRoot b squared minus 4 a c EndRoot right-parenthesis 2nd Column otherwise period EndLayout

<<Compute quadratic t values>>= 
double q; if (b < 0) q = -.5 * (b - rootDiscrim); else q = -.5 * (b + rootDiscrim); *t0 = q / a; *t1 = c / q; if (*t0 > *t1) std::swap(*t0, *t1); return true;

A.5.2 2 times 2 Linear Systems

There are a number of places throughout pbrt where we need to solve a 2 times 2 linear system upper A x equals upper B of the form

Start 2 By 2 Matrix 1st Row 1st Column a 00 2nd Column a 01 2nd Row 1st Column a 10 2nd Column a 11 EndMatrix StartBinomialOrMatrix x 0 Choose x 1 EndBinomialOrMatrix equals StartBinomialOrMatrix b 0 Choose b 1 EndBinomialOrMatrix

for values x 0 and x 1 . The SolveLinearSystem2x2() routine finds the closed-form solution to such a system. It returns true if it was successful and false if the determinant of upper A is very small, indicating that the system is numerically ill-conditioned and either not solvable or likely to have unacceptable floating-point errors. In this case, no solution is returned.

<<Matrix4x4 Method Definitions>>= 
bool SolveLinearSystem2x2(const Float A[2][2], const Float B[2], Float *x0, Float *x1) { Float det = A[0][0] * A[1][1] - A[0][1] * A[1][0]; if (std::abs(det) < 1e-10f) return false; *x0 = (A[1][1] * B[0] - A[0][1] * B[1]) / det; *x1 = (A[0][0] * B[1] - A[1][0] * B[0]) / det; if (std::isnan(*x0) || std::isnan(*x1)) return false; return true; }

A.5.3 4 times 4 Matrices

The Matrix4x4 structure provides a low-level representation of 4 times 4 matrices. It is an integral part of the Transform class.

<<Matrix4x4 Declarations>>= 
struct Matrix4x4 { <<Matrix4x4 Public Methods>> 
Matrix4x4() { m[0][0] = m[1][1] = m[2][2] = m[3][3] = 1.f; m[0][1] = m[0][2] = m[0][3] = m[1][0] = m[1][2] = m[1][3] = m[2][0] = m[2][1] = m[2][3] = m[3][0] = m[3][1] = m[3][2] = 0.f; } Matrix4x4(Float mat[4][4]); Matrix4x4(Float t00, Float t01, Float t02, Float t03, Float t10, Float t11, Float t12, Float t13, Float t20, Float t21, Float t22, Float t23, Float t30, Float t31, Float t32, Float t33); bool operator==(const Matrix4x4 &m2) const { for (int i = 0; i < 4; ++i) for (int j = 0; j < 4; ++j) if (m[i][j] != m2.m[i][j]) return false; return true; } bool operator!=(const Matrix4x4 &m2) const { for (int i = 0; i < 4; ++i) for (int j = 0; j < 4; ++j) if (m[i][j] != m2.m[i][j]) return true; return false; } friend Matrix4x4 Transpose(const Matrix4x4 &); void Print(FILE *f) const { fprintf(f, "[ "); for (int i = 0; i < 4; ++i) { fprintf(f, " [ "); for (int j = 0; j < 4; ++j) { fprintf(f, "%f", m[i][j]); if (j != 3) fprintf(f, ", "); } fprintf(f, " ]\n"); } fprintf(f, " ] "); } static Matrix4x4 Mul(const Matrix4x4 &m1, const Matrix4x4 &m2) { Matrix4x4 r; for (int i = 0; i < 4; ++i) for (int j = 0; j < 4; ++j) r.m[i][j] = m1.m[i][0] * m2.m[0][j] + m1.m[i][1] * m2.m[1][j] + m1.m[i][2] * m2.m[2][j] + m1.m[i][3] * m2.m[3][j]; return r; } friend Matrix4x4 Inverse(const Matrix4x4 &);
Float m[4][4]; };

The default constructor, not shown here, sets the matrix to the identity matrix. The Matrix4x4 implementation also provides constructors that allow the user to pass an array of floats or 16 individual floats to initialize a Matrix4x4:

<<Matrix4x4 Public Methods>>= 
Matrix4x4(Float mat[4][4]); Matrix4x4(Float t00, Float t01, Float t02, Float t03, Float t10, Float t11, Float t12, Float t13, Float t20, Float t21, Float t22, Float t23, Float t30, Float t31, Float t32, Float t33);

The implementations of operators that test for equality and inequality are straightforward and not included in the text here.

The Matrix4x4 class supports a few low-level matrix operations. For example, Transpose() returns a new matrix that is the transpose of the original matrix.

<<Matrix4x4 Method Definitions>>+= 
Matrix4x4 Transpose(const Matrix4x4 &m) { return Matrix4x4(m.m[0][0], m.m[1][0], m.m[2][0], m.m[3][0], m.m[0][1], m.m[1][1], m.m[2][1], m.m[3][1], m.m[0][2], m.m[1][2], m.m[2][2], m.m[3][2], m.m[0][3], m.m[1][3], m.m[2][3], m.m[3][3]); }

The product of two matrices bold upper M 1 and bold upper M 2 is computed by setting the left-parenthesis i comma j right-parenthesis th element of the result to the inner product of the i th row of bold upper M 1 with the j th column of bold upper M 2 .

<<Matrix4x4 Public Methods>>+=  
static Matrix4x4 Mul(const Matrix4x4 &m1, const Matrix4x4 &m2) { Matrix4x4 r; for (int i = 0; i < 4; ++i) for (int j = 0; j < 4; ++j) r.m[i][j] = m1.m[i][0] * m2.m[0][j] + m1.m[i][1] * m2.m[1][j] + m1.m[i][2] * m2.m[2][j] + m1.m[i][3] * m2.m[3][j]; return r; }

Finally, Inverse() returns the inverse of the matrix. The implementation (not shown here) uses a numerically stable Gauss–Jordan elimination routine to compute the inverse.

<<Matrix4x4 Public Methods>>+= 
friend Matrix4x4 Inverse(const Matrix4x4 &);