Mercurial > hg > orthanc-stone
annotate Framework/Toolbox/GeometryToolbox.cpp @ 40:7207a407bcd8
shared copyright with osimis
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Wed, 04 Jan 2017 16:37:42 +0100 |
parents | 517c46f527cd |
children | 28956ed68280 |
rev | line source |
---|---|
0 | 1 /** |
2 * Stone of Orthanc | |
3 * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics | |
4 * Department, University Hospital of Liege, Belgium | |
40
7207a407bcd8
shared copyright with osimis
Sebastien Jodogne <s.jodogne@gmail.com>
parents:
32
diff
changeset
|
5 * Copyright (C) 2017 Osimis, Belgium |
0 | 6 * |
7 * This program is free software: you can redistribute it and/or | |
8 * modify it under the terms of the GNU General Public License as | |
9 * published by the Free Software Foundation, either version 3 of the | |
10 * License, or (at your option) any later version. | |
11 * | |
12 * In addition, as a special exception, the copyright holders of this | |
13 * program give permission to link the code of its release with the | |
14 * OpenSSL project's "OpenSSL" library (or with modified versions of it | |
15 * that use the same license as the "OpenSSL" library), and distribute | |
16 * the linked executables. You must obey the GNU General Public License | |
17 * in all respects for all of the code used other than "OpenSSL". If you | |
18 * modify file(s) with this exception, you may extend this exception to | |
19 * your version of the file(s), but you are not obligated to do so. If | |
20 * you do not wish to do so, delete this exception statement from your | |
21 * version. If you delete this exception statement from all source files | |
22 * in the program, then also delete it here. | |
23 * | |
24 * This program is distributed in the hope that it will be useful, but | |
25 * WITHOUT ANY WARRANTY; without even the implied warranty of | |
26 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
27 * General Public License for more details. | |
28 * | |
29 * You should have received a copy of the GNU General Public License | |
30 * along with this program. If not, see <http://www.gnu.org/licenses/>. | |
31 **/ | |
32 | |
33 | |
34 #include "GeometryToolbox.h" | |
35 | |
32 | 36 #include "../../Resources/Orthanc/Core/Logging.h" |
16 | 37 #include "../../Resources/Orthanc/Core/OrthancException.h" |
38 #include "../../Resources/Orthanc/Core/Toolbox.h" | |
0 | 39 |
40 #include <stdio.h> | |
41 #include <boost/lexical_cast.hpp> | |
42 | |
43 namespace OrthancStone | |
44 { | |
45 namespace GeometryToolbox | |
46 { | |
47 void Print(const Vector& v) | |
48 { | |
49 for (size_t i = 0; i < v.size(); i++) | |
50 { | |
51 printf("%8.2f\n", v[i]); | |
52 } | |
53 printf("\n"); | |
54 } | |
55 | |
56 | |
57 bool ParseVector(Vector& target, | |
58 const std::string& value) | |
59 { | |
60 std::vector<std::string> items; | |
61 Orthanc::Toolbox::TokenizeString(items, value, '\\'); | |
62 | |
63 target.resize(items.size()); | |
64 | |
65 for (size_t i = 0; i < items.size(); i++) | |
66 { | |
67 try | |
68 { | |
69 target[i] = boost::lexical_cast<double>(Orthanc::Toolbox::StripSpaces(items[i])); | |
70 } | |
71 catch (boost::bad_lexical_cast&) | |
72 { | |
73 target.clear(); | |
74 return false; | |
75 } | |
76 } | |
77 | |
78 return true; | |
79 } | |
80 | |
81 | |
32 | 82 bool ParseVector(Vector& target, |
83 const OrthancPlugins::IDicomDataset& dataset, | |
84 const OrthancPlugins::DicomPath& tag) | |
85 { | |
86 std::string value; | |
87 return (dataset.GetStringValue(value, tag) && | |
88 ParseVector(target, value)); | |
89 } | |
90 | |
91 | |
0 | 92 void AssignVector(Vector& v, |
93 double v1, | |
94 double v2) | |
95 { | |
96 v.resize(2); | |
97 v[0] = v1; | |
98 v[1] = v2; | |
99 } | |
100 | |
101 | |
102 void AssignVector(Vector& v, | |
103 double v1, | |
104 double v2, | |
105 double v3) | |
106 { | |
107 v.resize(3); | |
108 v[0] = v1; | |
109 v[1] = v2; | |
110 v[2] = v3; | |
111 } | |
112 | |
113 | |
114 bool IsNear(double x, | |
115 double y) | |
116 { | |
117 // As most input is read as single-precision numbers, we take the | |
118 // epsilon machine for float32 into consideration to compare numbers | |
119 return IsNear(x, y, 10.0 * std::numeric_limits<float>::epsilon()); | |
120 } | |
121 | |
122 | |
123 void NormalizeVector(Vector& u) | |
124 { | |
125 double norm = boost::numeric::ublas::norm_2(u); | |
126 if (!IsCloseToZero(norm)) | |
127 { | |
128 u = u / norm; | |
129 } | |
130 } | |
131 | |
132 | |
133 void CrossProduct(Vector& result, | |
134 const Vector& u, | |
135 const Vector& v) | |
136 { | |
137 if (u.size() != 3 || | |
138 v.size() != 3) | |
139 { | |
140 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); | |
141 } | |
142 | |
143 result.resize(3); | |
144 | |
145 result[0] = u[1] * v[2] - u[2] * v[1]; | |
146 result[1] = u[2] * v[0] - u[0] * v[2]; | |
147 result[2] = u[0] * v[1] - u[1] * v[0]; | |
148 } | |
149 | |
150 | |
151 void ProjectPointOntoPlane(Vector& result, | |
152 const Vector& point, | |
153 const Vector& planeNormal, | |
154 const Vector& planeOrigin) | |
155 { | |
156 double norm = boost::numeric::ublas::norm_2(planeNormal); | |
157 if (IsCloseToZero(norm)) | |
158 { | |
159 // Division by zero | |
160 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); | |
161 } | |
162 | |
163 // Make sure the norm of the normal is 1 | |
164 Vector n; | |
165 n = planeNormal / norm; | |
166 | |
167 // Algebraic form of line–plane intersection, where the line passes | |
168 // through "point" along the direction "normal" (thus, l == n) | |
169 // https://en.wikipedia.org/wiki/Line%E2%80%93plane_intersection#Algebraic_form | |
170 result = boost::numeric::ublas::inner_prod(planeOrigin - point, n) * n + point; | |
171 } | |
172 | |
173 | |
174 | |
175 bool IsParallelOrOpposite(bool& isOpposite, | |
176 const Vector& u, | |
177 const Vector& v) | |
178 { | |
179 // The dot product of the two vectors gives the cosine of the angle | |
180 // between the vectors | |
181 // https://en.wikipedia.org/wiki/Dot_product | |
182 | |
183 double normU = boost::numeric::ublas::norm_2(u); | |
184 double normV = boost::numeric::ublas::norm_2(v); | |
185 | |
186 if (IsCloseToZero(normU) || | |
187 IsCloseToZero(normV)) | |
188 { | |
189 return false; | |
190 } | |
191 | |
192 double cosAngle = boost::numeric::ublas::inner_prod(u, v) / (normU * normV); | |
193 | |
194 // The angle must be zero, so the cosine must be almost equal to | |
195 // cos(0) == 1 (or cos(180) == -1 if allowOppositeDirection == true) | |
196 | |
197 if (IsCloseToZero(cosAngle - 1.0)) | |
198 { | |
199 isOpposite = false; | |
200 return true; | |
201 } | |
202 else if (IsCloseToZero(fabs(cosAngle) - 1.0)) | |
203 { | |
204 isOpposite = true; | |
205 return true; | |
206 } | |
207 else | |
208 { | |
209 return false; | |
210 } | |
211 } | |
212 | |
213 | |
214 bool IsParallel(const Vector& u, | |
215 const Vector& v) | |
216 { | |
217 bool isOpposite; | |
218 return (IsParallelOrOpposite(isOpposite, u, v) && | |
219 !isOpposite); | |
220 } | |
221 | |
222 | |
223 bool IntersectTwoPlanes(Vector& p, | |
224 Vector& direction, | |
225 const Vector& origin1, | |
226 const Vector& normal1, | |
227 const Vector& origin2, | |
228 const Vector& normal2) | |
229 { | |
230 // This is "Intersection of 2 Planes", possibility "(C) 3 Plane | |
231 // Intersect Point" of: | |
232 // http://geomalgorithms.com/a05-_intersect-1.html | |
233 | |
234 // The direction of the line of intersection is orthogonal to the | |
235 // normal of both planes | |
236 CrossProduct(direction, normal1, normal2); | |
237 | |
238 double norm = boost::numeric::ublas::norm_2(direction); | |
239 if (IsCloseToZero(norm)) | |
240 { | |
241 // The two planes are parallel or coincident | |
242 return false; | |
243 } | |
244 | |
245 double d1 = -boost::numeric::ublas::inner_prod(normal1, origin1); | |
246 double d2 = -boost::numeric::ublas::inner_prod(normal2, origin2); | |
247 Vector tmp = d2 * normal1 - d1 * normal2; | |
248 | |
249 CrossProduct(p, tmp, direction); | |
250 p /= norm; | |
251 | |
252 return true; | |
253 } | |
254 | |
255 | |
256 bool ClipLineToRectangle(double& x1, // Coordinates of the clipped line (out) | |
257 double& y1, | |
258 double& x2, | |
259 double& y2, | |
260 const double ax, // Two points defining the line (in) | |
261 const double ay, | |
262 const double bx, | |
263 const double by, | |
264 const double& xmin, // Coordinates of the rectangle (in) | |
265 const double& ymin, | |
266 const double& xmax, | |
267 const double& ymax) | |
268 { | |
269 // This is Skala algorithm for rectangles, "A new approach to line | |
270 // and line segment clipping in homogeneous coordinates" | |
271 // (2005). This is a direct, non-optimized translation of Algorithm | |
272 // 2 in the paper. | |
273 | |
274 static uint8_t tab1[16] = { 255 /* none */, | |
275 0, | |
276 0, | |
277 1, | |
278 1, | |
279 255 /* na */, | |
280 0, | |
281 2, | |
282 2, | |
283 0, | |
284 255 /* na */, | |
285 1, | |
286 1, | |
287 0, | |
288 0, | |
289 255 /* none */ }; | |
290 | |
291 | |
292 static uint8_t tab2[16] = { 255 /* none */, | |
293 3, | |
294 1, | |
295 3, | |
296 2, | |
297 255 /* na */, | |
298 2, | |
299 3, | |
300 3, | |
301 2, | |
302 255 /* na */, | |
303 2, | |
304 3, | |
305 1, | |
306 3, | |
307 255 /* none */ }; | |
308 | |
309 // Create the coordinates of the rectangle | |
310 Vector x[4]; | |
311 AssignVector(x[0], xmin, ymin, 1.0); | |
312 AssignVector(x[1], xmax, ymin, 1.0); | |
313 AssignVector(x[2], xmax, ymax, 1.0); | |
314 AssignVector(x[3], xmin, ymax, 1.0); | |
315 | |
316 // Move to homogoneous coordinates in 2D | |
317 Vector p; | |
318 | |
319 { | |
320 Vector a, b; | |
321 AssignVector(a, ax, ay, 1.0); | |
322 AssignVector(b, bx, by, 1.0); | |
323 CrossProduct(p, a, b); | |
324 } | |
325 | |
326 uint8_t c = 0; | |
327 | |
328 for (unsigned int k = 0; k < 4; k++) | |
329 { | |
330 if (boost::numeric::ublas::inner_prod(p, x[k]) >= 0) | |
331 { | |
332 c |= (1 << k); | |
333 } | |
334 } | |
335 | |
336 assert(c < 16); | |
337 | |
338 uint8_t i = tab1[c]; | |
339 uint8_t j = tab2[c]; | |
340 | |
341 if (i == 255 || j == 255) | |
342 { | |
343 return false; // No intersection | |
344 } | |
345 else | |
346 { | |
347 Vector a, b, e; | |
348 CrossProduct(e, x[i], x[(i + 1) % 4]); | |
349 CrossProduct(a, p, e); | |
350 CrossProduct(e, x[j], x[(j + 1) % 4]); | |
351 CrossProduct(b, p, e); | |
352 | |
353 // Go back to non-homogeneous coordinates | |
354 x1 = a[0] / a[2]; | |
355 y1 = a[1] / a[2]; | |
356 x2 = b[0] / b[2]; | |
357 y2 = b[1] / b[2]; | |
358 | |
359 return true; | |
360 } | |
361 } | |
32 | 362 |
363 | |
364 void GetPixelSpacing(double& spacingX, | |
365 double& spacingY, | |
366 const OrthancPlugins::IDicomDataset& dicom) | |
367 { | |
368 Vector v; | |
369 | |
370 if (ParseVector(v, dicom, OrthancPlugins::DICOM_TAG_PIXEL_SPACING)) | |
371 { | |
372 if (v.size() != 2 || | |
373 v[0] <= 0 || | |
374 v[1] <= 0) | |
375 { | |
376 LOG(ERROR) << "Bad value for PixelSpacing tag"; | |
377 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); | |
378 } | |
379 else | |
380 { | |
381 spacingX = v[0]; | |
382 spacingY = v[1]; | |
383 } | |
384 } | |
385 else | |
386 { | |
387 // The "PixelSpacing" is of type 1C: It could be absent, use | |
388 // default value in such a case | |
389 spacingX = 1; | |
390 spacingY = 1; | |
391 } | |
392 } | |
0 | 393 } |
394 } |