From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Change axes of ellipsoid
Date: Fri, 14 May 2010 22:14:05 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 32
Message-ID: <hskhvd$orn$>
References: <hsi8q2$jjo$>
Reply-To: <HIDDEN>
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: 1273875245 25463 (14 May 2010 22:14:05 GMT)
NNTP-Posting-Date: Fri, 14 May 2010 22:14:05 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1187260
Xref: comp.soft-sys.matlab:635992

"Els " <> wrote in message <hsi8q2$jjo$>...
> 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

  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