Mercurial > hg > orthanc-stl
changeset 35:ee3bc8f7df5b
reorganization
| author | Sebastien Jodogne <s.jodogne@gmail.com> |
|---|---|
| date | Thu, 04 Apr 2024 20:37:08 +0200 |
| parents | bee2017f3088 |
| children | 13698d34e059 |
| files | CMakeLists.txt Sources/Plugin.cpp Sources/StructureSet.cpp Sources/StructureSet.h Sources/StructureSetGeometry.cpp Sources/StructureSetGeometry.h |
| diffstat | 6 files changed, 611 insertions(+), 460 deletions(-) [+] |
line wrap: on
line diff
--- a/CMakeLists.txt Thu Apr 04 20:25:23 2024 +0200 +++ b/CMakeLists.txt Thu Apr 04 20:37:08 2024 +0200 @@ -207,6 +207,8 @@ Sources/Plugin.cpp Sources/STLToolbox.cpp Sources/StructurePolygon.cpp + Sources/StructureSet.cpp + Sources/StructureSetGeometry.cpp Sources/VTKToolbox.cpp Sources/Vector3D.cpp
--- a/Sources/Plugin.cpp Thu Apr 04 20:25:23 2024 +0200 +++ b/Sources/Plugin.cpp Thu Apr 04 20:37:08 2024 +0200 @@ -22,6 +22,8 @@ **/ +#include "StructureSetGeometry.h" +#include "StructureSet.h" #include "STLToolbox.h" #include "StructurePolygon.h" #include "VTKToolbox.h" @@ -156,466 +158,6 @@ #include <dcmtk/dcmdata/dcuid.h> -class StructureSet : public boost::noncopyable -{ -private: - std::vector<StructurePolygon*> polygons_; - std::string patientId_; - std::string studyInstanceUid_; - std::string seriesInstanceUid_; - std::string sopInstanceUid_; - bool hasFrameOfReferenceUid_; - std::string frameOfReferenceUid_; - -public: - explicit StructureSet(Orthanc::ParsedDicomFile& dicom) : - hasFrameOfReferenceUid_(false) - { - DcmDataset& dataset = *dicom.GetDcmtkObject().getDataset(); - patientId_ = STLToolbox::GetStringValue(dataset, DCM_PatientID); - studyInstanceUid_ = STLToolbox::GetStringValue(dataset, DCM_StudyInstanceUID); - seriesInstanceUid_ = STLToolbox::GetStringValue(dataset, DCM_SeriesInstanceUID); - sopInstanceUid_ = STLToolbox::GetStringValue(dataset, DCM_SOPInstanceUID); - - DcmSequenceOfItems* frame = NULL; - if (!dataset.findAndGetSequence(DCM_ReferencedFrameOfReferenceSequence, frame).good() || - frame == NULL) - { - throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); - } - - if (frame->card() == 1) - { - const char* v = NULL; - if (frame->getItem(0)->findAndGetString(DCM_FrameOfReferenceUID, v).good() && - v != NULL) - { - hasFrameOfReferenceUid_ = true; - frameOfReferenceUid_.assign(v); - } - } - - DcmSequenceOfItems* rois = NULL; - if (!dataset.findAndGetSequence(DCM_ROIContourSequence, rois).good() || - rois == NULL) - { - throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); - } - - std::vector<DcmSequenceOfItems*> contours(rois->card()); - size_t countPolygons = 0; - - for (unsigned long i = 0; i < rois->card(); i++) - { - DcmSequenceOfItems* contour = NULL; - if (!rois->getItem(i)->findAndGetSequence(DCM_ContourSequence, contour).good() || - contour == NULL) - { - throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); - } - else - { - contours[i] = contour; - countPolygons += contour->card(); - } - } - - polygons_.resize(countPolygons); - - size_t pos = 0; - for (unsigned long i = 0; i < contours.size(); i++) - { - for (unsigned long j = 0; j < contours[i]->card(); j++, pos++) - { - polygons_[pos] = new StructurePolygon(dicom, i, j); - } - } - - assert(pos == countPolygons); - } - - ~StructureSet() - { - for (size_t i = 0; i < polygons_.size(); i++) - { - assert(polygons_[i] != NULL); - delete polygons_[i]; - } - } - - const std::string& GetPatientId() const - { - return patientId_; - } - - const std::string& GetStudyInstanceUid() const - { - return studyInstanceUid_; - } - - const std::string& GetSeriesInstanceUid() const - { - return seriesInstanceUid_; - } - - const std::string& GetSopInstanceUid() const - { - return sopInstanceUid_; - } - - std::string HashStudy() const - { - Orthanc::DicomInstanceHasher hasher(patientId_, studyInstanceUid_, seriesInstanceUid_, sopInstanceUid_); - return hasher.HashStudy(); - } - - size_t GetPolygonsCount() const - { - return polygons_.size(); - } - - const StructurePolygon& GetPolygon(size_t i) const - { - if (i >= polygons_.size()) - { - throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); - } - else - { - assert(polygons_[i] != NULL); - return *polygons_[i]; - } - } - - bool HasFrameOfReferenceUid() const - { - return hasFrameOfReferenceUid_; - } - - const std::string& GetFrameOfReferenceUid() const - { - if (hasFrameOfReferenceUid_) - { - return frameOfReferenceUid_; - } - else - { - throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); - } - } - - // This static method is faster than constructing the full "StructureSet" object - static void ListStructuresNames(std::set<std::string>& target, - Orthanc::ParsedDicomFile& source) - { - target.clear(); - - DcmSequenceOfItems* sequence = NULL; - if (!source.GetDcmtkObject().getDataset()->findAndGetSequence(DCM_StructureSetROISequence, sequence).good() || - sequence == NULL) - { - throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); - } - - for (unsigned long i = 0; i < sequence->card(); i++) - { - DcmItem* item = sequence->getItem(i); - if (item == NULL) - { - throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); - } - else - { - target.insert(STLToolbox::GetStringValue(*item, DCM_ROIName)); - } - } - } -}; - - -class StructureSetGeometry : public boost::noncopyable -{ -private: - bool strict_; - Vector3D slicesNormal_; - double slicesSpacing_; - double minProjectionAlongNormal_; - double maxProjectionAlongNormal_; - size_t slicesCount_; - - bool 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; - } - } - -public: - StructureSetGeometry(const StructureSet& structures, - bool strict) : - strict_(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); - } - } - - const Vector3D& GetSlicesNormal() const - { - return slicesNormal_; - } - - double GetSlicesSpacing() const - { - return slicesSpacing_; - } - - double GetMinProjectionAlongNormal() const - { - return minProjectionAlongNormal_; - } - - double GetMaxProjectionAlongNormal() const - { - return maxProjectionAlongNormal_; - } - - bool 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; - } - } - - size_t GetSlicesCount() const - { - return slicesCount_; - } - - bool LookupSliceIndex(size_t& slice, - const StructurePolygon& polygon) const - { - double z; - return (ProjectAlongNormal(z, polygon) && - LookupProjectionIndex(slice, z)); - } -}; - - - - class XorFiller : public Orthanc::ImageProcessing::IPolygonFiller { private:
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Sources/StructureSet.cpp Thu Apr 04 20:37:08 2024 +0200 @@ -0,0 +1,170 @@ +/** + * 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 "StructureSet.h" + +#include "STLToolbox.h" + +#include <OrthancException.h> + +#include <dcmtk/dcmdata/dcdeftag.h> +#include <dcmtk/dcmdata/dcfilefo.h> + + +StructureSet::StructureSet(Orthanc::ParsedDicomFile& dicom) : + hasFrameOfReferenceUid_(false) +{ + DcmDataset& dataset = *dicom.GetDcmtkObject().getDataset(); + patientId_ = STLToolbox::GetStringValue(dataset, DCM_PatientID); + studyInstanceUid_ = STLToolbox::GetStringValue(dataset, DCM_StudyInstanceUID); + seriesInstanceUid_ = STLToolbox::GetStringValue(dataset, DCM_SeriesInstanceUID); + sopInstanceUid_ = STLToolbox::GetStringValue(dataset, DCM_SOPInstanceUID); + + DcmSequenceOfItems* frame = NULL; + if (!dataset.findAndGetSequence(DCM_ReferencedFrameOfReferenceSequence, frame).good() || + frame == NULL) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); + } + + if (frame->card() == 1) + { + const char* v = NULL; + if (frame->getItem(0)->findAndGetString(DCM_FrameOfReferenceUID, v).good() && + v != NULL) + { + hasFrameOfReferenceUid_ = true; + frameOfReferenceUid_.assign(v); + } + } + + DcmSequenceOfItems* rois = NULL; + if (!dataset.findAndGetSequence(DCM_ROIContourSequence, rois).good() || + rois == NULL) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); + } + + std::vector<DcmSequenceOfItems*> contours(rois->card()); + size_t countPolygons = 0; + + for (unsigned long i = 0; i < rois->card(); i++) + { + DcmSequenceOfItems* contour = NULL; + if (!rois->getItem(i)->findAndGetSequence(DCM_ContourSequence, contour).good() || + contour == NULL) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); + } + else + { + contours[i] = contour; + countPolygons += contour->card(); + } + } + + polygons_.resize(countPolygons); + + size_t pos = 0; + for (unsigned long i = 0; i < contours.size(); i++) + { + for (unsigned long j = 0; j < contours[i]->card(); j++, pos++) + { + polygons_[pos] = new StructurePolygon(dicom, i, j); + } + } + + assert(pos == countPolygons); +} + + +StructureSet::~StructureSet() +{ + for (size_t i = 0; i < polygons_.size(); i++) + { + assert(polygons_[i] != NULL); + delete polygons_[i]; + } +} + + +std::string StructureSet::HashStudy() const +{ + Orthanc::DicomInstanceHasher hasher(patientId_, studyInstanceUid_, seriesInstanceUid_, sopInstanceUid_); + return hasher.HashStudy(); +} + + +const StructurePolygon& StructureSet::GetPolygon(size_t i) const +{ + if (i >= polygons_.size()) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); + } + else + { + assert(polygons_[i] != NULL); + return *polygons_[i]; + } +} + + +const std::string& StructureSet::GetFrameOfReferenceUid() const +{ + if (hasFrameOfReferenceUid_) + { + return frameOfReferenceUid_; + } + else + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); + } +} + + +void StructureSet::ListStructuresNames(std::set<std::string>& target, + Orthanc::ParsedDicomFile& source) +{ + target.clear(); + + DcmSequenceOfItems* sequence = NULL; + if (!source.GetDcmtkObject().getDataset()->findAndGetSequence(DCM_StructureSetROISequence, sequence).good() || + sequence == NULL) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); + } + + for (unsigned long i = 0; i < sequence->card(); i++) + { + DcmItem* item = sequence->getItem(i); + if (item == NULL) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); + } + else + { + target.insert(STLToolbox::GetStringValue(*item, DCM_ROIName)); + } + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Sources/StructureSet.h Thu Apr 04 20:37:08 2024 +0200 @@ -0,0 +1,84 @@ +/** + * 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/>. + **/ + + +#pragma once + +#include "StructurePolygon.h" + +class StructureSet : public boost::noncopyable +{ +private: + std::vector<StructurePolygon*> polygons_; + std::string patientId_; + std::string studyInstanceUid_; + std::string seriesInstanceUid_; + std::string sopInstanceUid_; + bool hasFrameOfReferenceUid_; + std::string frameOfReferenceUid_; + +public: + explicit StructureSet(Orthanc::ParsedDicomFile& dicom); + + ~StructureSet(); + + const std::string& GetPatientId() const + { + return patientId_; + } + + const std::string& GetStudyInstanceUid() const + { + return studyInstanceUid_; + } + + const std::string& GetSeriesInstanceUid() const + { + return seriesInstanceUid_; + } + + const std::string& GetSopInstanceUid() const + { + return sopInstanceUid_; + } + + std::string HashStudy() const; + + size_t GetPolygonsCount() const + { + return polygons_.size(); + } + + const StructurePolygon& GetPolygon(size_t i) const; + + bool HasFrameOfReferenceUid() const + { + return hasFrameOfReferenceUid_; + } + + const std::string& GetFrameOfReferenceUid() const; + + // This static method is faster than constructing the full "StructureSet" object + static void ListStructuresNames(std::set<std::string>& target, + Orthanc::ParsedDicomFile& source); +};
--- /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)); +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Sources/StructureSetGeometry.h Thu Apr 04 20:37:08 2024 +0200 @@ -0,0 +1,76 @@ +/** + * 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/>. + **/ + + +#pragma once + +#include "StructureSet.h" + + +class StructureSetGeometry : public boost::noncopyable +{ +private: + Vector3D slicesNormal_; + double slicesSpacing_; + double minProjectionAlongNormal_; + double maxProjectionAlongNormal_; + size_t slicesCount_; + + bool LookupProjectionIndex(size_t& index, + double z) const; + +public: + StructureSetGeometry(const StructureSet& structures, + bool strict); + + const Vector3D& GetSlicesNormal() const + { + return slicesNormal_; + } + + double GetSlicesSpacing() const + { + return slicesSpacing_; + } + + double GetMinProjectionAlongNormal() const + { + return minProjectionAlongNormal_; + } + + double GetMaxProjectionAlongNormal() const + { + return maxProjectionAlongNormal_; + } + + bool ProjectAlongNormal(double& z, + const StructurePolygon& polygon) const; + + size_t GetSlicesCount() const + { + return slicesCount_; + } + + bool LookupSliceIndex(size_t& slice, + const StructurePolygon& polygon) const; +};