RDKit
Open-source cheminformatics and machine learning.
point.h
Go to the documentation of this file.
1 //
2 // Copyright (C) 2003-2008 Greg Landrum and Rational Discovery LLC
3 //
4 // @@ All Rights Reserved @@
5 // This file is part of the RDKit.
6 // The contents are covered by the terms of the BSD license
7 // which is included in the file license.txt, found at the root
8 // of the RDKit source tree.
9 //
10 
11 #ifndef __RD_POINT_H__
12 #define __RD_POINT_H__
13 #include <iostream>
14 #include <cmath>
15 #include <vector>
16 #include <map>
17 
18 #ifndef M_PI
19 #define M_PI 3.14159265358979323846
20 #endif
21 
22 #include <RDGeneral/Invariant.h>
23 #include <Numerics/Vector.h>
24 #include <boost/smart_ptr.hpp>
25 
26 namespace RDGeom {
27 
28 class Point {
29  // this is the virtual base class, mandating certain functions
30  public:
31  virtual ~Point(){};
32 
33  virtual double operator[](unsigned int i) const = 0;
34  virtual double &operator[](unsigned int i) = 0;
35 
36  virtual void normalize() = 0;
37  virtual double length() const = 0;
38  virtual double lengthSq() const = 0;
39  virtual unsigned int dimension() const = 0;
40 
41  virtual Point *copy() const = 0;
42 };
43 
44 // typedef class Point3D Point;
45 class Point3D : public Point {
46  public:
47  double x, y, z;
48 
49  Point3D() : x(0.0), y(0.0), z(0.0){};
50  Point3D(double xv, double yv, double zv) : x(xv), y(yv), z(zv){};
51 
52  ~Point3D(){};
53 
54  Point3D(const Point3D &other)
55  : Point(other), x(other.x), y(other.y), z(other.z) {}
56 
57  virtual Point *copy() const { return new Point3D(*this); }
58 
59  inline unsigned int dimension() const { return 3; }
60 
61  inline double operator[](unsigned int i) const {
62  PRECONDITION(i < 3, "Invalid index on Point3D");
63  if (i == 0) {
64  return x;
65  } else if (i == 1) {
66  return y;
67  } else {
68  return z;
69  }
70  }
71 
72  inline double &operator[](unsigned int i) {
73  PRECONDITION(i < 3, "Invalid index on Point3D");
74  if (i == 0) {
75  return x;
76  } else if (i == 1) {
77  return y;
78  } else {
79  return z;
80  }
81  }
82 
83  Point3D &operator=(const Point3D &other) {
84  x = other.x;
85  y = other.y;
86  z = other.z;
87  return *this;
88  };
89 
90  Point3D &operator+=(const Point3D &other) {
91  x += other.x;
92  y += other.y;
93  z += other.z;
94  return *this;
95  };
96 
97  Point3D &operator-=(const Point3D &other) {
98  x -= other.x;
99  y -= other.y;
100  z -= other.z;
101  return *this;
102  };
103 
104  Point3D &operator*=(double scale) {
105  x *= scale;
106  y *= scale;
107  z *= scale;
108  return *this;
109  };
110 
111  Point3D &operator/=(double scale) {
112  x /= scale;
113  y /= scale;
114  z /= scale;
115  return *this;
116  };
117 
118  Point3D operator-() const {
119  Point3D res(x, y, z);
120  res.x *= -1.0;
121  res.y *= -1.0;
122  res.z *= -1.0;
123  return res;
124  }
125 
126  void normalize() {
127  double l = this->length();
128  x /= l;
129  y /= l;
130  z /= l;
131  };
132 
133  double length() const {
134  double res = x * x + y * y + z * z;
135  return sqrt(res);
136  };
137 
138  double lengthSq() const {
139  // double res = pow(x,2) + pow(y,2) + pow(z,2);
140  double res = x * x + y * y + z * z;
141  return res;
142  };
143 
144  double dotProduct(const Point3D &other) const {
145  double res = x * (other.x) + y * (other.y) + z * (other.z);
146  return res;
147  };
148 
149  /*! \brief determines the angle between a vector to this point
150  * from the origin and a vector to the other point.
151  *
152  * The angle is unsigned: the results of this call will always
153  * be between 0 and M_PI
154  */
155  double angleTo(const Point3D &other) const {
156  Point3D t1, t2;
157  t1 = *this;
158  t2 = other;
159  t1.normalize();
160  t2.normalize();
161  double dotProd = t1.dotProduct(t2);
162  // watch for roundoff error:
163  if (dotProd < -1.0)
164  dotProd = -1.0;
165  else if (dotProd > 1.0)
166  dotProd = 1.0;
167  return acos(dotProd);
168  }
169 
170  /*! \brief determines the signed angle between a vector to this point
171  * from the origin and a vector to the other point.
172  *
173  * The results of this call will be between 0 and M_2_PI
174  */
175  double signedAngleTo(const Point3D &other) const {
176  double res = this->angleTo(other);
177  // check the sign of the z component of the cross product:
178  if ((this->x * other.y - this->y * other.x) < -1e-6) res = 2.0 * M_PI - res;
179  return res;
180  }
181 
182  /*! \brief Returns a normalized direction vector from this
183  * point to another.
184  *
185  */
186  Point3D directionVector(const Point3D &other) const {
187  Point3D res;
188  res.x = other.x - x;
189  res.y = other.y - y;
190  res.z = other.z - z;
191  res.normalize();
192  return res;
193  }
194 
195  /*! \brief Cross product of this point with the another point
196  *
197  * The order is important here
198  * The result is "this" cross with "other" not (other x this)
199  */
200  Point3D crossProduct(const Point3D &other) const {
201  Point3D res;
202  res.x = y * (other.z) - z * (other.y);
203  res.y = -x * (other.z) + z * (other.x);
204  res.z = x * (other.y) - y * (other.x);
205  return res;
206  };
207 
208  /*! \brief Get a unit perpendicular from this point (treating it as a vector):
209  *
210  */
212  Point3D res(0.0, 0.0, 0.0);
213  if (x) {
214  if (y) {
215  res.y = -1 * x;
216  res.x = y;
217  } else if (z) {
218  res.z = -1 * x;
219  res.x = z;
220  } else {
221  res.y = 1;
222  }
223  } else if (y) {
224  if (z) {
225  res.z = -1 * y;
226  res.y = z;
227  } else {
228  res.x = 1;
229  }
230  } else if (z) {
231  res.x = 1;
232  }
233  double l = res.length();
234  POSTCONDITION(l > 0.0, "zero perpendicular");
235  res /= l;
236  return res;
237  }
238 };
239 
240 // given a set of four pts in 3D compute the dihedral angle between the
241 // plane of the first three points (pt1, pt2, pt3) and the plane of the
242 // last three points (pt2, pt3, pt4)
243 // the computed angle is between 0 and PI
244 double computeDihedralAngle(const Point3D &pt1, const Point3D &pt2,
245  const Point3D &pt3, const Point3D &pt4);
246 
247 // given a set of four pts in 3D compute the signed dihedral angle between the
248 // plane of the first three points (pt1, pt2, pt3) and the plane of the
249 // last three points (pt2, pt3, pt4)
250 // the computed angle is between -PI and PI
251 double computeSignedDihedralAngle(const Point3D &pt1, const Point3D &pt2,
252  const Point3D &pt3, const Point3D &pt4);
253 
254 class Point2D : public Point {
255  public:
256  double x, y;
257 
258  Point2D() : x(0.0), y(0.0){};
259  Point2D(double xv, double yv) : x(xv), y(yv){};
260 
261  ~Point2D(){};
262 
263  Point2D(const Point2D &other) : Point(other), x(other.x), y(other.y) {}
264 
265  virtual Point *copy() const { return new Point2D(*this); }
266 
267  inline unsigned int dimension() const { return 2; }
268 
269  inline double operator[](unsigned int i) const {
270  PRECONDITION(i < 2, "Invalid index on Point2D");
271  if (i == 0) {
272  return x;
273  } else {
274  return y;
275  }
276  }
277 
278  inline double &operator[](unsigned int i) {
279  PRECONDITION(i < 2, "Invalid index on Point2D");
280  if (i == 0) {
281  return x;
282  } else {
283  return y;
284  }
285  }
286 
287  Point2D &operator=(const Point2D &other) {
288  x = other.x;
289  y = other.y;
290  return *this;
291  };
292 
293  Point2D &operator+=(const Point2D &other) {
294  x += other.x;
295  y += other.y;
296  return *this;
297  };
298 
299  Point2D &operator-=(const Point2D &other) {
300  x -= other.x;
301  y -= other.y;
302  return *this;
303  };
304 
305  Point2D &operator*=(double scale) {
306  x *= scale;
307  y *= scale;
308  return *this;
309  };
310 
311  Point2D &operator/=(double scale) {
312  x /= scale;
313  y /= scale;
314  return *this;
315  };
316 
317  Point2D operator-() const {
318  Point2D res(x, y);
319  res.x *= -1.0;
320  res.y *= -1.0;
321  return res;
322  }
323 
324  void normalize() {
325  double ln = this->length();
326  x /= ln;
327  y /= ln;
328  };
329 
330  void rotate90() {
331  double temp = x;
332  x = -y;
333  y = temp;
334  }
335 
336  double length() const {
337  // double res = pow(x,2) + pow(y,2);
338  double res = x * x + y * y;
339  return sqrt(res);
340  };
341 
342  double lengthSq() const {
343  double res = x * x + y * y;
344  return res;
345  };
346 
347  double dotProduct(const Point2D &other) const {
348  double res = x * (other.x) + y * (other.y);
349  return res;
350  };
351 
352  double angleTo(const Point2D &other) const {
353  Point2D t1, t2;
354  t1 = *this;
355  t2 = other;
356  t1.normalize();
357  t2.normalize();
358  double dotProd = t1.dotProduct(t2);
359  // watch for roundoff error:
360  if (dotProd < -1.0)
361  dotProd = -1.0;
362  else if (dotProd > 1.0)
363  dotProd = 1.0;
364  return acos(dotProd);
365  }
366 
367  double signedAngleTo(const Point2D &other) const {
368  double res = this->angleTo(other);
369  if ((this->x * other.y - this->y * other.x) < -1e-6) res = 2.0 * M_PI - res;
370  return res;
371  }
372 
373  Point2D directionVector(const Point2D &other) const {
374  Point2D res;
375  res.x = other.x - x;
376  res.y = other.y - y;
377  res.normalize();
378  return res;
379  }
380 };
381 
382 class PointND : public Point {
383  public:
384  typedef boost::shared_ptr<RDNumeric::Vector<double> > VECT_SH_PTR;
385 
386  PointND(unsigned int dim) {
388  dp_storage.reset(nvec);
389  };
390 
391  PointND(const PointND &other) : Point(other) {
393  new RDNumeric::Vector<double>(*other.getStorage());
394  dp_storage.reset(nvec);
395  }
396 
397  virtual Point *copy() const { return new PointND(*this); }
398 
399 #if 0
400  template <typename T>
401  PointND(const T &vals){
402  RDNumeric::Vector<double> *nvec = new RDNumeric::Vector<double>(vals.size(), 0.0);
403  dp_storage.reset(nvec);
404  unsigned int idx=0;
405  typename T::const_iterator it;
406  for(it=vals.begin();
407  it!=vals.end();
408  ++it){
409  nvec->setVal(idx,*it);
410  ++idx;
411  };
412  };
413 #endif
414 
415  ~PointND() {}
416 
417  inline double operator[](unsigned int i) const {
418  return dp_storage.get()->getVal(i);
419  }
420 
421  inline double &operator[](unsigned int i) { return (*dp_storage.get())[i]; }
422 
423  inline void normalize() { dp_storage.get()->normalize(); }
424 
425  inline double length() const { return dp_storage.get()->normL2(); }
426 
427  inline double lengthSq() const { return dp_storage.get()->normL2Sq(); }
428 
429  unsigned int dimension() const { return dp_storage.get()->size(); }
430 
431  PointND &operator=(const PointND &other) {
433  new RDNumeric::Vector<double>(*other.getStorage());
434  dp_storage.reset(nvec);
435  return *this;
436  }
437 
438  PointND &operator+=(const PointND &other) {
439  (*dp_storage.get()) += (*other.getStorage());
440  return *this;
441  }
442 
443  PointND &operator-=(const PointND &other) {
444  (*dp_storage.get()) -= (*other.getStorage());
445  return *this;
446  }
447 
448  PointND &operator*=(double scale) {
449  (*dp_storage.get()) *= scale;
450  return *this;
451  }
452 
453  PointND &operator/=(double scale) {
454  (*dp_storage.get()) /= scale;
455  return *this;
456  }
457 
459  PRECONDITION(this->dimension() == other.dimension(),
460  "Point dimensions do not match");
461  PointND np(other);
462  np -= (*this);
463  np.normalize();
464  return np;
465  }
466 
467  double dotProduct(const PointND &other) const {
468  return dp_storage.get()->dotProduct(*other.getStorage());
469  }
470 
471  double angleTo(const PointND &other) const {
472  double dp = this->dotProduct(other);
473  double n1 = this->length();
474  double n2 = other.length();
475  if ((n1 > 1.e-8) && (n2 > 1.e-8)) {
476  dp /= (n1 * n2);
477  }
478  if (dp < -1.0)
479  dp = -1.0;
480  else if (dp > 1.0)
481  dp = 1.0;
482  return acos(dp);
483  }
484 
485  private:
486  VECT_SH_PTR dp_storage;
487  inline const RDNumeric::Vector<double> *getStorage() const {
488  return dp_storage.get();
489  }
490 };
491 
492 typedef std::vector<RDGeom::Point *> PointPtrVect;
493 typedef PointPtrVect::iterator PointPtrVect_I;
494 typedef PointPtrVect::const_iterator PointPtrVect_CI;
495 
496 typedef std::vector<RDGeom::Point3D *> Point3DPtrVect;
497 typedef std::vector<RDGeom::Point2D *> Point2DPtrVect;
498 typedef Point3DPtrVect::iterator Point3DPtrVect_I;
499 typedef Point3DPtrVect::const_iterator Point3DPtrVect_CI;
500 typedef Point2DPtrVect::iterator Point2DPtrVect_I;
501 typedef Point2DPtrVect::const_iterator Point2DPtrVect_CI;
502 
503 typedef std::vector<const RDGeom::Point3D *> Point3DConstPtrVect;
504 typedef Point3DConstPtrVect::iterator Point3DConstPtrVect_I;
505 typedef Point3DConstPtrVect::const_iterator Point3DConstPtrVect_CI;
506 
507 typedef std::vector<Point3D> POINT3D_VECT;
508 typedef std::vector<Point3D>::iterator POINT3D_VECT_I;
509 typedef std::vector<Point3D>::const_iterator POINT3D_VECT_CI;
510 
511 typedef std::map<int, Point2D> INT_POINT2D_MAP;
512 typedef INT_POINT2D_MAP::iterator INT_POINT2D_MAP_I;
513 typedef INT_POINT2D_MAP::const_iterator INT_POINT2D_MAP_CI;
514 
515 std::ostream &operator<<(std::ostream &target, const RDGeom::Point &pt);
516 
519 RDGeom::Point3D operator*(const RDGeom::Point3D &p1, double v);
520 RDGeom::Point3D operator/(const RDGeom::Point3D &p1, double v);
521 
524 RDGeom::Point2D operator*(const RDGeom::Point2D &p1, double v);
525 RDGeom::Point2D operator/(const RDGeom::Point2D &p1, double v);
526 
529 RDGeom::PointND operator*(const RDGeom::PointND &p1, double v);
530 RDGeom::PointND operator/(const RDGeom::PointND &p1, double v);
531 }
532 
533 #endif
double & operator[](unsigned int i)
Definition: point.h:72
std::vector< RDGeom::Point3D * > Point3DPtrVect
Definition: point.h:496
unsigned int dimension() const
Definition: point.h:429
virtual double length() const =0
Point3D crossProduct(const Point3D &other) const
Cross product of this point with the another point.
Definition: point.h:200
Point3D & operator/=(double scale)
Definition: point.h:111
#define POSTCONDITION(expr, mess)
Definition: Invariant.h:115
Point3DPtrVect::iterator Point3DPtrVect_I
Definition: point.h:498
Point2D & operator*=(double scale)
Definition: point.h:305
unsigned int dimension() const
Definition: point.h:267
INT_POINT2D_MAP::iterator INT_POINT2D_MAP_I
Definition: point.h:512
std::vector< RDGeom::Point * > PointPtrVect
Definition: point.h:492
#define M_PI
Definition: point.h:19
RDGeom::Point3D operator+(const RDGeom::Point3D &p1, const RDGeom::Point3D &p2)
Point2D directionVector(const Point2D &other) const
Definition: point.h:373
std::vector< Point3D >::const_iterator POINT3D_VECT_CI
Definition: point.h:509
double angleTo(const Point3D &other) const
determines the angle between a vector to this point from the origin and a vector to the other point...
Definition: point.h:155
void rotate90()
Definition: point.h:330
Point3D(const Point3D &other)
Definition: point.h:54
double y
Definition: point.h:256
Point3D getPerpendicular() const
Get a unit perpendicular from this point (treating it as a vector):
Definition: point.h:211
std::ostream & operator<<(std::ostream &target, const RDGeom::Point &pt)
PointPtrVect::const_iterator PointPtrVect_CI
Definition: point.h:494
double length() const
Definition: point.h:425
unsigned int dimension() const
Definition: point.h:59
double dotProduct(const Point3D &other) const
Definition: point.h:144
double & operator[](unsigned int i)
Definition: point.h:421
double z
Definition: point.h:47
Point2D & operator=(const Point2D &other)
Definition: point.h:287
double signedAngleTo(const Point2D &other) const
Definition: point.h:367
boost::shared_ptr< RDNumeric::Vector< double > > VECT_SH_PTR
Definition: point.h:384
virtual double operator[](unsigned int i) const =0
Point2D & operator/=(double scale)
Definition: point.h:311
double & operator[](unsigned int i)
Definition: point.h:278
PointND directionVector(const PointND &other)
Definition: point.h:458
Point3DConstPtrVect::const_iterator Point3DConstPtrVect_CI
Definition: point.h:505
std::vector< Point3D > POINT3D_VECT
Definition: point.h:507
void normalize()
Definition: point.h:423
Point2D(const Point2D &other)
Definition: point.h:263
Point2DPtrVect::iterator Point2DPtrVect_I
Definition: point.h:500
void setVal(unsigned int i, TYPE val)
sets the index at a particular value
Definition: Vector.h:86
virtual Point * copy() const
Definition: point.h:57
PointND & operator-=(const PointND &other)
Definition: point.h:443
RDGeom::Point3D operator-(const RDGeom::Point3D &p1, const RDGeom::Point3D &p2)
double computeDihedralAngle(const Point3D &pt1, const Point3D &pt2, const Point3D &pt3, const Point3D &pt4)
std::vector< const RDGeom::Point3D * > Point3DConstPtrVect
Definition: point.h:503
PointND & operator*=(double scale)
Definition: point.h:448
std::map< int, Point2D > INT_POINT2D_MAP
Definition: point.h:511
RDGeom::Point3D operator*(const RDGeom::Point3D &p1, double v)
void normalize()
Definition: point.h:126
Point3D(double xv, double yv, double zv)
Definition: point.h:50
virtual void normalize()=0
Point3D & operator+=(const Point3D &other)
Definition: point.h:90
virtual Point * copy() const
Definition: point.h:265
double dotProduct(const Point2D &other) const
Definition: point.h:347
double signedAngleTo(const Point3D &other) const
determines the signed angle between a vector to this point from the origin and a vector to the other ...
Definition: point.h:175
double angleTo(const PointND &other) const
Definition: point.h:471
Point2DPtrVect::const_iterator Point2DPtrVect_CI
Definition: point.h:501
PointND & operator=(const PointND &other)
Definition: point.h:431
void normalize()
Definition: point.h:324
double dotProduct(const PointND &other) const
Definition: point.h:467
Point3DPtrVect::const_iterator Point3DPtrVect_CI
Definition: point.h:499
virtual ~Point()
Definition: point.h:31
double operator[](unsigned int i) const
Definition: point.h:417
virtual double lengthSq() const =0
PointND & operator+=(const PointND &other)
Definition: point.h:438
Point3D & operator*=(double scale)
Definition: point.h:104
virtual Point * copy() const =0
Point3D directionVector(const Point3D &other) const
Returns a normalized direction vector from this point to another.
Definition: point.h:186
double computeSignedDihedralAngle(const Point3D &pt1, const Point3D &pt2, const Point3D &pt3, const Point3D &pt4)
double length() const
Definition: point.h:336
Point2D(double xv, double yv)
Definition: point.h:259
double angleTo(const Point2D &other) const
Definition: point.h:352
Point2D & operator-=(const Point2D &other)
Definition: point.h:299
Point3D & operator-=(const Point3D &other)
Definition: point.h:97
double y
Definition: point.h:47
Point2D & operator+=(const Point2D &other)
Definition: point.h:293
PointND(unsigned int dim)
Definition: point.h:386
#define PRECONDITION(expr, mess)
Definition: Invariant.h:107
PointND & operator/=(double scale)
Definition: point.h:453
RDGeom::Point3D operator/(const RDGeom::Point3D &p1, double v)
double x
Definition: point.h:256
virtual unsigned int dimension() const =0
Point3D & operator=(const Point3D &other)
Definition: point.h:83
double lengthSq() const
Definition: point.h:138
INT_POINT2D_MAP::const_iterator INT_POINT2D_MAP_CI
Definition: point.h:513
double length() const
Definition: point.h:133
Point3D operator-() const
Definition: point.h:118
double x
Definition: point.h:47
double lengthSq() const
Definition: point.h:427
double operator[](unsigned int i) const
Definition: point.h:269
virtual Point * copy() const
Definition: point.h:397
std::vector< RDGeom::Point2D * > Point2DPtrVect
Definition: point.h:497
A class to represent vectors of numbers.
Definition: Vector.h:28
double operator[](unsigned int i) const
Definition: point.h:61
double lengthSq() const
Definition: point.h:342
Point3DConstPtrVect::iterator Point3DConstPtrVect_I
Definition: point.h:504
std::vector< Point3D >::iterator POINT3D_VECT_I
Definition: point.h:508
PointND(const PointND &other)
Definition: point.h:391
Point2D operator-() const
Definition: point.h:317
PointPtrVect::iterator PointPtrVect_I
Definition: point.h:493