Mercurial > hg > orthanc-stl
diff Sources/Plugin.cpp @ 35:ee3bc8f7df5b
reorganization
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Thu, 04 Apr 2024 20:37:08 +0200 |
parents | bee2017f3088 |
children | 13698d34e059 |
line wrap: on
line diff
--- 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: