Randomly Generate Orthogonal 3x3 Matrix

1.6k views Asked by At

I'm looking to do some complex part analysis within Seimens NX. I'm looking to implement the double caliper method of measuring a model in order to find the minimum possible box that it could possibly fit into(for machining purposes). I've got all of my measurement code in place, but I am completely baffled by the idea of a construct that can randomly output normalized 3x3 vectors for use as coordinate systems. The part is measured with respect to this coordinate system, so each coordinate system gives a unique "minimum part envelope". Once analyzed, the smallest envelope is selected and displayed.

this is the type of vector I am talking about:
1 0 0
0 1 0
0 0 1

numbers can be any value between -1 and 1, with decimals not only being accepted but pretty much required.

and no, this isn't my homework. More of an individual pursuit in my free time at work.

2

There are 2 answers

1
nsanders On BEST ANSWER

If you apply a rotation matrix to an already orthogonal matrix, then the result should also be orthogonal.

So you can redefine your problem as applying a random rotation matrix to the identity matrix.

Perhaps do one random rotation matrix for each axis (x,y,z) and then apply the matrices themselves in a random order?

0
Hanzhe Teng On

If you don't mind to consider only a special subset of the orthogonal matrices, there is an easier way to achieve this, which is to take advantage of the Rodrigues' rotation formula to generate rotation matrices (which has an additional constraint that its determinant is equal to 1).

Image of the Rodrigues' rotation formula

Image of the bracket notation

With this, you only need to generate a random 3x1 unit vector (as the rotation axis) and specify a rotation angle. This formula will transform them into a valid rotation matrix.

MATLAB example:

function R = rot(w, theta)
  bw = [0, -w(3), w(2); w(3), 0, -w(1); -w(2), w(1), 0];
  R = eye(3) + sin(theta)*bw + (1-cos(theta))*bw*bw;
end

w = rand(3,1)
w = w/norm(w)
R = rot(w, 3.14)

C++ example:

// w: the unit vector indicating the rotation axis
// theta: the rotation angle in radian
Eigen::Matrix3d MatrixExp3 (Eigen::Vector3d w, float theta){
  Eigen::Matrix3d bw, R;
  bw << 0, -w(2), w(1), w(2), 0, -w(0), -w(1), w(0), 0;
  R << Eigen::Matrix3d::Identity() + std::sin(theta)*bw + (1-std::cos(theta))*bw*bw;
  return R;
}

int main() {
  std::srand((unsigned int) time(0));
  Eigen::Vector3d w = Eigen::Vector3d::Random();
  Eigen::Matrix3d R = MatrixExp3(w.normalized(), 3.14f);
  std::cout << R << std::endl;
}