[comp.graphics] cartesian to spherical

mcooper@sunc5.cs.uiuc.edu (09/29/90)

The answer's staring me in the face for this one, but I seem to have gone blind.
Keep in mind that this is on a SGI PI, so the fhypot etc. are valid for
floats..

I'm trying to write a routine that will draw a cylinder in an arbitary 
position in an arbitrary orientation.  I have a cylinder object created, with
one end centered at the origin, and the other at (1.0, 0.0, 0.0).

the routine should accept two endpoints, rotate the cylinder to the
proper angle, and translate it to the proper position.

 
I've hit a dead end in my efforts.  

I pass the routine two endpoints, (x1,y1,z1) and (x2,y2,z2).

to get my phi and theta angles, (and rho, for length), I use the following 
equations with x = x2-x1, y = y2-y1, z = z2 - z1:  


        rho = fhypot(x, fhypot(y, z));
        theta = fatan2(z, x);
        phi = facos2(y, *rho);

I have also tried 
	phi = fatan2(y, fhypot(z,x));


But nothing seems to work quite right.

Incidentally, I'm performing the operations as follows:

pushmatrix();
	translate(x2,y2,z2);
	rot(theta,'y');
	rot(phi,'z');
	scale(rho, 1.0, 1.0);
popmatrix();

The scaling and translation seem to work fine. the rotations are a mess.

Any help is appreciated.  Note the different email address in the .sig..
Thanks.

+-------------------------------------+---------------------------------------+
|  "When the going gets weird,        |  Marc Cooper	marcc@ncsa.uiuc.edu   |
|              the weird turn pro."   |                                       |
|		  		      |  National Center for Supercomputing   |
|		 -Hunter S. Thompson  |		    Applications	      |
+-------------------------------------+---------------------------------------+
|  The above opinoins are my own, and do not reflect those of NCSA.   *Yawn*  |
+-----------------------------------------------------------------------------+

erich@eye.com ( Eric Haines) (10/01/90)

In article <23900019@sunc5> mcooper@sunc5.cs.uiuc.edu writes:
>I have a cylinder object created, with
>one end centered at the origin, and the other at (1.0, 0.0, 0.0).
>
>the routine should accept two endpoints, rotate the cylinder to the
>proper angle, and translate it to the proper position.
>
> 
>I've hit a dead end in my efforts.  
>
>I pass the routine two endpoints, (x1,y1,z1) and (x2,y2,z2).
>
>to get my phi and theta angles, (and rho, for length), I use the following 
>equations with x = x2-x1, y = y2-y1, z = z2 - z1:  
>
>
>        rho = fhypot(x, fhypot(y, z));
>        theta = fatan2(z, x);
>        phi = facos2(y, *rho);
>
>I have also tried 
>	phi = fatan2(y, fhypot(z,x));
>
>
>But nothing seems to work quite right.
>
>Incidentally, I'm performing the operations as follows:
>
>pushmatrix();
>	translate(x2,y2,z2);
>	rot(theta,'y');
>	rot(phi,'z');
>	scale(rho, 1.0, 1.0);
>popmatrix();

Let's see, without looking too carefully at what you've done, the problem is how
to rotate axis {1,0,0} to some {x,y,z}.  Since you don't really care about the
amount you rotate about the axis of the cylinder (i.e. you're not trying to
line up some point on the cylinder to some new position - all you want is for
the cylinder to be between the points), the problem's easier:

	Scale the cylinder.
	Take the cross product of {x2-x1,y2-y1,z2-z1} and {1,0,0} - this is
		the axis you must rotate about (or you could do it in two
		rotations, as you tried).  Take the dot product of these two
		axes to get an angle of rotation about this axis.
	Now do the translate.

However, you seem to have nifty functions to do rotations about an orthogonal
axis quick.  So, OK, you want to rotate around z, then around y (if I
understand SGI notation correctly, i.e.  read from bottom to top).  To compute
the angle between the axis you start with and the axis you want to get to by
rotation around the Z axis, we cast both vectors onto the X-Y plane (happily,
{1,0,0} is already cast upon this plane) and take the dot product between
them.  The arcosine of this dot product is the "raw" angle between the
vectors.  So, the rotation from +X to the projection of {x2-x1,y2-y1,z2-z1}
(call this "A") on to the X-Y plane is:

	if ( y > 0.0 ) {
	    phi = facos2(x,fhypot(x,y)) ;
	} else {
	    phi = -facos2(x,fhypot(x,y)) ;
	}

since if x == 1 you want to rotate 0, x == 0 you want to rotate 90 deg
(assuming Y is non-zero), either plus or minus depending on the sign of y, and
x == -1 you rotate 180 deg.

Anyway, the idea is to find the amount of rotation for each axis by finding
the arcosine of the dot product between {1,0,0} and the "A" axis vector cast
upon the plane you want to rotate on (& normalizing this cast vector).  That
is, the "A" vector is {x,y,0} in the X-Y plane, and you want to normalize
this.  Explicitly:

	phi = arcosine( ( {x,y,0} / sqrt(x*x,y*y) ) DOT { 1, 0, 0 } )
    reduces to:
	phi = arcosine( x / sqrt(x*x,y*y) )

You need the "if" tests so that you can get the full -180 to 180 degree
rotations (facos2 gives 0 to 180 only, for example).

Then the Y rotation (in the new frame of reference) should be the angle
between the A vector cast onto the X-Y, normalized (which is where the
previous rotation we just computed rotated the {1,0,0} axis), and the full A
vector (the final position that we want).  As before, this is simply the
arcosine of the dot product of {x,y,0} and {x,y,z}, normalizing {x,y,0} before
doing so.  This works out to be:

	if ( z > 0.0 ) {
	    theta = facos(fhypot(x,y)) ;
	} else {
	    theta = -facos(fhypot(x,y)) ;
	}

Undoubtedly, I have a sign wrong here somewhere, but you should be able to
figure this out quick enough.


There are also some special cases which I assume facos2() covers gracefully,
i.e. a 0 denominator.

Hopefully this is right!

"Anything free comes with no guarantee",

Eric Haines, erich@eye.com