utils.h

Go to the documentation of this file.

/* Author(s):     Francois Mauger <mauger@lpccaen.in2p3.fr>
 * Creation date: 2008-05-23
 * Last modified: 2015-02-18
 *
 * License:
 *
 * Description:
 *
 *  Utilities.
 *
 */

#ifndef GEOMTOOLS_UTILS_H
#define GEOMTOOLS_UTILS_H 1

// Standard library:
#include <string>
#include <iostream>
#include <list>
#include <vector>

// Third party:
// - Bayeux/datatools:
#include <datatools/bit_mask.h>
#include <datatools/reflection_interface.h>

// This project:
#include <geomtools/geomtools_config.h>
#include <geomtools/clhep.h>

namespace geomtools {

  bool position_is_in(double position_,
                      double start_position_, double delta_position_,
                      double tolerance_ = 0.0,
                      bool bounds_excluded_ = false);

  bool angle_is_in(double angle_,
                   double start_angle_, double delta_angle_,
                   double angular_tolerance_ = 0.0,
                   bool bounds_excluded_ = false);

  typedef std::list<vector_2d> basic_polyline_2d;

  typedef std::pair<vector_3d, vector_3d> basic_segment_3d;
  typedef basic_segment_3d segment_type;

  typedef std::list<vector_3d> basic_polyline_3d;
  typedef std::list<vector_3d> polyline_type;

  typedef std::vector<vector_3d> vertex_col_type;

  typedef std::list<polyline_type> wires_type;

  void print_xy(std::ostream & out_,
                const vector_2d & p_,
                bool endl_ = true);

  std::string to_xy(const vector_2d & p_);

  std::string vector_2d_to_xy(const vector_2d & p_);

  void print_xy_stdout(const vector_2d & p_);

  void print_xy_stderr(const vector_2d & p_);

  void print_xyz(std::ostream & out_,
                 const vector_3d & p_,
                 bool endl_ = true);

  void print(std::ostream & out_, const vector_3d & p_);

  std::string to_xyz(const vector_3d & p_);

  std::string vector_3d_to_xyz(const vector_3d & p_);

  void print_xyz_stdout(const vector_3d & p_);

  void print_xyz_stderr(const vector_3d & p_);

  void set_vector_2d(const std::vector<double> & source_, vector_2d & target_);

  void set_vector_3d(const std::vector<double> & source_, vector_3d & target_);

  bool parse(std::istream & in_, vector_3d & position_, bool nothing_more_ = false);

  bool parse(const std::string & token_, vector_3d & position_);

  void print_xy(std::ostream & out_,
                const basic_polyline_2d & p_,
                bool endl_ = true);

  std::string to_xy(const basic_polyline_2d & p_);

  std::string basic_polyline_2d_to_xy(const basic_polyline_2d & p_);

  void print_xy_stdout(const basic_polyline_2d & p_);

  void print_xy_stderr(const basic_polyline_2d & p_);

  void print_xyz(std::ostream & out_,
                 const basic_polyline_3d & p_,
                 bool endl_ = true);

  std::string to_xyz(const basic_polyline_3d & p_);

  std::string basic_polyline_3d_to_xyz(const basic_polyline_3d & p_);

  void print_xyz_stdout(const basic_polyline_3d & p_);

  void print_xyz_stderr(const basic_polyline_3d & p_);

  enum orientation_type {
    ORIENTATION_INVALID    = -1, 
    ORIENTATION_VERTICAL   =  0, 
    VERTICAL               =  0, 
    ORIENTATION_HORIZONTAL =  1, 
    HORIZONTAL             =  1  
  };

  enum direction_flags_type {
    DIRECTION_FLAGS_NONE   = 0x0,
    DIRECTION_FLAGS_BACK   = datatools::bit_mask::bit00, 
    DIRECTION_FLAGS_FRONT  = datatools::bit_mask::bit01, 
    DIRECTION_FLAGS_LEFT   = datatools::bit_mask::bit02, 
    DIRECTION_FLAGS_RIGHT  = datatools::bit_mask::bit03, 
    DIRECTION_FLAGS_BOTTOM = datatools::bit_mask::bit04, 
    DIRECTION_FLAGS_TOP    = datatools::bit_mask::bit05  
  };

  enum direction_type {
    DIRECTION_INVALID = -1,                
    DIRECTION_BACK    =  0,                
    DIRECTION_XMINUS  =  DIRECTION_BACK,   
    DIRECTION_FRONT   =  1,                
    DIRECTION_XPLUS   =  DIRECTION_FRONT,  
    DIRECTION_LEFT    =  2,                
    DIRECTION_YMINUS  =  DIRECTION_LEFT,   
    DIRECTION_RIGHT   =  3,                
    DIRECTION_YPLUS   =  DIRECTION_RIGHT,  
    DIRECTION_BOTTOM  =  4,                
    DIRECTION_ZMINUS  =  DIRECTION_BOTTOM, 
    DIRECTION_TOP     =  5,                
    DIRECTION_ZPLUS   =  DIRECTION_TOP     
  };

  direction_type get_direction_from_label(const std::string &);

  std::string get_direction_label(direction_type);

  enum vertex_1d_type {
    VERTEX_NONE     = 0x0,
    VERTEX_ALL_BITS = datatools::bit_mask::nbits31,
    VERTEX_ALL      = VERTEX_ALL_BITS
  };

  enum path_1d_type {
    PATH_NONE     = 0x0,
    PATH_ALL_BITS = datatools::bit_mask::nbits31,
    PATH_ALL      = PATH_ALL_BITS
  };

  enum side_2d_type {
    SIDE_NONE     = 0x0,
    EDGE_NONE     = SIDE_NONE,
    SIDE_ALL_BITS = datatools::bit_mask::nbits31,
    EDGE_ALL_BITS = SIDE_ALL_BITS,
    SIDE_ALL      = SIDE_ALL_BITS,
    EDGE_ALL      = SIDE_ALL
  };
  typedef side_2d_type edge_2d_type;

  enum shape_domain_flags_type {
    SHAPE_DOMAIN_NONE                = 0x0,
    SHAPE_DOMAIN_INSIDE              = datatools::bit_mask::bit00,
    SHAPE_DOMAIN_OUTSIDE             = datatools::bit_mask::bit01,
    SHAPE_DOMAIN_ON_SURFACE          = datatools::bit_mask::bit02,
    SHAPE_DOMAIN_INSIDE_DAUGHTER     = datatools::bit_mask::bit03,
    SHAPE_DOMAIN_ON_DAUGHTER_SURFACE = datatools::bit_mask::bit04
  };

  std::string shape_domain_flag_label(shape_domain_flags_type flag_);

  struct constants {
    static const int    NO_INTERCEPT = -1;
    static const double DEFAULT_TOLERANCE;
    static const double DEFAULT_ANGULAR_TOLERANCE;
    static const double ZERO_TOLERANCE;
    static const double USING_PROPER_TOLERANCE;

    static double get_default_tolerance();
    static double get_default_angular_tolerance();
    static double get_zero_tolerance();
    static double get_proper_tolerance();
  };

  class filled_utils
  {
  public:

    enum filled_type {
      FILLED_UNDEFINED    = -1,
      FILLED_NONE         =  0,
      FILLED_BY_ENVELOPE  =  1,
      FILLED_BY_EXTRUSION =  2
    };

    static const std::string & filled_none_label();
    static const std::string & filled_by_envelope_label();
    static const std::string & filled_by_extrusion_label();

    static std::string get_filled_mode_label(filled_utils::filled_type);

    static filled_type get_filled_mode(const std::string & filled_label);

  };

  /* Initialize a rotation matrix for
   * "World coordinates system->Local coordinates system":
   */

  enum euler_angles_type {
    EULER_ANGLES_INVALID = -1, 
    EULER_ANGLES_ZXZ     =  0, 
    EULER_ANGLES_XYX     =  1, 
    EULER_ANGLES_YZY     =  2, 
    EULER_ANGLES_ZYZ     =  3, 
    EULER_ANGLES_XZX     =  4, 
    EULER_ANGLES_YXY     =  5, 
    EULER_ANGLES_XYZ     =  6, 
    EULER_ANGLES_YZX     =  7, 
    EULER_ANGLES_ZXY     =  8, 
    EULER_ANGLES_XZY     =  9, 
    EULER_ANGLES_ZYX     = 10, 
    EULER_ANGLES_YXZ     = 11  
  };

  void create(rotation_3d & rot_,
              double angle0_,
              double angle1_,
              double angle2_,
              euler_angles_type et_);

  void create_zyz(rotation_3d & rot_,
                  double phi_,
                  double theta_,
                  double delta_);

  void create(rotation_3d & rot_,
              double phi_,
              double theta_,
              double delta_);

  void create_zxz(rotation_3d & rot_,
                  double phi_,
                  double theta_,
                  double psi_);

  void create_xyz(rotation_3d & rot_,
                  double phi_,
                  double theta_,
                  double psi_);

  void create_rotation_3d(rotation_3d & rot_,
                          double phi_,
                          double theta_,
                          double delta_);

  void create_rotation_from_zyz_euler_angles(rotation_3d & rot_,
                                             double phi_,
                                             double theta_,
                                             double delta_);

  void create_rotation(rotation_3d & rot_,
                       double phi_,
                       double theta_,
                       double delta_);

  bool is_identity(const rotation_3d & rot_);

