#include #include "GraphicsGems.h" /*=========================================================================* * R A N D _ R O T A T I O N Author: Jim Arvo, 1991 * * * * This routine maps three values (x[0], x[1], x[2]) in the range [0,1] * * into a 3x3 rotation matrix, M. Uniformly distributed random variables * * x0, x1, and x2 create uniformly distributed random rotation matrices. * * To create small uniformly distributed "perturbations", supply * * samples in the following ranges * * * * x[0] in [ 0, d ] * * x[1] in [ 0, 1 ] * * x[2] in [ 0, d ] * * * * where 0 < d < 1 controls the size of the perturbation. Any of the * * random variables may be stratified (or "jittered") for a slightly more * * even distribution. * * * *=========================================================================*/ void rand_rotation( float x[], Matrix3 *M ) { float theta = x[0] * PITIMES2; /* Rotation about the pole (Z). */ float phi = x[1] * PITIMES2; /* For direction of pole deflection. */ float z = x[2] * 2.0; /* For magnitude of pole deflection. */ /* Compute a vector V used for distributing points over the sphere */ /* via the reflection I - V Transpose(V). This formulation of V */ /* will guarantee that if x[1] and x[2] are uniformly distributed, */ /* the reflected points will be uniform on the sphere. Note that V */ /* has length sqrt(2) to eliminate the 2 in the Householder matrix. */ float r = sqrt( z ); float Vx = sin( phi ) * r; float Vy = cos( phi ) * r; float Vz = sqrt( 2.0 - z ); /* Compute the row vector S = Transpose(V) * R, where R is a simple */ /* rotation by theta about the z-axis. No need to compute Sz since */ /* it's just Vz. */ float st = sin( theta ); float ct = cos( theta ); float Sx = Vx * ct - Vy * st; float Sy = Vx * st + Vy * ct; /* Construct the rotation matrix ( V Transpose(V) - I ) R, which */ /* is equivalent to V S - R. */ M->element[0][0] = Vx * Sx - ct; M->element[0][1] = Vx * Sy - st; M->element[0][2] = Vx * Vz; M->element[1][0] = Vy * Sx + st; M->element[1][1] = Vy * Sy - ct; M->element[1][2] = Vy * Vz; M->element[2][0] = Vz * Sx; M->element[2][1] = Vz * Sy; M->element[2][2] = 1.0 - z; /* This equals Vz * Vz - 1.0 */ }