1 // $Id$
2 //
3 // Copyright (C) 2003-2006 Rational Discovery LLC
4 //
5 // @@ All Rights Reserved @@
6 // This file is part of the RDKit.
7 // The contents are covered by the terms of the BSD license
8 // which is included in the file license.txt, found at the root
9 // of the RDKit source tree.
10 //
11 #include <RDGeneral/Invariant.h>
12 #include <cstring>
13 #include "Transform.h"
14 #include "Transform2D.h"
15 #include <cmath>
16 #include "point.h"
17
18 namespace RDGeom {
19
setToIdentity()20 void Transform2D::setToIdentity() {
21 unsigned int i, id;
22 double *data = d_data.get();
23 memset(static_cast<void *>(data), 0, d_dataSize * sizeof(double));
24 for (i = 0; i < DIM_2D; i++) {
25 id = i * (DIM_2D + 1);
26 data[id] = 1.0;
27 }
28 }
29
TransformPoint(Point2D & pt) const30 void Transform2D::TransformPoint(Point2D &pt) const {
31 double x, y;
32 double *data = d_data.get();
33 x = data[0] * pt.x + data[1] * pt.y + data[2];
34 y = data[3] * pt.x + data[4] * pt.y + data[5];
35
36 pt.x = x;
37 pt.y = y;
38 }
39
SetTranslation(const Point2D & pt)40 void Transform2D::SetTranslation(const Point2D &pt) {
41 unsigned int i = DIM_2D - 1;
42 double *data = d_data.get();
43 data[i] = pt.x;
44 i += DIM_2D;
45 data[i] = pt.y;
46 i += DIM_2D;
47 data[i] = 1.0;
48 }
49
SetTransform(const Point2D & pt,double angle)50 void Transform2D::SetTransform(const Point2D &pt, double angle) {
51 this->setToIdentity();
52
53 Transform2D trans1;
54 trans1.SetTranslation(-pt);
55 double *data = d_data.get();
56 // set the rotation
57 data[0] = cos(angle);
58 data[1] = -sin(angle);
59 data[3] = sin(angle);
60 data[4] = cos(angle);
61
62 (*this) *= trans1;
63
64 // translation back to the original coordinate
65 Transform2D trans2;
66 trans2.SetTranslation(pt);
67 trans2 *= (*this);
68
69 // now combine them
70 this->assign(trans2);
71 }
72
SetTransform(const Point2D & ref1,const Point2D & ref2,const Point2D & pt1,const Point2D & pt2)73 void Transform2D::SetTransform(const Point2D &ref1, const Point2D &ref2,
74 const Point2D &pt1, const Point2D &pt2) {
75 // compute the angle between the two vectors
76 Point2D rvec = ref2 - ref1;
77 Point2D pvec = pt2 - pt1;
78
79 double dp = rvec.dotProduct(pvec);
80 double lp = (rvec.length()) * (pvec.length());
81 if (lp <= 0.0) {
82 this->setToIdentity();
83 return;
84 }
85
86 double cval = dp / lp;
87 if (cval < -1.0) {
88 cval = -1.0;
89 } else if (cval > 1.0) {
90 cval = 1.0;
91 }
92
93 double ang = acos(cval);
94
95 // figure out if we have to do clock wise or anti clock wise rotation
96 double cross = (pvec.x) * (rvec.y) - (pvec.y) * (rvec.x);
97
98 if (cross < 0.0) {
99 ang *= -1.0;
100 }
101 this->setToIdentity();
102 // set the rotation
103 double *data = d_data.get();
104 data[0] = cos(ang);
105 data[1] = -sin(ang);
106 data[3] = sin(ang);
107 data[4] = cos(ang);
108
109 // apply this rotation to pt1 and compute the translation
110 Point2D npt1 = pt1;
111 this->TransformPoint(npt1);
112 data[DIM_2D - 1] = ref1.x - npt1.x;
113 data[2 * DIM_2D - 1] = ref1.y - npt1.y;
114 }
115 } // namespace RDGeom
116
operator *(const RDGeom::Transform2D & t1,const RDGeom::Transform2D & t2)117 RDGeom::Transform2D operator*(const RDGeom::Transform2D &t1,
118 const RDGeom::Transform2D &t2) {
119 RDGeom::Transform2D res;
120 RDNumeric::multiply(t1, t2, res);
121 return res;
122 };
123