Mercurial > hg > orthanc-stone

changeset 2076:990f396484b1

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
fix rendering of RT-DOSE with negative GridFrameOffsetVector
author Sebastien Jodogne <s.jodogne@gmail.com>
date Tue, 11 Jul 2023 15:58:16 +0200
parents d84bdcbd8bf1
children 07964689cb0b
files Applications/StoneWebViewer/WebAssembly/CMakeLists.txt OrthancStone/Resources/WebAssemblyUnitTests/CMakeLists.txt OrthancStone/Sources/Loaders/DicomVolumeLoader.cpp OrthancStone/Sources/Loaders/OrthancMultiframeVolumeLoader.cpp OrthancStone/Sources/Loaders/OrthancMultiframeVolumeLoader.h OrthancStone/Sources/Toolbox/DicomInstanceParameters.cpp OrthancStone/Sources/Toolbox/DicomInstanceParameters.h OrthancStone/UnitTestsSources/CMakeLists.txt OrthancStone/UnitTestsSources/DicomTests.cpp
diffstat 9 files changed, 196 insertions(+), 64 deletions(-) [+]
line wrap: on
line diff
--- a/Applications/StoneWebViewer/WebAssembly/CMakeLists.txt	Tue Jul 11 13:20:20 2023 +0200
+++ b/Applications/StoneWebViewer/WebAssembly/CMakeLists.txt	Tue Jul 11 15:58:16 2023 +0200
@@ -134,6 +134,8 @@
   StoneWebViewer.cpp
   )
 
+DefineSourceBasenameForTarget(StoneWebViewer)
+
 set_target_properties(StoneWebViewer
   PROPERTIES
   COMPILE_FLAGS "${WASM_FLAGS}"
--- a/OrthancStone/Resources/WebAssemblyUnitTests/CMakeLists.txt	Tue Jul 11 13:20:20 2023 +0200
+++ b/OrthancStone/Resources/WebAssemblyUnitTests/CMakeLists.txt	Tue Jul 11 15:58:16 2023 +0200
@@ -118,6 +118,8 @@
   ${ORTHANC_STONE_SOURCES}
   )
 
+DefineSourceBasenameForTarget(UnitTests)
+
   
 # Declare installation files for the module
 # ---------------------------------------------------------------
--- a/OrthancStone/Sources/Loaders/DicomVolumeLoader.cpp	Tue Jul 11 13:20:20 2023 +0200
+++ b/OrthancStone/Sources/Loaders/DicomVolumeLoader.cpp	Tue Jul 11 15:58:16 2023 +0200
@@ -57,7 +57,7 @@
       double spacing;
       if (parameters.GetSopClassUid() == SopClassUid_RTDose)
       {
-        if (!parameters.ComputeRegularSpacing(spacing))
+        if (!parameters.ComputeFrameOffsetsSpacing(spacing))
         {
           LOG(WARNING) << "Unable to compute the spacing in a RT-DOSE instance";
           spacing = frames.GetSpacingBetweenSlices();
--- a/OrthancStone/Sources/Loaders/OrthancMultiframeVolumeLoader.cpp	Tue Jul 11 13:20:20 2023 +0200
+++ b/OrthancStone/Sources/Loaders/OrthancMultiframeVolumeLoader.cpp	Tue Jul 11 15:58:16 2023 +0200
@@ -215,7 +215,11 @@
     switch (parameters.GetSopClassUid())
     {
       case SopClassUid_RTDose:
-        spacingZ = parameters.GetSliceThickness();
+        if (!parameters.ComputeFrameOffsetsSpacing(spacingZ))
+        {
+          throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented,
+                                          "Cannot load RT-DOSE with incorrect GridFrameOffsetVector (3004,000C)");
+        }
         break;
 
       default:
@@ -224,6 +228,8 @@
           "No support for multiframe instances with SOP class UID: " + GetSopClassUid(dicom));
     }
 
+    isReversedFrameOffsets_ = parameters.IsReversedFrameOffsets();
+
     const unsigned int width = parameters.GetImageInformation().GetWidth();
     const unsigned int height = parameters.GetImageInformation().GetHeight();
     const unsigned int depth = parameters.GetImageInformation().GetNumberOfFrames();
@@ -231,7 +237,7 @@
     {
       VolumeImageGeometry geometry;
       geometry.SetSizeInVoxels(width, height, depth);
-      geometry.SetAxialGeometry(parameters.GetGeometry());
+      geometry.SetAxialGeometry(parameters.GetMultiFrameGeometry());
       geometry.SetVoxelDimensions(parameters.GetPixelSpacingX(),
                                   parameters.GetPixelSpacingY(), spacingZ);
       volume_->Initialize(geometry, format, true /* Do compute range */);
@@ -241,8 +247,6 @@
 
     ScheduleFrameDownloads();
 
-
-
     BroadcastMessage(DicomVolumeImage::GeometryReadyMessage(*volume_));
   }
 
@@ -322,7 +326,17 @@
 
       for (unsigned int z = 0; z < depth; z++)
       {
-        ImageBuffer3D::SliceWriter writer(target, VolumeProjection_Axial, z);
+        unsigned targetZ;
+        if (isReversedFrameOffsets_)
+        {
+          targetZ = depth - 1 - z;
+        }
+        else
+        {
+          targetZ = z;
+        }
+
+        ImageBuffer3D::SliceWriter writer(target, VolumeProjection_Axial, targetZ);
 
         assert(writer.GetAccessor().GetWidth() == width &&
           writer.GetAccessor().GetHeight() == height);
@@ -554,6 +568,7 @@
     float outliersHalfRejectionRate) 
     : LoaderStateMachine(loadersContext)
     , volume_(volume)
+    , isReversedFrameOffsets_(false)
     , pixelDataLoaded_(false)
     , outliersHalfRejectionRate_(outliersHalfRejectionRate)
     , distributionRawMin_(0)
--- a/OrthancStone/Sources/Loaders/OrthancMultiframeVolumeLoader.h	Tue Jul 11 13:20:20 2023 +0200
+++ b/OrthancStone/Sources/Loaders/OrthancMultiframeVolumeLoader.h	Tue Jul 11 15:58:16 2023 +0200
@@ -49,6 +49,7 @@
     };
 
     boost::shared_ptr<DicomVolumeImage>  volume_;
+    bool                                 isReversedFrameOffsets_;
     std::string                          instanceId_;
     std::string                          transferSyntaxUid_;
     bool                                 pixelDataLoaded_;
--- a/OrthancStone/Sources/Toolbox/DicomInstanceParameters.cpp	Tue Jul 11 13:20:20 2023 +0200
+++ b/OrthancStone/Sources/Toolbox/DicomInstanceParameters.cpp	Tue Jul 11 15:58:16 2023 +0200
@@ -114,16 +114,18 @@
     }
 
 
-    bool sliceThicknessPresent = true;
-    if (!dicom.ParseDouble(sliceThickness_, Orthanc::DICOM_TAG_SLICE_THICKNESS))
+    if (dicom.ParseDouble(sliceThickness_, Orthanc::DICOM_TAG_SLICE_THICKNESS))
     {
+      hasSliceThickness_ = true;
+    }
+    else
+    {
+      hasSliceThickness_ = false;
+
       if (numberOfFrames_ > 1)
       {
         LOG(INFO) << "The (non-madatory) slice thickness information is missing in a multiframe image";
       }
-      
-      sliceThickness_ = 100.0 * std::numeric_limits<double>::epsilon();
-      sliceThicknessPresent = false;
     }
 
     hasPixelSpacing_ = GeometryToolbox::GetPixelSpacing(pixelSpacingX_, pixelSpacingY_, dicom);
@@ -139,49 +141,6 @@
     if (numberOfFrames_ > 1)
     {
       ExtractFrameOffsets(frameOffsets_, dicom, numberOfFrames_);
-
-      // if the slice thickness is unknown, we try to infer it from the sequence of grid frame offsets
-      // this only works if:
-      // - the first offset is 0.0 (case (a) of http://dicom.nema.org/medical/Dicom/2017c/output/chtml/part03/sect_C.8.8.3.2.html)
-      // - the offsets are all equal, to some small tolerance
-      // - the offsets is positive (increasing throughout the frames)
-      if (!sliceThicknessPresent)
-      {
-        if (frameOffsets_.size() >= 2)
-        {
-          double sliceThickness = std::abs(frameOffsets_[1] - frameOffsets_[0]);
-          
-          if (sliceThickness > 0)
-          {
-            bool sameSized = true;
-            
-            for (size_t i = 2; i < frameOffsets_.size(); ++i)
-            {
-              double currentThickness = std::abs(frameOffsets_[i] - frameOffsets_[i-1]);
-              if (!LinearAlgebra::IsNear(sliceThickness, currentThickness))
-              {
-                LOG(ERROR) << "Unable to extract slice thickness from GridFrameOffsetVector (3004,000C) (reason: varying spacing)";
-                sameSized = false;
-                break;
-              }
-            }
-            
-            if (sameSized)
-            {
-              sliceThickness_ = sliceThickness;
-              LOG(INFO) << "SliceThickness was not specified in the Dicom but was inferred from GridFrameOffsetVector (3004,000C).";
-            }
-          }
-        }
-        else
-        {
-          LOG(ERROR) << "Unable to extract slice thickness from GridFrameOffsetVector (3004,000C) (reason: GridFrameOffsetVector not present or too small)";
-        }
-      }
-    }
-    else
-    {
-      frameOffsets_.resize(0);
     }
 
     if (sopClassUid_ == SopClassUid_RTDose)
@@ -274,6 +233,19 @@
   }
 
 
+  double DicomInstanceParameters::GetSliceThickness() const
+  {
+    if (data_.hasSliceThickness_)
+    {
+      return data_.sliceThickness_;
+    }
+    else
+    {
+      throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
+    }
+  }
+
+
   const Orthanc::DicomImageInformation& DicomInstanceParameters::GetImageInformation() const
   {
     assert(tags_.get() != NULL);
@@ -715,7 +687,7 @@
   }
 
 
-  bool DicomInstanceParameters::ComputeRegularSpacing(double& spacing) const
+  bool DicomInstanceParameters::ComputeFrameOffsetsSpacing(double& spacing) const
   {
     if (data_.frameOffsets_.size() == 0)  // Not a RT-DOSE
     {
@@ -728,18 +700,32 @@
     }
     else
     {
-      assert(data_.frameOffsets_.size() == GetNumberOfFrames());
+      static double THRESHOLD = 0.001;
       
-      spacing = std::abs(data_.frameOffsets_[1] - data_.frameOffsets_[0]);
+      if (data_.frameOffsets_.size() != GetNumberOfFrames())
+      {
+        throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
+      }
+
+      double a = data_.frameOffsets_[1] - data_.frameOffsets_[0];
 
       for (size_t i = 1; i + 1 < data_.frameOffsets_.size(); i++)
       {
-        double s = std::abs(data_.frameOffsets_[i + 1] - data_.frameOffsets_[i]);
-        if (!LinearAlgebra::IsNear(spacing, s, 0.001))
+        double b = data_.frameOffsets_[i + 1] - data_.frameOffsets_[i];
+        if (!LinearAlgebra::IsNear(a, b, THRESHOLD))
         {
+          LOG(ERROR) << "Unable to extract slice thickness from GridFrameOffsetVector (3004,000C) (reason: varying spacing)";
           return false;
         }
       }
+
+      spacing = std::abs(a);
+
+      if (HasSliceThickness() &&
+          !LinearAlgebra::IsNear(spacing, GetSliceThickness(), THRESHOLD))
+      {
+        LOG(WARNING) << "SliceThickness and GridFrameOffsetVector (3004,000C) do not match";
+      }
       
       return true;
     }
@@ -816,4 +802,41 @@
       }
     }
   }
