MeshKit  1.0
basic_math.hpp
Go to the documentation of this file.
00001 #ifndef BASIC_MATH_H
00002 #define BASIC_MATH_H
00003 
00004 #include <stdlib.h>
00005 #include <math.h>
00006 #include <assert.h>
00007 #include <iostream>
00008 #include <vector>
00009 using namespace std;
00010 
00011 #define ANGLE_IN_DEGREES  0
00012 #define ANGLE_IN_RADIANS  1
00013 
00014 template<class DataType, int n>
00015 class Array {
00016 public:
00017      DataType &operator[](int i) {
00018           return data[i];
00019      }
00020      const DataType &operator[](int i) const {
00021           return data[i];
00022      }
00023 private:
00024      DataType  data[n];
00025 };
00026 typedef Array<double,2> Point2D;
00027 typedef Array<double,3> Point3D;
00028 typedef Array<double,3> Array3D;
00029 typedef Array<double,4> Array4D;
00030 typedef Array<double,3> Vec3D;
00031 
00032 typedef Array<float,2> Point2F;
00033 typedef Array<float,3> Point3F;
00034 typedef Array<float,3> Array3F;
00035 typedef Array<float,4> Array4F;
00036 typedef Array<float,3> Vec3F;
00037 
00038 namespace Math {
00039 template<class T>
00040 T random_value(T minval, T maxval)
00041 {
00042 }
00043 
00044 template<>
00045 inline int random_value(int minval, int maxval)
00046 {
00047      return(minval + rand() % (int) (maxval - minval + 1));
00048 }
00049 
00050 template<>
00051 inline size_t random_value(size_t minval, size_t maxval)
00052 {
00053      return(size_t) (minval + rand() % (int) (maxval - minval + 1));
00054 }
00055 
00056 template<>
00057 inline double random_value(double minval, double maxval)
00058 {
00059      return minval + 0.5 * (1.0 + drand48())*(maxval - minval);
00060 }
00061 
00062 template<>
00063 inline float random_value(float minval, float maxval)
00064 {
00065      return minval + 0.5 * (1.0 + drand48())*(maxval - minval);
00066 }
00067 
00068 inline void create_vector( const Point3D &head, const Point3D &tail, Vec3D &xyz)
00069 {
00070      xyz[0] = head[0] - tail[0];
00071      xyz[1] = head[1] - tail[1];
00072      xyz[2] = head[2] - tail[2];
00073 }
00074 
00075 inline double length( const Point3D &A, const Point3D &B)
00076 {
00077      double dx = A[0] - B[0];
00078      double dy = A[1] - B[1];
00079      double dz = A[2] - B[2];
00080      return sqrt( dx*dx + dy*dy + dz*dz );
00081 }
00082 
00083 inline double length2( const Point3D &A, const Point3D &B)
00084 {
00085      double dx = A[0] - B[0];
00086      double dy = A[1] - B[1];
00087      double dz = A[2] - B[2];
00088      return dx*dx + dy*dy + dz*dz;
00089 }
00090 
00091 inline double magnitude( const Vec3D &A )
00092 {
00093      return sqrt( A[0]*A[0] + A[1]*A[1] + A[2]*A[2] );
00094 }
00095 
00096 inline double dot_product( const Vec3D &A, const Vec3D &B)
00097 {
00098      return A[0]*B[0] + A[1]*B[1] + A[2]*B[2];
00099 }
00100 
00101 inline void cross_product( const Vec3D &A, const Vec3D &B, Vec3D &C)
00102 {
00103      C[0] = A[1]*B[2] - A[2]*B[1];
00104      C[1] = A[2]*B[0] - A[0]*B[2];
00105      C[2] = A[0]*B[1] - A[1]*B[0];
00106 }
00107 
00108 inline double poly_area(const vector<Point2D> &p )
00109 {
00110      // Formula from wikipedia....
00111      double sum = 0.0;
00112      int nSize = p.size();
00113 
00114      assert( nSize >= 3);
00115 
00116      for( int i  = 0; i < nSize; i++)
00117           sum += p[i][0]*p[(i+1)%nSize][1] -p[(i+1)%nSize][0]*p[i][1];
00118      return 0.5*sum;
00119 }
00120 
00121 inline void poly_centroid( const vector<Point2D> &p, Point2D &c )
00122 {
00123      int nSize = p.size();
00124      assert( nSize >= 3);
00125 
00126      // Formula from wikipedia. Note that centroid does not depend on the
00127      // orientation of the polygon. The division by Area will take care
00128      // of correct value.
00129 
00130      // For convex bodies, centroid is always inside the region.
00131 
00132      double cx = 0.0;
00133      double cy = 0.0;
00134      double cf;
00135 
00136      for( int i  = 0; i < nSize; i++) {
00137           cf  = p[i][0]*p[(i+1)%nSize][1] - p[(i+1)%nSize][0]*p[i][1];
00138           cx +=  (p[i][0]+p[(i+1)%nSize][0])*cf;
00139           cy +=  (p[i][1]+p[(i+1)%nSize][1])*cf;
00140      }
00141 
00142      double A = poly_area( p );
00143 
00144      c[0] = cx/(6.0*A);
00145      c[1] = cy/(6.0*A);
00146 }
00147 
00148 inline void normal( const Point3D &p0, const Point3D &p1, const Point3D &p2,
00149                     Vec3D &normal)
00150 {
00151      Vec3D p1p0, p2p0;
00152      Math::create_vector( p2, p0, p2p0);
00153      Math::create_vector( p1, p0, p1p0);
00154      Math::cross_product( p2p0, p1p0, normal);
00155 
00156      double mag = Math::magnitude( normal );
00157      normal[0] /= mag;
00158      normal[1] /= mag;
00159      normal[2] /= mag;
00160 }
00161 
00162 inline Vec3D unit_vector( const Point3D &head, const Point3D &tail)
00163 {
00164      Vec3D uvec;
00165      create_vector( head, tail, uvec);
00166      double dl  = magnitude(uvec);
00167 
00168      uvec[0]  /=  dl;
00169      uvec[1]  /=  dl;
00170      uvec[2]  /=  dl;
00171 
00172      return uvec;
00173 }
00174 
00176 
00177 inline double getVectorAngle( const Vec3D &A, const Vec3D &B, int measure)
00178 {
00179      double AB = dot_product(A,B);
00180      double Am = magnitude(A);
00181      double Bm = magnitude(B);
00182 
00183      if( Am < 1.0E-15 || Bm < 1.0E-15) return 0.0;
00184 
00185      double x = AB/(Am*Bm);
00186 
00187      if( x > 1.0) x = 1.0;
00188      if( x < -1.0) x = -1.0;
00189 
00190      if( measure == ANGLE_IN_DEGREES ) return 180*acos(x)/M_PI;
00191      return acos(x);
00192 }
00193 
00195 
00196 template<class T>
00197 inline T max_value( const T &a, const T &b, const T &c)
00198 {
00199      return max(a,max(b, c));
00200 }
00201 
00202 template<class T>
00203 inline T min_value( const T &a, const T &b, const T &c)
00204 {
00205      return min(a,min(b, c));
00206 }
00207 
00208 template <class T, size_t n>
00209 inline double getAngle(const Array<T, n> &VecA, const Array<T, n> &VecB,
00210                        int unit_measure)
00211 {
00212      double Abar, Bbar, theta;
00213      Abar = magnitude(VecA);
00214      Bbar = magnitude(VecB);
00215 
00216      if (Abar < 1.0E-10 || Bbar < 1.0E-10) {
00217           cout << " Warning: Error in Angle calculation " << endl;
00218           cout << " Magnitude of Vector A is " << Abar << endl;
00219           cout << " Magnitude of Vector B is " << Bbar << endl;
00220           return 0.0;
00221      }
00222 
00223      double value = dot_product(VecA, VecB) / (Abar * Bbar);
00224 
00225      if (value > +1.0) value = +1.0;
00226      if (value < -1.0) value = -1.0;
00227 
00228      theta = acos(value);
00229 
00230      if (unit_measure == ANGLE_IN_DEGREES) theta *= (180.0 / M_PI);
00231 
00232      return theta;
00233 }
00234 
00236 
00237 template <class T, size_t n>
00238 inline T getAngle(const Array<T, n> &pa, const Array<T, n> &pb,
00239                   const Array<T, n> &pc, int unit_measure = 0)
00240 {
00241      Array<T, n> VecA = create_vector(pb, pa);
00242      Array<T, n> VecB = create_vector(pc, pa);
00243      return getAngle(VecA, VecB, unit_measure);
00244 }
00245 
00247 
00248 inline double getTriAngle(const Point3D &pa, const Point3D &pb, const Point3D &pc)
00249 {
00250      double a2   =  length2( pb, pc );
00251      double b2   =  length2( pc, pa );
00252      double c2   =  length2( pa, pb );
00253      double cosA =  (b2 + c2 - a2)/(2*sqrt(b2*c2) );
00254 
00255      if( cosA >  1.0) cosA =  1.0;
00256      if( cosA < -1.0) cosA = -1.0;
00257      return 180*acos(cosA)/M_PI;
00258 
00259 }
00261 
00262 inline void getTriAngles(const Point3D &pa, const Point3D &pb, const Point3D &pc, Point3D &angles)
00263 {
00264      double a2   =  length2( pb, pc );
00265      double b2   =  length2( pc, pa );
00266      double c2   =  length2( pa, pb );
00267      double cosA =  (b2 + c2 - a2)/(2*sqrt(b2*c2) );
00268      double cosB =  (a2 + c2 - b2)/(2*sqrt(a2*c2) );
00269      double cosC =  (a2 + b2 - c2)/(2*sqrt(a2*b2) );
00270 
00271      if( cosA >  1.0) cosA =  1.0;
00272      if( cosA < -1.0) cosA = -1.0;
00273      angles[0] = 180*acos(cosA)/M_PI;
00274 
00275      if( cosB >  1.0) cosB =  1.0;
00276      if( cosB < -1.0) cosB = -1.0;
00277      angles[1] = 180*acos(cosB)/M_PI;
00278 
00279      if( cosC >  1.0) cosC =  1.0;
00280      if( cosC < -1.0) cosC = -1.0;
00281      angles[2] = 180*acos(cosC)/M_PI;
00282 }
00283 
00285 }
00286 
00287 #include "FastArea.hpp"
00288 
00289 #endif
00290 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines