Generates the rototranslation matrix for the rotation around an arbitrary line in 3D. The line need not pass through the origin. Optionally, also, applies this transformation to a list of 3D coordinates.
SYNTAX 1:
M=AxelRot(deg,u,x0)
in:
u, x0: 3D vectors specifying the line in parametric form x(t)=x0+t*u
Default for x0 is [0,0,0] corresponding to pure rotation (no shift).
If x0=[] is passed as input, this is also equivalent to passing x0=[0,0,0].
deg: The counterclockwise rotation about the line in degrees. Counterclockwise is defined using the
right hand rule in reference to the direction of u.
out:
M: A 4x4 affine transformation matrix representing
the rototranslation. Namely, M will have the form
M=[R,t;0 0 0 1]
where R is a 3x3 rotation and t is a 3x1 translation vector.
SYNTAX 2:
[R,t]=AxelRot(deg,u,x0)
Same as Syntax 1 except that R and t are returned as separate arguments.
SYNTAX 3:
This syntax requires 4 input arguments be specified,
[XYZnew, R, t] = AxelRot(XYZold, deg, u, x0)
where the columns of the 3xN matrix XYZold specify a set of N point coordinates in 3D space. The output XYZnew is the transformation of the columns of XYZold, i.e., the original coordinates rotated appropriately about the axis. All other input/output arguments have the same meanings as before.
Matt J (2020). 3D Rotation about Shifted Axis (https://www.mathworks.com/matlabcentral/fileexchange/308643drotationaboutshiftedaxis), MATLAB Central File Exchange. Retrieved .
1.3.0.0  Added option to transform a set of coordinates. 

1.1.0.0  Clarified description and help doc only. No new code. 
Create scripts with code, output, and formatted text in a single executable document.
Zhao Li (view profile)
Matt J (view profile)
Hi Farzad,
You would use the 4th syntax
[XYZnew, R, t] = AxelRot(XYZold, deg, u, x0)
where u=cross(normal1,normal2) is a direction vector for the line of intersection and x0 is any point on that line.
Feri (view profile)
Sorry, the previous question should have been asked with my account, could you please answer it here.
How can I use your code for a 3D point cloud data, (XYZ), that should be rotated around the intersection of two planes in space?
ron garden (view profile)
Hi Matt,
How can I use your code for a 3D point cloud data, (XYZ), that should be rotated around the intersection of two planes in space?
Ali Morshedifard (view profile)
Ho Ha (view profile)
Got my answer!
Ho Ha (view profile)
Hi @Matt J,
I have a 3D data (x,y,z, and color). They initially have a cubeshape.
When I rotate the data around xaxis and point x0, I expect to see a cubeshape again, but the shape is different! For YZ plane, the shape is parallelogram (instead of rectangle)! Can you tell me what is going wrong?
I want to preserve the degree between planes.
D G (view profile)
Matt J (view profile)
Albert,
I'm afraid I don't follow your rationale for choosing deg=0. However, you must understand that deg specifies the amount of rotation. If you rotate everything by 0 degrees, then there is no rotation and nothing moves. So, it is no surprise that you always get M=eye(4). The argument u does not represent a target for the new zaxis. It is the axis you want to rotate everything around.
Albert (view profile)
@Matt J
I don't understand the 'deg' well. Let's say, I want to rotate coordinate system xyz such that z axis is parallel to line x=y=z, ie, u=(1,1,1). Actually xyz coordinate system doesn't rotate around the line x=y=z, thus I set deg to 0 degree.
deg=0;
u=[1, 1, 1];
x0=[0, 0, 0];
M=AxelRot(deg,u,x0)
>> M=
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
Similarly rotate xyz coordinate system such that z axis is parallel to line x=y=z, ie, u=(1,1,1).
deg=0;
u=[1, 1, 1];
x0=[0, 0, 0];
M=AxelRot(deg,u,x0)
>>M=
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
Obviously I shall not get same rotation matrix like this, right?
Please point out where I have done wrong.
Thank you very much for your code and time and energy to reply.
Matt J (view profile)
Hi tan,
That depends on which axis you want to rotate about. There are infinitely many choices, but for your first example one possible choice is u=cross([0 0 1],[1,1,1]) and similarly for your second example about u=cross([1 0 0],[2 3 0])
tan (view profile)
Hi Matt J,
How do I rotate a plane if it is initially pointing at the Z direction to point in the direction of X=1,Y=1 and Z=1 coordinate?
Alternatively, how do you rotate the plane from the existing Z direction to the coordinate at X=2 ,Y=3 and Z=0 coordinate. Thanks.
Sébastien Herbert (view profile)
Matt J (view profile)
Hi Tasnim,
You will need tformarray to rotate a 3D image, but you can compute the requisite tform using AxelRot
tasnim hasan (view profile)
Hi Matt J, Can I use this to rotate 3D MR data about an arbitrary and an angle
Matt J (view profile)
Hi Divya,
AxelRot does not rotate 2D or 3D images. It rotates coordinates. If using syntax 3, the Data must be 3xN whose columns form a list of 3D point coordinates.
Divya (view profile)
@Matt
The dimensions of "Data" are 1248x256x93 double
The error is the following:
Error using *
Inputs must be 2D, or at least one input must be scalar.
To compute elementwise TIMES, use TIMES (.*) instead.
Error in AxelRot (line 65)
XYZnew=bsxfun(@plus,R*XYZold,t);
Matt J (view profile)
@Divya
Sounds possible, but I'd need to see the error message and know the dimensions of "Data".
Divya (view profile)
I keep getting an error in line 65 about * usage when trying to run the code with a 3D matrix using syntax 3: AxelRot(Data, 5, [1 0 0], [ ]) . Is this an error due to my array sizes or something that I can fix. I am simply trying to rotate the dataset by a few degrees to be able to get particular features when I slice the data.
Fritz (view profile)
Matt J (view profile)
@Michael,
Ultimately, the problem you are having is not directly within the scope of what is offered here. This FEX submission is designed for coordinate rotation, rather than image rotation. Except that, of course, you could use AxelRot to generate input for maketform() which could then be used in imtransform(), tformarray(), etc...
The negative values you are getting from imrotate() is because cubic interpolation is an approximation of sinc interpolation. It therefore need not inherit positivity from the interpolated data points. It is typically a tradeoff you have to make in order to interpolate with C1 smoothness.
An alternative would be to use raised cosine interpolation, which is both C1continuous and positivitypreserving. However, it is an unconventional interpolation method and I can't be sure the rotated image would look the way you hope. Moreover, to implement this, you would have to use some advanced features of imtransfrom(). In particular, you would need to specify your own customized interpolation method using makeresampler().
Michael Young (view profile)
@Matt
My issue has been, for some reason, I am getting some negative values from an entirely positive matrix when using bicubic interpolation for imrotate(). I'm not sure why this is occuring, and I don't want to degrade the results by relying on bilinear.
Matt J (view profile)
@Michael,
It is not clear to me what distinction you draw between an NxN matrix and and NxN image.
Michael Young (view profile)
@Matt
I was under the impression you needed to upload an image to use imrotate(), instead of merely rotating a matrix. Is this incorrect? Thanks for the recommendation regardless!
Matt J (view profile)
@Michael,
Sounds like what you really need is imrotate() from the Image Processing Toolbox.
Michael Young (view profile)
Matt J:
I am working on a problem that I believe can utilize your technique, though I am not sure how exactly to do so. I wish to rotate an n x n matrix by x degrees, and then determine the appropriate values of the rotated matrix in regards to the original matrix. Visually, it will be a matrix overlapping with that same matrix rotated, and render the corresponding values on the original matrix from the rotated matrix. Thank you for any help you can provide!
Matt J (view profile)
Hi Evangelos,
You are correct. The purpose of AxelRot is to rotate about an axis that doesn't necessarily pass through the origin. As I think you'll find in the help doc, the inputs x0 and u define the desired rotation axis through the parametric equation x(t)=x0+t*u where x(t) is a point on the axis parametrized by t.
Evangelos Mazomenos (view profile)
Hi Matt,thanks for replying to my comment. After your reply and a little bit of more coding I think I know what is happening. Please verify my thoughts and correct me if I am wrong.
The AxelRot function you coded is function that rotates a point (or set of points) by "deg" about an axis. The axis is defined originally by "u" but the axis (of rotation) might be also shifted by "x0" so the term translation in your code comments refers to the shifting of the axis. My initial understanding of translation and the one that is implemented with the alternative code i wrote:
Mtrans = mkaff(eye(3), x0);
Mroto = mkaff(R3d(deg,u));
M = Mroto * Mtrans (considering translation first)
is that a point (or set of points) is initially translate by x0 (the point itself not the rotation axis) and then rotated by "deg" around u. This is also what MATLAB's inbuild function makehgttform implements. Basically this is independent transformations of points in 3D space (rotations, translations, scalings, etc ...)group together in a transformation matrix. On the other hand your intention is to provide a code only for rotation about an axis which might be shifted from his unit vector "u".
To support my claim I tried to rotate a point ([1 1 1]) for deg=360 and x0=[0 1 0] about axis u=[1 0 0].
With the initial AxelRot code the point since it is only rotated around an axis (albeit a shifted one by x0) for 360 it will return to its original location of [1 1 1]. On the other hand with the code (alternative code) in which x0 corresponds to a translation of the point in 3D the point will move to [1 2 1] as the result of the helical motion that corresponds to a rotation about an axis and a translation.
So I think both approaches (and codes) are correct but they implement a different operation. I got confused because I thought that "translation" in your comments corresponded to point translation and not axis shift. I must admit that you use the term shift in your code and I should be more careful.
Anyway let me know if you agree. Thanks for your time
Matt J (view profile)
Hi Evangelos, and thanks for your feedback. There is no specific literature reference that the code was built on. I derived the expressions on my own. However, I imagine most books on robotics, for example, would teach the principles behind deriving all the M matrices.
The alternative code that you have presented gives the wrong result, so I'm not sure what motivates it. You can see this by testing it with simple input data such as
x0=[0 1 0].'; u=[1 0 0 ].'; deg=90;
We know that the transformation specified by these values for x0,u,deg should rotate the origin to the point [0;1;1]. However, your code maps the origin to [0;0;1].
One thing I will mention though is that the choice of the value for AxisShift is not unique. I could also have chosen AxisShift=x0 or any other vector on the desired axis, u. You will see by testing them directly that all such choices produce the same result. However, x0(x0.'*u).*u is the shortest vector from the origin to the axis u and so this choice is somehow "cleaner" in my mind.
Hope that helps.
Evangelos Mazomenos (view profile)
Hi Matt J.
The code is really good thanks for sharing. Is it possible to provide a reference and briefly explain the way you construct the transformation matrix in your code and particularly the translation part. This is the part I do not understand why needs to be done:
AxisShift=x0(x0.'*u).*u; %%
Mshift=mkaff(eye(3),AxisShift); Mroto=mkaff(R3d(deg,u)); %% I get this
M=(Mshift\Mroto)*Mshift;
Why is the AxisShift variable calculated and why Mshift and M are defined like this??
As far I know the 4x4 transformation matrix M, for rotation about u for deg degrees and translation x0 (in 3D) should be defined as(given you code):
Mtrans = mkaff(eye(3), x0);
Mroto = mkaff(R3d(deg,u));
M = Mroto * Mtrans (considering translation first)
This is what matlab gives with the makehgttform command and I do not understand why you do it in a different way.
Thanks for your time.
Matt J (view profile)
Thanks, Karl. Yes, the code is from the days before I was better educated about backslash, but for such a small matrix, it is unlikely to matter.
Karl (view profile)
Hi Matt,
Very helpful script!
Matlab suggests replacing 'M=inv(Mshift)*Mroto*Mshift' with the likes of 'M=(Mshift\Mroto)*Mshift' for speed and accuracy (it will also remove the warning from the editor and code analysis).
Karl (view profile)
Matt J (view profile)
@Michael
Not sure I understand your question. You can rototranslate one or more 3D column vectors using SYNTAX 3 of the routine.
Michael (view profile)
How could you modify this code to rotate/translate a 3D velocity vector instead of just a coordinate frame?
Sagar (view profile)
Alan Jennings (view profile)
Good work. I'd recomend adding the copyright restrictions to the source code to keep them together.
Matt J (view profile)
Hi Florian,
I don't/can't know why you're getting a different result from your CAD software, but I'm pretty convinced that P3_CATIA is incorrect.
One check you can do is to find the perpendicular projection of the different points onto the rotation axis, as with the code below. They should all agree, but I get a significant discrepancy for prj_CATIA
axproj=@(z) u*(u\(zx0))+x0;
prj_old=axproj(P3_old),
prj_new=axproj(P3new_1),
prj_CATIA=axproj(P3_CATIA),
Once you've computed the perpendicular projection it's also easy to verify the angular separation between P3_old and P3new_1, as with the code below. I get 5 degrees to very high precision
anglesep=@(a,b)acosd(dot(a,b)/norm(a)/norm(b));
angle=anglesep(P3new_1prj_new,P3_oldprj_new),
Florian (view profile)
Hi Matt,
thanks for the code. But I've a little problem.
I tried to use your function with this points
P1 = [33.319 862.139 420.373]
P2 = [550.303 721.267 667.915]
P3 = [113 870.591 483.66]
%% Variablendeklaration
deg = 5;
k = [P1(1,1) P2(1,1)
P1(1,2) P2(1,2)
P1(1,3) P2(1,3) ];
x0 = [P1(1,1);P1(1,2);P1(1,3)];
x1 = [P2(1,1);P2(1,2);P3(1,3)];
%%
u = x1  x0;
%% test
P3_old = [P3(1,1);P3(1,2);P3(1,3)]
%%
[P3new_1, R, t] = AxelRot(P3_old, deg, u, x0)
the result is
P3new_1 =
111.4677
866.0073
485.9746
I'm making the same rotation in our CAD software (CATIA)> I'll get this result
P3_CATIA =
111.386
868.555
Do you have any idea about the differences?
485.872
Matt J (view profile)
@insa: I don't know what you mean exactly by a "circular node". Do you mean a single point in 3D space? Is each "node" defined simply by an (x,y,z) coordinate triple which you want to rotate to a different location? If so, this is a straightforward application of SYNTAX 3 of the tool and you shouldn't have any problem.
insa lyon (view profile)
Dear Mr.Matt J,
simply, we have a circular nodes lying on xy plane, i want to rotate theses nodes with a specific angle and get the rotated circular points and their coordinates ??what do you think?
Regards,
Matt J (view profile)
Hmmm. No, I guess not. It's sounding less and less like a rigid transformation.
insa lyon (view profile)
Dear Mr.Matt J,Thank you for your answer.
Actually, i have a series of '3D circular' coordinates( Nodes) lying on a straight line , this is the initial form.What i need is to propagate every single point(node)through the curved centerline points.
So, the final form i need will be a series of 3D circular coordinates( Nodes)lying on a curved centerline points , these 3D circular coordinates have to be perpendicular on this curved line ,so you think that i can transform those coordinates using SYNTAX 3?
My kind regards,
Matt J (view profile)
@insa: your question is unclear to me, I'm afraid. You mean the data you have is a series of 3D coordinates lying on some kind of spiral space curve? If so, you can transform those coordinates using SYNTAX 3.
insa lyon (view profile)
Hello,
i have the axis of a Cylinder straight,and the curved axis alongthat i want to calculate the new curved cylinder withe roto translation?
Do you thing that i can do it ??
Melissa (view profile)
Matt J (view profile)
@Melissa  finding the axis of your cylinder, if you don't already know it, is a data fitting problem. You should probably pose that question in the Newsgroup, since it doesn't concern this FEX submission directly.
Melissa (view profile)
How do I find a tilted axis so that I can use this?
Matt J (view profile)
@Melissa  You can use SYNTAX 3 to undo the tilt in your cylinder, assuming you already know its tilted zaxis, Ztilt. The cylinder becomes untilted when you rotate it about the rotation axis
u=cross(Ztilt/norm(Ztilt), [0 0 1])
by a rotation angle of deg=asind(u).
Melissa (view profile)
Can I use this to take a set of 3d cylinder cartesian coordinates and make a profile of polar coordinates by each slice of the cylinder?
Check this: http://www.mathworks.com/matlabcentral/answers/13398cartesiantopolar for more details of what I need.
Matt J (view profile)
@pink, you might want to explain to me what you're trying to do. I can't tell what elaboration of the help documentation you're looking for.
pink (view profile)
how to use?
thanks