Foundations of Game Engine Development Volume 1 Mathematics (Eric Lengyel 著)

http://www.foundationsofgameenginedev.com/

Chapter1 Vectors and Matrices (已看)

Chapter2 Transforms (已看)

Chapter3 Geometry

Chapter4 Adavanced Algebra

 

Chapter1 Vectors and Matrices

　　1.1 Vector Fundamentals

A scalar is a quantity such as distance,mass,or time that can be fully described using a single numerical value representing its size,or its magnitude.

A vector is a quantity that carries enough information to represent a direction in space in addition to a magnitude

An n-dimensional vector v can be written as　　 v = (v₀, v₁, ..., v_{n - 1});

Zero-based indices are a much better fit for the way in which computers access individual fields in data structures.

The meaning of a vector's components depends on the coordinate system in which those components are expressed.　　v = (v_x, v_y, v_z)

　　1.2 Basic Vecotr Operations

struct Vector3D {
    float x, y, z;

    Vector3D() = default;

    Vector3D(float a, float b, float c) {
        x = a;
        y = b;
        z = c;
    }

    float & operator [](int i) {
        return ((&x)[i]);
    }

    const float & operator [](int i) const {
        return ((&x)[i]);
    }
};

A vector can be used to represent a point that does have a location in space by thinking of the vector as a relative offset from a given origin.　　　　

　　　　1.2.1 Magnitude and Scalar Multiplication

The magnitude of an n-dimensional vector is calculated with the formula 

The vector whose components are all zero called the zero vector, and it is the only vector for which the magnitude is zero.

The magnitude of a vector can be changed by multiplying it by a scalar value.　　tv = (tv₀, tv₁, ..., tv_n-1)　　

A vector that has a magnitude of one is called a unit vector.Unit vector are particularly imortant because they are able to provide directional information without a magnitude when a meaningful size of some kind is not necessary.

The process of setting a vector's magnitude to one is called normalization, and a unit vector that has been produced by this process is often referred to as normalized.

struct Vector3D {
    float x, y, z;

private:
    Vector3D() = default;
public:
    Vector3D(float a, float b, float c) {
        x = a;
        y = b;
        z = c;
    }

    float & operator[](int i) {
        return ((&x)[i]);
    }

    const float & operator[](int i) const {
        return ((&x)[i]);
    }

    Vector3D & operator*=(float s) {
        x *= s;
        y *= s;
        z *= s;
        return (*this);
    }

    Vector3D & operator/=(float s) {
        s = 1.0f / s;
        x *= s;
        y *= s;
        z *= s;
        return (*this);
    }
};

inline Vector3D operator*(const Vector3D & v, float s) {
    return (Vector3D(v.x * s, v.y * s, v.z * s));
}

inline Vector3D operator/(const Vector3D & v, float s) {
    s = 1.0f / s;
    return (Vector3D(v.x * s, v.y * s, v.z  * s));
}

inline Vector3D operator-(const Vector3D & v) {
    return (Vector3D(-v.x, -v.y, -v.z));
}

inline float Magnitude(const Vector3D & v) {
    return (sqrt(v.x * v.x + v.y * v.y + v.z * v.z));
}

inline Vector3D Normalize(const Vector3D & v) {
    return (v / Magnitude(v));
}

　　　　1.2.2 Addtion and Subtraction

Vectors can be added and subtracted by applying these operations componentwise. That is, for two n-dimensional vectors a and b, we have  and 

#ifndef VECTOR3D_H_
#define VECTOR3D_H_

struct Vector3D {
    float x, y, z;

private:
    Vector3D() = default;
public:
    Vector3D(float a, float b, float c) {
        x = a;
        y = b;
        z = c;
    }

    float & operator[](int i) {
        return ((&x)[i]);
    }

    const float & operator[](int i) const {
        return ((&x)[i]);
    }

    Vector3D & operator*=(float s) {
        x *= s;
        y *= s;
        z *= s;
        return (*this);
    }

    Vector3D & operator/=(float s) {
        s = 1.0f / s;
        x *= s;
        y *= s;
        z *= s;
        return (*this);
    }

    Vector3D & operator+=(const Vector3D & v) {
        x += v.x;
        y += v.y;
        z += v.z;
        return (*this);
    }

    Vector3D & operator-=(const Vector3D & v) {
        x -= v.x;
        y -= v.y;
        z -= v.z;
        return (*this);
    }
};

inline Vector3D operator*(const Vector3D & v, float s) {
    return (Vector3D(v.x * s, v.y * s, v.z * s));
}

inline Vector3D operator/(const Vector3D & v, float s) {
    s = 1.0f / s;
    return (Vector3D(v.x * s, v.y * s, v.z  * s));
}

inline Vector3D operator-(const Vector3D & v) {
    return (Vector3D(-v.x, -v.y, -v.z));
}

inline float Magnitude(const Vector3D & v) {
    return (sqrt(v.x * v.x + v.y * v.y + v.z * v.z));
}

inline Vector3D Normalize(const Vector3D & v) {
    return (v / Magnitude(v));
}

inline Vector3D operator+(const Vector3D & a, const Vector3D & b) {
    return (Vector3D(a.x + b.x, a.y + b.y, a.z + b.z));
}

inline Vector3D operator-(const Vector3D & a, const Vector3D & b) {
    return (Vector3D(a.x - b.x, a.y - b.y, a.z - b.z));
}

#endif

　　1.3 Matrix Fundamentals

A matrix is a mathematical object composed of a set of numerical quantities arranged in a two-dimensional array of rows and columns. When a matrix has n rows and m columns, we say that its size is n x m,which is read "n by m". If it's the case that n = m,

then we say the matrix is "square matrix". The numbers that make up a matrix M are called its entries.

The symbol M_ij means the entry apearing int he i-th row and j-th column.Sometimes, a comma may be inserted between the indices for clarity.Using this notation, an n x m matrix M can be written as 

The entries M_ij , where the row index and the column index are equal to each other, are called the diagonal entries of the matrix M, and they are said to reside on the main diagonal of the matrix.

The main diagonal starts with the upper-left entry and continues downward and to the right.

The entries M_ij  for which i != j are called the off-diagonal entries of the matrix M.

Any matrix for which all of the off-diagonal entries are zero is called a diagnoal matrix. Note that some or all of the entries on the main diagonal itself could be zero,and the matrix would still be considered to be a diagonal matrix.

The transpose of a matrix M is the matrix denoted by M^T whose rows are equal to the columns of M, or equivalently, whose columns are equal to the rows of M.

If the matrix M has size n x m, then the matrix M^T has size m x n, and its entries are given by M_ij^T = M_ji.

The transpose of a matrix can be thought of as the reflection of its entries across the main diagonal.

If a matrix M is equal to its transpose M^T, meaning that it's always the case that M_ij = M_ji, then it is called a symmetric matrix because all of the entries above and right of the main diagonal are the same as the entries below and left of the main diagonal, 

with the row and column indices reversed.

A matrix must be square in order to be symmetric, and every diagoal matrix is automatically symmetric.

If the entries of a transpose M^T are equal to the negations of the same entries in the matrix M,meaning that it's always the case that M_ij^T = -M_ij , then the matrix M is called an antisymmetric matrix(反对称矩阵) or a skew-symmetric matrix(斜对称矩阵).

Note that for this to be the case, all of the entries on the main diagonal must be zero as in the example.

An n-dimensional vector can be regarded as an n x 1 matrix or as 1 x n matrix,and as such, is called a column vector or a row vector, respectively.

A comma-spearated list of n components is equivalent to an n x 1 matrix containing the same numbers in the same order.

If we apply the transpose operation to v, then we get the row vector  which despite its horizontal layout, is different from the comma-separated list of components.

It is frequently useful to treat a matrix as an array of column vectors or row vectors.For example,suppose that a,b, and c are column vectors. Then we can construct 3 x 3 matrix M by making those vectors the columns of the matrix and writing it as 

#ifndef MATRIX3D_H_
#define MATRIX3D_H_

#include "Vector3D.h"

struct Matrix3D {
private:
    float n[3][3];
    Matrix3D() = default;
public:
    // column major order
    Matrix3D(float n00, float n01, float n02,
        float n10, float n11, float n12,
        float n20, float n21, float n22) {
        n[0][0] = n00; n[0][1] = n10; n[0][2] = n20;
        n[1][0] = n01; n[1][1] = n11; n[1][2] = n21;
        n[2][0] = n02; n[2][1] = n12; n[2][2] = n22;
    }

