MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 #include <math.h> 00002 #include <vector> 00003 #include <stdexcept> 00004 #include <assert.h> 00005 00006 class geomObject 00007 { 00008 public: 00009 virtual ~geomObject() {} 00010 virtual void project_points2geom( int dim, double* oldcoords, double* newcoords, double* derivs ) const = 0; 00011 virtual double compute_projecterror( int dim, double* oldcoords ) const 00012 { 00013 double* newcoords = new double[dim](); 00014 project_points2geom( dim, oldcoords, newcoords, NULL ); 00015 double ans = 0; 00016 00017 for( int i = 0; i < dim; ++i ) 00018 { 00019 ans += ( newcoords[i] - oldcoords[i] ) * ( newcoords[i] - oldcoords[i] ); 00020 } 00021 00022 delete[] newcoords; 00023 return sqrt( ans ); 00024 } 00025 virtual void compute_projecterror( int dim, int nverts, double* oldcoords, double& l1err, double& l2err, 00026 double& linferr ) const 00027 { 00028 l1err = l2err = linferr = 0; 00029 00030 for( int i = 0; i < nverts; ++i ) 00031 { 00032 double err = compute_projecterror( dim, oldcoords + i * dim ); 00033 l1err += err; 00034 l2err += err * err; 00035 linferr = std::max( linferr, err ); 00036 } 00037 00038 l1err /= nverts; 00039 l2err = sqrt( l2err / nverts ); 00040 } 00041 inline double Twonorm( int dim, double* vec ) const 00042 { 00043 double ans = 0; 00044 00045 for( int i = 0; i < dim; ++i ) 00046 { 00047 ans += vec[i] * vec[i]; 00048 } 00049 00050 return sqrt( ans ); 00051 } 00052 }; 00053 00054 class sphere : public geomObject 00055 { 00056 public: 00057 sphere( double x = 0, double y = 0, double z = 0, double r = 1 ) 00058 : centerx( x ), centery( y ), centerz( z ), radius( r ) 00059 { 00060 assert( r > 0 ); 00061 } 00062 virtual ~sphere() {} 00063 void project_points2geom( int dim, double* oldcoords, double* newcoords, double* derivs ) const 00064 { 00065 if( oldcoords == NULL || newcoords == NULL ) { throw std::invalid_argument( "NULL pointer" ); } 00066 00067 std::vector< double > direction; 00068 direction.push_back( oldcoords[0] - centerx ); 00069 direction.push_back( oldcoords[1] - centery ); 00070 00071 if( dim == 3 ) { direction.push_back( oldcoords[2] - centerz ); } 00072 00073 double len = geomObject::Twonorm( dim, &( direction[0] ) ); 00074 assert( len > 0 ); 00075 00076 for( int i = 0; i < dim; ++i ) 00077 { 00078 direction[i] /= len; 00079 } 00080 00081 newcoords[0] = centerx + direction[0] * radius; 00082 newcoords[1] = centery + direction[1] * radius; 00083 00084 if( dim == 2 ) 00085 { 00086 if( derivs ) 00087 { 00088 derivs[0] = -direction[1]; 00089 derivs[1] = direction[0]; 00090 } 00091 00092 return; 00093 } 00094 else if( dim == 3 ) 00095 { 00096 newcoords[2] = centerz + direction[2] * radius; 00097 00098 if( derivs ) 00099 { 00100 derivs[0] = direction[0]; 00101 derivs[1] = direction[1]; 00102 derivs[2] = direction[2]; 00103 } 00104 } 00105 else 00106 { 00107 throw std::invalid_argument( "dim must be 2 or 3" ); 00108 } 00109 } 00110 00111 private: 00112 double centerx, centery, centerz, radius; 00113 }; 00114 00115 class torus : public geomObject 00116 { 00117 public: 00118 torus( double x = 0, double y = 0, double z = 0, double aa = 0.3, double cc = 1.0 ) 00119 : centerx( x ), centery( y ), centerz( z ), a( aa ), c( cc ) 00120 { 00121 assert( aa > 0 && cc > 0 ); 00122 } 00123 virtual ~torus() {} 00124 void project_points2geom( int dim, double* oldcoords, double* newcoords, double* derivs ) const 00125 { 00126 assert( dim == 3 && oldcoords && newcoords ); 00127 double transfer[3] = { oldcoords[0] - centerx, oldcoords[1] - centery, oldcoords[2] - centerz }; 00128 double twodnrm = geomObject::Twonorm( 2, transfer ); 00129 double tubecenter[3] = { c * transfer[0] / twodnrm + centerx, c * transfer[1] / twodnrm + centery, centerz }; 00130 double direction[3] = { oldcoords[0] - tubecenter[0], oldcoords[1] - tubecenter[1], 00131 oldcoords[2] - tubecenter[2] }; 00132 double len = geomObject::Twonorm( 3, direction ); 00133 assert( len > 0 ); 00134 direction[0] /= len; 00135 direction[1] /= len; 00136 direction[2] /= len; 00137 00138 if( derivs ) 00139 { 00140 derivs[0] = direction[0]; 00141 derivs[1] = direction[1]; 00142 derivs[2] = direction[2]; 00143 } 00144 00145 newcoords[0] = tubecenter[0] + a * direction[0]; 00146 newcoords[1] = tubecenter[1] + a * direction[1]; 00147 newcoords[2] = tubecenter[2] + a * direction[2]; 00148 } 00149 00150 private: 00151 double centerx, centery, centerz, a, c; 00152 };