1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkMatrix3x3.cxx
5
6 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 All rights reserved.
8 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9
10 This software is distributed WITHOUT ANY WARRANTY; without even
11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12 PURPOSE. See the above copyright notice for more information.
13
14 =========================================================================*/
15 #include "vtkMatrix3x3.h"
16 #include "vtkMath.h"
17 #include "vtkObjectFactory.h"
18
19 #include <cmath>
20 #include <cstdlib>
21
22 vtkStandardNewMacro(vtkMatrix3x3);
23
24 //------------------------------------------------------------------------------
vtkMatrix3x3()25 vtkMatrix3x3::vtkMatrix3x3()
26 {
27 vtkMatrix3x3::Identity(*this->Element);
28 }
29
30 //------------------------------------------------------------------------------
31 vtkMatrix3x3::~vtkMatrix3x3() = default;
32
33 //------------------------------------------------------------------------------
Zero(double elements[9])34 void vtkMatrix3x3::Zero(double elements[9])
35 {
36 for (int i = 0; i < 9; ++i)
37 {
38 elements[i] = 0.0;
39 }
40 }
41
42 //------------------------------------------------------------------------------
Identity(double elements[9])43 void vtkMatrix3x3::Identity(double elements[9])
44 {
45 elements[0] = elements[4] = elements[8] = 1.0;
46 elements[1] = elements[2] = elements[3] = elements[5] = elements[6] = elements[7] = 0.0;
47 }
48
49 //------------------------------------------------------------------------------
50 namespace
51 { // Enclose private helper function in anonymous namespace
52
53 template <class T1, class T2, class T3>
vtkMatrix3x3MultiplyPoint(T1 elem[9],T2 in[3],T3 out[3])54 void vtkMatrix3x3MultiplyPoint(T1 elem[9], T2 in[3], T3 out[3])
55 {
56 T3 v1 = in[0];
57 T3 v2 = in[1];
58 T3 v3 = in[2];
59
60 out[0] = v1 * elem[0] + v2 * elem[1] + v3 * elem[2];
61 out[1] = v1 * elem[3] + v2 * elem[4] + v3 * elem[5];
62 out[2] = v1 * elem[6] + v2 * elem[7] + v3 * elem[8];
63 }
64
65 } // End anonymous namespace
66
67 //------------------------------------------------------------------------------
68 // Multiply this matrix by a point (in homogeneous coordinates).
69 // and return the result in result. The in[3] and result[3]
70 // arrays must both be allocated but they can be the same array.
MultiplyPoint(const double elements[9],const float in[3],float result[3])71 void vtkMatrix3x3::MultiplyPoint(const double elements[9], const float in[3], float result[3])
72 {
73 vtkMatrix3x3MultiplyPoint(elements, in, result);
74 }
75
76 //------------------------------------------------------------------------------
MultiplyPoint(const double elements[9],const double in[3],double result[3])77 void vtkMatrix3x3::MultiplyPoint(const double elements[9], const double in[3], double result[3])
78 {
79 vtkMatrix3x3MultiplyPoint(elements, in, result);
80 }
81
82 //------------------------------------------------------------------------------
83 // Multiplies matrices a and b and stores the result in c.
Multiply3x3(const double a[9],const double b[9],double c[9])84 void vtkMatrix3x3::Multiply3x3(const double a[9], const double b[9], double c[9])
85 {
86 double accum[9];
87
88 for (int i = 0; i < 9; i += 3)
89 {
90 for (int k = 0; k < 3; k++)
91 {
92 accum[i + k] = a[i + 0] * b[k + 0] + a[i + 1] * b[k + 3] + a[i + 2] * b[k + 6];
93 }
94 }
95
96 // Copy to final dest
97 for (int j = 0; j < 9; j++)
98 {
99 c[j] = accum[j];
100 }
101 }
102
103 //------------------------------------------------------------------------------
104 // Matrix Inversion (adapted from Richard Carling in "Graphics Gems,"
105 // Academic Press, 1990).
Invert(const double inElements[9],double outElements[9])106 void vtkMatrix3x3::Invert(const double inElements[9], double outElements[9])
107 {
108 // inverse( original_matrix, inverse_matrix )
109 // calculate the inverse of a 3x3 matrix
110 //
111 // -1
112 // A = ___1__ adjoint A
113 // det A
114 //
115
116 // calculate the 3x3 determinent
117 // if the determinent is zero,
118 // then the inverse matrix is not unique.
119
120 double det = vtkMatrix3x3::Determinant(inElements);
121 if (det == 0.0)
122 {
123 return;
124 }
125
126 // calculate the adjoint matrix
127 vtkMatrix3x3::Adjoint(inElements, outElements);
128
129 // scale the adjoint matrix to get the inverse
130 for (int i = 0; i < 9; i++)
131 {
132 outElements[i] /= det;
133 }
134 }
135
136 //------------------------------------------------------------------------------
Determinant(const double elements[9])137 double vtkMatrix3x3::Determinant(const double elements[9])
138 {
139 return vtkMath::Determinant3x3(elements[0], elements[1], elements[2], elements[3], elements[4],
140 elements[5], elements[6], elements[7], elements[8]);
141 }
142
143 //------------------------------------------------------------------------------
Adjoint(const double inElements[9],double outElements[9])144 void vtkMatrix3x3::Adjoint(const double inElements[9], double outElements[9])
145 {
146 //
147 // adjoint( original_matrix, inverse_matrix )
148 //
149 // calculate the adjoint of a 3x3 matrix
150 //
151 // Let a denote the minor determinant of matrix A obtained by
152 // ij
153 //
154 // deleting the ith row and jth column from A.
155 //
156 // i+j
157 // Let b = (-1) a
158 // ij ji
159 //
160 // The matrix B = (b ) is the adjoint of A
161 // ij
162 //
163 double a1, a2, a3, b1, b2, b3, c1, c2, c3;
164
165 // assign to individual variable names to aid
166 // selecting correct values
167 a1 = inElements[0];
168 b1 = inElements[1];
169 c1 = inElements[2];
170 a2 = inElements[3];
171 b2 = inElements[4];
172 c2 = inElements[5];
173 a3 = inElements[6];
174 b3 = inElements[7];
175 c3 = inElements[8];
176
177 // row column labeling reversed since we transpose rows & columns
178
179 outElements[0] = vtkMath::Determinant2x2(b2, b3, c2, c3);
180 outElements[3] = -vtkMath::Determinant2x2(a2, a3, c2, c3);
181 outElements[6] = vtkMath::Determinant2x2(a2, a3, b2, b3);
182
183 outElements[1] = -vtkMath::Determinant2x2(b1, b3, c1, c3);
184 outElements[4] = vtkMath::Determinant2x2(a1, a3, c1, c3);
185 outElements[7] = -vtkMath::Determinant2x2(a1, a3, b1, b3);
186
187 outElements[2] = vtkMath::Determinant2x2(b1, b2, c1, c2);
188 outElements[5] = -vtkMath::Determinant2x2(a1, a2, c1, c2);
189 outElements[8] = vtkMath::Determinant2x2(a1, a2, b1, b2);
190 }
191
192 //------------------------------------------------------------------------------
DeepCopy(double elements[9],const double newElements[9])193 void vtkMatrix3x3::DeepCopy(double elements[9], const double newElements[9])
194 {
195 for (int i = 0; i < 9; ++i)
196 {
197 elements[i] = newElements[i];
198 }
199 }
200
201 //------------------------------------------------------------------------------
202 // Transpose the matrix and put it into out.
Transpose(const double inElements[9],double outElements[9])203 void vtkMatrix3x3::Transpose(const double inElements[9], double outElements[9])
204 {
205 for (int i = 0; i < 3; ++i)
206 {
207 for (int j = i; j < 3; ++j)
208 {
209 double temp = inElements[3 * i + j];
210 outElements[3 * i + j] = inElements[3 * j + i];
211 outElements[3 * j + i] = temp;
212 }
213 }
214 }
215
216 //------------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)217 void vtkMatrix3x3::PrintSelf(ostream& os, vtkIndent indent)
218 {
219 this->Superclass::PrintSelf(os, indent);
220
221 os << indent << "Elements:\n";
222 for (int i = 0; i < 3; ++i)
223 {
224 os << indent;
225 for (int j = 0; j < 3; ++j)
226 {
227 os << "\t" << this->Element[i][j];
228 }
229 os << "\n";
230 }
231 }
232