Mercurial > hg > orthanc-stone
changeset 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 | 050f01d7951b |
| children | 86a18d201e27 |
| files | Framework/Toolbox/CoordinateSystem3D.cpp Framework/Toolbox/CoordinateSystem3D.h Framework/Toolbox/DicomStructureSet.cpp Framework/Toolbox/GeometryToolbox.cpp Framework/Toolbox/GeometryToolbox.h |
| diffstat | 5 files changed, 103 insertions(+), 13 deletions(-) [+] |
line wrap: on
line diff
--- a/Framework/Toolbox/CoordinateSystem3D.cpp Fri Sep 20 13:59:36 2019 +0200 +++ b/Framework/Toolbox/CoordinateSystem3D.cpp Mon Sep 23 15:18:33 2019 +0200 @@ -168,6 +168,19 @@ return boost::numeric::ublas::inner_prod(point, normal_); } + void CoordinateSystem3D::ProjectPoint2(double& offsetX, double& offsetY, const Vector& point) const + { + // Project the point onto the slice + double projectionX,projectionY,projectionZ; + GeometryToolbox::ProjectPointOntoPlane2(projectionX, projectionY, projectionZ, point, normal_, origin_); + + // As the axes are orthonormal vectors thanks to + // CheckAndComputeNormal(), the following dot products give the + // offset of the origin of the slice wrt. the origin of the + // reference plane https://en.wikipedia.org/wiki/Vector_projection + offsetX = axisX_[0] * (projectionX - origin_[0]) + axisX_[1] * (projectionY - origin_[1]) + axisX_[2] * (projectionZ - origin_[2]); + offsetY = axisY_[0] * (projectionX - origin_[0]) + axisY_[1] * (projectionY - origin_[1]) + axisY_[2] * (projectionZ - origin_[2]); + } void CoordinateSystem3D::ProjectPoint(double& offsetX, double& offsetY,
--- a/Framework/Toolbox/CoordinateSystem3D.h Fri Sep 20 13:59:36 2019 +0200 +++ b/Framework/Toolbox/CoordinateSystem3D.h Mon Sep 23 15:18:33 2019 +0200 @@ -101,6 +101,13 @@ double& offsetY, const Vector& point) const; + /* + Alternated faster implementation (untested yet) + */ + void ProjectPoint2(double& offsetX, + double& offsetY, + const Vector& point) const; + bool IntersectSegment(Vector& p, const Vector& edgeFrom, const Vector& edgeTo) const;
--- a/Framework/Toolbox/DicomStructureSet.cpp Fri Sep 20 13:59:36 2019 +0200 +++ b/Framework/Toolbox/DicomStructureSet.cpp Mon Sep 23 15:18:33 2019 +0200 @@ -260,7 +260,7 @@ if (IsPointOnSliceIfAny(*it)) { double x, y; - geometry.ProjectPoint(x, y, *it); + geometry.ProjectPoint2(x, y, *it); extent_.AddPoint(x, y); } } @@ -301,7 +301,7 @@ } double x, y; - geometry_.ProjectPoint(x, y, slice.GetOrigin()); + geometry_.ProjectPoint2(x, y, slice.GetOrigin()); bool isOpposite; if (GeometryToolbox::IsParallelOrOpposite @@ -321,13 +321,13 @@ double xmax = -std::numeric_limits<double>::infinity(); double prevX, prevY; - geometry_.ProjectPoint(prevX, prevY, points_[points_.size() - 1]); + geometry_.ProjectPoint2(prevX, prevY, points_[points_.size() - 1]); for (size_t i = 0; i < points_.size(); i++) { // Reference: ../../Resources/Computations/IntersectSegmentAndHorizontalLine.py double curX, curY; - geometry_.ProjectPoint(curX, curY, points_[i]); + geometry_.ProjectPoint2(curX, curY, points_[i]); // if prev* and cur* are on opposite sides of y, this means that the // segment intersects the plane. @@ -365,8 +365,8 @@ sliceThickness_ / 2.0 * geometry_.GetNormal()); // then to the cutting plane geometry... - slice.ProjectPoint(x1, y1, p1); - slice.ProjectPoint(x2, y2, p2); + slice.ProjectPoint2(x1, y1, p1); + slice.ProjectPoint2(x2, y2, p2); return true; } } @@ -392,13 +392,13 @@ double ymax = -std::numeric_limits<double>::infinity(); double prevX, prevY; - geometry_.ProjectPoint(prevX, prevY, points_[points_.size() - 1]); + geometry_.ProjectPoint2(prevX, prevY, points_[points_.size() - 1]); for (size_t i = 0; i < points_.size(); i++) { // Reference: ../../Resources/Computations/IntersectSegmentAndVerticalLine.py double curX, curY; - geometry_.ProjectPoint(curX, curY, points_[i]); + geometry_.ProjectPoint2(curX, curY, points_[i]); if ((prevX < x && curX > x) || (prevX > x && curX < x)) @@ -424,8 +424,8 @@ Vector p2 = (geometry_.MapSliceToWorldCoordinates(x, ymax) - sliceThickness_ / 2.0 * geometry_.GetNormal()); - slice.ProjectPoint(x1, y1, p1); - slice.ProjectPoint(x2, y2, p2); + slice.ProjectPoint2(x1, y1, p1); + slice.ProjectPoint2(x2, y2, p2); // TODO WHY THIS??? y1 = -y1; @@ -838,7 +838,7 @@ p != polygon->GetPoints().end(); ++p) { double x, y; - slice.ProjectPoint(x, y, *p); + slice.ProjectPoint2(x, y, *p); polygons.back().push_back(Point2D(x, y)); } #else @@ -850,7 +850,7 @@ { Vector prev = points3D[0]; double prevX, prevY; - slice.ProjectPoint(prevX, prevY, prev); + slice.ProjectPoint2(prevX, prevY, prev); prev2D = Point2D(prevX, prevY); } @@ -859,7 +859,7 @@ { Vector next = points3D[ipt]; double nextX, nextY; - slice.ProjectPoint(nextX, nextY, next); + slice.ProjectPoint2(nextX, nextY, next); Point2D next2D(nextX, nextY); segments.push_back(std::pair<Point2D, Point2D>(prev2D, next2D)); prev2D = next2D;
--- 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,
--- a/Framework/Toolbox/GeometryToolbox.h Fri Sep 20 13:59:36 2019 +0200 +++ b/Framework/Toolbox/GeometryToolbox.h Mon Sep 23 15:18:33 2019 +0200 @@ -32,6 +32,16 @@ const Vector& planeNormal, const Vector& planeOrigin); + /* + Alternated faster implementation (untested yet) + */ + void ProjectPointOntoPlane2(double& resultX, + double& resultY, + double& resultZ, + const Vector& point, + const Vector& planeNormal, + const Vector& planeOrigin); + bool IsParallel(const Vector& u, const Vector& v);