Warning, /include/Geant4/tools/lina/geom3 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_geom3
0005 #define tools_geom3
0006
0007 #include "line"
0008 #include "plane"
0009 #include "vec3"
0010
0011 namespace tools {
0012
0013 template <class T>
0014 class cubic {
0015 public:
0016 cubic(const vec3<T>& a_p0,const vec3<T>& a_v0,
0017 const vec3<T>& a_p1,const vec3<T>& a_v1) {
0018 // Construct a cubic given 2 points and their tangents.
0019 initialize(a_p0.x(),a_p0.y(),a_p0.z(),
0020 a_v0.x(),a_v0.y(),a_v0.z(),
0021 a_p1.x(),a_p1.y(),a_p1.z(),
0022 a_v1.x(),a_v1.y(),a_v1.z());
0023 }
0024 cubic(const T& a_p0_x,const T& a_p0_y,const T& a_p0_z,
0025 const T& a_v0_x,const T& a_v0_y,const T& a_v0_z,
0026 const T& a_p1_x,const T& a_p1_y,const T& a_p1_z,
0027 const T& a_v1_x,const T& a_v1_y,const T& a_v1_z){
0028 initialize(a_p0_x,a_p0_y,a_p0_z,
0029 a_v0_x,a_v0_y,a_v0_z,
0030 a_p1_x,a_p1_y,a_p1_z,
0031 a_v1_x,a_v1_y,a_v1_z);
0032 }
0033
0034 virtual ~cubic() {}
0035 public:
0036 cubic(const cubic& a_from)
0037 :m_a(a_from.m_a)
0038 ,m_b(a_from.m_b)
0039 ,m_c(a_from.m_c)
0040 ,m_d(a_from.m_d)
0041 {}
0042 cubic& operator=(const cubic& a_from) {
0043 m_a = a_from.m_a;
0044 m_b = a_from.m_b;
0045 m_c = a_from.m_c;
0046 m_d = a_from.m_d;
0047 return *this;
0048 }
0049 public:
0050 void get_point(unsigned int a_index,unsigned int a_num,vec3<T>& a_p){
0051 //a_index = 0 is a_p0
0052 //a_index = a_num-1 is a_p1
0053 T _s = T(a_index)/T(a_num-1);
0054 T s2 = _s*_s;
0055 T s3 = s2*_s;
0056 a_p = m_a*s3 + m_b*s2 + m_c*_s + m_d;
0057 }
0058 void get_point(unsigned int a_index,unsigned int a_num,T& a_x,T& a_y,T& a_z){
0059 //a_index = 0 is a_p0
0060 //a_index = a_num-1 is a_p1
0061 T s = T(a_index)/T(a_num-1);
0062 T s2 = s*s;
0063 T s3 = s2*s;
0064 a_x = m_a.x()*s3 + m_b.x()*s2 + m_c.x()*s + m_d.x();
0065 a_y = m_a.y()*s3 + m_b.y()*s2 + m_c.y()*s + m_d.y();
0066 a_z = m_a.z()*s3 + m_b.z()*s2 + m_c.z()*s + m_d.z();
0067 }
0068 protected:
0069 void initialize(const T& a_p0_x,const T& a_p0_y,const T& a_p0_z,
0070 const T& a_v0_x,const T& a_v0_y,const T& a_v0_z,
0071 const T& a_p1_x,const T& a_p1_y,const T& a_p1_z,
0072 const T& a_v1_x,const T& a_v1_y,const T& a_v1_z){
0073 // Construct a cubic given 2 points and their tangents.
0074
0075 // f(s) = a s**3 + b s**2 + c s + d
0076 // f'(s) = 3 a s**2 + 2 b s + c
0077
0078 // f(0) = d = p0
0079 // f'(0) = c = v0
0080
0081 // f(1) = a + b + v0 + p0 = p1
0082 // f'(1) = 3 a + 2 b + v0 = v1
0083
0084 // f(1) = a + b = p1 - v0 - p0
0085 // f'(1) = 3 a + 2 b = v1 - v0
0086
0087 // b = 3(p1-v0-p0)-(v1-v0)
0088 // a = p1-v0-p0 - b = p1-v0-p0-3(p1-v0-p0)+(v1-v0)
0089 // a = -2p1 + v0 + 2p0 + v1
0090
0091 T a_x = -2*a_p1_x + a_v0_x + 2*a_p0_x + a_v1_x;
0092 T a_y = -2*a_p1_y + a_v0_y + 2*a_p0_y + a_v1_y;
0093 T a_z = -2*a_p1_z + a_v0_z + 2*a_p0_z + a_v1_z;
0094 m_a.set_value(a_x,a_y,a_z);
0095
0096 T b_x = 3*(a_p1_x - a_v0_x - a_p0_x) - (a_v1_x - a_v0_x);
0097 T b_y = 3*(a_p1_y - a_v0_y - a_p0_y) - (a_v1_y - a_v0_y);
0098 T b_z = 3*(a_p1_z - a_v0_z - a_p0_z) - (a_v1_z - a_v0_z);
0099 m_b.set_value(b_x,b_y,b_z);
0100
0101 m_c.set_value(a_v0_x,a_v0_y,a_v0_z);
0102 m_d.set_value(a_p0_x,a_p0_y,a_p0_z);
0103
0104 }
0105
0106 protected:
0107 vec3<T> m_a;
0108 vec3<T> m_b;
0109 vec3<T> m_c;
0110 vec3<T> m_d;
0111 };
0112
0113 }
0114
0115 #include <vector>
0116
0117 namespace tools {
0118
0119 template <class VEC3>
0120 class clip {
0121 protected:
0122 typedef typename VEC3::elem_t T;
0123 public:
0124 clip():m_cur(0){}
0125 virtual ~clip() {}
0126 private:
0127 clip(const clip&):m_cur(0){}
0128 clip& operator=(const clip&){return *this;}
0129 public:
0130 void reset() {
0131 m_data[0].clear();
0132 m_data[1].clear();
0133 m_cur = 0;
0134 }
0135
0136 void add(const VEC3& a_point) {m_data[m_cur].push_back(a_point);}
0137
0138 void execute(const plane<VEC3>& plane) {
0139 //Clip polygon against plane. This might change the number of
0140 //vertices in the polygon.
0141
0142 size_t n = m_data[m_cur].size();
0143 if(!n) return;
0144
0145 // create a loop :
0146 VEC3 dummy = m_data[m_cur][0];
0147 m_data[m_cur].push_back(dummy);
0148
0149 const VEC3& planeN = plane.normal();
0150
0151 for(size_t i = 0; i < n; i++) {
0152 VEC3 v0 = m_data[m_cur][i];
0153 VEC3 v1 = m_data[m_cur][i+1];
0154
0155 T d0 = plane.distance(v0);
0156 T d1 = plane.distance(v1);
0157
0158 if (d0 >= 0.0f && d1 < 0.0f) { // exit plane
0159 VEC3 dir = v1-v0;
0160 // we know that v0 != v1 since we got here
0161 dir.normalize();
0162 T dot = dir.dot(planeN);
0163 VEC3 newvertex = v0 - dir * (d0/dot);
0164 out_point(newvertex);
0165 } else if (d0 < 0.0f && d1 >= 0.0f) { // enter plane
0166 VEC3 dir = v1-v0;
0167 // we know that v0 != v1 since we got here
0168 dir.normalize();
0169 T dot = dir.dot(planeN);
0170 VEC3 newvertex = v0 - dir * (d0/dot);
0171 out_point(newvertex);
0172 out_point(v1);
0173 } else if (d0 >= 0.0f && d1 >= 0.0f) { // in plane
0174 out_point(v1);
0175 }
0176 }
0177 m_data[m_cur].clear();
0178 m_cur ^= 1;
0179 }
0180
0181 const std::vector<VEC3>& result() const {return m_data[m_cur];}
0182
0183 protected:
0184 void out_point(const VEC3& a_p) {m_data[m_cur ^ 1].push_back(a_p);}
0185
0186 protected:
0187 std::vector<VEC3> m_data[2];
0188 unsigned int m_cur;
0189 };
0190
0191 }
0192
0193 #endif
