changeset 1161:19b1c8caade4 broker

fix sagittal geometry
author Sebastien Jodogne <s.jodogne@gmail.com>
date Tue, 19 Nov 2019 17:30:05 +0100
parents 59906485896f
children 709aa65aca17
files Framework/Toolbox/CoordinateSystem3D.cpp Framework/Toolbox/CoordinateSystem3D.h Framework/Toolbox/DicomInstanceParameters.cpp Framework/Toolbox/DicomInstanceParameters.h Framework/Volumes/DicomVolumeImageMPRSlicer.cpp Framework/Volumes/VolumeImageGeometry.cpp
diffstat 6 files changed, 74 insertions(+), 76 deletions(-) [+]
line wrap: on
line diff
--- a/Framework/Toolbox/CoordinateSystem3D.cpp	Tue Nov 19 12:52:07 2019 +0100
+++ b/Framework/Toolbox/CoordinateSystem3D.cpp	Tue Nov 19 17:30:05 2019 +0100
@@ -241,4 +241,14 @@
     return s;
   }
 
+
+  CoordinateSystem3D CoordinateSystem3D::NormalizeCuttingPlane(const CoordinateSystem3D& plane)
+  {
+    double ox, oy;
+    plane.ProjectPoint(ox, oy, LinearAlgebra::CreateVector(0, 0, 0));
+
+    CoordinateSystem3D normalized(plane);
+    normalized.SetOrigin(plane.MapSliceToWorldCoordinates(ox, oy));
+    return normalized;
+  }
 }
--- a/Framework/Toolbox/CoordinateSystem3D.h	Tue Nov 19 12:52:07 2019 +0100
+++ b/Framework/Toolbox/CoordinateSystem3D.h	Tue Nov 19 17:30:05 2019 +0100
@@ -22,6 +22,7 @@
 #pragma once
 
 #include "LinearAlgebra.h"
+#include "../Scene2D/ScenePoint2D.h"
 
 #include <Plugins/Samples/Common/IDicomDataset.h>
 
@@ -95,12 +96,24 @@
     Vector MapSliceToWorldCoordinates(double x,
                                       double y) const;
     
+    Vector MapSliceToWorldCoordinates(const ScenePoint2D& p) const
+    {
+      return MapSliceToWorldCoordinates(p.GetX(), p.GetY());
+    }
+    
     double ProjectAlongNormal(const Vector& point) const;
 
     void ProjectPoint(double& offsetX,
                       double& offsetY,
                       const Vector& point) const;
 
+    ScenePoint2D ProjectPoint(const Vector& point) const
+    {
+      double x, y;
+      ProjectPoint(x, y, point);
+      return ScenePoint2D(x, y);
+    }
+
     /*
     Alternated faster implementation (untested yet)
     */
@@ -120,5 +133,9 @@
     static bool ComputeDistance(double& distance,
                                 const CoordinateSystem3D& a,
                                 const CoordinateSystem3D& b);
+
+    // Normalize a cutting plane so that the origin (0,0,0) of the 3D
+    // world is mapped to the origin of its (x,y) coordinate system
+    static CoordinateSystem3D NormalizeCuttingPlane(const CoordinateSystem3D& plane);
   };
 }
--- a/Framework/Toolbox/DicomInstanceParameters.cpp	Tue Nov 19 12:52:07 2019 +0100
+++ b/Framework/Toolbox/DicomInstanceParameters.cpp	Tue Nov 19 17:30:05 2019 +0100
@@ -470,4 +470,33 @@
 
     return (value * factor + offset);
   }
+
+
+  bool DicomInstanceParameters::Data::ComputeRegularSpacing(double& spacing) const
+  {
+    if (frameOffsets_.size() == 0)  // Not a RT-DOSE
+    {
+      return false;
+    }
+    else if (frameOffsets_.size() == 1)
+    {
+      spacing = 1;   // Edge case: RT-DOSE with one single frame
+      return true;
+    }
+    else
+    {
+      spacing = std::abs(frameOffsets_[1] - frameOffsets_[0]);
+
+      for (size_t i = 1; i + 1 < frameOffsets_.size(); i++)
+      {
+        double s = frameOffsets_[i + 1] - frameOffsets_[i];
+        if (!LinearAlgebra::IsNear(spacing, s, 0.001))
+        {
+          return false;
+        }
+      }
+      
+      return true;
+    }
+  }
 }
--- a/Framework/Toolbox/DicomInstanceParameters.h	Tue Nov 19 12:52:07 2019 +0100
+++ b/Framework/Toolbox/DicomInstanceParameters.h	Tue Nov 19 17:30:05 2019 +0100
@@ -75,6 +75,8 @@
                                       bool useDouble) const;
 
       double ApplyRescale(double value) const;
+
+      bool ComputeRegularSpacing(double& target) const;
     };
 
     
@@ -209,6 +211,11 @@
       return data_.doseUnits_;
     }
 
+    void SetDoseGridScaling(double value)
+    {
+      data_.doseGridScaling_ = value;
+    }
+
     double GetDoseGridScaling() const
     {
       return data_.doseGridScaling_;
@@ -218,5 +225,11 @@
     {
       return data_.ApplyRescale(value);
     }
+
+    // Required for RT-DOSE
+    bool ComputeRegularSpacing(double& target) const
+    {
+      return data_.ComputeRegularSpacing(target);
+    }
   };
 }
--- a/Framework/Volumes/DicomVolumeImageMPRSlicer.cpp	Tue Nov 19 12:52:07 2019 +0100
+++ b/Framework/Volumes/DicomVolumeImageMPRSlicer.cpp	Tue Nov 19 17:30:05 2019 +0100
@@ -45,8 +45,7 @@
     revision_(volume_.GetRevision())
   {
     valid_ = (volume_.HasDicomParameters() &&
-              volume_.GetGeometry().DetectSlice(projection_, sliceIndex_, 
-                                                cuttingPlane));
+              volume_.GetGeometry().DetectSlice(projection_, sliceIndex_, cuttingPlane));
   }
 
 
