JunkBox_Lib++ (for Windows): MathLib/Vector.h Source File

00001 
00002 #ifndef  __JBXL_CPP_VECTOR_H_
00003 #define  __JBXL_CPP_VECTOR_H_
00004 
00005 
00014 /*
00015 template <typename T=double> class DllExport Vector      // ３次元ベクトルの定義
00016 template <typename T=double> class DllExport Vector2D    // ２次元ベクトルの定義, パラメトリック曲面用
00017 template <typename T=int> class DllExport PCoordinate   // 多角形の座標の値を格納するクラス
00018 
00019 template <typename T=int> class DllExport RBound        // 境界構造体
00020 template <typename T> inline bool disJunctBounds(RBound<T> b1, RBound<T> b2);
00021 
00022 template <typename T> inline double VectorAngle(Vector<T> a, Vector<T> b)
00023 template <typename T> inline Vector<T> NewellMethod(Vector<T> v1, Vector<T> v2, Vector<T> v3)
00024 template <typename T> inline Vector<T> BSplineInterp4(Vector<T> p0, Vector<T> p1, double t)
00025 */
00026 
00027 
00028 
00029 #include "tools++.h"
00030 #include <math.h>
00031 #include "Tolerance.h"
00032 
00033 
00034 
00035 //
00036 namespace  jbxl {
00037 
00038 
00039 
00040 #define  VECTOR     Vector
00041 #define  UVMAP      UVMap
00042 
00043 
00044 #define  BOUNDARY_BLANK 5       
00045 
00046 
00048 // Vector
00049 //
00050 
00057 template <typename T=double> class DllExport Vector
00058 {
00059 public:
00060     T x;
00061     T y;
00062     T z;
00063 
00064     double n;       
00065     double c;       
00066     int    d;       
00067 
00068 public:
00069     Vector(T X=0, T Y=0, T Z=0, double N=0.0, double C=1.0, int D=0) { set(X, Y, Z, N, C, D);}
00070     virtual ~Vector(void) {}
00071 
00072     T     norm2(void) { return (x*x + y*y + z*z);}
00073     double norm(void) { n = (double)sqrt(x*x + y*y + z*z); return n;}
00074     Vector<T> normalize(void);
00075 
00076     void   init(double C=1.0) { x = y = z = (T)0; n = 0.0; c = C; d = 0;}
00077     void   set(T X, T Y=0, T Z=0, double N=0.0, double C=1.0, int D=0);
00078 
00079     // 別名
00080     T      element1(void) { return x;}
00081     T      element2(void) { return y;}
00082     T      element3(void) { return z;}
00083     T      element(int i) { if(i==1) return x; else if(i==2) return y; else if(i==3) return z; else return (T)0;} 
00084 
00085     template <typename R> Vector<T>& operator = (const Vector<R> a) { x=(T)a.x; y=(T)a.y; z=(T)a.z; n=a.n; c=a.c; d=a.d; return *this;}
00086 };
00087 
00088 
00089 
00090 template <typename T> Vector<T> Vector<T>::normalize(void)
00091 {
00092     double nrm = (double)sqrt(x*x+y*y+z*z);
00093 
00094     if (nrm>=Zero_Eps) {
00095         x  = (T)(x/nrm);
00096         y  = (T)(y/nrm);
00097         z  = (T)(z/nrm);
00098         n  = 1.0;
00099     }
00100     else {
00101         init();
00102     }
00103 
00104     return *this;
00105 }
00106 
00107 
00108 
00114 template <typename T> inline void Vector<T>::set(T X, T Y, T Z, double N, double C, int D) 
00115 {
00116     x  = X;
00117     y  = Y;
00118     z  = Z;
00119     n  = N;
00120     c  = C;
00121     d  = D;
00122 
00123     if (n<=0.0) {
00124         n = (double)sqrt(x*x + y*y + z*z);
00125     }
00126 }
00127 
00128 
00129 
00131 // Vector オペレータ
00132 
00133 template <typename T> inline Vector<T> operator - (const Vector<T> a)
00134 { return Vector<T>(-a.x, -a.y, -a.z, a.n, a.c, a.d); }
00135 
00136 
00137 template <typename T> inline Vector<T> operator + (const Vector<T> a, const Vector<T> b)
00138 { return Vector<T>(a.x+b.x, a.y+b.y, a.z+b.z, 0.0, Min(a.c, b.c), a.d); }
00139 
00140 
00141 template <typename T, typename R> inline Vector<T> operator + (const Vector<T> a, R c)
00142 { return Vector<T>(a.x+(T)c, a.y+(T)c, a.z+(T)c, 0.0, a.c, a.d); }
00143 
00144 
00145 template <typename T, typename R> inline Vector<T> operator + (const R c, Vector<T> a)
00146 { return Vector<T>((T)c+a.x, (T)c+a.y, (T)c+a.z, 0.0, a.c, a.d); }
00147 
00148 
00149 template <typename T> inline Vector<T> operator - (const Vector<T> a, const Vector<T> b)
00150 { return Vector<T>(a.x-b.x, a.y-b.y, a.z-b.z, 0.0, Min(a.c, b.c), a.d); }
00151 
00152 
00153 template <typename T, typename R> inline Vector<T> operator - (const Vector<T> a, R c)
00154 { return Vector<T>(a.x-(T)c, a.y-(T)c, a.z-(T)c, 0.0, a.c, a.d); }
00155 
00156 
00157 template <typename T, typename R> inline Vector<T> operator - (R c, const Vector<T> a)
00158 { return Vector<T>((T)c-a.x, (T)c-a.y, (T)c-a.z, 0.0, a.c, a.d); }
00159 
00160 
00161 template <typename T, typename R> inline Vector<T> operator * (const R d, const Vector<T> a)
00162 { return Vector<T>((T)d*a.x, (T)d*a.y, (T)d*a.z, (T)d*a.n, a.c, a.d); }
00163 
00164 
00165 template <typename T, typename R> inline Vector<T> operator * (const Vector<T> a, const R d)
00166 { return Vector<T>(a.x*(T)d, a.y*(T)d, a.z*(T)d, a.n*(T)d, a.c, a.d); }
00167 
00168 
00169 template <typename T, typename R> inline Vector<T> operator / (const Vector<T> a, const R d)
00170 { return Vector<T>(a.x/(T)d, a.y/(T)d, a.z/(T)d, a.n/(T)d, a.c, a.d); }
00171 
00172 
00173 template <typename T, typename R> inline Vector<T> operator / (const R d, const Vector<T> a)
00174 {
00175     Vector<T> v((T)d/a.x, (T)d/a.y, (T)d/a.z, 0.0, a.c, a.d);
00176     v.norm();
00177     return v;
00178 }
00179 
00180 
00182 template <typename T> inline Vector<T> operator ^ (const Vector<T> a, const Vector<T> b)
00183 { return Vector<T>(a.y*b.z-a.z*b.y, a.z*b.x-a.x*b.z, a.x*b.y-a.y*b.x, 0.0, Min(a.c, b.c)); }
00184 
00185 
00187 template <typename T> inline T operator * (const Vector<T> a, const Vector<T> b)
00188 { return (a.x*b.x + a.y*b.y + a.z*b.z); }
00189 
00190 
00191 template <typename T> inline bool operator == (const Vector<T> v1, const Vector<T> v2)
00192 { return (v1.x==v2.x && v1.y==v2.y && v1.z==v2.z); }
00193 
00194 
00195 template <typename T> inline bool operator != (const Vector<T> v1, const Vector<T> v2)
00196 { return (v1.x!=v2.x || v1.y!=v2.y || v1.z!=v2.z); }
00197 
00198 //
00199 template <typename T> inline Vector<T> MidPoint(const Vector<T> a, const Vector<T> b)
00200 { return 0.5*(a+b); }
00201 
00202 
00203 //
00205 template <typename T> inline double VectorDist(const Vector<T> a, const Vector<T> b)
00206 {
00207     //return (b-a).n;
00208     return (b-a).norm();
00209 }
00210 
00211 
00212 
00213 template <typename T> inline bool operator < (const Vector<T> v1, const Vector<T> v2)
00214 {
00215     if (v1.x != v2.x) return v1.x < v2.x;
00216     if (v1.y != v2.y) return v1.y < v2.y;
00217     if (v1.z != v2.z) return v1.z < v2.z;
00218     return false;
00219 }
00220 
00221 
00222 
00224 
00229 template <typename T> inline bool same_vector(Vector<T> v1, Vector<T> v2)
00230 {
00231     T dx = v1.x - v2.x;
00232     T dy = v1.y - v2.y;
00233     T dz = v1.z - v2.z;
00234     T d2 = dx*dx + dy*dy + dz*dz;
00235     double t2 = Vector_Tolerance*Vector_Tolerance;
00236 
00237     if ((double)d2>t2) return false;
00238     return true;
00239 }
00240 
00241 
00242 
00243 template <typename T> inline Vector<T>* dupVector(Vector<T>* a, int n)
00244 {
00245     Vector<T>* v = (Vector<T>*)malloc(sizeof(Vector<T>)*n);
00246     if (v==NULL) return NULL;
00247 
00248     for (int i=0; i<n; i++) v[i] = a[i];
00249 
00250     return v;
00251 }
00252 
00253 
00254 // angle between a and b
00255 template <typename T> inline double VectorAngle(Vector<T> a, Vector<T> b)
00256 {
00257     a.normalize();
00258     b.normalize();
00259     if (a.n<Zero_Eps || b.n<Zero_Eps) return 0.0;
00260 
00261     double cs = a*b;
00262     if      (cs>=1.0)  return 0.0;
00263     else if (cs<=-1.0) return PI;
00264     
00265     return acos(a*b);
00266 }
00267 
00268 
00269 template <typename T> inline double VectorAngle(Vector<T> a, Vector<T> b,  Vector<T> c)
00270 {
00271     return VectorAngle(b-a, c-b);
00272 }
00273 
00274 
00275 //
00277 template <typename T> inline Vector<T> NewellMethod(Vector<T> v1, Vector<T> v2, Vector<T> v3)
00278 {
00279     Vector<T> vect;
00280 
00281     vect.x = (v1.y-v2.y)*(v1.z+v2.z) + (v2.y-v3.y)*(v2.z+v3.z) + (v3.y-v1.y)*(v3.z+v1.z);
00282     vect.y = (v1.z-v2.z)*(v1.x+v2.x) + (v2.z-v3.z)*(v2.x+v3.x) + (v3.z-v1.z)*(v3.x+v1.x);
00283     vect.z = (v1.x-v2.x)*(v1.y+v2.y) + (v2.x-v3.x)*(v2.y+v3.y) + (v3.x-v1.x)*(v3.y+v1.y);
00284     vect.n = 0.0;
00285     vect.c = Min(v1.c, v2.c);
00286     vect.c = Min(v3.c, vect.c);
00287 
00288     return vect;
00289 }
00290 
00291 
00292 template <typename T> inline Vector<T> NewellMethod3(Vector<T> v1, Vector<T> v2, Vector<T> v3)
00293 {
00294     return NewellMethod(v1, v2, v3);
00295 }
00296 
00297 
00298 template <typename T> inline Vector<T> NewellMethod4(Vector<T> v1, Vector<T> v2, Vector<T> v3, Vector<T> v4)
00299 {
00300     Vector<T> vect;
00301 
00302     vect.x = (v1.y-v2.y)*(v1.z+v2.z) + (v2.y-v3.y)*(v2.z+v3.z) + (v3.y-v4.y)*(v3.z+v4.z) + (v4.y-v1.y)*(v4.z+v1.z);
00303     vect.y = (v1.z-v2.z)*(v1.x+v2.x) + (v2.z-v3.z)*(v2.x+v3.x) + (v3.z-v4.z)*(v3.x+v4.x) + (v4.z-v1.z)*(v4.x+v1.x);
00304     vect.z = (v1.x-v2.x)*(v1.y+v2.y) + (v2.x-v3.x)*(v2.y+v3.y) + (v3.x-v4.x)*(v3.y+v4.y) + (v4.x-v1.x)*(v4.y+v1.y);
00305     vect.n = 0.0;
00306     vect.c = Min(v1.c, v2.c);
00307     vect.c = Min(v3.c, vect.c);
00308     vect.c = Min(v4.c, vect.c);
00309     return vect;
00310 }
00311 
00312 
00313 
00324 template <typename T> inline Vector<T> BSplineInterp4(Vector<T> p0, Vector<T> p1, double t)
00325 {
00326     Vector<T> q0 = 2*p0 - p1;
00327 //  Vector<T> q1 = p0;
00328 //  Vector<T> q2 = p1;
00329     Vector<T> q3 = 2*p1 - p0;
00330 
00331     double t2 = t*t;
00332     double t3 = t2*t;
00333     
00334     double c0 = (t2 - t)*0.5 + (1.0 - t3)*0.1666666666666667;   // 1/6*(1-t)^3
00335     double c1 = t3*0.5 - t2 + 0.6666666666666667;               // 1/2*t3 - t2 + 2/3
00336     double c2 = (t + t2 - t3)*0.5 + 0.1666666666666667;         // -1/2*t3 + 1/2*t2 + 1/2*t + 1/6
00337     double c3 = t3*0.1666666666666667;                          // 1/6*t3
00338 
00339     Vector<T> vect = c0*q0 + c1*p0 + c2*p1 + c3*q3;
00340 
00341     return vect;
00342 }
00343 
00344 
00345 
00346 
00348 // 3D多角形座標
00349 //
00350 
00357 template <typename T=int> class DllExport  PCoordinate
00358 {
00359 public:
00360     int  dim;               // 角数
00361     Vector<T>** point;      // データ
00362 
00363 public:
00364     
00365     PCoordinate(void)  { dim = 0; point = NULL;}
00366     PCoordinate(int n) { init(n);}
00367 
00368     virtual ~PCoordinate(void) {}
00369 
00370     //
00371     void  init(int n) {
00372         dim   = 0;
00373         point = NULL;
00374 
00375         if (n>0) {
00376             point = (Vector<T>**)malloc(sizeof(Vector<T>*)*n);
00377             if (point!=NULL) {
00378                 dim = n;
00379                 for (int i=0; i<dim; i++) point[i] = new Vector<T>();
00380             }
00381         }
00382     }
00383 
00384     void  set(int m, T x, T y=0, T z=0, double n=0.0) {
00385         if (m>=0 && m<dim && point!=NULL) point[m]->set(x, y, z, n); 
00386     }       
00387     
00388     void  clear(void) {
00389         if (point!=NULL) {
00390             for (int i=0; i<dim; i++) point[i]->set((T)0);
00391         }
00392     }
00393 
00394     void  free(void) {
00395         if (point!=NULL) {
00396             for (int i=0; i<dim; i++) delete(point[i]);
00397             ::free(point);
00398             point = NULL;
00399         }
00400         dim = 0;
00401     }
00402 
00403 };
00404 
00405 
00406 
00407 
00409 // 境界構造体
00410 //
00411 
00417 template <typename T=int> class DllExport RBound 
00418 {
00419 public:
00420     T   xmin;       
00421     T   xmax;       
00422     T   ymin;       
00423     T   ymax;       
00424     T   zmin;       
00425     T   zmax;       
00426     T   tmin;       
00427     T   tmax;       
00428 
00429 public:
00430 
00431     RBound(T XMin=(T)0, T XMax=(T)0, T YMin=(T)0, T YMax=(T)0, T ZMin=(T)0, T ZMax=(T)0, T TMin=(T)0, T TMax=(T)0) {
00432         set(XMin, XMax, YMin, YMax, ZMin, ZMax, TMin, TMax);
00433     }
00434 
00435     virtual ~RBound() {}
00436 
00437 
00438     //
00439     void  set(T XMin=(T)0, T XMax=(T)0, T YMin=(T)0, T YMax=(T)0, T ZMin=(T)0, T ZMax=(T)0, T TMin=(T)0, T TMax=(T)0){
00440         xmin = XMin;
00441         ymin = YMin;
00442         zmin = ZMin;
00443         xmax = XMax;
00444         ymax = YMax;
00445         zmax = ZMax;
00446         tmin = TMin;
00447         tmax = TMax;
00448     }       
00449 
00450     void init() { xmin = xmax = ymin = ymax = zmin = zmax = tmin = tmax = (T)0;}
00451 
00453     // 拡大
00454 
00456     void  enlarge(T f) {
00457         xmin -= f;
00458         ymin -= f;
00459         zmin -= f;
00460         xmax += f;
00461         ymax += f;
00462         zmax += f;
00463     }
00464 
00465     void  multiple(T f) {
00466         xmin *= f;
00467         ymin *= f;
00468         zmin *= f;
00469         xmax *= f;
00470         ymax *= f;
00471         zmax *= f;
00472     }
00473 
00475     // 融合
00476 
00478     void  fusion(RBound<T> bound) {
00479         xmin = Min(xmin, bound.xmin);
00480         ymin = Min(ymin, bound.ymin);
00481         zmin = Min(zmin, bound.zmin);
00482         xmax = Max(xmax, bound.xmax);
00483         ymax = Max(ymax, bound.ymax);
00484         zmax = Max(zmax, bound.zmax);
00485     }
00486 
00488     void  fusion(Vector<T> vect){
00489         xmin = Min(xmin, vect.x);
00490         ymin = Min(ymin, vect.y);
00491         zmin = Min(zmin, vect.z);
00492         xmax = Max(xmax, vect.x);
00493         ymax = Max(ymax, vect.y);
00494         zmax = Max(zmax, vect.z);
00495     }
00496 
00498     void  fusion(T x, T y, T z){
00499         xmin = Min(xmin, x);
00500         ymin = Min(ymin, y);
00501         zmin = Min(zmin, z);
00502         xmax = Max(xmax, x);
00503         ymax = Max(ymax, y);
00504         zmax = Max(zmax, z);
00505     }
00506 
00508     // 共通領域
00509 
00511     void  commonarea(RBound<T> bound) {
00512         xmin = Max(xmin, bound.xmin);
00513         ymin = Max(ymin, bound.ymin);
00514         zmin = Max(zmin, bound.zmin);
00515         xmax = Min(xmax, bound.xmax);
00516         ymax = Min(ymax, bound.ymax);
00517         zmax = Min(zmax, bound.zmax);
00518     }
00519 
00522     void  cutdown(RBound<T> bound) {
00523         xmin += bound.xmin;
00524         ymin += bound.ymin;
00525         zmin += bound.zmin;
00526         xmax += bound.xmin;
00527         ymax += bound.ymin;
00528         zmax += bound.zmin;
00529     }
00530 
00531     //   
00532     void  cutdown(Vector<T> vect){
00533         xmin += vect.x;
00534         ymin += vect.y;
00535         zmin += vect.z;
00536         xmax += vect.x;
00537         ymax += vect.y;
00538         zmax += vect.z;
00539     }
00540 
00541     //  
00542     void  cutdown(T x, T y, T z){
00543         xmin += x;
00544         ymin += y;
00545         zmin += z;
00546         xmax += x;
00547         ymax += y;
00548         zmax += z;
00549     }
00550 
00552     bool  outofBounds(Vector<T> vect) {
00553         return  vect.x < xmin || vect.x > xmax ||
00554                 vect.y < ymin || vect.y > ymax ||
00555                 vect.z < zmin || vect.z > zmax;
00556     }
00557 
00559     bool  outofBounds(T x, T y, T z) {
00560         return  x < xmin || x > xmax ||
00561                 y < ymin || y > ymax ||
00562                 z < zmin || z > zmax;
00563     }
00564 };
00565 
00566 
00567 
00576 template <typename T> inline bool disJunctBounds(RBound<T> b1, RBound<T> b2)
00577 {
00578     return  (b1.xmin >= b2.xmax) || (b2.xmin >= b1.xmax) || 
00579             (b1.ymin >= b2.ymax) || (b2.ymin >= b1.ymax) || 
00580             (b1.zmin >= b2.zmax) || (b2.zmin >= b1.zmax);
00581 }
00582 
00583 
00584 
00585 
00587 // UVMap
00588 //
00589 
00596 template <typename T=double> class DllExport UVMap
00597 {
00598 public:
00599     T   u;
00600     T   v;
00601     int d;      
00602 
00603 public:
00604     UVMap(T U=0, T V=0, int d=0) { set(U, V, d);}
00605     virtual ~UVMap(void) {}
00606 
00607     void   init(void) { u = v = (T)0; d = 0;}
00608     void   set(T U, T V=0, int D=0) { u = U, v = V; d = D;}
00609 
00610     UVMap flip(void)  { u = (T)(1.0-u); v = (T)(1.0-v); return *this;}
00611     UVMap flipV(void) { v = (T)(1.0-v); return *this;}
00612     UVMap flipU(void) { u = (T)(1.0-u); return *this;}
00613 
00614     template <typename R> UVMap<T>& operator = (const UVMap<R> a) { u=(T)a.u; v=(T)a.v; d=a.d; return *this;}
00615 };
00616 
00617 
00618 
00619 
00621 // UVMap オペレータ
00622 
00623 template <typename T> inline UVMap<T> operator - (const UVMap<T> a)
00624 { return UVMap<T>(-a.u, -a.v, a.d); }
00625 
00626 
00627 template <typename T> inline UVMap<T> operator + (const UVMap<T> a, const UVMap<T> b)
00628 { return UVMap<T>(a.u+b.u, a.v+b.v, a.d); }
00629 
00630 
00631 template <typename T, typename R> inline UVMap<T> operator + (const UVMap<T> a, R c)
00632 { return UVMap<T>(a.u+(T)c, a.v+(T)c, a.d); }
00633 
00634 
00635 template <typename T, typename R> inline UVMap<T> operator + (const R c, UVMap<T> a)
00636 { return UVMap<T>((T)c+a.u, (T)c+a.v, a.d); }
00637 
00638 
00639 template <typename T> inline UVMap<T> operator - (const UVMap<T> a, const UVMap<T> b)
00640 { return UVMap<T>(a.u-b.u, a.v-b.v, a.d); }
00641 
00642 
00643 template <typename T, typename R> inline UVMap<T> operator - (const UVMap<T> a, R c)
00644 { return UVMap<T>(a.u-(T)c, a.v-(T)c, a.d); }
00645 
00646 
00647 template <typename T, typename R> inline UVMap<T> operator - (const R c, UVMap<T>a)
00648 { return UVMap<T>((T)c-a.u, (T)c-a.v, a.d); }
00649 
00650 
00651 template <typename T, typename R> inline UVMap<T> operator * (const R d, const UVMap<T> a)
00652 { return UVMap<T>((T)d*a.u, (T)d*a.v, a.d); }
00653 
00654 
00655 template <typename T, typename R> inline UVMap<T> operator * (const UVMap<T> a, const R d)
00656 { return UVMap<T>(a.u*(T)d, a.v*(T)d, a.d); }
00657 
00658 
00659 template <typename T, typename R> inline UVMap<T> operator / (const UVMap<T> a, const R d)
00660 { return UVMap<T>(a.u/(T)d, a.v/(T)d, a.d); }
00661 
00662 
00663 template <typename T, typename R> inline UVMap<T> operator / (const R d, const UVMap<T> a)
00664 { return UVMap<T>((T)d/a.u, (T)d/a.v, a.d); }
00665 
00666 
00667 template <typename T> inline bool operator == (const UVMap<T> a, const UVMap<T> b)
00668 { return (a.u==b.u && a.v==b.v); }
00669 
00670 
00671 template <typename T> inline bool operator != (const UVMap<T> a, const UVMap<T> b)
00672 { return (a.u!=b.u || a.v!=b.v); }
00673 
00674 
00675 
00676 }       // namespace
00677 
00678 #endif 
00679 
00680