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