  void extract_xyz_euler_angle_from_rotation(const rotation_3d & rot_,
                                             double & a_,
                                             double & b_,
                                             double & c_);

  // http://www.geometrictools.com/Documentation/EulerAngles.pdf
  void extract_zyz_euler_angle_from_rotation(const rotation_3d & rot_,
                                             double & a_,
                                             double & b_,
                                             double & c_);

  // http://www.cgafaq.info/wiki/Euler_angles_from_matrix
  // struct extract_euler_angle
  // {
  //   int i, neg, alt, rev;
  //   extract_euler_angle()  {}
  // };

  enum special_rotation_angle_type {
    ROTATION_ANGLE_INVALID = -1, 
    ROTATION_ANGLE_0   = 0, 
    ROTATION_ANGLE_90  = 1, 
    ROTATION_ANGLE_180 = 2, 
    ROTATION_ANGLE_270 = 3  
  };

  bool check_special_rotation_angle(int);

  double get_special_rotation_angle(int);

  int get_special_rotation_angle_from_label(const std::string & );

  std::string get_special_rotation_angle_label(int);

  enum axis_type {
    AXIS_INVALID = -1, 
    AXIS_X = 0, 
    AXIS_Y = 1, 
    AXIS_Z = 2  
  };

  enum rotation_axis_type {
    ROTATION_AXIS_INVALID = -1, 
    ROTATION_AXIS_X = AXIS_X, 
    ROTATION_AXIS_Y = AXIS_Y, 
    ROTATION_AXIS_Z = AXIS_Z  
  };

  bool check_rotation_axis(int);

  int get_rotation_axis_from_label(const std::string & );

  std::string get_rotation_label(int);

  void create_rotation_from_axis(rotation_3d & rot_,
                                 int axis_,
                                 double angle_);

  void create_rotation(rotation_3d & rot_,
                       int axis_,
                       double angle_);

  void create_rotation(rotation_3d & rot_,
                       int axis_,
                       int special_angle_);

  void create_rotation_from(rotation_3d & rot_,
                            const std::string &);

  void reset(rotation_3d & rot_);

  void reset_rotation_3d(rotation_3d & rot_);

  void tree_dump(const rotation_3d & rot_,
                 std::ostream & out_,
                 const std::string & title_ = "",
                 const std::string & indent_ = "");

  void invalidate(rotation_3d & rot_);

  void rectify(rotation_3d & rot_);

  void invalidate_rotation_3d(rotation_3d & rot_);

  bool is_valid(const rotation_3d & rot_);

  bool is_valid_rotation_3d(const rotation_3d & rot_);

  void create(vector_3d &,
              double x_,
              double y_,
              double z_);

  void create_xyz(vector_3d &,
                  double x_,
                  double y_,
                  double z_);

  void create_polar(vector_3d &,
                    double r_,
                    double theta_,
                    double z_);

  void create_spherical(vector_3d &,
                        double r_,
                        double phi_,
                        double theta_);

  void set(vector_3d & vec_, double x_, double y_, double z_);

  void set_r_theta_phi(vector_3d & vec_, double r_, double theta_, double phi_);

  void set_rho_phi_z(vector_3d & vec_, double rho_, double phi_, double z_);

  void zero(vector_3d & vec_);
  
  const vector_3d & invalid_vector_3d();

  void invalidate(vector_3d & vec_);

  void invalidate_vector_3d(vector_3d & vec_);

  bool is_valid(const vector_3d & vec_);

  bool is_valid_vector_3d(const vector_3d & vec_);

  bool are_near(const vector_3d & vec1_,
                const vector_3d & vec2_,
                double tolerance_ = constants::DEFAULT_TOLERANCE);

  bool are_near_vector_3d(const vector_3d & vec1_,
                          const vector_3d & vec2_,
                          double tolerance_);

  void set(vector_2d & vec_, double x_, double y_);

  void set_r_phi(vector_2d & vec_, double r_, double phi_);

  void zero(vector_2d & vec_);

  void invalidate(vector_2d & vec_);

  void invalidate_vector_2d(vector_2d & vec_);

  const vector_2d & invalid_vector_2d();

  bool is_valid(const vector_2d & vec_);

  bool is_valid_vector_2d(const vector_2d & vec_);

  void vector_2d_to_vector_3d(const vector_2d & v2d_, vector_3d & v3d_);

  void vector_3d_to_vector_2d(const vector_3d & v3d_, vector_2d & v2d_);

