# HG changeset patch
# User Sebastien Jodogne <s.jodogne@gmail.com>
# Date 1641927580 -3600
# Node ID b896f20d24cab0737acef138d412f778a5aaffe4
# Parent  08f2476e8f5e8714d5f7983f900c8176eb710283
added UnionOfRectangles algorithm

diff -r 08f2476e8f5e -r b896f20d24ca OrthancStone/Resources/CMake/OrthancStoneConfiguration.cmake
--- a/OrthancStone/Resources/CMake/OrthancStoneConfiguration.cmake	Tue Jan 11 18:58:37 2022 +0100
+++ b/OrthancStone/Resources/CMake/OrthancStoneConfiguration.cmake	Tue Jan 11 19:59:40 2022 +0100
@@ -445,6 +445,8 @@
   ${ORTHANC_STONE_ROOT}/Toolbox/TextRenderer.h
   ${ORTHANC_STONE_ROOT}/Toolbox/UndoRedoStack.cpp
   ${ORTHANC_STONE_ROOT}/Toolbox/UndoRedoStack.h
+  ${ORTHANC_STONE_ROOT}/Toolbox/UnionOfRectangles.cpp
+  ${ORTHANC_STONE_ROOT}/Toolbox/UnionOfRectangles.h
   
   ${ORTHANC_STONE_ROOT}/Viewport/IViewport.h
   ${ORTHANC_STONE_ROOT}/Viewport/DefaultViewportInteractor.cpp
diff -r 08f2476e8f5e -r b896f20d24ca OrthancStone/Sources/Scene2D/ScenePoint2D.h
--- a/OrthancStone/Sources/Scene2D/ScenePoint2D.h	Tue Jan 11 18:58:37 2022 +0100
+++ b/OrthancStone/Sources/Scene2D/ScenePoint2D.h	Tue Jan 11 19:59:40 2022 +0100
@@ -88,5 +88,10 @@
     static double SquaredDistancePtSegment(const ScenePoint2D& a,
                                            const ScenePoint2D& b,
                                            const ScenePoint2D& p);
+
+    bool IsEqual(const ScenePoint2D& other) const
+    {
+      return LinearAlgebra::IsCloseToZero(DistancePtPt(*this, other));
+    }
   };
 }
