hpoint_nd.h Source File

00001 /*============================================================================
00002         File: hpoint_nd.h
00003      Purpose: 
00004     Revision: $Id: hpoint_nd.h,v 1.3 2002/05/22 17:06:47 philosophil Exp $
00005   Created by: Philippe Lavoie          (26 January 1999)
00006  Modified by: Martin Schuerch
00007 
00008  Copyright notice:
00009           Copyright (C) 1996-1999 Philippe Lavoie
00010  
00011           This library is free software; you can redistribute it and/or
00012           modify it under the terms of the GNU Library General Public
00013           License as published by the Free Software Foundation; either
00014           version 2 of the License, or (at your option) any later version.
00015  
00016           This library is distributed in the hope that it will be useful,
00017           but WITHOUT ANY WARRANTY; without even the implied warranty of
00018           MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019           Library General Public License for more details.
00020  
00021           You should have received a copy of the GNU Library General Public
00022           License along with this library; if not, write to the Free
00023           Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
00024 =============================================================================*/
00025 #ifndef _Matrix_hpointnD_h_
00026 #define _Matrix_hpointnD_h_
00027 
00028 #include "matrix_global.h"
00029 #include "point_nd.h"
00030 
00031 
00032 namespace PLib {
00033 
00057   template <class T, int N>
00058   struct HPoint_nD {
00059     // only N=2,3 are defined for now
00060     //HPoint_nD() = 0 ; 
00061   };
00062 
00063 
00064 
00065   template <>
00066   struct HPoint_nD<float,3> {
00067     typedef float T;
00068     T *data ; 
00069     int created; // usefull to change the data pointer
00070 
00071     HPoint_nD(): data(new T[4]), created(1) { x() = y() = z() = w() = T(); }
00072     HPoint_nD(T *d,int c): data(d), created(c) { ; } 
00073     HPoint_nD(T a): data(new T[4]), created(1) { x() = y() = z() = w() = a; }
00074     HPoint_nD(T X, T Y, T Z, T W):  data(new T[4]), created(1) { x()=X ; y()=Y; z()=Z; w()=W;}
00075     HPoint_nD(const HPoint_nD& a): data(new T[4]), created(1) { memcpy((void*)data,(void*)a.data,4*sizeof(T)); }
00076     HPoint_nD(const Point_nD<T,3>& a): data(new T[4]), created(1) { memcpy((void*)data,(void*)a.data,(3)*sizeof(T));  data[3] = T(1); }
00077     ~HPoint_nD() { if(created) if(data) delete []data ; }
00078 
00079     inline T& x() { return data[0] ;}
00080     inline T x() const { return data[0] ; }
00081     inline T& y() { return data[1] ;}
00082     inline T y() const { return data[1] ;}
00083     inline T& z() { return data[2] ;}
00084     inline T z() const { return data[2] ;}
00085     inline T& w() { return data[3] ;}
00086     inline T w() const { return data[3] ;}
00087   
00088     HPoint_nD& operator=(const HPoint_nD& v) { data[0]=v.data[0]; data[1]=v.data[1]; data[2]=v.data[2]; data[3]=v.data[3]; return *this;}
00089     HPoint_nD& operator=(const Point_nD<T,3>& v) {x()=v.x() ; y()=v.y() ; z()=v.z() ; w()=1.0 ; return *this;}
00090     HPoint_nD& operator=(const T v) { x() = y() = z() = w() = v ; return *this ;} ;
00091     HPoint_nD& operator+=(const HPoint_nD& v) { data[0]+=v.data[0]; data[1]+=v.data[1]; data[2]+=v.data[2]; data[3]+=v.data[3]; return *this;}
00092     HPoint_nD& operator-=(const HPoint_nD& v) { data[0]-=v.data[0]; data[1]-=v.data[1]; data[2]-=v.data[2]; data[3]-=v.data[3]; return *this;}
00093 
00094     HPoint_nD& operator*=(T v) { data[0]*=v; data[1]*=v; data[2]*=v; data[3]*=v; return *this; }
00095     HPoint_nD& operator/=(T v) { data[0]/=v; data[1]/=v; data[2]/=v; data[3]/=v; return *this; }
00096 
00097     void move(const Point_nD<T,3>& m) {
00098       for(int i=3-1;i>=0;--i) {
00099         data[i] += m.data[i]*data[3];
00100       }
00101     }
00102     Point_nD<T,3> projectW() { return Point_nD<T,3>(x(),y(),z()) ; }
00103   };
00104 
00105 
00106   template <>
00107   struct HPoint_nD<double,3> {
00108     typedef double T;
00109     T *data ; 
00110     int created; // usefull to change the data pointer
00111 
00112     HPoint_nD(): data(new T[4]), created(1) { x() = y() = z() = w() = T(); }
00113     HPoint_nD(T *d,int c): data(d), created(c) { ; } 
00114     HPoint_nD(T a): data(new T[4]), created(1) { x() = y() = z() = w() = a; }
00115     HPoint_nD(T X, T Y, T Z, T W):  data(new T[4]), created(1) { x()=X ; y()=Y; z()=Z; w()=W;}
00116     HPoint_nD(const HPoint_nD& a): data(new T[4]), created(1) { memcpy((void*)data,(void*)a.data,4*sizeof(T)); }
00117     HPoint_nD(const Point_nD<T,3>& a): data(new T[4]), created(1) { memcpy((void*)data,(void*)a.data,(3)*sizeof(T));  data[3] = T(1); }
00118     ~HPoint_nD() { if(created) if(data) delete []data ; }
00119 
00120     inline T& x() { return data[0] ;}
00121     inline T x() const { return data[0] ; }
00122     inline T& y() { return data[1] ;}
00123     inline T y() const { return data[1] ;}
00124     inline T& z() { return data[2] ;}
00125     inline T z() const { return data[2] ;}
00126     inline T& w() { return data[3] ;}
00127     inline T w() const { return data[3] ;}
00128   
00129     HPoint_nD& operator=(const HPoint_nD& v) { data[0]=v.data[0]; data[1]=v.data[1]; data[2]=v.data[2]; data[3]=v.data[3]; return *this;}
00130     HPoint_nD& operator=(const Point_nD<T,3>& v) {x()=v.x() ; y()=v.y() ; z()=v.z() ; w()=1.0 ; return *this;}
00131     HPoint_nD& operator=(const T v) { x() = y() = z() = w() = v ; return *this ;} ;
00132     HPoint_nD& operator+=(const HPoint_nD& v) { data[0]+=v.data[0]; data[1]+=v.data[1]; data[2]+=v.data[2]; data[3]+=v.data[3]; return *this;}
00133     HPoint_nD& operator-=(const HPoint_nD& v) { data[0]-=v.data[0]; data[1]-=v.data[1]; data[2]-=v.data[2]; data[3]-=v.data[3]; return *this;}
00134 
00135     HPoint_nD& operator*=(T v) { data[0]*=v; data[1]*=v; data[2]*=v; data[3]*=v; return *this; }
00136     HPoint_nD& operator/=(T v) { data[0]/=v; data[1]/=v; data[2]/=v; data[3]/=v; return *this; }
00137 
00138     void move(const Point_nD<T,3>& m) {
00139       for(int i=3-1;i>=0;--i) {
00140         data[i] += m.data[i]*data[3];
00141       }
00142     }
00143     Point_nD<T,3> projectW() { return Point_nD<T,3>(x(),y(),z()) ; }
00144   };
00145 
00146 
00147 
00148 
00149 
00161   template <>
00162   struct HPoint_nD<float,2> {
00163     typedef float T;
00164     T *data ; 
00165     int created; // usefull to change the data pointer
00166 
00167     HPoint_nD(): data(new T[3]), created(1) { x() = y() = w() = T(); }
00168     HPoint_nD(T *d,int c): data(d), created(c) { ; } 
00169     HPoint_nD(T a): data(new T[3]), created(1) { x() = y() = w() = a; }
00170     HPoint_nD(T X, T Y, T W):  data(new T[3]), created(1) { x()=X ; y()=Y; w()=W;}
00171     HPoint_nD(T X, T Y, T Z, T W):  data(new T[3]), created(1) { x()=X ; y()=Y; w()=W;}
00172     HPoint_nD(const HPoint_nD& a): data(new T[3]), created(1) { memcpy((void*)data,(void*)a.data,3*sizeof(T)); }
00173     HPoint_nD(const Point_nD<T,2>& a): data(new T[3]), created(1) { memcpy((void*)data,(void*)a.data,(2)*sizeof(T));  data[2] = T(1); }
00174     ~HPoint_nD() { if(created) if(data) delete []data ; }
00175 
00176     inline T& x() { return data[0] ;}
00177     inline T x() const { return data[0] ; }
00178     inline T& y() { return data[1] ;}
00179     inline T y() const { return data[1] ;}
00180     inline T& z() { return  dumbVar;}
00181     inline T z() const { return T(0) ;}
00182     inline T& w() { return data[2] ;}
00183     inline T w() const { return data[2] ;}
00184   
00185 
00186     HPoint_nD<T,2>& operator=(const HPoint_nD<T,2>& v) { data[0]=v.data[0]; data[1]=v.data[1]; data[2]=v.data[2]; return *this;}
00187     HPoint_nD<T,2>& operator=(const Point_nD<T,2>& v) {x()=v.x() ; y()=v.y() ; w()=1.0 ; return *this;}
00188     HPoint_nD<T,2>& operator=(const T v) { x() = y() = w() = v ; return *this ;} ;
00189 
00190     HPoint_nD<T,2>& operator*=(T v) { data[0]*=v; data[1]*=v; data[2]*=v; return *this; }
00191     HPoint_nD<T,2>& operator/=(T v) { data[0]/=v; data[1]/=v; data[2]/=v; return *this; }
00192 
00193     HPoint_nD<T,2>& operator-=(const HPoint_nD<T,2>& v) { data[0]-=v.data[0]; data[1]-=v.data[1]; data[2]-=v.data[2]; return *this;}
00194     HPoint_nD<T,2>& operator+=(const HPoint_nD<T,2>& v) { data[0]+=v.data[0]; data[1]+=v.data[1]; data[2]+=v.data[2]; return *this;}
00195 
00196   void move(const Point_nD<T,2>& m) {
00197     for(int i=2-1;i>=0;--i) {
00198       data[i] += m.data[i]*data[2];
00199     }
00200   }
00201   Point_nD<T,2> projectW() { return Point_nD<T,2>(x(),y()) ; }
00202 
00203   protected:
00204     static T dumbVar ;
00205   };
00206 
00207 
00208   template <>
00209   struct HPoint_nD<double,2> {
00210     typedef double T;
00211     T *data ; 
00212     int created; // usefull to change the data pointer
00213 
00214     HPoint_nD(): data(new T[3]), created(1) { x() = y() = w() = T(); }
00215     HPoint_nD(T *d,int c): data(d), created(c) { ; } 
00216     HPoint_nD(T a): data(new T[3]), created(1) { x() = y() = w() = a; }
00217     HPoint_nD(T X, T Y, T W):  data(new T[3]), created(1) { x()=X ; y()=Y; w()=W;}
00218     HPoint_nD(T X, T Y, T Z, T W):  data(new T[3]), created(1) { x()=X ; y()=Y; w()=W;}
00219     HPoint_nD(const HPoint_nD& a): data(new T[3]), created(1) { memcpy((void*)data,(void*)a.data,3*sizeof(T)); }
00220     HPoint_nD(const Point_nD<T,2>& a): data(new T[3]), created(1) { memcpy((void*)data,(void*)a.data,(2)*sizeof(T));  data[2] = T(1); }
00221     ~HPoint_nD() { if(created) if(data) delete []data ; }
00222 
00223     inline T& x() { return data[0] ;}
00224     inline T x() const { return data[0] ; }
00225     inline T& y() { return data[1] ;}
00226     inline T y() const { return data[1] ;}
00227     inline T& z() { return  dumbVar;}
00228     inline T z() const { return T(0) ;}
00229     inline T& w() { return data[2] ;}
00230     inline T w() const { return data[2] ;}
00231   
00232 
00233     HPoint_nD<T,2>& operator=(const HPoint_nD<T,2>& v) { data[0]=v.data[0]; data[1]=v.data[1]; data[2]=v.data[2]; return *this;}
00234     HPoint_nD<T,2>& operator=(const Point_nD<T,2>& v) {x()=v.x() ; y()=v.y() ; w()=1.0 ; return *this;}
00235     HPoint_nD<T,2>& operator=(const T v) { x() = y() = w() = v ; return *this ;} ;
00236 
00237     HPoint_nD<T,2>& operator*=(T v) { data[0]*=v; data[1]*=v; data[2]*=v; return *this; }
00238     HPoint_nD<T,2>& operator/=(T v) { data[0]/=v; data[1]/=v; data[2]/=v; return *this; }
00239 
00240     HPoint_nD<T,2>& operator-=(const HPoint_nD<T,2>& v) { data[0]-=v.data[0]; data[1]-=v.data[1]; data[2]-=v.data[2]; return *this;}
00241     HPoint_nD<T,2>& operator+=(const HPoint_nD<T,2>& v) { data[0]+=v.data[0]; data[1]+=v.data[1]; data[2]+=v.data[2]; return *this;}
00242 
00243     void move(const Point_nD<T,2>& m) {
00244       for(int i=2-1;i>=0;--i) {
00245         data[i] += m.data[i]*data[2];
00246       }
00247     }
00248 
00249     Point_nD<T,2> projectW() { return Point_nD<T,2>(x(),y()) ; }
00250 
00251   protected:
00252     static T dumbVar ;
00253   };
00254 
00255 
00256 
00268   template <class T, const int D>
00269   class NoInitHPoint_nD : public HPoint_nD<T,D> {
00270   public:
00271     NoInitHPoint_nD() : HPoint_nD<T,D>(0,0) { ; }
00272   };
00273 
00274 
00286   template <class T>
00287   class NoInitHPoint_2D : public HPoint_nD<T,2> {
00288   public:
00289     NoInitHPoint_2D() : HPoint_nD<T,2>(0,0) { ; }
00290   };
00291 
00292 
00293 
00307   template <class T, int D>
00308   inline HPoint_nD<T,D> operator+(const HPoint_nD<T,D>&a,const HPoint_nD<T,D>&b) {
00309     HPoint_nD<T,D> sum(a) ;
00310     sum += b ;
00311     return sum ;
00312   }
00313 
00314 
00315   /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
00316         Routine: operator- --- the substraction operator for a HPoint_3D
00317      The substraction operator for a HPoint_3D
00318           Input: a --> the first vector
00319                  b --> the second vector
00320          Output: $a-b$
00321    Restrictions:
00322      author Philippe Lavoie (24 January 1997)
00323     Modified by:
00324    * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
00325   template <class T, int D>
00326   inline HPoint_nD<T,D> operator-(const HPoint_nD<T,D>&a,const HPoint_nD<T,D>&b) {
00327     HPoint_nD<T,D> diff(a) ;
00328     diff -= b ;
00329     return diff ;
00330   }
00331 
00332 
00333   /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
00334         Routine: dot --- The dot product between two points
00335      The dot product between two points in 4D.
00336           Input: a --> the first vector
00337                  b --> the second vector
00338          Output: $a.b$
00339    Restrictions:
00340      author Philippe Lavoie (24 January 1997)
00341     Modified by:
00342    * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
00343   template <class T>
00344   inline T dot(const HPoint_nD<T,3>&a,const HPoint_nD<T,3>&b) {
00345     return a.x()*b.x()+a.y()*b.y()+a.z()*b.z()+a.w()*b.w() ;
00346   }
00347 
00348   template <class T>
00349   inline T dot(const HPoint_nD<T,2>&a,const HPoint_nD<T,2>&b) {
00350     return a.x()*b.x()+a.y()*b.y()+a.w()*b.w() ;
00351   }
00352 
00353 
00354   /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
00355         Routine: operator* --- multiplying a point with a float
00356      Multiplying a point with a float
00357           Input: b --> the point to multiply
00358                  a --> the floating point value to multiply with
00359          Output: $ab$
00360    Restrictions:
00361      author Philippe Lavoie (24 January 1997)
00362     Modified by:
00363    * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
00364   template <class T>
00365   inline HPoint_nD<T,3> operator*(const HPoint_nD<T,3>&b, T a) {
00366     HPoint_nD<T,3> mul(b.x()*a,b.y()*a,b.z()*a,b.w()*a) ;
00367     return mul;
00368   }
00369 
00370   template <class T>
00371   inline HPoint_nD<T,2> operator*(const HPoint_nD<T,2>&b, T a) {
00372     HPoint_nD<T,3> mul(b.x()*a,b.y()*a,b.w()*a) ;
00373     return mul;
00374   }
00375 
00376   /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
00377         Routine: operator* --- multiplying a point with a float
00378      Multiplying a point with a float
00379           Input: a --> the floating point value to multiply with
00380                  b --> the point to multiply
00381          Output: $ab$
00382    Restrictions:
00383      author Philippe Lavoie (24 January 1997)
00384     Modified by:
00385    * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
00386   template <class T>
00387   inline HPoint_nD<T,3> operator*( T a, const HPoint_nD<T,3>&b) {
00388     HPoint_nD<T,3> mul(b.x()*a,b.y()*a,b.z()*a,b.w()*a) ;
00389     return mul;
00390     // the following code isn't as good because no convresion of double and float are made
00391     //HPoint_nD<T,N> mul(b) ;
00392     //return mul *= a;
00393   }
00394 
00395   inline HPoint_nD<float,3> operator*( double a, const HPoint_nD<float,3>&b) {
00396     HPoint_nD<float,3> mul(b.x()*a,b.y()*a,b.z()*a,b.w()*a) ;
00397     return mul;
00398   }
00399 
00400 
00401   template <class T>
00402   inline HPoint_nD<T,2> operator*(T a, const HPoint_nD<T,2>&b) {
00403     HPoint_nD<T,2> mul(b.x()*a,b.y()*a,b.w()*a) ;
00404     return mul;
00405   }
00406 
00407   inline HPoint_nD<float,2> operator*(double a, const HPoint_nD<float,2>&b) {
00408     HPoint_nD<float,2> mul(b.x()*a,b.y()*a,b.w()*a) ;
00409     return mul;
00410   }
00411   
00412 
00413   /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
00414         Routine: operator* --- The dot product between two points
00415      The dot product between two points in 4D.
00416           Input: a --> the first vector
00417                  b --> the second vector
00418          Output: $a.b$
00419    Restrictions:
00420      author Philippe Lavoie (24 January 1997)
00421     Modified by:
00422    * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
00423   template <class T>
00424   inline T operator*(const HPoint_nD<T,3>&a,const HPoint_nD<T,3>&b) {
00425     return a.x()*b.x()+a.y()*b.y()+a.z()*b.z()+a.w()*b.w() ;
00426   }
00427 
00428   template <class T>
00429   inline T operator*(const HPoint_nD<T,2>&a,const HPoint_nD<T,2>&b) {
00430     return a.x()*b.x()+a.y()*b.y()+a.w()*b.w() ;
00431   }
00432 
00433 
00434 
00435    /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
00436         Routine: operator/ --- Divide a point by a float
00437      Divide a point by a float
00438           Input: a --> the point in 4D
00439                  b --> the floating point value
00440          Output: $a/b$
00441    Restrictions:
00442      author Philippe Lavoie (24 January 1997)
00443     Modified by:
00444    * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
00445   template <class T>
00446   inline HPoint_nD<T,3> operator/(const HPoint_nD<T,3>&a,const T b) {
00447     HPoint_nD<T,3> div(a.x()/b,a.y()/b,a.z()/b,a.w()/b) ;
00448     return div ;
00449   }
00450 
00451   inline HPoint_nD<float,3> operator/(const HPoint_nD<float,3>&a,const double b) {
00452     HPoint_nD<float,3> div(float(a.x()/b),float(a.y()/b),float(a.z()/b),float(a.w()/b)) ;
00453     return div ;
00454   }
00455 
00456   template <class T>
00457   inline HPoint_nD<T,2> operator/(const HPoint_nD<T,2>&a,const T b) {
00458     HPoint_nD<T,2> div(a.x()/b,a.y()/b,a.w()/b) ;
00459     return div ;
00460   }
00461 
00462   inline HPoint_nD<float,2> operator/(const HPoint_nD<float,2>&a,const double b) {
00463     HPoint_nD<float,2> div(float(a.x()/b),float(a.y()/b),float(a.w()/b)) ;
00464     return div ;
00465   }
00466 
00467 
00468   /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
00469         Routine: operator== --- The equality operator with a float
00470      The equality operator with a float. Verifies if every elements
00471                  of $a$ are equal to $b$.
00472           Input: a --> the control point
00473                  b --> the floating point value
00474          Output: 1 if every elements of $a$ is equal to $b$, 0 otherwise
00475    Restrictions:
00476      author Philippe Lavoie (24 January 1997)
00477     Modified by:
00478    * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
00479   template <class T, int D>
00480   inline int operator==(const HPoint_nD<T,D>& a, T b) {
00481     int r ;
00482     r = 1 ;
00483     for(int i=3;i>=0;--i)
00484       r = r && a.data[i]==b ;
00485     return r ;
00486   }
00487 
00488   /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
00489         Routine: operator!= --- the inequality operator between two points
00490      The inequality operator between two points
00491           Input: a --> the first vector
00492                  b --> the second vector
00493          Output: 1 if they are different, 0 otherwise
00494    Restrictions: This might be of limited value since comparing two floating
00495                  point numbers doesn't usually give the proper result.
00496      author Philippe Lavoie (24 January 1997)
00497     Modified by:
00498    * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
00499   template <class T, int D>
00500   inline int operator!=(const HPoint_nD<T,D>& a, const HPoint_nD<T,D>& b){
00501     return !(a==b) ;
00502   }
00503 
00504   /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
00505         Routine: operator== --- the equality operator between two points
00506      The equality operator between two points
00507           Input: a --> the first vector
00508                  b --> the second vector
00509          Output: 1 if they are equal, 0 otherwise
00510    Restrictions: This might be of limited value since comparing two floating
00511                  point numbers doesn't usually give the proper result.
00512      author Philippe Lavoie (24 January 1997)
00513     Modified by:
00514    * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
00515   template <class T, int D>
00516   inline int operator==(const HPoint_nD<T,D>& a, const HPoint_nD<T,D>& b){
00517     int r = 1 ;
00518     for(int i=D;i>=0;--i)
00519       r = r && a.data[i]==b.data[i] ;
00520     return r ;
00521   }
00522 
00523 
00524 
00525   /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
00526         Routine: norm2 --- the sum of the square of all the elements of a point
00527      The sum of the square of all the elements of a point or the 
00528                  length of the vector in homogenous space.
00529           Input: a --> the point
00530          Output: $a_x^2+a_y^2+a_z^2+a_w^2$
00531    Restrictions:
00532      author Philippe Lavoie (24 January 1997)
00533     Modified by:
00534    * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
00535   template <class T,int N>
00536   inline T norm2(const HPoint_nD<T,N>& a){
00537     double temp = 0 ;
00538     for(int i=N-1;i>=0;--i)
00539       temp += a.data[i]*a.data[i] ;
00540     return temp ;
00541   }
00542 
00543   template <class T,int N> 
00544   inline double norm(const HPoint_nD<T,N>& a) { return sqrt(norm2(a)); }
00545 
00546 
00547 
00548   /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
00549         Routine: angle --- the angle between two control points
00550      The angle between two control points
00551           Input: a --> the first vector
00552                  b --> the second vector
00553          Output: the angle in radians between the two points
00554    Restrictions:
00555      author Philippe Lavoie (24 January 1997)
00556     Modified by:
00557    * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
00558   template <class T, int D>
00559   inline T angle(const HPoint_nD<T,D>& a,const HPoint_nD<T,D>& b) {
00560     if(b==0 || a==0 )
00561       return 0 ;
00562     return dot(a,b)/norm(a)/norm(b) ;
00563   }
00564 
00565 
00566 
00567   /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
00568          Member: move --- translates a point in 4D by a point in 3D
00569      Translates a point in 4D by a point in 3D. This will only
00570                  modify the $x,y$ and $z$ component of the point in 4D.
00571                  It will not modify the $w$.
00572           Input: m --> the point in 3D
00573          Output: $a'_x = a_x+m_x, a'_y = a_y+m_y, a'_z = a_z+m_z, a'_w = a_w$
00574    Restrictions:
00575      author Philippe Lavoie (24 January 1997)
00576     Modified by:
00577    * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
00578   // is now a in the class definition
00579 
00580   //inline void HPoint_nD<T,D>::move(const Point_nD<T,D>& m); {
00581 
00582 
00583   /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
00584         Routine: operator<< --- the output operator to an ostream
00585      The output operator to an ostream
00586           Input:  os --> the ostream
00587                point --> the point to output
00588          Output: the ostream with the point 
00589    Restrictions:
00590      author Philippe Lavoie (24 January 1997)
00591     Modified by:
00592    * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
00593   template <class T>
00594   inline ostream& operator<<(ostream& os,const HPoint_nD<T,3>& point)
00595   {
00596     os << point.x() << " " << point.y() << " " << point.z() << " " << point.w() << " " ;
00597     return os;  
00598   }
00599 
00600   template <class T>
00601   inline ostream& operator<<(ostream& os,const HPoint_nD<T,2>& point)
00602   {
00603     os << point.x() << " " << point.y() << " " << point.w() << " " ;
00604     return os;  
00605   }
00606 
00607   /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
00608         Routine: distance3D --- the 3D distance between 2 control points
00609      The 3D distance between 2 control points
00610           Input: a --> the first point
00611                  b --> the second point
00612          Output: the distance
00613    Restrictions:
00614      author Philippe Lavoie (24 January 1997)
00615     Modified by:
00616    * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
00617   template <class T, int D>
00618   inline T distance3D(const HPoint_nD<T,D>& a, const HPoint_nD<T,D>& b) { 
00619     return norm(project(a)-project(b)) ;
00620   }
00621 
00622   /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
00623         Routine: operator>> --- the input operator from an istream
00624      Initialize a point in 4D from the input stream.
00625           Input:  os --> the input stream
00626                point <-- the point to initialize
00627          Output: the istream without the point
00628    Restrictions:
00629      author Philippe Lavoie (24 January 1997)
00630     Modified by:
00631    * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
00632   template <class T>
00633   inline istream& operator>>(istream& os, HPoint_nD<T,3>& point){
00634     float x,y,z,w ;
00635     os >> x >> y >> z >> w;
00636     point.data[0] = x ;
00637     point.data[1] = y ; 
00638     point.data[2] = z ;
00639     point.data[3] = w ;
00640     return os ;
00641   }
00642 
00643   template <class T>
00644   inline istream& operator>>(istream& os, HPoint_nD<T,2>& point){
00645     float x,y,w ;
00646     os >> x >> y >> w;
00647     point.data[0] = x ;
00648     point.data[1] = y ; 
00649     point.data[2] = w ;
00650     return os ;
00651   }
00652 
00653 
00654 
00655   /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
00656         Routine: project --- projects a point in 4D to a point in 3D
00657      Projects a point in 4D to a point in 3D
00658           Input: a --> the point in 4D
00659          Output: the point in 3D $(a_x/a_w, a_y/a_w, a_z/a_w)$
00660    Restrictions:
00661      author Philippe Lavoie (24 January 1997)
00662     Modified by:
00663    * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
00664 
00665   template <class T>
00666   inline Point_nD<T,3> project(const HPoint_nD<T,3>&a) {
00667     Point_nD<T,3> p(a.data[0]/a.data[3],a.data[1]/a.data[3],a.data[2]/a.data[3]) ;
00668     return p ;
00669   }
00670 
00671   template <class T>
00672   inline Point_nD<T,2> project(const HPoint_nD<T,2>&a) {
00673     Point_nD<T,2> p(a.x()/a.w(),a.y()/a.w()) ;
00674     return p ;
00675   }
00676 
00677 
00678   template <class T, int N>
00679   inline int operator<(const HPoint_nD<T,N>& a, const HPoint_nD<T,N>& b) { 
00680     return norm2(a)<norm2(b) ? 1 : 0 ;}
00681   template <class T, int N>
00682   inline int operator>(const HPoint_nD<T,N>& a, const HPoint_nD<T,N>& b) { 
00683     return norm2(a)>norm2(b) ? 1 : 0 ;}
00684   template <class T, int N>
00685   inline int operator<=(const HPoint_nD<T,N>& a, const HPoint_nD<T,N>& b) { 
00686     return norm2(a)<=norm2(b) ? 1 : 0 ;}
00687   template <class T, int N>
00688   inline int operator>=(const HPoint_nD<T,N>& a, const HPoint_nD<T,N>& b) { 
00689     return norm2(a)>=norm2(b) ? 1 : 0 ;}
00690 
00691 
00692   template <class T>
00693   inline HPoint_nD<T,2> minimum(HPoint_nD<T,2> a, HPoint_nD<T,2> b){
00694     HPoint_nD<T,2> m ;
00695     m.x() = minimum(a.x(),b.x()) ;
00696     m.y() = minimum(a.y(),b.y()) ;
00697     m.w() = minimum(a.w(),b.w()) ;
00698     return m ; 
00699   }
00700 
00701   template <class T>
00702   inline HPoint_nD<T,3> minimum(HPoint_nD<T,3> a, HPoint_nD<T,3> b){
00703     HPoint_nD<T,3> m ;
00704     m.x() = minimum(a.x(),b.x()) ;
00705     m.y() = minimum(a.y(),b.y()) ;
00706     m.z() = minimum(a.z(),b.z()) ;
00707     m.w() = minimum(a.w(),b.w()) ;
00708     return m ; 
00709   }
00710 
00711 
00712   template <class T>
00713   inline HPoint_nD<T,2> maximum(HPoint_nD<T,2> a, HPoint_nD<T,2> b){
00714     HPoint_nD<T,2> m ;
00715     m.x() = maximum(a.x(),b.x()) ;
00716     m.y() = maximum(a.y(),b.y()) ;
00717     m.w() = maximum(a.w(),b.w()) ;
00718     return m ; 
00719   }
00720 
00721   template <class T>
00722   inline HPoint_nD<T,3> maximum(HPoint_nD<T,3> a, HPoint_nD<T,3> b){
00723     HPoint_nD<T,3> m ;
00724     m.x() = maximum(a.x(),b.x()) ;
00725     m.y() = maximum(a.y(),b.y()) ;
00726     m.z() = maximum(a.z(),b.z()) ;
00727     m.w() = maximum(a.w(),b.w()) ;
00728     return m ; 
00729   }
00730 
00731 
00732   template <class T>
00733   inline HPoint_nD<T,3> minimumByRef(const HPoint_nD<T,3> &a, const HPoint_nD<T,3> &b){
00734     HPoint_nD<T,3> m ;
00735     m.x() = minimum(a.x(),b.x()) ;
00736     m.y() = minimum(a.y(),b.y()) ;
00737     m.z() = minimum(a.z(),b.z()) ;
00738     m.w() = minimum(a.w(),b.w()) ;
00739     return m ; 
00740   }
00741 
00742 
00743   template <class T>
00744   inline HPoint_nD<T,2> minimumByRef(const HPoint_nD<T,2> &a, const HPoint_nD<T,2> &b){
00745     HPoint_nD<T,2> m ;
00746     m.x() = minimum(a.x(),b.x()) ;
00747     m.y() = minimum(a.y(),b.y()) ;
00748     m.w() = minimum(a.w(),b.w()) ;
00749     return m ; 
00750   }
00751 
00752 
00753   template <class T>
00754   inline HPoint_nD<T,3> maximumByRef(const HPoint_nD<T,3> &a, const HPoint_nD<T,3> &b){
00755     HPoint_nD<T,3> m ;
00756     m.x() = maximum(a.x(),b.x()) ;
00757     m.y() = maximum(a.y(),b.y()) ;
00758     m.z() = maximum(a.z(),b.z()) ;
00759     m.w() = maximum(a.w(),b.w()) ;
00760     return m ; 
00761   }
00762 
00763 
00764 
00765   typedef HPoint_nD<float,3> HPoint3Df ;
00766   typedef HPoint_nD<double,3> HPoint3Dd ;
00767 
00768   typedef HPoint_nD<float,2> HPoint2Df ;
00769   typedef HPoint_nD<double,2> HPoint2Dd ;
00770 
00771 
00772 
00773 } // end namespace
00774 
00775 typedef PLib::HPoint_nD<float,3> PlHPoint3Df ;
00776 typedef PLib::HPoint_nD<double,3> PlHPoint3Dd ;
00777 
00778 typedef PLib::HPoint_nD<float,2> PlHPoint2Df ;
00779 typedef PLib::HPoint_nD<double,2> PlHPoint2Dd ;
00780 
00781 #endif