cgma
Transform Class Reference

#include <geometry.hpp>

List of all members.

Public Types

enum  mat_format { C_STYLE, FORTRAN_STYLE }

Public Member Functions

 Transform ()
 Transform (const Vector3d &v)
 Transform (const std::vector< double > &inputs, bool degree_format_p=false, enum mat_format=FORTRAN_STYLE)
 Transform (double rot[9], const Vector3d &trans, enum mat_format=C_STYLE)
const Vector3dgetTranslation () const
bool hasRot () const
bool hasInversion () const
double getTheta () const
const Vector3dgetAxis () const
void print (std::ostream &str) const
Transform reverse () const

Protected Member Functions

void set_rots_from_matrix (double raw_matrix[9], enum mat_format)

Protected Attributes

Vector3d translation
bool has_rot
double theta
Vector3d axis
bool invert

Detailed Description

Definition at line 111 of file geometry.hpp.


Member Enumeration Documentation

Enumerator:
C_STYLE 
FORTRAN_STYLE 

Definition at line 114 of file geometry.hpp.


Constructor & Destructor Documentation

Transform::Transform ( ) [inline]

Definition at line 125 of file geometry.hpp.

:translation(),has_rot(false),invert(false){}
Transform::Transform ( const Vector3d v) [inline]

Definition at line 126 of file geometry.hpp.

:translation(v),has_rot(false),invert(false){}
Transform::Transform ( const std::vector< double > &  inputs,
bool  degree_format_p = false,
enum mat_format  f = FORTRAN_STYLE 
)

Definition at line 146 of file geometry.cpp.

                                                                                                  : 
  has_rot(false), invert(false)
{
  
  size_t num_inputs = inputs.size();
  
  // translation is always defined by first three inputs
  translation = Vector3d(inputs); 

  if( num_inputs == 9 ||                         // translation matrix with third vector missing
      num_inputs == 12 || num_inputs == 13 )  // translation matrix fully specified
    {
    
      has_rot = true;
      double raw_matrix[9];
    
    if( num_inputs == 9 ){
      for( int i = 3; i < 9; ++i){
        raw_matrix[i-3] = degree_format_p ? cos(inputs.at(i) * M_PI / 180.0 ) : inputs.at(i);
      }
      
      Vector3d v1( raw_matrix ); //v1 = v1.normalize();
      Vector3d v2( raw_matrix+3 ); //v2 = v2.normalize();
      Vector3d v3 = v1.cross(v2);//.normalize();
      raw_matrix[6] = v3.v[0];
      raw_matrix[7] = v3.v[1];
      raw_matrix[8] = v3.v[2];
    }
    else{
      for( int i = 3; i < 12; ++i){
        raw_matrix[i-3] = degree_format_p ? cos(inputs.at(i) * M_PI / 180.0 ) : inputs.at(i);
      }
      if( num_inputs == 13 && inputs.at(12) == -1.0 ){
        std::cout << "Notice: a transformation has M = -1.  Inverting the translation;" << std::endl;
        std::cout << " though this might not be what you wanted." << std::endl;
        translation = -translation;
      }
    }

    set_rots_from_matrix( raw_matrix, f );
    
  }
  else if( num_inputs != 3 ){
    // an unsupported number of transformation inputs
    std::cerr << "Warning: transformation with " << num_inputs << " input items is unsupported" << std::endl;
    std::cerr << "  (will pretend there's no rotation: expect incorrect geometry.)" << std::endl;
  }
  
}  
Transform::Transform ( double  rot[9],
const Vector3d trans,
enum mat_format  f = C_STYLE 
)

Definition at line 140 of file geometry.cpp.

                                                                              :
  translation(trans), has_rot(true), invert(false)
{
  set_rots_from_matrix( rot, f );
}

Member Function Documentation

const Vector3d& Transform::getAxis ( ) const [inline]

Definition at line 134 of file geometry.hpp.

{ return axis; }
double Transform::getTheta ( ) const [inline]

Definition at line 133 of file geometry.hpp.

{ return theta; }
const Vector3d& Transform::getTranslation ( ) const [inline]

Definition at line 130 of file geometry.hpp.

{ return translation; }
bool Transform::hasInversion ( ) const [inline]

Definition at line 132 of file geometry.hpp.

{ return invert; }
bool Transform::hasRot ( ) const [inline]

Definition at line 131 of file geometry.hpp.

{ return has_rot; }
void Transform::print ( std::ostream &  str) const

Definition at line 206 of file geometry.cpp.

                                            {
  str << "[trans " << translation;
  if(has_rot){
    str << "(" << theta << ":" << axis << ")";
  }
  if(invert){
    str << "(I)";
  }
  str << "]";
}
Transform Transform::reverse ( void  ) const

Definition at line 196 of file geometry.cpp.

                                   {
  Transform t;
  t.translation = -this->translation;
  t.has_rot = this->has_rot;
  t.invert = this->invert;
  t.axis = -this->axis;
  t.theta = this->theta;
  return t;
}
void Transform::set_rots_from_matrix ( double  raw_matrix[9],
enum mat_format  f 
) [protected]

Compute Euler axis/angle, given a rotation matix. See en.wikipedia.org/wiki/Rotation_representation_(mathematics)

