Mercurial > hg > orthanc-stone

comparison OrthancStone/Sources/Toolbox/CoordinateSystem3D.cpp @ 1648:4a43106bc122

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
cross-hair starts to work
author Sebastien Jodogne <s.jodogne@gmail.com>
date Thu, 12 Nov 2020 19:27:08 +0100
parents 52b8b96cb55f
children af312e145980
comparison
equal deleted inserted replaced
1647:adff3cd78775 1648:4a43106bc122
27 27
28 #include <Logging.h> 28 #include <Logging.h>
29 #include <Toolbox.h> 29 #include <Toolbox.h>
30 #include <OrthancException.h> 30 #include <OrthancException.h>
31 31
32 #include <boost/math/constants/constants.hpp>
33
32 namespace OrthancStone 34 namespace OrthancStone
33 { 35 {
34 void CoordinateSystem3D::CheckAndComputeNormal() 36 void CoordinateSystem3D::CheckAndComputeNormal()
35 { 37 {
36 // DICOM expects normal vectors to define the axes: "The row and 38 // DICOM expects normal vectors to define the axes: "The row and
59 { 61 {
60 LinearAlgebra::CrossProduct(normal_, axisX_, axisY_); 62 LinearAlgebra::CrossProduct(normal_, axisX_, axisY_);
61 63
62 d_ = -(normal_[0] * origin_[0] + normal_[1] * origin_[1] + normal_[2] * origin_[2]); 64 d_ = -(normal_[0] * origin_[0] + normal_[1] * origin_[1] + normal_[2] * origin_[2]);
63 65
64 // Just a sanity check, it should be useless by construction 66 // Just a sanity check, it should be useless by construction (*)
65 assert(LinearAlgebra::IsNear(boost::numeric::ublas::norm_2(normal_), 1.0)); 67 assert(LinearAlgebra::IsNear(boost::numeric::ublas::norm_2(normal_), 1.0));
66 } 68 }
67 } 69 }
68 70
69 71
227 { 229 {
228 return GeometryToolbox::IntersectPlaneAndLine(p, normal_, d_, origin, direction); 230 return GeometryToolbox::IntersectPlaneAndLine(p, normal_, d_, origin, direction);
229 } 231 }
230 232
231 233
234 double CoordinateSystem3D::ComputeDistance(const Vector& p) const
235 {
236 /**
237 * "normal_" is an unit vector (*) => sqrt(a_1^2+a_2^2+a_3^2) = 1,
238 * and the denominator equals 1 by construction.
239 * https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_plane#Closest_point_and_distance_for_a_hyperplane_and_arbitrary_point
240 **/
241
242 return std::abs(boost::numeric::ublas::inner_prod(p, normal_) + d_);
243 }
244
245
232 bool CoordinateSystem3D::ComputeDistance(double& distance, 246 bool CoordinateSystem3D::ComputeDistance(double& distance,
233 const CoordinateSystem3D& a, 247 const CoordinateSystem3D& a,
234 const CoordinateSystem3D& b) 248 const CoordinateSystem3D& b)
235 { 249 {
236 bool opposite = false; // Ignored 250 bool opposite = false; // Ignored
263 277
264 CoordinateSystem3D normalized(plane); 278 CoordinateSystem3D normalized(plane);
265 normalized.SetOrigin(plane.MapSliceToWorldCoordinates(ox, oy)); 279 normalized.SetOrigin(plane.MapSliceToWorldCoordinates(ox, oy));
266 return normalized; 280 return normalized;
267 } 281 }
282
283
284 CoordinateSystem3D CoordinateSystem3D::CreateFromPlaneGeneralForm(double a,
285 double b,
286 double c,
287 double d)
288 {
289 /**
290 * "a*x + b*y + c*z + d = 0" => The un-normalized normal is vector
291 * (a,b,c).
292 **/
293
294 Vector normal;
295 LinearAlgebra::AssignVector(normal, a, b, c);
296
297
298 /**
299 * Choose the origin of plane, as the point that is closest to the
300 * origin of the axes (0,0,0).
301 * https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_plane#Restatement_using_linear_algebra
302 **/
303
304 double squaredNorm = a * a + b * b + c * c;
305 if (LinearAlgebra::IsCloseToZero(squaredNorm))
306 {
307 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadGeometry, "Singular matrix");
308 }
309
310 Vector origin = -d * normal / squaredNorm;
311
312
313 /**
314 * Select the X axis by computing a vector that is perpendicular
315 * to the normal.
316 *
317 * "Exactly 1 and only 1 of the bools get set; b0/b1/b2 gets set
318 * if dimension "i" has magnitude strictly less than all
319 * subsequent dimensions and not greater than all previous
320 * dimensions. We then have a unit vector with a single non-zero
321 * dimension that corresponds to a dimension of minimum magnitude
322 * in "normal". The cross product of this with "normal" is
323 * orthogonal to "normal" by definition of cross product. Consider
324 * now that the cross product is numerically unstable only when
325 * the two vectors are very closely aligned. Consider that our
326 * unit vector is large in only a single dimension and that that
327 * dimension corresponds to the dimension where "normal" was
328 * small. It's thus guaranteed to be loosely orthogonal to
329 * "normal" before taking the cross product, with least
330 * orthogonality in the case where all dimensions of "normal" are
331 * equal. In this least-orthogonal case, we're still quite
332 * orthogonal given that our unit vector has all but one dimension
333 * 0 whereas "normal" has all equal. We thus avoid the unstable
334 * case of taking the cross product of two nearly-aligned
335 * vectors." https://stackoverflow.com/a/43454629/881731
336 **/
337
338 bool b0 = (normal[0] < normal[1]) && (normal[0] < normal[2]);
339 bool b1 = (normal[1] <= normal[0]) && (normal[1] < normal[2]);
340 bool b2 = (normal[2] <= normal[0]) && (normal[2] <= normal[1]);
341 Vector swap = LinearAlgebra::CreateVector(b0 ? 1 : 0, b1 ? 1 : 0, b2 ? 1 : 0);
342
343 Vector axisX;
344 LinearAlgebra::CrossProduct(axisX, normal, swap);
345 LinearAlgebra::NormalizeVector(axisX);
346
347
348 /**
349 * The Y axis follows as the cross-product of the normal and the X
350 * axis.
351 **/
352
353 Vector axisY;
354 LinearAlgebra::CrossProduct(axisY, axisX, normal);
355 LinearAlgebra::NormalizeVector(axisY);
356
357 return CoordinateSystem3D(origin, axisX, axisY);
358 }
268 } 359 }