1 // $Id$
2 //
3 // Copyright (C) 2005-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
12 #include <RDGeneral/test.h>
13 #include <RDGeneral/types.h>
14 #include <RDGeneral/Invariant.h>
15 #include <RDGeneral/utils.h>
16 #include "Transform2D.h"
17 #include "Transform3D.h"
18 #include "point.h"
19 #include <cstdlib>
20 #include <ctime>
21
22 using namespace RDGeom;
23 using namespace std;
24
ptEq(const Point3D pt1,const Point3D pt2,double val=1.e-8)25 bool ptEq(const Point3D pt1, const Point3D pt2, double val = 1.e-8) {
26 return ((fabs(pt1.x - pt2.x) < val) && (fabs(pt1.y - pt2.y) < val) &&
27 (fabs(pt1.z - pt2.z) < val));
28 }
29
ptEq(const Point2D pt1,const Point2D pt2,double val=1.e-8)30 bool ptEq(const Point2D pt1, const Point2D pt2, double val = 1.e-8) {
31 return ((fabs(pt1.x - pt2.x) < val) && (fabs(pt1.y - pt2.y) < val));
32 }
33
randNum(double x=5)34 double randNum(double x = 5) {
35 auto res = (double)rand();
36 res /= RAND_MAX;
37 res *= x;
38 return res;
39 }
40
testPointND()41 void testPointND() {
42 PointND pt(5);
43 TEST_ASSERT(pt.dimension() == 5);
44 unsigned int i;
45 for (i = 0; i < 5; ++i) {
46 pt[i] = i + 1.0;
47 }
48 pt.normalize();
49 double ep[5];
50 ep[0] = 0.13484;
51 ep[1] = 0.26968;
52 ep[2] = 0.40452;
53 ep[3] = 0.53936;
54 ep[4] = 0.6742;
55
56 for (i = 0; i < 5; ++i) {
57 TEST_ASSERT(fabs(pt[i] - ep[i]) < 1.e-4);
58 }
59
60 PointND pt2(pt);
61 for (i = 0; i < 5; ++i) {
62 TEST_ASSERT(fabs(pt2[i] - ep[i]) < 1.e-4);
63 }
64
65 pt2 += pt;
66 for (i = 0; i < 5; ++i) {
67 TEST_ASSERT(fabs(pt2[i] - 2 * ep[i]) < 1.e-4);
68 }
69
70 pt2 /= 2.0;
71 for (i = 0; i < 5; ++i) {
72 TEST_ASSERT(fabs(pt2[i] - ep[i]) < 1.e-4);
73 }
74
75 pt2 -= pt;
76 for (i = 0; i < 5; ++i) {
77 TEST_ASSERT(fabs(pt2[i] - 0.0) < 1.e-4);
78 }
79 pt2 = pt;
80 pt2 *= 2.;
81 for (i = 0; i < 5; ++i) {
82 TEST_ASSERT(fabs(pt2[i] - 2 * ep[i]) < 1.e-4);
83 }
84
85 double dp = pt.dotProduct(pt2);
86 TEST_ASSERT(fabs(dp - 2.0) < 1.e-4);
87
88 double angle = pt.angleTo(pt2);
89 TEST_ASSERT(fabs(angle - 0.0) < 1.e-4);
90 }
91
testPointOps3D()92 void testPointOps3D() {
93 Point3D pt0(1, 0, 0);
94 Point3D pt1(0, 1, 0);
95 Point3D pt2(-1, 0, 0);
96 Point3D pt3(0, -1, 0);
97
98 TEST_ASSERT(fabs(pt0.angleTo(pt0)) < 1e-4);
99 TEST_ASSERT(fabs(pt0.angleTo(pt1) - M_PI / 2.) < 1e-4);
100 TEST_ASSERT(fabs(pt0.angleTo(pt2) - M_PI) < 1e-4);
101 TEST_ASSERT(fabs(pt0.angleTo(pt3) - M_PI / 2.) < 1e-4);
102
103 TEST_ASSERT(fabs(pt1.angleTo(pt0) - M_PI / 2.) < 1e-4);
104 TEST_ASSERT(fabs(pt1.angleTo(pt1)) < 1e-4);
105 TEST_ASSERT(fabs(pt1.angleTo(pt2) - M_PI / 2.) < 1e-4);
106 TEST_ASSERT(fabs(pt1.angleTo(pt3) - M_PI) < 1e-4);
107
108 TEST_ASSERT(fabs(pt0.signedAngleTo(pt0)) < 1e-4);
109 TEST_ASSERT(fabs(pt0.signedAngleTo(pt1) - M_PI / 2.) < 1e-4);
110 TEST_ASSERT(fabs(pt0.signedAngleTo(pt2) - M_PI) < 1e-4);
111 TEST_ASSERT(fabs(pt0.signedAngleTo(pt3) - 3. * M_PI / 2.) < 1e-4);
112
113 Point3D diffPt = pt0.directionVector(pt1);
114 Point3D ref(-sqrt(2.) / 2., sqrt(2.) / 2., 0);
115 TEST_ASSERT(ptEq(diffPt, ref));
116 }
testPointOps2D()117 void testPointOps2D() {
118 Point2D pt0(1, 0);
119 Point2D pt1(0, 1);
120 Point2D pt2(-1, 0);
121 Point2D pt3(0, -1);
122
123 TEST_ASSERT(fabs(pt0.angleTo(pt0)) < 1e-4);
124 TEST_ASSERT(fabs(pt0.angleTo(pt1) - M_PI / 2.) < 1e-4);
125 TEST_ASSERT(fabs(pt0.angleTo(pt2) - M_PI) < 1e-4);
126 TEST_ASSERT(fabs(pt0.angleTo(pt3) - M_PI / 2.) < 1e-4);
127
128 TEST_ASSERT(fabs(pt1.angleTo(pt0) - M_PI / 2.) < 1e-4);
129 TEST_ASSERT(fabs(pt1.angleTo(pt1)) < 1e-4);
130 TEST_ASSERT(fabs(pt1.angleTo(pt2) - M_PI / 2.) < 1e-4);
131 TEST_ASSERT(fabs(pt1.angleTo(pt3) - M_PI) < 1e-4);
132
133 TEST_ASSERT(fabs(pt0.signedAngleTo(pt0)) < 1e-4);
134 TEST_ASSERT(fabs(pt0.signedAngleTo(pt1) - M_PI / 2.) < 1e-4);
135 TEST_ASSERT(fabs(pt0.signedAngleTo(pt2) - M_PI) < 1e-4);
136 TEST_ASSERT(fabs(pt0.signedAngleTo(pt3) - 3. * M_PI / 2.) < 1e-4);
137
138 Point2D diffPt = pt0.directionVector(pt1);
139 Point2D ref(-sqrt(2.) / 2., sqrt(2.) / 2.);
140 TEST_ASSERT(ptEq(diffPt, ref));
141 }
142
test12D()143 void test12D() {
144 Point2D pt(1.0, 2.0);
145 Transform2D trans;
146 trans.TransformPoint(pt);
147
148 CHECK_INVARIANT(fabs(pt.x - 1.0) < 1.e-8, "");
149 CHECK_INVARIANT(fabs(pt.y - 2.0) < 1.e-8, "");
150
151 Point2D ref1(randNum(), randNum());
152 Point2D ref2(randNum(), randNum());
153
154 std::cout << "ref1: " << ref1 << " ref2: " << ref2 << "\n";
155
156 Point2D pt1(randNum(), randNum());
157 Point2D pt2(randNum(), randNum());
158 Point2D pt1o = pt1;
159 Point2D pt2o = pt2;
160 std::cout << "pt1: " << pt1 << " pt2: " << pt2 << "\n";
161
162 Transform2D t2d;
163 t2d.SetTransform(ref1, ref2, pt1, pt2);
164 t2d.TransformPoint(pt1);
165 t2d.TransformPoint(pt2);
166
167 // make sure pt1 overlaps ref1
168 Point2D dif1 = pt1 - ref1;
169 CHECK_INVARIANT(fabs(dif1.x) < 1.e-8, "");
170 CHECK_INVARIANT(fabs(dif1.y) < 1.e-8, "");
171
172 // now check that the angle between the two vectors (ref2 - ref1) and
173 // (pt2 - pt1) is zero
174 Point2D rvec = ref2 - ref1;
175 Point2D pvec = pt2 - pt1;
176 rvec.normalize();
177 pvec.normalize();
178 double pdot = rvec.dotProduct(pvec);
179 CHECK_INVARIANT(fabs(pdot - 1.0) < 1.e-8, "");
180
181 // compute the reverse transform and make sure we are basically getting the
182 // identity
183 Transform2D tdi;
184 tdi.SetTransform(pt1o, pt2o, pt1, pt2);
185 tdi.TransformPoint(pt1);
186 tdi.TransformPoint(pt2);
187
188 CHECK_INVARIANT(ptEq(pt1, pt1o), "");
189 CHECK_INVARIANT(ptEq(pt2, pt2o), "");
190
191 // the following product should result in an identity matrix
192 tdi *= t2d;
193
194 tdi.TransformPoint(pt1);
195 tdi.TransformPoint(pt2);
196
197 CHECK_INVARIANT(ptEq(pt1, pt1o), "");
198 CHECK_INVARIANT(ptEq(pt2, pt2o), "");
199
200 Point2D npt1(1.0, 0.0);
201 Point2D npt2(5.0, 0.0);
202 Point2D opt1 = npt1;
203 Point2D opt2(1.0, 4.0);
204 Transform2D ntd;
205 ntd.SetTransform(npt1, M_PI / 2);
206 ntd.TransformPoint(npt1);
207 ntd.TransformPoint(npt2);
208
209 CHECK_INVARIANT(ptEq(npt1, opt1), "");
210 CHECK_INVARIANT(ptEq(npt2, opt2), "");
211 }
212
test23D()213 void test23D() {
214 Point3D pt(1.0, 0.0, 0.0);
215 Point3D tpt = pt;
216 Transform3D trans;
217 trans.SetRotation(M_PI / 2., X_Axis);
218 trans.TransformPoint(pt);
219 CHECK_INVARIANT(ptEq(tpt, pt), "");
220
221 Point3D pt2(0.0, 1.0, 0.0);
222 Point3D tpt2(0.0, 0.0, 1.0);
223 trans.TransformPoint(pt2);
224 CHECK_INVARIANT(ptEq(tpt2, pt2), "");
225
226 Point3D pt3(0.0, 0.0, 1.0);
227 Point3D tpt3(0.0, -1.0, 0.0);
228 trans.TransformPoint(pt3);
229 CHECK_INVARIANT(ptEq(tpt3, pt3), "");
230
231 // rotate around y-axis
232 Transform3D transy;
233 transy.SetRotation(M_PI / 2., Y_Axis);
234 transy.TransformPoint(pt);
235 Point3D tpt4(0.0, 0.0, -1.0);
236 CHECK_INVARIANT(ptEq(tpt4, pt), "");
237
238 Point3D pt5(0.0, 1.0, 0.0);
239 Point3D tpt5(0.0, 1.0, 0.0);
240 transy.TransformPoint(pt5);
241 CHECK_INVARIANT(ptEq(tpt5, pt5), "");
242
243 Point3D pt6(0.0, 0.0, 1.0);
244 Point3D tpt6(1.0, 0.0, 0.0);
245 transy.TransformPoint(pt6);
246 CHECK_INVARIANT(ptEq(tpt6, pt6), "");
247
248 // z-axis
249 Transform3D transz;
250 transz.SetRotation(M_PI / 2., Z_Axis);
251 Point3D pt7(1.0, 0.0, 0.0);
252 Point3D tpt7(0.0, 1.0, 0.0);
253 transz.TransformPoint(pt7);
254 CHECK_INVARIANT(ptEq(tpt7, pt7), "");
255
256 Point3D pt8(0.0, 1.0, 0.0);
257 Point3D tpt8(-1.0, 0.0, 0.0);
258 transz.TransformPoint(pt8);
259 CHECK_INVARIANT(ptEq(tpt8, pt8), "");
260
261 Point3D pt9(0.0, 0.0, 1.0);
262 Point3D tpt9(0.0, 0.0, 1.0);
263 transz.TransformPoint(pt9);
264 CHECK_INVARIANT(ptEq(tpt9, pt9), "");
265 }
266
test3MatMultiply()267 void test3MatMultiply() {
268 // start with line on the axis starting at 1.0,
269 // transform it into a line on z-axis starting at 3.0
270 Point3D pt1(1.0, 0.0, 0.0);
271 Point3D pt2(2.0, 0.0, 0.0);
272
273 std::cout << "Pt1: " << pt1 << " Pt2: " << pt2 << "\n";
274 std::cout << "-Pt1: " << (-pt1) << "\n";
275 // move to origin
276 Transform3D t1;
277 t1.SetTranslation(-pt1);
278
279 Point3D tp1 = pt1;
280 Point3D tp2 = pt1;
281 t1.TransformPoint(tp1);
282 t1.TransformPoint(tp2);
283 std::cout << "tp1: " << tp1 << " tp2: " << tp2 << "\n";
284
285 // rotation around origin
286 Transform3D t2;
287 t2.SetRotation(-M_PI / 2.0, Y_Axis);
288
289 t2.TransformPoint(tp1);
290 t2.TransformPoint(tp2);
291 std::cout << "tp1: " << tp1 << " tp2: " << tp2 << "\n";
292
293 // move on z-axis
294 Transform3D t3;
295 Point3D npt1(0.0, 0.0, 3.0);
296 t3.SetTranslation(npt1);
297 Point3D npt2(0.0, 0.0, 4.0);
298
299 t3.TransformPoint(tp1);
300 t3.TransformPoint(tp2);
301 std::cout << "tp1: " << tp1 << " tp2: " << tp2 << "\n";
302
303 std::cout << "npt1: " << npt1 << " npt2: " << npt2 << "\n";
304
305 // combine the transform;
306
307 Transform3D t4 = t3 * t2 * t1;
308 t2 *= t1;
309 t3 *= t2;
310
311 Point3D opt1 = pt1;
312 Point3D opt2 = pt2;
313
314 t3.TransformPoint(pt1);
315 t3.TransformPoint(pt2);
316 std::cout << "Pt1: " << pt1 << " Pt2: " << pt2 << "\n";
317 // check the transformed points align with the new points on z-axis
318 CHECK_INVARIANT(ptEq(pt1, npt1), "");
319 CHECK_INVARIANT(ptEq(pt2, npt2), "");
320
321 t4.TransformPoint(opt1);
322 t4.TransformPoint(opt2);
323 CHECK_INVARIANT(ptEq(opt1, npt1), "");
324 CHECK_INVARIANT(ptEq(opt2, npt2), "");
325 }
326
testFromQuaternion()327 void testFromQuaternion() {
328 double qt[4];
329 qt[0] = cos(M_PI / 6);
330 qt[1] = -sin(M_PI / 6);
331 qt[2] = 0.0;
332 qt[3] = 0.0;
333 Transform3D trans;
334 trans.SetRotationFromQuaternion(qt);
335
336 Transform3D ntrans;
337 ntrans.SetRotation(M_PI / 3, Point3D(1.0, 0.0, 0.0));
338
339 unsigned int i, j;
340 for (i = 0; i < 4; i++) {
341 for (j = 0; j < 4; j++) {
342 CHECK_INVARIANT(RDKit::feq(trans.getVal(i, j), ntrans.getVal(i, j)), "");
343 }
344 }
345 }
346
main()347 int main() {
348 srand(time(nullptr));
349
350 std::cout << "****************************************\n";
351 std::cout << "testPointND\n";
352 testPointND();
353
354 std::cout << "****************************************\n";
355 std::cout << "testPointOps3D\n";
356 testPointOps3D();
357
358 std::cout << "****************************************\n";
359 std::cout << "testPointOps2D\n";
360 testPointOps2D();
361
362 std::cout << "****************************************\n";
363 std::cout << "test12D\n";
364 test12D();
365
366 std::cout << "****************************************\n";
367 std::cout << "test23D\n";
368 test23D();
369
370 std::cout << "****************************************\n";
371 std::cout << "test3MatMultiply\n";
372 test3MatMultiply();
373
374 std::cout << "****************************************\n";
375 std::cout << "testFromQuaternion\n";
376 testFromQuaternion();
377 std::cout << "****************************************\n";
378 return 0;
379 }
380