Mercurial > hg > orthanc-stone
comparison Framework/Toolbox/GeometryToolbox.cpp @ 1013:53cc787bd7bc toa2019092301
- 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 { |