Mercurial > hg > orthanc-stone
comparison OrthancStone/Sources/Toolbox/LinearAlgebra.h @ 1512:244ad1e4e76a
reorganization of folders
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Tue, 07 Jul 2020 16:21:02 +0200 |
parents | Framework/Toolbox/LinearAlgebra.h@30deba7bc8e2 |
children | 4fb8fdf03314 |
comparison
equal
deleted
inserted
replaced
1511:9dfeee74c1e6 | 1512:244ad1e4e76a |
---|---|
1 /** | |
2 * Stone of Orthanc | |
3 * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics | |
4 * Department, University Hospital of Liege, Belgium | |
5 * Copyright (C) 2017-2020 Osimis S.A., Belgium | |
6 * | |
7 * This program is free software: you can redistribute it and/or | |
8 * modify it under the terms of the GNU Affero General Public License | |
9 * as published by the Free Software Foundation, either version 3 of | |
10 * the License, or (at your option) any later version. | |
11 * | |
12 * This program is distributed in the hope that it will be useful, but | |
13 * WITHOUT ANY WARRANTY; without even the implied warranty of | |
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
15 * Affero General Public License for more details. | |
16 * | |
17 * You should have received a copy of the GNU Affero General Public License | |
18 * along with this program. If not, see <http://www.gnu.org/licenses/>. | |
19 **/ | |
20 | |
21 | |
22 #pragma once | |
23 | |
24 // Patch for ublas in Boost 1.64.0 | |
25 // https://github.com/dealii/dealii/issues/4302 | |
26 #include <boost/version.hpp> | |
27 #if BOOST_VERSION >= 106300 // or 64, need to check | |
28 # include <boost/serialization/array_wrapper.hpp> | |
29 #endif | |
30 | |
31 #include <DicomFormat/DicomMap.h> | |
32 | |
33 #include <boost/numeric/ublas/matrix.hpp> | |
34 #include <boost/numeric/ublas/vector.hpp> | |
35 | |
36 namespace OrthancStone | |
37 { | |
38 typedef boost::numeric::ublas::matrix<double> Matrix; | |
39 typedef boost::numeric::ublas::vector<double> Vector; | |
40 | |
41 // logs, debugging... | |
42 std::ostream& operator<<(std::ostream& s, const Vector& vec); | |
43 std::ostream& operator<<(std::ostream& s, const Matrix& m); | |
44 | |
45 namespace LinearAlgebra | |
46 { | |
47 void Print(const Vector& v); | |
48 | |
49 void Print(const Matrix& m); | |
50 | |
51 bool ParseVector(Vector& target, | |
52 const std::string& s); | |
53 | |
54 bool ParseVector(Vector& target, | |
55 const Orthanc::DicomMap& dataset, | |
56 const Orthanc::DicomTag& tag); | |
57 | |
58 inline void AssignVector(Vector& v, | |
59 double v1, | |
60 double v2) | |
61 { | |
62 v.resize(2); | |
63 v[0] = v1; | |
64 v[1] = v2; | |
65 } | |
66 | |
67 | |
68 inline void AssignVector(Vector& v, | |
69 double v1) | |
70 { | |
71 v.resize(1); | |
72 v[0] = v1; | |
73 } | |
74 | |
75 | |
76 inline void AssignVector(Vector& v, | |
77 double v1, | |
78 double v2, | |
79 double v3) | |
80 { | |
81 v.resize(3); | |
82 v[0] = v1; | |
83 v[1] = v2; | |
84 v[2] = v3; | |
85 } | |
86 | |
87 | |
88 inline void AssignVector(Vector& v, | |
89 double v1, | |
90 double v2, | |
91 double v3, | |
92 double v4) | |
93 { | |
94 v.resize(4); | |
95 v[0] = v1; | |
96 v[1] = v2; | |
97 v[2] = v3; | |
98 v[3] = v4; | |
99 } | |
100 | |
101 | |
102 inline Vector CreateVector(double v1) | |
103 { | |
104 Vector v; | |
105 AssignVector(v, v1); | |
106 return v; | |
107 } | |
108 | |
109 | |
110 inline Vector CreateVector(double v1, | |
111 double v2) | |
112 { | |
113 Vector v; | |
114 AssignVector(v, v1, v2); | |
115 return v; | |
116 } | |
117 | |
118 | |
119 inline Vector CreateVector(double v1, | |
120 double v2, | |
121 double v3) | |
122 { | |
123 Vector v; | |
124 AssignVector(v, v1, v2, v3); | |
125 return v; | |
126 } | |
127 | |
128 | |
129 inline Vector CreateVector(double v1, | |
130 double v2, | |
131 double v3, | |
132 double v4) | |
133 { | |
134 Vector v; | |
135 AssignVector(v, v1, v2, v3, v4); | |
136 return v; | |
137 } | |
138 | |
139 | |
140 inline bool IsNear(double x, | |
141 double y, | |
142 double threshold) | |
143 { | |
144 return fabs(x - y) <= threshold; | |
145 } | |
146 | |
147 inline bool IsNear(double x, | |
148 double y) | |
149 { | |
150 // As most input is read as single-precision numbers, we take the | |
151 // epsilon machine for float32 into consideration to compare numbers | |
152 return IsNear(x, y, 10.0 * std::numeric_limits<float>::epsilon()); | |
153 } | |
154 | |
155 inline bool IsCloseToZero(double x) | |
156 { | |
157 return IsNear(x, 0.0); | |
158 } | |
159 | |
160 void NormalizeVector(Vector& u); | |
161 | |
162 void CrossProduct(Vector& result, | |
163 const Vector& u, | |
164 const Vector& v); | |
165 | |
166 double DotProduct(const Vector& u, const Vector& v); | |
167 | |
168 void FillMatrix(Matrix& target, | |
169 size_t rows, | |
170 size_t columns, | |
171 const double values[]); | |
172 | |
173 void FillVector(Vector& target, | |
174 size_t size, | |
175 const double values[]); | |
176 | |
177 void Convert(Matrix& target, | |
178 const Vector& source); | |
179 | |
180 inline Matrix Transpose(const Matrix& a) | |
181 { | |
182 return boost::numeric::ublas::trans(a); | |
183 } | |
184 | |
185 | |
186 inline Matrix IdentityMatrix(size_t size) | |
187 { | |
188 return boost::numeric::ublas::identity_matrix<double>(size); | |
189 } | |
190 | |
191 | |
192 inline Matrix ZeroMatrix(size_t size1, | |
193 size_t size2) | |
194 { | |
195 return boost::numeric::ublas::zero_matrix<double>(size1, size2); | |
196 } | |
197 | |
198 | |
199 inline Matrix Product(const Matrix& a, | |
200 const Matrix& b) | |
201 { | |
202 return boost::numeric::ublas::prod(a, b); | |
203 } | |
204 | |
205 | |
206 inline Vector Product(const Matrix& a, | |
207 const Vector& b) | |
208 { | |
209 return boost::numeric::ublas::prod(a, b); | |
210 } | |
211 | |
212 | |
213 inline Matrix Product(const Matrix& a, | |
214 const Matrix& b, | |
215 const Matrix& c) | |
216 { | |
217 return Product(a, Product(b, c)); | |
218 } | |
219 | |
220 | |
221 inline Matrix Product(const Matrix& a, | |
222 const Matrix& b, | |
223 const Matrix& c, | |
224 const Matrix& d) | |
225 { | |
226 return Product(a, Product(b, c, d)); | |
227 } | |
228 | |
229 | |
230 inline Matrix Product(const Matrix& a, | |
231 const Matrix& b, | |
232 const Matrix& c, | |
233 const Matrix& d, | |
234 const Matrix& e) | |
235 { | |
236 return Product(a, Product(b, c, d, e)); | |
237 } | |
238 | |
239 | |
240 inline Vector Product(const Matrix& a, | |
241 const Matrix& b, | |
242 const Vector& c) | |
243 { | |
244 return Product(Product(a, b), c); | |
245 } | |
246 | |
247 | |
248 inline Vector Product(const Matrix& a, | |
249 const Matrix& b, | |
250 const Matrix& c, | |
251 const Vector& d) | |
252 { | |
253 return Product(Product(a, b, c), d); | |
254 } | |
255 | |
256 | |
257 double ComputeDeterminant(const Matrix& a); | |
258 | |
259 bool IsOrthogonalMatrix(const Matrix& q, | |
260 double threshold); | |
261 | |
262 bool IsOrthogonalMatrix(const Matrix& q); | |
263 | |
264 bool IsRotationMatrix(const Matrix& r, | |
265 double threshold); | |
266 | |
267 bool IsRotationMatrix(const Matrix& r); | |
268 | |
269 void InvertUpperTriangularMatrix(Matrix& output, | |
270 const Matrix& k); | |
271 | |
272 /** | |
273 * This function computes the RQ decomposition of a 3x3 matrix, | |
274 * using Givens rotations. Reference: Algorithm A4.1 (page 579) of | |
275 * "Multiple View Geometry in Computer Vision" (2nd edition). The | |
276 * output matrix "Q" is a rotation matrix, and "R" is upper | |
277 * triangular. | |
278 **/ | |
279 void RQDecomposition3x3(Matrix& r, | |
280 Matrix& q, | |
281 const Matrix& a); | |
282 | |
283 void InvertMatrix(Matrix& target, | |
284 const Matrix& source); | |
285 | |
286 // This is the same as "InvertMatrix()", but without exception | |
287 bool InvertMatrixUnsafe(Matrix& target, | |
288 const Matrix& source); | |
289 | |
290 void CreateSkewSymmetric(Matrix& s, | |
291 const Vector& v); | |
292 | |
293 Matrix InvertScalingTranslationMatrix(const Matrix& t); | |
294 | |
295 bool IsShearMatrix(const Matrix& shear); | |
296 | |
297 Matrix InvertShearMatrix(const Matrix& shear); | |
298 }; | |
299 } |