--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/OrthancStone/Resources/Graveyard/RTStructTentativeReimplementation-BGO/DicomStructurePolygon2.cpp	Tue Feb 01 08:38:32 2022 +0100
@@ -0,0 +1,307 @@
+/**
+ * Stone of Orthanc
+ * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics
+ * Department, University Hospital of Liege, Belgium
+ * Copyright (C) 2017-2022 Osimis S.A., Belgium
+ * Copyright (C) 2021-2022 Sebastien Jodogne, ICTEAM UCLouvain, Belgium
+ *
+ * This program is free software: you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation, either version 3 of
+ * the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this program. If not, see
+ * <http://www.gnu.org/licenses/>.
+ **/
+
+#ifdef BGO_ENABLE_DICOMSTRUCTURESETLOADER2
+
+#include "DicomStructurePolygon2.h"
+
+#include "../Toolbox/LinearAlgebra.h"
+
+#include <Logging.h>
+
+namespace OrthancStone
+{
+  void DicomStructurePolygon2::ComputeDependentProperties()
+  {
+    ORTHANC_ASSERT(state_ == Building);
+
+    for (size_t j = 0; j < points_.size(); ++j)
+    {
+      // TODO: move to AddPoint!
+      const Vector& p = points_[j];
+      if (p[0] < minX_)
+        minX_ = p[0];
+      if (p[0] > maxX_)
+        maxX_ = p[0];
+
+      if (p[1] < minY_)
+        minY_ = p[1];
+      if (p[1] > maxY_)
+        maxY_ = p[1];
+
+      if (p[2] < minZ_)
+        minZ_ = p[2];
+      if (p[2] > maxZ_)
+        maxZ_ = p[2];
+    }
+
+    if (LinearAlgebra::IsNear(minX_, maxX_))
+    {
+      LinearAlgebra::AssignVector(normal_, 1, 0, 0);
+      //ORTHANC_ASSERT(!LinearAlgebra::IsNear(minX, maxX));
+      ORTHANC_ASSERT(!LinearAlgebra::IsNear(minY_, maxY_));
+      ORTHANC_ASSERT(!LinearAlgebra::IsNear(minZ_, maxZ_));
+    }
+    else if (LinearAlgebra::IsNear(minY_, maxY_))
+    {
+      LinearAlgebra::AssignVector(normal_, 0, 1, 0);
+      ORTHANC_ASSERT(!LinearAlgebra::IsNear(minX_, maxX_));
+      ORTHANC_ASSERT(!LinearAlgebra::IsNear(minZ_, maxZ_));
+    }
+    else if (LinearAlgebra::IsNear(minZ_, maxZ_))
+    {
+      LinearAlgebra::AssignVector(normal_, 0, 0, 1);
+      ORTHANC_ASSERT(!LinearAlgebra::IsNear(minX_, maxX_));
+      ORTHANC_ASSERT(!LinearAlgebra::IsNear(minY_, maxY_));
+    }
+    else
+    {
+      LOG(ERROR) << "The contour is not coplanar and not parallel to any axis.";
+      throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
+    }
+    state_ = Valid;
+  }
+
+
+  void DicomStructurePolygon2::ProjectOnConstantPlane(
+    std::vector<ScenePoint2D>& intersections, const CoordinateSystem3D& plane) const
+  {
+    // the plane can either have constant X, or constant Y.
+    // - for constant Z planes, use the ProjectOnParallelPlane method
+    // - other type of planes are not supported
+
+    // V is the coordinate that is constant in the plane
+    double planeV = 0.0;
+
+    // if true, then "u" in the code is "x" and "v" is "y". 
+    // (v is constant in the plane)
+    bool uvxy = false;
+
+    size_t uindex = static_cast<size_t>(-1);
+    size_t vindex = static_cast<size_t>(-1);
+
+    ORTHANC_ASSERT(LinearAlgebra::IsNear(plane.GetNormal()[2], 0.0));
+
+    if (LinearAlgebra::IsNear(plane.GetNormal()[1], 0.0))
+    {
+      // normal is 1,0,0 (or -1,0,0). 
+      // plane is constant X
+      uindex = 1;
+      vindex = 0;
+
+      uvxy = false;
+      planeV = plane.GetOrigin()[0];
+      if (planeV < minX_)
+        return;
+      if (planeV > maxX_)
+        return;
+    }
+    else if (LinearAlgebra::IsNear(plane.GetNormal()[0], 0.0))
+    {
+      // normal is 0,1,0 (or 0,-1,0). 
+      // plane is constant Y
+      uindex = 0;
+      vindex = 1;
+
+      uvxy = true;
+      planeV = plane.GetOrigin()[1];
+      if (planeV < minY_)
+        return;
+      if (planeV > maxY_)
+        return;
+    }
+    else
+    {
+      // if the following assertion(s) fail(s), it means the plane is NOT a constant-X or constant-Y plane
+      LOG(ERROR) << "Plane normal must be (a,0,0) or (0,a,0), with a == -1 or a == 1";
+      throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
+    }
+
+    size_t pointCount = GetPointCount();
+    if (pointCount >= 3)
+    {
+      // this vector will contain the coordinates of the intersection points
+      // between the plane and the polygon.
+      // these are expressed in the U coordinate, that is either X or Y, 
+      // depending upon the plane orientation
+      std::vector<double> uIntersections;
+
+      // we loop on the segments of the polygon (TODO: optimize)
+      // and we compute the intersection between each segment and the cut
+      // cutting plane (slice) has a constant X
+
+      for (size_t iPoint = 0; iPoint < pointCount; ++iPoint)
+      {
+        double u1 = points_[iPoint][uindex];
+        double v1 = points_[iPoint][vindex];
+
+        double u2 = 0;
+        double v2 = 0;
+
+        if (iPoint < pointCount - 1)
+        {
+          u2 = points_[iPoint + 1][uindex];
+          v2 = points_[iPoint + 1][vindex];
+        }
+        else
+        {
+          u2 = points_[0][uindex];
+          v2 = points_[0][vindex];
+        }
+
+        // Check if the segment intersects the plane
+        if ((std::min(v1, v2) <= planeV) && (std::max(v1, v2) >= planeV))
+        {
+          // special case: the segment is parallel to the plane but close to it
+          if (LinearAlgebra::IsNear(v1, v2))
+          {
+            // in that case, we choose to label both points as an intersection
+            double x, y;
+            plane.ProjectPoint(x, y, points_[iPoint]);
+            intersections.push_back(ScenePoint2D(x, y));
+
+            plane.ProjectPoint(x, y, points_[iPoint + 1]);
+            intersections.push_back(ScenePoint2D(x, y));
+          }
+          else
+          {
+            // we are looking for u so that (u,planeV) belongs to the segment
+            // let's define alpha = (u-u2)/(u1-u2) --> u = alpha*(u1-u2) + u2
+            // alpha = (v2-planeV)/(v2-v1)
+            // because the following two triangles are similar
+            // [ (planeY,x)  , (y2,x2), (planeY,x2) ] or
+            // [ (planeX,y)  , (x2,y2), (planeX,y2) ]
+            // and
+            // [ (y1    ,x1) , (y2,x2), (y1    ,x2) ] or
+            // [ (x1    ,y1) , (x2,y2), (x1    ,y2) ]
+
+            /*
+              void CoordinateSystem3D::ProjectPoint(double& offsetX,
+              double& offsetY,
+              const Vector& point) const
+            */
+            double alpha = (v2 - planeV) / (v2 - v1);
+
+            // get rid of numerical oddities
+            if (alpha < 0.0)
+              alpha = 0.0;
+            if (alpha > 1.0)
+              alpha = 1.0;
+            double u = alpha * (u1 - u2) + u2;
+
+            // here is the intersection in world coordinates
+            Vector intersection;
+            if(uvxy)
+              LinearAlgebra::AssignVector(intersection, u, planeV, minZ_);
+            else
+              LinearAlgebra::AssignVector(intersection, planeV, u, minZ_);
+
+            // and we convert it to plane coordinates
+            {
+              double xi, yi;
+              plane.ProjectPoint(xi, yi, intersection);
+
+              // we consider that the x axis is always parallel to the polygons
+              // TODO: is this hypothesis safe??????
+              uIntersections.insert(std::lower_bound(uIntersections.begin(), uIntersections.end(), xi), xi);
+            }
+          }
+        }
+      } // end of for (size_t iPoint = 0; iPoint < pointCount; ++iPoint)
+    
+      // now we convert the intersections to plane points
+      // we consider that the x axis is always parallel to the polygons
+      // TODO: same hypothesis as above: plane is perpendicular to polygons, 
+      // plane is parallel to the XZ (constant Y) or YZ (constant X) 3D planes
+      for (size_t i = 0; i < uIntersections.size(); ++i)
+      {
+        double x = uIntersections[i];
+        intersections.push_back(ScenePoint2D(x, minZ_));
+      }
+    } // end of if (pointCount >= 3)
+    else
+    {
+      LOG(ERROR) << "This polygon has " << pointCount << " vertices, which is less than 3 --> skipping";
+    }
+  } 
+
+  void DicomStructurePolygon2::ProjectOnParallelPlane(
+    std::vector< std::pair<ScenePoint2D, ScenePoint2D> >& segments, 
+    const CoordinateSystem3D& plane) const
+  {
+    if (points_.size() < 3)
+      return;
+
+    // the plane is horizontal
+    ORTHANC_ASSERT(LinearAlgebra::IsNear(plane.GetNormal()[0], 0.0));
+    ORTHANC_ASSERT(LinearAlgebra::IsNear(plane.GetNormal()[1], 0.0));
+
+    segments.clear();
+    segments.reserve(points_.size());
+    // since the returned values need to be expressed in the supplied coordinate
+    // system, we need to subtract origin_ from the returned points
+
+    double planeOriginX = plane.GetOrigin()[0];
+    double planeOriginY = plane.GetOrigin()[1];
+
+    // precondition: points_.size() >= 3
+    for (size_t j = 0; j < points_.size()-1; ++j)
+    {
+      // segment between point j and j+1
+
+      const Vector& point0 = GetPoint(j);
+      // subtract plane origin x and y
+      ScenePoint2D p0(point0[0] - planeOriginX, point0[1] - planeOriginY);
+    
+      const Vector& point1 = GetPoint(j+1);
+      // subtract plane origin x and y
+      ScenePoint2D p1(point1[0] - planeOriginX, point1[1] - planeOriginY);
+
+      segments.push_back(std::pair<ScenePoint2D, ScenePoint2D>(p0,p1));
+    }
+
+
+    // final segment 
+
+    const Vector& point0 = GetPoint(points_.size() - 1);
+    // subtract plane origin x and y
+    ScenePoint2D p0(point0[0] - planeOriginX, point0[1] - planeOriginY);
+
+    const Vector& point1 = GetPoint(0);
+    // subtract plane origin x and y
+    ScenePoint2D p1(point1[0] - planeOriginX, point1[1] - planeOriginY);
+
+    segments.push_back(std::pair<ScenePoint2D, ScenePoint2D>(p0, p1));
+  }
+
+  double DicomStructurePolygon2::GetZ() const
+  {
+    ORTHANC_ASSERT(LinearAlgebra::IsNear(normal_[0], 0.0));
+    ORTHANC_ASSERT(LinearAlgebra::IsNear(normal_[1], 0.0));
+    ORTHANC_ASSERT(LinearAlgebra::IsNear(minZ_, maxZ_));
+    return minZ_;
+  }
+}
+
+#endif 
+// BGO_ENABLE_DICOMSTRUCTURESETLOADER2
+