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