---

math_tools.h

Go to the documentation of this file.

#ifndef __MATH_TOOLS__
#define __MATH_TOOLS__




#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <ctime>

#ifndef WITHOUT_LIMITS
#include <limits>
#endif

#include "msg_stream.h"




namespace Math_tools{
  
  template<typename Vec2,typename Real>
  inline void barycentric_coordinates(const Vec2& A,
                                      const Vec2& B,
                                      const Vec2& C,
                                      const Vec2& M,
                                      Real* const a,
                                      Real* const b,
                                      Real* const c);
  
  inline void init_random() {std::srand(std::time(NULL));}

  template<typename Real>
  inline Real random(const Real min_value_included,
                     const Real max_value_included);

  
  template<typename ValueIterator>
  inline double entropy(ValueIterator begin,
                        ValueIterator end);


  
  template<typename Real>
  inline Real clamp(const Real min_value,
                    const Real max_value,
                    const Real x);

  
  template<typename Array,typename Real>
  inline
  typename Array::value_type
  bilinear_interpolation(const Array& array,
                         const Real x,
                         const Real y);


  inline double degree_to_radian(const double d);
  
  inline double radian_to_degree(const double r);

  template<typename Real>
  inline Real sign(const Real r);


#ifndef WITHOUT_LIMITS    

  template<typename T> inline bool is_quiet_NaN(T x);

  template<typename T> inline bool is_signaling_NaN(T x);

  template<typename T> inline bool is_NaN(T x);
  
#endif

/*
  
  #############################################
  #############################################
  #############################################
  ######                                 ######
  ######   I M P L E M E N T A T I O N   ######
  ######                                 ######
  #############################################
  #############################################
  #############################################
  
*/


  
  
  template<typename Vec2,typename Real>
  void barycentric_coordinates(const Vec2& A,
                               const Vec2& B,
                               const Vec2& C,
                               const Vec2& M,
                               Real* const a,
                               Real* const b,
                               Real* const c){

    const Real Xa = A[0] - C[0];
    const Real Ya = A[1] - C[1];

    const Real Xb = B[0] - C[0];
    const Real Yb = B[1] - C[1];

    const Real Xm = M[0] - C[0];
    const Real Ym = M[1] - C[1];

    const Real det = Xa*Yb - Ya*Xb;

    if (det != 0){
      
      const Real inv_det = 1.0/det;

      *a = inv_det * (Yb*Xm  - Xb*Ym);
      *b = inv_det * (-Ya*Xm + Ym*Xa);
      *c = 1.0 - (*a) - (*b);
    }
    else{

      Message::error
        <<"barycentric_coordinates: colinear points not yet implemented"
        <<Message::done;
      
    }
  }


  

  template<typename Real>
  Real random(const Real min_value_included,
              const Real max_value_included){
    
    const double delta = static_cast<double>(max_value_included - min_value_included);
    
    return min_value_included + static_cast<Real>(delta*rand()/(RAND_MAX));
  }



  
  template<typename ValueIterator>
  double entropy(ValueIterator begin,
                 ValueIterator end){
    
    double sum = 0;
    for(ValueIterator i=begin;i!=end;i++){
      if (*i>=0){
        sum+=*i;
      }
      else{
        Message::error<<"entropy: non-positive value"<<Message::done;
      }
    }

    if (sum==0){
      Message::error<<"entropy: sum==0"<<Message::done;
    }

    double result = 0;
    for(ValueIterator i=begin;i!=end;i++){
      const double p = (*i)/sum;
      result -= (p==0) ? 0 : p*log(p);
    }

    return result;
  }

  

  template<typename Real>
  inline Real clamp(const Real min_value,
                    const Real max_value,
                    const Real x){

    return std::max(std::min(x,max_value),min_value);
  }


  template<typename Array,typename Real>
  inline
  typename Array::value_type
  bilinear_interpolation(const Array& array,
                         const Real x,
                         const Real y){

    typedef unsigned int size_type;
    typedef float        real_type;

    const size_type x_size = array.x_size();
    const size_type y_size = array.y_size();
    
    const size_type x_index  =
      clamp<size_type>(0,x_size-1,static_cast<size_type>(x));
    const size_type xx_index =
      clamp<size_type>(0,x_size-1,x_index+1);
    
    const size_type y_index  =
      clamp<size_type>(0,y_size-1,static_cast<size_type>(y));
    const size_type yy_index =
      clamp<size_type>(0,y_size-1,y_index+1);

    const real_type x_alpha = x - x_index;
    const real_type y_alpha = y - y_index;

    return
      (1.0f-x_alpha) * (1.0f-y_alpha) * array(x_index, y_index) +
      x_alpha        * (1.0f-y_alpha) * array(xx_index,y_index) +
      (1.0f-x_alpha) * y_alpha        * array(x_index, yy_index) +
      x_alpha        * y_alpha        * array(xx_index,yy_index);
      
  }


  inline double degree_to_radian(const double d){
    return d*M_PI/180;
  }
  

  inline double radian_to_degree(const double r){
    return r*180/M_PI;
  }


  template<typename Real>
  inline Real sign(const Real r){
    return (r==0) ? 0 : ((r>0) ? 1 : -1);
  }


  
#ifndef WITHOUT_LIMITS    

  template<typename T> bool is_quiet_NaN(T x){

    const unsigned int size = sizeof(T)/sizeof(int);
    
    T ref = std::numeric_limits<T>::quiet_NaN();
    
    int* index_ref = reinterpret_cast<int*>(&ref);
    int* index_cmp = reinterpret_cast<int*>(&x);

    for(unsigned int i=0;i<size;i++){
      if (index_ref[i]!=index_cmp[i]){
        return false;
      }
    }

    return true;
  }

  
  template<typename T> bool is_signaling_NaN(T x){
    
    const unsigned int size = sizeof(T)/sizeof(int);
    
    T ref = std::numeric_limits<T>::signaling_NaN();
    
    int* index_ref = reinterpret_cast<int*>(&ref);
    int* index_cmp = reinterpret_cast<int*>(&x);

    for(unsigned int i=0;i<size;i++){
      if (index_ref[i]!=index_cmp[i]){
        return false;
      }
    }

    return true;
  }


  template<typename T> bool is_NaN(T x){
    return (is_quiet_NaN(x)||is_signaling_NaN(x));
  }
  
#endif

  
} // end of namespace


#endif

---