MOAB: Mesh Oriented datABase  (version 5.4.0)
HomXform.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 #include "moab/HomXform.hpp"
00017 #include <cassert>
00018 
00019 namespace moab
00020 {
00021 
00022 HomCoord HomCoord::unitv[3] = { HomCoord( 1, 0, 0 ), HomCoord( 0, 1, 0 ), HomCoord( 0, 0, 1 ) };
00023 HomCoord HomCoord::IDENTITY( 1, 1, 1 );
00024 
00025 HomCoord& HomCoord::getUnitv( int c )
00026 {
00027     return unitv[c];
00028 }
00029 
00030 int dum[] = { 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1 };
00031 HomXform HomXform::IDENTITY( dum );
00032 
00033 void HomXform::three_pt_xform( const HomCoord& p1,
00034                                const HomCoord& q1,
00035                                const HomCoord& p2,
00036                                const HomCoord& q2,
00037                                const HomCoord& p3,
00038                                const HomCoord& q3 )
00039 {
00040     // pmin and pmax are min and max bounding box corners which are mapped to
00041     // qmin and qmax, resp.  qmin and qmax are not necessarily min/max corners,
00042     // since the mapping can change the orientation of the box in the q reference
00043     // frame.  Re-interpreting the min/max bounding box corners does not change
00044     // the mapping.
00045 
00046     // change that: base on three points for now (figure out whether we can
00047     // just use two later); three points are assumed to define an orthogonal
00048     // system such that (p2-p1)%(p3-p1) = 0
00049 
00050     // use the three point rule to compute the mapping, from Mortensen,
00051     // "Geometric Modeling".  If p1, p2, p3 and q1, q2, q3 are three points in
00052     // the two coordinate systems, the three pt rule is:
00053     //
00054     // v1 = p2 - p1
00055     // v2 = p3 - p1
00056     // v3 = v1 x v2
00057     // w1-w3 similar, with q1-q3
00058     // V = matrix with v1-v3 as rows
00059     // W similar, with w1-w3
00060     // R = V^-1 * W
00061     // t = q1 - p1 * W
00062     // Form transform matrix M from R, t
00063 
00064     // check to see whether unity transform applies
00065     if( p1 == q1 && p2 == q2 && p3 == q3 )
00066     {
00067         *this = HomXform::IDENTITY;
00068         return;
00069     }
00070 
00071     // first, construct 3 pts from input
00072     HomCoord v1 = p2 - p1;
00073     assert( v1.i() != 0 || v1.j() != 0 || v1.k() != 0 );
00074     HomCoord v2 = p3 - p1;
00075     HomCoord v3 = v1 * v2;
00076 
00077     if( v3.length_squared() == 0 )
00078     {
00079         // 1d coordinate system; set one of v2's coordinates such that
00080         // it's orthogonal to v1
00081         if( v1.i() == 0 )
00082             v2.set( 1, 0, 0 );
00083         else if( v1.j() == 0 )
00084             v2.set( 0, 1, 0 );
00085         else if( v1.k() == 0 )
00086             v2.set( 0, 0, 1 );
00087         else
00088             assert( false );
00089         v3 = v1 * v2;
00090         assert( v3.length_squared() != 0 );
00091     }
00092     // assert to make sure they're each orthogonal
00093     assert( v1 % v2 == 0 && v1 % v3 == 0 && v2 % v3 == 0 );
00094     v1.normalize();
00095     v2.normalize();
00096     v3.normalize();
00097     // Make sure h is set to zero here, since it'll mess up things if it's one
00098     v1.homCoord[3] = v2.homCoord[3] = v3.homCoord[3] = 0;
00099 
00100     HomCoord w1 = q2 - q1;
00101     assert( w1.i() != 0 || w1.j() != 0 || w1.k() != 0 );
00102     HomCoord w2 = q3 - q1;
00103     HomCoord w3 = w1 * w2;
00104 
00105     if( w3.length_squared() == 0 )
00106     {
00107         // 1d coordinate system; set one of w2's coordinates such that
00108         // it's orthogonal to w1
00109         if( w1.i() == 0 )
00110             w2.set( 1, 0, 0 );
00111         else if( w1.j() == 0 )
00112             w2.set( 0, 1, 0 );
00113         else if( w1.k() == 0 )
00114             w2.set( 0, 0, 1 );
00115         else
00116             assert( false );
00117         w3 = w1 * w2;
00118         assert( w3.length_squared() != 0 );
00119     }
00120     // assert to make sure they're each orthogonal
00121     assert( w1 % w2 == 0 && w1 % w3 == 0 && w2 % w3 == 0 );
00122     w1.normalize();
00123     w2.normalize();
00124     w3.normalize();
00125     // Make sure h is set to zero here, since it'll mess up things if it's one
00126     w1.homCoord[3] = w2.homCoord[3] = w3.homCoord[3] = 0;
00127 
00128     // form v^-1 as transpose (ok for orthogonal vectors); put directly into
00129     // transform matrix, since it's eventually going to become R
00130     *this = HomXform( v1.i(), v2.i(), v3.i(), 0, v1.j(), v2.j(), v3.j(), 0, v1.k(), v2.k(), v3.k(), 0, 0, 0, 0, 1 );
00131 
00132     // multiply by w to get R
00133     *this *= HomXform( w1.i(), w1.j(), w1.k(), 0, w2.i(), w2.j(), w2.k(), 0, w3.i(), w3.j(), w3.k(), 0, 0, 0, 0, 1 );
00134 
00135     // compute t and put into last row
00136     HomCoord t              = q1 - p1 * *this;
00137     ( *this ).XFORM( 3, 0 ) = t.i();
00138     ( *this ).XFORM( 3, 1 ) = t.j();
00139     ( *this ).XFORM( 3, 2 ) = t.k();
00140 
00141     // M should transform p to q
00142     assert( ( q1 == p1 * *this ) && ( q2 == p2 * *this ) && ( q3 == p3 * *this ) );
00143 }
00144 
00145 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines