MOAB: Mesh Oriented datABase  (version 5.5.0)
affinexform_test.cpp
Go to the documentation of this file.
1 /**
2  * MOAB, a Mesh-Oriented datABase, is a software component for creating,
3  * storing and accessing finite element mesh data.
4  *
5  * Copyright 2004 Sandia Corporation. Under the terms of Contract
6  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
7  * retains certain rights in this software.
8  *
9  * This library is free software; you can redistribute it and/or
10  * modify it under the terms of the GNU Lesser General Public
11  * License as published by the Free Software Foundation; either
12  * version 2.1 of the License, or (at your option) any later version.
13  *
14  */
15 
16 /**
17  * \class AffineXform
18  * \brief Define an affine transformatino
19  * \author Jason Kraftcheck ([email protected])
20  * \date August, 2006
21  */
22 
23 #ifdef _WIN32
24 #define _USE_MATH_DEFINES
25 #endif
26 
27 #include "AffineXform.hpp"
28 #include "moab/Interface.hpp"
29 #include <cassert>
30 
31 /*
32 namespace moab {
33 
34 AffineXform AffineXform::rotation( const double* from_vec, const double* to_vec )
35 {
36  CartVect from(from_vec);
37  CartVect to(to_vec);
38  CartVect a = from * to;
39  double len = a.length();
40 
41  // If input vectors are not parallel (the normal case)
42  if (len >= std::numeric_limits<double>::epsilon()) {
43  from.normalize();
44  to.normalize();
45  return rotation( from % to, (from * to).length(), a/len );
46  }
47 
48  // Vectors are parallel:
49  //
50  // If vectors are in same direction then rotation is identity (no transform)
51  if (from % to >= 0.0)
52  return AffineXform();
53 
54  // Parallel vectors in opposite directions:
55  //
56  // NOTE: This case is ill-defined. There are infinitely
57  // many rotations that can align the two vectors. The angle
58  // of rotation is 180 degrees, but the axis of rotation may
59  // be any unit vector orthogonal to the input vectors.
60  //
61  from.normalize();
62  double lenxy = std::sqrt( from[0]*from[0] + from[1]*from[1] );
63  CartVect axis( -from[0]*from[2]/lenxy,
64  -from[1]*from[2]/lenxy,
65  lenxy );
66  return rotation( -1, 0, axis );
67 }
68 
69 
70 } // namespace moab
71 */
72 
73 using namespace moab;
74 
75 #include <iostream>
76 #define ASSERT_VECTORS_EQUAL( A, B ) assert_vectors_equal( ( A ), ( B ), #A, #B, __LINE__ )
77 #define ASSERT_DOUBLES_EQUAL( A, B ) assert_doubles_equal( ( A ), ( B ), #A, #B, __LINE__ )
78 #define ASSERT( B ) assert_bool( ( B ), #B, __LINE__ )
79 
80 const double TOL = 1e-6;
81 
82 int error_count = 0;
83 
84 void assert_vectors_equal( const double* a, const double* b, const char* sa, const char* sb, int lineno )
85 {
86  if( fabs( a[0] - b[0] ) > TOL || fabs( a[1] - b[1] ) > TOL || fabs( a[2] - b[2] ) > TOL )
87  {
88  std::cerr << "Assertion failed at line " << lineno << std::endl
89  << "\t" << sa << " == " << sb << std::endl
90  << "\t[" << a[0] << ", " << a[1] << ", " << a[2] << "] == [" << b[0] << ", " << b[1] << ", " << b[2]
91  << "]" << std::endl;
92  ++error_count;
93  }
94 }
95 
96 void assert_vectors_equal( const CartVect& a, const CartVect& b, const char* sa, const char* sb, int lineno )
97 {
98  assert_vectors_equal( a.array(), b.array(), sa, sb, lineno );
99 }
100 
101 void assert_doubles_equal( double a, double b, const char* sa, const char* sb, int lineno )
102 {
103  if( fabs( a - b ) > TOL )
104  {
105  std::cerr << "Assertion failed at line " << lineno << std::endl
106  << "\t" << sa << " == " << sb << std::endl
107  << "\t" << a << " == " << b << std::endl;
108  ++error_count;
109  }
110 }
111 
112 void assert_bool( bool b, const char* sb, int lineno )
113 {
114  if( !b )
115  {
116  std::cerr << "Assertion failed at line " << lineno << std::endl << "\t" << sb << std::endl;
117  ++error_count;
118  }
119 }
120 
121 const CartVect point1( 0.0, 0.0, 0.0 ), point2( 3.5, 1000, -200 );
122 const CartVect vect1( 0.0, 0.0, -100.0 ), vect2( 1.0, 0.0, 1.0 );
123 
124 void test_none()
125 {
126  // default xform should do nothing.
127  CartVect output;
128  AffineXform none;
129  none.xform_point( point1.array(), output.array() );
130  ASSERT_VECTORS_EQUAL( output, point1 );
131  none.xform_point( point2.array(), output.array() );
132  ASSERT_VECTORS_EQUAL( output, point2 );
133  none.xform_vector( vect1.array(), output.array() );
134  ASSERT_VECTORS_EQUAL( output, vect1 );
135  none.xform_vector( vect2.array(), output.array() );
136  ASSERT_VECTORS_EQUAL( output, vect2 );
137 }
138 
139 void test_translation()
140 {
141  CartVect offset( 1.0, 2.0, 3.0 );
142  CartVect output;
143 
144  AffineXform move = AffineXform::translation( offset.array() );
145 
146  // test that points are moved by offset
147  move.xform_point( point1.array(), output.array() );
148  ASSERT_VECTORS_EQUAL( output, point1 + offset );
149  move.xform_point( point2.array(), output.array() );
150  ASSERT_VECTORS_EQUAL( output, point2 + offset );
151 
152  // vectors should not be changed by a translation
153  move.xform_vector( vect1.array(), output.array() );
154  ASSERT_VECTORS_EQUAL( output, vect1 );
155  move.xform_vector( vect2.array(), output.array() );
156  ASSERT_VECTORS_EQUAL( output, vect2 );
157 }
158 
159 void test_rotation()
160 {
161  CartVect output;
162 
163  // rotate 90 degress about Z axis
164 
165  AffineXform rot = AffineXform::rotation( M_PI / 2.0, CartVect( 0, 0, 1 ).array() );
166  ASSERT_DOUBLES_EQUAL( rot.matrix().determinant(), 1.0 );
167 
168  rot.xform_point( point1.array(), output.array() );
169  ASSERT_VECTORS_EQUAL( output, point1 ); // origin not affected by transform
170 
171  CartVect expectedz( -point2[1], point2[0], point2[2] ); // in first quadrant
172  rot.xform_point( point2.array(), output.array() );
173  ASSERT_DOUBLES_EQUAL( output.length(), point2.length() );
174  ASSERT_VECTORS_EQUAL( output, expectedz );
175 
176  rot.xform_vector( vect1.array(), output.array() );
177  ASSERT_DOUBLES_EQUAL( output.length(), vect1.length() );
178  ASSERT_VECTORS_EQUAL( output, vect1 );
179 
180  rot.xform_vector( vect2.array(), output.array() );
181  ASSERT_DOUBLES_EQUAL( output.length(), vect2.length() );
182  ASSERT_VECTORS_EQUAL( output, CartVect( 0, 1, 1 ) );
183 
184  // rotate 90 degress about Y axis
185 
186  rot = AffineXform::rotation( M_PI / 2.0, CartVect( 0, 1, 0 ).array() );
187  ASSERT_DOUBLES_EQUAL( rot.matrix().determinant(), 1.0 );
188 
189  rot.xform_point( point1.array(), output.array() );
190  ASSERT_VECTORS_EQUAL( output, point1 ); // origin not affected by transform
191 
192  CartVect expectedy( point2[2], point2[1], -point2[0] ); // in second quadrant
193  rot.xform_point( point2.array(), output.array() );
194  ASSERT_DOUBLES_EQUAL( output.length(), point2.length() );
195  ASSERT_VECTORS_EQUAL( output, expectedy );
196 
197  rot.xform_vector( vect1.array(), output.array() );
198  ASSERT_DOUBLES_EQUAL( output.length(), vect1.length() );
199  ASSERT_VECTORS_EQUAL( output, CartVect( -100, 0, 0 ) );
200 
201  rot.xform_vector( vect2.array(), output.array() );
202  ASSERT_DOUBLES_EQUAL( output.length(), vect2.length() );
203  ASSERT_VECTORS_EQUAL( output, CartVect( 1, 0, -1 ) );
204 
205  // rotate 90 degress about X axis
206 
207  rot = AffineXform::rotation( M_PI / 2.0, CartVect( 1, 0, 0 ).array() );
208  ASSERT_DOUBLES_EQUAL( rot.matrix().determinant(), 1.0 );
209 
210  rot.xform_point( point1.array(), output.array() );
211  ASSERT_VECTORS_EQUAL( output, point1 ); // origin not affected by transform
212 
213  CartVect expectedx( point2[0], -point2[2], point2[1] ); // in third quadrant
214  rot.xform_point( point2.array(), output.array() );
215  ASSERT_DOUBLES_EQUAL( output.length(), point2.length() );
216  ASSERT_VECTORS_EQUAL( output, expectedx );
217 
218  rot.xform_vector( vect1.array(), output.array() );
219  ASSERT_DOUBLES_EQUAL( output.length(), vect1.length() );
220  ASSERT_VECTORS_EQUAL( output, CartVect( 0, 100, 0 ) );
221 
222  rot.xform_vector( vect2.array(), output.array() );
223  ASSERT_DOUBLES_EQUAL( output.length(), vect2.length() );
224  ASSERT_VECTORS_EQUAL( output, CartVect( 1, -1, 0 ) );
225 
226  // rotate 180 degrees about vector in XY plane
227 
228  rot = AffineXform::rotation( M_PI, CartVect( 1, 1, 0 ).array() );
229  ASSERT_DOUBLES_EQUAL( rot.matrix().determinant(), 1.0 );
230 
231  rot.xform_point( point1.array(), output.array() );
232  ASSERT_VECTORS_EQUAL( output, point1 ); // origin not affected by transform
233 
234  rot.xform_point( point2.array(), output.array() );
235  ASSERT_DOUBLES_EQUAL( output.length(), point2.length() );
236  ASSERT_VECTORS_EQUAL( output, CartVect( point2[1], point2[0], -point2[2] ) );
237 
238  rot.xform_vector( vect1.array(), output.array() );
239  ASSERT_DOUBLES_EQUAL( output.length(), vect1.length() );
240  ASSERT_VECTORS_EQUAL( output, -vect1 ); // vector is in xy plane
241 
242  rot.xform_vector( vect2.array(), output.array() );
243  ASSERT_DOUBLES_EQUAL( output.length(), vect2.length() );
244  ASSERT_VECTORS_EQUAL( output, CartVect( 0, 1, -1 ) );
245 }
246 
247 void test_rotation_from_vec()
248 {
249  CartVect v1( 1, 1, 1 );
250  CartVect v2( 1, 0, 0 );
251  AffineXform rot = AffineXform::rotation( v1.array(), v2.array() );
252  CartVect result;
253  rot.xform_vector( v1.array(), result.array() );
254  // vectors should be parallel, but not same length
255  ASSERT_DOUBLES_EQUAL( result.length(), v1.length() );
256  result.normalize();
257  ASSERT_VECTORS_EQUAL( result, v2 );
258 
259  double v3[] = { -1, 0, 0 };
260  rot = AffineXform::rotation( v3, v2.array() );
261  rot.xform_vector( v3, result.array() );
262  ASSERT_VECTORS_EQUAL( result, v2 );
263 }
264 
265 CartVect refl( const CartVect& vect, const CartVect& norm )
266 {
267  CartVect n( norm );
268  n.normalize();
269  double d = vect % n;
270  return vect - 2 * d * n;
271 }
272 
273 void test_reflection()
274 {
275  CartVect output;
276 
277  // reflect about XY plane
278  AffineXform ref = AffineXform::reflection( CartVect( 0, 0, 1 ).array() );
279  ASSERT_DOUBLES_EQUAL( ref.matrix().determinant(), -1.0 );
280  ref.xform_point( point1.array(), output.array() );
281  ASSERT_VECTORS_EQUAL( output, point1 );
282  ref.xform_point( point2.array(), output.array() );
283  ASSERT_DOUBLES_EQUAL( output.length(), point2.length() );
284  ASSERT_VECTORS_EQUAL( output, CartVect( point2[0], point2[1], -point2[2] ) );
285  ref.xform_vector( vect1.array(), output.array() );
286  ASSERT_DOUBLES_EQUAL( output.length(), vect1.length() );
287  ASSERT_VECTORS_EQUAL( output, -vect1 );
288  ref.xform_vector( vect2.array(), output.array() );
289  ASSERT_DOUBLES_EQUAL( output.length(), vect2.length() );
290  ASSERT_VECTORS_EQUAL( output, CartVect( 1, 0, -1 ) );
291 
292  // reflect about arbitrary palne
293  CartVect norm( 3, 2, 1 );
294  ref = AffineXform::reflection( norm.array() );
295  ASSERT_DOUBLES_EQUAL( ref.matrix().determinant(), -1.0 );
296  ref.xform_point( point1.array(), output.array() );
297  ASSERT_VECTORS_EQUAL( output, point1 );
298  ref.xform_point( point2.array(), output.array() );
299  ASSERT_DOUBLES_EQUAL( output.length(), point2.length() );
300  ASSERT_VECTORS_EQUAL( output, refl( point2, norm ) );
301  ref.xform_vector( vect1.array(), output.array() );
302  ASSERT_DOUBLES_EQUAL( output.length(), vect1.length() );
303  ASSERT_VECTORS_EQUAL( output, refl( vect1, norm ) );
304  ref.xform_vector( vect2.array(), output.array() );
305  ASSERT_DOUBLES_EQUAL( output.length(), vect2.length() );
306  ASSERT_VECTORS_EQUAL( output, refl( vect2, norm ) );
307 }
308 
309 void test_scale()
310 {
311  CartVect output;
312 
313  AffineXform scale = AffineXform::scale( 1.0 );
314  ASSERT( !scale.scale() );
315  scale = AffineXform::scale( -1.0 );
316  ASSERT( !scale.scale() );
317 
318  // scale in X only
319  scale = AffineXform::scale( CartVect( 2, 1, 1 ).array() );
320  ASSERT( scale.scale() );
321  scale.xform_point( point1.array(), output.array() );
322  ASSERT_VECTORS_EQUAL( output, CartVect( 2 * point1[0], point1[1], point1[2] ) );
323  scale.xform_point( point2.array(), output.array() );
324  ASSERT_VECTORS_EQUAL( output, CartVect( 2 * point2[0], point2[1], point2[2] ) );
325  scale.xform_vector( vect1.array(), output.array() );
326  ASSERT_VECTORS_EQUAL( output, CartVect( 2 * vect1[0], vect1[1], vect1[2] ) );
327  scale.xform_vector( vect2.array(), output.array() );
328  ASSERT_VECTORS_EQUAL( output, CartVect( 2 * vect2[0], vect2[1], vect2[2] ) );
329 
330  // scale in all
331  scale = AffineXform::scale( CartVect( 0.5, 0.5, 0.5 ).array() );
332  ASSERT( scale.scale() );
333  scale.xform_point( point1.array(), output.array() );
334  ASSERT_VECTORS_EQUAL( output, 0.5 * point1 );
335  scale.xform_point( point2.array(), output.array() );
336  ASSERT_VECTORS_EQUAL( output, 0.5 * point2 );
337  scale.xform_vector( vect1.array(), output.array() );
338  ASSERT_VECTORS_EQUAL( output, 0.5 * vect1 );
339  scale.xform_vector( vect2.array(), output.array() );
340  ASSERT_VECTORS_EQUAL( output, 0.5 * vect2 );
341 }
342 
343 void test_scale_point()
344 {
345  const double point[] = { 2, 3, 4 };
346  const double f[] = { 0.2, 0.1, 0.3 };
347  double result[3];
348  AffineXform scale = AffineXform::scale( f, point );
349  scale.xform_point( point, result );
350  ASSERT_VECTORS_EQUAL( result, point );
351 
352  const double delta[3] = { 1, 0, 2 };
353  const double pt2[] = { point[0] + delta[0], point[1] + delta[1], point[2] + delta[2] };
354  scale = AffineXform::scale( f, point );
355  scale.xform_point( pt2, result );
356 
357  const double expected[] = { point[0] + f[0] * delta[0], point[1] + f[1] * delta[1], point[2] + f[2] * delta[2] };
358  ASSERT_VECTORS_EQUAL( result, expected );
359 }
360 
361 void test_accumulate()
362 {
363  CartVect indiv, accum;
364 
365  // build an group of transforms. make sure translation is somewhere in the middle
366  AffineXform move, scal, rot1, rot2, refl;
367  move = AffineXform::translation( CartVect( 5, -5, 1 ).array() );
368  scal = AffineXform::scale( CartVect( 1, 0.5, 2 ).array() );
369  rot1 = AffineXform::rotation( M_PI / 3, CartVect( 0.5, 0.5, 1 ).array() );
370  rot2 = AffineXform::rotation( M_PI / 4, CartVect( 1.0, 0.0, 0.0 ).array() );
371  refl = AffineXform::reflection( CartVect( -1, -1, 0 ).array() );
372  AffineXform accu;
373  accu.accumulate( scal );
374  accu.accumulate( rot1 );
375  accu.accumulate( move );
376  accu.accumulate( refl );
377  accu.accumulate( rot2 );
378 
379  accu.xform_point( point1.array(), accum.array() );
380  scal.xform_point( point1.array(), indiv.array() );
381  rot1.xform_point( indiv.array() );
382  move.xform_point( indiv.array() );
383  refl.xform_point( indiv.array() );
384  rot2.xform_point( indiv.array() );
385  ASSERT_VECTORS_EQUAL( accum, indiv );
386 
387  accu.xform_point( point2.array(), accum.array() );
388  scal.xform_point( point2.array(), indiv.array() );
389  rot1.xform_point( indiv.array() );
390  move.xform_point( indiv.array() );
391  refl.xform_point( indiv.array() );
392  rot2.xform_point( indiv.array() );
393  ASSERT_VECTORS_EQUAL( accum, indiv );
394 
395  accu.xform_vector( vect1.array(), accum.array() );
396  scal.xform_vector( vect1.array(), indiv.array() );
397  rot1.xform_vector( indiv.array() );
398  move.xform_vector( indiv.array() );
399  refl.xform_vector( indiv.array() );
400  rot2.xform_vector( indiv.array() );
401  ASSERT_VECTORS_EQUAL( accum, indiv );
402 
403  accu.xform_vector( vect2.array(), accum.array() );
404  scal.xform_vector( vect2.array(), indiv.array() );
405  rot1.xform_vector( indiv.array() );
406  move.xform_vector( indiv.array() );
407  refl.xform_vector( indiv.array() );
408  rot2.xform_vector( indiv.array() );
409  ASSERT_VECTORS_EQUAL( accum, indiv );
410 }
411 
412 void test_inversion()
413 {
414  CartVect result;
415 
416  // build an group of transforms. make sure translation is somewhere in the middle
417  AffineXform move, scal, rot1, rot2, refl;
418  move = AffineXform::translation( CartVect( 5, -5, 1 ).array() );
419  scal = AffineXform::scale( CartVect( 1, 0.5, 2 ).array() );
420  rot1 = AffineXform::rotation( M_PI / 3, CartVect( 0.5, 0.5, 1 ).array() );
421  rot2 = AffineXform::rotation( M_PI / 4, CartVect( 1.0, 0.0, 0.0 ).array() );
422  refl = AffineXform::reflection( CartVect( -1, -1, 0 ).array() );
423  AffineXform acc;
424  acc.accumulate( scal );
425  acc.accumulate( rot1 );
426  acc.accumulate( move );
427  acc.accumulate( refl );
428  acc.accumulate( rot2 );
429 
430  AffineXform inv = acc.inverse();
431 
432  acc.xform_point( point1.array(), result.array() );
433  inv.xform_point( result.array() );
434  ASSERT_VECTORS_EQUAL( point1, result );
435 
436  acc.xform_point( point2.array(), result.array() );
437  inv.xform_point( result.array() );
438  ASSERT_VECTORS_EQUAL( point2, result );
439 
440  acc.xform_vector( vect1.array(), result.array() );
441  inv.xform_vector( result.array() );
442  ASSERT_VECTORS_EQUAL( vect1, result );
443 
444  acc.xform_vector( vect2.array(), result.array() );
445  inv.xform_vector( result.array() );
446  ASSERT_VECTORS_EQUAL( vect2, result );
447 }
448 
449 void test_is_reflection()
450 {
451  AffineXform refl1, refl2, scale;
452  refl1 = AffineXform::reflection( CartVect( -1, -1, 0 ).array() );
453  refl2 = AffineXform::reflection( CartVect( 1, 0, 0 ).array() );
454  scale = AffineXform::scale( CartVect( -1, 1, 1 ).array() );
455 
456  ASSERT( refl1.reflection() );
457  ASSERT( refl2.reflection() );
458  ASSERT( scale.reflection() );
459 
460  AffineXform inv1, inv2, inv3;
461  inv1 = refl1.inverse();
462  inv2 = refl2.inverse();
463  inv3 = scale.inverse();
464 
465  ASSERT( inv1.reflection() );
466  ASSERT( inv2.reflection() );
467  ASSERT( inv3.reflection() );
468 
469  refl1.accumulate( refl2 );
470  refl2.accumulate( scale );
471  ASSERT( !refl1.reflection() );
472  ASSERT( !refl2.reflection() );
473 
474  AffineXform rot, mov;
475  rot = AffineXform::rotation( M_PI / 4, CartVect( 1, 1, 1 ).array() );
476  mov = AffineXform::translation( CartVect( -5, 6, 7 ).array() );
477  ASSERT( !rot.reflection() );
478  ASSERT( !mov.reflection() );
479 }
480 
481 int main()
482 {
483  test_none();
484  test_translation();
485  test_rotation();
486  test_reflection();
487  test_rotation_from_vec();
488  test_scale();
489  test_scale_point();
490  test_accumulate();
491  test_inversion();
492  test_is_reflection();
493  return error_count;
494 }