#include <geometry.hpp>

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 Vector3d &	getTranslation () const
bool	hasRot () const
bool	hasInversion () const
double	getTheta () const
const Vector3d &	getAxis () 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

enum Transform::mat_format

Enumerator:

C_STYLE
FORTRAN_STYLE

Definition at line 114 of file geometry.hpp.

{ C_STYLE, FORTRAN_STYLE };

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

Vector3d Transform::axis [protected]

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.

Vector3d Transform::translation [protected]

Definition at line 117 of file geometry.hpp.

The documentation for this class was generated from the following files:

Public Types

Public Member Functions

Protected Member Functions

Protected Attributes

Detailed Description

Member Enumeration Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation