Mercurial > hg > orthanc-stone
diff 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 |
line wrap: on
line diff
--- a/Framework/Toolbox/GeometryToolbox.cpp Fri Sep 20 13:59:36 2019 +0200 +++ b/Framework/Toolbox/GeometryToolbox.cpp Mon Sep 23 15:18:33 2019 +0200 @@ -52,7 +52,67 @@ result = boost::numeric::ublas::inner_prod(planeOrigin - point, n) * n + point; } + /* + undefined results if vector are not 3D + */ + void ProjectPointOntoPlane2( + double& resultX, + double& resultY, + double& resultZ, + const Vector& point, + const Vector& planeNormal, + const Vector& planeOrigin) + { + double pointX = point[0]; + double pointY = point[1]; + double pointZ = point[2]; + double planeNormalX = planeNormal[0]; + double planeNormalY = planeNormal[1]; + double planeNormalZ = planeNormal[2]; + + double planeOriginX = planeOrigin[0]; + double planeOriginY = planeOrigin[1]; + double planeOriginZ = planeOrigin[2]; + + double normSq = (planeNormalX * planeNormalX) + (planeNormalY * planeNormalY) + (planeNormalZ * planeNormalZ); + + // Algebraic form of line–plane intersection, where the line passes + // through "point" along the direction "normal" (thus, l == n) + // https://en.wikipedia.org/wiki/Line%E2%80%93plane_intersection#Algebraic_form + + if (LinearAlgebra::IsNear(1.0, normSq)) + { + double nX = planeNormalX; + double nY = planeNormalY; + double nZ = planeNormalZ; + + double prod = (planeOriginX - pointX) * nX + (planeOriginY - pointY) * nY + (planeOriginZ - pointZ) * nZ; + + resultX = prod * nX + pointX; + resultY = prod * nY + pointY; + resultZ = prod * nZ + pointZ; + } + else + { + double norm = sqrt(normSq); + if (LinearAlgebra::IsCloseToZero(norm)) + { + // Division by zero + throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); + } + double invNorm = 1.0 / norm; + double nX = planeNormalX * invNorm; + double nY = planeNormalY * invNorm; + double nZ = planeNormalZ * invNorm; + + double prod = (planeOriginX - pointX) * nX + (planeOriginY - pointY) * nY + (planeOriginZ - pointZ) * nZ; + + resultX = prod * nX + pointX; + resultY = prod * nY + pointY; + resultZ = prod * nZ + pointZ; + } + } bool IsParallelOrOpposite(bool& isOpposite, const Vector& u,