VerdictVector.hpp

Go to the documentation of this file.

/*=========================================================================

  Module:    $RCSfile: VerdictVector.hpp,v $

  Copyright (c) 2006 Sandia Corporation.
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/

/*
 *
 * VerdictVector.hpp contains declarations of vector operations
 *
 * This file is part of VERDICT
 *
 */

#ifndef VERDICTVECTOR_HPP
#define VERDICTVECTOR_HPP

#include "moab/verdict.h"
#include <cassert>
#include <cmath>

class VerdictVector;
typedef void ( VerdictVector::*transform_function )( double gamma, double gamma2 );
// a pointer to some function that transforms the point,
// taking a double parameter.  e.g. blow_out, rotate, and scale_angle

class VerdictVector
{
  public:
    //- Heading: Constructors and Destructor
    VerdictVector();  //- Default constructor.

    VerdictVector( const double x, const double y, const double z );
    //- Constructor: create vector from three components

    VerdictVector( const double xyz[3] );
    //- Constructor: create vector from tuple

    VerdictVector( const VerdictVector& tail, const VerdictVector& head );
    //- Constructor for a VerdictVector starting at tail and pointing
    //- to head.

    VerdictVector( const VerdictVector& copy_from );  //- Copy Constructor

    //- Heading: Set and Inquire Functions
    void set( const double xv, const double yv, const double zv );
    //- Change vector components to {x}, {y}, and {z}

    void set( const double xyz[3] );
    //- Change vector components to xyz[0], xyz[1], xyz[2]

    void set( const VerdictVector& tail, const VerdictVector& head );
    //- Change vector to go from tail to head.

    void set( const VerdictVector& to_copy );
    //- Same as operator=(const VerdictVector&)

    double x() const;  //- Return x component of vector

    double y() const;  //- Return y component of vector

    double z() const;  //- Return z component of vector

    void get_xyz( double& x, double& y, double& z );  //- Get x, y, z components
    void get_xyz( double xyz[3] );                    //- Get xyz tuple

    double& r();  //- Return r component of vector, if (r,theta) format

    double& theta();  //- Return theta component of vector, if (r,theta) format

    void x( const double xv );  //- Set x component of vector

    void y( const double yv );  //- Set y component of vector

    void z( const double zv );  //- Set z component of vector

    void r( const double xv );  //- Set r component of vector, if (r,theta) format

    void theta( const double yv );  //- Set theta component of vector, if (r,theta) format

    void xy_to_rtheta();
    //- convert from cartesian to polar coordinates, just 2d for now
    //- theta is in [0,2 PI)

    void rtheta_to_xy();
    //- convert from  polar to cartesian coordinates, just 2d for now

    void scale_angle( double gamma, double );
    //- tranform_function.
    //- transform  (x,y) to (r,theta) to (r,gamma*theta) to (x',y')
    //- plus some additional scaling so long chords won't cross short ones

    void blow_out( double gamma, double gamma2 = 0.0 );
    //- transform_function
    //- blow points radially away from the origin,

    void rotate( double angle, double );
    //- transform function.
    //- transform  (x,y) to (r,theta) to (r,theta+angle) to (x',y')

    void reflect_about_xaxis( double dummy, double );
    //- dummy argument to make it a transform function

    double normalize();
    //- Normalize (set magnitude equal to 1) vector - return the magnitude

    VerdictVector& length( const double new_length );
    //- Change length of vector to {new_length}. Can be used to move a
    //- location a specified distance from the origin in the current
    //- orientation.

    double length() const;
    //- Calculate the length of the vector.
    //- Use {length_squared()} if only comparing lengths, not adding.

    double distance_between( const VerdictVector& test_vector );
    //- Calculate the distance from the head of one vector
    //  to the head of the test_vector.

    double length_squared() const;
    //- Calculate the squared length of the vector.
    //- Faster than {length()} since it eliminates the square root if
    //- only comparing other lengths.