diff -r 08f2476e8f5e -r b896f20d24ca OrthancStone/Sources/Toolbox/UnionOfRectangles.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/OrthancStone/Sources/Toolbox/UnionOfRectangles.cpp	Tue Jan 11 19:59:40 2022 +0100
@@ -0,0 +1,552 @@
+/**
+ * Stone of Orthanc
+ * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics
+ * Department, University Hospital of Liege, Belgium
+ * Copyright (C) 2017-2022 Osimis S.A., Belgium
+ * Copyright (C) 2021-2022 Sebastien Jodogne, ICTEAM UCLouvain, Belgium
+ *
+ * This program is free software: you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation, either version 3 of
+ * the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this program. If not, see
+ * <http://www.gnu.org/licenses/>.
+ **/
+
+
+#include "UnionOfRectangles.h"
+
+#include "Internals/OrientedIntegerLine2D.h"
+#include "Internals/RectanglesIntegerProjection.h"
+#include "SegmentTree.h"
+
+#include <OrthancException.h>
+
+#include <stack>
+
+
+namespace OrthancStone
+{
+  class UnionOfRectangles::Payload : public Orthanc::IDynamicObject
+  {    
+  private:
+    int     counter_;
+    Status  status_;
+  
+  public:
+    Payload() :
+      counter_(0),
+      status_(Status_Empty)
+    {
+    }
+
+    int GetCounter() const
+    {
+      return counter_;
+    }
+
+    Status GetStatus() const
+    {
+      return status_;
+    }
+
+    void SetStatus(Status status)
+    {
+      status_ = status;
+    }
+
+    void Increment()
+    {
+      counter_ ++;
+    }
+
+    void Decrement()
+    {
+      if (counter_ == 0)
+      {
+        throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
+      }
+      else
+      {
+        counter_ --;
+      }
+    }
+  };
+  
+  
+  class UnionOfRectangles::Factory : public SegmentTree::IPayloadFactory
+  {
+  public:
+    virtual Orthanc::IDynamicObject* Create()
+    {
+      return new Payload;
+    }
+  };
+
+
+  class UnionOfRectangles::Visitor : public SegmentTree::IVisitor
+  {
+  private:
+    Operation  operation_;
+  
+  public:
+    Visitor(Operation operation) :
+      operation_(operation)
+    {
+    }
+  
+    virtual void Visit(const SegmentTree& node,
+                       bool fullyInside)
+    {
+      Payload& payload = node.GetTypedPayload<Payload>();
+
+      if (fullyInside)
+      {
+        switch (operation_)
+        {
+          case Operation_Insert:
+            payload.Increment();
+            break;
+            
+          case Operation_Delete:
+            payload.Decrement();
+            break;
+
+          default:
+            throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
+        }
+      }
+
+      if (payload.GetCounter() > 0)
+      {
+        payload.SetStatus(Status_Full);
+      }
+      else if (node.IsLeaf())
+      {
+        payload.SetStatus(Status_Empty);
+      }
+      else if (node.GetLeftChild().GetTypedPayload<Payload>().GetStatus() == Status_Empty &&
+               node.GetRightChild().GetTypedPayload<Payload>().GetStatus() == Status_Empty)
+      {
+        payload.SetStatus(Status_Empty);
+      }
+      else
+      {
+        payload.SetStatus(Status_Partial);
+      }
+    }
+
+
+    // This is the "CONTR()" function from the textbook
+    static void IntersectComplement(std::stack<size_t>& stack,
+                                    size_t low,
+                                    size_t high,
+                                    const SegmentTree& node)
+    {
+      if (low >= high)
+      {
+        throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
+      }
+    
+      Status status = node.GetTypedPayload<Payload>().GetStatus();
+    
+      if (status != Status_Full)
+      {
+        assert(status == Status_Partial ||
+               status == Status_Empty);
+
+        // Aliases to use the same variable names as in the textbook
+        const size_t& b = low;
+        const size_t& e = high;
+        const size_t& bv = node.GetLowBound();
+        const size_t& ev = node.GetHighBound();
+      
+        if (b <= bv &&
+            ev <= e &&
+            status == Status_Empty)
+        {
+          // [B[v], E[v]] is contributed
+          if (!stack.empty() &&
+              stack.top() == bv)
+          {
+            stack.pop();     // Merge continuous segments
+          }
+          else
+          {
+            stack.push(bv);  // Beginning of edge
+          }
+
+          stack.push(ev);    // Current termination of edge
+        }
+        else
+        {
+          size_t middle = (bv + ev) / 2;
+
+          if (b < middle)
+          {
+            IntersectComplement(stack, b, e, node.GetLeftChild());
+          }
+
+          if (middle < e)
+          {
+            IntersectComplement(stack, b, e, node.GetRightChild());
+          }
+        }
+      }
+    }
+  };
+
+
+  static void AddVerticalEdges(std::list<Internals::OrientedIntegerLine2D>& edges,
+                               std::stack<size_t>& stack,
+                               size_t x,
+                               bool isLeft)
+  {
+    if (stack.size() % 2 != 0)
+    {
+      throw std::runtime_error("internal error");
+    }
+
+    typedef std::pair<size_t, size_t>  Interval;
+
+    std::list<Interval> intervals;
+
+    while (!stack.empty())
+    {
+      size_t high = stack.top();
+      stack.pop();
+      size_t low = stack.top();
+      stack.pop();
+    
+      if (!intervals.empty() &&
+          intervals.back().second == low)
+      {
+        // Extend the previous interval
+        intervals.back().second = high;
+      }
+      else
+      {
+        intervals.push_back(std::make_pair(low, high));
+      }
+    }
+  
+    for (std::list<Interval>::const_iterator
+           it = intervals.begin(); it != intervals.end(); ++it)
+    {
+      if (it->first == it->second)
+      {
+        throw std::runtime_error("parameter out of range");
+      }
+    
+      // By convention, the left sides go downward, and the right go upward
+      if (isLeft)
+      {
+        edges.push_back(Internals::OrientedIntegerLine2D(x, it->first, x, it->second));
+      }
+      else
+      {
+        edges.push_back(Internals::OrientedIntegerLine2D(x, it->second, x, it->first));
+      }
+    }
+  }
+  
+
+  class UnionOfRectangles::VerticalSide
+  {
+  private:
+    size_t  x_;
+    bool    isLeft_;
+    size_t  y1_;
+    size_t  y2_;
+
+  public:
+    VerticalSide(size_t x,
+                 bool isLeft,
+                 size_t y1,
+                 size_t y2) :
+      x_(x),
+      isLeft_(isLeft),
+      y1_(y1),
+      y2_(y2)
+    {
+      assert(y1 < y2);
+    }
+
+    size_t GetX() const
+    {
+      return x_;
+    }
+
+    bool IsLeft() const
+    {
+      return isLeft_;
+    }
+
+    size_t GetY1() const
+    {
+      return y1_;
+    }
+
+    size_t GetY2() const
+    {
+      return y2_;
+    }
+
+    bool operator< (const VerticalSide& other) const
+    {
+      if (x_ < other.x_)
+      {
+        return true;
+      }
+      else if (x_ == other.x_)
+      {
+        return isLeft_;
+      }
+      else
+      {
+        return false;
+      }
+    }
+
+    bool Equals(const VerticalSide& other) const
+    {
+      return (x_ == other.x_ &&
+              isLeft_ == other.isLeft_);
+    }
+  };
+
+
+  class UnionOfRectangles::HorizontalJunction
+  {
+  private:
+    size_t  x_;
+    size_t  y_;
+    size_t  ybis_;
+    bool    downward_;
+
+  public:
+    HorizontalJunction(size_t x,
+                       size_t y,
+                       size_t ybis,
+                       bool downward) :
+      x_(x),
+      y_(y),
+      ybis_(ybis),
+      downward_(downward)
+    {
+    }
+
+    size_t GetX() const
+    {
+      return x_;
+    }
+
+    size_t GetY() const
+    {
+      return y_;
+    }
+
+    size_t GetYBis() const
+    {
+      return ybis_;
+    }
+
+    bool IsDownward() const
+    {
+      return downward_;
+    }
+
+    bool operator< (const HorizontalJunction& other) const
+    {
+      if (y_ > other.y_)
+      {
+        return true;
+      }
+      else if (y_ == other.y_)
+      {
+        return x_ < other.x_;
+      }
+      else
+      {
+        return false;
+      }
+    }
+  };
+  
+  
+  void UnionOfRectangles::Apply(std::list< std::vector<ScenePoint2D> >& contours,
+                                const std::list<Extent2D>& rectangles)
+  {
+    /**
+     * STEP 1
+     **/
+    Internals::RectanglesIntegerProjection horizontalProjection(rectangles, true);
+    Internals::RectanglesIntegerProjection verticalProjection(rectangles, false);
+
+    assert(horizontalProjection.GetProjectedRectanglesCount() == verticalProjection.GetProjectedRectanglesCount());
+
+    
+    /**
+     * STEP 2
+     **/
+    Factory factory;
+    SegmentTree tree(0, verticalProjection.GetEndpointsCount() - 1, factory);
+
+    
+    /**
+     * STEP 3
+     **/
+    std::vector<VerticalSide> verticalSides;
+
+    const size_t countRectangles = horizontalProjection.GetProjectedRectanglesCount();
+    
+    verticalSides.reserve(2 * countRectangles);
+    
+    for (size_t i = 0; i < countRectangles; i++)
+    {
+      size_t horizontalLow = horizontalProjection.GetProjectedRectangleLow(i);
+      size_t horizontalHigh = horizontalProjection.GetProjectedRectangleHigh(i);
+      size_t verticalLow = verticalProjection.GetProjectedRectangleLow(i);
+      size_t verticalHigh = verticalProjection.GetProjectedRectangleHigh(i);
+      
+      verticalSides.push_back(VerticalSide(horizontalLow, true, verticalLow, verticalHigh));
+      verticalSides.push_back(VerticalSide(horizontalHigh, false, verticalLow, verticalHigh));
+    }
+    
+    assert(verticalSides.size() == 2 * countRectangles);
+    
+    std::sort(verticalSides.begin(), verticalSides.end());
+
+    
+    /**
+     * STEP 4
+     **/
+    std::list<Internals::OrientedIntegerLine2D> verticalEdges;
+    std::stack<size_t> stack;
+    
+    for (size_t i = 0; i < verticalSides.size(); i++)
+    {
+      if (i > 0 &&
+          !verticalSides[i].Equals(verticalSides[i - 1]))
+      {
+        // Output the stack
+        AddVerticalEdges(verticalEdges, stack, verticalSides[i - 1].GetX(),
+                         verticalSides[i - 1].IsLeft());
+      }
+
+      size_t y1 = verticalSides[i].GetY1();
+      size_t y2 = verticalSides[i].GetY2();
+      
+      if (verticalSides[i].IsLeft())
+      {
+        Visitor::IntersectComplement(stack, y1, y2, tree);
+
+        Visitor visitor(Operation_Insert);
+        tree.VisitSegment(y1, y2, visitor);
+      }
+      else
+      {
+        Visitor visitor(Operation_Delete);
+        tree.VisitSegment(y1, y2, visitor);
+
+        Visitor::IntersectComplement(stack, y1, y2, tree);
+      }
+    }
+
+    if (!verticalSides.empty() &&
+        !stack.empty())
+    {
+      AddVerticalEdges(verticalEdges, stack, verticalSides.back().GetX(),
+                       verticalSides.back().IsLeft());
+    }
+
+
+    /**
+     * STEP 5: Horizontal edges
+     **/
+    std::vector<HorizontalJunction> horizontalJunctions;
+    horizontalJunctions.reserve(2 * verticalEdges.size());
+
+    for (std::list<Internals::OrientedIntegerLine2D>::const_iterator
+           it = verticalEdges.begin(); it != verticalEdges.end(); it++)
+    {
+      assert(it->GetX1() == it->GetX2());
+      horizontalJunctions.push_back(HorizontalJunction(it->GetX1(), it->GetY1(), it->GetY2(), it->IsDownward()));
+      horizontalJunctions.push_back(HorizontalJunction(it->GetX1(), it->GetY2(), it->GetY1(), it->IsDownward()));
+    }  
+
+    assert(horizontalJunctions.size() == 2 * verticalEdges.size());
+    std::sort(horizontalJunctions.begin(), horizontalJunctions.end());
+
+    std::list<Internals::OrientedIntegerLine2D> horizontalEdges;
+    for (size_t i = 0; i < horizontalJunctions.size(); i += 2)
+    {
+      size_t x1 = horizontalJunctions[i].GetX();
+      size_t x2 = horizontalJunctions[i + 1].GetX();
+      size_t y = horizontalJunctions[i].GetY();
+
+      if ((horizontalJunctions[i].IsDownward() && y > horizontalJunctions[i].GetYBis()) ||
+          (!horizontalJunctions[i].IsDownward() && y < horizontalJunctions[i].GetYBis()))
+      {
+        horizontalEdges.push_back(Internals::OrientedIntegerLine2D(x1, y, x2, y));
+      }
+      else
+      {
+        horizontalEdges.push_back(Internals::OrientedIntegerLine2D(x2, y, x1, y));
+      }
+    }
+
+
+    /**
+     * POST-PROCESSING: Combine the separate sets of horizontal and
+     * vertical edges, into a set of 2D chains
+     **/
+    std::vector<Internals::OrientedIntegerLine2D> allEdges;
+    allEdges.reserve(horizontalEdges.size() + verticalEdges.size());
+    
+    for (std::list<Internals::OrientedIntegerLine2D>::const_iterator
+           it = horizontalEdges.begin(); it != horizontalEdges.end(); it++)
+    {
+      allEdges.push_back(*it);
+    }
+    
+    for (std::list<Internals::OrientedIntegerLine2D>::const_iterator
+           it = verticalEdges.begin(); it != verticalEdges.end(); it++)
+    {
+      allEdges.push_back(*it);
+    }
+    
+    assert(allEdges.size() == horizontalEdges.size() + verticalEdges.size());
+
+    std::list<Internals::OrientedIntegerLine2D::Chain> chains;
+    Internals::OrientedIntegerLine2D::ExtractChains(chains, allEdges);
+
+    contours.clear();
+
+    for (std::list<Internals::OrientedIntegerLine2D::Chain>::const_iterator
+           it = chains.begin(); it != chains.end(); ++it)
+    {
+      assert(!it->empty());
+      
+      std::vector<ScenePoint2D> chain;
+      chain.reserve(it->size());
+      
+      for (Internals::OrientedIntegerLine2D::Chain::const_iterator
+             it2 = it->begin(); it2 != it->end(); ++it2)
+      {
+        chain.push_back(ScenePoint2D(horizontalProjection.GetEndpointCoordinate(it2->first),
+                                     verticalProjection.GetEndpointCoordinate(it2->second)));
+      }
+
+      assert(!chain.empty());
+      contours.push_back(chain);
+    }
+  }
+}
diff -r 08f2476e8f5e -r b896f20d24ca OrthancStone/Sources/Toolbox/UnionOfRectangles.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/OrthancStone/Sources/Toolbox/UnionOfRectangles.h	Tue Jan 11 19:59:40 2022 +0100
@@ -0,0 +1,69 @@
+/**
+ * Stone of Orthanc
+ * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics
+ * Department, University Hospital of Liege, Belgium
+ * Copyright (C) 2017-2022 Osimis S.A., Belgium
+ * Copyright (C) 2021-2022 Sebastien Jodogne, ICTEAM UCLouvain, Belgium
+ *
+ * This program is free software: you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation, either version 3 of
+ * the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this program. If not, see
+ * <http://www.gnu.org/licenses/>.
+ **/
+
+
+#pragma once
+
+#include "Extent2D.h"
+#include "../Scene2D/ScenePoint2D.h"
+
+#include <list>
+
+
+namespace OrthancStone
+{
+  /**
+   * This implementation closely follows "Finding the Contour of a
+   * Union of Iso-Oriented Rectangles" by Lipski and Preparata (1980),
+   * as well as Section 8.5 (pages 340-348) of "Computational Geometry
+   * - An Introduction" by Preparata and Ian Shamos (1985).
+   **/
+  class UnionOfRectangles : public boost::noncopyable
+  {
+  private:
+    enum Operation
+    {
+      Operation_Insert,
+      Operation_Delete
+    };
+
+    enum Status
+    {
+      Status_Full,
+      Status_Partial,
+      Status_Empty
+    };
+
+    class Payload;
+    class Factory;
+    class Visitor;
+
+    class VerticalSide;
+    class HorizontalJunction;
+    
+    UnionOfRectangles();  // Pure static class
+
+  public:
+    static void Apply(std::list< std::vector<ScenePoint2D> >& contours,
+                      const std::list<Extent2D>& rectangles);
+  };
+}
diff -r 08f2476e8f5e -r b896f20d24ca UnitTestsSources/ComputationalGeometryTests.cpp
--- a/UnitTestsSources/ComputationalGeometryTests.cpp	Tue Jan 11 18:58:37 2022 +0100
+++ b/UnitTestsSources/ComputationalGeometryTests.cpp	Tue Jan 11 19:59:40 2022 +0100
@@ -25,6 +25,7 @@
 #include "../OrthancStone/Sources/Toolbox/Internals/OrientedIntegerLine2D.h"
 #include "../OrthancStone/Sources/Toolbox/Internals/RectanglesIntegerProjection.h"
 #include "../OrthancStone/Sources/Toolbox/SegmentTree.h"
