Mercurial > hg > orthanc-stone

comparison OrthancStone/Sources/Toolbox/LinearAlgebra.h @ 1512:244ad1e4e76a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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 }