/* /usr/local/gismo/repo/geometry/Volume.cc,v 1.36 1994/08/11 22:42:54 alanb Exp */ // File: Volume.cc // Author: Alan Breakstone // Feb-06-1995 :: Paul Fuchs - fixes in the calculation of roll,pitch,yaw angles // Contents ---------------------------------------------------------- // // Volume::Volume() // Volume::~Volume() // Volume::Volume( const Volume& v ) // operator<<( ostream& os, const Volume& v ) // Volume::printOn( ostream& os ) const // Volume::operator==( const Volume& v ) // Volume::distanceToLeave( const Ray& r, // ThreeVec& p, const Surface*& s ) const // Volume::distanceToEnter( const Ray& r, // ThreeVec& p, const Surface*& s ) const // Volume::distanceToLeave( const Ray& r, // ThreeVec& p, const Surface*& s, double testDist ) const // Volume::distanceToEnter( const Ray& r, // ThreeVec& p, const Surface*& s, double testDist ) const // Volume::inside( const ThreeVec& x ) const // Volume::move( const ThreeVec& p ) // Volume::rotate( double alpha, double beta, double gamma, int inverse ) // Volume::rotateLocalToGlobal( const ThreeVec& v ) // Volume::rotateGlobalToLocal( const ThreeVec& v ) // Volume::localToGlobal( const ThreeVec& v ) // Volume::globalToLocal( const ThreeVec& v ) // Volume::localToGlobal( Ray* r ); // Volume::globalToLocal( Ray* r ); // Volume::localToGlobal( Surface* s ) // Volume::globalToLocal( Surface* s ) // Volume::localToGlobal( Volume* v ) // Volume::globalToLocal( Volume* v ) // Volume::removeSurfaces() // Volume::createSurfaces() // Volume::setDefaults() // Volume::read( ipstream& str ) // Volume::write( opstream& str ) // Volume::numberOfParameters() // Volume::getParameterName( int i ) // Volume::getParameterType( int ) // Volume::getParameter( int i ) // Volume::setParameter( int i, void* value ) // Volume::checkParameters() // Volume::drawYourSelf( View3D* v ) // Volume::createPolyLine( ThreeVec v[], int n ) // Volume::createPolyLine( float x[],float y[],float z[], int n ) // Volume::makeA( const char* name ) // // End --------------------------------------------------------------- #ifdef __GNUC__ #pragma implementation #endif #include "Volume.h" //static StreamableLink vfl( new Volume() ); // create link to factory Volume::Volume() { // default constructor // creates an empty list and gives a default title to the Volume surface_list = new GmsList(); setTitle( "myVolume" ); // center point and orientation angles set to zero // to change these use the move and rotate member functions center = ThreeVec( 0.0, 0.0, 0.0 ); roll = 0.0; pitch = 0.0; yaw = 0.0; // Set the maximum dimension to a large number. calcMaxDimension(); } Volume::~Volume() { // destructor removes all Surfaces from its list and deletes them removeSurfaces(); delete surface_list; } Volume::Volume( const Volume& v ) { // copy constructor // must make deep copy (i.e., use the Surface copy constructor to make // real copies of the Surfaces to add to the GmsList) surface_list = new GmsList(); setTitle( v.getTitle() ); center = v.center; roll = v.roll; pitch = v.pitch; yaw = v.yaw; max_dimension = v.max_dimension; Surface* s = (Surface *)v.surface_list->first(); while ( s ) { Surface* sc = new Surface( *s ); surface_list->append( sc ); s = (Surface *)s->next(); } } ostream& operator<<( ostream& os, const Volume& v ) { // overwrite output operator << to print out Volume objects // using the printOn function defined below v.printOn( os ); return os; } void Volume::printOn( ostream& os ) const { // printing function using C++ ostream class os << "Title of this Volume: " << title << "\n"; os << "Volume centered on " << center << "\n" << " roll angle " << roll << "\n" << " pitch angle " << pitch << "\n" << " yaw angle " << yaw << "\n" << " maximum dimension " << max_dimension << "\n" << "This Volume is bounded by the following Surfaces:\n"; Surface *s = (Surface *)surface_list->first(); while ( s ) { os << *s; s = (Surface *)s->next(); } } int Volume::operator==( const Volume& v ) { return ( center == v.center && roll == v.roll && pitch == v.pitch && yaw == v.yaw ); } double Volume::distanceToLeave( const Ray& r, ThreeVec& p, const Surface *&sf ) const { // Distance along the Ray r or any of its derived classes // to the nearest Surface to leave the Volume. // Gives the point p of intercept and a pointer to Surface which is hit. // If there is no intersection, the Surface pointer is set to 0. // This point is not guaranteed to be inside (or outside) the Volume. double d = 0.0, t = FLT_MAX; ThreeVec temp ( t, t, t ); p = temp; sf = 0; Surface *s = (Surface *)surface_list->first(); // Loop over each Surface in the Volume's list, // find the distance to that Surface, and pick the smallest positive // (or zero) distance. // The distanceToLeaveSurface member function of Ray will make sure // that the Ray or its derived class is leaving the Surface (defined by // the direction of the normal to the Surface) at the point of intersection. while ( s ) { d = r.distanceToLeaveSurface( s, temp ); if ( ( t > d ) && ( d >= 0.0 ) ) { t = d; p = temp; sf = s; } s = (Surface *)s->next(); } return t; } double Volume::distanceToEnter( const Ray& r, ThreeVec& p, const Surface *&sf ) const { // Distance along the Ray r or any of its derived classes // to the nearest Surface to enter the Volume. // Gives the point p of intercept and a pointer to Surface which is hit. // If there is no intersection, the Surface pointer is set to 0. // This point is not guaranteed to be inside (or outside) the Volume. double d = 0.0, t = FLT_MAX; ThreeVec temp ( t, t, t ); p = temp; sf = 0; Surface *s = (Surface *)surface_list->first(); // Loop over each Surface in the Volume's list, // find the distance to that Surface, and pick the smallest positive // (or zero) distance. // The distanceToEnterSurface member function of Ray will make sure // that the Ray or its derived class is entering the Surface (defined by // the direction of the normal to the Surface) at the point of intersection. while ( s ) { d = r.distanceToEnterSurface( s, temp ); if ( ( t > d ) && ( d >= 0.0 ) ) { t = d; p = temp; sf = s; } s = (Surface *)s->next(); } return t; } double Volume::distanceToLeave( const Ray& r, ThreeVec& p, const Surface *&sf, double testDist ) const { // Distance along the Ray r or any of its derived classes // to the nearest Surface to leave the Volume, constrained to be // less than testDist. // Gives the point p of intercept and a pointer to Surface which is hit. // If there is no intersection, the Surface pointer is set to 0. // This point is not guaranteed to be inside (or outside) the Volume. // If no intersection point is found, returns a large number (FLT_MAX). double d = 0.0, thn = FLT_MAX, t = FLT_MAX; ThreeVec temp ( t, t, t ); p = temp; sf = 0; Surface *s = (Surface *)surface_list->first(); // Loop over each Surface in the Volume's list. // Use the howNear function to find the closest distance between the // origin of the Ray (or its derived class) and the Surface. // Only if that distance is non-negative, less than the test distance, // and less than the distance along the Ray found for the previous Surface, // find the distance along the Ray (or its derived class) to that Surface. // This will get the shortest distance along the Ray to leave the Volume. // The distanceToLeaveSurface member function of Ray will make sure // that the Ray or its derived class is leaving the Surface (defined by // the direction of the normal to the Surface) at the point of intersection. while ( s ) { d = s->howNear( r.position() ); if ( fabs( d ) < thn ) thn = fabs( d ); if ( d >= 0.0 ) { if ( ( d < testDist ) && ( d < t ) ) { d = r.distanceToLeaveSurface( s, temp ); if ( ( t > d ) && ( d >= 0.0 ) ) { t = d; p = temp; sf = s; } } } s = (Surface *)s->next(); } // Return a negative distance in certain cases as a flag to the calling // routine in the medium subproject. if ( thn < FLT_MAX && t == FLT_MAX && testDist < FLT_MAX / 2 ) return -thn; else return t; } double Volume::distanceToEnter( const Ray& r, ThreeVec& p, const Surface *&sf, double testDist ) const { // Distance along the Ray r or any of its derived classes // to the nearest Surface to enter the Volume, constrained to be // less than testDist. // Gives the point p of intercept and a pointer to Surface which is hit. // If there is no intersection, the Surface pointer is set to 0. // This point is not guaranteed to be inside (or outside) the Volume. // If no intersection point is found, returns a large number (FLT_MAX). double d = 0.0, t = FLT_MAX; ThreeVec temp ( t, t, t ); p = temp; sf = 0; // Quick test to see if the Ray can possibly enter this Volume. // If the distance from the origin of the Ray to the center of the Volume // minus the maximum dimension (which is the radius of a sphere which // circumscribes the Volume) is greater than testDist, then the Ray cannot // enter the Volume. ThreeVec oc = center - r.position(); if ( ( oc.magnitude() - max_dimension ) > testDist ) return t; // Surface *s = (Surface *)surface_list->first(); // Loop over each Surface in the Volume's list. // Use the howNear function to find the closest distance between the // origin of the Ray (or its derived class) and the Surface. // For some types of Surfaces, howNear returns a negative distance for // points outside the Surface, so one must take the absolute value. // Only if that distance is non-negative, less than the test distance, // and less than the distance along the Ray found for the previous Surface, // find the distance along the Ray (or its derived class) to that Surface. // This will get the shortest distance along the Ray to enter the Volume. // The distanceToEnterSurface member function of Ray will make sure // that the Ray or its derived class is entering the Surface (defined by // the direction of the normal to the Surface) at the point of intersection. while ( s ) { d = fabs( s->howNear( r.position() ) ); if ( ( d >= 0.0 ) && ( d <= testDist ) && ( d < t ) ) { d = r.distanceToEnterSurface( s, temp ); if ( ( t > d ) && ( d >= 0.0 ) ) { t = d; p = temp; sf = s; } } s = (Surface *)s->next(); } return t; } int Volume::inside( const ThreeVec& x ) const { // Return 0 if a point x is outside the Volume, 1 if inside. // The point is inside the Volume if it is inside all the Surfaces // which define the Volume. // Surfaces provide an outward-pointing unit normal for this purpose. int in = 0; Surface *s = (Surface *)surface_list->first(); while ( s ) { in = s->inside( x ); if ( !in ) break; s = (Surface *)s->next(); } return in; } void Volume::move( const ThreeVec& p ) { // translate Volume by p // do this by translating all Surfaces in the Volume's list // and translate the center point center += p; Surface *s = (Surface *)surface_list->first(); while ( s ) { s->move( p ); s = (Surface *)s->next(); } // delete the representation of this Volume for the drawing package // to force a new representation to be created deleteRepresentation(); } void Volume::rotate( double alpha, double beta, double gamma, int inverse ) { // rotate Volume first about global x-axis by angle alpha, // second about global y-axis by angle beta, // and third about global z-axis by angle gamma // by rotating all the Surfaces in the Volume list // angles are assumed to be given in radians // if inverse is non-zero, the order of rotations is reversed ThreeMat m; Surface *s = (Surface *)surface_list->first(); while ( s ) { s->rotate( alpha, beta, gamma, m, inverse ); s = (Surface *)s->next(); } // rotate the center of the volume center.rotate( alpha, beta, gamma, inverse ); // calculate new orientation angles appropriate for all classes // unit vectors in x and y directions ThreeVec n1 ( 1.0, 0.0, 0.0 ); ThreeVec n2 ( 0.0, 1.0, 0.0 ); // rotate the unit vectors first by the old orientation angles n1.rotate( roll, pitch, yaw, 0 ); n2.rotate( roll, pitch, yaw, 0 ); // now rotate them by the new angles n1.rotate( alpha, beta, gamma, inverse ); n2.rotate( alpha, beta, gamma, inverse ); // now calculate new roll, pitch, and yaw angles based on the axes // rotate first about global x-axis = roll angle // second about global y-axis = pitch angle // third about global z-axis = yaw angle pitch = -asin( n1.get_z() ); // pi/2 ... -pi/2 double cosp = cos( pitch ); if ( fabs( cosp ) < FLT_EPSILON ) { // special case for pitch angle = + or - PI/2 // for this case, only the sum (for -PI/2) or difference (for +PI/2) of // the yaw and roll angles is calculable, thus we set the yaw angle // arbitrarily to zero, and calculate only the roll angle // yaw = 0.0; // if ( n1.get_z() < 0.0 ) // roll = asin( n2.get_x() ); // else // roll = -asin( n2.get_x() ); // } // --- Feb 06, 1995 Paul Fuchs -- there is a parity error here for angles // larger than |pi/2|. Do it the other way around : Set the roll angle // arbitrarily to zero and calculate the yaw angle of the final rotation. // n1.z = 1.0 n3 has become -n1, n2 remains n2 // n1.z = -1.0 n3 has become n1, n2 remains n2 // ===> yaw ranges from -pi to pi (like the return value of phi()). roll = 0.0; ThreeVec nx = ThreeVec(n2.get_y(),-n2.get_x(),0.0); yaw = nx.phi(); } // general case, calculate both yaw and roll angles else { // There are two possible solutions, check which one is correct with // the x-component of the rotated n1 unit vector. // yaw = asin( n1.get_y() / cosp ); // roll = asin( n2.get_z() / cosp ); // double ax = cosp * cos(yaw); // if ( fabs( ax - n1.get_x() ) > FLT_EPSILON ) { // pitch = M_PI - pitch; // cosp = cos( pitch ); // yaw = asin( n1.get_y() / cosp ); // roll = asin( n2.get_z() / cosp ); // } // --- Feb 06, 1995 Paul Fuchs -- this is unneccesary. The goal is to define // a set of angles which are consistent with globalToLocal and localToGlobal // operations. // ===> roll ranges from -pi to pi // ===> pitch ranges from -pi/2 to pi/2 // ===> yaw ranges from -pi to pi // This solves a parity problem in the transformation between local and // global coordinate systems. yaw = n1.phi(); // selfconistent phi n2.rotate( 0.0, -pitch, -yaw, 1 ); // undo pitch and yaw ThreeVec nx = ThreeVec( n2.get_y(),n2.get_z(),0.0 ); roll = nx.phi(); // roll is left. } // delete the representation of this Volume for the drawing package // to force a new representation to be created deleteRepresentation(); } ThreeVec Volume::rotateLocalToGlobal( const ThreeVec& v ) { // Given the direction vector v whose components are given in the local // coordinate system of the Volume (defined by roll, pitch, and yaw // angles), calculate and return its components in the global coordinate // system. // First do rotations by roll about the global x-axis, // second by pitch about the global y-axis, and // third by yaw about the global z-axis. // Do not translate the vector. ThreeVec q = v; q.rotate_x( roll ); q.rotate_y( pitch ); q.rotate_z( yaw ); return q; } ThreeVec Volume::rotateGlobalToLocal( const ThreeVec& v ) { // Given the direction vector v whose components are given in the global // coordinate system, calculate and return its components in the local // coordinate system of the Volume (defined by roll, pitch, and yaw // angles). // Do rotations by -yaw about the global z-axis, followed // by -pitch about the global y-axis, and finally // by -roll about the global x-axis. // Do not translate the vector. ThreeVec q = v; q.rotate_z( -yaw ); q.rotate_y( -pitch ); q.rotate_x( -roll ); return q; } ThreeVec Volume::localToGlobal( const ThreeVec& v ) { // Given the vector v whose components are given in the local // coordinate system of the Volume (defined by center, roll, pitch, // and yaw), calculate and return its components in the global system. // First do rotations by roll about the global x-axis, // second by pitch about the global y-axis, and // third by yaw about the global z-axis. // Then translate the vector by center. ThreeVec q = v; q.rotate_x( roll ); q.rotate_y( pitch ); q.rotate_z( yaw ); q += center; return q; } ThreeVec Volume::globalToLocal( const ThreeVec& v ) { // Given the vector v whose components are given in the global system, // calculate and return its components in the local coordinate system // of the Volume (defined by center, roll, pitch, and yaw). // First translate the vector by -center, // then do rotations by -yaw about the global z-axis, followed // by -pitch about the global y-axis, and finally // by -roll about the global x-axis. ThreeVec q = v; q -= center; q.rotate_z( -yaw ); q.rotate_y( -pitch ); q.rotate_x( -roll ); return q; } void Volume::localToGlobal( Ray* r ) { // Transforms a Ray or Helix from the local coordinate system of the // Volume (defined by center, roll, pitch, and yaw) to the global system. // The origin of the Ray or Helix is both rotated and translated, but the // direction vectors are only rotated. // First do rotations by roll about the global x-axis, // second by pitch about the global y-axis, and // third by yaw about the global z-axis, // then do the translation. r->transform( center, roll, pitch, yaw, 0 ); } void Volume::globalToLocal( Ray* r ) { // Transforms a Ray or Helix to the local coordinate system of the // Volume (defined by center, roll, pitch, and yaw) from the global // system. The origin of the Ray or Helix is both rotated and translated, // but the direction vectors are only rotated. // First translate by -center, then // do rotations by -yaw about the global z-axis, then // by -pitch about the global y-axis, and finally // by -roll about the global x-axis. ThreeVec c = -center; r->transform( c, -roll, -pitch, -yaw, 1 ); } void Volume::localToGlobal( Surface* s ) { // Transforms a Surface from the local coordinate system of the Volume // (defined by center, roll, pitch, and yaw) to the global system. // What this means is that any vector data member of the Surface will // have its coordinates transformed from the local to global system. // First do rotations by roll about the global x-axis, // second by pitch about the global y-axis, and // third by yaw about the global z-axis, using the Surface::rotate // function. Then translate using the Surface::move function. ThreeMat m; s->rotate( roll, pitch, yaw, m, 0 ); s->move( center ); } void Volume::globalToLocal( Surface* s ) { // Transforms a Surface to the local coordinate system of the Volume // (defined by center, roll, pitch, and yaw) from the global system. // What this means is that any vector data member of the Surface will // have its coordinates transformed from the global to local system. // First translate by -center using the Surface::move function, then // do rotations by -yaw about the global z-axis, then // by -pitch about the global y-axis, and finally // by -roll about the global x-axis, using the Surface::rotate // function with the inverse flag set to 1. ThreeMat m; s->move( -center ); s->rotate( -roll, -pitch, -yaw, m, 1 ); } void Volume::localToGlobal( Volume* v ) { // Transforms a Volume from the local coordinate system of this Volume // (defined by center, roll, pitch, and yaw) to the global system. // What this means is that any vector data member of the Volume will // have its coordinates transformed from the local to global system. // First do rotations by roll about the global x-axis, // second by pitch about the global y-axis, and // third by yaw about the global z-axis, using the Volume::rotate // function. Then translate using the Volume::move function. v->rotate( roll, pitch, yaw, 0 ); v->move( center ); } void Volume::globalToLocal( Volume* v ) { // Transforms a Volume to the local coordinate system of this Volume // (defined by center, roll, pitch, and yaw) from the global system. // What this means is that any vector data member of the Volume will // have its coordinates transformed from the global to local system. // First translate by -center using the Volume::move function, then // do rotations by -yaw about the global z-axis, then // by -pitch about the global y-axis, and finally // by -roll about the global x-axis, using the Volume::rotate // function with the inverse flag set to 1. v->move( -center ); v->rotate( -roll, -pitch, -yaw, 1 ); } void Volume::removeSurfaces() { // remove surfaces from list and delete them // typically followed by polymorphic createSurfaces() Surface *s = (Surface *)surface_list->first(); while ( s ) { surface_list->remove( s ); delete s; s = (Surface *)surface_list->first(); } } void Volume::createSurfaces() { // use to make surface list consistent with position and orientation // following reading in or change of parameters // remove any surfaces that may have existed in Volume's list if ( surface_list->count() > 0 ) removeSurfaces(); // make sure this is a subclass of Volume, if not, quit if ( !strcmp( nameOf(),"Volume" ) )return; // call the particular surface creator createSurfaces(); // check if the Surfaces created need to be rotated or translated int needToRotate = ( roll != 0.0 || pitch != 0.0 || yaw != 0.0 ); int needToTranslate = ( center != ThreeVec(0.0,0.0,0.0) ); if ( ! needToRotate && ! needToTranslate ) return; // if so, go through the list and rotate and translate them as necessary ThreeMat m; Surface *s = (Surface *)surface_list->first(); while ( s ) { if ( needToRotate ) s->rotate( roll, pitch, yaw, m, 0 ); if ( needToTranslate ) s->move(center); s = (Surface *)s->next(); } // also must rotate vector data members if ( needToRotate ) rotateVectors( m ); } void Volume::setDefaults() { // set default parameters for a Volume // this is equivalent to destroying a Volume and creating one with // the null constructor center = ThreeVec( 0.0, 0.0, 0.0 ); roll = pitch = yaw = 0.0; removeSurfaces(); } // ------------------------------------------------------------------ // I/O functions /* void* Volume::read( ipstream& str ) { // read in title, center coordinates, and orientation angles from // input stream double x,y,z; // if anyone ever makes title a char *, enable the following #if 0 str >> title; #else char * t; str >> t; strncpy(title,t,sizeof(title)-1); delete t; #endif str >> x >> y >> z >> roll >> pitch >> yaw; center = ThreeVec(x,y,z); return this; } void Volume::write( opstream& str ) { // write out title, center coordinates, and orientation angles to // output stream str << title << center.get_x() << center.get_y() << center.get_z() << roll << pitch << yaw ; } */ // ------------------------------------------------------------ // Parameter functions int Volume::numberOfParameters() { return 6; } const char* Volume::getParameterName( int i ) { static char* parname[] = { "x", "y", "z", "roll", "pitch", "yaw" }; return parname[i]; } ParameterType Volume::getParameterType( int ) { // all of Volume's parameters are double precision return PAR_TYPE_DOUBLE; } const void* Volume::getParameter( int i ) { static double temp; switch(i) { case 0: return &(temp= center.get_x()); case 1: return &(temp= center.get_y()); case 2: return &(temp= center.get_z()); case 3: return &roll; case 4: return &pitch; case 5: return &yaw; default: return 0; } } void Volume::setParameter( int i, void* value ) { double v = *(double*)value; switch(i) { case 0: center = ThreeVec( v, center.get_y(), center.get_z() ); break; case 1: center = ThreeVec( center.get_x(), v, center.get_z() ); break; case 2: center = ThreeVec( center.get_x(), center.get_y(), v ); break; case 3: roll = v; break; case 4: pitch = v; break; case 5: yaw = v; break; default: return; // should not happen } } const char* Volume::checkParameters() { // this is meant to generate an error message if parameters are // illegal removeSurfaces(); return 0; } // ------------------------------------------------------------ // Drawing functions void Volume::drawYourSelf( View3D* v ) { // override the drawYourSelf from Viewable in order to call the // polymorphic createRepresentation for the first display if( !hasRepresentation() ) createRepresentation(); Viewable::drawYourSelf( v ); } void Volume::createPolyLine( ThreeVec v[], int n ) { // for use by the createRepresentation methods moveTo( v[0] ); for ( int i = 1; i < n; i++ ) lineTo( v[i] ); flush(); setGrayLevel( 0.0 ); } void Volume::createPolyLine( float x[],float y[],float z[], int n ) { // for use by the createRepresentation methods moveTo( ThreeVec( x[0], y[0], z[0] ) ); for ( int i = 1; i < n; i++ ) lineTo( ThreeVec( x[i], y[i], z[i] ) ); flush(); setGrayLevel( 0.0 ); } /* // // static function for object persistence Volume* Volume::makeA( const char* name ) { // Interface to the Streamable class data base Volume* v = (Volume*)Streamable::makeA(name); if (!v || strcmp(v->nameOfBase(),"Volume")) return 0; return v; } // // Overloaded friend functions opstream& operator << ( opstream& os, Volume& v ) { v.write(os); return os; } opstream& operator << ( opstream& os, Volume* v ) { v->write(os); return os; } ipstream& operator >> ( ipstream& is, Volume& v ) { v.read(is); return is; } ipstream& operator >> ( ipstream& is, Volume* v ) { v->read(is); return is; } */