# HG changeset patch
# User Sebastien Jodogne <s.jodogne@gmail.com>
# Date 1510158458 -3600
# Node ID c9e88e7935a47237b21dd0aebce3d41cb85f83a9
# Parent  44fc253d487663c2026c53a1d3628f026a5482e9
rt-struct projection

diff -r 44fc253d4876 -r c9e88e7935a4 Applications/Samples/SingleVolumeApplication.h
--- a/Applications/Samples/SingleVolumeApplication.h	Wed Nov 08 16:13:55 2017 +0100
+++ b/Applications/Samples/SingleVolumeApplication.h	Wed Nov 08 17:27:38 2017 +0100
@@ -203,7 +203,7 @@
         //ct->ScheduleLoadSeries("dd069910-4f090474-7d2bba07-e5c10783-f9e4fb1d");
         ct->ScheduleLoadSeries("a04ecf01-79b2fc33-58239f7e-ad9db983-28e81afa");  // IBA
         //ct->ScheduleLoadSeries("03677739-1d8bca40-db1daf59-d74ff548-7f6fc9c0");  // 0522c0001 TCIA
-  
+        
         std::auto_ptr<OrthancVolumeImage> pet(new OrthancVolumeImage(context.GetWebService(), true));
         //pet->ScheduleLoadSeries("aabad2e7-80702b5d-e599d26c-4f13398e-38d58a9e");
         pet->ScheduleLoadInstance("830a69ff-8e4b5ee3-b7f966c8-bccc20fb-d322dceb"); // IBA 1
@@ -220,6 +220,7 @@
         
         context.AddInteractor(new Interactor(*pet, *widget, projection, 1));
         //context.AddInteractor(new VolumeImageInteractor(*ct, *widget, projection));
+
         context.AddVolume(ct.release());
         context.AddVolume(pet.release());
 
diff -r 44fc253d4876 -r c9e88e7935a4 Framework/Layers/DicomStructureSetRendererFactory.cpp
--- a/Framework/Layers/DicomStructureSetRendererFactory.cpp	Wed Nov 08 16:13:55 2017 +0100
+++ b/Framework/Layers/DicomStructureSetRendererFactory.cpp	Wed Nov 08 17:27:38 2017 +0100
@@ -54,6 +54,11 @@
 
         for (size_t k = 0; k < structureSet_.GetStructureCount(); k++)
         {
+          /*if (structureSet_.GetStructureName(k) != "CORD")
+          {
+            continue;
+            }*/
+          
           std::vector< std::vector<DicomStructureSet::PolygonPoint> >  polygons;
 
           if (structureSet_.ProjectStructure(polygons, k, slice_))
diff -r 44fc253d4876 -r c9e88e7935a4 Framework/Toolbox/DicomStructureSet.cpp
--- a/Framework/Toolbox/DicomStructureSet.cpp	Wed Nov 08 16:13:55 2017 +0100
+++ b/Framework/Toolbox/DicomStructureSet.cpp	Wed Nov 08 17:27:38 2017 +0100
@@ -28,8 +28,61 @@
 #include <Plugins/Samples/Common/FullOrthancDataset.h>
 #include <Plugins/Samples/Common/DicomDatasetReader.h>
 
+#include <limits>
 #include <stdio.h>
 #include <boost/lexical_cast.hpp>
+#include <boost/geometry.hpp>
+#include <boost/geometry/geometries/point_xy.hpp>
+#include <boost/geometry/geometries/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
 {
@@ -449,7 +502,7 @@
 
   Vector DicomStructureSet::GetNormal() const
   {
-    if (!referencedSlices_.empty())
+    if (referencedSlices_.empty())
     {
       Vector v;
       GeometryToolbox::AssignVector(v, 0, 0, 1);
@@ -540,9 +593,81 @@
     else if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisX()) ||
              GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisY()))
     {
-      // Sagittal or coronal projections
+      // Sagittal or coronal projection
+
+      std::vector<BoostPolygon> projected;
+  
+      for (Polygons::iterator polygon = structure.polygons_.begin();
+           polygon != structure.polygons_.end(); ++polygon)
+      {
+        polygon->UpdateReferencedSlice(referencedSlices_);
+
+        // First, check that the polygon intersects the input slice
+        double zmin = std::numeric_limits<double>::infinity();
+        double zmax = -std::numeric_limits<double>::infinity();
+
+        double zref = slice.ProjectAlongNormal(slice.GetOrigin());
+        
+        for (Points::const_iterator p = polygon->GetPoints().begin();
+             p != polygon->GetPoints().end(); ++p)
+        {
+          double z = slice.ProjectAlongNormal(*p);
+          zmin = std::min(zmin, z);
+          zmax = std::max(zmax, z);
+        }
+
+        if (zmin < zref && zref < zmax)
+        {
+          // The polygon intersect the input slice, project the polygon
+          double xmin = std::numeric_limits<double>::infinity();
+          double xmax = -std::numeric_limits<double>::infinity();
+          double ymin = std::numeric_limits<double>::infinity();
+          double ymax = -std::numeric_limits<double>::infinity();
+
+          for (Points::const_iterator p = polygon->GetPoints().begin();
+               p != polygon->GetPoints().end(); ++p)
+          {
+            double x, y;
+            slice.ProjectPoint(x, y, *p);
+            xmin = std::min(xmin, x);
+            xmax = std::max(xmax, x);
+            ymin = std::min(ymin, y);
+            ymax = std::max(ymax, y);
+          }
+
+          if (GeometryToolbox::IsNear(xmin, xmax))
+          {
+            projected.push_back(CreateRectangle(xmin - polygon->GetSliceThickness() / 2.0, ymin,
+                                                xmax + polygon->GetSliceThickness() / 2.0, ymax));
+          }
+          else if (GeometryToolbox::IsNear(ymin, ymax))
+          {
+            projected.push_back(CreateRectangle(xmin, ymin - polygon->GetSliceThickness() / 2.0,
+                                                xmax, ymax + polygon->GetSliceThickness() / 2.0));
+          }
+          else
+          {
+            // Should not happen
+          }
+        }
+      }
+
+      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());
+        }
+      }  
       
-      return false;
+      return true;
     }
     else
     {
diff -r 44fc253d4876 -r c9e88e7935a4 Framework/Toolbox/DicomStructureSet.h
--- a/Framework/Toolbox/DicomStructureSet.h	Wed Nov 08 16:13:55 2017 +0100
+++ b/Framework/Toolbox/DicomStructureSet.h	Wed Nov 08 17:27:38 2017 +0100
@@ -93,6 +93,11 @@
       {
         return points_;
       }
+
+      double GetSliceThickness() const
+      {
+        return sliceThickness_;
+      }
     };
 
     typedef std::list<Polygon>  Polygons;