+
+
+  CoordinateSystem3D DicomInstanceParameters::GetMultiFrameGeometry() const
+  {
+    if (data_.frameOffsets_.empty())
+    {
+      return data_.geometry_;
+    }
+    else
+    {
+      assert(data_.frameOffsets_.size() == data_.numberOfFrames_);
+
+      double lowest = data_.frameOffsets_[0];
+      for (size_t i = 1; i < data_.frameOffsets_.size(); i++)
+      {
+        lowest = std::min(lowest, data_.frameOffsets_[i]);
+      }
+
+      return CoordinateSystem3D(
+        data_.geometry_.GetOrigin() + lowest * data_.geometry_.GetNormal(),
+        data_.geometry_.GetAxisX(),
+        data_.geometry_.GetAxisY());
+    }
+  }
+
+
+  bool DicomInstanceParameters::IsReversedFrameOffsets() const
+  {
+    if (data_.frameOffsets_.size() < 2)
+    {
+      return false;
+    }
+    else
+    {
+      return (data_.frameOffsets_[0] > data_.frameOffsets_[1]);
+    }
+  }
 }
--- a/OrthancStone/Sources/Toolbox/DicomInstanceParameters.h	Tue Jul 11 13:20:20 2023 +0200
+++ b/OrthancStone/Sources/Toolbox/DicomInstanceParameters.h	Tue Jul 11 15:58:16 2023 +0200
@@ -48,6 +48,7 @@
       unsigned int        numberOfFrames_;
       unsigned int        width_;
       unsigned int        height_;
+      bool                hasSliceThickness_;
       double              sliceThickness_;
       double              pixelSpacingX_;
       double              pixelSpacingY_;
@@ -143,11 +144,13 @@
       return data_.height_;
     }
 
-    double GetSliceThickness() const
+    bool HasSliceThickness() const
     {
-      return data_.sliceThickness_;
+      return data_.hasSliceThickness_;
     }
 
+    double GetSliceThickness() const;
+
     double GetPixelSpacingX() const
     {
       return data_.pixelSpacingX_;
@@ -234,7 +237,7 @@
     double ApplyRescale(double value) const;
 
     // Required for RT-DOSE
-    bool ComputeRegularSpacing(double& target) const;
+    bool ComputeFrameOffsetsSpacing(double& target) const;
 
     const std::string& GetFrameOfReferenceUid() const
     {
@@ -260,5 +263,9 @@
     {
       return data_.instanceNumber_;
     }
+
+    CoordinateSystem3D GetMultiFrameGeometry() const;
+
+    bool IsReversedFrameOffsets() const;
   };
 }
--- a/OrthancStone/UnitTestsSources/CMakeLists.txt	Tue Jul 11 13:20:20 2023 +0200
+++ b/OrthancStone/UnitTestsSources/CMakeLists.txt	Tue Jul 11 15:58:16 2023 +0200
@@ -59,6 +59,8 @@
   ${ORTHANC_STONE_SOURCES}
   )
 
+DefineSourceBasenameForTarget(UnitTests)
+
 add_custom_command(
   TARGET UnitTests
   POST_BUILD
--- a/OrthancStone/UnitTestsSources/DicomTests.cpp	Tue Jul 11 13:20:20 2023 +0200
+++ b/OrthancStone/UnitTestsSources/DicomTests.cpp	Tue Jul 11 15:58:16 2023 +0200
@@ -54,7 +54,8 @@
   ASSERT_EQ(1u, p->GetNumberOfFrames());
   ASSERT_EQ(0u, p->GetWidth());
   ASSERT_EQ(0u, p->GetHeight());
-  ASSERT_TRUE(OrthancStone::LinearAlgebra::IsCloseToZero(p->GetSliceThickness()));
+  ASSERT_FALSE(p->HasSliceThickness());
+  ASSERT_THROW(p->GetSliceThickness(), Orthanc::OrthancException);
   ASSERT_FLOAT_EQ(1, p->GetPixelSpacingX());
   ASSERT_FLOAT_EQ(1, p->GetPixelSpacingY());
   ASSERT_FALSE(p->GetGeometry().IsValid());
@@ -81,7 +82,7 @@
   ASSERT_DOUBLE_EQ(1.0, p->ApplyRescale(1.0));
 
   double s;
-  ASSERT_FALSE(p->ComputeRegularSpacing(s));
+  ASSERT_FALSE(p->ComputeFrameOffsetsSpacing(s));
   ASSERT_TRUE(p->GetFrameOfReferenceUid().empty());
 }
 
