Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Change axes of ellipsoid

Subject: Change axes of ellipsoid

From: Els

Date: 14 May, 2010 01:25:22

Message: 1 of 9

I want to rotate a standard ellipsoid, constructed by the next code:
-------------
% construct of center ellipsoid
a= rand(1,3)*10
cx=a(1)
cy=a(2)
cz=a(3)

% construct the three different radii
b= rand(1,3)*10
rx=b(1)
ry=b(2)
rz=b(3)

% construct ellipsoid
[x,y,z] = ellipsoid(cx,cy,cz,rx,ry,rz,25)
h=surf(x,y,z)
alpha(0.05)
hold on
-------------------------
The ellipsoid has 3 axes peripendicular to cx, cy, and cz.
But now I have a line, how doi I rotate my ellipsoid so that the old x axes is rotated to the line. And furthermore, the surface points (x,y,z) will change as well. How do I get the rotated surface points (xRot,yRot,zRot)

Subject: Change axes of ellipsoid

From: Roger Stafford

Date: 14 May, 2010 04:30:22

Message: 2 of 9

"Els " <y.e.t.reeuwijk@student.utwente.nl> wrote in message <hsi8q2$jjo$1@fred.mathworks.com>...
> I want to rotate a standard ellipsoid, constructed by the next code:
> -------------
> % construct of center ellipsoid
> a= rand(1,3)*10
> cx=a(1)
> cy=a(2)
> cz=a(3)
>
> % construct the three different radii
> b= rand(1,3)*10
> rx=b(1)
> ry=b(2)
> rz=b(3)
>
> % construct ellipsoid
> [x,y,z] = ellipsoid(cx,cy,cz,rx,ry,rz,25)
> h=surf(x,y,z)
> alpha(0.05)
> hold on
> -------------------------
> The ellipsoid has 3 axes peripendicular to cx, cy, and cz.
> But now I have a line, how doi I rotate my ellipsoid so that the old x axes is rotated to the line. And furthermore, the surface points (x,y,z) will change as well. How do I get the rotated surface points (xRot,yRot,zRot)
- - - - - - - - - -
  I assume your "line" runs through the origin. Let t be a 3D row vector pointing along this line in the direction toward which the positive x-axis will rotate. The arrays x, y, and z here are your arrays. Then do this:

% Calculate necessary rotation data
t = t/norm(t); % Make t a unit vector
u = [1,0,0]; % Unit vector along the positive x-axis
w = cross(u,t); % The axis of rotation
cosa = dot(u,t); % Cosine of the angle of rotation
sina = norm(w); % Sine of the angle of rotation
w = w/sina; % Make w a unit vector
v = cross(w,u); % Unit vector orthogonal to u and w

% Proceed with the rotation transformation
n = numel(x); % Number of points in surface
p = [x(:),y(:),z(:)]; % Create n x 3 array of ellipsoid surface points
u = repmat(u,n,1); % Match sizes of u, v, and w to that of p
v = repmat(v,n,1);
w = repmat(w,n,1);
wp = cross(w,p,2);
q = dot(p,w,2)*w + cosa*cross(wp,w,2) + sina*wp; % Rotate p
xr = q(:,1); yr = q(:,2); zr = q(:,3); % Extract components
xr = reshape(xr,size(x)); % These are the rotated coordinates
yr = reshape(yr,size(x)); % in mesh format again
zr = reshape(zr,size(x));

% Now plot the rotated surface
surf(xr,yr,zr)

  I'm sure this can be made a lot more streamlined, but it should get you there. Unless I have made some blunder - I haven't had time to test it. You can do that for me.

  Note: If your "line" doesn't run through the origin, you can first translate everything so that it does. Afterward the rotation you can translate back again by the same amount. I was too lazy to do that for you.

Roger Stafford

Subject: Change axes of ellipsoid

From: Els

Date: 14 May, 2010 05:50:37

Message: 3 of 9

Roger, thanks a lot!! This will get me there for sure. And I will fix the possible translation as well.

Subject: Change axes of ellipsoid

From: Roger Stafford

Date: 14 May, 2010 06:46:37

Message: 4 of 9

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <hsijku$ltv$1@fred.mathworks.com>...
> ........
> q = dot(p,w,2)*w + cosa*cross(wp,w,2) + sina*wp; % Rotate p
> ........

  I thought of one blunder in that code. In generating q the intent was that each row of w was to be multiplied by the single value in the corresponding element of the single column of dot(p,w,2), but the '*' operation won't do the job there. It ought to read

