diff Framework/Toolbox/DicomStructureSet.cpp @ 122:e3433dabfb8d wasm

refactoring DicomStructureSet
author Sebastien Jodogne <s.jodogne@gmail.com>
date Fri, 06 Oct 2017 17:25:08 +0200
parents e66b2c757790
children 44fc253d4876
line wrap: on
line diff
--- a/Framework/Toolbox/DicomStructureSet.cpp	Wed Oct 04 17:53:47 2017 +0200
+++ b/Framework/Toolbox/DicomStructureSet.cpp	Fri Oct 06 17:25:08 2017 +0200
@@ -36,6 +36,7 @@
   static const OrthancPlugins::DicomTag DICOM_TAG_CONTOUR_GEOMETRIC_TYPE(0x3006, 0x0042);
   static const OrthancPlugins::DicomTag DICOM_TAG_CONTOUR_IMAGE_SEQUENCE(0x3006, 0x0016);
   static const OrthancPlugins::DicomTag DICOM_TAG_CONTOUR_SEQUENCE(0x3006, 0x0040);
+  static const OrthancPlugins::DicomTag DICOM_TAG_CONTOUR_DATA(0x3006, 0x0050);
   static const OrthancPlugins::DicomTag DICOM_TAG_NUMBER_OF_CONTOUR_POINTS(0x3006, 0x0046);
   static const OrthancPlugins::DicomTag DICOM_TAG_REFERENCED_SOP_INSTANCE_UID(0x0008, 0x1155);
   static const OrthancPlugins::DicomTag DICOM_TAG_ROI_CONTOUR_SEQUENCE(0x3006, 0x0039);
@@ -72,96 +73,81 @@
             GeometryToolbox::ParseVector(target, value));
   }
 
-  CoordinateSystem3D DicomStructureSet::
-  ExtractSliceGeometry(double& sliceThickness,
-                       OrthancPlugins::IOrthancConnection& orthanc,
-                       const OrthancPlugins::IDicomDataset& tags,
-                       size_t contourIndex,
-                       size_t sliceIndex)
-  {
-    using namespace OrthancPlugins;
-
-    size_t size;
-    if (!tags.GetSequenceSize(size, DicomPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, contourIndex,
-                                              DICOM_TAG_CONTOUR_SEQUENCE, sliceIndex,
-                                              DICOM_TAG_CONTOUR_IMAGE_SEQUENCE)) ||
-        size != 1)
-    {
-      throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);          
-    }
-
-    DicomDatasetReader reader(tags);
-    std::string parentUid = reader.GetMandatoryStringValue
-      (DicomPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, contourIndex,
-                 DICOM_TAG_CONTOUR_SEQUENCE, sliceIndex,
-                 DICOM_TAG_CONTOUR_IMAGE_SEQUENCE, 0,
-                 DICOM_TAG_REFERENCED_SOP_INSTANCE_UID));
-
-    Json::Value parentLookup;
-    MessagingToolbox::RestApiPost(parentLookup, orthanc, "/tools/lookup", parentUid);
 
-    if (parentLookup.type() != Json::arrayValue ||
-        parentLookup.size() != 1 ||
-        !parentLookup[0].isMember("Type") ||
-        !parentLookup[0].isMember("Path") ||
-        parentLookup[0]["Type"].type() != Json::stringValue ||
-        parentLookup[0]["ID"].type() != Json::stringValue ||
-        parentLookup[0]["Type"].asString() != "Instance")
+  void DicomStructureSet::Polygon::CheckPoint(const Vector& v)
+  {
+    if (hasSlice_)
     {
-      throw Orthanc::OrthancException(Orthanc::ErrorCode_UnknownResource);          
+      if (!GeometryToolbox::IsNear(GeometryToolbox::ProjectAlongNormal(v, normal_),
+                                   projectionAlongNormal_, 
+                                   sliceThickness_ / 2.0 /* in mm */))
+      {
+        LOG(ERROR) << "This RT-STRUCT contains a point that is off the slice of its instance";
+        throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);          
+      }
     }
-
-    Json::Value parentInstance;
-    MessagingToolbox::RestApiGet(parentInstance, orthanc, "/instances/" + parentLookup[0]["ID"].asString());
-
-    if (parentInstance.type() != Json::objectValue ||
-        !parentInstance.isMember("ParentSeries") ||
-        parentInstance["ParentSeries"].type() != Json::stringValue)
-    {
-      throw Orthanc::OrthancException(Orthanc::ErrorCode_NetworkProtocol);          
-    }
+  }
 
