--- a/Framework/Toolbox/DicomStructureSet2.cpp	Sat Sep 14 17:27:41 2019 +0200
+++ b/Framework/Toolbox/DicomStructureSet2.cpp	Fri Sep 20 11:58:00 2019 +0200
@@ -18,73 +18,19 @@
  * along with this program. If not, see <http://www.gnu.org/licenses/>.
  **/
 
-
 #include "DicomStructureSet2.h"
 
-#include "../Toolbox/GeometryToolbox.h"
+#include "../Toolbox/LinearAlgebra.h"
+#include "../StoneException.h"
 
 #include <Core/Logging.h>
 #include <Core/OrthancException.h>
 #include <Core/Toolbox.h>
+#include <Core/DicomFormat/DicomTag.h>
+
 #include <Plugins/Samples/Common/FullOrthancDataset.h>
 #include <Plugins/Samples/Common/DicomDatasetReader.h>
 
-#include <limits>
-#include <stdio.h>
-#include <boost/geometry.hpp>
-#include <boost/geometry/geometries/point_xy.hpp>
-#include <boost/geometry/geometries/polygon.hpp>
-#include <boost/geometry/multi/geometries/multi_polygon.hpp>
-
-typedef boost::geometry::model::d2::point_xy<double> BoostPoint;
-typedef boost::geometry::model::polygon<BoostPoint> BoostPolygon;
-typedef boost::geometry::model::multi_polygon<BoostPolygon>  BoostMultiPolygon;
-
-
-static void Union(BoostMultiPolygon& output,
-                  std::vector<BoostPolygon>& input)
-{
-  for (size_t i = 0; i < input.size(); i++)
-  {
-    boost::geometry::correct(input[i]);
-  }
-  
-  if (input.size() == 0)
-  {
-    output.clear();
-  }
-  else if (input.size() == 1)
-  {
-    output.resize(1);
-    output[0] = input[0];
-  }
-  else
-  {
-    boost::geometry::union_(input[0], input[1], output);
-
-    for (size_t i = 0; i < input.size(); i++)
-    {
-      BoostMultiPolygon tmp;
-      boost::geometry::union_(output, input[i], tmp);
-      output = tmp;
-    }
-  }
-}
-
-
-static BoostPolygon CreateRectangle(float x1, float y1,
-                                    float x2, float y2)
-{
-  BoostPolygon r;
-  boost::geometry::append(r, BoostPoint(x1, y1));
-  boost::geometry::append(r, BoostPoint(x1, y2));
-  boost::geometry::append(r, BoostPoint(x2, y2));
-  boost::geometry::append(r, BoostPoint(x2, y1));
-  return r;
-}
-
-
-
 namespace OrthancStone
 {
   static const OrthancPlugins::DicomTag DICOM_TAG_CONTOUR_GEOMETRIC_TYPE(0x3006, 0x0042);
@@ -100,8 +46,7 @@
   static const OrthancPlugins::DicomTag DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE(0x3006, 0x0080);
   static const OrthancPlugins::DicomTag DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE(0x3006, 0x0020);
 
-
-  static uint8_t ConvertColor(double v)
+  static inline uint8_t ConvertAndClipToByte(double v)
   {
     if (v < 0)
     {
@@ -117,325 +62,98 @@
     }
   }
 
-
-  static bool ParseVector(Vector& target,
-                          const OrthancPlugins::IDicomDataset& dataset,
-                          const OrthancPlugins::DicomPath& tag)
+  static bool ReadDicomToVector(Vector& target,
+    const OrthancPlugins::IDicomDataset& dataset,
+    const OrthancPlugins::DicomPath& tag)
   {
     std::string value;
     return (dataset.GetStringValue(value, tag) &&
-            LinearAlgebra::ParseVector(target, value));
+      LinearAlgebra::ParseVector(target, value));
   }
-
-  void DicomStructureSet2::Polygon::CheckPointIsOnSlice(const Vector& v) const
+  
+  std::ostream& operator<<(std::ostream& s, const OrthancPlugins::DicomPath& dicomPath)
   {
-    if (hasSlice_)
-    {
-      double magnitude =
-        GeometryToolbox::ProjectAlongNormal(v, geometry_.GetNormal());
-      if(!LinearAlgebra::IsNear(
-        magnitude,
-        projectionAlongNormal_,
-        sliceThickness_ / 2.0 /* in mm */ ))
-      {
-        LOG(ERROR) << "This RT-STRUCT contains a point that is off the "
-          << "slice of its instance | "
-          << "magnitude = " << magnitude << " | "
-          << "projectionAlongNormal_ = " << projectionAlongNormal_ << " | "
-          << "tolerance (sliceThickness_ / 2.0) = " << (sliceThickness_ / 2.0);
-
-        throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
-      }
-    }
-  }
-
-  bool DicomStructureSet2::Polygon::IsPointOnSliceIfAny(const Vector& v) const
-  {
-    if (hasSlice_)
+    std::stringstream tmp;
+    for (size_t i = 0; i < dicomPath.GetPrefixLength(); ++i)
     {
-      double magnitude = 
-        GeometryToolbox::ProjectAlongNormal(v, geometry_.GetNormal());
-      bool onSlice = LinearAlgebra::IsNear(
-        magnitude,
-        projectionAlongNormal_,
-        sliceThickness_ / 2.0 /* in mm */);
-      if (!onSlice)
-      {
-        LOG(WARNING) << "This RT-STRUCT contains a point that is off the "
-          << "slice of its instance | "
-          << "magnitude = " << magnitude << " | "
-          << "projectionAlongNormal_ = " << projectionAlongNormal_ << " | "
-          << "tolerance (sliceThickness_ / 2.0) = " << (sliceThickness_ / 2.0);
-      }
-      return onSlice;
+      OrthancPlugins::DicomTag tag = dicomPath.GetPrefixTag(i);
+      
+      // We use this other object to be able to use GetMainTagsName
+      // and Format
+      Orthanc::DicomTag tag2(tag.GetGroup(), tag.GetElement());
+      size_t index = dicomPath.GetPrefixIndex(i);
+      tmp << tag2.GetMainTagsName() << " (" << tag2.Format() << ") [" << index << "] / ";
     }
-    else
-    {
-      return false;
-    }
-  }
-
-  void DicomStructureSet2::Polygon::AddPoint(const Vector& v)
-  {
-#if 1
-    // BGO 2019-09-03
-    if (IsPointOnSliceIfAny(v))
-    {
-      points_.push_back(v);
-    }
-#else
-    CheckPoint(v);
-    points_.push_back(v);
-#endif 
+    const OrthancPlugins::DicomTag& tag = dicomPath.GetFinalTag();
+    Orthanc::DicomTag tag2(tag.GetGroup(), tag.GetElement());
+    tmp << tag2.GetMainTagsName() << " (" << tag2.Format() << ")";
+    s << tmp.str();
+    return s;
   }
 
 
-  bool DicomStructureSet2::Polygon::UpdateReferencedSlice(const ReferencedSlices& slices)
+  void DicomStructureSet2::SetContents(const OrthancPlugins::FullOrthancDataset& tags)
   {
-    if (hasSlice_)
-    {
-      return true;
-    }
-    else
+    FillStructuresFromDataset(tags);
+    ComputeDependentProperties();
+  }
+
+  void DicomStructureSet2::ComputeDependentProperties()
+  {
+    for (size_t i = 0; i < structures_.size(); ++i)
     {
-      ReferencedSlices::const_iterator it = slices.find(sopInstanceUid_);
-      
-      if (it == slices.end())
-      {
-        return false;
-      }
-      else
-      {
-        const CoordinateSystem3D& geometry = it->second.geometry_;
-        
-        hasSlice_ = true;
-        geometry_ = geometry;
-        projectionAlongNormal_ = GeometryToolbox::ProjectAlongNormal(geometry.GetOrigin(), geometry.GetNormal());
-        sliceThickness_ = it->second.thickness_;
-
-        extent_.Reset();
-        
-        for (Points::const_iterator it = points_.begin(); it != points_.end(); ++it)
-        {
-          if (IsPointOnSliceIfAny(*it))
-          {
-            double x, y;
-            geometry.ProjectPoint(x, y, *it);
-            extent_.AddPoint(x, y);
-          }
-        }
-        return true;
-      }
+      structures_[i].ComputeDependentProperties();
     }
   }
 
-  bool DicomStructureSet2::Polygon::IsOnSlice(const CoordinateSystem3D& slice) const
-  {
-    bool isOpposite = false;
-
-    if (points_.empty() ||
-        !hasSlice_ ||
-        !GeometryToolbox::IsParallelOrOpposite(isOpposite, slice.GetNormal(), geometry_.GetNormal()))
-    {
-      return false;
-    }
-    
-    double d = GeometryToolbox::ProjectAlongNormal(slice.GetOrigin(), geometry_.GetNormal());
-
-    return (LinearAlgebra::IsNear(d, projectionAlongNormal_,
-                                  sliceThickness_ / 2.0));
-  }
-    
-  bool DicomStructureSet2::Polygon::Project(double& x1,
-                                           double& y1,
-                                           double& x2,
-                                           double& y2,
-                                           const CoordinateSystem3D& slice) const
-  {
-    // TODO: optimize this method using a sweep-line algorithm for polygons
-    
-    if (!hasSlice_ ||
-        points_.size() <= 1)
-    {
-      return false;
-    }
-
-    double x, y;
-    geometry_.ProjectPoint(x, y, slice.GetOrigin());
-      
-    bool isOpposite;
-    if (GeometryToolbox::IsParallelOrOpposite
-        (isOpposite, slice.GetNormal(), geometry_.GetAxisY()))
-    {
-      if (y < extent_.GetY1() ||
-          y > extent_.GetY2())
-      {
-        // The polygon does not intersect the input slice
-        return false;
-      }
-        
-      bool isFirst = true;
-      double xmin = std::numeric_limits<double>::infinity();
-      double xmax = -std::numeric_limits<double>::infinity();
-
-      double prevX, prevY;
-      geometry_.ProjectPoint(prevX, prevY, points_[points_.size() - 1]);
-        
-      for (size_t i = 0; i < points_.size(); i++)
-      {
-        // Reference: ../../Resources/Computations/IntersectSegmentAndHorizontalLine.py
-        double curX, curY;
-        geometry_.ProjectPoint(curX, curY, points_[i]);
-
-        if ((prevY < y && curY > y) ||
-            (prevY > y && curY < y))
-        {
-          double p = (curX * prevY - curY * prevX + y * (prevX - curX)) / (prevY - curY);
-          xmin = std::min(xmin, p);
-          xmax = std::max(xmax, p);
-          isFirst = false;
-        }
-
-        prevX = curX;
-        prevY = curY;
-      }
-        
-      if (isFirst)
-      {
-        return false;
-      }
-      else
-      {
-        Vector p1 = (geometry_.MapSliceToWorldCoordinates(xmin, y) +
-                     sliceThickness_ / 2.0 * geometry_.GetNormal());
-        Vector p2 = (geometry_.MapSliceToWorldCoordinates(xmax, y) -
-                     sliceThickness_ / 2.0 * geometry_.GetNormal());
-          
-        slice.ProjectPoint(x1, y1, p1);
-        slice.ProjectPoint(x2, y2, p2);
-        return true;
-      }
-    }
-    else if (GeometryToolbox::IsParallelOrOpposite
-             (isOpposite, slice.GetNormal(), geometry_.GetAxisX()))
-    {
-      if (x < extent_.GetX1() ||
-          x > extent_.GetX2())
-      {
-        return false;
-      }
-
-      bool isFirst = true;
-      double ymin = std::numeric_limits<double>::infinity();
-      double ymax = -std::numeric_limits<double>::infinity();
-
-      double prevX, prevY;
-      geometry_.ProjectPoint(prevX, prevY, points_[points_.size() - 1]);
-        
-      for (size_t i = 0; i < points_.size(); i++)
-      {
-        // Reference: ../../Resources/Computations/IntersectSegmentAndVerticalLine.py
-        double curX, curY;
-        geometry_.ProjectPoint(curX, curY, points_[i]);
-
-        if ((prevX < x && curX > x) ||
-            (prevX > x && curX < x))
-        {
-          double p = (curX * prevY - curY * prevX + x * (curY - prevY)) / (curX - prevX);
-          ymin = std::min(ymin, p);
-          ymax = std::max(ymax, p);
-          isFirst = false;
-        }
-
-        prevX = curX;
-        prevY = curY;
-      }
-        
-      if (isFirst)
-      {
-        return false;
-      }
-      else
-      {
-        Vector p1 = (geometry_.MapSliceToWorldCoordinates(x, ymin) +
-                     sliceThickness_ / 2.0 * geometry_.GetNormal());
-        Vector p2 = (geometry_.MapSliceToWorldCoordinates(x, ymax) -
-                     sliceThickness_ / 2.0 * geometry_.GetNormal());
-
-        slice.ProjectPoint(x1, y1, p1);
-        slice.ProjectPoint(x2, y2, p2);
-
-        // TODO WHY THIS???
-        y1 = -y1;
-        y2 = -y2;
-
-        return true;
-      }
-    }
-    else
-    {
-      // Should not happen
-      return false;
-    }
-  }
-
-  
-  const DicomStructureSet2::Structure& DicomStructureSet2::GetStructure(size_t index) const
-  {
-    if (index >= structures_.size())
-    {
-      throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
-    }
-
-    return structures_[index];
-  }
-
-
-  DicomStructureSet2::Structure& DicomStructureSet2::GetStructure(size_t index)
-  {
-    if (index >= structures_.size())
-    {
-      throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
-    }
-
-    return structures_[index];
-  }
-
-  DicomStructureSet2::DicomStructureSet2(const OrthancPlugins::FullOrthancDataset& tags)
+  void DicomStructureSet2::FillStructuresFromDataset(const OrthancPlugins::FullOrthancDataset& tags)
   {
     OrthancPlugins::DicomDatasetReader reader(tags);
-    
-    size_t count, tmp;
+
+    // a few sanity checks
+    size_t count = 0, tmp = 0;
+
+    //  DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE (0x3006, 0x0080);
+    //  DICOM_TAG_ROI_CONTOUR_SEQUENCE         (0x3006, 0x0039);
+    //  DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE   (0x3006, 0x0020);
     if (!tags.GetSequenceSize(count, DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE) ||
-        !tags.GetSequenceSize(tmp, DICOM_TAG_ROI_CONTOUR_SEQUENCE) ||
-        tmp != count ||
-        !tags.GetSequenceSize(tmp, DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE) ||
-        tmp != count)
+      !tags.GetSequenceSize(tmp, DICOM_TAG_ROI_CONTOUR_SEQUENCE) ||
+      tmp != count ||
+      !tags.GetSequenceSize(tmp, DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE) ||
+      tmp != count)
     {
       throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
     }
 