    Matrix3D(const Vector3D & a, const Vector3D & b, const Vector3D & c) {
        n[0][0] = a.x; n[0][1] = a.y; n[0][2] = a.z;
        n[1][0] = b.x; n[1][1] = b.y; n[1][2] = b.z;
        n[2][0] = c.x; n[2][1] = c.y; n[2][2] = c.z;
    }

    float & operator()(int i, int j) {
        return (n[j][i]);
    }

    const float & operator()(int i, int j) const {
        return (n[j][i]);
    }

    Vector3D & operator[](int j) {
        return (*reinterpret_cast<Vector3D *>(n[j]));
    }

    const Vector3D & operator[](int j) const {
        return (*reinterpret_cast<const Vector3D *>(n[j]));
    }
};

#endif

　　1.4 Basic Matrix Operations

　　　　1.4.1 Addition,Subtraction,and Scalar Multiplication

Two matrices of the same size can be added or subtracted by simply adding or subtracting corresponding entries.

A matrix M is multiplied by a scalar t by applying the multiplication to every entry of the matrix.

　　　　1.4.2 Matrix Multiplication

Two matrices can be multiplied together if and only if the number of columns in the first matrix is equal to the number of rows in the second matrix.

The result is a new matrix having the same number of rows as the first matrix and the same number of columns as the second matrix.

When an n x p matrix A is multiplied by a p x m matrix B, the (i,j) entry of the matrix product AB is given by the formula

#ifndef MATRIX3D_H_
#define MATRIX3D_H_

#include "Vector3D.h"

struct Matrix3D {
private:
    float n[3][3];
    Matrix3D() = default;
public:
    // column major order
    Matrix3D(float n00, float n01, float n02,
        float n10, float n11, float n12,
        float n20, float n21, float n22) {
        n[0][0] = n00; n[0][1] = n10; n[0][2] = n20;
        n[1][0] = n01; n[1][1] = n11; n[1][2] = n21;
        n[2][0] = n02; n[2][1] = n12; n[2][2] = n22;
    }

    Matrix3D(const Vector3D & a, const Vector3D & b, const Vector3D & c) {
        n[0][0] = a.x; n[0][1] = a.y; n[0][2] = a.z;
        n[1][0] = b.x; n[1][1] = b.y; n[1][2] = b.z;
        n[2][0] = c.x; n[2][1] = c.y; n[2][2] = c.z;
    }

    float & operator()(int i, int j) {
        return (n[j][i]);
    }

    const float & operator()(int i, int j) const {
        return (n[j][i]);
    }

    Vector3D & operator[](int j) {
        return (*reinterpret_cast<Vector3D *>(n[j]));
    }

    const Vector3D & operator[](int j) const {
        return (*reinterpret_cast<const Vector3D *>(n[j]));
    }
};

Matrix3D operator*(const Matrix3D & A, const Matrix3D & B) {
    return (Matrix3D(
        A(0, 0) * B(0, 0) + A(0, 1) * B(1, 0) + A(0, 2) * B(2, 0),
        A(0, 0) * B(0, 1) + A(0, 1) * B(1, 1) + A(0, 2) * B(2, 1),
        A(0, 0) * B(0, 2) + A(0, 1) * B(1, 2) + A(0, 2) * B(2, 2),

        A(1, 0) * B(0, 0) + A(1, 1) * B(1, 0) + A(1, 2) * B(2, 0),
        A(1, 0) * B(0, 1) + A(1, 1) * B(1, 1) + A(1, 2) * B(2, 1),
        A(1, 0) * B(0, 2) + A(1, 1) * B(1, 2) + A(1, 2) * B(2, 2),

        A(2, 0) * B(0, 0) + A(2, 1) * B(1, 0) + A(2, 2) * B(2, 0),
        A(2, 0) * B(0, 1) + A(2, 1) * B(1, 1) + A(2, 2) * B(2, 1),
        A(2, 0) * B(0, 2) + A(2, 1) * B(1, 2) + A(2, 2) * B(2, 2)
    ));
}

Vector3D operator*(const Matrix3D & M, const Vector3D & v) {
    return (Vector3D(
        M(0, 0) * v.x + M(0, 1) * v.y + M(0, 2) * v.z,
        M(1, 0) * v.x + M(1, 1) * v.y + M(1, 2) * v.z,
        M(2, 0) * v.x + M(2, 1) * v.y + M(2, 2) * v.z
    ));

}

#endif

 

　　1.5 Vector Multiplication

　　　　1.5.1 Dot Product

The dot product between two n-dimensional vectors a and b is a scalar quantity given by the formula 

If the vectors a and b are regarded as column matrices, then the dot product be also be expressed as which produces a 1 x 1 matrix having a single entry.

Although not at all obvious from its definition, the dot product between two vectors a and b statisfies the equality  where θ is the planar angle between the directions of a and b if the were to be drawn as arrows starting at the same location

This equality represent the main application of the dot product, and it provides a computationally cheap way to determine how much two vectors are parallel to each other or perpendicular to each other.

If a and b are both unit vectors,then a · b is always in the range [-1,1] because in this case a · b = cosθ, and the range of the cosine function is [-1, 1]

Assuming the magnitude of a and b remain the same, the dot product a · b attains its largest positive value when a and b are parallel and point in the same direction.

When a and b are parallel and point in opposite directions, a · b attains its largest negative value. If a and b are perpendicular, then a · b is zero no matter what the magnitudes of a and b are.

In general, the dot product is positive when the angle between the vectors is less then 90 degrees and negative when the angle is greater than 90 degrees.Loosely speaking, the dot product provides a measure of how much one vector is like another.

When a · b = 0, the vectors a and b are said to be orthogonal, and this term is used even if one of the vectors being multiplied together is the zero vector.

Orthogonality is a concept that includes the geometric state of two vectors being perpendicular to each other, but it also has a more abstract meaning in different settings that are not explored in this book

　　　　1.5.2 Cross Product

#ifndef VECTOR3D_H_
#define VECTOR3D_H_

struct Vector3D {
    float x, y, z;

private:
    Vector3D() = default;
public:
    Vector3D(float a, float b, float c) {
        x = a;
        y = b;
        z = c;
    }

    float & operator[](int i) {
        return ((&x)[i]);
    }

    const float & operator[](int i) const {
        return ((&x)[i]);
    }

    Vector3D & operator*=(float s) {
        x *= s;
        y *= s;
        z *= s;
        return (*this);
    }

    Vector3D & operator/=(float s) {
        s = 1.0f / s;
        x *= s;
        y *= s;
        z *= s;
        return (*this);
    }

    Vector3D & operator+=(const Vector3D & v) {
        x += v.x;
        y += v.y;
        z += v.z;
        return (*this);
    }

    Vector3D & operator-=(const Vector3D & v) {
        x -= v.x;
        y -= v.y;
        z -= v.z;
        return (*this);
    }
};

inline Vector3D operator*(const Vector3D & v, float s) {
    return (Vector3D(v.x * s, v.y * s, v.z * s));
}

inline Vector3D operator/(const Vector3D & v, float s) {
    s = 1.0f / s;
    return (Vector3D(v.x * s, v.y * s, v.z  * s));
}

inline Vector3D operator-(const Vector3D & v) {
    return (Vector3D(-v.x, -v.y, -v.z));
}

inline float Magnitude(const Vector3D & v) {
    return (sqrt(v.x * v.x + v.y * v.y + v.z * v.z));
}

inline Vector3D Normalize(const Vector3D & v) {
    return (v / Magnitude(v));
}

inline Vector3D operator+(const Vector3D & a, const Vector3D & b) {
    return (Vector3D(a.x + b.x, a.y + b.y, a.z + b.z));
}

inline Vector3D operator-(const Vector3D & a, const Vector3D & b) {
    return (Vector3D(a.x - b.x, a.y - b.y, a.z - b.z));
}

inline float Dot(const Vector3D & a, const Vector3D & b) {
    return (a.x * b.x + a.y * b.y + a.z * b.z);
}

inline Vector3D Cross(const Vector3D & a, const Vector3D & b) {
    return (Vector3D(
        a.y * b.z - a.z * b.y,
        a.z * b.x - a.x * b.z,
        a.x * b.y - a.y * b.x
    ));
}

#endif

