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
246  const Point3D &pt2,
247  const Point3D &pt3,
248  const Point3D &pt4);
249 
250 // given a set of four pts in 3D compute the signed dihedral angle between the
251 // plane of the first three points (pt1, pt2, pt3) and the plane of the
252 // last three points (pt2, pt3, pt4)
253 // the computed angle is between -PI and PI
255  const Point3D &pt1, const Point3D &pt2, const Point3D &pt3,
256  const Point3D &pt4);
257 
259  public:
260  double x, y;
261 
262  Point2D() : x(0.0), y(0.0){};
263  Point2D(double xv, double yv) : x(xv), y(yv){};
264  ~Point2D(){};
265 
266  Point2D(const Point2D &other) : Point(other), x(other.x), y(other.y) {}
267  //! construct from a Point3D (ignoring the z coordinate)
268  Point2D(const Point3D &p3d) : Point(p3d), x(p3d.x), y(p3d.y){};
269 
270  virtual Point *copy() const { return new Point2D(*this); }
271 
272  inline unsigned int dimension() const { return 2; }
273 
274  inline double operator[](unsigned int i) const {
275  PRECONDITION(i < 2, "Invalid index on Point2D");
276  if (i == 0) {
277  return x;
278  } else {
279  return y;
280  }
281  }
282 
283  inline double &operator[](unsigned int i) {
284  PRECONDITION(i < 2, "Invalid index on Point2D");
285  if (i == 0) {
286  return x;
287  } else {
288  return y;
289  }
290  }
291 
292  Point2D &operator=(const Point2D &other) {
293  x = other.x;
294  y = other.y;
295  return *this;
296  };
297 
298  Point2D &operator+=(const Point2D &other) {
299  x += other.x;
300  y += other.y;
301  return *this;
302  };
303 
304  Point2D &operator-=(const Point2D &other) {
305  x -= other.x;
306  y -= other.y;
307  return *this;
308  };
309 
310  Point2D &operator*=(double scale) {
311  x *= scale;
312  y *= scale;
313  return *this;
314  };
315 
316  Point2D &operator/=(double scale) {
317  x /= scale;
318  y /= scale;
319  return *this;
320  };
321 
322  Point2D operator-() const {
323  Point2D res(x, y);
324  res.x *= -1.0;
325  res.y *= -1.0;
326  return res;
327  }
328 
329  void normalize() {
330  double ln = this->length();
331  x /= ln;
332  y /= ln;
333  };
334 
335  void rotate90() {
336  double temp = x;
337  x = -y;
338  y = temp;
339  }
340 
341  double length() const {
342  // double res = pow(x,2) + pow(y,2);
343  double res = x * x + y * y;
344  return sqrt(res);
345  };
346 
347  double lengthSq() const {
348  double res = x * x + y * y;
349  return res;
350  };
351 
352  double dotProduct(const Point2D &other) const {
353  double res = x * (other.x) + y * (other.y);
354  return res;
355  };
356 
357  double angleTo(const Point2D &other) const {
358  Point2D t1, t2;
359  t1 = *this;
360  t2 = other;
361  t1.normalize();
362  t2.normalize();
363  double dotProd = t1.dotProduct(t2);
364  // watch for roundoff error:
365  if (dotProd < -1.0)
366  dotProd = -1.0;
367  else if (dotProd > 1.0)
368  dotProd = 1.0;
369  return acos(dotProd);
370  }
371 
372  double signedAngleTo(const Point2D &other) const {
373  double res = this->angleTo(other);
374  if ((this->x * other.y - this->y * other.x) < -1e-6) res = 2.0 * M_PI - res;
375  return res;
376  }
377 
378  Point2D directionVector(const Point2D &other) const {
379  Point2D res;
380  res.x = other.x - x;
381  res.y = other.y - y;
382  res.normalize();
383  return res;
384  }
385 };
386 
388  public:
389  typedef boost::shared_ptr<RDNumeric::Vector<double>> VECT_SH_PTR;
390 
391  PointND(unsigned int dim) {
393  dp_storage.reset(nvec);
394  };
395 
396  PointND(const PointND &other) : Point(other) {
398  new RDNumeric::Vector<double>(*other.getStorage());
399  dp_storage.reset(nvec);
400  }
401 
402  virtual Point *copy() const { return new PointND(*this); }
403 
404 #if 0
405  template <typename T>
406  PointND(const T &vals){
407  RDNumeric::Vector<double> *nvec = new RDNumeric::Vector<double>(vals.size(), 0.0);
408  dp_storage.reset(nvec);
409  unsigned int idx=0;
410  typename T::const_iterator it;
411  for(it=vals.begin();
412  it!=vals.end();
413  ++it){
414  nvec->setVal(idx,*it);
415  ++idx;
416  };
417  };
418 #endif
419 
420  ~PointND() {}
421 
422  inline double operator[](unsigned int i) const {
423  return dp_storage.get()->getVal(i);
424  }
425 
426  inline double &operator[](unsigned int i) { return (*dp_storage.get())[i]; }
427 
428  inline void normalize() { dp_storage.get()->normalize(); }
429 
430  inline double length() const { return dp_storage.get()->normL2(); }
431 
432  inline double lengthSq() const { return dp_storage.get()->normL2Sq(); }
433 
434  unsigned int dimension() const { return dp_storage.get()->size(); }
435 
436  PointND &operator=(const PointND &other) {
437  if (this == &other) return *this;
438 
440  new RDNumeric::Vector<double>(*other.getStorage());
441  dp_storage.reset(nvec);
442  return *this;
443  }
444 
445  PointND &operator+=(const PointND &other) {
446  (*dp_storage.get()) += (*other.getStorage());
447  return *this;
448  }
449 
450  PointND &operator-=(const PointND &other) {
451  (*dp_storage.get()) -= (*other.getStorage());
452  return *this;
453  }
454 
455  PointND &operator*=(double scale) {
456  (*dp_storage.get()) *= scale;
457  return *this;
458  }
459 
460  PointND &operator/=(double scale) {
461  (*dp_storage.get()) /= scale;
462  return *this;
463  }
464 
466  PRECONDITION(this->dimension() == other.dimension(),
467  "Point dimensions do not match");
468  PointND np(other);
469  np -= (*this);
470  np.normalize();
471  return np;
472  }
473 
474  double dotProduct(const PointND &other) const {
475  return dp_storage.get()->dotProduct(*other.getStorage());
476  }
477 
478  double angleTo(const PointND &other) const {
479  double dp = this->dotProduct(other);
480  double n1 = this->length();
481  double n2 = other.length();
482  if ((n1 > 1.e-8) && (n2 > 1.e-8)) {
483  dp /= (n1 * n2);
484  }
485  if (dp < -1.0)
486  dp = -1.0;
487  else if (dp > 1.0)
488  dp = 1.0;
489  return acos(dp);
490  }
491 
492  private:
493  VECT_SH_PTR dp_storage;
494  inline const RDNumeric::Vector<double> *getStorage() const {
495  return dp_storage.get();
496  }
497 };
498 
499 typedef std::vector<RDGeom::Point *> PointPtrVect;
500 typedef PointPtrVect::iterator PointPtrVect_I;
501 typedef PointPtrVect::const_iterator PointPtrVect_CI;
502 
503 typedef std::vector<RDGeom::Point3D *> Point3DPtrVect;
504 typedef std::vector<RDGeom::Point2D *> Point2DPtrVect;
505 typedef Point3DPtrVect::iterator Point3DPtrVect_I;
506 typedef Point3DPtrVect::const_iterator Point3DPtrVect_CI;
507 typedef Point2DPtrVect::iterator Point2DPtrVect_I;
508 typedef Point2DPtrVect::const_iterator Point2DPtrVect_CI;
509 
510 typedef std::vector<const RDGeom::Point3D *> Point3DConstPtrVect;
511 typedef Point3DConstPtrVect::iterator Point3DConstPtrVect_I;
512 typedef Point3DConstPtrVect::const_iterator Point3DConstPtrVect_CI;
513 
514 typedef std::vector<Point3D> POINT3D_VECT;
515 typedef std::vector<Point3D>::iterator POINT3D_VECT_I;
516 typedef std::vector<Point3D>::const_iterator POINT3D_VECT_CI;
517 
518 typedef std::map<int, Point2D> INT_POINT2D_MAP;
519 typedef INT_POINT2D_MAP::iterator INT_POINT2D_MAP_I;
520 typedef INT_POINT2D_MAP::const_iterator INT_POINT2D_MAP_CI;
521 
522 RDKIT_RDGEOMETRYLIB_EXPORT std::ostream &operator<<(std::ostream &target,
523  const RDGeom::Point &pt);
524 
526  const RDGeom::Point3D &p2);
528  const RDGeom::Point3D &p2);
530  double v);
532  double v);
533 
535  const RDGeom::Point2D &p2);
537  const RDGeom::Point2D &p2);
539  double v);
541  double v);
542 
544  const RDGeom::PointND &p2);
546  const RDGeom::PointND &p2);
548  double v);
550  double v);
551 } // namespace RDGeom
552 
553 #endif
double & operator[](unsigned int i)
Definition: point.h:73
std::vector< RDGeom::Point3D * > Point3DPtrVect
Definition: point.h:503
Point3D & operator/=(double scale)
Definition: point.h:112
#define RDKIT_RDGEOMETRYLIB_EXPORT
Definition: export.h:515
#define POSTCONDITION(expr, mess)
Definition: Invariant.h:117
Point3DPtrVect::iterator Point3DPtrVect_I
Definition: point.h:505
Point2D & operator*=(double scale)
Definition: point.h:310
double length() const
Definition: point.h:430
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:519
std::vector< RDGeom::Point * > PointPtrVect
Definition: point.h:499
#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:516
void rotate90()
Definition: point.h:335
Point3D(const Point3D &other)
Definition: point.h:55
double y
Definition: point.h:260
Point2D(const Point3D &p3d)
construct from a Point3D (ignoring the z coordinate)
Definition: point.h:268
PointPtrVect::const_iterator PointPtrVect_CI
Definition: point.h:501
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:426
double z
Definition: point.h:48
Point2D & operator=(const Point2D &other)
Definition: point.h:292
double angleTo(const Point2D &other) const
Definition: point.h:357
boost::shared_ptr< RDNumeric::Vector< double > > VECT_SH_PTR
Definition: point.h:389
Point2D & operator/=(double scale)
Definition: point.h:316
double & operator[](unsigned int i)
Definition: point.h:283
PointND directionVector(const PointND &other)
Definition: point.h:465
virtual Point * copy() const
Definition: point.h:402
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:512
std::vector< Point3D > POINT3D_VECT
Definition: point.h:514
void normalize()
Definition: point.h:428
Point2D(const Point2D &other)
Definition: point.h:266
Point2DPtrVect::iterator Point2DPtrVect_I
Definition: point.h:507
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:450
std::vector< const RDGeom::Point3D * > Point3DConstPtrVect
Definition: point.h:510
RDKIT_RDGEOMETRYLIB_EXPORT std::ostream & operator<<(std::ostream &target, const RDGeom::Point &pt)
PointND & operator*=(double scale)
Definition: point.h:455
std::map< int, Point2D > INT_POINT2D_MAP
Definition: point.h:518
void normalize()
Definition: point.h:127
double angleTo(const PointND &other) const
Definition: point.h:478
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:352
Point3D & operator+=(const Point3D &other)
Definition: point.h:91
Point2DPtrVect::const_iterator Point2DPtrVect_CI
Definition: point.h:508
double dotProduct(const PointND &other) const
Definition: point.h:474
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:436
unsigned int dimension() const
Definition: point.h:272
void normalize()
Definition: point.h:329
Point3DPtrVect::const_iterator Point3DPtrVect_CI
Definition: point.h:506
virtual Point * copy() const
Definition: point.h:58
virtual ~Point()
Definition: point.h:32
PointND & operator+=(const PointND &other)
Definition: point.h:445
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:263
Point2D & operator-=(const Point2D &other)
Definition: point.h:304
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:270
Point3D & operator-=(const Point3D &other)
Definition: point.h:98
unsigned int dimension() const
Definition: point.h:434
double y
Definition: point.h:48
Point2D & operator+=(const Point2D &other)
Definition: point.h:298
Point2D directionVector(const Point2D &other) const
Definition: point.h:378
double operator[](unsigned int i) const
Definition: point.h:274
Point2D operator-() const
Definition: point.h:322
PointND(unsigned int dim)
Definition: point.h:391
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:109
PointND & operator/=(double scale)
Definition: point.h:460
double x
Definition: point.h:260
double lengthSq() const
Definition: point.h:432
Point3D & operator=(const Point3D &other)
Definition: point.h:84
INT_POINT2D_MAP::const_iterator INT_POINT2D_MAP_CI
Definition: point.h:520
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:347
std::vector< RDGeom::Point2D * > Point2DPtrVect
Definition: point.h:504
double length() const
Definition: point.h:341
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:422
double signedAngleTo(const Point2D &other) const
Definition: point.h:372
Point3DConstPtrVect::iterator Point3DConstPtrVect_I
Definition: point.h:511
double dotProduct(const Point3D &other) const
Definition: point.h:145
std::vector< Point3D >::iterator POINT3D_VECT_I
Definition: point.h:515
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:396
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:500
RDKIT_RDGEOMETRYLIB_EXPORT double computeDihedralAngle(const Point3D &pt1, const Point3D &pt2, const Point3D &pt3, const Point3D &pt4)