Saturday, March 31, 2012

3D Rotation Algorithm about arbitrary axis with C/C++ code

When an object is to be rotated about an axis that is not parallel to one of the coordinate axes, we need to perform some additional transformations. In this case, we also need rotations to align the axis with a selected coordinate axis and to bring the axis back to its original orientation. Given the specifications for the rotation axis and the rotation angle, we can accomplish the required rotation in five steps.
  1. Translate the object so that the rotation axis passes through the coordinate origin.
  2. Rotate the object so that the axis of rotation coincides with one of the coordinate axes.
  3. Perform the specified rotation about that coordinate axis.
  4. Apply inverse rotations to bring the rotation axis back to its original orientation.
  5. Apply the inverse translation to bring the rotation axis back to its original position.
Let (u, v, w) be a vector that specify axis about which the object is to be rotated. Suppose the vector passes through a point (a, b, c). Then the general composite matrix for the rotation can be written as:
rotationMatrix 
Here L is magnitude of the axis vector. C source code for the arbitrary rotation is given below. Here (a, b, c) is assumed to be origin for simplicity.
Source Code
#include <iostream>
#include <cmath> 
using namespace std;
typedef struct {     float x;     float y;     float z; }Point; Point points;
float rotationMatrix[4][4]; float inputMatrix[4][1] = {0.0, 0.0, 0.0, 0.0}; float outputMatrix[4][1] = {0.0, 0.0, 0.0, 0.0};
void showPoint(){     cout<<"("<<outputMatrix[0][0]<<","<<outputMatrix[1][0]<<","<<outputMatrix[2][0]<<")"<<endl; }
void multiplyMatrix() {     for(int i = 0; i < 4; i++ ){         for(int j = 0; j < 1; j++){             outputMatrix[i][j] = 0;             for(int k = 0; k < 4; k++){                 outputMatrix[i][j] += rotationMatrix[i][k] * inputMatrix[k][j];             }         }     } } void setUpRotationMatrix(float angle, float u, float v, float w) {     float L = (u*u + v * v + w * w);     angle = angle * M_PI / 180.0; //converting to radian value     float u2 = u * u;     float v2 = v * v;     float w2 = w * w;
    rotationMatrix[0][0] = (u2 + (v2 + w2) * cos(angle)) / L;     rotationMatrix[0][1] = (u * v * (1 - cos(angle)) - w * sqrt(L) * sin(angle)) / L;     rotationMatrix[0][2] = (u * w * (1 - cos(angle)) + v * sqrt(L) * sin(angle)) / L;     rotationMatrix[0][3] = 0.0;
    rotationMatrix[1][0] = (u * v * (1 - cos(angle)) + w * sqrt(L) * sin(angle)) / L;     rotationMatrix[1][1] = (v2 + (u2 + w2) * cos(angle)) / L;     rotationMatrix[1][2] = (v * w * (1 - cos(angle)) - u * sqrt(L) * sin(angle)) / L;     rotationMatrix[1][3] = 0.0;
    rotationMatrix[2][0] = (u * w * (1 - cos(angle)) - v * sqrt(L) * sin(angle)) / L;     rotationMatrix[2][1] = (v * w * (1 - cos(angle)) + u * sqrt(L) * sin(angle)) / L;     rotationMatrix[2][2] = (w2 + (u2 + v2) * cos(angle)) / L;     rotationMatrix[2][3] = 0.0;
    rotationMatrix[3][0] = 0.0;     rotationMatrix[3][1] = 0.0;     rotationMatrix[3][2] = 0.0;     rotationMatrix[3][3] = 1.0; }
int main() {     float angle;     float u, v, w;     cout<<"Enter the initial point you want to transform:";     cin>>points.x>>points.y>>points.z;     inputMatrix[0][0] = points.x;     inputMatrix[1][0] = points.y;     inputMatrix[2][0] = points.z;     inputMatrix[3][0] = 1.0;
    cout<<"Enter axis vector: ";     cin>>u>>v>>w;
    cout<<"Enter the rotating angle in degree: ";     cin>>angle;
    setUpRotationMatrix(angle, u, v, w);     multiplyMatrix();     showPoint();
    return 0; }

11 comments:

  1. Thanx a lot dude. Was having one hell of a time programming with all the matrices. Cheers!

    ReplyDelete
  2. Thanx a lot!

    You have helped save my master's thesis from vanishing into oblivion!

    Regards,
    -Bala.
    IIT-Bombay.

    ReplyDelete
  3. "Here L is magnitude of the axis vector"
    Magnitude means length of vector.

    But in the source : float L = (u*u + v * v + w * w);
    And that is the square of the magnitude

    Example.
    Vector (0, 0, 2).
    It is trivial to compute it's length, it is 2.
    But
    L = 0*0 + 0*0 + 2*2 = 4

    The source code is correct, and where magnitude is required sqrt(L) is used. Only the explanation needs amends to properly describe L as the square of the magnitude.

    ReplyDelete
    Replies
    1. Thank you for your quick catch up and reading the post in detail. It is a simple logical mistake. I appreciate that I will make the correction

      Delete
  4. angle = angle * M_PI / 180.0; //converting to radian value

    ^^ that line saved me so much time. I was getting so frustrated because i knew my matrix multiplication was right, but my rotations were completely broken. thanks man.

    ReplyDelete
  5. Hello,

    Can you please give a reference to this, how did you obtain the matrix? Thank you!

    ReplyDelete
  6. Thanks, this was eye opening

    ReplyDelete
  7. Thanks, that was really helpful! But I wanted to ask you a question.
    I used your code to rotate an object (an ellipsoid for example) but almost have of the voxels went missing, what should I do in this case?
    Thanks

    ReplyDelete
  8. Hi! Thanks first of all for your awesome work here.
    I have a question about the rotation matrix, specifically about the 4th column.
    Shouldn't the first three entries of the 4th column be permutation of each other (in both a,b,c and u,v,w)? It would seem that there is an equal term in the 1st and 3rd entry: u*(b*v + c*w).
    I was wondering if maybe that term should be w*(a*u + b*v) in the 3rd entry of the 4th column (and thus the second permutation of the term in the 1st entry and the first permutation of the term in the 2nd entry.)

    Best regards,

    Peter

    ReplyDelete