Vectors

2.2 Vectors

pbrt provides both 2D and 3D vector classes. Both are parameterized by the type of the underlying vector element, thus making it easy to instantiate vectors of both integer and floating-point types.

<<Vector Declarations>>=

template <typename T> class Vector2 { public: <<Vector2 Public Methods>>

Vector2() { x = y = 0; } Vector2(T xx, T yy) : x(xx), y(yy) { Assert(!HasNaNs()); } bool HasNaNs() const { return std::isnan(x) || std::isnan(y); } explicit Vector2(const Point2<T> &p); explicit Vector2(const Point3<T> &p); #ifndef NDEBUG // The default versions of these are fine for release builds; for debug // we define them so that we can add the Assert checks. Vector2(const Vector2<T> &v) { Assert(!v.HasNaNs()); x = v.x; y = v.y; } Vector2<T> &operator=(const Vector2<T> &v) { Assert(!v.HasNaNs()); x = v.x; y = v.y; return *this; } #endif // !NDEBUG friend std::ostream& operator<<(std::ostream& os, const Vector2<T> &v) { os << "[" << v.x << ", " << v.y << "]"; return os; } Vector2<T> operator+(const Vector2<T> &v) const { Assert(!v.HasNaNs()); return Vector2(x + v.x, y + v.y); } Vector2<T>& operator+=(const Vector2<T> &v) { Assert(!v.HasNaNs()); x += v.x; y += v.y; return *this; } Vector2<T> operator-(const Vector2<T> &v) const { Assert(!v.HasNaNs()); return Vector2(x - v.x, y - v.y); } Vector2<T>& operator-=(const Vector2<T> &v) { Assert(!v.HasNaNs()); x -= v.x; y -= v.y; return *this; } bool operator==(const Vector2<T> &v) const { return x == v.x && y == v.y; } bool operator!=(const Vector2<T> &v) const { return x != v.x || y != v.y; } Vector2<T> operator*(T f) const { return Vector2<T>(f*x, f*y); } Vector2<T> &operator*=(T f) { Assert(!std::isnan(f)); x *= f; y *= f; return *this; } Vector2<T> operator/(T f) const { Assert(f != 0); Float inv = (Float)1 / f; return Vector2<T>(x * inv, y * inv); } Vector2<T> &operator/=(T f) { Assert(f != 0); Float inv = (Float)1 / f; x *= inv; y *= inv; return *this; } Vector2<T> operator-() const { return Vector2<T>(-x, -y); } T operator[](int i) const { Assert(i >= 0 && i <= 1); if (i == 0) return x; return y; } T &operator[](int i) { Assert(i >= 0 && i <= 1); if (i == 0) return x; return y; } Float LengthSquared() const { return x*x + y*y; } Float Length() const { return std::sqrt(LengthSquared()); }

<<Vector2 Public Data>>

T x, y;

};

<<Vector Declarations>>+=

template <typename T> class Vector3 { public: <<Vector3 Public Methods>>

T operator[](int i) const { Assert(i >= 0 && i <= 2); if (i == 0) return x; if (i == 1) return y; return z; } T &operator[](int i) { Assert(i >= 0 && i <= 2); if (i == 0) return x; if (i == 1) return y; return z; } Vector3() { x = y = z = 0; } Vector3(T x, T y, T z) : x(x), y(y), z(z) { Assert(!HasNaNs()); } bool HasNaNs() const { return std::isnan(x) || std::isnan(y) || std::isnan(z); } explicit Vector3(const Point3<T> &p); #ifndef NDEBUG // The default versions of these are fine for release builds; for debug // we define them so that we can add the Assert checks. Vector3(const Vector3<T> &v) { Assert(!v.HasNaNs()); x = v.x; y = v.y; z = v.z; } Vector3<T> &operator=(const Vector3<T> &v) { Assert(!v.HasNaNs()); x = v.x; y = v.y; z = v.z; return *this; } #endif // !NDEBUG friend std::ostream& operator<<(std::ostream& os, const Vector3<T> &v) { os << "[" << v.x << ", " << v.y << ", " << v.z << "]"; return os; } Vector3<T> operator+(const Vector3<T> &v) const { return Vector3(x + v.x, y + v.y, z + v.z); } Vector3<T>& operator+=(const Vector3<T> &v) { x += v.x; y += v.y; z += v.z; return *this; } Vector3<T> operator-(const Vector3<T> &v) const { return Vector3(x - v.x, y - v.y, z - v.z); } Vector3<T>& operator-=(const Vector3<T> &v) { x -= v.x; y -= v.y; z -= v.z; return *this; } bool operator==(const Vector3<T> &v) const { return x == v.x && y == v.y && z == v.z; } bool operator!=(const Vector3<T> &v) const { return x != v.x || y != v.y || z != v.z; } Vector3<T> operator*(T s) const { return Vector3<T>(s*x, s*y, s*z); } Vector3<T> &operator*=(T s) { x *= s; y *= s; z *= s; return *this; } Vector3<T> operator/(T f) const { Assert(f != 0); Float inv = (Float)1 / f; return Vector3<T>(x * inv, y * inv, z * inv); } Vector3<T> &operator/=(T f) { Assert(f != 0); Float inv = (Float)1 / f; x *= inv; y *= inv; z *= inv; return *this; } Vector3<T> operator-() const { return Vector3<T>(-x, -y, -z); } Float LengthSquared() const { return x * x + y * y + z * z; } Float Length() const { return std::sqrt(LengthSquared()); } explicit Vector3(const Normal3<T> &n);

<<Vector3 Public Data>>

T x, y, z;

};

In the following, we will generally only include implementations of Vector3 methods; all have Vector2 parallels that have straightforward implementation differences.

A vector is represented with a tuple of components that gives its representation in terms of the , , (in 3D) axes of the space it is defined in. The individual components of a 3D vector will be written , , and .

<<Vector2 Public Data>>=

T x, y;

<<Vector3 Public Data>>=

T x, y, z;

An alternate implementation would be to have a single template class that is also parameterized with an integer number of dimensions and to represent the coordinates with an array of that many T values. While this approach would reduce the total amount of code, individual components of the vector couldn’t be accessed as v.x and so forth. We believe that in this case, a bit more code in the vector implementations is worthwhile in return for more transparent access to elements.

However, some routines do find it useful to be able to easily loop over the components of vectors; the vector classes also provide a C++ operator to index into the components so that, given a vector v, v[0] == v.x and so forth.

<<Vector3 Public Methods>>=

For convenience, a number of widely used types of vectors are given a typedef, so that they have more concise names in code elsewhere.

<<Vector Declarations>>+=

typedef Vector2<Float> Vector2f; typedef Vector2<int> Vector2i; typedef Vector3<Float> Vector3f; typedef Vector3<int> Vector3i;

Readers who have been exposed to object-oriented design may question our decision to make the vector element data publicly accessible. Typically, data members are only accessible inside their class, and external code that wishes to access or modify the contents of a class must do so through a well-defined API of selector and mutator functions. Although we generally agree with this design principle (though see the discussion of data-oriented design in the “Further Reading” section of Chapter 1), it is not appropriate here. The purpose of selector and mutator functions is to hide the class’s internal implementation details. In the case of vectors, hiding this basic part of their design gains nothing and adds bulk to code that uses them.

By default, the values are set to zero, although the user of the class can optionally supply values for each of the components. If the user does supply values, we check that none of them has the floating-point “not a number” (NaN) value using the Assert() macro. When compiled in optimized mode, this macro disappears from the compiled code, saving the expense of verifying this case. NaNs almost certainly indicate a bug in the system; if a NaN is generated by some computation, we’d like to catch it as soon as possible in order to make isolating its source easier. (See Section 3.9.1 for more discussion of NaN values.)

<<Vector3 Public Methods>>+=

Vector3() { x = y = z = 0; } Vector3(T x, T y, T z) : x(x), y(y), z(z) { Assert(!HasNaNs()); }

The code to check for NaNs just calls the std::isnan() function on each of the , , and components.

<<Vector3 Public Methods>>+=

bool HasNaNs() const { return std::isnan(x) || std::isnan(y) || std::isnan(z); }

Addition and subtraction of vectors are done component-wise. The usual geometric interpretation of vector addition and subtraction is shown in Figures 2.3 and 2.4.

Figure 2.3: (left) Vector addition: . (right) Notice that the sum forms the diagonal of the parallelogram formed by and , which shows the commutativity of vector addition: .

Figure 2.4: (left) Vector subtraction. (right) If we consider the parallelogram formed by two vectors, the diagonals are given by (dashed line) and (not shown).

<<Vector3 Public Methods>>+=

Vector3<T> operator+(const Vector3<T> &v) const { return Vector3(x + v.x, y + v.y, z + v.z); } Vector3<T>& operator+=(const Vector3<T> &v) { x += v.x; y += v.y; z += v.z; return *this; }