    double interior_angle( const VerdictVector& otherVector );
    //- Calculate the interior angle: acos((a%b)/(|a||b|))
    //- Returns angle in degrees.

    double vector_angle_quick( const VerdictVector& vec1, const VerdictVector& vec2 );
    //- Calculate the interior angle between the projections of
    //- {vec1} and {vec2} onto the plane defined by the {this} vector.
    //- The angle returned is the right-handed angle around the {this}
    //- vector from {vec1} to {vec2}. Angle is in radians.

    double vector_angle( const VerdictVector& vector1, const VerdictVector& vector2 ) const;
    //- Compute the angle between the projections of {vector1} and {vector2}
    //- onto the plane defined by *this. The angle is the
    //- right-hand angle, in radians, about *this from {vector1} to {vector2}.

    void perpendicular_z();
    //- Transform this vector to a perpendicular one, leaving
    //- z-component alone. Rotates clockwise about the z-axis by pi/2.

    void print_me();
    //- Prints out the coordinates of this vector.

    void orthogonal_vectors( VerdictVector& vector2, VerdictVector& vector3 );
    //- Finds 2 (arbitrary) vectors that are orthogonal to this one

    void next_point( const VerdictVector& direction, double distance, VerdictVector& out_point );
    //- Finds the next point in space based on *this* point (starting point),
    //- a direction and the distance to extend in the direction. The
    //- direction vector need not be a unit vector.  The out_point can be
    //- "*this" (i.e., modify point in place).

    bool within_tolerance( const VerdictVector& vectorPtr2, double tolerance ) const;
    //- Compare two vectors to see if they are spatially equal.  They
    //- compare if x, y, and z are all within tolerance.

    //- Heading: Operator Overloads  *****************************
    VerdictVector& operator+=( const VerdictVector& vec );
    //- Compound Assignment: addition: {this = this + vec}

    VerdictVector& operator-=( const VerdictVector& vec );
    //- Compound Assignment: subtraction: {this = this - vec}

    VerdictVector& operator*=( const VerdictVector& vec );
    //- Compound Assignment: cross product: {this = this * vec},
    //- non-commutative

    VerdictVector& operator*=( const double scalar );
    //- Compound Assignment: multiplication: {this = this * scalar}

    VerdictVector& operator/=( const double scalar );
    //- Compound Assignment: division: {this = this / scalar}

    VerdictVector operator-() const;
    //- unary negation.

    friend VerdictVector operator~( const VerdictVector& vec );
    //- normalize. Returns a new vector which is a copy of {vec},
    //- scaled such that {|vec|=1}. Uses overloaded bitwise NOT operator.

    friend VerdictVector operator+( const VerdictVector& v1, const VerdictVector& v2 );
    //- vector addition

    friend VerdictVector operator-( const VerdictVector& v1, const VerdictVector& v2 );
    //- vector subtraction

    friend VerdictVector operator*( const VerdictVector& v1, const VerdictVector& v2 );
    //- vector cross product, non-commutative

    friend VerdictVector operator*( const VerdictVector& v1, const double sclr );
    //- vector * scalar

    friend VerdictVector operator*( const double sclr, const VerdictVector& v1 );
    //- scalar * vector

    friend double operator%( const VerdictVector& v1, const VerdictVector& v2 );
    //- dot product

    friend VerdictVector operator/( const VerdictVector& v1, const double sclr );
    //- vector / scalar

    friend int operator==( const VerdictVector& v1, const VerdictVector& v2 );
    //- Equality operator

    friend int operator!=( const VerdictVector& v1, const VerdictVector& v2 );
    //- Inequality operator

    friend VerdictVector interpolate( const double param, const VerdictVector& v1, const VerdictVector& v2 );
    //- Interpolate between two vectors. Returns (1-param)*v1 + param*v2

    VerdictVector& operator=( const VerdictVector& from );

  private:
    double xVal;  //- x component of vector.
    double yVal;  //- y component of vector.
    double zVal;  //- z component of vector.
};

VerdictVector vectorRotate( const double angle, const VerdictVector& normalAxis, const VerdictVector& referenceAxis );
//- A new coordinate system is created with the xy plane corresponding
//- to the plane normal to {normalAxis}, and the x axis corresponding to
//- the projection of {referenceAxis} onto the normal plane.  The normal
//- plane is the tangent plane at the root point.  A unit vector is
//- constructed along the local x axis and then rotated by the given
//- ccw angle to form the new point.  The new point, then is a unit
//- distance from the global origin in the tangent plane.
//- {angle} is in radians.

inline double VerdictVector::x() const
{
    return xVal;
}
inline double VerdictVector::y() const
{
    return yVal;
}
inline double VerdictVector::z() const
{
    return zVal;
}
inline void VerdictVector::get_xyz( double xyz[3] )
{
    xyz[0] = xVal;
    xyz[1] = yVal;
    xyz[2] = zVal;
}
inline void VerdictVector::get_xyz( double& xv, double& yv, double& zv )
{
    xv = xVal;
    yv = yVal;
    zv = zVal;
}
inline double& VerdictVector::r()
{
    return xVal;
}
inline double& VerdictVector::theta()
{
    return yVal;
}
inline void VerdictVector::x( const double xv )
{
    xVal = xv;
}
inline void VerdictVector::y( const double yv )
{
    yVal = yv;
}
inline void VerdictVector::z( const double zv )
{
    zVal = zv;
}
inline void VerdictVector::r( const double xv )
{
    xVal = xv;
}
inline void VerdictVector::theta( const double yv )
{
    yVal = yv;
}
inline VerdictVector& VerdictVector::operator+=( const VerdictVector& vector )
{
    xVal += vector.x();
    yVal += vector.y();
    zVal += vector.z();
    return *this;
}

inline VerdictVector& VerdictVector::operator-=( const VerdictVector& vector )
{
    xVal -= vector.x();
    yVal -= vector.y();
    zVal -= vector.z();
    return *this;
}

inline VerdictVector& VerdictVector::operator*=( const VerdictVector& vector )
{
    double xcross, ycross, zcross;
    xcross = yVal * vector.z() - zVal * vector.y();
    ycross = zVal * vector.x() - xVal * vector.z();
    zcross = xVal * vector.y() - yVal * vector.x();
    xVal   = xcross;
    yVal   = ycross;
    zVal   = zcross;
    return *this;
}

inline VerdictVector::VerdictVector( const VerdictVector& copy_from )
    : xVal( copy_from.xVal ), yVal( copy_from.yVal ), zVal( copy_from.zVal )
{
}

inline VerdictVector::VerdictVector() : xVal( 0 ), yVal( 0 ), zVal( 0 ) {}

inline VerdictVector::VerdictVector( const VerdictVector& tail, const VerdictVector& head )
    : xVal( head.xVal - tail.xVal ), yVal( head.yVal - tail.yVal ), zVal( head.zVal - tail.zVal )
{
}

inline VerdictVector::VerdictVector( const double xIn, const double yIn, const double zIn )
    : xVal( xIn ), yVal( yIn ), zVal( zIn )
{
}

// This sets the vector to be perpendicular to it's current direction.
// NOTE:
//      This is a 2D function.  It only works in the XY Plane.
inline void VerdictVector::perpendicular_z()
{
    double temp = x();
    x( y() );
    y( -temp );
}

inline void VerdictVector::set( const double xv, const double yv, const double zv )
{
    xVal = xv;
    yVal = yv;
    zVal = zv;
}

inline void VerdictVector::set( const double xyz[3] )
{
    xVal = xyz[0];
    yVal = xyz[1];
    zVal = xyz[2];
}

inline void VerdictVector::set( const VerdictVector& tail, const VerdictVector& head )
{
    xVal = head.xVal - tail.xVal;
    yVal = head.yVal - tail.yVal;
    zVal = head.zVal - tail.zVal;
}