+#include "../OrthancStone/Sources/Toolbox/UnionOfRectangles.h"
 
 #include <Logging.h>
 #include <OrthancException.h>
@@ -479,3 +480,224 @@
   ASSERT_EQ(5u, v[4]);
   ASSERT_EQ(5u, v[5]);
 }
+
+
+TEST(UnionOfRectangles, Textbook)
+{
+  // This is Figure 8.12 from textbook
+
+  std::list<OrthancStone::Extent2D>  rectangles;
+  rectangles.push_back(OrthancStone::Extent2D(1, 3, 13, 5));
+  rectangles.push_back(OrthancStone::Extent2D(3, 1, 7, 12));
+  rectangles.push_back(OrthancStone::Extent2D(5, 7, 11, 10));
+  rectangles.push_back(OrthancStone::Extent2D(10, 2, 14, 8));
+  rectangles.push_back(OrthancStone::Extent2D(3, 3, 4, 3));  // empty rectangle
+
+  for (unsigned int fillHole = 0; fillHole < 2; fillHole++)
+  {
+    if (fillHole)
+    {
+      rectangles.push_back(OrthancStone::Extent2D(6.5, 4.5, 10.5, 7.5));
+    }
+    
+    std::list< std::vector<OrthancStone::ScenePoint2D> > contours;
+    OrthancStone::UnionOfRectangles::Apply(contours, rectangles);
+
+    ASSERT_EQ(fillHole ? 1u : 2u, contours.size());
+    ASSERT_EQ(17u, contours.front().size());
+    ASSERT_TRUE(contours.front()[0].IsEqual(OrthancStone::ScenePoint2D(3, 12)));
+    ASSERT_TRUE(contours.front()[1].IsEqual(OrthancStone::ScenePoint2D(7, 12)));
+    ASSERT_TRUE(contours.front()[2].IsEqual(OrthancStone::ScenePoint2D(7, 10)));
+    ASSERT_TRUE(contours.front()[3].IsEqual(OrthancStone::ScenePoint2D(11, 10)));
+    ASSERT_TRUE(contours.front()[4].IsEqual(OrthancStone::ScenePoint2D(11, 8)));
+    ASSERT_TRUE(contours.front()[5].IsEqual(OrthancStone::ScenePoint2D(14, 8)));
+    ASSERT_TRUE(contours.front()[6].IsEqual(OrthancStone::ScenePoint2D(14, 2)));
+    ASSERT_TRUE(contours.front()[7].IsEqual(OrthancStone::ScenePoint2D(10, 2)));
+    ASSERT_TRUE(contours.front()[8].IsEqual(OrthancStone::ScenePoint2D(10, 3)));
+    ASSERT_TRUE(contours.front()[9].IsEqual(OrthancStone::ScenePoint2D(7, 3)));
+    ASSERT_TRUE(contours.front()[10].IsEqual(OrthancStone::ScenePoint2D(7, 1)));
+    ASSERT_TRUE(contours.front()[11].IsEqual(OrthancStone::ScenePoint2D(3, 1)));
+    ASSERT_TRUE(contours.front()[12].IsEqual(OrthancStone::ScenePoint2D(3, 3)));
+    ASSERT_TRUE(contours.front()[13].IsEqual(OrthancStone::ScenePoint2D(1, 3)));
+    ASSERT_TRUE(contours.front()[14].IsEqual(OrthancStone::ScenePoint2D(1, 5)));
+    ASSERT_TRUE(contours.front()[15].IsEqual(OrthancStone::ScenePoint2D(3, 5)));
+    ASSERT_TRUE(contours.front()[16].IsEqual(OrthancStone::ScenePoint2D(3, 12)));
+
+    if (!fillHole)
+    {
+      ASSERT_EQ(5u, contours.back().size());
+      ASSERT_TRUE(contours.back()[0].IsEqual(OrthancStone::ScenePoint2D(10, 7)));
+      ASSERT_TRUE(contours.back()[1].IsEqual(OrthancStone::ScenePoint2D(7, 7)));
+      ASSERT_TRUE(contours.back()[2].IsEqual(OrthancStone::ScenePoint2D(7, 5)));
+      ASSERT_TRUE(contours.back()[3].IsEqual(OrthancStone::ScenePoint2D(10, 5)));
+      ASSERT_TRUE(contours.back()[4].IsEqual(OrthancStone::ScenePoint2D(10, 7)));
+    }
+  }
+}
+
+
+#if 0
+static void SaveSvg(const std::list< std::vector<OrthancStone::ScenePoint2D> >& contours)
+{
+  float ww = 15;
+  float hh = 13;
+
+  FILE* fp = fopen("test.svg", "w");
+  fprintf(fp, "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n");
+  fprintf(fp, "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\" \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n");
+  fprintf(fp, "<svg width=\"%f\" height=\"%f\" viewBox=\"0 0 %f %f\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n", 100.0f*ww, 100.0f*hh, ww, hh);
+
+  // http://thenewcode.com/1068/Making-Arrows-in-SVG
+  fprintf(fp, "<defs>\n");
+  fprintf(fp, "<marker id=\"arrowhead\" markerWidth=\"2\" markerHeight=\"3\" \n");
+  fprintf(fp, "refX=\"2\" refY=\"1.5\" orient=\"auto\">\n");
+  fprintf(fp, "<polygon points=\"0 0, 2 1.5, 0 3\" />\n");
+  fprintf(fp, "</marker>\n");
+  fprintf(fp, "</defs>\n");
+
+  fprintf(fp, "<rect fill=\"#fff\" stroke=\"#000\" x=\"0\" y=\"0\" width=\"%f\" height=\"%f\"/>\n", ww, hh);
+
+  for (std::list< std::vector<OrthancStone::ScenePoint2D> >::const_iterator
+         it = contours.begin(); it != contours.end(); ++it)
+  {
+    for (size_t i = 0; i + 1 < it->size(); i++)
+    {
+      float x1 = (*it)[i].GetX();
+      float x2 = (*it)[i + 1].GetX();
+      float y1 = (*it)[i].GetY();
+      float y2 = (*it)[i + 1].GetY();
+      
+      fprintf(fp, "<line x1=\"%f\" y1=\"%f\" x2=\"%f\" y2=\"%f\" stroke=\"blue\" stroke-width=\"0.05\" marker-end=\"url(#arrowhead)\"/>\n", x1, y1, x2, y2);
+    }
+  }
+  fprintf(fp, "</svg>\n");
+  
+  fclose(fp);
+}
+#endif
+
+
+TEST(UnionOfRectangles, Complex)
+{
+  {
+    std::list<OrthancStone::Extent2D>  rectangles;
+    rectangles.push_back(OrthancStone::Extent2D(1, 4, 4, 6));
+    rectangles.push_back(OrthancStone::Extent2D(4, 6, 7, 8));
+    rectangles.push_back(OrthancStone::Extent2D(4, 2, 7, 4));
+    rectangles.push_back(OrthancStone::Extent2D(7, 4, 10, 6));
+  
+    std::list< std::vector<OrthancStone::ScenePoint2D> > contours;
+    OrthancStone::UnionOfRectangles::Apply(contours, rectangles);
+
+    ASSERT_EQ(1u, contours.size());
+    ASSERT_EQ(17u, contours.front().size());
+    ASSERT_TRUE(contours.front()[0].IsEqual(OrthancStone::ScenePoint2D(4, 8)));
+    ASSERT_TRUE(contours.front()[1].IsEqual(OrthancStone::ScenePoint2D(7, 8)));
+    ASSERT_TRUE(contours.front()[2].IsEqual(OrthancStone::ScenePoint2D(7, 6)));
+    ASSERT_TRUE(contours.front()[3].IsEqual(OrthancStone::ScenePoint2D(10, 6)));
+    ASSERT_TRUE(contours.front()[4].IsEqual(OrthancStone::ScenePoint2D(10, 4)));
+    ASSERT_TRUE(contours.front()[5].IsEqual(OrthancStone::ScenePoint2D(7, 4)));
+    ASSERT_TRUE(contours.front()[6].IsEqual(OrthancStone::ScenePoint2D(7, 2)));
+    ASSERT_TRUE(contours.front()[7].IsEqual(OrthancStone::ScenePoint2D(4, 2)));
+    ASSERT_TRUE(contours.front()[8].IsEqual(OrthancStone::ScenePoint2D(4, 4)));
+    ASSERT_TRUE(contours.front()[9].IsEqual(OrthancStone::ScenePoint2D(7, 4)));
+    ASSERT_TRUE(contours.front()[10].IsEqual(OrthancStone::ScenePoint2D(7, 6)));
+    ASSERT_TRUE(contours.front()[11].IsEqual(OrthancStone::ScenePoint2D(4, 6)));
+    ASSERT_TRUE(contours.front()[12].IsEqual(OrthancStone::ScenePoint2D(4, 4)));
+    ASSERT_TRUE(contours.front()[13].IsEqual(OrthancStone::ScenePoint2D(1, 4)));
+    ASSERT_TRUE(contours.front()[14].IsEqual(OrthancStone::ScenePoint2D(1, 6)));
+    ASSERT_TRUE(contours.front()[15].IsEqual(OrthancStone::ScenePoint2D(4, 6)));
+    ASSERT_TRUE(contours.front()[16].IsEqual(OrthancStone::ScenePoint2D(4, 8)));
+  }
+  
+  {
+    std::list<OrthancStone::Extent2D>  rectangles;
+    rectangles.push_back(OrthancStone::Extent2D(1, 4, 4, 6));
+    rectangles.push_back(OrthancStone::Extent2D(4, 4, 7, 6));
+  
+    std::list< std::vector<OrthancStone::ScenePoint2D> > contours;
+    OrthancStone::UnionOfRectangles::Apply(contours, rectangles);
+
+    ASSERT_EQ(1u, contours.size());
+    ASSERT_EQ(5u, contours.front().size());
+    ASSERT_TRUE(contours.front()[0].IsEqual(OrthancStone::ScenePoint2D(1, 6)));
+    ASSERT_TRUE(contours.front()[1].IsEqual(OrthancStone::ScenePoint2D(7, 6)));
+    ASSERT_TRUE(contours.front()[2].IsEqual(OrthancStone::ScenePoint2D(7, 4)));
+    ASSERT_TRUE(contours.front()[3].IsEqual(OrthancStone::ScenePoint2D(1, 4)));
+    ASSERT_TRUE(contours.front()[4].IsEqual(OrthancStone::ScenePoint2D(1, 6)));
+  }
+  
+  {
+    std::list<OrthancStone::Extent2D>  rectangles;
+    rectangles.push_back(OrthancStone::Extent2D(1, 4, 4, 6));
+    rectangles.push_back(OrthancStone::Extent2D(1, 6, 4, 8));
+  
+    std::list< std::vector<OrthancStone::ScenePoint2D> > contours;
+    OrthancStone::UnionOfRectangles::Apply(contours, rectangles);
+
+    ASSERT_EQ(1u, contours.size());
+    ASSERT_EQ(5u, contours.front().size());
+    ASSERT_TRUE(contours.front()[0].IsEqual(OrthancStone::ScenePoint2D(1, 8)));
+    ASSERT_TRUE(contours.front()[1].IsEqual(OrthancStone::ScenePoint2D(4, 8)));
+    ASSERT_TRUE(contours.front()[2].IsEqual(OrthancStone::ScenePoint2D(4, 4)));
+    ASSERT_TRUE(contours.front()[3].IsEqual(OrthancStone::ScenePoint2D(1, 4)));
+    ASSERT_TRUE(contours.front()[4].IsEqual(OrthancStone::ScenePoint2D(1, 8)));
+  }
+  
+  {
+    std::list<OrthancStone::Extent2D>  rectangles;
+    rectangles.push_back(OrthancStone::Extent2D(1, 1, 2, 2));
+    rectangles.push_back(OrthancStone::Extent2D(4, 4, 7, 6));
+    rectangles.push_back(OrthancStone::Extent2D(4, 6, 7, 8));
+
+    std::list< std::vector<OrthancStone::ScenePoint2D> > contours;
+    OrthancStone::UnionOfRectangles::Apply(contours, rectangles);
+    
+    ASSERT_EQ(2u, contours.size());
+    
+    ASSERT_EQ(5u, contours.front().size());
+    ASSERT_TRUE(contours.front()[0].IsEqual(OrthancStone::ScenePoint2D(4, 8)));
+    ASSERT_TRUE(contours.front()[1].IsEqual(OrthancStone::ScenePoint2D(7, 8)));
+    ASSERT_TRUE(contours.front()[2].IsEqual(OrthancStone::ScenePoint2D(7, 4)));
+    ASSERT_TRUE(contours.front()[3].IsEqual(OrthancStone::ScenePoint2D(4, 4)));
+    ASSERT_TRUE(contours.front()[4].IsEqual(OrthancStone::ScenePoint2D(4, 8)));
+    
+    ASSERT_EQ(5u, contours.back().size());
+    ASSERT_TRUE(contours.back()[0].IsEqual(OrthancStone::ScenePoint2D(1, 2)));
+    ASSERT_TRUE(contours.back()[1].IsEqual(OrthancStone::ScenePoint2D(2, 2)));
+    ASSERT_TRUE(contours.back()[2].IsEqual(OrthancStone::ScenePoint2D(2, 1)));
+    ASSERT_TRUE(contours.back()[3].IsEqual(OrthancStone::ScenePoint2D(1, 1)));
+    ASSERT_TRUE(contours.back()[4].IsEqual(OrthancStone::ScenePoint2D(1, 2)));
+  }
+
+#if 0
+  {
+    std::list<OrthancStone::Extent2D>  rectangles;
+    rectangles.push_back(OrthancStone::Extent2D(1, 4, 4, 6));
+    rectangles.push_back(OrthancStone::Extent2D(6, 4, 9, 6));
+    rectangles.push_back(OrthancStone::Extent2D(4, 6, 7, 8));
+    rectangles.push_back(OrthancStone::Extent2D(4, 2, 7, 6));
+
+    std::list< std::vector<OrthancStone::ScenePoint2D> > contours;
+    OrthancStone::UnionOfRectangles::Apply(contours, rectangles);
+
+    SaveSvg(contours);
+    
+    ASSERT_EQ(1u, contours.size());
+    ASSERT_EQ(13u, contours.front().size());
+    ASSERT_TRUE(contours.front()[0].IsEqual(OrthancStone::ScenePoint2D(4, 8)));
+    ASSERT_TRUE(contours.front()[1].IsEqual(OrthancStone::ScenePoint2D(7, 8)));
+    ASSERT_TRUE(contours.front()[2].IsEqual(OrthancStone::ScenePoint2D(7, 6)));
+    ASSERT_TRUE(contours.front()[3].IsEqual(OrthancStone::ScenePoint2D(9, 6)));
+    ASSERT_TRUE(contours.front()[4].IsEqual(OrthancStone::ScenePoint2D(7, 6)));
+    ASSERT_TRUE(contours.front()[5].IsEqual(OrthancStone::ScenePoint2D(4, 4)));
+    ASSERT_TRUE(contours.front()[6].IsEqual(OrthancStone::ScenePoint2D(1, 4)));
+    ASSERT_TRUE(contours.front()[7].IsEqual(OrthancStone::ScenePoint2D()));
+    ASSERT_TRUE(contours.front()[8].IsEqual(OrthancStone::ScenePoint2D()));
+    ASSERT_TRUE(contours.front()[9].IsEqual(OrthancStone::ScenePoint2D()));
+    ASSERT_TRUE(contours.front()[10].IsEqual(OrthancStone::ScenePoint2D()));
+    ASSERT_TRUE(contours.front()[11].IsEqual(OrthancStone::ScenePoint2D()));
+    ASSERT_TRUE(contours.front()[12].IsEqual(OrthancStone::ScenePoint2D()));
+  }
+#endif
+}