HomXform.cpp

Go to the documentation of this file.

/**
 * MOAB, a Mesh-Oriented datABase, is a software component for creating,
 * storing and accessing finite element mesh data.
 *
 * Copyright 2004 Sandia Corporation. Under the terms of Contract
 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
 * retains certain rights in this software.
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 */
 
#include "moab/HomXform.hpp"
#include <cassert>
 
namespace moab
{
 
HomCoord HomCoord::unitv[3] = { HomCoord( 1, 0, 0 ), HomCoord( 0, 1, 0 ), HomCoord( 0, 0, 1 ) };
HomCoord HomCoord::IDENTITY( 1, 1, 1 );
 
HomCoord& HomCoord::getUnitv( int c )
{
 return unitv[c];
}
 
int dum[] = { 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1 };
HomXform HomXform::IDENTITY( dum );
 
void HomXform::three_pt_xform( const HomCoord& p1,
 const HomCoord& q1,
 const HomCoord& p2,
 const HomCoord& q2,
 const HomCoord& p3,
 const HomCoord& q3 )
{
 // pmin and pmax are min and max bounding box corners which are mapped to
 // qmin and qmax, resp. qmin and qmax are not necessarily min/max corners,
 // since the mapping can change the orientation of the box in the q reference
 // frame. Re-interpreting the min/max bounding box corners does not change
 // the mapping.
 
 // change that: base on three points for now (figure out whether we can
 // just use two later); three points are assumed to define an orthogonal
 // system such that (p2-p1)%(p3-p1) = 0
 
 // use the three point rule to compute the mapping, from Mortensen,
 // "Geometric Modeling". If p1, p2, p3 and q1, q2, q3 are three points in
 // the two coordinate systems, the three pt rule is:
 //
 // v1 = p2 - p1
 // v2 = p3 - p1
 // v3 = v1 x v2
 // w1-w3 similar, with q1-q3
 // V = matrix with v1-v3 as rows
 // W similar, with w1-w3
 // R = V^-1 * W
 // t = q1 - p1 * W
 // Form transform matrix M from R, t
 
 // check to see whether unity transform applies
 if( p1 == q1 && p2 == q2 && p3 == q3 )
 {
 *this = HomXform::IDENTITY;
 return;
 }
 
 // first, construct 3 pts from input
 HomCoord v1 = p2 - p1;
 assert( v1.i() != 0 || v1.j() != 0 || v1.k() != 0 );
 HomCoord v2 = p3 - p1;
 HomCoord v3 = v1 * v2;
 
 if( v3.length_squared() == 0 )
 {
 // 1d coordinate system; set one of v2's coordinates such that
 // it's orthogonal to v1
 if( v1.i() == 0 )
 v2.set( 1, 0, 0 );
 else if( v1.j() == 0 )
 v2.set( 0, 1, 0 );
 else if( v1.k() == 0 )
 v2.set( 0, 0, 1 );
 else
 assert( false );
 v3 = v1 * v2;
 assert( v3.length_squared() != 0 );
 }
 // assert to make sure they're each orthogonal
 assert( v1 % v2 == 0 && v1 % v3 == 0 && v2 % v3 == 0 );
 v1.normalize();
 v2.normalize();
 v3.normalize();
 // Make sure h is set to zero here, since it'll mess up things if it's one
 v1.homCoord[3] = v2.homCoord[3] = v3.homCoord[3] = 0;
 
 HomCoord w1 = q2 - q1;
 assert( w1.i() != 0 || w1.j() != 0 || w1.k() != 0 );
 HomCoord w2 = q3 - q1;
 HomCoord w3 = w1 * w2;
 
 if( w3.length_squared() == 0 )
 {
 // 1d coordinate system; set one of w2's coordinates such that
 // it's orthogonal to w1
 if( w1.i() == 0 )
 w2.set( 1, 0, 0 );
 else if( w1.j() == 0 )
 w2.set( 0, 1, 0 );
 else if( w1.k() == 0 )
 w2.set( 0, 0, 1 );
 else
 assert( false );
 w3 = w1 * w2;
 assert( w3.length_squared() != 0 );
 }
 // assert to make sure they're each orthogonal
 assert( w1 % w2 == 0 && w1 % w3 == 0 && w2 % w3 == 0 );
 w1.normalize();
 w2.normalize();
 w3.normalize();
 // Make sure h is set to zero here, since it'll mess up things if it's one
 w1.homCoord[3] = w2.homCoord[3] = w3.homCoord[3] = 0;
 
 // form v^-1 as transpose (ok for orthogonal vectors); put directly into
 // transform matrix, since it's eventually going to become R
 *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 );
 
 // multiply by w to get R
 *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 );
 
 // compute t and put into last row
 HomCoord t = q1 - p1 * *this;
 ( *this ).XFORM( 3, 0 ) = t.i();
 ( *this ).XFORM( 3, 1 ) = t.j();
 ( *this ).XFORM( 3, 2 ) = t.k();
 
 // M should transform p to q
 assert( ( q1 == p1 * *this ) && ( q2 == p2 * *this ) && ( q3 == p3 * *this ) );
}
 
} // namespace moab