//=========================================================================== // GoTools Core - SINTEF Geometry Tools Core library, version 2.0.1 // // Copyright (C) 2000-2007, 2010 SINTEF ICT, Applied Mathematics, Norway. // // This program is free software; you can redistribute it and/or // modify it under the terms of the GNU General Public License // as published by the Free Software Foundation version 2 of the License. // // This program is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU General Public License for more details. // // You should have received a copy of the GNU General Public License // along with this program; if not, write to the Free Software // Foundation, Inc., // 59 Temple Place - Suite 330, // Boston, MA 02111-1307, USA. // // Contact information: E-mail: tor.dokken@sintef.no // SINTEF ICT, Department of Applied Mathematics, // P.O. Box 124 Blindern, // 0314 Oslo, Norway. // // Other licenses are also available for this software, notably licenses // for: // - Building commercial software. // - Building software whose source code you wish to keep private. //=========================================================================== 00013 #ifndef _POINT_H 00014 #define _POINT_H 00015 00016 00017 #include "GoTools/utils/Array.h" 00018 #include "GoTools/utils/errormacros.h" 00019 #include <algorithm> 00020 #include <math.h> 00021 #include "GoTools/utils/config.h" 00022 00023 00024 namespace Go { 00025 00026 00035 class GO_API Point 00036 { 00037 private: 00038 double* pstart_; 00039 int n_; 00040 bool owns_; 00041 00042 public: 00047 Point() 00048 : pstart_(0), n_(0), owns_(true) 00049 {} 00053 explicit Point(int dim) 00054 : pstart_(new double[dim]), n_(dim), owns_(true) 00055 {} 00058 Point(double x, double y) 00059 : pstart_(new double[2]), n_(2), owns_(true) 00060 { 00061 pstart_[0] = x; 00062 pstart_[1] = y; 00063 } 00066 Point(double x, double y, double z) 00067 : pstart_(new double[3]), n_(3), owns_(true) 00068 { 00069 pstart_[0] = x; 00070 pstart_[1] = y; 00071 pstart_[2] = z; 00072 } 00073 00076 // @@@ 'typename T' must be present in order to compile on MS 00077 // v6.0. (Using just 'int Dim' and 'Array<double, Dim>' will cause 00078 // "internal compiler error" !) @jbt 00079 template <typename T, int Dim> 00080 explicit Point(const Array<T, Dim>& v) 00081 : pstart_(0), n_(Dim), owns_(true) 00082 { 00083 pstart_ = new double[n_]; 00084 #ifndef _MSC_VER // Getting rid of warning C4996 on Windows 00085 std::copy(v.begin(), v.end(), pstart_); 00086 #else 00087 stdext::unchecked_copy(v.begin(), v.end(), pstart_); 00088 #endif // _MSC_VER 00089 } 00090 00091 // #ifndef _MSC_VER 00093 template <typename RandomAccessIterator> 00094 Point(RandomAccessIterator first, RandomAccessIterator last) 00095 : pstart_(0), n_(last - first), owns_(true) 00096 { 00097 pstart_ = new double[n_]; 00098 #ifndef _MSC_VER // Getting rid of warning C4996 on Windows 00099 std::copy(first, last, pstart_); 00100 #else 00101 stdext::unchecked_copy(first, last, pstart_); 00102 #endif // _MSC_VER 00103 } 00104 // #endif 00105 // #ifdef _MSC_VER 00106 // #if _MSC_VER > 1200 00107 // template <typename ForwardIterator> 00108 // Point(ForwardIterator first, ForwardIterator last) 00109 // : pstart_(0), n_(last - first), owns_(true) 00110 // { 00111 // pstart_ = new double[n_]; 00112 // std::copy(first, last, pstart_); 00113 // } 00114 // #else // Old version 00115 // Point(const double* first, const double* last) 00116 // : pstart_(0), n_(last - first), owns_(true) 00117 // { 00118 // pstart_ = new double[n_]; 00119 // std::copy(first, last, pstart_); 00120 // } 00121 // Point(const float* first, const float* last) 00122 // : pstart_(0), n_(last - first), owns_(true) 00123 // { 00124 // pstart_ = new double[n_]; 00125 // std::copy(first, last, pstart_); 00126 // } 00127 // #endif // _MSC_VER > 1200 00128 // #endif 00129 00134 Point(double* begin, double* end, bool own) 00135 : pstart_(0), n_(end-begin), owns_(own) 00136 { 00137 if (owns_) { 00138 pstart_ = new double[n_]; 00139 #ifndef _MSC_VER // Getting rid of warning C4996 on Windows 00140 std::copy(begin, end, pstart_); 00141 #else 00142 stdext::unchecked_copy(begin, end, pstart_); 00143 #endif // _MSC_VER 00144 } else { 00145 pstart_ = begin; 00146 } 00147 } 00148 00150 Point(const Point& v) 00151 : pstart_(0), n_(v.n_), owns_(true) 00152 { 00153 pstart_ = new double[n_]; 00154 #ifndef _MSC_VER // Getting rid of warning C4996 on Windows 00155 std::copy(v.pstart_, v.pstart_ + n_, pstart_); 00156 #else 00157 stdext::unchecked_copy(v.pstart_, v.pstart_ + n_, pstart_); 00158 #endif // _MSC_VER 00159 } 00160 00162 Point& operator = (const Point &v) 00163 { 00164 Point temp(v); 00165 swap(temp); 00166 return *this; 00167 } 00168 00170 ~Point() 00171 { 00172 if (owns_) delete [] pstart_; 00173 } 00174 00176 void swap(Point& other) 00177 { 00178 std::swap(pstart_, other.pstart_); 00179 std::swap(n_, other.n_); 00180 std::swap(owns_, other.owns_); 00181 } 00182 00186 void read(std::istream& is); 00187 00192 void write(std::ostream& os) const; 00193 00194 00196 const double& operator [] (int i) const { return pstart_[i]; } 00198 double& operator [] (int i) { return pstart_[i]; } 00200 const double* begin() const { return pstart_; } 00202 double* begin() { return pstart_; } 00204 const double* end() const { return pstart_ + n_; } 00206 double* end() { return pstart_ + n_; } 00208 int size() const { return n_; } 00210 int dimension() const { return n_; } 00211 00213 void resize(int d) 00214 { 00215 if (n_ < d) { 00216 Point temp(d); 00217 swap(temp); 00218 } else { 00219 n_ = d; 00220 } 00221 } 00222 00224 void setValue(double x, double y) 00225 { 00226 resize(2); 00227 pstart_[0] = x; 00228 pstart_[1] = y; 00229 } 00231 void setValue(double x, double y, double z) 00232 { 00233 resize(3); 00234 pstart_[0] = x; 00235 pstart_[1] = y; 00236 pstart_[2] = z; 00237 } 00241 void setValue(const double* array) 00242 { 00243 for (int i = 0; i < n_; ++i) { 00244 pstart_[i] = array[i]; 00245 } 00246 } 00250 void setValue(double val) 00251 { 00252 for (int i = 0; i < n_; ++i) { 00253 pstart_[i] = val; 00254 } 00255 } 00256 00257 void resetValue(int idx, double val) 00258 { 00259 DEBUG_ERROR_IF(idx < 0 || idx >= n_, 00260 "Dimension mismatch."); 00261 pstart_[idx] = val; 00262 } 00263 00266 double length2() const 00267 { 00268 double l2 = 0; 00269 for (int i = 0; i < n_; ++i) 00270 l2 += pstart_[i]*pstart_[i]; 00271 return l2; 00272 } 00275 double length() const 00276 { 00277 return sqrt(length2()); 00278 } 00281 double lengthInf() const 00282 { 00283 double linf = 0; 00284 for (int i = 0; i < n_; ++i) { 00285 linf = std::max(linf, fabs(pstart_[i])); 00286 } 00287 return linf; 00288 } 00292 double dist2(const Point &v) const 00293 { 00294 double l2 = 0; 00295 double d; 00296 for (int i = 0; i < n_; ++i) { 00297 d = pstart_[i] - v.pstart_[i]; 00298 l2 += d*d; 00299 } 00300 return l2; 00301 } 00305 double dist(const Point &v) const 00306 { 00307 return sqrt(dist2(v)); 00308 } 00312 double distInf(const Point &v) const 00313 { 00314 double linf = 0; 00315 double d; 00316 for (int i = 0; i < n_; ++i) { 00317 d = pstart_[i] - v.pstart_[i]; 00318 linf = std::max(linf, fabs(d)); 00319 } 00320 return linf; 00321 } 00324 void normalize() 00325 { 00326 double tl = length(); 00327 DEBUG_ERROR_IF(tl == 0.0, "Cannot normalize vector of zero length"); 00328 (*this) /= tl; 00329 } 00330 00331 double normalize_checked() 00332 { 00333 double tl = length(); 00334 if (tl<1.0e-12) 00335 return 0; 00336 00337 (*this) /= tl; 00338 return tl; 00339 } 00340 00342 Point operator + (const Point &v) const 00343 { 00344 Point res(*this); 00345 res += v; 00346 return res; 00347 } 00349 void operator +=(const Point &v) 00350 { 00351 DEBUG_ERROR_IF(n_!=v.n_, 00352 "Dimension mismatch."); 00353 for (int i = 0; i < n_; ++i) 00354 pstart_[i] += v.pstart_[i]; 00355 } 00356 00358 Point operator - (const Point &v) const 00359 { 00360 Point res(*this); 00361 res -= v; 00362 return res; 00363 } 00365 void operator -=(const Point &v) 00366 { 00367 DEBUG_ERROR_IF(n_!=v.n_, 00368 "Dimension mismatch."); 00369 for (int i = 0; i < n_; ++i) 00370 pstart_[i] -= v.pstart_[i]; 00371 } 00373 Point operator * (double d) const 00374 { 00375 Point res(*this); 00376 res *= d; 00377 return res; 00378 } 00380 void operator *=(double d) 00381 { 00382 for (int i = 0; i < n_; ++i) 00383 pstart_[i] *= d; 00384 } 00385 00387 Point operator / (double d) const 00388 { 00389 Point res(*this); 00390 res /= d; 00391 return res; 00392 } 00394 void operator /=(double d) 00395 { 00396 for (int i = 0; i < n_; ++i) 00397 pstart_[i] /= d; 00398 } 00399 00401 Point operator - () const 00402 { 00403 Point res(*this); 00404 for (int i = 0; i < n_; ++i) 00405 res.pstart_[i] = - pstart_[i]; 00406 return res; 00407 } 00408 00411 double operator * (const Point &v) const 00412 { 00413 DEBUG_ERROR_IF(n_!=v.n_, 00414 "Dimension mismatch."); 00415 double res = 0; 00416 for (int i = 0; i < n_; ++i) 00417 res += pstart_[i]*v.pstart_[i]; 00418 return res; 00419 } 00420 00423 Point operator % (const Point &v) const 00424 { 00425 DEBUG_ERROR_IF(n_ != 3 || v.n_ != 3, 00426 "Dimension mismatch."); 00427 return Point(pstart_[1]*v.pstart_[2] - pstart_[2]*v.pstart_[1], 00428 pstart_[2]*v.pstart_[0] - pstart_[0]*v.pstart_[2], 00429 pstart_[0]*v.pstart_[1] - pstart_[1]*v.pstart_[0]); 00430 } 00431 00434 Point cross(const Point &v) const 00435 { 00436 return operator%(v); 00437 } 00438 00440 void setToCrossProd(const Point &u, const Point &v) 00441 { 00442 DEBUG_ERROR_IF(u.n_!=v.n_, 00443 "Dimension mismatch."); 00444 DEBUG_ERROR_IF(u.n_!=3, 00445 "Dimension must be 3."); 00446 00447 bool have_already = owns_ && (n_ >= v.n_); 00448 if (!have_already) { 00449 Point temp(3); 00450 swap(temp); 00451 } else { 00452 n_ = 3; 00453 } 00454 pstart_[0]=u.pstart_[1]*v.pstart_[2] - u.pstart_[2]*v.pstart_[1]; 00455 pstart_[1]=u.pstart_[2]*v.pstart_[0] - u.pstart_[0]*v.pstart_[2]; 00456 pstart_[2]=u.pstart_[0]*v.pstart_[1] - u.pstart_[1]*v.pstart_[0]; 00457 } 00458 00461 double cosAngle(const Point& v) const 00462 { 00463 double tl1 = length(); 00464 double tl2 = v.length(); 00465 //(DEBUG_ERROR_IF(tl1*tl2 == 0.0, "Vector of zero length"); 00466 if (tl1*tl2 == 0.0) 00467 { 00468 //MESSAGE("Vector of zero length in angle compuation"); 00469 return 0.0; 00470 } 00471 double res = ((*this)*v)/(tl1*tl2); 00472 res = std::min(1.0, std::max(-1.0, res)); 00473 return res; 00474 } 00475 00477 double angle(const Point& v) const 00478 { 00479 return acos(cosAngle(v)); 00480 } 00481 00484 double angle_smallest(const Point& v) const 00485 { 00486 double tl1 = length(); 00487 double tl2 = v.length(); 00488 //DEBUG_ERROR_IF(tl1*tl2 == 0.0, "Vector of zero length"); 00489 double tcos = fabs(((*this)*v)/(tl1*tl2)); 00490 tcos = std::min(1.0, tcos); 00491 return acos(tcos); 00492 } 00493 00494 }; 00495 00497 inline Point operator * (double d, const Point& p) 00498 { return p*d; } 00499 00501 inline std::istream& operator>>(std::istream& is, Go::Point& v) 00502 { v.read(is); return is; } 00503 00505 inline std::ostream& operator<<(std::ostream& os, const Go::Point& v) 00506 { v.write(os); return os; } 00507 00508 // comparison operator 00509 inline bool operator<(const Point& p1, const Point& p2) 00510 { 00511 const int dim = p1.dimension(); 00512 DEBUG_ERROR_IF(p2.dimension() != dim, "Dimension Mismatch"); 00513 for (int i = dim-1; i >= 0; --i) { 00514 if (p1[i] > p2[i]) return false; 00515 if (p1[i] < p2[i]) return true; 00516 } 00517 return false; 00518 } 00519 00520 } // End of namespace Go 00521 00522 00523 00524 #endif // _POINT_H 00525 00526 00527 00528 00529