The cross product between two 3D vectors a and b is another 3D vector given by the formula 

The cross product can also be expressed as a matrix product by forming a special 3 x 3 antisymmetric matrix denoted by [a]_x and multiplying by the column vector b. The matrix [a]_x is defined as 

and when multiplied by b, it gives us 

It's important to emphasize that the cross product is defined only for three dimensions, whereas the dot product is defined for all numbers of dimensions.

This limitation is a consequence of the fact that the cross product is actually a subtle misinterpretation of a more general and more algebraically sound operation called the wedge product(楔积). 

If two vectors a and b are parallel, either because they point the same direction or they point in opposite directions,then the cross product a x bis zero no matter what the magnitudes of a and b are. 

Because the are parallel, one of the vectors can be written as a scalar multiple of the other, so we can say that b = ta for some scalar t.

The fact that the cross product is zero then becomes obvious when b is replaced by ta in the definition to get 

When two vectors a and b are not parallel, the cross product a x b is a new vector that is perpendicular to both a and b.This is evident if we calculate the dot products a · (a x b) and b · (a x b) because both of them are always zero

The vectors a and b, not being parallel, establish a plane to which the cross product a x b is perpendicular, and we have two choices as to what direction the vector a x b actually points.

If we are looking at the plane from a position not lying in the plane itself, then a x b could point toward us, or it could point away from us.

The correct direction is determined by the handeness of the coordinates system.

The magnitude of the cross product between two vectors a and b satisfies the equality  where θ is the planar angle between the directions of a and b if they were to be drawn as arrows starting at the same location.

　　　　1.5.3 Scalar Triple Product

The scalar triple product of three vectors a,b and c is the scalar quantity produced by multiplying two of the vectors together using the cross product and then multiplying by the third vector using the dot product, as in (a x b) · c.

It turns out that it doesn't matter where the cross product goes and where the dot product gose, and the scalar triple product gives the same result as long as the vector multiplied together in the same order with wraparound.

As such, the special notation [a, b, c], without any specific multiplication symbols, is often used to represent the scalar triple product, and we can write 

If the order of the input vectors is reversed, then the scalr tripple product is negated to give us 

This accounts for all six possible permutation of the vectors a, b and c.

　　1.6 Vector Projection

Given a particular vector, we might want to find two or more other vectors with specific alignments that add up to our original vector. This is a process called decomposing a vector into its separate components.

The most straightforward decomposition invovles breaking a vector into pieces that are aligned to the coordinates axes.

The letters i, j and k are commonly used to represent unit vector aligned to the positive x,y and z axes,and they are thus defined as 

If we wanted to decompose a 3D vector v into components parallel to the coordinate axes, then we could write 

In general, we can use the dot product to project any vector a onto any other nonzero vector b using the formula 

The notation a_||b indicates the component of the vector a that is parallel to the vector b,and equation gives us the projection of a onto b(The alternative notation proj_ba is sometimes used in other texts for the projection of a onto b)

The projection of a onto b can be expressed as the matrix product

The product bb^T yields a symmetric matrix that can be multiplied by the vector a to perform the projection. In three dimensions, we have 

The matrix in this equation is an example of an outer product. Ingeneral, the outer product between two vectors u and v is written as , and it produces a matrix for which the (i,j) entry is equal to u_iv_j as in 

The perpendicular part of the decomposition is called the rejection of a from b and is written

#ifndef VECTOR3D_H_
#define VECTOR3D_H_

struct Vector3D {
    float x, y, z;

private:
    Vector3D() = default;
public:
    Vector3D(float a, float b, float c) {
        x = a;
        y = b;
        z = c;
    }

    float & operator[](int i) {
        return ((&x)[i]);
    }

    const float & operator[](int i) const {
        return ((&x)[i]);
    }

    Vector3D & operator*=(float s) {
        x *= s;
        y *= s;
        z *= s;
        return (*this);
    }

    Vector3D & operator/=(float s) {
        s = 1.0f / s;
        x *= s;
        y *= s;
        z *= s;
        return (*this);
    }

    Vector3D & operator+=(const Vector3D & v) {
        x += v.x;
        y += v.y;
        z += v.z;
        return (*this);
    }

    Vector3D & operator-=(const Vector3D & v) {
        x -= v.x;
        y -= v.y;
        z -= v.z;
        return (*this);
    }
};

inline Vector3D operator*(const Vector3D & v, float s) {
    return (Vector3D(v.x * s, v.y * s, v.z * s));
}

inline Vector3D operator/(const Vector3D & v, float s) {
    s = 1.0f / s;
    return (Vector3D(v.x * s, v.y * s, v.z  * s));
}

inline Vector3D operator-(const Vector3D & v) {
    return (Vector3D(-v.x, -v.y, -v.z));
}

inline float Magnitude(const Vector3D & v) {
    return (sqrt(v.x * v.x + v.y * v.y + v.z * v.z));
}

inline Vector3D Normalize(const Vector3D & v) {
    return (v / Magnitude(v));
}

inline Vector3D operator+(const Vector3D & a, const Vector3D & b) {
    return (Vector3D(a.x + b.x, a.y + b.y, a.z + b.z));
}

inline Vector3D operator-(const Vector3D & a, const Vector3D & b) {
    return (Vector3D(a.x - b.x, a.y - b.y, a.z - b.z));
}

inline float Dot(const Vector3D & a, const Vector3D & b) {
    return (a.x * b.x + a.y * b.y + a.z * b.z);
}

inline Vector3D Cross(const Vector3D & a, const Vector3D & b) {
    return (Vector3D(
        a.y * b.z - a.z * b.y,
        a.z * b.x - a.x * b.z,
        a.x * b.y - a.y * b.x
    ));
}

inline Vector3D Project(const Vector3D & a, const Vector3D & b) {
    return (b * (Dot(a, b) / Dot(b, b)));
}

inline Vector3D Reject(const Vector3D & a, const Vector3D & b) {
    return (a - b * (Dot(a, b) / Dot(b, b)));
}

#endif

　　1.7 Matrix Inversion

One of the main reasons that matrices appear so frequently in game engines is because they represent coordinate system transformations.

That is, a matrix M describes how a vector, point, line, plane, or even another transformation can be moved from a coordinate system A with its own origin and set of axes to another coordinate system B with a different origin and set of axes.

It is often necessary to be able to perform the reverse transformation from coordinate system B back to coordinate system A, and to accomplish this, we must be able to find a matrix M^-1, called  the inverse of M, that undoes the transformation performed by the matrix M

　　　　1.7.1 Identity Matrices

The identity matrix I_n is the n x nmatrix whose entries on the main diagonal are all ones and whose entries everywhere else are all zero, which can be written as 

Inverse are defined onlly for square matrices, and the inverse of an n x n matrix M is a matrix denoted by M^-1 having the property that MM^-1 = I_n and M^-1M = I_n. An inverse does not always exist.

　　　　1.7.2 Determinants

The determinant of an n x n matrix M is a scalar value that can be thought of as a sort of magnitude M.It is written as det(M) or |M|.

A matrix has an inverse if and only if its determinant is non zero

The calculation of the determinant for an n x n matrix M requires that we sum over n! terms corresponding to all of the possible permutations of the sequence {0, 1, 2, ..., n - 1}

The determinant of  3 x 3 matrix M is given by 

The determinant of a 1 x 1 matrix is simply the value of its single entry.

The determinant of a 2 x 2 matrix 

#ifndef MATRIX3D_H_
#define MATRIX3D_H_

#include "Vector3D.h"

struct Matrix3D {
private:
    float n[3][3];
    Matrix3D() = default;
public:
    // column major order
    Matrix3D(float n00, float n01, float n02,
        float n10, float n11, float n12,
        float n20, float n21, float n22) {
        n[0][0] = n00; n[0][1] = n10; n[0][2] = n20;
        n[1][0] = n01; n[1][1] = n11; n[1][2] = n21;
        n[2][0] = n02; n[2][1] = n12; n[2][2] = n22;
    }

    Matrix3D(const Vector3D & a, const Vector3D & b, const Vector3D & c) {
        n[0][0] = a.x; n[0][1] = a.y; n[0][2] = a.z;
        n[1][0] = b.x; n[1][1] = b.y; n[1][2] = b.z;
        n[2][0] = c.x; n[2][1] = c.y; n[2][2] = c.z;
    }

    float & operator()(int i, int j) {
        return (n[j][i]);
    }

