00001
00014 #ifndef _VECTOR_H
00015 #define _VECTOR_H
00016
00017
00018 #include <iostream>
00019 #include <iomanip>
00020 #include <vector>
00021 #include <numeric>
00022 #include <set>
00023 #include <xmath.H>
00024 #include <error.H>
00025 #include <ListAssignment.H>
00026
00027 using namespace std;
00028
00029 namespace xchen
00030 {
00031
00033 template<class T> class BaseVector
00034 {
00035 public:
00036 virtual size_t size() const { return 0; };
00037 };
00038
00040 template <class T, int Sz>
00041 class Vector : public BaseVector<T>
00042 {
00043 T D[Sz];
00044 typedef const T& R;
00045
00046 public:
00048 typedef T value_type;
00049 typedef T& reference;
00050 typedef const T& const_reference;
00051 typedef T* iterator;
00052 typedef const T* const_iterator;
00053 typedef std::reverse_iterator<iterator> reverse_iterator;
00054 typedef std::reverse_iterator<const_iterator> const_reverse_iterator;
00055 iterator begin() { return D; }
00056 iterator end() { return D+Sz; }
00057 const_iterator begin() const { return D; }
00058 const_iterator end() const { return D+Sz; }
00059 reverse_iterator rbegin() { return reverse_iterator( end() ); }
00060 const_reverse_iterator rbegin() const { return const_reverse_iterator( end() ); }
00061 reverse_iterator rend() { return reverse_iterator( begin() ); }
00062 const_reverse_iterator rend() const { return const_reverse_iterator( begin() ); }
00063 T front() { return D[0]; }
00064 const_reference front() const { return D[0]; }
00065 T back() { return D[Sz-1]; }
00066 const_reference back() const { return D[Sz-1]; }
00067 static bool empty() { return false; }
00068 static size_t max_size() { return Sz; }
00069
00070 T& operator[](size_t idx) { return D[idx]; }
00071 const T& operator[](size_t idx) const { return D[idx]; }
00072 bool operator<(Vector const& rhs) const;
00073 bool operator==(const Vector& rhs) const { return equal(D, D+Sz, rhs.D); }
00075
00076 public:
00078 Vector(const_iterator itr) { copy(itr, itr+Sz, D); }
00079
00081 template <class T2, int Sz2>
00082 Vector(const Vector<T2, Sz2>& rhs) { init_from(rhs); }
00083
00085 template <class T2>
00086 Vector(const BaseVector<T2>& rhs);
00087
00088
00089 Vector() { }
00090
00092 Vector(R e0) { ERR(1>Sz,Args);D[0] = e0; for(int i=1; i<Sz; i++) D[i] = T(.0); }
00093 Vector(R e0,R e1) { ERR(2>Sz,Args);D[0] = e0,D[1] = e1; for(int i=2; i<Sz; i++) D[i] = T(.0); }
00094 Vector(R e0,R e1,R e2) { ERR(3>Sz,Args);D[0] = e0,D[1] = e1,D[2] = e2; for(int i=3; i<Sz; i++) D[i] = T(.0); }
00095 Vector(R e0,R e1,R e2,R e3) { ERR(4>Sz,Args);D[0] = e0,D[1] = e1,D[2] = e2,D[3] = e3; for(int i=4; i<Sz; i++) D[i] = T(.0); }
00096 Vector(R e0,R e1,R e2,R e3,R e4) { ERR(5>Sz,Args);D[0]=e0,D[1]=e1,D[2]=e2,D[3]=e3,D[4]=e4; for(int i=5; i<Sz; i++) D[i]=T(.0); }
00097 Vector(R e0,R e1,R e2,R e3,R e4,R e5) { ERR(6>Sz,Args);D[0]=e0,D[1]=e1,D[2]=e2,D[3]=e3,D[4]=e4,D[5]=e5; for(int i=6; i<Sz; i++) D[i]=T(.0); }
00099
00100 ListAssignmentCheckSzCt<T, T*, Sz-1> operator=(T x)
00101 {
00102 D[0] = x;
00103 return ListAssignmentCheckSzCt<T, T*, Sz-1>(D + 1);
00104 }
00105
00106 size_t size() const { return Sz; }
00107
00108 public:
00110 Vector operator- ();
00111 Vector& operator+= (const Vector &r);
00112 Vector& operator-= (const Vector &r);
00113 Vector& operator*= (double s) { transform(begin(),end(),begin(),scale_fun(s)); return *this; }
00114 Vector& operator/= (double s) { return (*this) *= 1.0/s; }
00115 Vector operator+ (const Vector &r) const { Vector tmp(*this); tmp += r; return tmp; }
00116 Vector operator- (const Vector &r) const { Vector tmp(*this); tmp -= r; return tmp; }
00117 Vector operator* (double s) const { Vector tmp(*this); tmp *= s; return tmp; }
00118 Vector operator/ (double s) const { Vector tmp(*this); tmp /= s; return tmp; }
00119
00120 Vector& Sxw(T s) { assert(Sz==4);D[0]*=s;D[3]*=s; return *this; }
00121 Vector& Syw(T s) { assert(Sz==4);D[1]*=s;D[3]*=s; return *this; }
00122 Vector& Szw(T s) { assert(Sz==4);D[2]*=s;D[3]*=s; return *this; }
00123
00124 T* GetData() { return D; }
00125 const T* GetData() const { return D; }
00126
00127 T operator* (const Vector &r) const { return inner_product(begin(),end(),r.begin(),T(.0)); }
00128 Vector operator^(const Vector &r) const;
00129 Vector& operator^=(const Vector &r);
00130
00131 double GetNorm() const { return sqrt( double( inner_product(begin(),end(),begin(),0.0) ) ); }
00132 T GetSquareNorm() const { return inner_product(begin(),end(),begin(),.0); }
00133 void Normalize() { operator/= (GetNorm()); }
00134 Vector GetNormalizedVector() const { return (*this) / GetNorm(); }
00135
00136 Vector GetProject2Sphere() const { return GetNormalizedVector(); }
00137 Vector GetProject2AffinePlane() const { return (*this) / D[Sz-1]; }
00138 Vector GetProjection(const Vector &r) const { return ( r/r.GetSquareNorm() * (*this * r) ); }
00139 Vector<T,Sz-1> GetProjection() const { assert(Sz>1); return Vector<T,Sz-1>(*this); }
00140 Vector<T,Sz+1> GetPromotion() const { return Vector<T,Sz+1>(*this); }
00141
00142 double Distance(const Vector& v) const { Vector tmp(*this); tmp -= v; return tmp.GetNorm(); }
00143 double Angle(const Vector &r) const;
00145
00146
00147
00149 T LinearCombine(const Vector<double,Sz> lambda) const;
00150 T AffineCombine(const Vector<double,Sz> lambda) const;
00151 T VectorCombine(const Vector<double,Sz> lambda) const;
00153
00154
00155 public:
00156 friend ostream &operator<<<T,Sz>(ostream &os,const Vector &v);
00157 friend istream &operator>><T,Sz>(istream &is,Vector &v);
00158
00159 private:
00160 template <class T2, int Sz2>
00161 void init_from(const Vector<T2, Sz2>& rhs) { bzero(D, sizeof(T)*Sz); copy(rhs.begin(), rhs.size()<Sz ? rhs.end() : rhs.begin()+Sz, D); }
00162
00163 struct less_mag : public binary_function<const T&,const T&,bool> {
00164 bool operator()(const T& x,const T& y) { return fabs(x) < fabs(y); }
00165 };
00166
00167 struct scale_fun : public unary_function<const T&,T>
00168 {
00169 scale_fun(double& s) : scale(s) { }
00170 double scale;
00171 T operator()(const T& x) { return x*scale; }
00172 };
00173 };
00174
00175
00176
00177
00178
00179 template<class T, int Sz>
00180 template<class T2>
00181 Vector<T, Sz> :: Vector<T, Sz>(const BaseVector<T2>& rhs)
00182 {
00183 switch(rhs.size())
00184 {
00185 case 1: init_from( (Vector<T2, 1>&)(rhs) ); break;
00186 case 2: init_from( (Vector<T2, 2>&)(rhs) ); break;
00187 case 3: init_from( (Vector<T2, 3>&)(rhs) ); break;
00188 case 4: init_from( (Vector<T2, 4>&)(rhs) ); break;
00189 case 5: init_from( (Vector<T2, 5>&)(rhs) ); break;
00190 case 6: init_from( (Vector<T2, 6>&)(rhs) ); break;
00191 default: cerr << "Vector size support up to only 6, but given " << rhs.size() << "\n"; exit(0);
00192 }
00193 }
00194
00195
00196 template <class T,int Sz>
00197 Vector<T,Sz> operator*(const T& t,const Vector<T,Sz> &v)
00198 {
00199 return v * t;
00200 }
00201
00202
00203
00204 template <class T,int Sz>
00205 istream& operator>>(istream &is,Vector<T,Sz> &v)
00206 {
00207 typename Vector<T,Sz> :: iterator iter = v.begin();
00208 while( iter != v.end() )
00209 is >> (*iter++);
00210 return is;
00211 }
00212
00213
00214 template <class T,int Sz>
00215 ostream &operator<<(ostream &os,const Vector<T,Sz> &v)
00216 {
00217 os << "[";
00218 for(typename Vector<T,Sz>::const_iterator itr = v.begin(); itr != v.end()-1; ++itr)
00219 {
00220 os << setw(10) << setprecision(4) << *itr << ", ";
00221 }
00222 os << setw(10) << setprecision(4) << v.back() << "]";
00223 return os;
00224 }
00225
00226
00227 template <class T,int Sz>
00228 Vector<T,Sz> Vector<T,Sz> :: operator- ()
00229 {
00230 Vector tmp(*this);
00231 transform(tmp.begin(),tmp.end(),tmp.begin(),negate<T>());
00232 return tmp;
00233 }
00234
00235
00236 template <class T,int Sz>
00237 Vector<T,Sz>& Vector<T,Sz> :: operator+= (const Vector &r)
00238 {
00239 transform(begin(),end(),r.begin(),begin(),plus<T>());
00240 return *this;
00241 }
00242
00243
00244 template <class T,int Sz>
00245 Vector<T,Sz>& Vector<T,Sz> :: operator-= (const Vector &r)
00246 {
00247 transform(begin(),end(),r.begin(),begin(),minus<T>());
00248 return *this;
00249 }
00250
00251
00252 template <class T,int Sz>
00253 double Vector<T,Sz> :: Angle(const Vector &r) const
00254 {
00255 double t = (*this * r) / (GetNorm()*r.GetNorm());
00256 if(approx_eq(t,1.0))
00257 return 0;
00258 if(approx_eq(t,-1.0))
00259 return pi;
00260 return acos( t ) * 180/pi;
00261 }
00262
00263
00264 template <class T,int Sz>
00265 Vector<T,Sz> Vector<T,Sz> :: operator^(const Vector &r) const
00266 {
00267 return Vector(((*this)[1]*r[2] - (*this)[2]*r[1]),
00268 ((*this)[2]*r[0] - (*this)[0]*r[2]),
00269 ((*this)[0]*r[1] - (*this)[1]*r[0]));
00270 }
00271
00272
00273 template <class T,int Sz>
00274 Vector<T,Sz>& Vector<T,Sz> :: operator^=(const Vector &r)
00275 {
00276 #if Sz == 3
00277 T x=(*this)[0],y=(*this)[1],z=(*this)[2];
00278 (*this)[0]=y*r[2] - z*r[1];
00279 (*this)[1]=z*r[0] - x*r[2];
00280 (*this)[2]=x*r[1] - y*r[0];
00281 return *this;
00282 #else
00283 cerr << "cross product is only defined for Vector<T,3>.\n";
00284 exit(0);
00285 #endif
00286 }
00287
00288 template <class T>
00289 inline Vector<T,3> operator^(Vector<T,3>const &v1, Vector<T,3>const &v2)
00290 {
00291 Vector<T,3> ret(v1);
00292 ret ^= v2;
00293 return ret;
00294 }
00295
00296 template <class T,int Sz>
00297 T Vector<T,Sz> :: LinearCombine(const Vector<double,Sz> lambda) const
00298 {
00299 const_iterator itr1 = begin();
00300 typename Vector<double,Sz>::const_iterator itr2 = lambda.begin();
00301 T ret(0.0);
00302 for(; itr1 != end(); itr1++,itr2++)
00303 ret += (*itr1) * (*itr2);
00304 return ret;
00305 }
00306
00307 template <class T,int Sz>
00308 T Vector<T,Sz> :: AffineCombine(const Vector<double,Sz> lambda) const
00309 {
00310 double chksum;
00311 assert((chksum = 0.0,1));
00312
00313 iterator itr1 = begin();
00314 typename Vector<double,Sz>::const_iterator itr2 = lambda.begin();
00315 T ret(0.0);
00316 for(; itr1 != end(); itr1++,itr2++)
00317 {
00318 ret += (*itr1) * (*itr2);
00319 assert((chksum += *itr2,1));
00320 }
00321 assert(approx_eq(chksum,1.0));
00322 return ret;
00323 }
00324
00325
00326 template <class T,int Sz>
00327 T Vector<T,Sz> :: VectorCombine(const Vector<double,Sz> lambda) const
00328 {
00329 double chksum;
00330 assert((chksum = 0.0,1));
00331
00332 iterator itr1 = begin();
00333 typename Vector<double,Sz>::const_iterator itr2 = lambda.begin();
00334 T ret = 0;
00335 for(; itr1 != end(); itr1++,itr2++)
00336 {
00337 ret += (*itr1) * (*itr2);
00338 assert((chksum += *itr2,1));
00339 }
00340 assert(approx_eq(chksum,0.0));
00341 return ret;
00342 }
00343
00344
00345
00346 template <class T,int Sz>
00347 inline bool Vector<T,Sz> :: operator<(Vector<T,Sz> const& rhs) const
00348 {
00349 for(int i=Sz-1; i>=0; i--)
00350 {
00351 if(D[i]<rhs[i])
00352 return true;
00353 else if(D[i]>rhs[i])
00354 return false;
00355 }
00356 return false;
00357 }
00358
00359
00360
00361
00362 #define NormalizedVector Vector
00363
00364 template<class T, int sz>
00365 inline Vector<T,sz> interpolate(const NormalizedVector<T,sz>& t, const Vector<T,sz>& p1, const Vector<T,sz>& p2)
00366 {
00367 Vector<T,sz> p;
00368 for(int i=0; i<sz; i++)
00369 p[i] = interpolate( t[i], p1[i], p2[i] );
00370 return p;
00371 }
00372 template<class T, int sz>
00373 inline NormalizedVector<T,sz> ratio(const Vector<T,sz>& p, const Vector<T,sz>& p1, const Vector<T,sz>& p2)
00374 {
00375 NormalizedVector<T,sz> r;
00376 for(int i=0; i<sz; i++)
00377 r[i] = ratio( p[i], p1[i], p2[i] );
00378 return r;
00379 }
00380
00381
00382
00383
00384 #define Point Vector
00385
00386 typedef BaseVector<double> dVector;
00387 typedef BaseVector<float> fVector;
00388 typedef BaseVector<int> iVector;
00389 typedef BaseVector<char> cVector;
00390
00391 typedef BaseVector<double> dPoint;
00392 typedef BaseVector<float> fPoint;
00393 typedef BaseVector<int> iPoint;
00394 typedef BaseVector<char> cPoint;
00395
00396 typedef Vector<double,5> dVector5D;
00397 typedef Vector<double,4> dVector4D;
00398 typedef Vector<double,3> dVector3D;
00399 typedef Vector<double,2> dVector2D;
00400 typedef Vector<double,1> dVector1D;
00401
00402
00403 typedef Vector<float,5> fVector5D;
00404 typedef Vector<float,4> fVector4D;
00405 typedef Vector<float,3> fVector3D;
00406 typedef Vector<float,2> fVector2D;
00407 typedef Vector<float,1> fVector1D;
00408
00409
00410 typedef Vector<int,5> iVector5D;
00411 typedef Vector<int,4> iVector4D;
00412 typedef Vector<int,3> iVector3D;
00413 typedef Vector<int,2> iVector2D;
00414 typedef Vector<int,1> iVector1D;
00415
00416
00417 typedef Vector<char,5> cVector5D;
00418 typedef Vector<char,4> cVector4D;
00419 typedef Vector<char,3> cVector3D;
00420 typedef Vector<char,2> cVector2D;
00421 typedef Vector<char,1> cVector1D;
00422
00423
00424 typedef iVector5D iPoint5D;
00425 typedef iVector4D iPoint4D;
00426 typedef iVector3D iPoint3D;
00427 typedef iVector2D iPoint2D;
00428 typedef iVector1D iPoint1D;
00429
00430
00431 typedef fVector5D fPoint5D;
00432 typedef fVector4D fPoint4D;
00433 typedef fVector3D fPoint3D;
00434 typedef fVector2D fPoint2D;
00435 typedef fVector1D fPoint1D;
00436
00437
00438 typedef dVector5D dPoint5D;
00439 typedef dVector4D dPoint4D;
00440 typedef dVector3D dPoint3D;
00441 typedef dVector2D dPoint2D;
00442 typedef dVector1D dPoint1D;
00443
00444 typedef cVector5D cPoint5D;
00445 typedef cVector4D cPoint4D;
00446 typedef cVector3D cPoint3D;
00447 typedef cVector2D cPoint2D;
00448 typedef cVector1D cPoint1D;
00449
00450
00451 typedef vector<dVector2D> dVector2Dvector;
00452 typedef vector<dVector3D> dVector3Dvector;
00453 typedef vector<iVector2D> iVector2Dvector;
00454 typedef vector<iVector3D> iVector3Dvector;
00455 typedef vector<fVector2D> fVector2Dvector;
00456 typedef vector<fVector3D> fVector3Dvector;
00457
00458 typedef set<dVector2D, less<dVector2D> > dVector2Dset;
00459 typedef set<dVector3D, less<dVector3D> > dVector3Dset;
00460 typedef set<iVector2D, less<iVector2D> > iVector2Dset;
00461 typedef set<iVector3D, less<iVector3D> > iVector3Dset;
00462 typedef set<fVector2D, less<fVector2D> > fVector2Dset;
00463 typedef set<fVector3D, less<fVector3D> > fVector3Dset;
00464
00465 typedef dVector2Dset dPoint2Dset;
00466 typedef dVector3Dset dPoint3Dset;
00467 typedef iVector2Dset iPoint2Dset;
00468 typedef iVector3Dset iPoint3Dset;
00469 typedef fVector2Dset fPoint2Dset;
00470 typedef fVector3Dset fPoint3Dset;
00471
00472 typedef dVector2Dvector dPoint2Dvector;
00473 typedef dVector3Dvector dPoint3Dvector;
00474 typedef iVector2Dvector iPoint2Dvector;
00475 typedef iVector3Dvector iPoint3Dvector;
00476 typedef fVector2Dvector fPoint2Dvector;
00477 typedef fVector3Dvector fPoint3Dvector;
00478
00479
00480 typedef vector<dVector2D> dVector2Dvector;
00481 typedef vector<dVector3D> dVector3Dvector;
00482 typedef vector<iVector2D> iVector2Dvector;
00483 typedef vector<iVector3D> iVector3Dvector;
00484 typedef vector<fVector2D> fVector2Dvector;
00485 typedef vector<fVector3D> fVector3Dvector;
00486
00487 typedef set<dVector2D, less<dVector2D> > dVector2Dset;
00488 typedef set<dVector3D, less<dVector3D> > dVector3Dset;
00489 typedef set<iVector2D, less<iVector2D> > iVector2Dset;
00490 typedef set<iVector3D, less<iVector3D> > iVector3Dset;
00491 typedef set<fVector2D, less<fVector2D> > fVector2Dset;
00492 typedef set<fVector3D, less<fVector3D> > fVector3Dset;
00493
00494 typedef dVector2Dset dPoint2Dset;
00495 typedef dVector3Dset dPoint3Dset;
00496 typedef iVector2Dset iPoint2Dset;
00497 typedef iVector3Dset iPoint3Dset;
00498 typedef fVector2Dset fPoint2Dset;
00499 typedef fVector3Dset fPoint3Dset;
00500
00501 typedef dVector2Dvector dPoint2Dvector;
00502 typedef dVector3Dvector dPoint3Dvector;
00503 typedef iVector2Dvector iPoint2Dvector;
00504 typedef iVector3Dvector iPoint3Dvector;
00505 typedef fVector2Dvector fPoint2Dvector;
00506 typedef fVector3Dvector fPoint3Dvector;
00507
00508 typedef dPoint5D E5;
00509 typedef dPoint4D E4;
00510 typedef dPoint3D E3;
00511 typedef dPoint2D E2;
00512 typedef dPoint1D E1;
00513
00514 typedef fPoint5D fE5;
00515 typedef fPoint4D fE4;
00516 typedef fPoint3D fE3;
00517 typedef fPoint2D fE2;
00518 typedef fPoint1D fE1;
00519
00520 typedef E5 P4;
00521 typedef E4 P3;
00522 typedef E3 P2;
00523 typedef E2 P1;
00524
00525 typedef Vector<int,5> iE5;
00526 typedef Vector<int,4> iE4;
00527 typedef Vector<int,3> iE3;
00528 typedef Vector<int,2> iE2;
00529 }
00530
00531
00532 #endif
00533
00534
