MOAB: Mesh Oriented datABase  (version 5.3.0)
affinexform_test.cpp
Go to the documentation of this file.
00001 /**
00002  * MOAB, a Mesh-Oriented datABase, is a software component for creating,
00003  * storing and accessing finite element mesh data.
00004  *
00005  * Copyright 2004 Sandia Corporation.  Under the terms of Contract
00006  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
00007  * retains certain rights in this software.
00008  *
00009  * This library is free software; you can redistribute it and/or
00010  * modify it under the terms of the GNU Lesser General Public
00011  * License as published by the Free Software Foundation; either
00012  * version 2.1 of the License, or (at your option) any later version.
00013  *
00014  */
00015 
00016 /**
00017  * \class AffineXform
00018  * \brief Define an affine transformatino
00019  * \author Jason Kraftcheck (kraftche@cae.wisc.edu)
00020  * \date August, 2006
00021  */
00022 
00023 #ifdef _WIN32
00024 #define _USE_MATH_DEFINES
00025 #endif
00026 
00027 #include "AffineXform.hpp"
00028 #include "moab/Interface.hpp"
00029 #include <cassert>
00030 
00031 /*
00032 namespace moab {
00033 
00034 AffineXform AffineXform::rotation( const double* from_vec, const double* to_vec )
00035 {
00036   CartVect from(from_vec);
00037   CartVect to(to_vec);
00038   CartVect a = from * to;
00039   double len = a.length();
00040 
00041   // If input vectors are not parallel (the normal case)
00042   if (len >= std::numeric_limits<double>::epsilon()) {
00043     from.normalize();
00044     to.normalize();
00045     return rotation( from % to, (from * to).length(), a/len );
00046   }
00047 
00048   // Vectors are parallel:
00049   //
00050   // If vectors are in same direction then rotation is identity (no transform)
00051   if (from % to >= 0.0)
00052     return AffineXform();
00053 
00054   // Parallel vectors in opposite directions:
00055   //
00056   // NOTE:  This case is ill-defined.  There are infinitely
00057   // many rotations that can align the two vectors.  The angle
00058   // of rotation is 180 degrees, but the axis of rotation may
00059   // be any unit vector orthogonal to the input vectors.
00060   //
00061   from.normalize();
00062   double lenxy = std::sqrt( from[0]*from[0] + from[1]*from[1] );
00063   CartVect axis( -from[0]*from[2]/lenxy,
00064                  -from[1]*from[2]/lenxy,
00065                                   lenxy );
00066   return rotation( -1, 0, axis );
00067 }
00068 
00069 
00070 } // namespace moab
00071 */
00072 
00073 using namespace moab;
00074 
00075 #include <iostream>
00076 #define ASSERT_VECTORS_EQUAL( A, B ) assert_vectors_equal( ( A ), ( B ), #A, #B, __LINE__ )
00077 #define ASSERT_DOUBLES_EQUAL( A, B ) assert_doubles_equal( ( A ), ( B ), #A, #B, __LINE__ )
00078 #define ASSERT( B )                  assert_bool( ( B ), #B, __LINE__ )
00079 
00080 const double TOL = 1e-6;
00081 
00082 int error_count = 0;
00083 
00084 void assert_vectors_equal( const double* a, const double* b, const char* sa, const char* sb, int lineno )
00085 {
00086     if( fabs( a[0] - b[0] ) > TOL || fabs( a[1] - b[1] ) > TOL || fabs( a[2] - b[2] ) > TOL )
00087     {
00088         std::cerr << "Assertion failed at line " << lineno << std::endl
00089                   << "\t" << sa << " == " << sb << std::endl
00090                   << "\t[" << a[0] << ", " << a[1] << ", " << a[2] << "] == [" << b[0] << ", " << b[1] << ", " << b[2]
00091                   << "]" << std::endl;
00092         ++error_count;
00093     }
00094 }
00095 
00096 void assert_vectors_equal( const CartVect& a, const CartVect& b, const char* sa, const char* sb, int lineno )
00097 {
00098     assert_vectors_equal( a.array(), b.array(), sa, sb, lineno );
00099 }
00100 
00101 void assert_doubles_equal( double a, double b, const char* sa, const char* sb, int lineno )
00102 {
00103     if( fabs( a - b ) > TOL )
00104     {
00105         std::cerr << "Assertion failed at line " << lineno << std::endl
00106                   << "\t" << sa << " == " << sb << std::endl
00107                   << "\t" << a << " == " << b << std::endl;
00108         ++error_count;
00109     }
00110 }
00111 
00112 void assert_bool( bool b, const char* sb, int lineno )
00113 {
00114     if( !b )
00115     {
00116         std::cerr << "Assertion failed at line " << lineno << std::endl << "\t" << sb << std::endl;
00117         ++error_count;
00118     }
00119 }
00120 
00121 const CartVect point1( 0.0, 0.0, 0.0 ), point2( 3.5, 1000, -200 );
00122 const CartVect vect1( 0.0, 0.0, -100.0 ), vect2( 1.0, 0.0, 1.0 );
00123 
00124 void test_none()
00125 {
00126     // default xform should do nothing.
00127     CartVect output;
00128     AffineXform none;
00129     none.xform_point( point1.array(), output.array() );
00130     ASSERT_VECTORS_EQUAL( output, point1 );
00131     none.xform_point( point2.array(), output.array() );
00132     ASSERT_VECTORS_EQUAL( output, point2 );
00133     none.xform_vector( vect1.array(), output.array() );
00134     ASSERT_VECTORS_EQUAL( output, vect1 );
00135     none.xform_vector( vect2.array(), output.array() );
00136     ASSERT_VECTORS_EQUAL( output, vect2 );
00137 }
00138 
00139 void test_translation()
00140 {
00141     CartVect offset( 1.0, 2.0, 3.0 );
00142     CartVect output;
00143 
00144     AffineXform move = AffineXform::translation( offset.array() );
00145 
00146     // test that points are moved by offset
00147     move.xform_point( point1.array(), output.array() );
00148     ASSERT_VECTORS_EQUAL( output, point1 + offset );
00149     move.xform_point( point2.array(), output.array() );
00150     ASSERT_VECTORS_EQUAL( output, point2 + offset );
00151 
00152     // vectors should not be changed by a translation
00153     move.xform_vector( vect1.array(), output.array() );
00154     ASSERT_VECTORS_EQUAL( output, vect1 );
00155     move.xform_vector( vect2.array(), output.array() );
00156     ASSERT_VECTORS_EQUAL( output, vect2 );
00157 }
00158 
00159 void test_rotation()
00160 {
00161     CartVect output;
00162 
00163     // rotate 90 degress about Z axis
00164 
00165     AffineXform rot = AffineXform::rotation( M_PI / 2.0, CartVect( 0, 0, 1 ).array() );
00166     ASSERT_DOUBLES_EQUAL( rot.matrix().determinant(), 1.0 );
00167 
00168     rot.xform_point( point1.array(), output.array() );
00169     ASSERT_VECTORS_EQUAL( output, point1 );  // origin not affected by transform
00170 
00171     CartVect expectedz( -point2[1], point2[0], point2[2] );  // in first quadrant
00172     rot.xform_point( point2.array(), output.array() );
00173     ASSERT_DOUBLES_EQUAL( output.length(), point2.length() );
00174     ASSERT_VECTORS_EQUAL( output, expectedz );
00175 
00176     rot.xform_vector( vect1.array(), output.array() );
00177     ASSERT_DOUBLES_EQUAL( output.length(), vect1.length() );
00178     ASSERT_VECTORS_EQUAL( output, vect1 );
00179 
00180     rot.xform_vector( vect2.array(), output.array() );
00181     ASSERT_DOUBLES_EQUAL( output.length(), vect2.length() );
00182     ASSERT_VECTORS_EQUAL( output, CartVect( 0, 1, 1 ) );
00183 
00184     // rotate 90 degress about Y axis
00185 
00186     rot = AffineXform::rotation( M_PI / 2.0, CartVect( 0, 1, 0 ).array() );
00187     ASSERT_DOUBLES_EQUAL( rot.matrix().determinant(), 1.0 );
00188 
00189     rot.xform_point( point1.array(), output.array() );
00190     ASSERT_VECTORS_EQUAL( output, point1 );  // origin not affected by transform
00191 
00192     CartVect expectedy( point2[2], point2[1], -point2[0] );  // in second quadrant
00193     rot.xform_point( point2.array(), output.array() );
00194     ASSERT_DOUBLES_EQUAL( output.length(), point2.length() );
00195     ASSERT_VECTORS_EQUAL( output, expectedy );
00196 
00197     rot.xform_vector( vect1.array(), output.array() );
00198     ASSERT_DOUBLES_EQUAL( output.length(), vect1.length() );
00199     ASSERT_VECTORS_EQUAL( output, CartVect( -100, 0, 0 ) );
00200 
00201     rot.xform_vector( vect2.array(), output.array() );
00202     ASSERT_DOUBLES_EQUAL( output.length(), vect2.length() );
00203     ASSERT_VECTORS_EQUAL( output, CartVect( 1, 0, -1 ) );
00204 
00205     // rotate 90 degress about X axis
00206 
00207     rot = AffineXform::rotation( M_PI / 2.0, CartVect( 1, 0, 0 ).array() );
00208     ASSERT_DOUBLES_EQUAL( rot.matrix().determinant(), 1.0 );
00209 
00210     rot.xform_point( point1.array(), output.array() );
00211     ASSERT_VECTORS_EQUAL( output, point1 );  // origin not affected by transform
00212 
00213     CartVect expectedx( point2[0], -point2[2], point2[1] );  // in third quadrant
00214     rot.xform_point( point2.array(), output.array() );
00215     ASSERT_DOUBLES_EQUAL( output.length(), point2.length() );
00216     ASSERT_VECTORS_EQUAL( output, expectedx );
00217 
00218     rot.xform_vector( vect1.array(), output.array() );
00219     ASSERT_DOUBLES_EQUAL( output.length(), vect1.length() );
00220     ASSERT_VECTORS_EQUAL( output, CartVect( 0, 100, 0 ) );
00221 
00222     rot.xform_vector( vect2.array(), output.array() );
00223     ASSERT_DOUBLES_EQUAL( output.length(), vect2.length() );
00224     ASSERT_VECTORS_EQUAL( output, CartVect( 1, -1, 0 ) );
00225 
00226     // rotate 180 degrees about vector in XY plane
00227 
00228     rot = AffineXform::rotation( M_PI, CartVect( 1, 1, 0 ).array() );
00229     ASSERT_DOUBLES_EQUAL( rot.matrix().determinant(), 1.0 );
00230 
00231     rot.xform_point( point1.array(), output.array() );
00232     ASSERT_VECTORS_EQUAL( output, point1 );  // origin not affected by transform
00233 
00234     rot.xform_point( point2.array(), output.array() );
00235     ASSERT_DOUBLES_EQUAL( output.length(), point2.length() );
00236     ASSERT_VECTORS_EQUAL( output, CartVect( point2[1], point2[0], -point2[2] ) );
00237 
00238     rot.xform_vector( vect1.array(), output.array() );
00239     ASSERT_DOUBLES_EQUAL( output.length(), vect1.length() );
00240     ASSERT_VECTORS_EQUAL( output, -vect1 );  // vector is in xy plane
00241 
00242     rot.xform_vector( vect2.array(), output.array() );
00243     ASSERT_DOUBLES_EQUAL( output.length(), vect2.length() );
00244     ASSERT_VECTORS_EQUAL( output, CartVect( 0, 1, -1 ) );
00245 }
00246 
00247 void test_rotation_from_vec()
00248 {
00249     CartVect v1( 1, 1, 1 );
00250     CartVect v2( 1, 0, 0 );
00251     AffineXform rot = AffineXform::rotation( v1.array(), v2.array() );
00252     CartVect result;
00253     rot.xform_vector( v1.array(), result.array() );
00254     // vectors should be parallel, but not same length
00255     ASSERT_DOUBLES_EQUAL( result.length(), v1.length() );
00256     result.normalize();
00257     ASSERT_VECTORS_EQUAL( result, v2 );
00258 
00259     double v3[] = { -1, 0, 0 };
00260     rot         = AffineXform::rotation( v3, v2.array() );
00261     rot.xform_vector( v3, result.array() );
00262     ASSERT_VECTORS_EQUAL( result, v2 );
00263 }
00264 
00265 CartVect refl( const CartVect& vect, const CartVect& norm )
00266 {
00267     CartVect n( norm );
00268     n.normalize();
00269     double d = vect % n;
00270     return vect - 2 * d * n;
00271 }
00272 
00273 void test_reflection()
00274 {
00275     CartVect output;
00276 
00277     // reflect about XY plane
00278     AffineXform ref = AffineXform::reflection( CartVect( 0, 0, 1 ).array() );
00279     ASSERT_DOUBLES_EQUAL( ref.matrix().determinant(), -1.0 );
00280     ref.xform_point( point1.array(), output.array() );
00281     ASSERT_VECTORS_EQUAL( output, point1 );
00282     ref.xform_point( point2.array(), output.array() );
00283     ASSERT_DOUBLES_EQUAL( output.length(), point2.length() );
00284     ASSERT_VECTORS_EQUAL( output, CartVect( point2[0], point2[1], -point2[2] ) );
00285     ref.xform_vector( vect1.array(), output.array() );
00286     ASSERT_DOUBLES_EQUAL( output.length(), vect1.length() );
00287     ASSERT_VECTORS_EQUAL( output, -vect1 );
00288     ref.xform_vector( vect2.array(), output.array() );
00289     ASSERT_DOUBLES_EQUAL( output.length(), vect2.length() );
00290     ASSERT_VECTORS_EQUAL( output, CartVect( 1, 0, -1 ) );
00291 
00292     // reflect about arbitrary palne
00293     CartVect norm( 3, 2, 1 );
00294     ref = AffineXform::reflection( norm.array() );
00295     ASSERT_DOUBLES_EQUAL( ref.matrix().determinant(), -1.0 );
00296     ref.xform_point( point1.array(), output.array() );
00297     ASSERT_VECTORS_EQUAL( output, point1 );
00298     ref.xform_point( point2.array(), output.array() );
00299     ASSERT_DOUBLES_EQUAL( output.length(), point2.length() );
00300     ASSERT_VECTORS_EQUAL( output, refl( point2, norm ) );
00301     ref.xform_vector( vect1.array(), output.array() );
00302     ASSERT_DOUBLES_EQUAL( output.length(), vect1.length() );
00303     ASSERT_VECTORS_EQUAL( output, refl( vect1, norm ) );
00304     ref.xform_vector( vect2.array(), output.array() );
00305     ASSERT_DOUBLES_EQUAL( output.length(), vect2.length() );
00306     ASSERT_VECTORS_EQUAL( output, refl( vect2, norm ) );
00307 }
00308 
00309 void test_scale()
00310 {
00311     CartVect output;
00312 
00313     AffineXform scale = AffineXform::scale( 1.0 );
00314     ASSERT( !scale.scale() );
00315     scale = AffineXform::scale( -1.0 );
00316     ASSERT( !scale.scale() );
00317 
00318     // scale in X only
00319     scale = AffineXform::scale( CartVect( 2, 1, 1 ).array() );
00320     ASSERT( scale.scale() );
00321     scale.xform_point( point1.array(), output.array() );
00322     ASSERT_VECTORS_EQUAL( output, CartVect( 2 * point1[0], point1[1], point1[2] ) );
00323     scale.xform_point( point2.array(), output.array() );
00324     ASSERT_VECTORS_EQUAL( output, CartVect( 2 * point2[0], point2[1], point2[2] ) );
00325     scale.xform_vector( vect1.array(), output.array() );
00326     ASSERT_VECTORS_EQUAL( output, CartVect( 2 * vect1[0], vect1[1], vect1[2] ) );
00327     scale.xform_vector( vect2.array(), output.array() );
00328     ASSERT_VECTORS_EQUAL( output, CartVect( 2 * vect2[0], vect2[1], vect2[2] ) );
00329 
00330     // scale in all
00331     scale = AffineXform::scale( CartVect( 0.5, 0.5, 0.5 ).array() );
00332     ASSERT( scale.scale() );
00333     scale.xform_point( point1.array(), output.array() );
00334     ASSERT_VECTORS_EQUAL( output, 0.5 * point1 );
00335     scale.xform_point( point2.array(), output.array() );
00336     ASSERT_VECTORS_EQUAL( output, 0.5 * point2 );
00337     scale.xform_vector( vect1.array(), output.array() );
00338     ASSERT_VECTORS_EQUAL( output, 0.5 * vect1 );
00339     scale.xform_vector( vect2.array(), output.array() );
00340     ASSERT_VECTORS_EQUAL( output, 0.5 * vect2 );
00341 }
00342 
00343 void test_scale_point()
00344 {
00345     const double point[] = { 2, 3, 4 };
00346     const double f[]     = { 0.2, 0.1, 0.3 };
00347     double result[3];
00348     AffineXform scale = AffineXform::scale( f, point );
00349     scale.xform_point( point, result );
00350     ASSERT_VECTORS_EQUAL( result, point );
00351 
00352     const double delta[3] = { 1, 0, 2 };
00353     const double pt2[]    = { point[0] + delta[0], point[1] + delta[1], point[2] + delta[2] };
00354     scale                 = AffineXform::scale( f, point );
00355     scale.xform_point( pt2, result );
00356 
00357     const double expected[] = { point[0] + f[0] * delta[0], point[1] + f[1] * delta[1], point[2] + f[2] * delta[2] };
00358     ASSERT_VECTORS_EQUAL( result, expected );
00359 }
00360 
00361 void test_accumulate()
00362 {
00363     CartVect indiv, accum;
00364 
00365     // build an group of transforms.  make sure translation is somewhere in the middle
00366     AffineXform move, scal, rot1, rot2, refl;
00367     move = AffineXform::translation( CartVect( 5, -5, 1 ).array() );
00368     scal = AffineXform::scale( CartVect( 1, 0.5, 2 ).array() );
00369     rot1 = AffineXform::rotation( M_PI / 3, CartVect( 0.5, 0.5, 1 ).array() );
00370     rot2 = AffineXform::rotation( M_PI / 4, CartVect( 1.0, 0.0, 0.0 ).array() );
00371     refl = AffineXform::reflection( CartVect( -1, -1, 0 ).array() );
00372     AffineXform accu;
00373     accu.accumulate( scal );
00374     accu.accumulate( rot1 );
00375     accu.accumulate( move );
00376     accu.accumulate( refl );
00377     accu.accumulate( rot2 );
00378 
00379     accu.xform_point( point1.array(), accum.array() );
00380     scal.xform_point( point1.array(), indiv.array() );
00381     rot1.xform_point( indiv.array() );
00382     move.xform_point( indiv.array() );
00383     refl.xform_point( indiv.array() );
00384     rot2.xform_point( indiv.array() );
00385     ASSERT_VECTORS_EQUAL( accum, indiv );
00386 
00387     accu.xform_point( point2.array(), accum.array() );
00388     scal.xform_point( point2.array(), indiv.array() );
00389     rot1.xform_point( indiv.array() );
00390     move.xform_point( indiv.array() );
00391     refl.xform_point( indiv.array() );
00392     rot2.xform_point( indiv.array() );
00393     ASSERT_VECTORS_EQUAL( accum, indiv );
00394 
00395     accu.xform_vector( vect1.array(), accum.array() );
00396     scal.xform_vector( vect1.array(), indiv.array() );
00397     rot1.xform_vector( indiv.array() );
00398     move.xform_vector( indiv.array() );
00399     refl.xform_vector( indiv.array() );
00400     rot2.xform_vector( indiv.array() );
00401     ASSERT_VECTORS_EQUAL( accum, indiv );
00402 
00403     accu.xform_vector( vect2.array(), accum.array() );
00404     scal.xform_vector( vect2.array(), indiv.array() );
00405     rot1.xform_vector( indiv.array() );
00406     move.xform_vector( indiv.array() );
00407     refl.xform_vector( indiv.array() );
00408     rot2.xform_vector( indiv.array() );
00409     ASSERT_VECTORS_EQUAL( accum, indiv );
00410 }
00411 
00412 void test_inversion()
00413 {
00414     CartVect result;
00415 
00416     // build an group of transforms.  make sure translation is somewhere in the middle
00417     AffineXform move, scal, rot1, rot2, refl;
00418     move = AffineXform::translation( CartVect( 5, -5, 1 ).array() );
00419     scal = AffineXform::scale( CartVect( 1, 0.5, 2 ).array() );
00420     rot1 = AffineXform::rotation( M_PI / 3, CartVect( 0.5, 0.5, 1 ).array() );
00421     rot2 = AffineXform::rotation( M_PI / 4, CartVect( 1.0, 0.0, 0.0 ).array() );
00422     refl = AffineXform::reflection( CartVect( -1, -1, 0 ).array() );
00423     AffineXform acc;
00424     acc.accumulate( scal );
00425     acc.accumulate( rot1 );
00426     acc.accumulate( move );
00427     acc.accumulate( refl );
00428     acc.accumulate( rot2 );
00429 
00430     AffineXform inv = acc.inverse();
00431 
00432     acc.xform_point( point1.array(), result.array() );
00433     inv.xform_point( result.array() );
00434     ASSERT_VECTORS_EQUAL( point1, result );
00435 
00436     acc.xform_point( point2.array(), result.array() );
00437     inv.xform_point( result.array() );
00438     ASSERT_VECTORS_EQUAL( point2, result );
00439 
00440     acc.xform_vector( vect1.array(), result.array() );
00441     inv.xform_vector( result.array() );
00442     ASSERT_VECTORS_EQUAL( vect1, result );
00443 
00444     acc.xform_vector( vect2.array(), result.array() );
00445     inv.xform_vector( result.array() );
00446     ASSERT_VECTORS_EQUAL( vect2, result );
00447 }
00448 
00449 void test_is_reflection()
00450 {
00451     AffineXform refl1, refl2, scale;
00452     refl1 = AffineXform::reflection( CartVect( -1, -1, 0 ).array() );
00453     refl2 = AffineXform::reflection( CartVect( 1, 0, 0 ).array() );
00454     scale = AffineXform::scale( CartVect( -1, 1, 1 ).array() );
00455 
00456     ASSERT( refl1.reflection() );
00457     ASSERT( refl2.reflection() );
00458     ASSERT( scale.reflection() );
00459 
00460     AffineXform inv1, inv2, inv3;
00461     inv1 = refl1.inverse();
00462     inv2 = refl2.inverse();
00463     inv3 = scale.inverse();
00464 
00465     ASSERT( inv1.reflection() );
00466     ASSERT( inv2.reflection() );
00467     ASSERT( inv3.reflection() );
00468 
00469     refl1.accumulate( refl2 );
00470     refl2.accumulate( scale );
00471     ASSERT( !refl1.reflection() );
00472     ASSERT( !refl2.reflection() );
00473 
00474     AffineXform rot, mov;
00475     rot = AffineXform::rotation( M_PI / 4, CartVect( 1, 1, 1 ).array() );
00476     mov = AffineXform::translation( CartVect( -5, 6, 7 ).array() );
00477     ASSERT( !rot.reflection() );
00478     ASSERT( !mov.reflection() );
00479 }
00480 
00481 int main()
00482 {
00483     test_none();
00484     test_translation();
00485     test_rotation();
00486     test_reflection();
00487     test_rotation_from_vec();
00488     test_scale();
00489     test_scale_point();
00490     test_accumulate();
00491     test_inversion();
00492     test_is_reflection();
00493     return error_count;
00494 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines