D:/Programmierung/C/TRTK/include/TRTK/Coordinate.hpp

Go to the documentation of this file.

/*
    Copyright (C) 2010 - 2017 Christoph Haenisch

    Chair of Medical Engineering (mediTEC)
    RWTH Aachen University
    Pauwelsstr. 20
    52074 Aachen
    Germany

    See license.txt for more information.

    Version 0.8.2 (2017-12-08)
*/

#ifndef COORDINATE_HPP_1023074891
#define COORDINATE_HPP_1023074891


#include <assert.h>
#include <cmath>
#include <stdexcept>
#include <string>

#include <Eigen/Core>
#include <Eigen/Dense>

#include "Tools.hpp"
#include "ErrorObj.hpp"

#ifdef TRTK_SUPPORT_CVECTOR
#include <Mathbib/Vector.h>
#endif // TRTK_SUPPORT_CVECTOR


namespace TRTK
{

// Forward declarations

namespace Tools
{
    template <class T>
    bool isZero(const T value);

    template <class T>
    T rand();

    template <class T>
    T rand(T a, T b);

    template <class T>
    T randn();

    template <class T>
    T randn(T mu, T sigma);
}

template <class T>
class Coordinate
{
public:
    typedef T value_type;
    typedef Eigen::Array<T, Eigen::Dynamic, 1> data_type;                   
    typedef Eigen::Matrix<T, Eigen::Dynamic, 1> matrix_type;                

    Coordinate();               
    Coordinate(T);              
    Coordinate(T, T);           
    Coordinate(T, T, T);        
    Coordinate(T, T, T, T);     

    Coordinate(const Coordinate &);

    template <class U>
    Coordinate(const Coordinate<U> & coordinate);

    #ifdef TRTK_SUPPORT_CVECTOR
    Coordinate(const CVector &);
    #endif // TRTK_SUPPORT_CVECTOR

    template <typename Derived>
    explicit Coordinate(const Eigen::DenseBase<Derived> & vector); 

    virtual ~Coordinate();

    T & x();                
    T & y();                
    T & z();                
    T & w();                

    const T & x() const;    
    const T & y() const;    
    const T & z() const;    
    const T & w() const;    

    T & operator[](int);
    const T & operator[](int) const;

    bool operator==(const Coordinate &) const;

    Coordinate & operator+=(T);
    Coordinate & operator-=(T);
    Coordinate & operator*=(T);
    Coordinate & operator/=(T);

    Coordinate & operator+=(const Coordinate &);        
    Coordinate & operator-=(const Coordinate &);        
    Coordinate & operator*=(const Coordinate &);        
    Coordinate & operator/=(const Coordinate &);        

    Coordinate operator+(T) const;
    Coordinate operator-(T) const;
    Coordinate operator*(T) const;
    Coordinate operator/(T) const;

    Coordinate operator+(const Coordinate &) const;     
    Coordinate operator-(const Coordinate &) const;     
    Coordinate operator*(const Coordinate &) const;     
    Coordinate operator/(const Coordinate &) const;     

    Coordinate operator,(T) const;                      

    operator data_type();                               
    operator data_type() const;                         

    #ifdef TRTK_SUPPORT_CVECTOR
    operator CVector() const;                           
    #endif // TRTK_SUPPORT_CVECTOR

    Coordinate cross(const Coordinate &) const;         
    double dot(const Coordinate &) const;               
    Coordinate & fill(T);
    double norm() const;                                
    Coordinate & normalize();                           
    Coordinate normalized() const;                      
    Coordinate orthogonal() const;                      
    Coordinate orthonormal() const;                     
    Coordinate & reserve(unsigned int size);            
    Coordinate & resize(unsigned int size, T value = T());
    unsigned int size() const;
    double squaredNorm() const;                         

    data_type & toArray();                              
    const data_type & toArray() const;                  

    matrix_type toMatrix() const;                       

    std::string toString() const;                       

