Branch data Line data Source code
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 : : #include "moab/HomXform.hpp"
17 : : #include <assert.h>
18 : :
19 : : namespace moab
20 : : {
21 : :
22 : 58 : HomCoord HomCoord::unitv[3] = { HomCoord( 1, 0, 0 ), HomCoord( 0, 1, 0 ), HomCoord( 0, 0, 1 ) };
23 : 58 : HomCoord HomCoord::IDENTITY( 1, 1, 1 );
24 : :
25 : 132 : HomCoord& HomCoord::getUnitv( int c )
26 : : {
27 : 132 : return unitv[c];
28 : : }
29 : :
30 : : int dum[] = { 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1 };
31 : 58 : HomXform HomXform::IDENTITY( dum );
32 : :
33 : 41 : void HomXform::three_pt_xform( const HomCoord& p1, const HomCoord& q1, const HomCoord& p2, const HomCoord& q2,
34 : : const HomCoord& p3, const HomCoord& q3 )
35 : : {
36 : : // pmin and pmax are min and max bounding box corners which are mapped to
37 : : // qmin and qmax, resp. qmin and qmax are not necessarily min/max corners,
38 : : // since the mapping can change the orientation of the box in the q reference
39 : : // frame. Re-interpreting the min/max bounding box corners does not change
40 : : // the mapping.
41 : :
42 : : // change that: base on three points for now (figure out whether we can
43 : : // just use two later); three points are assumed to define an orthogonal
44 : : // system such that (p2-p1)%(p3-p1) = 0
45 : :
46 : : // use the three point rule to compute the mapping, from Mortensen,
47 : : // "Geometric Modeling". If p1, p2, p3 and q1, q2, q3 are three points in
48 : : // the two coordinate systems, the three pt rule is:
49 : : //
50 : : // v1 = p2 - p1
51 : : // v2 = p3 - p1
52 : : // v3 = v1 x v2
53 : : // w1-w3 similar, with q1-q3
54 : : // V = matrix with v1-v3 as rows
55 : : // W similar, with w1-w3
56 : : // R = V^-1 * W
57 : : // t = q1 - p1 * W
58 : : // Form transform matrix M from R, t
59 : :
60 : : // check to see whether unity transform applies
61 [ + - ][ + + ]: 41 : if( p1 == q1 && p2 == q2 && p3 == q3 )
[ + - ][ + + ]
[ + - ][ + - ]
[ + + ]
62 : : {
63 [ + - ]: 22 : *this = HomXform::IDENTITY;
64 : 41 : return;
65 : : }
66 : :
67 : : // first, construct 3 pts from input
68 [ + - ]: 19 : HomCoord v1 = p2 - p1;
69 [ + - ][ + + ]: 19 : assert( v1.i() != 0 || v1.j() != 0 || v1.k() != 0 );
[ + - ][ - + ]
[ # # ][ # # ]
70 [ + - ]: 19 : HomCoord v2 = p3 - p1;
71 [ + - ]: 19 : HomCoord v3 = v1 * v2;
72 : :
73 [ + - ][ + + ]: 19 : if( v3.length_squared() == 0 )
74 : : {
75 : : // 1d coordinate system; set one of v2's coordinates such that
76 : : // it's orthogonal to v1
77 [ + - ][ + + ]: 7 : if( v1.i() == 0 )
78 [ + - ]: 1 : v2.set( 1, 0, 0 );
79 [ + - ][ + - ]: 6 : else if( v1.j() == 0 )
80 [ + - ]: 6 : v2.set( 0, 1, 0 );
81 [ # # ][ # # ]: 0 : else if( v1.k() == 0 )
82 [ # # ]: 0 : v2.set( 0, 0, 1 );
83 : : else
84 : 0 : assert( false );
85 [ + - ][ + - ]: 7 : v3 = v1 * v2;
86 [ + - ][ - + ]: 7 : assert( v3.length_squared() != 0 );
87 : : }
88 : : // assert to make sure they're each orthogonal
89 [ + - ][ + - ]: 19 : assert( v1 % v2 == 0 && v1 % v3 == 0 && v2 % v3 == 0 );
[ + - ][ + - ]
[ + - ][ - + ]
90 [ + - ]: 19 : v1.normalize();
91 [ + - ]: 19 : v2.normalize();
92 [ + - ]: 19 : v3.normalize();
93 : : // Make sure h is set to zero here, since it'll mess up things if it's one
94 : 19 : v1.homCoord[3] = v2.homCoord[3] = v3.homCoord[3] = 0;
95 : :
96 [ + - ]: 19 : HomCoord w1 = q2 - q1;
97 [ + - ][ + + ]: 19 : assert( w1.i() != 0 || w1.j() != 0 || w1.k() != 0 );
[ + - ][ - + ]
[ # # ][ # # ]
98 [ + - ]: 19 : HomCoord w2 = q3 - q1;
99 [ + - ]: 19 : HomCoord w3 = w1 * w2;
100 : :
101 [ + - ][ + + ]: 19 : if( w3.length_squared() == 0 )
102 : : {
103 : : // 1d coordinate system; set one of w2's coordinates such that
104 : : // it's orthogonal to w1
105 [ + - ][ + + ]: 7 : if( w1.i() == 0 )
106 [ + - ]: 3 : w2.set( 1, 0, 0 );
107 [ + - ][ + - ]: 4 : else if( w1.j() == 0 )
108 [ + - ]: 4 : w2.set( 0, 1, 0 );
109 [ # # ][ # # ]: 0 : else if( w1.k() == 0 )
110 [ # # ]: 0 : w2.set( 0, 0, 1 );
111 : : else
112 : 0 : assert( false );
113 [ + - ][ + - ]: 7 : w3 = w1 * w2;
114 [ + - ][ - + ]: 7 : assert( w3.length_squared() != 0 );
115 : : }
116 : : // assert to make sure they're each orthogonal
117 [ + - ][ + - ]: 19 : assert( w1 % w2 == 0 && w1 % w3 == 0 && w2 % w3 == 0 );
[ + - ][ + - ]
[ + - ][ - + ]
118 [ + - ]: 19 : w1.normalize();
119 [ + - ]: 19 : w2.normalize();
120 [ + - ]: 19 : w3.normalize();
121 : : // Make sure h is set to zero here, since it'll mess up things if it's one
122 : 19 : w1.homCoord[3] = w2.homCoord[3] = w3.homCoord[3] = 0;
123 : :
124 : : // form v^-1 as transpose (ok for orthogonal vectors); put directly into
125 : : // transform matrix, since it's eventually going to become R
126 [ + - ][ + - ]: 19 : *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 );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
127 : :
128 : : // multiply by w to get R
129 [ + - ][ + - ]: 19 : *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 );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
130 : :
131 : : // compute t and put into last row
132 [ + - ][ + - ]: 19 : HomCoord t = q1 - p1 * *this;
133 [ + - ]: 19 : ( *this ).XFORM( 3, 0 ) = t.i();
134 [ + - ]: 19 : ( *this ).XFORM( 3, 1 ) = t.j();
135 [ + - ]: 19 : ( *this ).XFORM( 3, 2 ) = t.k();
136 : :
137 : : // M should transform p to q
138 [ + - ][ + - ]: 19 : assert( ( q1 == p1 * *this ) && ( q2 == p2 * *this ) && ( q3 == p3 * *this ) );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ - + ][ + - ]
[ + - ][ + - ]
[ # # # #
# # ]
139 : : }
140 : :
141 [ + - ][ + - ]: 232 : } // namespace moab
|