The code for subtracting two vectors is similar and therefore not shown here.

A vector can be multiplied component-wise by a scalar, thereby changing its length. Three functions are needed in order to cover all of the different ways that this operation may be written in source code (i.e., v*s, s*v, and v *= s):

<<Vector3 Public Methods>>+=

Vector3<T> operator*(T s) const { return Vector3<T>(s*x, s*y, s*z); } Vector3<T> &operator*=(T s) { x *= s; y *= s; z *= s; return *this; }

<<Geometry Inline Functions>>=

template <typename T> inline Vector3<T> operator*(T s, const Vector3<T> &v) { return v * s; }

Similarly, a vector can be divided component-wise by a scalar. The code for scalar division is similar to scalar multiplication, although division of a scalar by a vector is not well defined and so is not permitted.

In the implementation of these methods, we use a single division to compute the scalar’s reciprocal and then perform three component-wise multiplications. This is a useful trick for avoiding division operations, which are generally much slower than multiplies on modern CPUs.

We use the Assert() macro to make sure that the provided divisor is not zero; this should never happen and would indicate a bug elsewhere in the system.

<<Vector3 Public Methods>>+=

Vector3<T> operator/(T f) const { Assert(f != 0); Float inv = (Float)1 / f; return Vector3<T>(x * inv, y * inv, z * inv); } Vector3<T> &operator/=(T f) { Assert(f != 0); Float inv = (Float)1 / f; x *= inv; y *= inv; z *= inv; return *this; }

The Vector3 class also provides a unary negation operator that returns a new vector pointing in the opposite direction of the original one:

<<Vector3 Public Methods>>+=

Vector3<T> operator-() const { return Vector3<T>(-x, -y, -z); }

Finally, Abs() returns a vector with the absolute value operation applied to its components.

<<Geometry Inline Functions>>+=

template <typename T> Vector3<T> Abs(const Vector3<T> &v) { return Vector3<T>(std::abs(v.x), std::abs(v.y), std::abs(v.z)); }

2.2.1 Dot and Cross Product

Two useful operations on vectors are the dot product (also known as the scalar or inner product) and the cross product. For two vectors and , their dot product is defined as:

<<Geometry Inline Functions>>+=

template <typename T> inline T Dot(const Vector3<T> &v1, const Vector3<T> &v2) { return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z; }

The dot product has a simple relationship to the angle between the two vectors:

(2.1)

where is the angle between and , and denotes the length of the vector . It follows from this that is zero if and only if and are perpendicular, provided that neither nor is degenerate—equal to . A set of two or more mutually perpendicular vectors is said to be orthogonal. An orthogonal set of unit vectors is called orthonormal.

It immediately follows from Equation (2.1) that if and are unit vectors, their dot product is the cosine of the angle between them. As the cosine of the angle between two vectors often needs to be computed for rendering, we will frequently make use of this property. A few basic properties directly follow from the definition. For example, if , , and are vectors and is a scalar value, then:

We will frequently need to compute the absolute value of the dot product as well. The AbsDot() function does this for us so that a separate call to std::abs() isn’t necessary.

<<Geometry Inline Functions>>+=

template <typename T> inline T AbsDot(const Vector3<T> &v1, const Vector3<T> &v2) { return std::abs(Dot(v1, v2)); }

The cross product is another useful operation for 3D vectors. Given two vectors in 3D, the cross product is a vector that is perpendicular to both of them. Given orthogonal vectors and , then is defined to be a vector such that form an orthogonal coordinate system.

The cross product is defined as:

A way to remember this is to compute the determinant of the matrix:

where , , and represent the axes , , and , respectively. Note that this equation is merely a memory aid and not a rigorous mathematical construction, since the matrix entries are a mix of scalars and vectors.

In the implementation here, the vector elements are converted to double-precision (regardless of the type of Float) before the subtractions in the Cross() function. Using extra precision for 32-bit floating-point values here protects against error from catastrophic cancellation, a type of floating-point error that can happen when subtracting two values that are very close together. This isn’t a theoretical concern: this change was necessary to fix bugs that came up from this issue previously. See Section 3.9 for more information on floating-point rounding error.

<<Geometry Inline Functions>>+=

template <typename T> inline Vector3<T> Cross(const Vector3<T> &v1, const Vector3<T> &v2) { double v1x = v1.x, v1y = v1.y, v1z = v1.z; double v2x = v2.x, v2y = v2.y, v2z = v2.z; return Vector3<T>((v1y * v2z) - (v1z * v2y), (v1z * v2x) - (v1x * v2z), (v1x * v2y) - (v1y * v2x)); }

From the definition of the cross product, we can derive

(2.2)

where is the angle between and . An important implication of this is that the cross product of two perpendicular unit vectors is itself a unit vector. Note also that the result of the cross product is a degenerate vector if and are parallel.

This definition also shows a convenient way to compute the area of a parallelogram (Figure 2.5). If the two edges of the parallelogram are given by vectors and , and it has height , the area is given by . Since , we can use Equation (2.2) to see that the area is .

Figure 2.5: The area of a parallelogram with edges given by vectors and is equal to . From Equation (2.2), the length of the cross product of and is equal to the product of the two vector lengths times the sine of the angle between them—the parallelogram area.

2.2.2 Normalization

It is often necessary to normalize a vector—that is, to compute a new vector pointing in the same direction but with unit length. A normalized vector is often called a unit vector. The notation used in this book for normalized vectors is that is the normalized version of . To normalize a vector, it’s first useful to be able to compute its length.

<<Vector3 Public Methods>>+=

Float LengthSquared() const { return x * x + y * y + z * z; } Float Length() const { return std::sqrt(LengthSquared()); }

Normalize() normalizes a vector. It divides each component by the length of the vector, . It returns a new vector; it does not normalize the vector in place:

<<Geometry Inline Functions>>+=

template <typename T> inline Vector3<T> Normalize(const Vector3<T> &v) { return v / v.Length(); }

2.2.3 Miscellaneous Operations

A few additional operations are useful when working with vectors. The MinComponent() and MaxComponent() methods return the smallest and largest coordinate value, respectively.

<<Geometry Inline Functions>>+=

template <typename T> T MinComponent(const Vector3<T> &v) { return std::min(v.x, std::min(v.y, v.z)); } template <typename T> T MaxComponent(const Vector3<T> &v) { return std::max(v.x, std::max(v.y, v.z)); }

Related, MaxDimension() returns the index of the component with the largest value.

<<Geometry Inline Functions>>+=

template <typename T> int MaxDimension(const Vector3<T> &v) { return (v.x > v.y) ? ((v.x > v.z) ? 0 : 2) : ((v.y > v.z) ? 1 : 2); }

Component-wise minimum and maximum operations are also available.

<<Geometry Inline Functions>>+=

template <typename T> Vector3<T> Min(const Vector3<T> &p1, const Vector3<T> &p2) { return Vector3<T>(std::min(p1.x, p2.x), std::min(p1.y, p2.y), std::min(p1.z, p2.z)); } template <typename T> Vector3<T> Max(const Vector3<T> &p1, const Vector3<T> &p2) { return Vector3<T>(std::max(p1.x, p2.x), std::max(p1.y, p2.y), std::max(p1.z, p2.z)); }

Finally, Permute() permutes the coordinate values according to the index values provided.

<<Geometry Inline Functions>>+=

template <typename T> Vector3<T> Permute(const Vector3<T> &v, int x, int y, int z) { return Vector3<T>(v[x], v[y], v[z]); }

2.2.4 Coordinate System from a Vector

We will frequently want to construct a local coordinate system given only a single 3D vector. Because the cross product of two vectors is orthogonal to both, we can apply the cross product two times to get a set of three orthogonal vectors for the coordinate system. Note that the two vectors generated by this technique are unique only up to a rotation about the given vector.

The implementation of this function assumes that the vector passed in, v1, has already been normalized. It first constructs a perpendicular vector by zeroing one of the components of the original vector, swapping the remaining two, and negating one of them. Inspection of the two cases should make clear that v2 will be normalized and that the dot product must be equal to zero. Given these two perpendicular vectors, a single cross product gives the third, which by definition will be perpendicular to the first two.

<<Geometry Inline Functions>>+=

template <typename T> inline void CoordinateSystem(const Vector3<T> &v1, Vector3<T> *v2, Vector3<T> *v3) { if (std::abs(v1.x) > std::abs(v1.y)) *v2 = Vector3<T>(-v1.z, 0, v1.x) / std::sqrt(v1.x * v1.x + v1.z * v1.z); else *v2 = Vector3<T>(0, v1.z, -v1.y) / std::sqrt(v1.y * v1.y + v1.z * v1.z); *v3 = Cross(v1, *v2); }