-    std::string parentSeriesId = parentInstance["ParentSeries"].asString();
-    bool isFirst = parentSeriesId_.empty();
 
-    if (isFirst)
-    {
-      parentSeriesId_ = parentSeriesId;
-    }
-    else if (parentSeriesId_ != parentSeriesId)
+  void DicomStructureSet::Polygon::AddPoint(const Vector& v)
+  {
+    CheckPoint(v);
+    points_.push_back(v);
+  }
+
+
+  bool DicomStructureSet::Polygon::UpdateReferencedSlice(const ReferencedSlices& slices)
+  {
+    if (hasSlice_)
     {
-      LOG(ERROR) << "This RT-STRUCT refers to several different series";
-      throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
-    }
-
-    FullOrthancDataset parentTags(orthanc, "/instances/" + parentLookup[0]["ID"].asString() + "/tags");
-    CoordinateSystem3D slice(parentTags);
-
-    Vector v;
-    if (ParseVector(v, parentTags, DICOM_TAG_SLICE_THICKNESS) &&
-        v.size() > 0)
-    {
-      sliceThickness = v[0];
+      return true;
     }
     else
     {
-      sliceThickness = 1;  // 1 mm by default
-    }
+      ReferencedSlices::const_iterator it = slices.find(sopInstanceUid_);
+      
+      if (it == slices.end())
+      {
+        return false;
+      }
+      else
+      {
+        const CoordinateSystem3D& geometry = it->second.geometry_;
+        
+        hasSlice_ = true;
+        normal_ = geometry.GetNormal();
+        projectionAlongNormal_ = GeometryToolbox::ProjectAlongNormal(geometry.GetOrigin(), normal_);
+        sliceThickness_ = it->second.thickness_;
 
-    if (isFirst)
-    {
-      normal_ = slice.GetNormal();
+        for (Points::const_iterator it = points_.begin(); it != points_.end(); ++it)
+        {
+          CheckPoint(*it);
+        }
+
+        return true;
+      }
     }
-    else if (!GeometryToolbox::IsParallel(normal_, slice.GetNormal()))
-    {
-      LOG(ERROR) << "Incompatible orientation of slices in this RT-STRUCT";
-      throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
-    }
-
-    return slice;
   }
 
 
+  bool DicomStructureSet::Polygon::IsOnSlice(const CoordinateSystem3D& slice) const
+  {
+    bool isOpposite;
+
+    if (points_.empty() ||
+        !hasSlice_ ||
+        !GeometryToolbox::IsParallelOrOpposite(isOpposite, slice.GetNormal(), normal_))
+    {
+      return false;
+    }
+    
+    double d = GeometryToolbox::ProjectAlongNormal(slice.GetOrigin(), normal_);
+
+    return (GeometryToolbox::IsNear(d, projectionAlongNormal_,
+                                    sliceThickness_ / 2.0));
+  }
+
+  
   const DicomStructureSet::Structure& DicomStructureSet::GetStructure(size_t index) const
   {
     if (index >= structures_.size())
@@ -173,22 +159,10 @@
   }
 
 