Definition at line 26 of file geometry.cpp.

                                                                             {
    
  double mat[3][3]  = {{ raw_matrix[0], raw_matrix[3], raw_matrix[6] },
                       { raw_matrix[1], raw_matrix[4], raw_matrix[7] },
                       { raw_matrix[2], raw_matrix[5], raw_matrix[8] } };


  if( f == C_STYLE ){

    /* row-major ordering: rewrite mat as
       mat = { { raw_matrix[0], raw_matrix[1], raw_matrix[2] }; 
               { raw_matrix[3], raw_matrix[4], raw_matrix[5] },
               { raw_matrix[6], raw_matrix[7], raw_matrix[8] } };
    */

    mat[0][1] = raw_matrix [1];
    mat[0][2] = raw_matrix [2];
    mat[1][0] = raw_matrix [3];
    mat[1][2] = raw_matrix [5];
    mat[2][0] = raw_matrix [6];
    mat[2][1] = raw_matrix [7];

  }

  double det = matrix_det(raw_matrix); // determinant is the same regardless of ordering

  if( OPT_DEBUG ){
    std::cout << "Constructing rotation: " << std::endl;
    for( int i = 0; i < 3; i++ ){
      std::cout << "  [ ";
      for ( int j = 0; j < 3; j++ ){
        std::cout << mat[i][j] << " ";
      }
      std::cout << "]" << std::endl;
    }
    std::cout << "  det = " << det << std::endl;
  }

  if( det < 0.0 ){
    // negative determinant-> this transformation contains a reflection.
    invert = true;
    det *= -1;
    for( int i = 0; i < 9; i++){  mat[i/3][i%3] = -mat[i/3][i%3]; } 
    if( OPT_DEBUG ) std::cout << "  negative determinant => improper rotation (adding inversion)" << std::endl;
  }
  
  if( fabs( det - 1.0 ) > DBL_EPSILON ){
    std::cout << "Warning: determinant of rotation matrix " << det << " != +-1" << std::endl;
  }

  /* Older, more straightforward approach:
    theta = acos( (mat[0][0] + mat[1][1] + mat[2][2] - 1) / 2 );
    double twoSinTheta = 2 * sin(theta);
    axis.v[0] = ( mat[2][1]-mat[1][2]) / twoSinTheta;
    axis.v[1] = ( mat[0][2]-mat[2][0]) / twoSinTheta;
    axis.v[2] = ( mat[1][0]-mat[0][1]) / twoSinTheta;
  */

  /* I have switched from the simple implementation above to the more robust and complex
   * approach below.  It has better numerical stability and (what is more important)
   * handles correctly the cases where theta is a multiple of 180 degrees.
   * See also
   * en.wikipedia.org/wiki/Rotation_matrix#Axis_and_angle (for why to use atan2 instead of acos)
   * www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToAngle/index.htm (for handling 180-degree case)
   */

  double x = mat[2][1]-mat[1][2];
  double y = mat[0][2]-mat[2][0];
  double z = mat[1][0]-mat[0][1];
  double r = hypot( x, hypot( y,z ));
  double t = mat[0][0] + mat[1][1] + mat[2][2];
  theta = atan2(r,t-1);
  
  if( OPT_DEBUG ){
    std:: cout << "  r = " << r << " t = " << t << " theta = " << theta << std::endl;
    std:: cout << "  x = " << x << " y = " << y << " z = " << z << std::endl;
  }

  if( std::fabs(theta) <= DBL_EPSILON ){ 
    // theta is 0 or extremely close to it, so let's say there's no rotation after all.
    has_rot = false; 
    axis = Vector3d(); // zero vector
    if( OPT_DEBUG ) std::cout << "  (0) "; // std::endl comes below
  }
  else if( std::fabs( theta - M_PI ) <= DBL_EPSILON ){
    // theta is pi (180 degrees) or extremely close to it
    // find the column of mat with the largest diagonal
    int col = 0;
    if( mat[1][1] > mat[col][col] ){ col = 1; }
    if( mat[2][2] > mat[col][col] ){ col = 2; }

    axis.v[col] = sqrt( (mat[col][col]+1)/2 ); 
    double denom = 2*axis.v[col];
    axis.v[(col+1)%3] = mat[col][(col+1)%3] / denom;
    axis.v[(col+2)%3] = mat[col][(col+2)%3] / denom;

    if( OPT_DEBUG ) std::cout << "  (180) "; // std::endl comes below

  }
  else{ 
    // standard case: theta isn't 0 or 180
    axis.v[0] = x / r;
    axis.v[1] = y / r;
    axis.v[2] = z / r;
  }


  // convert theta back to degrees, as used by iGeom
  theta *= 180.0 / M_PI;
  
  if( OPT_DEBUG ) std::cerr << "computed rotation: " << *this << std::endl;

}

Member Data Documentation

Definition at line 119 of file geometry.hpp.

bool Transform::has_rot [protected]

Definition at line 118 of file geometry.hpp.

bool Transform::invert [protected]

Definition at line 120 of file geometry.hpp.

double Transform::theta [protected]

Definition at line 119 of file geometry.hpp.

Definition at line 117 of file geometry.hpp.


The documentation for this class was generated from the following files:
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines