3D Rotation Algorithm about arbitrary axis with C/C++ code
- Translate the object so that the rotation axis passes through the coordinate origin.
- Rotate the object so that the axis of rotation coincides with one of the coordinate axes.
- Perform the specified rotation about that coordinate axis.
- Apply inverse rotations to bring the rotation axis back to its original orientation.
Let (u, v, w) be a vector that specifies the 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:
Here L is the magnitude of the axis vector. C source code for the arbitrary rotation is given below. Here (a, b, c) is assumed to be the 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; }
Thanx a lot dude. Was having one hell of a time programming with all the matrices. Cheers!
Thanx a lot!
You have helped save my master's thesis from vanishing into oblivion!
Regards,
-Bala.
IIT-Bombay.
I am so glad about it Sir..
"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.
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
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.
Hello,
Can you please give a reference to this, how did you obtain the matrix? Thank you!
Thanks, this was eye opening
I have derived it myself
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
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
Hi,
good job. Everything in a single place 🙂
I suggest fixing two errors in the matrix (above the source code): diagonal elements, namely (2,2) and (3,3), contain terms u+w^2 and u+v^2 respectively. I guess it should be u^2 in both cases, i.e., u^2+w^2 and u^2 + v^2.
Adam