inline VerdictVector& VerdictVector::operator=( const VerdictVector& from )
{
    xVal = from.xVal;
    yVal = from.yVal;
    zVal = from.zVal;
    return *this;
}

inline void VerdictVector::set( const VerdictVector& to_copy )
{
    *this = to_copy;
}

// Scale all values by scalar.
inline VerdictVector& VerdictVector::operator*=( const double scalar )
{
    xVal *= scalar;
    yVal *= scalar;
    zVal *= scalar;
    return *this;
}

// Scales all values by 1/scalar
inline VerdictVector& VerdictVector::operator/=( const double scalar )
{
    assert( scalar != 0 );
    xVal /= scalar;
    yVal /= scalar;
    zVal /= scalar;
    return *this;
}

// Returns the normalized 'this'.
inline VerdictVector operator~( const VerdictVector& vec )
{
    double mag = sqrt( vec.xVal * vec.xVal + vec.yVal * vec.yVal + vec.zVal * vec.zVal );

    VerdictVector temp = vec;
    if( mag != 0.0 )
    {
        temp /= mag;
    }
    return temp;
}

// Unary minus.  Negates all values in vector.
inline VerdictVector VerdictVector::operator-() const
{
    return VerdictVector( -xVal, -yVal, -zVal );
}

inline VerdictVector operator+( const VerdictVector& vector1, const VerdictVector& vector2 )
{
    double xv = vector1.x() + vector2.x();
    double yv = vector1.y() + vector2.y();
    double zv = vector1.z() + vector2.z();
    return VerdictVector( xv, yv, zv );
    //  return VerdictVector(vector1) += vector2;
}

inline VerdictVector operator-( const VerdictVector& vector1, const VerdictVector& vector2 )
{
    double xv = vector1.x() - vector2.x();
    double yv = vector1.y() - vector2.y();
    double zv = vector1.z() - vector2.z();
    return VerdictVector( xv, yv, zv );
    //  return VerdictVector(vector1) -= vector2;
}

// Cross products.
// vector1 cross vector2
inline VerdictVector operator*( const VerdictVector& vector1, const VerdictVector& vector2 )
{
    return VerdictVector( vector1 ) *= vector2;
}

// Returns a scaled vector.
inline VerdictVector operator*( const VerdictVector& vector1, const double scalar )
{
    return VerdictVector( vector1 ) *= scalar;
}

// Returns a scaled vector
inline VerdictVector operator*( const double scalar, const VerdictVector& vector1 )
{
    return VerdictVector( vector1 ) *= scalar;
}

// Returns a vector scaled by 1/scalar
inline VerdictVector operator/( const VerdictVector& vector1, const double scalar )
{
    return VerdictVector( vector1 ) /= scalar;
}

inline int operator==( const VerdictVector& v1, const VerdictVector& v2 )
{
    return ( v1.xVal == v2.xVal && v1.yVal == v2.yVal && v1.zVal == v2.zVal );
}

inline int operator!=( const VerdictVector& v1, const VerdictVector& v2 )
{
    return ( v1.xVal != v2.xVal || v1.yVal != v2.yVal || v1.zVal != v2.zVal );
}

inline double VerdictVector::length_squared() const
{
    return ( xVal * xVal + yVal * yVal + zVal * zVal );
}

inline double VerdictVector::length() const
{
    return ( sqrt( xVal * xVal + yVal * yVal + zVal * zVal ) );
}

inline double VerdictVector::normalize()
{
    double mag = length();
    if( mag != 0 )
    {
        xVal = xVal / mag;
        yVal = yVal / mag;
        zVal = zVal / mag;
    }
    return mag;
}

// Dot Product.
inline double operator%( const VerdictVector& vector1, const VerdictVector& vector2 )
{
    return ( vector1.x() * vector2.x() + vector1.y() * vector2.y() + vector1.z() * vector2.z() );
}

#endif