@@ -223,3 +224,82 @@
     }
   }
 }
+
+
+TEST(DicomInstanceParameters, ReverseFrameOffsetsGrid)
+{
+  Orthanc::DicomMap m;
+  SetupUids(m);
+
+  m.SetValue(Orthanc::DICOM_TAG_IMAGE_POSITION_PATIENT, "-276.611\\-274.463\\100", false);
+  m.SetValue(Orthanc::DICOM_TAG_IMAGE_ORIENTATION_PATIENT, "1\\0\\0\\0\\1\\0", false);
+  m.SetValue(Orthanc::DICOM_TAG_NUMBER_OF_FRAMES, "126", false);
+  m.SetValue(Orthanc::DICOM_TAG_GRID_FRAME_OFFSET_VECTOR, "0\\-2\\-4\\-6\\-8\\-10\\-12\\-14\\-16\\-18\\-20\\-22\\-24\\-26\\-28\\-30\\-32\\-34\\-36\\-38\\-40\\-42\\-44\\-46\\-48\\-50\\-52\\-54\\-56\\-58\\-60\\-62\\-64\\-66\\-68\\-70\\-72\\-74\\-76\\-78\\-80\\-82\\-84\\-86\\-88\\-90\\-92\\-94\\-96\\-98\\-100\\-102\\-104\\-106\\-108\\-110\\-112\\-114\\-116\\-118\\-120\\-122\\-124\\-126\\-128\\-130\\-132\\-134\\-136\\-138\\-140\\-142\\-144\\-146\\-148\\-150\\-152\\-154\\-156\\-158\\-160\\-162\\-164\\-166\\-168\\-170\\-172\\-174\\-176\\-178\\-180\\-182\\-184\\-186\\-188\\-190\\-192\\-194\\-196\\-198\\-200\\-202\\-204\\-206\\-208\\-210\\-212\\-214\\-216\\-218\\-220\\-222\\-224\\-226\\-228\\-230\\-232\\-234\\-236\\-238\\-240\\-242\\-244\\-246\\-248\\-250", false);
+
+  std::unique_ptr<OrthancStone::DicomInstanceParameters> p;
+  p.reset(OrthancStone::DicomInstanceParameters(m).Clone());
+
+  ASSERT_FALSE(p->HasSliceThickness());
+  ASSERT_THROW(p->GetSliceThickness(), Orthanc::OrthancException);
+
+  double s;
+  ASSERT_TRUE(p->ComputeFrameOffsetsSpacing(s));
+  ASSERT_DOUBLE_EQ(s, 2.0);
+  ASSERT_TRUE(p->IsReversedFrameOffsets());
+
+  ASSERT_DOUBLE_EQ(1, p->GetMultiFrameGeometry().GetAxisX() [0]);
+  ASSERT_DOUBLE_EQ(0, p->GetMultiFrameGeometry().GetAxisX() [1]);
+  ASSERT_DOUBLE_EQ(0, p->GetMultiFrameGeometry().GetAxisX() [2]);
+
+  ASSERT_DOUBLE_EQ(0, p->GetMultiFrameGeometry().GetAxisY() [0]);
+  ASSERT_DOUBLE_EQ(1, p->GetMultiFrameGeometry().GetAxisY() [1]);
+  ASSERT_DOUBLE_EQ(0, p->GetMultiFrameGeometry().GetAxisY() [2]);
+
+  ASSERT_DOUBLE_EQ(0, p->GetMultiFrameGeometry().GetNormal() [0]);
+  ASSERT_DOUBLE_EQ(0, p->GetMultiFrameGeometry().GetNormal() [1]);
+  ASSERT_DOUBLE_EQ(1, p->GetMultiFrameGeometry().GetNormal() [2]);
+
+  ASSERT_DOUBLE_EQ(-276.611, p->GetMultiFrameGeometry().GetOrigin() [0]);
+  ASSERT_DOUBLE_EQ(-274.463, p->GetMultiFrameGeometry().GetOrigin() [1]);
+  ASSERT_DOUBLE_EQ(100.0 - 250.0, p->GetMultiFrameGeometry().GetOrigin() [2]);
+}
+
+
+TEST(DicomInstanceParameters, StandardFrameOffsetsGrid)
+{
+  Orthanc::DicomMap m;
+  SetupUids(m);
+
+  m.SetValue(Orthanc::DICOM_TAG_SLICE_THICKNESS, "2", false);
+  m.SetValue(Orthanc::DICOM_TAG_IMAGE_POSITION_PATIENT, "-205.2157\\-388.4679\\-120", false);
+  m.SetValue(Orthanc::DICOM_TAG_IMAGE_ORIENTATION_PATIENT, "1\\0\\0\\0\\1\\0", false);
+  m.SetValue(Orthanc::DICOM_TAG_NUMBER_OF_FRAMES, "155", false);
+  m.SetValue(Orthanc::DICOM_TAG_GRID_FRAME_OFFSET_VECTOR, "0\\2\\4\\6\\8\\10\\12\\14\\16\\18\\20\\22\\24\\26\\28\\30\\32\\34\\36\\38\\40\\42\\44\\46\\48\\50\\52\\54\\56\\58\\60\\62\\64\\66\\68\\70\\72\\74\\76\\78\\80\\82\\84\\86\\88\\90\\92\\94\\96\\98\\100\\102\\104\\106\\108\\110\\112\\114\\116\\118\\120\\122\\124\\126\\128\\130\\132\\134\\136\\138\\140\\142\\144\\146\\148\\150\\152\\154\\156\\158\\160\\162\\164\\166\\168\\170\\172\\174\\176\\178\\180\\182\\184\\186\\188\\190\\192\\194\\196\\198\\200\\202\\204\\206\\208\\210\\212\\214\\216\\218\\220\\222\\224\\226\\228\\230\\232\\234\\236\\238\\240\\242\\244\\246\\248\\250\\252\\254\\256\\258\\260\\262\\264\\266\\268\\270\\272\\274\\276\\278\\280\\282\\284\\286\\288\\290\\292\\294\\296\\298\\300\\302\\304\\306\\308", false);
+
+  std::unique_ptr<OrthancStone::DicomInstanceParameters> p;
+  p.reset(OrthancStone::DicomInstanceParameters(m).Clone());
+
+  ASSERT_TRUE(p->HasSliceThickness());
+  ASSERT_DOUBLE_EQ(2.0, p->GetSliceThickness());
+
+  double s;
+  ASSERT_TRUE(p->ComputeFrameOffsetsSpacing(s));
+  ASSERT_DOUBLE_EQ(s, 2.0);
+  ASSERT_FALSE(p->IsReversedFrameOffsets());
+
+  ASSERT_DOUBLE_EQ(1, p->GetMultiFrameGeometry().GetAxisX() [0]);
+  ASSERT_DOUBLE_EQ(0, p->GetMultiFrameGeometry().GetAxisX() [1]);
+  ASSERT_DOUBLE_EQ(0, p->GetMultiFrameGeometry().GetAxisX() [2]);
+
+  ASSERT_DOUBLE_EQ(0, p->GetMultiFrameGeometry().GetAxisY() [0]);
+  ASSERT_DOUBLE_EQ(1, p->GetMultiFrameGeometry().GetAxisY() [1]);
+  ASSERT_DOUBLE_EQ(0, p->GetMultiFrameGeometry().GetAxisY() [2]);
+
+  ASSERT_DOUBLE_EQ(0, p->GetMultiFrameGeometry().GetNormal() [0]);
+  ASSERT_DOUBLE_EQ(0, p->GetMultiFrameGeometry().GetNormal() [1]);
+  ASSERT_DOUBLE_EQ(1, p->GetMultiFrameGeometry().GetNormal() [2]);
+
+  ASSERT_DOUBLE_EQ(-205.2157, p->GetMultiFrameGeometry().GetOrigin() [0]);
+  ASSERT_DOUBLE_EQ(-388.4679, p->GetMultiFrameGeometry().GetOrigin() [1]);
+  ASSERT_DOUBLE_EQ(-120.0, p->GetMultiFrameGeometry().GetOrigin() [2]);
+}