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