--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Sources/StructureSetGeometry.cpp	Thu Apr 04 20:37:08 2024 +0200
@@ -0,0 +1,277 @@
+/**
+ * SPDX-FileCopyrightText: 2023-2024 Sebastien Jodogne, UCLouvain, Belgium
+ * SPDX-License-Identifier: GPL-3.0-or-later
+ */
+
+/**
+ * STL plugin for Orthanc
+ * Copyright (C) 2023-2024 Sebastien Jodogne, UCLouvain, Belgium
+ *
+ * This program is free software: you can redistribute it and/or
+ * modify it under the terms of the GNU 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
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ **/
+
+
+#include "StructureSetGeometry.h"
+
+#include "STLToolbox.h"
+
+#include <OrthancException.h>
+
+#include <list>
+
+
+bool StructureSetGeometry::LookupProjectionIndex(size_t& index,
+                                                 double z) const
+{
+  if (slicesCount_ == 0)
+  {
+    throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
+  }
+
+  if (z < minProjectionAlongNormal_ ||
+      z > maxProjectionAlongNormal_)
+  {
+    throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
+  }
+
+  assert(slicesSpacing_ > 0 &&
+         minProjectionAlongNormal_ < maxProjectionAlongNormal_);
+
+  double d = (z - minProjectionAlongNormal_) / slicesSpacing_;
+
+  if (STLToolbox::IsNear(d, round(d)))
+  {
+    if (d < 0.0 ||
+        d > static_cast<double>(slicesCount_) - 1.0)
+    {
+      throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
+    }
+    else
+    {
+      index = static_cast<size_t>(round(d));
+      return true;
+    }
+  }
+  else
+  {
+    return false;
+  }
+}
+
+
+StructureSetGeometry::StructureSetGeometry(const StructureSet& structures,
+                                           bool strict)
+{
+  bool isValid = false;
+
+  std::vector<double> projections;
+  projections.reserve(structures.GetPolygonsCount());
+
+  for (size_t i = 0; i < structures.GetPolygonsCount(); i++)
+  {
+    Vector3D normal;
+    if (structures.GetPolygon(i).IsCoplanar(normal))
+    {
+      // Initialize the normal of the whole volume, if need be
+      if (!isValid)
+      {
+        isValid = true;
+        slicesNormal_ = normal;
+      }
+
+      if (Vector3D::AreParallel(normal, slicesNormal_))
+      {
+        // This is a valid slice (it is parallel to the normal)
+        const Vector3D& point = structures.GetPolygon(i).GetPoint(0);
+        projections.push_back(Vector3D::DotProduct(point, slicesNormal_));
+      }
+      else
+      {
+        // RT-STRUCT with non-parallel slices
+        throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
+      }
+    }
+    else
+    {
+      // Ignore slices that are not coplanar
+    }
+  }
+
+  if (projections.empty())
+  {
+    throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented,
+                                    "Structure set without a valid geometry");
+  }
+
+  // Only keep unique projections
+
+  std::sort(projections.begin(), projections.end());
+  STLToolbox::RemoveDuplicateValues(projections);
+  assert(!projections.empty());
+
+  if (projections.size() == 1)
+  {
+    // Volume with one single slice
+    minProjectionAlongNormal_ = projections[0];
+    maxProjectionAlongNormal_ = projections[0];
+    slicesSpacing_ = 1;   // Arbitrary value
+    slicesCount_ = 1;
+    return;
+  }
+
+
+  // Compute the most probable spacing between the slices
+
+  {
+    std::vector<double> spacings;
+    spacings.resize(projections.size() - 1);
+
+    for (size_t i = 0; i < spacings.size(); i++)
+    {
+      spacings[i] = projections[i + 1] - projections[i];
+      assert(spacings[i] > 0);
+    }
+
+    std::sort(spacings.begin(), spacings.end());
+    STLToolbox::RemoveDuplicateValues(spacings);
+
+    if (spacings.empty())
+    {
+      throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
+    }
+
+    slicesSpacing_ = spacings[spacings.size() / 10];  // Take the 90% percentile of smallest spacings
+    assert(slicesSpacing_ > 0);
+  }
+
+
+  // Find the projection along the normal with the largest support
+
+  bool first = true;
+  size_t bestSupport;
+  double bestProjection;
+
+  std::list<size_t> candidates;
+  for (size_t i = 0; i < projections.size(); i++)
+  {
+    candidates.push_back(i);
+  }
+
+  while (!candidates.empty())
+  {
+    std::list<size_t> next;
+
+    size_t countSupport = 0;
+
+    std::list<size_t>::const_iterator it = candidates.begin();
+    size_t reference = *it;
+    it++;
+
+    while (it != candidates.end())
+    {
+      double d = (projections[*it] - projections[reference]) / slicesSpacing_;
+      if (STLToolbox::IsNear(d, round(d)))
+      {
+        countSupport ++;
+      }
+      else
+      {
+        next.push_back(*it);
+      }
+
+      it++;
+    }
+
+    if (first ||
+        countSupport > bestSupport)
+    {
+      first = false;
+      bestSupport = countSupport;
+      bestProjection = projections[reference];
+    }
+
+    if (strict &&
+        !next.empty())
+    {
+      throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat,
+                                      "Structure set with multiple support, which is not allowed in Strict mode");
+    }
+
+    candidates.swap(next);
+  }
+
+
+  // Compute the range of the projections
+
+  minProjectionAlongNormal_ = bestProjection;
+  maxProjectionAlongNormal_ = bestProjection;
+
+  for (size_t i = 0; i < projections.size(); i++)
+  {
+    double d = (projections[i] - bestProjection) / slicesSpacing_;
+    if (STLToolbox::IsNear(d, round(d)))
+    {
+      minProjectionAlongNormal_ = std::min(minProjectionAlongNormal_, projections[i]);
+      maxProjectionAlongNormal_ = std::max(maxProjectionAlongNormal_, projections[i]);
+    }
+  }
+
+  double d = (maxProjectionAlongNormal_ - minProjectionAlongNormal_) / slicesSpacing_;
+  if (STLToolbox::IsNear(d, round(d)))
+  {
+    slicesCount_ = static_cast<size_t>(round(d)) + 1;
+  }
+  else
+  {
+    throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
+  }
+
+
+  // Sanity check
+
+  size_t a, b;
+  if (!LookupProjectionIndex(a, minProjectionAlongNormal_) ||
+      !LookupProjectionIndex(b, maxProjectionAlongNormal_) ||
+      a != 0 ||
+      b + 1 != slicesCount_)
+  {
+    throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
+  }
+}
+
+
+bool StructureSetGeometry::ProjectAlongNormal(double& z,
+                                              const StructurePolygon& polygon) const
+{
+  Vector3D normal;
+  if (polygon.IsCoplanar(normal) &&
+      Vector3D::AreParallel(normal, slicesNormal_))
+  {
+    z = Vector3D::DotProduct(polygon.GetPoint(0), slicesNormal_);
+    return true;
+  }
+  else
+  {
+    return false;
+  }
+}
+
+
+bool StructureSetGeometry::LookupSliceIndex(size_t& slice,
+                                            const StructurePolygon& polygon) const
+{
+  double z;
+  return (ProjectAlongNormal(z, polygon) &&
+          LookupProjectionIndex(slice, z));
+}