    const float & operator()(int i, int j) const {
        return (n[j][i]);
    }

    Vector3D & operator[](int j) {
        return (*reinterpret_cast<Vector3D *>(n[j]));
    }

    const Vector3D & operator[](int j) const {
        return (*reinterpret_cast<const Vector3D *>(n[j]));
    }
};

Matrix3D operator*(const Matrix3D & A, const Matrix3D & B) {
    return (Matrix3D(
        A(0, 0) * B(0, 0) + A(0, 1) * B(1, 0) + A(0, 2) * B(2, 0),
        A(0, 0) * B(0, 1) + A(0, 1) * B(1, 1) + A(0, 2) * B(2, 1),
        A(0, 0) * B(0, 2) + A(0, 1) * B(1, 2) + A(0, 2) * B(2, 2),

        A(1, 0) * B(0, 0) + A(1, 1) * B(1, 0) + A(1, 2) * B(2, 0),
        A(1, 0) * B(0, 1) + A(1, 1) * B(1, 1) + A(1, 2) * B(2, 1),
        A(1, 0) * B(0, 2) + A(1, 1) * B(1, 2) + A(1, 2) * B(2, 2),

        A(2, 0) * B(0, 0) + A(2, 1) * B(1, 0) + A(2, 2) * B(2, 0),
        A(2, 0) * B(0, 1) + A(2, 1) * B(1, 1) + A(2, 2) * B(2, 1),
        A(2, 0) * B(0, 2) + A(2, 1) * B(1, 2) + A(2, 2) * B(2, 2)
    ));
}

Vector3D operator*(const Matrix3D & M, const Vector3D & v) {
    return (Vector3D(
        M(0, 0) * v.x + M(0, 1) * v.y + M(0, 2) * v.z,
        M(1, 0) * v.x + M(1, 1) * v.y + M(1, 2) * v.z,
        M(2, 0) * v.x + M(2, 1) * v.y + M(2, 2) * v.z
    ));

}

float Determinant(const Matrix3D & M) {
    return (
        M(0, 0) * M(1, 1) * M(2, 2) - M(0, 2) * M(1, 1) * M(2, 0) +
        M(0, 1) * M(1, 2) * M(2, 0) - M(0, 1) * M(1, 0) * M(2, 2) +
        M(0, 2) * M(1, 0) * M(2, 1) - M(0, 0) * M(1, 2) * M(2, 1));
}

#endif

The determinant of any n xn matrix M can be expressed as 

Using expansion by minors, the determinant  of an n x n matrix M is given by 

　　　　1.7.3 Elementary Matrices 

The inverse of a matrix M can be found by systematically performing a sequence of basic operations on M until it has been transformed into the identity matrix.

Each of these basic operations can be represented by an elementary matrix(初等矩阵), and the product of all the elementary matrices corresponding to all the basic operations used to transform M into the identity provides us with the inverse of M

There are three elementary row operations(初等行运算), named as such because they each affect one or two entiere rows of a matrix, and the particular types perform the following modifications to a matrix M.

　　(a) Multiply one row of M by a nonzero scalar value

　　(b) Exchange two rows of M

　　(c) Add a scalar multiple of one row of M to another row of M

　　　　1.7.4 Inverse Calculation

　　　　1.7.5 Inverse of Small Matrices

#ifndef MATRIX3D_H_
#define MATRIX3D_H_

#include "Vector3D.h"

struct Matrix3D {
private:
    float n[3][3];
    Matrix3D() = default;
public:
    // column major order
    Matrix3D(float n00, float n01, float n02,
        float n10, float n11, float n12,
        float n20, float n21, float n22) {
        n[0][0] = n00; n[0][1] = n10; n[0][2] = n20;
        n[1][0] = n01; n[1][1] = n11; n[1][2] = n21;
        n[2][0] = n02; n[2][1] = n12; n[2][2] = n22;
    }

    Matrix3D(const Vector3D & a, const Vector3D & b, const Vector3D & c) {
        n[0][0] = a.x; n[0][1] = a.y; n[0][2] = a.z;
        n[1][0] = b.x; n[1][1] = b.y; n[1][2] = b.z;
        n[2][0] = c.x; n[2][1] = c.y; n[2][2] = c.z;
    }

    float & operator()(int i, int j) {
        return (n[j][i]);
    }

    const float & operator()(int i, int j) const {
        return (n[j][i]);
    }

    Vector3D & operator[](int j) {
        return (*reinterpret_cast<Vector3D *>(n[j]));
    }

    const Vector3D & operator[](int j) const {
        return (*reinterpret_cast<const Vector3D *>(n[j]));
    }
};

Matrix3D operator*(const Matrix3D & A, const Matrix3D & B) {
    return (Matrix3D(
        A(0, 0) * B(0, 0) + A(0, 1) * B(1, 0) + A(0, 2) * B(2, 0),
        A(0, 0) * B(0, 1) + A(0, 1) * B(1, 1) + A(0, 2) * B(2, 1),
        A(0, 0) * B(0, 2) + A(0, 1) * B(1, 2) + A(0, 2) * B(2, 2),

        A(1, 0) * B(0, 0) + A(1, 1) * B(1, 0) + A(1, 2) * B(2, 0),
        A(1, 0) * B(0, 1) + A(1, 1) * B(1, 1) + A(1, 2) * B(2, 1),
        A(1, 0) * B(0, 2) + A(1, 1) * B(1, 2) + A(1, 2) * B(2, 2),

        A(2, 0) * B(0, 0) + A(2, 1) * B(1, 0) + A(2, 2) * B(2, 0),
        A(2, 0) * B(0, 1) + A(2, 1) * B(1, 1) + A(2, 2) * B(2, 1),
        A(2, 0) * B(0, 2) + A(2, 1) * B(1, 2) + A(2, 2) * B(2, 2)
    ));
}

Vector3D operator*(const Matrix3D & M, const Vector3D & v) {
    return (Vector3D(
        M(0, 0) * v.x + M(0, 1) * v.y + M(0, 2) * v.z,
        M(1, 0) * v.x + M(1, 1) * v.y + M(1, 2) * v.z,
        M(2, 0) * v.x + M(2, 1) * v.y + M(2, 2) * v.z
    ));

}

float Determinant(const Matrix3D & M) {
    return (
        M(0, 0) * M(1, 1) * M(2, 2) - M(0, 2) * M(1, 1) * M(2, 0) +
        M(0, 1) * M(1, 2) * M(2, 0) - M(0, 1) * M(1, 0) * M(2, 2) +
        M(0, 2) * M(1, 0) * M(2, 1) - M(0, 0) * M(1, 2) * M(2, 1));
}

Matrix3D Inverse(const Matrix3D & M) {
    const Vector3D & a = M[0];
    const Vector3D & b = M[1];
    const Vector3D & c = M[2];

    Vector3D r0 = Cross(b, c);
    Vector3D r1 = Cross(c, a);
    Vector3D r2 = Cross(a, b);

    float invDet = 1.0f / Dot(r2, c);
    return (Matrix3D(
        r0.x * invDet, r0.y * invDet, r0.z * invDet,
        r1.x * invDet, r1.y * invDet, r1.z * invDet,
        r2.x * invDet, r2.y * invDet, r2.z * invDet
    ));
}



#endif

　　Exercises for Chapter1

Chapter2 Transforms

A dynamic object in a game may need to move from one point to another, and it may need to rotate itself to different orientations.

A model may be composed of a collection of objects arranged in a tree structure,and each part may move in a manner that is relative to another part above it in the hierarchy.

It may be necessary to express positions and orientations in may different local coordinate systems used by various components of a rendering system.

All of these are examples of cases that require the application of transforms.

　　2.1 Coordinate Spaces

It is typical for a game engine to define a number of different coordinate systems.

There is usually a coordinate system called world space or global space that serves as a fixed background relative to which other coordinate systems can be established.]

Various objects in a game, which can include things like models, light sources, and cameras, often have their own independent coordinate systems called object space or local space.

When an interaction occurs between two objects using different coordinate systems, either one object need to be transformed into the coordinate system used by the other object

or both objects need to be transformed into some other common coordinate system.

　　　　2.1.1 Transformation Matrices

The transformation from a position P_A in coordinate system A to the position P_B in coordinate system B can be expressed as 

where M is a 3 x 3 martrix that reorients the coordinate axes, and t is a 3D translation vector that moves the origin of the coordinate system.