+    // let's now parse the structures stored in the dicom file
+    //  DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE  (0x3006, 0x0080)
+    //  DICOM_TAG_RT_ROI_INTERPRETED_TYPE       (0x3006, 0x00a4)
+    //  DICOM_TAG_ROI_DISPLAY_COLOR             (0x3006, 0x002a)
+    //  DICOM_TAG_ROI_NAME                      (0x3006, 0x0026)
     structures_.resize(count);
     for (size_t i = 0; i < count; i++)
     {
+      // (0x3006, 0x0080)[i]/(0x3006, 0x00a4)
       structures_[i].interpretation_ = reader.GetStringValue
-        (OrthancPlugins::DicomPath(DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE, i,
-                                   DICOM_TAG_RT_ROI_INTERPRETED_TYPE),
-         "No interpretation");
+      (OrthancPlugins::DicomPath(DICOM_TAG_RT_ROI_OBSERVATIONS_SEQUENCE, i,
+        DICOM_TAG_RT_ROI_INTERPRETED_TYPE),
+        "No interpretation");
 
+      // (0x3006, 0x0020)[i]/(0x3006, 0x0026)
       structures_[i].name_ = reader.GetStringValue
-        (OrthancPlugins::DicomPath(DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE, i,
-                                   DICOM_TAG_ROI_NAME),
-         "No name");
+      (OrthancPlugins::DicomPath(DICOM_TAG_STRUCTURE_SET_ROI_SEQUENCE, i,
+        DICOM_TAG_ROI_NAME),
+        "No name");
 
       Vector color;
-      if (ParseVector(color, tags, OrthancPlugins::DicomPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i,
-                                                             DICOM_TAG_ROI_DISPLAY_COLOR)) &&
-          color.size() == 3)
+      // (0x3006, 0x0039)[i]/(0x3006, 0x002a)
+      if (ReadDicomToVector(color, tags, OrthancPlugins::DicomPath(
+        DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, DICOM_TAG_ROI_DISPLAY_COLOR)) 
+        && color.size() == 3)
       {
-        structures_[i].red_ = ConvertColor(color[0]);
-        structures_[i].green_ = ConvertColor(color[1]);
-        structures_[i].blue_ = ConvertColor(color[2]);
+        structures_[i].red_   = ConvertAndClipToByte(color[0]);
+        structures_[i].green_ = ConvertAndClipToByte(color[1]);
+        structures_[i].blue_  = ConvertAndClipToByte(color[2]);
       }
       else
       {
@@ -445,91 +163,103 @@
       }
 
       size_t countSlices;