q = repmat(dot(p,w,2),1,3).*w + cosa*cross(wp,w,2) + sina*wp;

Roger Stafford

Subject: Change axes of ellipsoid

From: Roger Stafford

Date: 14 May, 2010 22:14:05

Message: 5 of 9

"Els " <y.e.t.reeuwijk@student.utwente.nl> wrote in message <hsi8q2$jjo$1@fred.mathworks.com>...
> I want to rotate a standard ellipsoid, constructed by the next code:
> ............
> But now I have a line, how do i I rotate my ellipsoid so that the old x axes is rotated to the line. And furthermore, the surface points (x,y,z) will change as well. How do I get the rotated surface points (xRot,yRot,zRot)
- - - - - - - - -
  I was not satisfied with that code I sent you, because it is too inefficient for large sets of points in your meshes, and because it does not handle cases well where the line and the x-axis are pointing in almost the same direction. Here is revised code which I think you will like much better. Also I have subjected it to a fair amount of testing. I generalize it here to a transformation which rotates any vector, u, by way of a "great circle" path, to any other vector, t. In your case the first of these, u, was to be the x-axis, but the following works for any u.

  Let u be a vector which is to be rotated via a great circle path to another vector t. Let x, y, and z be your mesh point coordinate arrays. Then do the following.

 % Prepare a 3 x 3 rotation transformation matrix, R
 t = t(:).'/norm(t); u = u(:).'/norm(u) ; % Make u & t unit row vectors
 w = cross(u,t); % w points along the axis of rotation
 cosa = u*t.'; % This gives scalar cosine of rotation angle
 I = eye(3); % The 3 x 3 identity matrix
 R = cosa*I + 1/(1+cosa)*w.'*w + cross(repmat(w,3,1),I,2);

 % Perform the rotation transformation on the points
 p = [x(:),y(:),z(:)]; % Rows of p are points' vectors
 q = p*R; % This rotates all the points
 xr = q(:,1); yr = q(:,2); zr = q(:,3); % Separate the coordinates
 xr = reshape(xr,size(x)); % Restore to original array sizes
 yr = reshape(yr,size(y));
 zr = reshape(zr,size(z));

 % Display rotated surface
 surf(xr,yr,zr)

  You will note that when the angle of rotation would be zero, the previous method would encounter a zero-divided-by-zero situation and would yield a NaN, but that is not true of the above algorithm. A NaN can be produced with this new method only when u and t are in completely opposite directions where cosa = -1, but that is an inherently indeterminate situation because the "great circle" route is not uniquely defined and the R matrix is therefore not unique.

  If timing is of any concern to you, note also that the step q = p*R should execute much faster than the corresponding step in the previous method.

Roger Stafford

Subject: Change axes of ellipsoid

From: Els

Date: 15 May, 2010 23:50:16

Message: 6 of 9

Roger, the previous script did work good. But you are so dedicated, this code does even work better. Thanks a lot for all your effort.

Subject: Change axes of ellipsoid

From: Els

Date: 25 May, 2010 12:15:21

Message: 7 of 9

Dear Roger,

Now, two weeks later, I detected a bug. And I hope you can help me with this.
The center of the original ellipsoid, should be the same as the rotated ellipsoid.
When I use your code, I get that the center of the rotated ellipsoid is translated.

Is there an easy way to correct this?

Subject: Change axes of ellipsoid

From: Roger Stafford

Date: 25 May, 2010 22:02:08

Message: 8 of 9

"Els " <y.e.t.reeuwijk@student.utwente.nl> wrote in message <htgf0p$s56$1@fred.mathworks.com>...
> Dear Roger,
>
> Now, two weeks later, I detected a bug. And I hope you can help me with this.
> The center of the original ellipsoid, should be the same as the rotated ellipsoid.
> When I use your code, I get that the center of the rotated ellipsoid is translated.
>
> Is there an easy way to correct this?
- - - - - - - -
  The code I sent you only does a rotation about an axis which runs through the origin. If your ellipsoid's center is not at the origin, it will very likely be displaced. You should first do a translation of your ellipsoid to center it at the origin. Then do the rotate. Then do a reverse translation. That should carry the center back to its original position, but the ellipsoid points around it will have been rotated about that center.

  I wouldn't call that a "bug". It's a capability assumed by the user (you) that isn't actually present. :-)

Roger Stafford

Subject: Change axes of ellipsoid

From: Els

Date: 26 May, 2010 05:23:04

Message: 9 of 9

Thanks, solved it!

Best wishes,

Els

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us