-  bool DicomStructureSet::IsPolygonOnSlice(const Polygon& polygon,
-                                           const CoordinateSystem3D& geometry) const
-  {
-    double d = boost::numeric::ublas::inner_prod(geometry.GetOrigin(), normal_);
-
-    return (GeometryToolbox::IsNear(d, polygon.projectionAlongNormal_, polygon.sliceThickness_ / 2.0) &&
-            !polygon.points_.empty());
-  }
-
-
-  DicomStructureSet::DicomStructureSet(OrthancPlugins::IOrthancConnection& orthanc,
-                                       const std::string& instanceId)
+  DicomStructureSet::DicomStructureSet(const OrthancPlugins::FullOrthancDataset& tags)
   {
     using namespace OrthancPlugins;
 
-    FullOrthancDataset tags(orthanc, "/instances/" + instanceId + "/tags");
     DicomDatasetReader reader(tags);
     
     size_t count, tmp;
@@ -204,13 +178,15 @@
     structures_.resize(count);
     for (size_t i = 0; i < count; i++)
     {
-      structures_[i].interpretation_ = reader.GetStringValue(DicomPath(DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE, i,
-                                                                       DICOM_TAG_RT_ROI_INTERPRETED_TYPE),
-                                                             "No interpretation");
+      structures_[i].interpretation_ = reader.GetStringValue
+        (DicomPath(DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE, i,
+                   DICOM_TAG_RT_ROI_INTERPRETED_TYPE),
+         "No interpretation");
 
-      structures_[i].name_ = reader.GetStringValue(DicomPath(DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE, i,
-                                                             DICOM_TAG_ROI_NAME),
-                                                   "No interpretation");
+      structures_[i].name_ = reader.GetStringValue
+        (DicomPath(DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE, i,
+                   DICOM_TAG_ROI_NAME),
+         "No interpretation");
 
       Vector color;
       if (ParseVector(color, tags, DicomPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i,
@@ -246,31 +222,45 @@
       {
         unsigned int countPoints;
 
-        if (!reader.GetUnsignedIntegerValue(countPoints, DicomPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i,
-                                                                   DICOM_TAG_CONTOUR_SEQUENCE, j,
-                                                                   DICOM_TAG_NUMBER_OF_CONTOUR_POINTS)))
+        if (!reader.GetUnsignedIntegerValue
+            (countPoints, DicomPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i,
+                                    DICOM_TAG_CONTOUR_SEQUENCE, j,
+                                    DICOM_TAG_NUMBER_OF_CONTOUR_POINTS)))
         {
           throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
         }
             
         //LOG(INFO) << "Parsing slice containing " << countPoints << " vertices";
 
-        std::string type = reader.GetMandatoryStringValue(DicomPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i,
-                                                                    DICOM_TAG_CONTOUR_SEQUENCE, j,
-                                                                    DICOM_TAG_CONTOUR_GEOMETRIC_TYPE));
+        std::string type = reader.GetMandatoryStringValue
+          (DicomPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i,
+                     DICOM_TAG_CONTOUR_SEQUENCE, j,
+                     DICOM_TAG_CONTOUR_GEOMETRIC_TYPE));
         if (type != "CLOSED_PLANAR")
         {
           LOG(ERROR) << "Cannot handle contour with geometry type: " << type;
           throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);          
         }
 
-        // The "CountourData" tag (3006,0050) is too large to be
-        // returned by the "/instances/{id}/tags" URI: Access it using
-        // the raw "/instances/{id}/content/{...}" endpoint
-        std::string slicesData;
-        orthanc.RestApiGet(slicesData, "/instances/" + instanceId + "/content/3006-0039/" +
-                           boost::lexical_cast<std::string>(i) + "/3006-0040/" +
-                           boost::lexical_cast<std::string>(j) + "/3006-0050");
+        size_t size;
+        if (!tags.GetSequenceSize(size, DicomPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i,
+                                                  DICOM_TAG_CONTOUR_SEQUENCE, j,
+                                                  DICOM_TAG_CONTOUR_IMAGE_SEQUENCE)) ||
+            size != 1)
+        {
+          throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);          
+        }
+
+        std::string sopInstanceUid = reader.GetMandatoryStringValue
+          (DicomPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i,
+                     DICOM_TAG_CONTOUR_SEQUENCE, j,
+                     DICOM_TAG_CONTOUR_IMAGE_SEQUENCE, 0,
+                     DICOM_TAG_REFERENCED_SOP_INSTANCE_UID));
+        
+        std::string slicesData = reader.GetMandatoryStringValue
+          (DicomPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i,
+                     DICOM_TAG_CONTOUR_SEQUENCE, j,
+                     DICOM_TAG_CONTOUR_DATA));
 
         Vector points;
         if (!GeometryToolbox::ParseVector(points, slicesData) ||
@@ -279,9 +269,7 @@
           throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);          
         }
 
