ANIMATING TRANSFORMATION
pbrt supports keyframe matrix animation for cameras and geometric primitive in the scene. Tanther than just supplying a single transformation to place the corresponding object in the scene, the user may supply a number of keyframe transformations, each one associated with a particular point in time. This makes it possible for the camera to move and for objects in scene to be moving during the time the simulated camera's shutter is open. In genenral, the problem of interpolating between kwyframe matrices is under-defined.
Keyframe matrix interpolation is an important problem in computer animation, where a number of different approaches have been developed. Fortunately, the problem of matrix interpolation in a renderer is generally less chanlenging than it is for animation systems for two reasons.
SRE实战 互联网时代守护先锋，助力企业售后服务体系运筹帷幄！一键直达领取阿里云限量特价优惠。First, in a renderer like pbrt, we generally have a keyframe matrix at the camera shutter open time and another at the shutter close time; we only need to interpolate between the two of them across the time of a single image. In animation systems, the matrices are generally available at a lower time frequency, so that there are many frames between pairs of keyframe matrices; as such, there is more opportunity to notice shortcomings in the interpolation.
Second, in a physically based renderer, the longer the period of time over which we need to interpolate the pair of matrices the longer the virtual camera shutter is open and the more motion blur there will be in the final image; the increased amount of motion blur often hides sins of the interpolation.
The most straightforward approach to interpolate transformations defined bykeyframe matrices - directly interpolating the individual components of the matrices - is not a good one, as it will generally lead to unexpected and undesirable results. For example, if the transformations apply different rotations, then even if we have a rigid-body motion, then even if we have a rigid-body motion, the intermediate matrices may scale the object, which is clearly undersirable.
The approach used for transformation interpolation in pbrt is based on matrix decomposition - given an arbitrary transformation matrix M, we decompose it into a concatentation of scale (S), rotation (R), and translation (T) transformations, $ M = SRT $, where each of those components is independently interpolated and then the composite interpolated marix is found by multiplying the three interpolated matrices together.
Interpolation of translation and scale can be performed easily and accurately with linear interpolation of components of their matrices; interpolating rotations is more difficult. Before describing the matrix decomposition implementation in pbrt, we will first introduce the quaternions, an elegant represnentation of rotations that leads to effective methods for interpolating them.
QUATERNIONS: Quaternions were originally invented by Sir William Rowan Hamilton in 1843 as a generalization of complex numbers. He determined that just as in two dimensions (x, y), where complex numbers could be defined as a sum of a real and an imaginary part $x+iy$, with $i^2 = -1$, a generaliztion could be made to four dimensions, giving quaternions. A quaternion is a four-tuple, $q = (x, y, z, w) = w + xi +yj + zk$ where i, j, and k are defined so that $i^2 = j^2 = k^2 = ijk = -1$. Other important relationships between the components are that $ij = k$ and $ji = -k$. This implies that quaternion multiplication is generally not communitative.
A quaternion can be represented as a quadruple $\mathbf{q} = (q_x, q_y, q_z, q_w)$ or as $\mathbf{q} = (\mathbf{q}_{xyz}, q_w)$, where $\mathbf{q}_{xyz}$ is an imaginary 3-vector and $q_w$ is the real part. We will use both representations. An expression for the product of two arbitrary quaternions can be found by expanding their definition in terms of real and imaginary components: $\mathbf{qq^{'}} = (q_w + q_xi + q_yj + q_zk)(q_w^{'} +q_x^{'}i + q_y^{'}j + q_z^{'}k)$. Collecting terms and using identities among the components like those listed above (e.g., $i^2 = -1$). the result can be expressed concisely using vector cross and dot products: $(\mathbf{qq^{'}})_{xyz} = \mathbf{q}_{xyz} \times \mathbf{q}^{'}_{xyz} + q_w \mathbf{q}^{'}_{xyz} + q_w^{'} \mathbf{q}_{xyz}$. There is a useful relationship between unit quaternion ($x^2 + y^2 +z^2 +w^2 = 1$) and the space of rotations in $R^3$: specifically, a rotation of angle $2\theta$ about a unit axis $\mathbf{v}$ can be mapped to a unit quaternion $(\mathbf{v}\cos\theta, \cos\theta)$, in which case the following quaternion product is equivalent to applying the rotation to a point p expressed in homogeneous coordinate form: $p^{'} = qpq^{-1}$.
It is useful to be able to compute the transformation matrix that represents the same rotation as a quaternion. In particular, after interpolating rotations with quaternions in the AnimatedTransform class, we will need to convert the interpolated rotation back to a transformation matrix to compute the final composite interpolated transformation.
QUATERNION INTERPOLATION: The last quaternion function we will define, Slerp(), interpolates between two quaternions using spherical linear interpolation. Spherical linear interpolation gives constant speed motion along great circle arcs on the surface of a sphere and consequently has two desirable properties for interpolation rotations: (1) The interpolated rotation path exhibits torque minimization: the path to get between two rotations is the shortest possible path in rotation space. (2) The interpolation has constant angular velocity: the relationship between change in the animation parameter t and the change in the resulting rotation is constant over the course of interpolation.
Spherical linear interpolation for quaternions was originally presented as follows, where two quaternions $\mathbf{q}_1$ and $\mathbf{q}_2$ are given and t belongs to $[0, 1]$: $slerp(\mathbf{q}_1, \mathbf{q}_2, t) = \frac{q_1\sin(1-t)\theta + q_2 \sin t \theta}{\sin\theta}$. Now, we pewsent an intuitive way to understand Slerp(). As context, given the quaternions to interpolate between $\mathbf{q}_1 and \mathbf{q}_2$, denote by $\theta$ the angle between them. Then, given a parameter value $t \in [0, 1]$, we'd like to find the intermediate quaternion $\mathbf{q}^{'}$ that makes angle ${\theta}^{'} = \theta t$ between it and $\mathbf{q}_1$, along the path from $\mathbf{q}_1$ to $\mathbf{q}_2$.
An easy way to compute $\mathbf{q}^{'}$ is to first compute an orthogonal coordinate system in he space of quaternions where one axis is $\mathbf{q}_1$ and the other is a quaternion orthogonal to $\mathbf{q}_1$ such that the two axes form a basis that spans $\mathbf{q}_1$ and $\mathbf{q}_2$. Given such a coordinate system, we can compute rotations with respect to $\mathbf{q}_1$. An orthogonal vector $\mathbf{q}_{\perp}$ can be found by projecting $\mathbf{q}_1$ onto $\mathbf{q}_2$ and then subtracting the orthogonal rpojection from $\mathbf{q}_2$; the remainder is guanranteed to be orthogonal to $\mathbf{q}_1$: $\mathbf{q}_{\perp} = \mathbf{q}_2 - (\mathbf{q}_1 * \mathbf{q}_2) \mathbf{q}_1$. Given the coordinate system, quaternions along the animation path are given by $\mathbf{q}^{'} = \mathbf{q}_{1} \cos (\theta t) + \mathbf{q}_{\perp} \sin(\theta t)$.
The implemetation of the Slerp() function checks to see if the two quaternions are nearly parallel, in which case it uses regular linear interpolation of quaternion components in order to avoid numerical instability.
1 Quaternion Slerp(Float t, const Quaternion &q1, const Quaternion &q2) 2 { 3 Float cosTheta = Dot(q1, q2); 4 if (cosTheta > .9995f) 5 return Normalize((1-t)*q1 + t*q2); 6 else { 7 Float theta = std::acos(Clamp(cosTheta, -1, 1)); 8 Float thetap = theta * t; 9 Quaternion qperp = Normalize(q2-q1 * cosTheta); 10 return q1 * std::cos(thetap) + qperp * std::sin(thetap); 11 } 12 }
ANIMATEDTRANSFORM IMPLEMENTATION
Given the foundations of the quaternion infrastracture, we can now implement the AnimatedTransform class, which implements keyframe transformation interpolation in pbrt. Its constructor takes two transformations and time values they are associated with. As mentioned earlier, AnimatedTransform decomposes the given composite transformation matrices into a scaling, rotation, and translation components. The decomposition is performed by the AnimatedTransform:Decompose() method.
1 class AnimatedTransform { 2 public: 3 // AnimatedTransform Public Methods 4 AnimatedTransform(const Transform *startTransform, Float startTime, const Transform *endTransform, Float endTime); 5 static void Decompose(const Matrix4x4 &m, Vector3f *T, Quaternion *R, Matrix4x4 *S); 6 void Interpolate(Float time, Transform *t) const; 7 Ray operator()(const Ray &r) const; 8 RayDifferential operator()(const RayDifferential &r) const; 9 Point3f operator()(Float time, const Point3f &p) const; 10 Vector3f operator()(Float time, Vector3f &v) const; 11 bool HasScale() const { 12 return startTransform->HasScale() || endTransform->HasScale(); 13 } 14 Bounds3f MotionBounds(const Bounds3f &b) const; 15 Bounds3f BoundPointMotion(const Point3f &p) const; 16 17 private: 18 // AnimatedTransform Private Data 19 const Transform *startTransform, *endTransform; 20 const Float startTime, endTime; 21 const bool actuallyAnimated; 22 Vector3f T[2]; 23 Quaternion R[2]; 24 Matrix4x4 S[2]; 25 bool HasRotation; 26 struct DerivativeTerm { 27 DerivativeTerm() {} 28 DerivativeTerm(Float c, Float x, Float y, Float z): kc(c), kx(x), ky(y), kz(z) {} 29 Float kc, kx, ky, kz; 30 Float Eval(const Point3f &p) const { 31 return kc + kx * p.x + ky * p.y + kz * p.z; 32 } 33 }; 34 DrerivativeTerm c1[3], c2[3], c3[3], c4[3], c5[3]; 35 };
1 AnimatedTransform::AnimatedTransform(const Transform *startTransform, Float startTime, const Transform *endTransform, Float endTime): startTransform(startTransform), endTransform(endTransform), startTime(startTime), endTime(endTime), actuallyAnimated(*startTransform != *endTransform) { 2 Decompose(startTransform->m, &T[0], &R[0], &S[0]); 3 Decompose(endTransform->m, &T[1], &R[1], &S[1]); 4 <Flip R[1] if needed to select shortest path> 5 hasRotation = Dot(R[0], R[1]) < 0.9995f; 6 <Compute terms of motion derivative function> 7 }
Given the composite matrix for a transformation, information has been lost about any individual transformations that were composed to compute it. We need to choose a canonical sequence of transformations for the decomposition. We will handle only affine transformations here, which is needed for animating cameras and geometric primitives in a rendering system; perspective transformations are not generally relevant to animation of objects like these.
The transformation decomposition we will use is the following: $M = TRS$, where M is the given transformation, T is a translation, R is a rotation, and S is a scale. S is actually a generalized scale that represents a scale in some ordinate system, just not necessarily the current one. In any case, it can still be correctly interpolated with linear interpolation of its components. The Decompose() method computes the decomposition given a Matrix4x4.
1 void AnimatedTransform::Decompose(const Matrix4x4 &m, Vector3f *T, Quaternion *Rquat, Matrix4x4 *S) { 2 <Extract translation T from transformation matrix> 3 <Compute new transformation matrix M without translation> 4 <Extract rotation R from transformation matrix> 5 <Compute scale S using rotation and original matrix> 6 }
1 <Extract translation T from transformation matrix> = 2 T->x = m.m[0][3]; 3 T->y = m.m[1][3]; 4 T->z = m.m[2][3];
Since we assuming an affine transformation, after we remove the translation, what is left is the upper 3*3 matrix that represents scaling and rotation together. This matrix is copied into a new matrix M for further processing.
1 <Compute new transformation matrix M without translation> = 2 Matrix4x4 M = m; 3 for (int i = 0; i < 3; ++i) 4 M.m[i][3] = M.m[3][i] = 0.f; 5 M.m[3][3] = 1.f;
Next we'd like to extrac the pure rotation component of M. We'll use a technique called polar decomposition to do this. It can be shown that the polar decomposition of a mtrix M into rotation R and scale S can be computed by successively averaging M with its inverse transpose $M_{i+1} = 0.5 (M_i + (M_i^T)^{-1})$ until convergence, at which point $M_i = R$. In practice, this series generally converges quickly.