  void make_phi_theta(vector_3d & vec_, double phi_, double theta_);

#if __GNUC__ >= 7
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wnoexcept-type"
#endif
  template <class ran_func>
  void randomize_direction(ran_func & ran_, vector_3d & ran_dir_) 
  {
    double phi = 2. * M_PI * ran_();
    double cos_theta = -1 + 2 * ran_();
    double sin_theta = std::sqrt(1. - cos_theta * cos_theta);
    double x = sin_theta * std::cos(phi);
    double y = sin_theta * std::sin(phi);
    double z = cos_theta;
    ran_dir_.set(x, y, z);
  }
#if __GNUC__ >= 7
#pragma GCC diagnostic pop
#endif

  template <class ran_func>
  vector_3d randomize_direction(ran_func & ran_)
  {
    vector_3d dir;
    randomize_direction(ran_, dir);
    return dir;
  }

  template <class ran_func>
  void randomize_orthogonal_direction(ran_func & ran_,
                                      const vector_3d & dir_,
                                      vector_3d & ran_dir_)
  {
    double theta = 2. * M_PI * ran_();
    double dx = cos(theta);
    double dy = sin(theta);
    double dz = 0.0;
    double dir_theta = dir_.theta();
    double dir_phi = dir_.phi();
    rotation_3d dir_rot;
    create_rotation_3d(dir_rot, dir_phi, dir_theta, 0.0);
    rotation_3d dir_inverse_rot;
    dir_inverse_rot = dir_rot.inverse();
    vector_3d v(dx, dy, dz);
    ran_dir_ = v.transform(dir_inverse_rot);
    return;
  }

  template <class ran_func>
  vector_3d randomize_orthogonal_direction(ran_func & ran_,
                                           const vector_3d & ref_dir_)
  {
    vector_3d dir;
    randomize_direction(ran_, ref_dir_, dir);
    return dir;
  }

  void compute_barycenter(const std::vector<vector_3d> & points_,
                          vector_3d & barycenter_);

  vector_3d compute_barycenter(const std::vector<vector_3d> & points_);

  void compute_weighted_barycenter(const std::vector<vector_3d> & points_,
                                   const std::vector<double> & weights_,
                                   vector_3d & weighted_barycenter_);

  vector_3d compute_weighted_barycenter(const std::vector<vector_3d> & points_,
                                        const std::vector<double> & weights_);

  class rotation_wrapper_type : public rotation_3d
  {
  public:
    rotation_wrapper_type(double mxx_, double mxy_, double mxz_,
                          double myx_, double myy_, double myz_,
                          double mzx_, double mzy_, double mzz_)
      : rotation_3d(mxx_, mxy_, mxz_,
                    myx_, myy_, myz_,
                    mzx_, mzy_, mzz_) {}
  };

  typedef rotation_wrapper_type rotation_wrapper_t;

  struct io
  {
    static const std::string & vector_2d_serial_tag();
    static const std::string & vector_3d_serial_tag();
    static const std::string & rotation_3d_serial_tag();
    static const std::string & position_suffix();
    static const std::string & rotation_suffix();
  };

} // end namespace geomtools

namespace boost {

  namespace serialization {

    template<class Archive>
    void save(Archive & a_ar ,
              const geomtools::vector_3d & v_,
              const unsigned int a_version);

    template<class Archive>
    void load(Archive & a_ar ,
              geomtools::vector_3d & v_,
              const unsigned int a_version);

    template<class Archive>
    void serialize(Archive & a_ar,
                   geomtools::vector_3d  & v_,
                   const unsigned int a_version);

  } // namespace serialization

} // namespace boost

namespace boost {

  namespace serialization {

    template<class Archive>
    void save(Archive & a_ar ,
              const geomtools::vector_2d & v_,
              const unsigned int a_version);

    template<class Archive>
    void load(Archive & a_ar ,
              geomtools::vector_2d & v_,
              const unsigned int a_version);

    template<class Archive>
    void serialize(Archive & a_ar,
                   geomtools::vector_2d  & v_,
                   const unsigned int a_version);

  } // namespace serialization

} // namespace boost

namespace boost {

  namespace serialization {

    template<class Archive>
    void save(Archive & a_ar ,
              const geomtools::rotation_3d & r_,
              const unsigned int a_version);

    template<class Archive>
    void load(Archive & a_ar ,
              geomtools::rotation_3d & a_rotation,
              const unsigned int a_version);

    template<class Archive>
    void serialize(Archive & a_ar,
                   geomtools::rotation_3d  & a_rotation,
                   const unsigned int a_version);

  } // namespace serialization

} // namespace boost

// Activate reflection layer for some 'geomtools::XXX' enums:
DR_TYPE_INIT(::geomtools::orientation_type)
DR_TYPE_INIT(::geomtools::direction_type)
DR_TYPE_INIT(::geomtools::axis_type)
DR_TYPE_INIT(::geomtools::rotation_axis_type)

#endif // GEOMTOOLS_UTILS_H

/*
** Local Variables: --
** mode: c++ --
** c-file-style: "gnu" --
** tab-width: 2 --
** End: --
*/