Mercurial > hg > orthanc-stone

comparison Framework/Toolbox/GeometryToolbox.cpp @ 1013:53cc787bd7bc toa2019092301

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
- Added an optimized ProjectPoint2 to CoordinateSystem3D. It has *not* replaced the ProjectPoint method because more tests need to be written. - ProjectPointOntoPlane2 is a faster version of GeometryToolbox::ProjectPointOntoPlane. Same remark as above. - DicomStructureSet.cpp now uses this optimized call.
author Benjamin Golinvaux <bgo@osimis.io>
date Mon, 23 Sep 2019 15:18:33 +0200
parents b70e9be013e4
children 2d8ab34c8c91
comparison
equal deleted inserted replaced
1012:050f01d7951b 1013:53cc787bd7bc
50 // through "point" along the direction "normal" (thus, l == n) 50 // through "point" along the direction "normal" (thus, l == n)
51 // https://en.wikipedia.org/wiki/Line%E2%80%93plane_intersection#Algebraic_form 51 // https://en.wikipedia.org/wiki/Line%E2%80%93plane_intersection#Algebraic_form
52 result = boost::numeric::ublas::inner_prod(planeOrigin - point, n) * n + point; 52 result = boost::numeric::ublas::inner_prod(planeOrigin - point, n) * n + point;
53 } 53 }
54 54
55 55 /*
56 undefined results if vector are not 3D
57 */
58 void ProjectPointOntoPlane2(
59 double& resultX,
60 double& resultY,
61 double& resultZ,
62 const Vector& point,
63 const Vector& planeNormal,
64 const Vector& planeOrigin)
65 {
66 double pointX = point[0];
67 double pointY = point[1];
68 double pointZ = point[2];
69
70 double planeNormalX = planeNormal[0];
71 double planeNormalY = planeNormal[1];
72 double planeNormalZ = planeNormal[2];
73
74 double planeOriginX = planeOrigin[0];
75 double planeOriginY = planeOrigin[1];
76 double planeOriginZ = planeOrigin[2];
77
78 double normSq = (planeNormalX * planeNormalX) + (planeNormalY * planeNormalY) + (planeNormalZ * planeNormalZ);
79
80 // Algebraic form of line–plane intersection, where the line passes
81 // through "point" along the direction "normal" (thus, l == n)
82 // https://en.wikipedia.org/wiki/Line%E2%80%93plane_intersection#Algebraic_form
83
84 if (LinearAlgebra::IsNear(1.0, normSq))
85 {
86 double nX = planeNormalX;
87 double nY = planeNormalY;
88 double nZ = planeNormalZ;
89
90 double prod = (planeOriginX - pointX) * nX + (planeOriginY - pointY) * nY + (planeOriginZ - pointZ) * nZ;
91
92 resultX = prod * nX + pointX;
93 resultY = prod * nY + pointY;
94 resultZ = prod * nZ + pointZ;
95 }
96 else
97 {
98 double norm = sqrt(normSq);
99 if (LinearAlgebra::IsCloseToZero(norm))
100 {
101 // Division by zero
102 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
103 }
104 double invNorm = 1.0 / norm;
105 double nX = planeNormalX * invNorm;
106 double nY = planeNormalY * invNorm;
107 double nZ = planeNormalZ * invNorm;
108
109 double prod = (planeOriginX - pointX) * nX + (planeOriginY - pointY) * nY + (planeOriginZ - pointZ) * nZ;
110
111 resultX = prod * nX + pointX;
112 resultY = prod * nY + pointY;
113 resultZ = prod * nZ + pointZ;
114 }
115 }
56 116
57 bool IsParallelOrOpposite(bool& isOpposite, 117 bool IsParallelOrOpposite(bool& isOpposite,
58 const Vector& u, 118 const Vector& u,
59 const Vector& v) 119 const Vector& v)
60 { 120 {