Mercurial > hg > orthanc-stone

comparison Framework/Toolbox/LinearAlgebra.cpp @ 158:a053ca7fa5c6 wasm

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
LinearAlgebra toolbox
author Sebastien Jodogne <s.jodogne@gmail.com>
date Wed, 14 Feb 2018 08:58:31 +0100
parents
children 0a73d76333db
comparison
equal deleted inserted replaced
157:2309e8d86efe 158:a053ca7fa5c6
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-2018 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 #include "LinearAlgebra.h"
23
24 #include <Core/Logging.h>
25 #include <Core/OrthancException.h>
26 #include <Core/Toolbox.h>
27
28 #include <stdio.h>
29 #include <boost/lexical_cast.hpp>
30
31 namespace OrthancStone
32 {
33 namespace LinearAlgebra
34 {
35 void Print(const Vector& v)
36 {
37 for (size_t i = 0; i < v.size(); i++)
38 {
39 printf("%g\n", v[i]);
40 //printf("%8.2f\n", v[i]);
41 }
42 printf("\n");
43 }
44
45
46 void Print(const Matrix& m)
47 {
48 for (size_t i = 0; i < m.size1(); i++)
49 {
50 for (size_t j = 0; j < m.size2(); j++)
51 {
52 printf("%g ", m(i,j));
53 //printf("%8.2f ", m(i,j));
54 }
55 printf("\n");
56 }
57 printf("\n");
58 }
59
60
61 bool ParseVector(Vector& target,
62 const std::string& value)
63 {
64 std::vector<std::string> items;
65 Orthanc::Toolbox::TokenizeString(items, value, '\\');
66
67 target.resize(items.size());
68
69 for (size_t i = 0; i < items.size(); i++)
70 {
71 try
72 {
73 target[i] = boost::lexical_cast<double>(Orthanc::Toolbox::StripSpaces(items[i]));
74 }
75 catch (boost::bad_lexical_cast&)
76 {
77 target.clear();
78 return false;
79 }
80 }
81
82 return true;
83 }
84
85
86 bool ParseVector(Vector& target,
87 const Orthanc::DicomMap& dataset,
88 const Orthanc::DicomTag& tag)
89 {
90 std::string value;
91 return (dataset.CopyToString(value, tag, false) &&
92 ParseVector(target, value));
93 }
94
95
96 void AssignVector(Vector& v,
97 double v1,
98 double v2)
99 {
100 v.resize(2);
101 v[0] = v1;
102 v[1] = v2;
103 }
104
105
106 void AssignVector(Vector& v,
107 double v1,
108 double v2,
109 double v3)
110 {
111 v.resize(3);
112 v[0] = v1;
113 v[1] = v2;
114 v[2] = v3;
115 }
116
117
118 bool IsNear(double x,
119 double y)
120 {
121 // As most input is read as single-precision numbers, we take the
122 // epsilon machine for float32 into consideration to compare numbers
123 return IsNear(x, y, 10.0 * std::numeric_limits<float>::epsilon());
124 }
125
126
127 void NormalizeVector(Vector& u)
128 {
129 double norm = boost::numeric::ublas::norm_2(u);
130 if (!IsCloseToZero(norm))
131 {
132 u = u / norm;
133 }
134 }
135
136
137 void CrossProduct(Vector& result,
138 const Vector& u,
139 const Vector& v)
140 {
141 if (u.size() != 3 ||
142 v.size() != 3)
143 {
144 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
145 }
146
147 result.resize(3);
148
149 result[0] = u[1] * v[2] - u[2] * v[1];
150 result[1] = u[2] * v[0] - u[0] * v[2];
151 result[2] = u[0] * v[1] - u[1] * v[0];
152 }
153
154
155 void FillMatrix(Matrix& target,
156 size_t rows,
157 size_t columns,
158 const double values[])
159 {
160 target.resize(rows, columns);
161
162 size_t index = 0;
163
164 for (size_t y = 0; y < rows; y++)
165 {
166 for (size_t x = 0; x < columns; x++, index++)
167 {
168 target(y, x) = values[index];
169 }
170 }
171 }
172
173
174 void FillVector(Vector& target,
175 size_t size,
176 const double values[])
177 {
178 target.resize(size);
179
180 for (size_t i = 0; i < size; i++)
181 {
182 target[i] = values[i];
183 }
184 }
185
186
187 void Convert(Matrix& target,
188 const Vector& source)
189 {
190 const size_t n = source.size();
191
192 target.resize(n, 1);
193
194 for (size_t i = 0; i < n; i++)
195 {
196 target(i, 0) = source[i];
197 }
198 }
199 }
200 }