00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027 #ifndef WFMATH_POINT_FUNCS_H
00028 #define WFMATH_POINT_FUNCS_H
00029
00030 #include <wfmath/const.h>
00031 #include <wfmath/vector.h>
00032 #include <wfmath/point.h>
00033
00034 namespace WFMath {
00035
00036 template<const int dim>
00037 Point<dim>::Point(const Point<dim>& p) : m_valid(p.m_valid)
00038 {
00039 for(int i = 0; i < dim; ++i)
00040 m_elem[i] = p.m_elem[i];
00041 }
00042
00043 template<const int dim>
00044 Point<dim>& Point<dim>::setToOrigin()
00045 {
00046 for(int i = 0; i < dim; ++i)
00047 m_elem[i] = 0;
00048
00049 m_valid = true;
00050
00051 return *this;
00052 }
00053
00054 template<const int dim>
00055 bool Point<dim>::isEqualTo(const Point<dim> &p, double epsilon) const
00056 {
00057 CoordType delta = (CoordType) _ScaleEpsilon(m_elem, p.m_elem, dim, epsilon);
00058
00059 for(int i = 0; i < dim; ++i)
00060 if(fabs(m_elem[i] - p.m_elem[i]) > delta)
00061 return false;
00062
00063 return true;
00064 }
00065
00066 template<const int dim>
00067 Vector<dim> operator-(const Point<dim>& c1, const Point<dim>& c2)
00068 {
00069 Vector<dim> out;
00070
00071 for(int i = 0; i < dim; ++i)
00072 out.m_elem[i] = c1.m_elem[i] - c2.m_elem[i];
00073
00074 out.m_valid = c1.m_valid && c2.m_valid;
00075
00076 return out;
00077 }
00078
00079 template<const int dim>
00080 Point<dim> operator+(const Point<dim>& c, const Vector<dim>& v)
00081 {
00082 Point<dim> out;
00083
00084 for(int i = 0; i < dim; ++i)
00085 out.m_elem[i] = c.m_elem[i] + v.m_elem[i];
00086
00087 out.m_valid = c.m_valid && v.m_valid;
00088
00089 return out;
00090 }
00091
00092 template<const int dim>
00093 Point<dim> operator-(const Point<dim>& c, const Vector<dim>& v)
00094 {
00095 Point<dim> out;
00096
00097 for(int i = 0; i < dim; ++i)
00098 out.m_elem[i] = c.m_elem[i] - v.m_elem[i];
00099
00100 out.m_valid = c.m_valid && v.m_valid;
00101
00102 return out;
00103 }
00104
00105 template<const int dim>
00106 Point<dim> operator+(const Vector<dim>& v, const Point<dim>& c)
00107 {
00108 Point<dim> out;
00109
00110 for(int i = 0; i < dim; ++i)
00111 out.m_elem[i] = c.m_elem[i] + v.m_elem[i];
00112
00113 out.m_valid = c.m_valid && v.m_valid;
00114
00115 return out;
00116 }
00117
00118 template<const int dim>
00119 Point<dim>& Point<dim>::operator=(const Point<dim>& rhs)
00120 {
00121
00122
00123 if (this == &rhs)
00124 return *this;
00125
00126 for(int i = 0; i < dim; ++i)
00127 m_elem[i] = rhs.m_elem[i];
00128
00129 m_valid = rhs.m_valid;
00130
00131 return *this;
00132 }
00133
00134 template<const int dim>
00135 Point<dim>& operator+=(Point<dim>& p, const Vector<dim> &rhs)
00136 {
00137 for(int i = 0; i < dim; ++i)
00138 p.m_elem[i] += rhs.m_elem[i];
00139
00140 p.m_valid = p.m_valid && rhs.m_valid;
00141
00142 return p;
00143 }
00144
00145 template<const int dim>
00146 Point<dim>& operator-=(Point<dim>& p, const Vector<dim> &rhs)
00147 {
00148 for(int i = 0; i < dim; ++i)
00149 p.m_elem[i] -= rhs.m_elem[i];
00150
00151 p.m_valid = p.m_valid && rhs.m_valid;
00152
00153 return p;
00154 }
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164 template<const int dim>
00165 bool Point<dim>::operator< (const Point<dim>& rhs) const
00166 {
00167 if(operator==(rhs))
00168 return false;
00169
00170 for(int i = 0; i < dim; ++i)
00171 if(m_elem[i] != rhs.m_elem[i])
00172 return m_elem[i] < rhs.m_elem[i];
00173
00174 assert(false);
00175 return false;
00176 }
00177
00178 template<const int dim>
00179 CoordType SquaredDistance(const Point<dim>& p1, const Point<dim>& p2)
00180 {
00181 CoordType ans = 0;
00182
00183 for(int i = 0; i < dim; ++i) {
00184 CoordType diff = p1.m_elem[i] - p2.m_elem[i];
00185 ans += diff * diff;
00186 }
00187
00188 return (fabs(ans) >= _ScaleEpsilon(p1.m_elem, p2.m_elem, dim)) ? ans : 0;
00189 }
00190
00191 #ifndef WFMATH_NO_TEMPLATES_AS_TEMPLATE_PARAMETERS
00192 template<const int dim, template<class> class container,
00193 template<class> class container2>
00194 Point<dim> Barycenter(const container<Point<dim> >& c,
00195 const container2<CoordType>& weights)
00196 {
00197
00198
00199 typename container<Point<dim> >::const_iterator c_i = c.begin(), c_end = c.end();
00200 typename container2<CoordType>::const_iterator w_i = weights.begin(),
00201 w_end = weights.end();
00202
00203 assert("nonempty list of points" && c_i != c_end);
00204 assert("nonempty list of weights" && w_i != w_end);
00205
00206 bool valid = c_i->isValid();
00207
00208 CoordType tot_weight = *w_i, max_weight = fabs(*w_i);
00209 Point<dim> out;
00210 for(int j = 0; j < dim; ++j)
00211 out[j] = (*c_i)[j] * *w_i;
00212
00213 while(++c_i != c_end && ++w_i != w_end) {
00214 tot_weight += *w_i;
00215 CoordType val = fabs(*w_i);
00216 if(val > max_weight)
00217 max_weight = val;
00218 if(!c_i->isValid())
00219 valid = false;
00220 for(int j = 0; j < dim; ++j)
00221 out[j] += (*c_i)[j] * *w_i;
00222 }
00223
00224
00225 assert("sum of weights must be nonzero" && max_weight > 0
00226 && fabs(tot_weight) > max_weight * WFMATH_EPSILON);
00227
00228 for(int j = 0; j < dim; ++j)
00229 out[j] /= tot_weight;
00230
00231 out.setValid(valid);
00232
00233 return out;
00234 }
00235
00236 template<const int dim, template<class> class container>
00237 Point<dim> Barycenter(const container<Point<dim> >& c)
00238 {
00239
00240
00241 typename container<Point<dim> >::const_iterator i = c.begin(), end = c.end();
00242
00243 assert("nonempty list of points" && i != end);
00244
00245 Point<dim> out = *i;
00246 int num_points = 1;
00247
00248 bool valid = i->isValid();
00249
00250 while(++i != end) {
00251 ++num_points;
00252 if(!i->isValid())
00253 valid = false;
00254 for(int j = 0; j < dim; ++j)
00255 out[j] += (*i)[j];
00256 }
00257
00258 for(int j = 0; j < dim; ++j)
00259 out[j] /= num_points;
00260
00261 out.setValid(valid);
00262
00263 return out;
00264 }
00265 #endif
00266
00267 template<const int dim>
00268 Point<dim> Midpoint(const Point<dim>& p1, const Point<dim>& p2, CoordType dist)
00269 {
00270 Point<dim> out;
00271 CoordType conj_dist = 1 - dist;
00272
00273 for(int i = 0; i < dim; ++i)
00274 out.m_elem[i] = p1.m_elem[i] * conj_dist + p2.m_elem[i] * dist;
00275
00276 out.m_valid = p1.m_valid && p2.m_valid;
00277
00278 return out;
00279 }
00280
00281 template<> inline Point<2>::Point(CoordType x, CoordType y) : m_valid(true)
00282 {
00283 m_elem[0] = x;
00284 m_elem[1] = y;
00285 }
00286
00287 template<> inline Point<3>::Point(CoordType x, CoordType y, CoordType z) : m_valid(true)
00288 {
00289 m_elem[0] = x;
00290 m_elem[1] = y;
00291 m_elem[2] = z;
00292 }
00293
00294 #ifndef WFMATH_NO_CLASS_FUNCTION_SPECIALIZATION
00295 template<> Point<2>& Point<2>::polar(CoordType r, CoordType theta);
00296 template<> void Point<2>::asPolar(CoordType& r, CoordType& theta) const;
00297
00298 template<> Point<3>& Point<3>::polar(CoordType r, CoordType theta,
00299 CoordType z);
00300 template<> void Point<3>::asPolar(CoordType& r, CoordType& theta,
00301 CoordType& z) const;
00302 template<> Point<3>& Point<3>::spherical(CoordType r, CoordType theta,
00303 CoordType phi);
00304 template<> void Point<3>::asSpherical(CoordType& r, CoordType& theta,
00305 CoordType& phi) const;
00306 #else
00307 void _NCFS_Point2_polar(CoordType *m_elem, CoordType r, CoordType theta);
00308 void _NCFS_Point2_asPolar(CoordType *m_elem, CoordType& r, CoordType& theta);
00309
00310 void _NCFS_Point3_polar(CoordType *m_elem, CoordType r, CoordType theta,
00311 CoordType z);
00312 void _NCFS_Point3_asPolar(CoordType *m_elem, CoordType& r, CoordType& theta,
00313 CoordType& z);
00314 void _NCFS_Point3_spherical(CoordType *m_elem, CoordType r, CoordType theta,
00315 CoordType phi);
00316 void _NCFS_Point3_asSpherical(CoordType *m_elem, CoordType& r, CoordType& theta,
00317 CoordType& phi);
00318
00319 template<>
00320 inline Point<2>& Point<2>::polar(CoordType r, CoordType theta)
00321 {
00322 _NCFS_Point2_polar((CoordType*) m_elem, r, theta);
00323 m_valid = true;
00324 return *this;
00325 }
00326
00327 template<>
00328 inline void Point<2>::asPolar(CoordType& r, CoordType& theta) const
00329 {
00330 _NCFS_Point2_asPolar((CoordType*) m_elem, r, theta);
00331 }
00332
00333 template<>
00334 inline Point<3>& Point<3>::polar(CoordType r, CoordType theta, CoordType z)
00335 {
00336 _NCFS_Point3_polar((CoordType*) m_elem, r, theta, z);
00337 m_valid = true;
00338 return *this;
00339 }
00340
00341 template<>
00342 inline void Point<3>::asPolar(CoordType& r, CoordType& theta, CoordType& z) const
00343 {
00344 _NCFS_Point3_asPolar((CoordType*) m_elem, r, theta, z);
00345 }
00346
00347 template<>
00348 inline Point<3>& Point<3>::spherical(CoordType r, CoordType theta, CoordType phi)
00349 {
00350 _NCFS_Point3_spherical((CoordType*) m_elem, r, theta, phi);
00351 m_valid = true;
00352 return *this;
00353 }
00354
00355 template<>
00356 inline void Point<3>::asSpherical(CoordType& r, CoordType& theta, CoordType& phi) const
00357 {
00358 _NCFS_Point3_asSpherical((CoordType*) m_elem, r, theta, phi);
00359 }
00360 #endif
00361
00362 }
00363
00364 #endif // WFMATH_POINT_FUNCS_H