--- a/Framework/Layers/DicomStructureSetRendererFactory.cpp	Thu Nov 16 15:00:45 2017 +0100
+++ b/Framework/Layers/DicomStructureSetRendererFactory.cpp	Fri Nov 17 18:01:31 2017 +0100
@@ -26,57 +26,85 @@
   class DicomStructureSetRendererFactory::Renderer : public ILayerRenderer
   {
   private:
-    DicomStructureSet&  structureSet_;
+    class Structure
+    {
+    private:
+      bool                                                         visible_;
+      uint8_t                                                      red_;
+      uint8_t                                                      green_;
+      uint8_t                                                      blue_;
+      std::string                                                  name_;
+      std::vector< std::vector<DicomStructureSet::PolygonPoint> >  polygons_;
+
+    public:
+      Structure(DicomStructureSet& structureSet,
+                const CoordinateSystem3D& slice,
+                size_t index) :
+        name_(structureSet.GetStructureName(index))
+      {
+        structureSet.GetStructureColor(red_, green_, blue_, index);
+        visible_ = structureSet.ProjectStructure(polygons_, index, slice);
+      }
+
+      void Render(CairoContext& context)
+      {
+        if (visible_)
+        {
+          cairo_t* cr = context.GetObject();
+        
+          context.SetSourceColor(red_, green_, blue_);
+
+          for (size_t i = 0; i < polygons_.size(); i++)
+          {
+            cairo_move_to(cr, polygons_[i][0].first, polygons_[i][0].second);
+
+            for (size_t j = 1; j < polygons_[i].size(); j++)
+            {
+              cairo_line_to(cr, polygons_[i][j].first, polygons_[i][j].second);
+            }
+
+            cairo_line_to(cr, polygons_[i][0].first, polygons_[i][0].second);
+            cairo_stroke(cr);
+          }
+        }
+      }
+    };
+
+    typedef std::list<Structure*>  Structures;
+    
     CoordinateSystem3D  slice_;
-    bool                visible_;
-
+    Structures          structures_;
+    
   public:
     Renderer(DicomStructureSet& structureSet,
              const CoordinateSystem3D& slice) :
-      structureSet_(structureSet),
-      slice_(slice),
-      visible_(true)
+      slice_(slice)
     {
+      for (size_t k = 0; k < structureSet.GetStructureCount(); k++)
+      {
+        structures_.push_back(new Structure(structureSet, slice, k));
+      }
+    }
+
+    virtual ~Renderer()
+    {
+      for (Structures::iterator it = structures_.begin();
+           it != structures_.end(); ++it)
+      {
+        delete *it;
+      }
     }
 
     virtual bool RenderLayer(CairoContext& context,
                              const ViewportGeometry& view)
     {
-      if (visible_)
-      {
-        cairo_set_line_width(context.GetObject(), 3.0f / view.GetZoom());
-
-        cairo_t* cr = context.GetObject();
-
-        for (size_t k = 0; k < structureSet_.GetStructureCount(); k++)
-        {
-          /*if (structureSet_.GetStructureName(k) != "CORD")
-          {
-            continue;
-            }*/
-          
-          std::vector< std::vector<DicomStructureSet::PolygonPoint> >  polygons;
+      cairo_set_line_width(context.GetObject(), 2.0f / view.GetZoom());
 
-          if (structureSet_.ProjectStructure(polygons, k, slice_))
-          {
-            uint8_t red, green, blue;
-            structureSet_.GetStructureColor(red, green, blue, k);
-            context.SetSourceColor(red, green, blue);
-
-            for (size_t i = 0; i < polygons.size(); i++)
-            {
-              cairo_move_to(cr, polygons[i][0].first, polygons[i][0].second);
-
-              for (size_t j = 1; j < polygons[i].size(); j++)
-              {
-                cairo_line_to(cr, polygons[i][j].first, polygons[i][j].second);
-              }
-
-              cairo_line_to(cr, polygons[i][0].first, polygons[i][0].second);
-              cairo_stroke(cr);
-            }
-          }
-        }
+      for (Structures::const_iterator it = structures_.begin();
+           it != structures_.end(); ++it)
+      {
+        assert(*it != NULL);
+        (*it)->Render(context);
       }
 
       return true;
@@ -89,9 +117,8 @@
 
     virtual void SetLayerStyle(const RenderStyle& style)
     {
-      visible_ = style.visible_;
     }
-
+    
     virtual bool IsFullQuality()
     {
       return true;

--- a/Framework/Toolbox/DicomStructureSet.cpp	Thu Nov 16 15:00:45 2017 +0100
+++ b/Framework/Toolbox/DicomStructureSet.cpp	Fri Nov 17 18:01:31 2017 +0100
@@ -213,8 +213,10 @@
                                            double& y2,
                                            const CoordinateSystem3D& slice) const
   {
+    // TODO Optimize this method using a sweep-line algorithm for polygons
+    
     if (!hasSlice_ ||
-        points_.empty())
+        points_.size() <= 1)
     {
       return false;
     }
@@ -226,42 +228,110 @@
     if (GeometryToolbox::IsParallelOrOpposite
         (isOpposite, slice.GetNormal(), geometry_.GetAxisY()))
     {
-      if (y >= extent_.GetY1() &&
-          y <= extent_.GetY2())
+      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++)
       {
-        Vector p1 = (geometry_.MapSliceToWorldCoordinates(extent_.GetX1(), y) +
+        // 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(extent_.GetX2(), y) -
+        Vector p2 = (geometry_.MapSliceToWorldCoordinates(xmax, y) -
                      sliceThickness_ / 2.0 * geometry_.GetNormal());
-
+          
         slice.ProjectPoint(x1, y1, p1);
         slice.ProjectPoint(x2, y2, p2);
         return true;
       }
-      else
-      {
-        return false;
-      }
     }
     else if (GeometryToolbox::IsParallelOrOpposite
              (isOpposite, slice.GetNormal(), geometry_.GetAxisX()))
     {
-      if (x >= extent_.GetX1() &&
-          x <= extent_.GetX2())
+      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++)
       {
-        Vector p1 = (geometry_.MapSliceToWorldCoordinates(x, extent_.GetY1()) +
+        // 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, extent_.GetY2()) -
+        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
-      {
-        return false;
-      }
     }
     else
     {
@@ -404,6 +474,7 @@
         }
 
         Polygon polygon(sopInstanceUid);
+        polygon.Reserve(countPoints);
 
         for (size_t k = 0; k < countPoints; k++)
         {
@@ -671,6 +742,7 @@
     else if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisX()) ||
              GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisY()))
     {
+#if 1
       // Sagittal or coronal projection
       std::vector<BoostPolygon> projected;
   
@@ -698,6 +770,22 @@
           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<PolygonPoint> 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;
     }

--- a/Framework/Toolbox/DicomStructureSet.h	Thu Nov 16 15:00:45 2017 +0100
+++ b/Framework/Toolbox/DicomStructureSet.h	Fri Nov 17 18:01:31 2017 +0100
@@ -57,7 +57,7 @@
 
     typedef std::map<std::string, ReferencedSlice>  ReferencedSlices;
     
-    typedef std::list<Vector>  Points;
+    typedef std::vector<Vector>  Points;
 
     class Polygon
     {
@@ -79,6 +79,11 @@
       {
       }
 
+      void Reserve(size_t n)
+      {
+        points_.reserve(n);
+      }
+
       void AddPoint(const Vector& v);
 
       bool UpdateReferencedSlice(const ReferencedSlices& slices);

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Resources/Computations/IntersectSegmentAndHorizontalLine.py	Fri Nov 17 18:01:31 2017 +0100
@@ -0,0 +1,21 @@
+#!/usr/bin/env python
+
+from sympy import *
+
+# Intersection between the 2D line segment (prevX,prevY)-(curX,curY) and the
+# horizontal line "y = y0" using homogeneous coordinates
+
+prevX, prevY, curX, curY, y0 = symbols('prevX prevY curX curY y0')
+
+p1 = Matrix([prevX, prevY, 1])
+p2 = Matrix([curX, curY, 1])
+l1 = p1.cross(p2)
+
+h1 = Matrix([0, y0, 1])
+h2 = Matrix([1, y0, 1])
+l2 = h1.cross(h2)
+
+a = l1.cross(l2)
+
+#pprint(cse(a/a[2], symbols = symbols('a b')))
+pprint(a / a[2])

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Resources/Computations/IntersectSegmentAndVerticalLine.py	Fri Nov 17 18:01:31 2017 +0100
@@ -0,0 +1,20 @@
+#!/usr/bin/env python
+
+from sympy import *
+
+# Intersection between the 2D line segment (prevX,prevY)-(curX,curY) and the
+# vertical line "x = x0" using homogeneous coordinates
+
+prevX, prevY, curX, curY, x0 = symbols('prevX prevY curX curY x0')
+
+p1 = Matrix([prevX, prevY, 1])
+p2 = Matrix([curX, curY, 1])
+l1 = p1.cross(p2)
+
+h1 = Matrix([x0, 0, 1])
+h2 = Matrix([x0, 1, 1])
+l2 = h1.cross(h2)
+
+a = l1.cross(l2)
+
+pprint(a / a[2])