-      if (!tags.GetSequenceSize(countSlices, OrthancPlugins::DicomPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i,
-                                                                       DICOM_TAG_CONTOUR_SEQUENCE)))
+      //  DICOM_TAG_ROI_CONTOUR_SEQUENCE          (0x3006, 0x0039);
+      //  DICOM_TAG_CONTOUR_SEQUENCE              (0x3006, 0x0040);
+      if (!tags.GetSequenceSize(countSlices, OrthancPlugins::DicomPath(
+        DICOM_TAG_ROI_CONTOUR_SEQUENCE, i, DICOM_TAG_CONTOUR_SEQUENCE)))
       {
+        LOG(WARNING) << "DicomStructureSet2::SetContents | structure \"" << structures_[i].name_ << "\" has no slices!";
         countSlices = 0;
       }
 
-      LOG(INFO) << "New RT structure: \"" << structures_[i].name_ 
-                   << "\" with interpretation \"" << structures_[i].interpretation_
-                   << "\" containing " << countSlices << " slices (color: " 
-                   << static_cast<int>(structures_[i].red_) << "," 
-                   << static_cast<int>(structures_[i].green_) << ","
-                   << static_cast<int>(structures_[i].blue_) << ")";
+      LOG(INFO) << "New RT structure: \"" << structures_[i].name_
+        << "\" with interpretation \"" << structures_[i].interpretation_
+        << "\" containing " << countSlices << " slices (color: "
+        << static_cast<int>(structures_[i].red_) << ","
+        << static_cast<int>(structures_[i].green_) << ","
+        << static_cast<int>(structures_[i].blue_) << ")";
 
       // These temporary variables avoid allocating many vectors in the loop below
-      OrthancPlugins::DicomPath countPointsPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i,
-                                                DICOM_TAG_CONTOUR_SEQUENCE, 0,
-                                                DICOM_TAG_NUMBER_OF_CONTOUR_POINTS);
+      
+      // (0x3006, 0x0039)[i]/(0x3006, 0x0040)[0]/(0x3006, 0x0046)
+      OrthancPlugins::DicomPath countPointsPath(
+        DICOM_TAG_ROI_CONTOUR_SEQUENCE, i,
+        DICOM_TAG_CONTOUR_SEQUENCE, 0,
+        DICOM_TAG_NUMBER_OF_CONTOUR_POINTS);
 
