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