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