-      OrthancPlugins::DicomPath geometricTypePath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i,
-                                                  DICOM_TAG_CONTOUR_SEQUENCE, 0,
-                                                  DICOM_TAG_CONTOUR_GEOMETRIC_TYPE);
-      
-      OrthancPlugins::DicomPath imageSequencePath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i,
-                                                  DICOM_TAG_CONTOUR_SEQUENCE, 0,
-                                                  DICOM_TAG_CONTOUR_IMAGE_SEQUENCE);
+      OrthancPlugins::DicomPath geometricTypePath(
+        DICOM_TAG_ROI_CONTOUR_SEQUENCE, i,
+        DICOM_TAG_CONTOUR_SEQUENCE, 0,
+        DICOM_TAG_CONTOUR_GEOMETRIC_TYPE);
+
+      OrthancPlugins::DicomPath imageSequencePath(
+        DICOM_TAG_ROI_CONTOUR_SEQUENCE, i,
+        DICOM_TAG_CONTOUR_SEQUENCE, 0,
+        DICOM_TAG_CONTOUR_IMAGE_SEQUENCE);
 
       // (3006,0039)[i] / (0x3006, 0x0040)[0] / (0x3006, 0x0016)[0] / (0x0008, 0x1155)
-      OrthancPlugins::DicomPath referencedInstancePath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i,
-                                                       DICOM_TAG_CONTOUR_SEQUENCE, 0,
-                                                       DICOM_TAG_CONTOUR_IMAGE_SEQUENCE, 0,
-                                                       DICOM_TAG_REFERENCED_SOP_INSTANCE_UID);
+      OrthancPlugins::DicomPath referencedInstancePath(
+        DICOM_TAG_ROI_CONTOUR_SEQUENCE, i,
+        DICOM_TAG_CONTOUR_SEQUENCE, 0,
+        DICOM_TAG_CONTOUR_IMAGE_SEQUENCE, 0,
+        DICOM_TAG_REFERENCED_SOP_INSTANCE_UID);
 
-      OrthancPlugins::DicomPath contourDataPath(DICOM_TAG_ROI_CONTOUR_SEQUENCE, i,
-                                                DICOM_TAG_CONTOUR_SEQUENCE, 0,
-                                                DICOM_TAG_CONTOUR_DATA);
+      OrthancPlugins::DicomPath contourDataPath(
+        DICOM_TAG_ROI_CONTOUR_SEQUENCE, i,
+        DICOM_TAG_CONTOUR_SEQUENCE, 0,
+        DICOM_TAG_CONTOUR_DATA);
 
       for (size_t j = 0; j < countSlices; j++)
       {
-        unsigned int countPoints;
+        unsigned int countPoints = 0;
 
         countPointsPath.SetPrefixIndex(1, j);
         if (!reader.GetUnsignedIntegerValue(countPoints, countPointsPath))
         {
+          LOG(ERROR) << "Dicom path " << countPointsPath << " is not valid (should contain an unsigned integer)";
           throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
         }
-            
+
         //LOG(INFO) << "Parsing slice containing " << countPoints << " vertices";
 
         geometricTypePath.SetPrefixIndex(1, j);
         std::string type = reader.GetMandatoryStringValue(geometricTypePath);
         if (type != "CLOSED_PLANAR")
         {
+          // TODO: support points!!
           LOG(WARNING) << "Ignoring contour with geometry type: " << type;
           continue;
         }
 
-        size_t size;
+        size_t size = 0;
 
         imageSequencePath.SetPrefixIndex(1, j);
         if (!tags.GetSequenceSize(size, imageSequencePath) || size != 1)
         {
           LOG(ERROR) << "The ContourImageSequence sequence (tag 3006,0016) must be present and contain one entry.";
-          throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);          
+          throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
         }
 
         referencedInstancePath.SetPrefixIndex(1, j);
         std::string sopInstanceUid = reader.GetMandatoryStringValue(referencedInstancePath);
 
-        contourDataPath.SetPrefixIndex(1, j);        
+        contourDataPath.SetPrefixIndex(1, j);
         std::string slicesData = reader.GetMandatoryStringValue(contourDataPath);
 
         Vector points;
         if (!LinearAlgebra::ParseVector(points, slicesData) ||
-            points.size() != 3 * countPoints)
+          points.size() != 3 * countPoints)
         {
-          throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);          
+          throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
         }
 
         // seen in real world
