--- a/Framework/Toolbox/GeometryToolbox.cpp	Tue Jan 16 17:44:16 2018 +0100
+++ b/Framework/Toolbox/GeometryToolbox.cpp	Thu Jan 18 14:42:33 2018 +0100
@@ -259,40 +259,40 @@
       // (2005). This is a direct, non-optimized translation of Algorithm
       // 2 in the paper.
 
-      static uint8_t tab1[16] = { 255 /* none */,
-                                  0,
-                                  0,
-                                  1,
-                                  1,
-                                  255 /* na */,
-                                  0,
-                                  2,
-                                  2,
-                                  0,
-                                  255 /* na */,
-                                  1,
-                                  1,
-                                  0,
-                                  0,
-                                  255 /* none */ };
+      static const uint8_t tab1[16] = { 255 /* none */,
+                                        0,
+                                        0,
+                                        1,
+                                        1,
+                                        255 /* na */,
+                                        0,
+                                        2,
+                                        2,
+                                        0,
+                                        255 /* na */,
+                                        1,
+                                        1,
+                                        0,
+                                        0,
+                                        255 /* none */ };
 
 
-      static uint8_t tab2[16] = { 255 /* none */,
-                                  3,
-                                  1,
-                                  3,
-                                  2,
-                                  255 /* na */,
-                                  2,
-                                  3,
-                                  3,
-                                  2,
-                                  255 /* na */,
-                                  2,
-                                  3,
-                                  1,
-                                  3,
-                                  255 /* none */ };
+      static const uint8_t tab2[16] = { 255 /* none */,
+                                        3,
+                                        1,
+                                        3,
+                                        2,
+                                        255 /* na */,
+                                        2,
+                                        3,
+                                        3,
+                                        2,
+                                        255 /* na */,
+                                        2,
+                                        3,
+                                        1,
+                                        3,
+                                        255 /* none */ };
 
       // Create the coordinates of the rectangle
       Vector x[4];
@@ -378,5 +378,59 @@
         spacingY = 1;
       }
     }
+
+    
+    Matrix CreateRotationMatrixAlongX(double a)
+    {
+      // Rotate along X axis (R_x)
+      // https://en.wikipedia.org/wiki/Rotation_matrix#Basic_rotations
+      Matrix r(3, 3);
+      r(0,0) = 1;
+      r(0,1) = 0;
+      r(0,2) = 0;
+      r(1,0) = 0;
+      r(1,1) = cos(a);
+      r(1,2) = -sin(a);
+      r(2,0) = 0;
+      r(2,1) = sin(a);
+      r(2,2) = cos(a);
+      return r;
+    }
+    
+
+    Matrix CreateRotationMatrixAlongY(double a)
+    {
+      // Rotate along Y axis (R_y)
+      // https://en.wikipedia.org/wiki/Rotation_matrix#Basic_rotations
+      Matrix r(3, 3);
+      r(0,0) = cos(a);
+      r(0,1) = 0;
+      r(0,2) = sin(a);
+      r(1,0) = 0;
+      r(1,1) = 1;
+      r(1,2) = 0;
+      r(2,0) = -sin(a);
+      r(2,1) = 0;
+      r(2,2) = cos(a);
+      return r;
+    }
+
+
+    Matrix CreateRotationMatrixAlongZ(double a)
+    {
+      // Rotate along Z axis (R_z)
+      // https://en.wikipedia.org/wiki/Rotation_matrix#Basic_rotations
+      Matrix r(3, 3);
+      r(0,0) = cos(a);
+      r(0,1) = -sin(a);
+      r(0,2) = 0;
+      r(1,0) = sin(a);
+      r(1,1) = cos(a);
+      r(1,2) = 0;
+      r(2,0) = 0;
+      r(2,1) = 0;
+      r(2,2) = 1;
+      return r;
+    }    
   }
 }

--- a/Framework/Toolbox/GeometryToolbox.h	Tue Jan 16 17:44:16 2018 +0100
+++ b/Framework/Toolbox/GeometryToolbox.h	Thu Jan 18 14:42:33 2018 +0100
@@ -28,12 +28,14 @@
 #  include <boost/serialization/array_wrapper.hpp>
 #endif
 
