Warning, /include/Geant4/tools/lina/vec3 is written in an unsupported language. File is not indexed.
0001 // Copyright (C) 2010, Guy Barrand. All rights reserved.
0002 // See the file tools.license for terms.
0003
0004 #ifndef tools_vec3
0005 #define tools_vec3
0006
0007 #include <cstddef> //size_t
0008
0009 #ifdef TOOLS_MEM
0010 #include "../mem"
0011 #endif
0012
0013 namespace tools {
0014
0015 template <class T>
0016 class vec3 {
0017 #ifdef TOOLS_MEM
0018 static const std::string& s_class() {
0019 static const std::string s_v("tools::vec3");
0020 return s_v;
0021 }
0022 #endif
0023 public:
0024 typedef T elem_t;
0025 unsigned int dimension() const {return 3;}
0026 public:
0027 vec3(){
0028 #ifdef TOOLS_MEM
0029 mem::increment(s_class().c_str());
0030 #endif
0031 m_data[0] = T();
0032 m_data[1] = T();
0033 m_data[2] = T();
0034 }
0035 vec3(const T a_vec[3]) {
0036 #ifdef TOOLS_MEM
0037 mem::increment(s_class().c_str());
0038 #endif
0039 m_data[0] = a_vec[0];
0040 m_data[1] = a_vec[1];
0041 m_data[2] = a_vec[2];
0042 }
0043 vec3(const T& a0,const T& a1,const T& a2
0044 #ifdef TOOLS_MEM
0045 ,bool a_inc = true
0046 #endif
0047 ) {
0048 #ifdef TOOLS_MEM
0049 if(a_inc) mem::increment(s_class().c_str());
0050 #endif
0051 m_data[0] = a0;
0052 m_data[1] = a1;
0053 m_data[2] = a2;
0054 }
0055 virtual ~vec3() {
0056 #ifdef TOOLS_MEM
0057 mem::decrement(s_class().c_str());
0058 #endif
0059 }
0060 public:
0061 vec3(const vec3& a_from){
0062 #ifdef TOOLS_MEM
0063 mem::increment(s_class().c_str());
0064 #endif
0065 m_data[0] = a_from.m_data[0];
0066 m_data[1] = a_from.m_data[1];
0067 m_data[2] = a_from.m_data[2];
0068 }
0069 vec3& operator=(const vec3& a_from) {
0070 m_data[0] = a_from.m_data[0];
0071 m_data[1] = a_from.m_data[1];
0072 m_data[2] = a_from.m_data[2];
0073 return *this;
0074 }
0075 public:
0076 const T& v0() const { return m_data[0];}
0077 const T& v1() const { return m_data[1];}
0078 const T& v2() const { return m_data[2];}
0079
0080 void v0(const T& a_value) { m_data[0] = a_value;}
0081 void v1(const T& a_value) { m_data[1] = a_value;}
0082 void v2(const T& a_value) { m_data[2] = a_value;}
0083
0084 const T& x() const {return m_data[0];}
0085 const T& y() const {return m_data[1];}
0086 const T& z() const {return m_data[2];}
0087
0088 T& x() {return m_data[0];}
0089 T& y() {return m_data[1];}
0090 T& z() {return m_data[2];}
0091
0092 void set_value(const T& a0,const T& a1,const T& a2) {
0093 m_data[0] = a0;
0094 m_data[1] = a1;
0095 m_data[2] = a2;
0096 }
0097 void set_value(const T aV[3]) {
0098 m_data[0] = aV[0];
0099 m_data[1] = aV[1];
0100 m_data[2] = aV[2];
0101 }
0102 void value(T& a0,T& a1,T& a2) const {
0103 a0 = m_data[0];
0104 a1 = m_data[1];
0105 a2 = m_data[2];
0106 }
0107
0108 //bool set_value(unsigned int a_index,const T& a_value) {
0109 // if(a_index>=3) return false;
0110 // m_[a_index] = a_value;
0111 // return true;
0112 //}
0113
0114 T length(T(*a_sqrt)(T)) const {
0115 return a_sqrt(m_data[0]*m_data[0]+m_data[1]*m_data[1]+m_data[2]*m_data[2]);
0116 }
0117
0118 T normalize(T(*a_sqrt)(T)) {
0119 T norme = length(a_sqrt);
0120 if(norme==T()) return T();
0121 divide(norme);
0122 return norme;
0123 }
0124
0125 T dot(const vec3& aV) const {
0126 return (m_data[0] * aV.m_data[0] +
0127 m_data[1] * aV.m_data[1] +
0128 m_data[2] * aV.m_data[2]);
0129 }
0130
0131 void cross(const vec3<T>& aV,vec3<T>& a_value) const {
0132 a_value.set_value(m_data[1] * aV.m_data[2] - m_data[2] * aV.m_data[1],
0133 m_data[2] * aV.m_data[0] - m_data[0] * aV.m_data[2],
0134 m_data[0] * aV.m_data[1] - m_data[1] * aV.m_data[0]);
0135 }
0136
0137 bool equal(const vec3& aV) const {
0138 if(m_data[0]!=aV.m_data[0]) return false;
0139 if(m_data[1]!=aV.m_data[1]) return false;
0140 if(m_data[2]!=aV.m_data[2]) return false;
0141 return true;
0142 }
0143
0144 template <class PREC>
0145 bool equal_prec(const vec3& a_v,PREC a_prec,PREC(*a_fabs)(const T&)) const {
0146 if(&a_v==this) return true;
0147 for(unsigned int index=0;index<3;index++) {
0148 T diff = m_data[index]-a_v.m_data[index];
0149 if(a_fabs(diff)>=a_prec) return false;
0150 }
0151 return true;
0152 }
0153
0154 vec3<T> _cross(const vec3<T>& aV) const {
0155 //not effective.
0156 return vec3<T>(m_data[1] * aV.m_data[2] - m_data[2] * aV.m_data[1],
0157 m_data[2] * aV.m_data[0] - m_data[0] * aV.m_data[2],
0158 m_data[0] * aV.m_data[1] - m_data[1] * aV.m_data[0]);
0159 }
0160
0161 bool divide(const T& a_T) {
0162 if(a_T==T()) return false;
0163 m_data[0] /= a_T;
0164 m_data[1] /= a_T;
0165 m_data[2] /= a_T;
0166 return true;
0167 }
0168
0169 void multiply(const T& a_T) {
0170 m_data[0] *= a_T;
0171 m_data[1] *= a_T;
0172 m_data[2] *= a_T;
0173 }
0174
0175 void add(const vec3& a_v) {
0176 m_data[0] += a_v.m_data[0];
0177 m_data[1] += a_v.m_data[1];
0178 m_data[2] += a_v.m_data[2];
0179 }
0180
0181 void add(const T& a0,const T& a1,const T& a2) {
0182 m_data[0] += a0;
0183 m_data[1] += a1;
0184 m_data[2] += a2;
0185 }
0186
0187 void subtract(const vec3& a_v) {
0188 m_data[0] -= a_v.m_data[0];
0189 m_data[1] -= a_v.m_data[1];
0190 m_data[2] -= a_v.m_data[2];
0191 }
0192
0193 void subtract(const T& a0,const T& a1,const T& a2) {
0194 m_data[0] -= a0;
0195 m_data[1] -= a1;
0196 m_data[2] -= a2;
0197 }
0198
0199 bool cos_angle(const vec3& a_v,T& a_cos,T(*a_sqrt)(T)) const {
0200 //WARNING : if ret false, a_cos is not set.
0201 T this_length = length(a_sqrt);
0202 if(this_length==T()) return false;
0203 T a_v_length = a_v.length(a_sqrt);
0204 if(a_v_length==T()) return false;
0205 a_cos = dot(a_v)/(this_length*a_v_length);
0206 return true;
0207 }
0208
0209 bool theta_phi(T& a_theta,T& a_phi,T(*a_sqrt)(T),T(*a_atan2)(T,T)) const {
0210 //WARNING : if ret false, a_theta, a_phi are not set.
0211 if(length(a_sqrt)==T()) return false;
0212 a_phi = a_atan2(m_data[1],m_data[0]);
0213 T xy = a_sqrt(m_data[0]*m_data[0]+m_data[1]*m_data[1]);
0214 a_theta = a_atan2(xy,m_data[2]);
0215 return true;
0216 }
0217
0218 public: //operators
0219 T& operator[](size_t a_index) {
0220 //WARNING : no check on a_index.
0221 return m_data[a_index];
0222 }
0223 const T& operator[](size_t a_index) const {
0224 //WARNING : no check on a_index.
0225 return m_data[a_index];
0226 }
0227
0228 vec3& operator*=(const T& a_v) {
0229 m_data[0] *= a_v;
0230 m_data[1] *= a_v;
0231 m_data[2] *= a_v;
0232 return *this;
0233 }
0234
0235 vec3 operator+(const vec3& a_v) const {
0236 return vec3(m_data[0]+a_v.m_data[0],
0237 m_data[1]+a_v.m_data[1],
0238 m_data[2]+a_v.m_data[2]);
0239 }
0240
0241 vec3 operator-(const vec3& a_v) const {
0242 return vec3(m_data[0]-a_v.m_data[0],
0243 m_data[1]-a_v.m_data[1],
0244 m_data[2]-a_v.m_data[2]);
0245 }
0246
0247 vec3 operator*(const T& a_v) const {
0248 return vec3(m_data[0]*a_v,
0249 m_data[1]*a_v,
0250 m_data[2]*a_v);
0251 }
0252
0253 vec3 operator/(const T& a_v) const {
0254 if(a_v==T()) return vec3();
0255 return vec3(m_data[0]/a_v,
0256 m_data[1]/a_v,
0257 m_data[2]/a_v);
0258 }
0259
0260 bool operator==(const vec3& a_v) const {return equal(a_v);}
0261 bool operator!=(const vec3& a_v) const {return !operator==(a_v);}
0262
0263 public: //for tools/sg/sf_vec
0264 typedef unsigned int size_type;
0265 size_type size() const {return 3;}
0266 const T* data() const {return m_data;}
0267 size_type data_size() const {return 3;} //for eqT.
0268 public: //for iv2sg
0269 const T* getValue() const {return m_data;}
0270 void setValue(const T& a0,const T& a1,const T& a2) {
0271 m_data[0] = a0;
0272 m_data[1] = a1;
0273 m_data[2] = a2;
0274 }
0275 void getValue(T& a0,T& a1,T& a2) const {
0276 a0 = m_data[0];
0277 a1 = m_data[1];
0278 a2 = m_data[2];
0279 }
0280 void setValue(const vec3& a_v) {
0281 m_data[0] = a_v.m_data[0];
0282 m_data[1] = a_v.m_data[1];
0283 m_data[2] = a_v.m_data[2];
0284 }
0285 void setValue(const T aV[3]) {
0286 m_data[0] = aV[0];
0287 m_data[1] = aV[1];
0288 m_data[2] = aV[2];
0289 }
0290
0291 vec3& setValue(const vec3& a_bary,
0292 const vec3& a_v0,const vec3& a_v1,const vec3& a_v2) {
0293 m_data[0] = a_bary[0]*a_v0[0]+a_bary[1]*a_v1[0]+a_bary[2]*a_v2[0];
0294 m_data[1] = a_bary[0]*a_v0[1]+a_bary[1]*a_v1[1]+a_bary[2]*a_v2[1];
0295 m_data[2] = a_bary[0]*a_v0[2]+a_bary[1]*a_v1[2]+a_bary[2]*a_v2[2];
0296 return *this;
0297 }
0298
0299 public:
0300 #if defined(TOOLS_MEM) && !defined(TOOLS_MEM_ATEXIT)
0301 static const vec3<T>& s_x() {static const vec3<T> s_v(1,0,0,false);return s_v;}
0302 static const vec3<T>& s_y() {static const vec3<T> s_v(0,1,0,false);return s_v;}
0303 static const vec3<T>& s_z() {static const vec3<T> s_v(0,0,1,false);return s_v;}
0304 #else
0305 static const vec3<T>& s_x() {static const vec3<T> s_v(1,0,0);return s_v;}
0306 static const vec3<T>& s_y() {static const vec3<T> s_v(0,1,0);return s_v;}
0307 static const vec3<T>& s_z() {static const vec3<T> s_v(0,0,1);return s_v;}
0308 #endif
0309 protected:
0310 T m_data[3];
0311
0312 private:static void check_instantiation() {vec3<float> v;}
0313 };
0314
0315 //for sf, mf :
0316 template <class T>
0317 inline const T* get_data(const vec3<T>& a_v) {return a_v.data();}
0318
0319 template <class T>
0320 inline void get_normal(const vec3<T>& a_p0,const vec3<T>& a_p1,const vec3<T>& a_p2,vec3<T>& a_nm,
0321 vec3<T>& a_tmp_1,vec3<T>& a_tmp_2,T(*a_sqrt)(T)) {
0322 // Used to optimize sg::bin().
0323 //(a_p1-a_p0).cross(a_p2-a_p1,a_nm);
0324
0325 a_tmp_1 = a_p1;
0326 a_tmp_1.subtract(a_p0);
0327
0328 a_tmp_2 = a_p2;
0329 a_tmp_2.subtract(a_p1);
0330
0331 a_tmp_1.cross(a_tmp_2,a_nm);
0332
0333 a_nm.normalize(a_sqrt);
0334 }
0335
0336 /*
0337 template <class VEC3>
0338 inline void get_normal(const VEC3& a_p0,const VEC3& a_p1,const VEC3& a_p2,VEC3& a_nm) {
0339 VEC3 tmp1,tmp2;
0340 get_normal(a_p0,a_p1,a_p2,a_nm,tmp1,tmp2);
0341 }
0342 */
0343 template <class VEC3>
0344 inline void direction(const VEC3& a_p0,const VEC3& a_p1,const VEC3& a_p2,VEC3& a_value) {
0345 // Orientation is computed by taking (p1 - p0) x (p2 - p0)
0346 VEC3 P = a_p1;
0347 P.subtract(a_p0);
0348 VEC3 P2 = a_p2;
0349 P2.subtract(a_p0);
0350 P.cross(P2,a_value);
0351 }
0352
0353 template <class VEC3>
0354 inline void area(const VEC3& a_p0,const VEC3& a_p1,const VEC3& a_p2,typename VEC3::elem_t& a_value,
0355 VEC3& a_tmp_1,VEC3& a_tmp_2,VEC3& a_tmp_3) {
0356 // area of the triangle (a_p0,a_p1,a_p2)
0357 typedef typename VEC3::elem_t T;
0358
0359 a_tmp_1 = a_p1;
0360 a_tmp_1.subtract(a_p0);
0361
0362 a_tmp_2 = a_p2;
0363 a_tmp_2.subtract(a_p1);
0364
0365 a_tmp_1.cross(a_tmp_2,a_tmp_3);
0366
0367 a_value = a_tmp_3.length()/T(2);
0368 }
0369
0370 template <class VEC3>
0371 inline void area(const VEC3& a_p0,const VEC3& a_p1,const VEC3& a_p2,typename VEC3::elem_t& a_value) {
0372 VEC3 tmp1,tmp2,tmp3;
0373 area(a_p0,a_p1,a_p2,a_value,tmp1,tmp2,tmp3);
0374 }
0375
0376 template <class T>
0377 inline void direction(const T& a_0_x,const T& a_0_y,const T& a_0_z,
0378 const T& a_1_x,const T& a_1_y,const T& a_1_z,
0379 const T& a_2_x,const T& a_2_y,const T& a_2_z,vec3<T>& a_value) {
0380 direction(vec3<T>(a_0_x,a_0_y,a_0_z),
0381 vec3<T>(a_1_x,a_1_y,a_1_z),
0382 vec3<T>(a_2_x,a_2_y,a_2_z),a_value);
0383 }
0384
0385 }
0386
0387 #include <ostream>
0388
0389 namespace tools {
0390
0391 // for sf_vec::dump().
0392 template <class T>
0393 inline std::ostream& operator<<(std::ostream& a_out,const vec3<T>& a_this){
0394 a_out << "x = " << a_this.v0()
0395 << ",y = " << a_this.v1()
0396 << ",z = " << a_this.v2();
0397 return a_out;
0398 }
0399
0400 }
0401
0402 #endif