The transformation is called an affine transformation

Assuming that M is invertible, we can solve this equation for P_A to obtain the inverse transformation from coordinate system B to coordinate system A as 

In general, the linear transformation V_B = MV_A replaces the axes in coordinate system A with the columns of the matrix M in coordinate system B.

The vector V_B is then a combination of the new axes in the same way that V_A was a combination of the old axes.

Suppose that the columns of M are given by a, b, and c so that we can write M = [a b c].

Then and for an arbitrary vector v,

　　　　2.1.2 Orthogonal Transforms

Many transformation matrices appearing in game engine do have mutually perpendicular unit-length columns, and such matrices are called orthogonal matrices(正交矩阵).

Orthogonal matrices have a number of interesting properties, the first of which is that the inverse of an orthogonal matrix is equal its transpose.

Assuming that a,b, and c all have unit length and are mutually perpendicular, we can calculate the product M^TM as 

Since a,b, and c each have unit length, the entries along the main diagonal are all ones.

Since a,b, and c are mutually perpendicular, all of the other entries are zero.

Euqation also demonstrates that the reverse implication is true.

If we assume that M^-1 = M^T, then the right side of the equation must be the identity matrix, so a,b, and c must be mutually perpendicular and have unit length.

Additionally, if M is orthogonal, the M^T must also be orthogonal because its inverse is equal to its transpose, which is just M.

This implies that the columns of M^T form a set of mutually perpendicular unit-length vectors, and that is equivalent to making the same statement about the rows of M.

Thus, the following list of statements all have the same meaning, and each one implies all of the others.

　　M is an othogonal matrix.

　　The inverse of M is equal to M^T

　　The columns of M are mutually perpendicular unit-length vectors

　　The rows of M are mutually perpendicular unit-length vectors

When an object is transformed by an orthogonal matrix, it may be reoriented  in space and/or reflected in a mirror, but it still has the exact same shape that it had before the transformation.

Orthogonal matrices preserve the dot product between any two vectors a and b, and this is easily proven by considering the dot product of the two vectors after the are transformed by an orthogonal matrix M as in 

Since a·a gives the squared  magnitude of a, this also proves that the length of a vector is not changed by an orthogonal matrix.

It must therefore be true that the angle θ between a and b is preserved as well because the dot product would otherwise change due to its relationship to the cosine of θ

The transform performed by an orthogonal matrix is always a rotation, a reflection, or a combination of the two.

The determinant of an orthogonal matrix is always +-1,positive for a pure rotation, and negative for a transform that includes a reflection.

Transforms that include a reflection reverse the handeness of the coordinate system, changing a right-handed coordinate system into a left-handed one, and vice versa.

　　　　2.1.3 Transform Composition

It is often the case that multiple transforms need to be applied to an object to either changes its shape in several steps or to express its position and orientation in a common coordinate system that might be several levels away in a hierarchical model.

Whenever a vector v needs to be transformed first by a matrix M₁ and then by another matrix M₂, we calculate the result v' as 

Since matrix multiplication is associative, we can multiply the two transformation matrices together first to create a single combined matrix N = M₂M₁,

and then we can calculate v' with the simpler equation v' = Nv.

This can be continued indefinitely with n transformation matrices, and their product MnMn-1...M2M1 can be precomputed and stored as a single transform

There are times when a transform is expressed in one coordinate system but needs to be applied to an object using a different coordinate system.

Fore example, it's convenient to express a scale transform in a coordinate system where the scale factors apply in directions parallel to the coordinate axes.

If a particular transform A can be applied in coordinate system A, but we have an object using coordinate system B, then the equivalent transform B in cooridnate system B is given by  where the matrix M transforms vectors from coordinate system A to coordinate system B.

When the transform B is applied to a vector v in coordinate system B, you can think of it as first transforming v from coordinate system B to coordinate system A with the matrix M^-1,

applying the transform A in that setting, and then transforming back into coordinate system B with the matrix M.

　　2.2 Rotations

Rotations often occur in a local coordinate system in which the axis of rotation is aligned to one of the coordinate axes, but the can also be applied about an arbitrary axis specified by a direction vector.

When a vector v is rotated about an axis, we can consider whta happens to the components of v that are parallel and perpendicular to the axis.

The component that is parallel to the axis is unaffected while the component that is perpendicular to the axis changes.

In general, a rotation always occurs among the directions parallel to an oriented plane in space, leaving every direction orthogonal to the plane alone.

In three dimensions, there is only one remaining direction, and that is what we call the axis of the rotation, but this concept does not extend to other numbers of dimensions.

　　　　2.2.1 Rotation About a Coordinate Axis

A rotation about the x,y,z axis occurs in the plane formed by the other two axes.

Suppose that we would like to rotate the vector v through an angle θ about the z axis.

Using the unit vectors i,j,k parallel to the coordinate axes, we can write v as 

Since k is parallel to the z axis, the component v_z does not change during the rotation,

It is only the component v_xi + v_yj that changes, and it does so in such a way that the vector v_xi is rotated through the angle θ, the result is a new vector having the same length v_x that can be written as a sum of vectors parallel to i and j, forming a right triangle.

As shown in Figure 2.1, when the vector v_xi is rotated through the angle θ, the result is a new vector having the same length v_x that can be written as a sum of vectors parallel to i and j, forming a right triangle.

Basic trignometry tells us that their lengths are v_xcosθ and v_xsinθ, respectively, so we can write the rotated vector as 

This represents the rotation only of the vector v_xi. 

We must also account for the vector v_yj, and the result of its rotation is similar in a setting that's rotated 90 degrees counterclockwise, so it is written in terms of j and -i instead of i and j.

The rotation of v_yj through the angle θ is 

When we combine equations and include the unaltered component v_zk, we can express the complete formula for the rotated vector v' as 

The transforms M_{rot x}(θ), M_{rot y}(θ),M_{rot z}(θ) that rotate a vector through an angle θ about the x,y,z axes are given by the matrices