+#include <boost/numeric/ublas/matrix.hpp>
 #include <boost/numeric/ublas/vector.hpp>
 
 #include <Core/DicomFormat/DicomMap.h>
 
 namespace OrthancStone
 {
+  typedef boost::numeric::ublas::matrix<double>   Matrix;
   typedef boost::numeric::ublas::vector<double>   Vector;
 
   namespace GeometryToolbox
@@ -118,5 +120,11 @@
     {
       return boost::numeric::ublas::inner_prod(point, normal);
     }
+
+    Matrix CreateRotationMatrixAlongX(double a);
+
+    Matrix CreateRotationMatrixAlongY(double a);
+
+    Matrix CreateRotationMatrixAlongZ(double a);    
   };
 }

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Framework/Toolbox/OrientedBoundingBox.cpp	Thu Jan 18 14:42:33 2018 +0100
@@ -0,0 +1,267 @@
+/**
+ * Stone of Orthanc
+ * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics
+ * Department, University Hospital of Liege, Belgium
+ * Copyright (C) 2017-2018 Osimis S.A., Belgium
+ *
+ * This program is free software: you can redistribute it and/or
+ * modify it under the terms of the GNU Affero 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
+ * Affero General Public License for more details.
+ * 
+ * You should have received a copy of the GNU Affero General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ **/
+
+
+#include "OrientedBoundingBox.h"
+
+#include <Core/OrthancException.h>
+
+#include <cassert>
+
+namespace OrthancStone
+{
+  static bool IntersectPlaneAndSegment(Vector& p,
+                                       const Vector& normal,
+                                       double d,
+                                       const Vector& edgeFrom,
+                                       const Vector& edgeTo)
+  {
+    // http://geomalgorithms.com/a05-_intersect-1.html#Line-Plane-Intersection
+
+    // Check for parallel line and plane
+    Vector direction = edgeTo - edgeFrom;
+    double denominator = boost::numeric::ublas::inner_prod(direction, normal);
+
+    if (fabs(denominator) < 100.0 * std::numeric_limits<double>::epsilon())
+    {
+      return false;
+    }
+    else
+    {
+      // Compute intersection
+      double t = -(normal[0] * edgeFrom[0] + 
+                   normal[1] * edgeFrom[1] + 
+                   normal[2] * edgeFrom[2] + d) / denominator;
+
+      if (t >= 0 && t <= 1)
+      {
+        // The intersection lies inside edge segment
+        p = edgeFrom + t * direction;
+        return true;
+      }
+      else
+      {
+        return false;
+      }
+    }
+  }
+
+
+  OrientedBoundingBox::OrientedBoundingBox(const ImageBuffer3D& image)
+  {
+    unsigned int n = image.GetDepth();
+    if (n < 1)
+    {
+      throw Orthanc::OrthancException(Orthanc::ErrorCode_IncompatibleImageSize);      
+    }
+
+    const CoordinateSystem3D& geometry = image.GetAxialGeometry();
+    Vector dim = image.GetVoxelDimensions(VolumeProjection_Axial);
+
+    u_ = geometry.GetAxisX();
+    v_ = geometry.GetAxisY();
+    w_ = geometry.GetNormal();
+
+    hu_ = static_cast<double>(image.GetWidth() * dim[0] / 2.0);
+    hv_ = static_cast<double>(image.GetHeight() * dim[1] / 2.0);
+    hw_ = static_cast<double>(image.GetDepth() * dim[2] / 2.0);
+      
+    c_ = (geometry.GetOrigin() + 
+          (hu_ - dim[0] / 2.0) * u_ +
+          (hv_ - dim[1] / 2.0) * v_ +
+          (hw_ - dim[2] / 2.0) * w_);
+  }
+
+
+  bool OrientedBoundingBox::HasIntersectionWithPlane(std::vector<Vector>& points,
+                                                     const Vector& normal,
+                                                     double d) const
+  {
+    assert(normal.size() == 3);
+
+    double r = (hu_ * fabs(boost::numeric::ublas::inner_prod(normal, u_)) +
+                hv_ * fabs(boost::numeric::ublas::inner_prod(normal, v_)) +
+                hw_ * fabs(boost::numeric::ublas::inner_prod(normal, w_)));
+
+    double s = boost::numeric::ublas::inner_prod(normal, c_) + d;
+
+    if (fabs(s) >= r)
+    {
+      // No intersection, or intersection is reduced to a single point
+      return false;
+    }
+    else
+    {
+      Vector p;
+
+      // Loop over all the 12 edges (segments) of the oriented
+      // bounding box, and check whether they intersect the plane
+        
+      // X-aligned edges
+      if (IntersectPlaneAndSegment(p, normal, d,
+                                   c_ - u_ * hu_ - v_ * hv_ - w_ * hw_,
+                                   c_ + u_ * hu_ - v_ * hv_ - w_ * hw_))
+      {
+        points.push_back(p);
+      }
+
+      if (IntersectPlaneAndSegment(p, normal, d,
+                                   c_ - u_ * hu_ + v_ * hv_ - w_ * hw_,
+                                   c_ + u_ * hu_ + v_ * hv_ - w_ * hw_))
+      {
+        points.push_back(p);
+      }
+
+      if (IntersectPlaneAndSegment(p, normal, d,
+                                   c_ - u_ * hu_ - v_ * hv_ + w_ * hw_,
+                                   c_ + u_ * hu_ - v_ * hv_ + w_ * hw_))
+      {
+        points.push_back(p);
+      }
+
+      if (IntersectPlaneAndSegment(p, normal, d,
+                                   c_ - u_ * hu_ + v_ * hv_ + w_ * hw_,
+                                   c_ + u_ * hu_ + v_ * hv_ + w_ * hw_))
+      {
+        points.push_back(p);
+      }
+
+      // Y-aligned edges
+      if (IntersectPlaneAndSegment(p, normal, d,
+                                   c_ - u_ * hu_ - v_ * hv_ - w_ * hw_,
+                                   c_ - u_ * hu_ + v_ * hv_ - w_ * hw_))
+      {
+        points.push_back(p);
+      }
+
+      if (IntersectPlaneAndSegment(p, normal, d,
+                                   c_ + u_ * hu_ - v_ * hv_ - w_ * hw_,
+                                   c_ + u_ * hu_ + v_ * hv_ - w_ * hw_))
+      {
+        points.push_back(p);
+      }
+
+      if (IntersectPlaneAndSegment(p, normal, d,
+                                   c_ - u_ * hu_ - v_ * hv_ + w_ * hw_,
+                                   c_ - u_ * hu_ + v_ * hv_ + w_ * hw_))
+      {
+        points.push_back(p);
+      }
+
+      if (IntersectPlaneAndSegment(p, normal, d,
+                                   c_ + u_ * hu_ - v_ * hv_ + w_ * hw_,
+                                   c_ + u_ * hu_ + v_ * hv_ + w_ * hw_))
+      {
+        points.push_back(p);
+      }
+
+      // Z-aligned edges
+      if (IntersectPlaneAndSegment(p, normal, d,
+                                   c_ - u_ * hu_ - v_ * hv_ - w_ * hw_,
+                                   c_ - u_ * hu_ - v_ * hv_ + w_ * hw_))
+      {
+        points.push_back(p);
+      }
+
+      if (IntersectPlaneAndSegment(p, normal, d,
+                                   c_ + u_ * hu_ - v_ * hv_ - w_ * hw_,
+                                   c_ + u_ * hu_ - v_ * hv_ + w_ * hw_))
+      {
+        points.push_back(p);
+      }
+
+      if (IntersectPlaneAndSegment(p, normal, d,
+                                   c_ - u_ * hu_ + v_ * hv_ - w_ * hw_,
+                                   c_ - u_ * hu_ + v_ * hv_ + w_ * hw_))
+      {
+        points.push_back(p);
+      }
+
+      if (IntersectPlaneAndSegment(p, normal, d,
+                                   c_ + u_ * hu_ + v_ * hv_ - w_ * hw_,
+                                   c_ + u_ * hu_ + v_ * hv_ + w_ * hw_))
+      {
+        points.push_back(p);
+      }
+
+      return true;
+    }
+  }
+
+
+  bool OrientedBoundingBox::HasIntersection(std::vector<Vector>& points,
+                                            const CoordinateSystem3D& plane) const
+  {
+    // From the vector equation of a 3D plane (specified by origin
+    // and normal), to the general equation of a 3D plane (which
+    // looses information about the origin of the coordinate system)
+    const Vector& normal = plane.GetNormal();
+    const Vector& origin = plane.GetOrigin();
+    double d = -(normal[0] * origin[0] + normal[1] * origin[1] + normal[2] * origin[2]);
+
+    return HasIntersectionWithPlane(points, normal, d);
+  }
+  
+
+  bool OrientedBoundingBox::Contains(const Vector& p) const
+  {
+    assert(p.size() == 3);
+
+    const Vector q = p - c_;
+
+    return (fabs(boost::numeric::ublas::inner_prod(q, u_)) <= hu_ &&
+            fabs(boost::numeric::ublas::inner_prod(q, v_)) <= hv_ &&
+            fabs(boost::numeric::ublas::inner_prod(q, w_)) <= hw_);
+  }
+
+  
+  void OrientedBoundingBox::FromInternalCoordinates(Vector& target,
+                                                    double x,
+                                                    double y,
+                                                    double z) const
+  {
+    target = (c_ +
+              u_ * 2.0 * hu_ * (x - 0.5) +
+              v_ * 2.0 * hv_ * (y - 0.5) +
+              w_ * 2.0 * hw_ * (z - 0.5));
+  }
+
+  
+  void OrientedBoundingBox::FromInternalCoordinates(Vector& target,
+                                                    const Vector& source) const
+  {
+    assert(source.size() == 3);
+    FromInternalCoordinates(target, source[0], source[1], source[2]);
+  }
+
+  
+  void OrientedBoundingBox::ToInternalCoordinates(Vector& target,
+                                                  const Vector& source) const
+  {
+    assert(source.size() == 3);
+    const Vector q = source - c_;
+
+    double x = boost::numeric::ublas::inner_prod(q, u_) / (2.0 * hu_) + 0.5;
+    double y = boost::numeric::ublas::inner_prod(q, v_) / (2.0 * hv_) + 0.5;
+    double z = boost::numeric::ublas::inner_prod(q, w_) / (2.0 * hw_) + 0.5;
+
+    GeometryToolbox::AssignVector(target, x, y, z);
+  }
+}

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Framework/Toolbox/OrientedBoundingBox.h	Thu Jan 18 14:42:33 2018 +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-2018 Osimis S.A., Belgium
+ *
+ * This program is free software: you can redistribute it and/or
+ * modify it under the terms of the GNU Affero 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
+ * Affero General Public License for more details.
+ * 
+ * You should have received a copy of the GNU Affero General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ **/
+
+
+#pragma once
+
+#include "GeometryToolbox.h"
+#include "../Volumes/ImageBuffer3D.h"
+
+namespace OrthancStone
+{
+  class OrientedBoundingBox : public boost::noncopyable
+  {
+  private:
+    Vector  c_;   // center
+    Vector  u_;   // normalized width vector
+    Vector  v_;   // normalized height vector
+    Vector  w_;   // normalized depth vector
+    double  hu_;  // half width
+    double  hv_;  // half height
+    double  hw_;  // half depth
+
+  public:
+    OrientedBoundingBox(const ImageBuffer3D& image);
+
+    bool HasIntersectionWithPlane(std::vector<Vector>& points,
+                                  const Vector& normal,
+                                  double d) const;
+
+    bool HasIntersection(std::vector<Vector>& points,
+                         const CoordinateSystem3D& plane) const;
+
+    bool Contains(const Vector& p) const;
+
+    void FromInternalCoordinates(Vector& target,
+                                 double x,
+                                 double y,
+                                 double z) const;
+
+    void FromInternalCoordinates(Vector& target,
+                                 const Vector& source) const;
+
+    void ToInternalCoordinates(Vector& target,
+                               const Vector& source) const;
+
+    const Vector& GetCenter() const
+    {
+      return c_;
+    }
+  };
+}
+

--- a/Resources/CMake/OrthancStoneConfiguration.cmake	Tue Jan 16 17:44:16 2018 +0100
+++ b/Resources/CMake/OrthancStoneConfiguration.cmake	Thu Jan 18 14:42:33 2018 +0100
@@ -182,6 +182,7 @@
   ${ORTHANC_STONE_DIR}/Framework/Toolbox/Extent2D.cpp
   ${ORTHANC_STONE_DIR}/Framework/Toolbox/GeometryToolbox.cpp
   ${ORTHANC_STONE_DIR}/Framework/Toolbox/MessagingToolbox.cpp
+  ${ORTHANC_STONE_DIR}/Framework/Toolbox/OrientedBoundingBox.cpp
   ${ORTHANC_STONE_DIR}/Framework/Toolbox/OrthancSlicesLoader.cpp
   ${ORTHANC_STONE_DIR}/Framework/Toolbox/ParallelSlices.cpp
   ${ORTHANC_STONE_DIR}/Framework/Toolbox/ParallelSlicesCursor.cpp