# HG changeset patch
# User Benjamin Golinvaux <bgo@osimis.io>
# Date 1569244713 -7200
# Node ID 53cc787bd7bce69639a6003ee20c8a2b304030e4
# Parent  050f01d7951ba5240fe2beddcf60ee2b3f520bda
- 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.

diff -r 050f01d7951b -r 53cc787bd7bc Framework/Toolbox/CoordinateSystem3D.cpp
--- 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,
diff -r 050f01d7951b -r 53cc787bd7bc Framework/Toolbox/CoordinateSystem3D.h
--- 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;
diff -r 050f01d7951b -r 53cc787bd7bc Framework/Toolbox/DicomStructureSet.cpp
--- 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;
diff -r 050f01d7951b -r 53cc787bd7bc Framework/Toolbox/GeometryToolbox.cpp
--- 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,
diff -r 050f01d7951b -r 53cc787bd7bc Framework/Toolbox/GeometryToolbox.h
--- 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);