Matrix3D MakeRotationX(float t) {
    float c = cos(t);
    float s = sin(t);
    return (Matrix3D(1.0f, 0.0f, 0.0f,
                             0.0f, c, -s,
                             0.0f, s, c));
Matrix3D MakeRotationY(float t) {
    float c = cos(t);
    float s = sin(t);
    return (Matrix3D(c, 0.0f, s,
                             0.0f, 1.0f, 0.0f,
                             -s, 0.0f, c));

Matrix3D MakeRotationZ(float t) {
    float c = cos(t);
    float s = sin(t);
    return (Matrix3D(c, -s, 0.0f,
                             s, c, 0.0f,
                             0.0f, 0.0f, 1.0f));

　　　　2.2.2 Rotation About an Arbitrary Axis

Let α be the angle between the vectors a and v.

The vector (the rejection of v from a) is perpendicular to a, and it has a length equal to ||v||sinα because it forms the side opposite the angle α in the right triangle.

It's this component that we need to rotate in the plane perpendicular  to the axis a.

As before, we can perform this rotation by expressing the result as a linear combination of the original vector and another vector in the plane that is the 90-degree counterclockwise rotation of the original vector.

Fortunately, this second vector is easily obtained as a x v, and it just happens to have the same length, ||v||sinα, as  does.

This means that we can express the rotated vector v' as , where θ is the angle of rotation about the axis a.

The component  v_||a parallel to the axis of rotation does not change, and the component  is replaced by a linear combination of  and a x v.

When we expand the definitions of the projection and rejection,Equation takes the form  where we have omitted the divisions by a² because a is a unit vector. This can be simplified a little bit to obtain

The projection (v·a)a and the cross product a x v can each be expressed as a 3 x 3 matrix multiplying the vector v.

Making these replacements gives uswhere we have also inserted an identity matrix in front of the first term so that all three terms contain a 3 x 3 matrix.

When we combine everything into a single matrix, we get the formula for the transform M_rot (θ,a) that rotates a vector through an angle θ about the axis a,where we have used the abbreviations c = cosθ and s = sinθ.

Matrix3D MakeRotation(float t, const Vector3D & a) {
    float c = cos(t);
    float s = sin(t);
    float d = 1.0f - c;

    float x = a.x * d;
    float y = a.y * d;
    float z = a.z * d;
    
    float axay = x * a.y;
    float axaz = x * a.z;
    float ayaz = y * a.z;
    
    return (Matrix3D(c + x * a.x, axay - s * a.z, axaz + s *a.y,
                             axay + s * a.z, c + y * a.y, ayaz - s * a.x,
                             axaz - s * a.y, ayaz + s * a.x, c + z * a.z));

　　2.3 Reflections

A vector v can be reflected through a plane perpendicular to vector a by decomposing v into its components perpendicular to a and parallel to a and then simply negating the parallel componenet.

The original vector can be written as , and the reflected vector v' is then given by 

When we replace each component by its matrix representation, we get  where we have assumed that a has unit length. By combining these terms into a single matrix, we arrive at the formula for the transform M_reflect(a) that reflects a vector through the plane perpendicular to the unit vector a.

The matrix M_reflect(a) reflects through a plane by negating the component of a vector parallel to a one-dimensional subspace represented by the direction a.

We can also construct a transform that instead negates the perpendicular component and leaves the parallel component unchanged.

Constructing this transform is a simple matter of negating equation to get 

Since the component  represents everything that is perpendicular to a, we are actually negating a two-dimensional subspace by performing two reflections aligned to vectors that are both orthogonal to a and each other.

The composition of two reflection is a rotation, so equation really represents a transform belonging to a larger set of transforms called involutions.

An involution(对合矩阵) is a matrix that, when multiplied by itself, produces the identity matrix.

Generalizing to n dimensions, the matrix M_invol(a) represents a composition of n-1 reflections through a set of n - 1 orthogonal planes containing the vector a, whereas the matrix M_reflect(a) represents a single reflection through one plane perpendicular to the vector a.

Since M_invol  = -M_reflect, the determinant can be calculated as .

The determinant of M_reflect is always -1, so the determinant  of M_invol is -1 when n is even and +1 when n is odd.

Matrix3D MakeInvolution(const Vector3D & a) {
    float x = a.x * 2.0f;
    float y = a.y * 2.0f;
    float z = a.z * 2.0f;
    
    float axay = x * a.y;
    float axaz = x * a.y;
    float ayaz = y * a.z;

    return (Matrix3D(x * a.x - 1.0f, axay, axaz,
                             axay, y * a.y - 1.0f, ayaz,
                             axaz, ayaz, z * a.z - 1.0f));

　　2.4 Scales

A scale transform is used to enlarger or shrink an object to a new overall size.

If the scale changes the size equally in every direction, then it is called a uniform scale.

If the scale expands or compresses an object more in some directions than it does in other directions, then it is called a nonuniform scale.

A uniform scale simply multiplies all vectors v by a scale factor s so that the transformed vector v' is given by 

A nonuniform scale aligned to the coordinate axes is similar, but the scale factors appearing on the main diagonal of the transformation matrix are not all equal to the same value.

Matrix3D MakeScale(float sx, float sy, float sz) {
    return (Matrix3D(sx, 0.0f, 0.0f,
                             0.0f, sy, 0.0f,
                             0.0f, 0.0f, sz));

We can scale an object along a single arbitrary direction a while perserving the object's size in every direction orthogonal to a by decomposing a vector v into its components   

and and scaling only the parallel component.

The scaled vector v' is then given by  for a scale factor s.

When we combine the matrix representation of the components, we obtain as the transform M_scale(s,a) that scales by a factor of s along the direction a.

Matrix3D MakeScale(float s, const Vector3D & a) {
    s -= 1.0f;
    float x = a.x * s;
    float y = a.y * s;
    float z = a.z * s;

    float axay = x * a.y;
    float axaz = x * a.z;
    float ayaz = y * a.z;
    
    return (Matrix3D(x * a.x + 1.0f, axay, axaz,
                             axay, y * a.y + 1.0f, ayaz,
                             axaz, ayaz, z * a.z + 1.0f));

　　2.5 Skews

A skew transform is used to shear an object along one direction by an angle made with a perpendicular direction.

A vector v is skewed by adding a vector parallel to a of length (b · a)tanθ to it, where the factor b · v represents the length of projection of v onto b.

The skewed vector v' is given by  where the factor tanθ is the ratio of the skew distance along a to the projected length of v along b, which can be interpreted as the amount by which b · v needs to be multiplied to obtaina skew by the angle θ.

Matrix3D MakeSkew(float t, const Vector3D & a, const Vector3D & b) {
    t = tan(t);
    float x = a.x * t;
    float y = a.y * t;
    float z = a.z * t;
    
    return (Matrix3D(x * b.x + 1.0f, x * b.y, x * b.z,
                             y * b.x, y * b.y + 1.0f, y * b.z,
                             z * b.x, z * b.y, z * b.z + 1.0f));
}

　　2.6 Homogeneous Coordinates

Game engines and other types of computer graphics applications integrate location into their transforms by using a four-dimensional projective space called homogeneous coordinates.

In homogeneous coordinates, we append a fourth number called the w cooridnate to every vector so that an arbitrary vector is written as (v_x,v_y,v_z,v_w).

A point in 3D space is associated with each 4D vector v by considering a line of infinite extent that passes through the origin in 4D space and is parallel to v.

The 3D point corresponding to v is given by the x, y and z coordinates at the unique location where a point on the associated line has a w coordinate equal to one.

Because all scalar multiplies of v correspond to offsets from the origin to points on the line parallel to v, we can simply divide all of the components of v by the coordinate v_w to find the location where the line intersects the subspace for which w = 1.

Homogeneous coordinates are so named because any nonzero scalar multiple of a 4D vector v produces the same 3D point after dividing by the w coordinate.

This is a projection of an intrinsically one-dimensional object, a line, to an intrinsically zero-dimensional object, a point, accomplished by viewing only one 3D slice of 4D space.

A 3D vector v is converted to a 4D vector by appending a w coordinate equal to zero, and a 3D point p is converted to a 4D homogeneous vector by appending a w coordinate equal to one, v = (v_x, v_y, v_z, 0), p = (p_x, p_y, p_z, 1)

One of the main advantages to using homogeneous coordinates is the ability to incorporate translations into our transforms by using 4 x 4 matrices.

Recall that a general affine transformation from coordinate system A to coordinate system B is given by ,  where M is a 3 x 3 transformation matrix and t is a 3D translation vector.

These can be combined into a single 4 x 4 transformation matrix H having the form 

struct Point3D : Vector3D {
    Point3D() = default;
    Point3D(float a, float b, float c) : Vector3D(a, b, c) {}
};

inline Point3D operator+(const Point3D & a, const Vector3D & b) {
    return (Point3D(a.x + b.x, a.y + b.y, a.z + b.z));
}

inline Vector3D operator-(const Point3D & a, const Point3D & b) {
    return (Vector3D(a.x - b.x, a.y - b.y, a.z - b.z));
}

struct Transform4D : Matrix4D
{
    Transform4D() = default;

    Transform4D(float n00, float n01, float n02, float n03,
                float n10, float n11, float n12, float n13,
                float n20, float n21, float n22, float n23)
    {
        n[0][0] = n00; n[0][1] = n10; n[0][2] = n20;
        n[1][0] = n01; n[1][1] = n11; n[1][2] = n21;
        n[2][0] = n02; n[2][1] = n12; n[2][2] = n22;
        n[3][0] = n03; n[3][1] = n13; n[3][2] = n23;

        n[0][3] = n[1][3] = n[2][3] = 0.0F;
        n[3][3] = 1.0F;
    }

    Transform4D(const Vector3D& a, const Vector3D& b,
                const Vector3D& c, const Point3D& p)
    {
        n[0][0] = a.x; n[0][1] = a.y; n[0][2] = a.z;
        n[1][0] = b.x; n[1][1] = b.y; n[1][2] = b.z;
        n[2][0] = c.x; n[2][1] = c.y; n[2][2] = c.z;
        n[3][0] = p.x; n[3][1] = p.y; n[3][2] = p.z;

        n[0][3] = n[1][3] = n[2][3] = 0.0F;
        n[3][3] = 1.0F;
    }

    Vector3D& operator [](int j)
    {
        return (*reinterpret_cast<Vector3D *>(n[j]));
    }

    const Vector3D& operator [](int j) const
    {
        return (*reinterpret_cast<const Vector3D *>(n[j]));
    }

    const Point3D& GetTranslation(void) const
    {
        return (*reinterpret_cast<const Point3D *>(n[3]));
    }

    void SetTranslation(const Point3D& p)
    {
        n[3][0] = p.x;
        n[3][1] = p.y;
        n[3][2] = p.z;
    }
};

Transform4D Inverse(const Transform4D& H)
{
    const Vector3D& a = H[0];
    const Vector3D& b = H[1];
    const Vector3D& c = H[2];
    const Vector3D& d = H[3];

    Vector3D s = Cross(a, b);
    Vector3D t = Cross(c, d);

    float invDet = 1.0F / Dot(s, c);

    s *= invDet;
    t *= invDet;
    Vector3D v = c * invDet;

    Vector3D r0 = Cross(b, v);
    Vector3D r1 = Cross(v, a);

    return (Transform4D(r0.x, r0.y, r0.z, -Dot(b, t),
                        r1.x, r1.y, r1.z,  Dot(a, t),
                         s.x,  s.y,  s.z, -Dot(d, s)));
}

Transform4D operator *(const Transform4D& A, const Transform4D& B)
{
    return (Transform4D(
        A(0,0) * B(0,0) + A(0,1) * B(1,0) + A(0,2) * B(2,0),
        A(0,0) * B(0,1) + A(0,1) * B(1,1) + A(0,2) * B(2,1),
        A(0,0) * B(0,2) + A(0,1) * B(1,2) + A(0,2) * B(2,2),
        A(0,0) * B(0,3) + A(0,1) * B(1,3) + A(0,2) * B(2,3) + A(0,3),
        A(1,0) * B(0,0) + A(1,1) * B(1,0) + A(1,2) * B(2,0),
        A(1,0) * B(0,1) + A(1,1) * B(1,1) + A(1,2) * B(2,1),
        A(1,0) * B(0,2) + A(1,1) * B(1,2) + A(1,2) * B(2,2),
        A(1,0) * B(0,3) + A(1,1) * B(1,3) + A(1,2) * B(2,3) + A(1,3),
        A(2,0) * B(0,0) + A(2,1) * B(1,0) + A(2,2) * B(2,0),
        A(2,0) * B(0,1) + A(2,1) * B(1,1) + A(2,2) * B(2,1),
        A(2,0) * B(0,2) + A(2,1) * B(1,2) + A(2,2) * B(2,2),
        A(2,0) * B(0,3) + A(2,1) * B(1,3) + A(2,2) * B(2,3) + A(2,3)));
}

Vector3D operator *(const Transform4D& H, const Vector3D& v)
{
    return (Vector3D(H(0,0) * v.x + H(0,1) * v.y + H(0,2) * v.z,
                     H(1,0) * v.x + H(1,1) * v.y + H(1,2) * v.z,
                     H(2,0) * v.x + H(2,1) * v.y + H(2,2) * v.z));
}

Point3D operator *(const Transform4D& H, const Point3D& p)
{
    return (Point3D(H(0,0) * p.x + H(0,1) * p.y + H(0,2) * p.z + H(0,3),
                    H(1,0) * p.x + H(1,1) * p.y + H(1,2) * p.z + H(1,3),
                    H(2,0) * p.x + H(2,1) * p.y + H(2,2) * p.z + H(2,3)));
}

　　2.7 Quaternions

　　　　2.7.1 Quaternion Fundamentals

A typical quaternion q has four components that can be written as  where x,y,z and w are real number.

It does'nt matter what order these components are written in because multiplication by i,j and k provide all the necessary identification for the imaginary terms.

It is often convenient to write a quaternion in the form q = v + s, where v, called the vector part, corresponds to the imaginary triplet(x, y, z), and s called the scalar part,corresponds to the real component w.

As with ordinary vectors and complex numbers, quaternion addition is performed componentwise.

Multiplication, however, follows the rule given by Hamilton, which can also be expressed in the more explicit form

If we represent the quaternions by q₁ = v₁ + s₁ and q₂ = v₂ + s₂ instead, then the product can be written as 

The only noncommutative piece appearing in Equation is the cross product, a fact from which we can quickly deduce that reversing the order of the factors in quaternion multiplication changes the product by twice the cross product between the vector parts, as stated by 

This expose the fact that two quaternions commute only if their vector parts are parallel because when that is the case, the cross product v₁ x v₂ is zero.

A quaternion q has a conjugate(共轭) denoted by q^* that is similar to the complex conjugate except that we are now negating three imaginary components instead of just one.

That is, the conjugate of quaternion q = v + s is given by 

The product of a quaternion and its conjugate gives us, which is a real number that we equate to the squared magnitude of the quaternion.

We denote the magnitude of a quaternion using two vertical bars, as with ordinary vectors, and define it as 

As with vectors, multiplying a quaternion q by a scalar value t has the effect of multiplying the magnitude of q by |t|.

Quaternions also have the property that the magnitude of the product of two quaternions q₁ and q₂ is equal to the product of their individual magnitudes,

struct Quaternion {
    float x, y, z, w;
    Quaternion() = default;
    Quaternion(float a, float b, float c, float s) {
        x = a;
        y = b;
        z = c;
        w = s;
    }

    Quaternion(const Vector3D& v, float s) {
        x = v.x;
        y = v.y;
        z = v.z;
        w = s;
    }
    
    const Vector3D& GetVectorPart() const {
        return (reinterpret_cast<const Vector3D&>(x));
    }

    Matrix3D GetRotationMatrix();
    void SetRotationMatrix(const Matrix3D & m);
};

Quaternion operator*(const Quaternion & q1, const Quaternion & q2) {
    return (Quaternion(
        q1.w * q2.x + q1.x * q2.w + q1.y * q2.z - q1.z * q2.y,
        q1.w * q2.y - q1.x * q2.z + q1.y * q2.w + q1.z * q2.x,
        q1.w * q2.z + q1.x * q2.y - q1.y * q2.x + q1.z * q2.w,
        q1.w * q2.w - q1.x * q2.x - q1.y * q2.y - q1.z * q2.z));
}

　　　　2.7.2 Rotations With Quaternions

Given a quaternion q = xi + yj + zk + w and a vector v = (v_x, v_y, v_z), a rotation is performed by considering the vector to be the quaternion v_xi + v_yj + v_zk and calculating a new vector v' with the product 

To be clear, each of the products in this equation is a quaternion multiplied by a quaternion.

This is sometimes called the sandwich product because the quaternion v is sandwiched between the quaternion q and its inverse

The quaternion v is known as a pure quaternion, which is any quaternion that has a zero scalar component and is thus made up of only imaginary terms.

When v is a pure quaternion, the sandwich product qvq^-1 always yields another pure quaternion.

Since we have established an equivalence between vectors and pure quaternion, we can say that the sandwich product transform a vector v into another vector v'.

The magnitude of q in Equation doesn't matter, as long as it's nonzero, because if ||q|| = m, then m can be factored out of q, and 1/m can be factored out of q^-1.

These cancel each other out and leave quaternions with magnitudes of one behind.

A quaternion q having a magnitude of one is called a unit quaternion, and it has the special property that its inverse is simply equal to its conjugate because qq^* = 1.

In the case that q is a unit quaternion, Equation simpifies to 

The set of unit quaternions form a multiplicatively closed subset of H because the product of any two unit quaternions is another unit quaternion

For this reason and the fact that vector transforms become simpler, only unit quaternions are typically used to represent rotations in practice.

The product qv is given by 

When we multiply this by q^* = -b + c, we get  after some simplification that includes an application of the vector triple product identity

-b x v x b = (b · v)b - b²v. 

If we set b = sa, where s = ||b|| and a is a unit vector, then we can write qvq^* as 

The right side of this equation has the same three terms that appear in the formula for rotation about an arbitrary axis a given by equation

except that the scalar coefficients are written in a different way.

In order for Equation  to perform a rotation through an angle θ, the values of c and s must satisfy the equalities

All three of these requirements are satisfied when we choose  because these values produce valid trigonometric identities.

We conclude that the quaternion represents a rotation through the angle θ about the unit-length axis a that can be applied to a vector v using the sandwich product qvq^*.

As with all the rotations previously described in this chapter, a quaternion rotation through a positive angle is counterclockwise when the axis points toward the viewer.

// This code rotate the vector v using the quaternion q by calculating  the sandwich product.It is assumed that q is a unit quaternion
Vector3D Transform(const Vector3D & v, const Quaternion & q) {
    const Vector3D & b = q.GetVectorPart();
    float b2 = b.x * b.x + b.y * b.y + b.z * b.z;
    return (v * (q.w * q.w - b2) + b * (Dot(v, b) * 2.0f) + Cross(b, v) * (q.w * 2.0f));
}

One advantage to using quaternions is that multiple rotations can easily be composed.

To first rotate a vector v using a quaternion q1 and then rotate the result using another  quaternion q2, we calculate the sandwich product of a sandwich product as in 

By reassociating the factors, this can be written as  showing that the two successive rotations are equivalent to a single rotation using the quaternion given by q2q1.

The product of two quaternions can be calculated with 16 multiplies and 12 adds, and that has a significantly lower cost than the 27 multiplies and 18 adds required to calculate the product of two 3 x 3 matrices.

A quaternion also has the advantage that it has much lower storage requirements because it comprises only four floating-point components compared to the nine floating-point entries needed by an equivalent 3 x 3 rotation matrix.

// This function creates a 3 x 3 matrix that corresponds to the Quaternion data structure for which it's called. It is assumed that the  quaternion has a magnitude of one.
Matrix3D Quaternion::GetRotationMatrix(void)
{
    float x2 = x * x;
    float y2 = y * y;
    float z2 = z * z;
    float xy = x * y;
    float xz = x * z;
    float yz = y * z;
    float wx = w * x;
    float wy = w * y;
    float wz = w * z;

    return (Matrix3D(
        1.0F - 2.0F * (y2 + z2), 2.0F * (xy - wz), 2.0F * (xz + wy),
        2.0F * (xy + wz), 1.0F - 2.0F * (x2 + z2), 2.0F * (yz - wx),
        2.0F * (xz - wy), 2.0F * (yz + wx), 1.0F - 2.0F * (x2 + y2)));
}

// This function sets the members of the Quaternion data structure for which it's called to the values corresponding to a quaternion equivalent to the 3 x 3 rotation matrix m.void Quaternion::SetRotationMatrix(const Matrix3D& m)
{
    float m00 = m(0,0);
    float m11 = m(1,1);
    float m22 = m(2,2);
    float sum = m00 + m11 + m22;

    if (sum > 0.0F)
    {
        w = sqrt(sum + 1.0F) * 0.5F;
        float f = 0.25F / w;

        x = (m(2,1) - m(1,2)) * f;
        y = (m(0,2) - m(2,0)) * f;
        z = (m(1,0) - m(0,1)) * f;
    }
    else if ((m00 > m11) && (m00 > m22))
    {
        x = sqrt(m00 - m11 - m22 + 1.0F) * 0.5F;
        float f = 0.25F / x;

        y = (m(1,0) + m(0,1)) * f;
        z = (m(0,2) + m(2,0)) * f;
        w = (m(2,1) - m(1,2)) * f;
    }
    else if (m11 > m22)
    {
        y = sqrt(m11 - m00 - m22 + 1.0F) * 0.5F;
        float f = 0.25F / y;

        x = (m(1,0) + m(0,1)) * f;
        z = (m(2,1) + m(1,2)) * f;
        w = (m(0,2) - m(2,0)) * f;
    }
    else
    {
        z = sqrt(m22 - m00 - m11 + 1.0F) * 0.5F;
        float f = 0.25F / z;

        x = (m(0,2) + m(2,0)) * f;
        y = (m(2,1) + m(1,2)) * f;
        w = (m(1,0) - m(0,1)) * f;
    }
}

　　Exercises for Chapter2

Chapter3 Geometry

Most of the computation performed bny a game engine involves some kind of geometry.

Geometry defines the world in which a game takes palce, geometry describes the characters in a game and their movements, geometry tells the graphics haradware how to render a scene, geometry allows an engine to determine what's visible to the camera, and geometry is necessary for detecting collisions between various objects.　　

　　3.1 Triangle Meshes

A triangle mesh is a collection of triangles that fit together to model the surface of a solid volume.

At a minimum, the data associated with a triangle mesh includes a list of vertex positions stored as 3D points with floating-point coordinates.

In most cases, the data also includes an index list that contains a triplet of integer indices for each triangle in the mesh specifying which three vertices define the triangle's boundary.

 

There are always multiple ways to break a polygon having at least four vertices into triangles.

The particular choice of triangles composing an entire mesh is called its triangulation.

A triangle mesh is called closed if it is the case that every edge is used by excatly two triangles.

That is, for any pair of vertices used by one triangle, there must be one more triangle, and no others, that also uses the same pair of vertices for one of its edges but does not share the third vertex with the first triangle.

A closed triangle mesh satisfies the Euler formula, which states  where V is the number of vertices, E is the number of edges, and F is the number of faces.

An important property of a triangle mesh is the winding direction used by its triangles.

In Figure 3.1, the lower-left triangle is wound in the counterclockwise direction when its vertices are referenced in the order (0, 1, 4), (1, 4, 0), or (4, 0, 1);

　　3.2 Normal Vectors

A normal vector, or just normal for short, is a vector that is perpendicular to a surface, and the direction in which it points is said to be normal to the surface.

A flat plane has only one normal direction, but most surfaces aren't so simple and thus have normal vectors that vary from point to point.

　　　　3.2.1 Calculating Normal Vectors

In the case that a surface is defined implicity by a scalar function f(p), the normal vector at p = (x, y, z) is given by the gradient  because it is perpendicular to every direction tangent to the level surface(等高面) of f at p.

Suppose we have the ellipsoid defined by the equation 

The point lies on the surface of this ellipsoid, and the normal vector n at that point is given by 

Calculating normal vector with the gradient is something that's usually done only in the process of constructing a triangle mesh to approximate an ideal surface described by some mathematical formula.

The normal vectors are typically scaled to unit length and stored with the vertex coordinates that they're associated with.

Most of the time, a game engine is working with a triangle mesh having an arbitrary shape that was created in a modeling program, and the only information available is the set of vertex coordinates and the list of of indices that tell how vertices are grouped into triangles.

An outward-facing normal vector is then given by 

Any permutation of the subscript that keeps them in the same cyclic order produces the same normal vector.

It doesn't matter which vertex is chosen to be subtracted from the other two as long as the first factor in the cross product involves the next vertex in the counterclockwise order.

If the order is reversed, then the calculated normal vector still lies along the same line, but it points in the opposite direction.

To calculate per-vertex normal vectors for a triangle mesh, it is typical for a game engine's model processing pipeline to calculate all of the per-face normal vectors and then take an average at each vertex over all of the faces that use that vertex.

The average may be weighted based on triangle area or other factors to create a smooth field of normal vectors over a curved surface.

In case in which a hard edge is desired, such as for a cube or the pyramid, vertex positions are typically duplicated, and different normal vectors corresponding to different faces are associated with the various copies of the vertex coordinates.

　　　　3.2.2 Transforming Normal Vectors

 

　　3.3 Lines and Rays

　　　　3.3.1 Parametric Lines

　　　　3.3.2 Distance Between a Point and a Line

　　　　3.3.3 Distance Between Two Lines

　　3.4 Planes

　　　　3.4.1 Implicit Planes

　　　　3.4.2 Distance Between a Point and a Plane

　　　　3.4.3 Reflection Through a Plane

　　　　3.4.4 Intersection of a Line and a Plane

　　　　3.4.5 Intersection of Three Planes

　　　　3.4.6 Intersection of Two Planes

　　　　3.4.7 Transforming Planes

　　3.5 Plucker Coordinates

　　　　3.5.1 Implicit Lines

　　　　3.5.2 Homogeneous Formulas

　　　　3.5.3 Transforming Lines

　　Exercises for Chapter3

Chapter4 Adavanced Algebra

　　4.1 Grassmann Algebra

　　　　4.1.1 Wedge Product

　　　　4.1.2 Bivectors

　　　　4.1.3 Trivectors

　　　　4.1.4 Algebraic Structure

　　　　4.1.5 Complements

　　　　4.1.6 Antivectors

　　　　4.1.7 Antiwedge Product

　　4.2 Projective Geometry

　　　　4.2.1 Lines

　　　　4.2.2 Planes

　　　　4.2.3 Join and Meet

　　　　4.2.4 Line Crossing

　　　　4.2.5 Plane Distance

　　　　4.2.6 Summary and Implementation

　　4.3 Matrix Inverses

　　4.4 Geometric Algebra

　　　　4.4.1 Geometric Product

　　　　4.4.2 Vector Division

　　　　4.4.3 Rotors

　　4.5 Conclusion

　　Exercises for Chapter4