-        Polygon polygon;
-        CoordinateSystem3D geometry = ExtractSliceGeometry(polygon.sliceThickness_, orthanc, tags, i, j);
-        polygon.projectionAlongNormal_ = geometry.ProjectAlongNormal(geometry.GetOrigin());
+        Polygon polygon(sopInstanceUid);
 
         for (size_t k = 0; k < countPoints; k++)
         {
@@ -289,27 +277,12 @@
           v[0] = points[3 * k];
           v[1] = points[3 * k + 1];
           v[2] = points[3 * k + 2];
-
-          if (!GeometryToolbox::IsNear(geometry.ProjectAlongNormal(v), 
-                                       polygon.projectionAlongNormal_, 
-                                       polygon.sliceThickness_ / 2.0 /* in mm */))
-          {
-            LOG(ERROR) << "This RT-STRUCT contains a point that is off the slice of its instance";
-            throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);          
-          }
-
-          polygon.points_.push_back(v);
+          polygon.AddPoint(v);
         }
 
         structures_[i].polygons_.push_back(polygon);
       }
     }
-
-    if (parentSeriesId_.empty() ||
-        normal_.size() != 3)
-    {
-      throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);          
-    }
   }
 
 
@@ -329,9 +302,9 @@
     for (Polygons::const_iterator polygon = structure.polygons_.begin();
          polygon != structure.polygons_.end(); ++polygon)
     {
-      if (!polygon->points_.empty())
+      if (!polygon->GetPoints().empty())
       {
-        center += polygon->points_.front() / n;
+        center += polygon->GetPoints().front() / n;
       }
     }
 