@@ -82,21 +81,6 @@
       ImageBuffer3D::SliceReader reader(volume_.GetPixelData(), projection_, sliceIndex_);
       texture.reset(dynamic_cast<TextureBaseSceneLayer*>
                     (configurator->CreateTextureFromDicom(reader.GetAccessor(), parameters)));
-
-      // <DEBUG-BLOCK>
-#if 0
-      Orthanc::JpegWriter writer;
-      writer.SetQuality(60);
-      static int index = 0;
-      std::string filePath = "C:\\temp\\sliceReader_P";
-      filePath += boost::lexical_cast<std::string>(projection_);
-      filePath += "_I";
-      filePath += boost::lexical_cast<std::string>(index);
-      filePath += ".jpg";
-      index++;
-      writer.WriteToFile(filePath, reader.GetAccessor());
-#endif
-      // <END-OF-DEBUG-BLOCK>
     }
     
     const CoordinateSystem3D& system = volume_.GetGeometry().GetProjectionGeometry(projection_);
@@ -105,53 +89,11 @@
     cuttingPlane.ProjectPoint(x0, y0, system.GetOrigin());
     cuttingPlane.ProjectPoint(x1, y1, system.GetOrigin() + system.GetAxisX());
 
-    // <DEBUG-BLOCK>
-#if 0
     {
-      LOG(ERROR) << "+----------------------------------------------------+";
-      LOG(ERROR) << "| DicomVolumeImageMPRSlicer::Slice::CreateSceneLayer |";
-      LOG(ERROR) << "+----------------------------------------------------+";
-      std::string projectionString;
-      switch (projection_)
-      {
-      case VolumeProjection_Coronal:
-        projectionString = "CORONAL";
-        break;
-      case VolumeProjection_Axial:
-        projectionString = "CORONAL";
-        break;
-      case VolumeProjection_Sagittal:
-        projectionString = "SAGITTAL";
-        break;
-      default:
-        ORTHANC_ASSERT(false);
-      }
-      if(volume_.GetGeometry().GetDepth() == 200)
-        LOG(ERROR) << "| CT     IMAGE 512x512 with projection " << projectionString;
-      else
-        LOG(ERROR) << "| RTDOSE IMAGE NNNxNNN with projection " << projectionString;
-      LOG(ERROR) << "+----------------------------------------------------+";
-      LOG(ERROR) << "| cuttingPlane = " << cuttingPlane;
-      LOG(ERROR) << "| point to project = " << system.GetOrigin();
-      LOG(ERROR) << "| result = x0: " << x0 << " y0: " << y0;
-      LOG(ERROR) << "+----------------------- END ------------------------+";
+      double xz, yz;
+      cuttingPlane.ProjectPoint(xz, yz, LinearAlgebra::CreateVector(0, 0, 0));
+      texture->SetOrigin(x0 - xz, y0 - yz);
     }
-#endif
-    // <END-OF-DEBUG-BLOCK>
-
-#if 1 // BGO 2019-08-13
-    // The sagittal coordinate system has a Y vector going down. The displayed
-    // image (scene coords) has a Y vector pointing upwards (towards the patient 
-    // coord Z index)
-    // we need to flip the Y local coordinates to get the scene-coord offset.
-    // TODO: this is quite ugly. Isn't there a better way?
-    if(projection_ == VolumeProjection_Sagittal)
-      texture->SetOrigin(x0, -y0);
-    else
-      texture->SetOrigin(x0, y0);
-#else
-    texture->SetOrigin(x0, y0);
-#endif
 
     double dx = x1 - x0;
     double dy = y1 - y0;
@@ -164,19 +106,6 @@
     Vector tmp = volume_.GetGeometry().GetVoxelDimensions(projection_);
     texture->SetPixelSpacing(tmp[0], tmp[1]);
 
-    // <DEBUG-BLOCK>
-    {
-      //using std::endl;
-      //std::stringstream ss;
-      //ss << "DicomVolumeImageMPRSlicer::Slice::CreateSceneLayer | cuttingPlane = " << cuttingPlane << " | projection_ = " << projection_ << endl;
-      //ss << "volume_.GetGeometry().GetProjectionGeometry(projection_) = " << system << endl;
-      //ss << "cuttingPlane.ProjectPoint(x0, y0, system.GetOrigin()); --> | x0 = " << x0 << " | y0 = " << y0 << "| x1 = " << x1 << " | y1 = " << y1 << endl;
-      //ss << "volume_.GetGeometry() = " << volume_.GetGeometry() << endl;
-      //ss << "volume_.GetGeometry() = " << volume_.GetGeometry() << endl;
-      //LOG(ERROR) << ss.str();
-    }
-    // <END-OF-DEBUG-BLOCK>
-
     return texture.release();
   }
 
--- a/Framework/Volumes/VolumeImageGeometry.cpp	Tue Nov 19 12:52:07 2019 +0100
+++ b/Framework/Volumes/VolumeImageGeometry.cpp	Tue Nov 19 17:30:05 2019 +0100
@@ -39,7 +39,7 @@
     
     sagittalGeometry_ = CoordinateSystem3D(p,
                                            axialGeometry_.GetAxisY(),
-                                           axialGeometry_.GetNormal());
+                                           -axialGeometry_.GetNormal());
 
     Vector origin = (
       axialGeometry_.MapSliceToWorldCoordinates(-0.5 * voxelDimensions_[0],