    enum Error {DIVISION_BY_ZERO,                       
                UNKNOWN_ERROR                           
    };

    static Coordinate rand(unsigned int size, T a = T(0), T b = T(1));              
    static Coordinate randn(unsigned int size, T mu = T(0), T sigma =  T(1));       

private:
    data_type m_data;
};


// The subsequent functions allow expressions, where the scalar is in front
// of the coordinate. See Corrdinate.cpp for their implementations

Coordinate<double> operator+(int, const Coordinate<double> &);


Coordinate<double> operator-(int, const Coordinate<double> &);


Coordinate<double> operator*(int, const Coordinate<double> &);


Coordinate<double> operator/(int, const Coordinate<double> &);


Coordinate<int> operator+(double, const Coordinate<int> &);


Coordinate<int> operator-(double, const Coordinate<int> &);


Coordinate<int> operator*(double, const Coordinate<int> &);


Coordinate<int> operator/(double, const Coordinate<int> &);


template <class T>
Coordinate<T> operator+(T, const Coordinate<T> & coordinate);


template <class T>
Coordinate<T> operator-(T, const Coordinate<T> & coordinate);


template <class T>
Coordinate<T> operator*(T, const Coordinate<T> & coordinate);


template <class T>
Coordinate<T> operator/(T, const Coordinate<T> & coordinate);


// Stream operators

template <class T>
std::istream & operator >> (std::istream &, Coordinate<T> &);           


template <class T>
std::ostream & operator << (std::ostream &, const Coordinate<T> &);     



// Arithmetical functions

template <class T>
TRTK::Coordinate<T> sqrt(const TRTK::Coordinate<T> & coordinate);




template <class T>
Coordinate<T>::Coordinate()
{
}


template <class T>
Coordinate<T>::Coordinate(T x) : m_data(1)
{
    m_data(0) = x;
}


template <class T>
Coordinate<T>::Coordinate(T x, T y) : m_data(2)
{
    m_data(0) = x;
    m_data(1) = y;
}


template <class T>
Coordinate<T>::Coordinate(T x, T y, T z) : m_data(3)
{
    m_data(0) = x;
    m_data(1) = y;
    m_data(2) = z;
}


template <class T>
Coordinate<T>::Coordinate(T x, T y, T z, T w) : m_data(4)
{
    m_data(0) = x;
    m_data(1) = y;
    m_data(2) = z;
    m_data(3) = w;
}


template <class T>
Coordinate<T>::Coordinate(const Coordinate & coordinate)
{
    m_data = coordinate.m_data;
}


template <class T>
template <class U>
Coordinate<T>::Coordinate(const Coordinate<U> & coordinate)
{
    const unsigned int size = coordinate.size();

    m_data.resize(size);

    for (unsigned i = 0; i < size; ++i)
    {
        m_data(i) = static_cast<T>(coordinate[i]);
    }
}


#ifdef TRTK_SUPPORT_CVECTOR

template <class T>
Coordinate<T>::Coordinate(const CVector & cVector)
{
    const unsigned int size = cVector.GetSize();

    this->resize(size);

    for (int i = 0; i < size; ++i)
    {
        m_data[i] = cVector(i);
    }
}

#endif // TRTK_SUPPORT_CVECTOR


template <class T>
template <typename Derived>
Coordinate<T>::Coordinate(const Eigen::DenseBase<Derived> & vector)
{
    assert(vector.rows() == 1 || vector.cols() == 1);
    m_data = vector; // no conversion from row vector to column vector necessary
}


template <class T>
Coordinate<T>::~Coordinate()
{
}


template <class T>
inline T & Coordinate<T>::x()
{
    assert(size() >= 1);

    return m_data(0);
}


template <class T>
inline T & Coordinate<T>::y()
{
    assert(size() >= 2);

    return m_data(1);
}


template <class T>
inline T & Coordinate<T>::z()
{
    assert(size() >= 3);

    return m_data(2);
}


template <class T>
inline T & Coordinate<T>::w()
{
    assert(size() >= 4);

    return m_data(3);
}


template <class T>
inline const T & Coordinate<T>::x() const
{
    assert(size() >= 1);

    return m_data(0);
}


template <class T>
inline const T & Coordinate<T>::y() const
{
    assert(size() >= 2);

    return m_data(1);
}


template <class T>
inline const T & Coordinate<T>::z() const
{
   assert(size() >= 3);

   return m_data(2);
}


template <class T>
inline const T & Coordinate<T>::w() const
{
    assert(size() >= 4);

    return m_data(3);
}


template <class T>
inline T & Coordinate<T>::operator[](int index)
{
    assert(0 <= index && index < m_data.size());

    return m_data(index);
}


template <class T>
inline const T & Coordinate<T>::operator[](int index) const
{
    assert(0 <= index && index < m_data.size());

    return m_data(index);
}


template <class T>
inline bool Coordinate<T>::operator==(const Coordinate & coordinate) const
{
    assert(this->size() == coordinate.size());

    return (m_data == coordinate.m_data).all();
}


template <class T>
inline Coordinate<T> & Coordinate<T>::operator+=(T value)
{
    m_data += value;
    return *this;
}


template <class T>
inline Coordinate<T> & Coordinate<T>::operator-=(T value)
{
    m_data -= value;
    return *this;
}


template <class T>
inline Coordinate<T> & Coordinate<T>::operator*=(T value)
{
    m_data *= value;
    return *this;
}


template <class T>
Coordinate<T> & Coordinate<T>::operator/=(T value)
{
    if (value == 0)
    {
        ErrorObj error;
        error.setClassName("Coordinate");
        error.setFunctionName("operator/=");
        error.setErrorMessage("Division by zero.");
        error.setErrorCode(DIVISION_BY_ZERO);
        throw error;
    }

    m_data /= value;
    return *this;
}


template <class T>
inline Coordinate<T> & Coordinate<T>::operator+=(const Coordinate & coordinate)
{
    assert(this->size() == coordinate.size());

    m_data += coordinate.m_data;
    return *this;
}


template <class T>
inline Coordinate<T> & Coordinate<T>::operator-=(const Coordinate & coordinate)
{
    assert(this->size() == coordinate.size());

    m_data -= coordinate.m_data;
    return *this;
}


template <class T>
inline Coordinate<T> & Coordinate<T>::operator*=(const Coordinate & coordinate)
{
    assert(this->size() == coordinate.size());

    m_data *= coordinate.m_data;
    return *this;
}


template <class T>
Coordinate<T> & Coordinate<T>::operator/=(const Coordinate & coordinate)
{
    assert(this->size() == coordinate.size());

    if (!coordinate.m_data.all())
    {
        ErrorObj error;
        error.setClassName("Coordinate");
        error.setFunctionName("operator/=");
        error.setErrorMessage("Division by zero.");
        error.setErrorCode(DIVISION_BY_ZERO);
        throw error;
    }

    m_data /= coordinate.m_data;

    return *this;
}


template <class T>
Coordinate<T> Coordinate<T>::operator+(T value) const
{
    Coordinate<T> coordinate;
    coordinate.m_data = m_data + value;

    return coordinate;
}


template <class T>
Coordinate<T> Coordinate<T>::operator-(T value) const
{
    Coordinate<T> coordinate;
    coordinate.m_data = m_data - value;

    return coordinate;
}


template <class T>
Coordinate<T> Coordinate<T>::operator*(T value) const
{
    Coordinate<T> coordinate;
    coordinate.m_data = m_data * value;

    return coordinate;
}


template <class T>
Coordinate<T> Coordinate<T>::operator/(T value) const
{
    if (value == 0)
    {
        ErrorObj error;
        error.setClassName("Coordinate");
        error.setFunctionName("operator/");
        error.setErrorMessage("Division by zero.");
        error.setErrorCode(DIVISION_BY_ZERO);
        throw error;
    }

    Coordinate<T> coordinate;
    coordinate.m_data = m_data / value;

    return coordinate;
}


template <class T>
Coordinate<T> Coordinate<T>::operator+(const Coordinate & coordinate) const
{
    assert(this->size() == coordinate.size());

    Coordinate<T> tmp;
    tmp.m_data = m_data + coordinate.m_data;

    return tmp;
}


template <class T>
Coordinate<T> Coordinate<T>::operator-(const Coordinate & coordinate) const
{
    assert(this->size() == coordinate.size());

    Coordinate<T> tmp;
    tmp.m_data = m_data - coordinate.m_data;

    return tmp;
}


template <class T>
Coordinate<T> Coordinate<T>::operator*(const Coordinate & coordinate) const
{
    assert(this->size() == coordinate.size());

    Coordinate<T> tmp;
    tmp.m_data = m_data * coordinate.m_data;

    return tmp;
}


template <class T>
Coordinate<T> Coordinate<T>::operator/(const Coordinate & coordinate) const
{
    assert(this->size() == coordinate.size());

    if (!coordinate.m_data.all())
    {
        ErrorObj error;
        error.setClassName("Coordinate");
        error.setFunctionName("operator/");
        error.setErrorMessage("Division by zero.");
        error.setErrorCode(DIVISION_BY_ZERO);
        throw error;
    }

    Coordinate<T> tmp;
    tmp.m_data = m_data / coordinate.m_data;

    return tmp;
}


template <class T>
Coordinate<T> Coordinate<T>::operator,(T value) const
{
    const unsigned int size = m_data.size();
    data_type data(size + 1);
    data.head(size) = m_data;
    data(size) = value;

    Coordinate coordinate;
    coordinate.m_data = data;

    return coordinate;
}


template <class T>
Coordinate<T>::operator data_type()
{
    return m_data;
}


template <class T>
Coordinate<T>::operator data_type() const
{
    return m_data;
}


#ifdef TRTK_SUPPORT_CVECTOR

template <class T>
Coordinate<T>::operator CVector() const
{
    const unsigned int size = m_data.size();

    CVector cVector(size);

    for (unsigned int i = 0; i < size; ++i)
    {
        cVector(i) = m_data(i);
    }

    return cVector;
}

#endif // TRTK_SUPPORT_CVECTOR


template <class T>
Coordinate<T> Coordinate<T>::cross(const Coordinate & other) const
{
    assert((this->size() == 3) && (other.size() == 3));

    const T x1 = m_data(0);
    const T x2 = m_data(1);
    const T x3 = m_data(2);

    const T y1 = other.m_data(0);
    const T y2 = other.m_data(1);
    const T y3 = other.m_data(2);

    return Coordinate<T>(x2 * y3 - x3 * y2,
                         x3 * y1 - x1 * y3,
                         x1 * y2 - x2 * y1);
}


template <class T>
double Coordinate<T>::dot(const Coordinate & other) const
{
    assert(this->size() == other.size());

    return m_data.matrix().dot(other.m_data.matrix());
}


template <class T>
inline Coordinate<T> & Coordinate<T>::fill(T value)
{
    m_data.setConstant(value);
    return *this;
}


template <class T>
inline double Coordinate<T>::norm() const
{
    return m_data.matrix().norm();
}


template <class T>
Coordinate<T> & Coordinate<T>::normalize()
{
    m_data.matrix().normalize();
    return *this;
}


template <class T>
Coordinate<T> Coordinate<T>::normalized() const
{
    return Coordinate(m_data.matrix().normalized());
}


template <class T>
Coordinate<T> Coordinate<T>::orthogonal() const
{
    assert(this->size() >= 2);
    assert(!Tools::isZero(this->norm()));


    // Find largest and second largest element of this coordinate
    // and their according positions.

    T max[2] = {0};
    int pos[2] = {0};

    for (unsigned int i = 0; i < this->size(); ++i)
    {
        using std::abs;

        if (abs(m_data(i)) > abs(max[0]))
        {
            max[1] = max[0];
            max[0] = m_data(i);
            pos[1] = pos[0];
            pos[0] = i;
        }
        else if (abs(m_data(i)) >= abs(max[1]))
        {
            max[1] = m_data(i);
            pos[1] = i;
        }
    }

    // Construct the orthogonal vector by switching the positions of
    // the two largest elements and negating the second largest element.
    // All other components are set to zero.
    // The dot product of this constructed vector with *this always equals zero.

    Coordinate<T> result;
    result.resize(this->size(), 0);
    result[pos[0]] = - max[1];
    result[pos[1]] =   max[0];

    return result;
}


template <class T>
Coordinate<T> Coordinate<T>::orthonormal() const
{
    assert(this->size() >= 2);
    assert(!Tools::isZero(this->norm()));

    // Find largest and second largest element of this coordinate
    // and their according positions.

    T max[2] = {0};
    int pos[2] = {0};

    for (unsigned int i = 0; i < this->size(); ++i)
    {
        using std::abs;

        if (abs(m_data(i)) > abs(max[0]))
        {
            max[1] = max[0];
            max[0] = m_data(i);
            pos[1] = pos[0];
            pos[0] = i;
        }
        else if (abs(m_data(i)) >= abs(max[1]))
        {
            max[1] = m_data(i);
            pos[1] = i;
        }
    }

    // Construct the orthonormal vector by switching the positions of
    // the two largest elements and negating the second largest element.
    // All other components are set to zero.

    using std::sqrt;
    T norm = sqrt(max[0] * max[0] + max[1] * max[1]);

    Coordinate<T> result;
    result.resize(this->size(), 0);
    result[pos[0]] = - max[1] / norm;
    result[pos[1]] =   max[0] / norm;

    return result;
}


template <class T>
Coordinate<T> Coordinate<T>::rand(unsigned int size, T a, T b)
{
    Coordinate<T> coordinate;
    coordinate.reserve(size);

    for (unsigned i = 0; i < size; ++i)
    {
        coordinate[i] = Tools::rand(a, b);
    }

    return coordinate;
}


template <class T>
Coordinate<T> Coordinate<T>::randn(unsigned int size, T mu, T sigma)
{
    Coordinate<T> coordinate;
    coordinate.reserve(size);

    for (unsigned i = 0; i < size; ++i)
    {
        coordinate[i] = Tools::randn(mu, sigma);
    }

    return coordinate;
}


template <class T>
Coordinate<T> & Coordinate<T>::reserve(unsigned int size)
{
    m_data.resize(size);
    return *this;
}


template <class T>
Coordinate<T> & Coordinate<T>::resize(unsigned int size, T value)
{
    const int former_size = signed(this->size());
    const int difference_in_size = signed(size) - signed(former_size);

    // resizing allocates new memory and thus all data are lost
    // --> save data to tempory and restore it afterwards

    data_type data = m_data;

    m_data.resize(size);

    if (difference_in_size <= 0)
    {
        m_data = data.head(size);
    }
    else
    {
        m_data.head(former_size) = data;
        m_data.tail(difference_in_size).setConstant(value);
    }

    return *this;
}


template <class T>
inline unsigned int Coordinate<T>::size() const
{
    return static_cast<int>(m_data.size());
}


template <class T>
inline double Coordinate<T>::squaredNorm() const
{
    return m_data.matrix().squaredNorm();
}


template <class T>
inline typename Coordinate<T>::data_type & Coordinate<T>::toArray()
{
    return m_data;
}


template <class T>
inline const typename Coordinate<T>::data_type & Coordinate<T>::toArray() const
{
    return m_data;
}


template <class T>
inline typename Coordinate<T>::matrix_type Coordinate<T>::toMatrix() const
{
    return m_data;
}


template <class T>
std::string Coordinate<T>::toString() const
{
    std::stringstream ss;

    // throw an exception if the output operation fails for any reason

    ss.exceptions(std::ios_base::badbit | std::ios_base::failbit | std::ios_base::eofbit);

    // put out formatted coordinate

    ss // << std::fixed // or std::scientific
       // << std::setprecision(4)
       // << std::showpos
       // << std::setfill('0')
       // << std::internal
       << "(";

    int i = 0;

    while (i < signed(m_data.size()) - 1)
    {
        ss // << std::setw(10)
           << m_data(i++)
           << ", ";
    }

    if (m_data.size() != 0)
    {
        ss // << std::setw(10)
           << m_data(i);
    }

    ss <<  ")";

    return ss.str();
}




// Operators of the form T op Coordinate<class T>

template <class T>
Coordinate<T> operator+(T value, const Coordinate<T> & coordinate)
{
    return coordinate + value;
}


template <class T>
Coordinate<T> operator-(T value, const Coordinate<T> & coordinate)
{
    return coordinate * (-1) + value;
}


template <class T>
Coordinate<T> operator*(T value, const Coordinate<T> & coordinate)
{
    return coordinate * value;
}


template <class T>
Coordinate<T> operator/(T value, const Coordinate<T> & coordinate)
{
    Coordinate<T> coordinate2;
    coordinate2.resize(coordinate.size(), value);

    return coordinate2.operator/(coordinate);
}


// Stream operators

template <class T>
std::istream & operator >> (std::istream & input, Coordinate<T> & coordinate)
{
    // The input format is a comma separated list enclosed by round brackets.
    // Anything diverging from it, leads to an exception.
    //
    // Example: (1, 2, 3)

    char input_character;
    const std::string error_msg("Wrong input format or coordinate size.");

    // throw an exception if the input character differs from what was expected

    std::ios_base::iostate old_iostate = input.exceptions();

    input.exceptions(std::ios_base::badbit | std::ios_base::failbit);

    // read in coordinate

    input >> input_character;
    if (input_character != '(') throw std::invalid_argument(error_msg);

    int i = 0;

    while (i < signed(coordinate.size()) - 1)
    {
        input >> coordinate[i++];

        input >> input_character;
        if (input_character != ',') throw std::invalid_argument(error_msg);
    }

    if (coordinate.size() != 0) // allow empty coordinates like '( )'
    {
        input >> coordinate[i];
    }

    input >> input_character;
    if (input_character != ')') throw std::invalid_argument(error_msg);

    // restore old io_state

    input.exceptions(old_iostate);

    return input;
}


template <class T>
std::ostream & operator << (std::ostream & output, const Coordinate<T> & coordinate)
{
    // save stream state

    std::ios_base::fmtflags saved_flags = output.flags();
    // std::streamsize saved_precision = output.precision();

    // throw an exception if the output operation fails for any reason

    output.exceptions(std::ios_base::badbit | std::ios_base::failbit | std::ios_base::eofbit);

    // put out formatted coordinate

    output.clear();

    output // << std::fixed // or std::scientific
           // << std::setprecision(4)
           // << std::showpos
           // << std::setfill('0')
           // << std::internal
           << "(";

    int i = 0;

    while (i < signed(coordinate.size()) - 1)
    {
        output // << std::setw(10)
               << coordinate[i++]
               << ", ";
    }

    if (coordinate.size() != 0) // allow empty coordinates like '( )'
    {
        output // << std::setw(10)
               << coordinate[i];
    }

    output <<  ")";

    // restore old stream state

    // output.precision(saved_precision);
    output.flags(saved_flags);

    return output;
}


template <class T>
TRTK::Coordinate<T> sqrt(const TRTK::Coordinate<T> & coordinate)
{
    TRTK::Coordinate<T> result;
    result.resize(coordinate.size());

    for (unsigned i = 0; i < coordinate.size(); ++i)
    {
        using std::sqrt;
        result[i] = ::sqrt(coordinate[i]);
    }

    return result;
}


} // end of namespace TRTK

#endif // COORDINATE_HPP_1023074891