@@ -364,48 +337,37 @@
 
 
   void DicomStructureSet::Render(CairoContext& context,
-                                 const CoordinateSystem3D& slice) const
+                                 const CoordinateSystem3D& slice)
   {
     cairo_t* cr = context.GetObject();
 
-    for (Structures::const_iterator structure = structures_.begin();
+    for (Structures::iterator structure = structures_.begin();
          structure != structures_.end(); ++structure)
     {
-      if (structure->name_ != "SKIN" &&
-          structure->name_ != "HEART" &&
-          //structure->name_ != "CORD" &&
-          structure->name_ != "ESOPHAGUS" &&
-          structure->name_ != "LUNG_LT" &&
-          structure->name_ != "LUNG_RT" &&
-          structure->name_ != "GTV_EXH_PRIMARY" &&
-          structure->name_ != "GTV_INH_PRIMARY" &&
-          structure->name_ != "GTV_PRIMARY")
-      {
-        continue;
-      }
-      
-      for (Polygons::const_iterator polygon = structure->polygons_.begin();
+      for (Polygons::iterator polygon = structure->polygons_.begin();
            polygon != structure->polygons_.end(); ++polygon)
       {
-        if (IsPolygonOnSlice(*polygon, slice))
+        polygon->UpdateReferencedSlice(referencedSlices_);
+          
+        if (polygon->IsOnSlice(slice))
         {
           context.SetSourceColor(structure->red_, structure->green_, structure->blue_);
 
-          Points::const_iterator p = polygon->points_.begin();
+          Points::const_iterator p = polygon->GetPoints().begin();
 
           double x, y;
           slice.ProjectPoint(x, y, *p);
           cairo_move_to(cr, x, y);
           ++p;
             
-          while (p != polygon->points_.end())
+          while (p != polygon->GetPoints().end())
           {
             slice.ProjectPoint(x, y, *p);
             cairo_line_to(cr, x, y);
             ++p;
           }
 
-          slice.ProjectPoint(x, y, *polygon->points_.begin());
+          slice.ProjectPoint(x, y, *polygon->GetPoints().begin());
           cairo_line_to(cr, x, y);
 
           cairo_stroke(cr);
@@ -413,4 +375,161 @@
       }
     }
   }
+
+
+  void DicomStructureSet::GetReferencedInstances(std::set<std::string>& instances)
+  {
+    for (Structures::const_iterator structure = structures_.begin();
+         structure != structures_.end(); ++structure)
+    {
+      for (Polygons::const_iterator polygon = structure->polygons_.begin();
+           polygon != structure->polygons_.end(); ++polygon)
+      {
+        instances.insert(polygon->GetSopInstanceUid());
+      }
+    }
+  }
+
+
+  void DicomStructureSet::AddReferencedSlice(const std::string& sopInstanceUid,
+                                             const std::string& seriesInstanceUid,
+                                             const CoordinateSystem3D& geometry,
+                                             double thickness)
+  {
+    if (referencedSlices_.find(sopInstanceUid) != referencedSlices_.end())
+    {
+      // This geometry is already known
+      throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
+    }
+    else
+    {
+      if (thickness < 0)
+      {
+        throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
+      }
+        
+      if (!referencedSlices_.empty())
+      {
+        const ReferencedSlice& reference = referencedSlices_.begin()->second;
+
+        if (reference.seriesInstanceUid_ != seriesInstanceUid)
+        {
+          LOG(ERROR) << "This RT-STRUCT refers to several different series";
+          throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
+        }
+
+        if (!GeometryToolbox::IsParallel(reference.geometry_.GetNormal(), geometry.GetNormal()))
+        {
+          LOG(ERROR) << "The slices in this RT-STRUCT are not parallel";
+          throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
+        }
+      }
+        
+      referencedSlices_[sopInstanceUid] = ReferencedSlice(seriesInstanceUid, geometry, thickness);
+    }
+  }
+
+
+  void DicomStructureSet::AddReferencedSlice(const Orthanc::DicomMap& dataset)
+  {
+    CoordinateSystem3D slice(dataset);
+
+    double thickness = 1;  // 1 mm by default
+
+    std::string s;
+    Vector v;
+    if (dataset.CopyToString(s, Orthanc::DICOM_TAG_SLICE_THICKNESS, false) &&
+        GeometryToolbox::ParseVector(v, s) &&
+        v.size() > 0)
+    {
+      thickness = v[0];
+    }
+
+    std::string instance, series;
+    if (dataset.CopyToString(instance, Orthanc::DICOM_TAG_SOP_INSTANCE_UID, false) &&
+        dataset.CopyToString(series, Orthanc::DICOM_TAG_SERIES_INSTANCE_UID, false))
+    {
+      AddReferencedSlice(instance, series, slice, thickness);
+    }
+    else
+    {
+      throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
+    }
+  }
+
+
+  void DicomStructureSet::CheckReferencedSlices()
+  {
+    for (Structures::iterator structure = structures_.begin();
+         structure != structures_.end(); ++structure)
+    {
+      for (Polygons::iterator polygon = structure->polygons_.begin();
+           polygon != structure->polygons_.end(); ++polygon)
+      {
+        if (!polygon->UpdateReferencedSlice(referencedSlices_))
+        {
+          LOG(ERROR) << "Missing information about referenced instance: "
+                     << polygon->GetSopInstanceUid();
+          throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
+        }
+      }
+    }
+  }
+
+
+  Vector DicomStructureSet::GetNormal() const
+  {
+    if (!referencedSlices_.empty())
+    {
+      Vector v;
+      GeometryToolbox::AssignVector(v, 0, 0, 1);
+      return v;
+    }
+    else
+    {
+      return referencedSlices_.begin()->second.geometry_.GetNormal();
+    }
+  }
+
+  
+  DicomStructureSet* DicomStructureSet::SynchronousLoad(OrthancPlugins::IOrthancConnection& orthanc,
+                                                        const std::string& instanceId)
+  {
+    const std::string uri = "/instances/" + instanceId + "/tags?ignore-length=3006-0050";
+    OrthancPlugins::FullOrthancDataset dataset(orthanc, uri);
+
+    std::auto_ptr<DicomStructureSet> result(new DicomStructureSet(dataset));
+
+    std::set<std::string> instances;
+    result->GetReferencedInstances(instances);
+
+    for (std::set<std::string>::const_iterator it = instances.begin();
+         it != instances.end(); ++it)
+    {
+      Json::Value lookup;
+      MessagingToolbox::RestApiPost(lookup, orthanc, "/tools/lookup", *it);
+
+      if (lookup.type() != Json::arrayValue ||
+          lookup.size() != 1 ||
+          !lookup[0].isMember("Type") ||
+          !lookup[0].isMember("Path") ||
+          lookup[0]["Type"].type() != Json::stringValue ||
+          lookup[0]["ID"].type() != Json::stringValue ||
+          lookup[0]["Type"].asString() != "Instance")
+      {
+        throw Orthanc::OrthancException(Orthanc::ErrorCode_UnknownResource);          
+      }
+
+      OrthancPlugins::FullOrthancDataset slice
+        (orthanc, "/instances/" + lookup[0]["ID"].asString() + "/tags");
+      Orthanc::DicomMap m;
+      MessagingToolbox::ConvertDataset(m, slice);
+      result->AddReferencedSlice(m);
+    }
+
+    result->CheckReferencedSlices();
+
+    return result.release();
+  }
+
 }