-        if(Orthanc::Toolbox::StripSpaces(sopInstanceUid) == "") 
+        if (Orthanc::Toolbox::StripSpaces(sopInstanceUid) == "")
         {
           LOG(ERROR) << "WARNING. The following Dicom tag (Referenced SOP Instance UID) contains an empty value : // (3006,0039)[" << i << "] / (0x3006, 0x0040)[0] / (0x3006, 0x0016)[0] / (0x0008, 0x1155)";
         }
 
-        Polygon polygon(sopInstanceUid);
+        DicomStructurePolygon2 polygon(sopInstanceUid,type);
         polygon.Reserve(countPoints);
 
         for (size_t k = 0; k < countPoints; k++)
@@ -540,297 +270,15 @@
           v[2] = points[3 * k + 2];
           polygon.AddPoint(v);
         }
-
-        structures_[i].polygons_.push_back(polygon);
-      }
-    }
-  }
-
-
-  Vector DicomStructureSet2::GetStructureCenter(size_t index) const
-  {
-    const Structure& structure = GetStructure(index);
-
-    Vector center;
-    LinearAlgebra::AssignVector(center, 0, 0, 0);
-    if (structure.polygons_.empty())
-    {
-      return center;
-    }
-
-    double n = static_cast<double>(structure.polygons_.size());
-
-    for (Polygons::const_iterator polygon = structure.polygons_.begin();
-         polygon != structure.polygons_.end(); ++polygon)
-    {
-      if (!polygon->GetPoints().empty())
-      {
-        center += polygon->GetPoints().front() / n;
-      }
-    }
-
-    return center;
-  }
-
-
-  const std::string& DicomStructureSet2::GetStructureName(size_t index) const
-  {
-    return GetStructure(index).name_;
-  }
-
-
-  const std::string& DicomStructureSet2::GetStructureInterpretation(size_t index) const
-  {
-    return GetStructure(index).interpretation_;
-  }
-
-
-  Color DicomStructureSet2::GetStructureColor(size_t index) const
-  {
-    const Structure& s = GetStructure(index);
-    return Color(s.red_, s.green_, s.blue_);
-  }
-  
-    
-  void DicomStructureSet2::GetStructureColor(uint8_t& red,
-                                            uint8_t& green,
-                                            uint8_t& blue,
-                                            size_t index) const
-  {
-    const Structure& s = GetStructure(index);
-    red = s.red_;
-    green = s.green_;
-    blue = s.blue_;
-  }
-
-
-  void DicomStructureSet2::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 DicomStructureSet2::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
-      LOG(ERROR) << "DicomStructureSet2::AddReferencedSlice(): (referencedSlices_.find(sopInstanceUid) != referencedSlices_.end()). sopInstanceUid = " << sopInstanceUid;
-     
-      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);
-
-      for (Structures::iterator structure = structures_.begin();
-           structure != structures_.end(); ++structure)
-      {
-        for (Polygons::iterator polygon = structure->polygons_.begin();
-             polygon != structure->polygons_.end(); ++polygon)
-        {
-          polygon->UpdateReferencedSlice(referencedSlices_);
-        }
+        structures_[i].AddPolygon(polygon);
       }
     }
   }
 
 
-  void DicomStructureSet2::AddReferencedSlice(const Orthanc::DicomMap& dataset)
-  {
-    CoordinateSystem3D slice(dataset);
-
-    double thickness = 1;  // 1 mm by default
-
-    std::string s;
-    Vector v;
-    if (dataset.LookupStringValue(s, Orthanc::DICOM_TAG_SLICE_THICKNESS, false) &&
-        LinearAlgebra::ParseVector(v, s) &&
-        v.size() > 0)
-    {
-      thickness = v[0];
-    }
-
-    std::string instance, series;
-    if (dataset.LookupStringValue(instance, Orthanc::DICOM_TAG_SOP_INSTANCE_UID, false) &&
-        dataset.LookupStringValue(series, Orthanc::DICOM_TAG_SERIES_INSTANCE_UID, false))
-    {
-      AddReferencedSlice(instance, series, slice, thickness);
-    }
-    else
-    {
-      throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
-    }
-  }
-
-
-  void DicomStructureSet2::CheckReferencedSlices()
+  void DicomStructureSet2::Clear()
   {
-    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_))
-        {
-          std::string sopInstanceUid = polygon->GetSopInstanceUid();
-          if (Orthanc::Toolbox::StripSpaces(sopInstanceUid) == "")
-          {
-            LOG(ERROR) << "DicomStructureSet2::CheckReferencedSlices(): "
-              << " missing information about referenced instance "
-              << "(sopInstanceUid is empty!)";
-          }
-          else
-          {
-            LOG(ERROR) << "DicomStructureSet2::CheckReferencedSlices(): "
-              << " missing information about referenced instance "
-              << "(sopInstanceUid = " << sopInstanceUid << ")";
-          }
-          //throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
-        }
-      }
-    }
-  }
-
-
-  Vector DicomStructureSet2::GetNormal() const
-  {
-    if (referencedSlices_.empty())
-    {
-      Vector v;
-      LinearAlgebra::AssignVector(v, 0, 0, 1);
-      return v;
-    }
-    else
-    {
-      return referencedSlices_.begin()->second.geometry_.GetNormal();
-    }
+    structures_.clear();
   }
 
-  
-  bool DicomStructureSet2::ProjectStructure(std::vector< std::vector<PolygonPoint2D> >& polygons,
-                                           const Structure& structure,
-                                           const CoordinateSystem3D& slice) const
-  {
-    polygons.clear();
-
-    Vector normal = GetNormal();
-    
-    bool isOpposite;    
-    if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetNormal()))
-    {
-      // This is an axial projection
-
-      for (Polygons::const_iterator polygon = structure.polygons_.begin();
-           polygon != structure.polygons_.end(); ++polygon)
-      {
-        if (polygon->IsOnSlice(slice))
-        {
-          polygons.push_back(std::vector<PolygonPoint2D>());
-          
-          for (Points::const_iterator p = polygon->GetPoints().begin();
-               p != polygon->GetPoints().end(); ++p)
-          {
-            double x, y;
-            slice.ProjectPoint(x, y, *p);
-            polygons.back().push_back(std::make_pair(x, y));
-          }
-        }
-      }
-
-      return true;
-    }
-    else if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisX()) ||
-             GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisY()))
-    {
-#if 1
-      // Sagittal or coronal projection
-      std::vector<BoostPolygon> projected;
-  
-      for (Polygons::const_iterator polygon = structure.polygons_.begin();
-           polygon != structure.polygons_.end(); ++polygon)
-      {
-        double x1, y1, x2, y2;
-        if (polygon->Project(x1, y1, x2, y2, slice))
-        {
-          projected.push_back(CreateRectangle(
-            static_cast<float>(x1), 
-            static_cast<float>(y1), 
-            static_cast<float>(x2), 
-            static_cast<float>(y2)));
-        }
-      }
-
-      BoostMultiPolygon merged;
-      Union(merged, projected);
-
-      polygons.resize(merged.size());
-      for (size_t i = 0; i < merged.size(); i++)
-      {
-        const std::vector<BoostPoint>& outer = merged[i].outer();
-
-        polygons[i].resize(outer.size());
-        for (size_t j = 0; j < outer.size(); j++)
-        {
-          polygons[i][j] = std::make_pair(outer[j].x(), outer[j].y());
-        }
-      }  
-#else
-      for (Polygons::iterator polygon = structure.polygons_.begin();
-           polygon != structure.polygons_.end(); ++polygon)
-      {
-        double x1, y1, x2, y2;
-        if (polygon->Project(x1, y1, x2, y2, slice))
-        {
-          std::vector<PolygonPoint2D> p(4);
-          p[0] = std::make_pair(x1, y1);
-          p[1] = std::make_pair(x2, y1);
-          p[2] = std::make_pair(x2, y2);
-          p[3] = std::make_pair(x1, y2);
-          polygons.push_back(p);
-        }
-      }
-#endif
-      
-      return true;
-    }
-    else
-    {
-      return false;
-    }
-  }
 }