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