This is tm3rotate.c in view mode; [Download] [Up]
/* Copyright (c) 1992 The Geometry Center; University of Minnesota
1300 South Second Street; Minneapolis, MN 55454, USA;
This file is part of geomview/OOGL. geomview/OOGL is free software;
you can redistribute it and/or modify it only under the terms given in
the file COPYING, which you should have received along with this file.
This and other related software may be obtained via anonymous ftp from
geom.umn.edu; email: software@geom.umn.edu. */
/* Authors: Charlie Gunn, Pat Hanrahan, Stuart Levy, Tamara Munzner,
Mark Phillips, Nathaniel Thurston */
#include <math.h>
#include "radians.h"
#include "transform3.h"
HPoint3 TM3_XAXIS = { 1.0, 0.0, 0.0, 0.0 };
HPoint3 TM3_YAXIS = { 0.0, 1.0, 0.0, 0.0 };
HPoint3 TM3_ZAXIS = { 0.0, 0.0, 1.0, 0.0 };
void
Tm3RotateX( T, angle )
Transform3 T;
float angle;
{
Tm3Identity( T );
Ctm3RotateX( T, angle );
}
void
Tm3RotateY( T, angle )
Transform3 T;
float angle;
{
Tm3Identity( T );
Ctm3RotateY( T, angle );
}
void
Tm3RotateZ( T, angle )
Transform3 T;
float angle;
{
Tm3Identity( T );
Ctm3RotateZ( T, angle );
}
/*
* Construct the matrix for a rotation about an axis.
*/
void
Tm3Rotate( T, angle, axis)
Transform3 T;
float angle;
HPoint3 *axis;
{
HPoint3 Vu;
float sinA, cosA, versA;
if( axis == &TM3_XAXIS ) Tm3RotateX( T, angle );
else if( axis == &TM3_YAXIS ) Tm3RotateY( T, angle );
else if( axis == &TM3_ZAXIS ) Tm3RotateZ( T, angle );
else {
Pt3Copy( (Point3 *)axis, (Point3 *)&Vu );
Pt3Unit( (Point3 *)&Vu );
sinA = sin(angle);
cosA = cos(angle);
versA = 1 - cosA;
Tm3Identity( T );
T[X][X] = Vu.x*Vu.x*versA + cosA;
T[Y][X] = Vu.x*Vu.y*versA - Vu.z*sinA;
T[Z][X] = Vu.x*Vu.z*versA + Vu.y*sinA;
T[X][Y] = Vu.y*Vu.x*versA + Vu.z*sinA;
T[Y][Y] = Vu.y*Vu.y*versA + cosA;
T[Z][Y] = Vu.y*Vu.z*versA - Vu.x*sinA;
T[X][Z] = Vu.x*Vu.z*versA - Vu.y*sinA;
T[Y][Z] = Vu.y*Vu.z*versA + Vu.x*sinA;
T[Z][Z] = Vu.z*Vu.z*versA + cosA;
}
}
/*
* Construct the matrix for the geodesic rotation between two vectors.
*/
void
Tm3RotateBetween( T, vfrom, vto )
Transform3 T;
Point3 *vfrom, *vto;
{
Point3 Vu;
float len, sinA, cosA, versA;
Tm3Identity( T );
len = sqrt(Pt3Dot(vfrom,vfrom) * Pt3Dot(vto,vto));
if(len == 0) return;
cosA = Pt3Dot(vfrom,vto) / len;
versA = 1 - cosA;
Pt3Cross( vfrom, vto, &Vu );
sinA = Pt3Length( &Vu ) / len;
if(sinA == 0) return;
Pt3Mul( 1/(len*sinA), &Vu, &Vu ); /* Normalize Vu */
T[X][X] = Vu.x*Vu.x*versA + cosA;
T[Y][X] = Vu.x*Vu.y*versA - Vu.z*sinA;
T[Z][X] = Vu.x*Vu.z*versA + Vu.y*sinA;
T[X][Y] = Vu.y*Vu.x*versA + Vu.z*sinA;
T[Y][Y] = Vu.y*Vu.y*versA + cosA;
T[Z][Y] = Vu.y*Vu.z*versA - Vu.x*sinA;
T[X][Z] = Vu.x*Vu.z*versA - Vu.y*sinA;
T[Y][Z] = Vu.y*Vu.z*versA + Vu.x*sinA;
T[Z][Z] = Vu.z*Vu.z*versA + cosA;
}
/*-----------------------------------------------------------------------
* Function: Tm3RotateTowardZ
* Description: Rotation of 3-space moving pt to positive z-axis
* Args: T: output matrix
* pt: input point
* Returns: nonthing
* Author: njt (comments by mbp)
* Date: (before) Tue Aug 18 23:31:58 1992
*/
void Tm3RotateTowardZ(T, pt)
Transform3 T;
HPoint3 *pt;
{
Transform3 S;
float r = pt->z;
/* Construct T = rotation about x-axis moving pt into x-z plane */
Tm3Identity(T);
r = sqrt(pt->y*pt->y + r*r);
if (r > 0) {
T[2][1] = -(T[1][2] = pt->y/r);
T[2][2] = T[1][1] = pt->z/r;
}
/* Construct S = rotation about y axis moving T(pt) into y-z plane */
Tm3Identity(S);
r = sqrt(pt->x*pt->x + r*r);
if (r > 0) {
S[2][0] = -(S[0][2] = pt->x/r);
S[2][2] = S[0][0] = sqrt(pt->z*pt->z + pt->y*pt->y)/r;
}
/* Desired transform is then S * T */
Tm3Concat(T, S, T);
}
/* Tm3CarefulRotateTowardZ gives a matrix which rotates the world
* about an axis perpinducular to both pos and the z axis, which moves
* pos to the z axis. Unlike Tm3RotateTowardZ, the "twist" is well-defined
* provided that it is not asked to rotate the negative z axis toward
* the positive z axis.
*/
void Tm3CarefulRotateTowardZ(T, pos)
Transform3 T;
HPoint3 *pos;
{
Transform3 S;
Transform3 Sinv;
HPoint3 perp, axis;
double dist, c, s;
/* first, find a rotation takes both pos and the z axis into the xy plane */
perp.x = -pos->y; perp.y = pos->x; perp.z = 0; perp.w = 1;
Tm3RotateTowardZ(S, &perp);
/* now, rotate pos and the z axis to this plane */
axis.x = 0; axis.y = 0; axis.z = -1; axis.w = 1;
HPt3Transform(S, &axis, &axis);
HPt3Transform(S, pos, pos);
/* find the rotation matrix for the transformed points */
c = axis.x * pos->x + axis.y * pos->y;
s = axis.x * pos->y - axis.y * pos->x;
dist = sqrt(c*c+s*s);
Tm3Identity(T);
if (dist == 0) return;
c /= dist;
s /= dist;
T[0][0] = c; T[0][1] = s;
T[1][0] = -s; T[1][1] = c;
/* Finally, conjugate the result back to the original coordinate system */
Tm3Invert(S, Sinv);
Tm3Concat(S, T, T);
Tm3Concat(T, Sinv, T);
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.