--- a/Framework/Layers/CircleMeasureTracker.cpp	Fri Feb 09 16:00:29 2018 +0100
+++ b/Framework/Layers/CircleMeasureTracker.cpp	Wed Feb 14 08:58:31 2018 +0100
@@ -59,7 +59,7 @@
     double y = (y1_ + y2_) / 2.0;
 
     Vector tmp;
-    GeometryToolbox::AssignVector(tmp, x2_ - x1_, y2_ - y1_);
+    LinearAlgebra::AssignVector(tmp, x2_ - x1_, y2_ - y1_);
     double r = boost::numeric::ublas::norm_2(tmp) / 2.0;
 
     context.SetSourceColor(color_[0], color_[1], color_[2]);

--- a/Framework/Toolbox/CoordinateSystem3D.cpp	Fri Feb 09 16:00:29 2018 +0100
+++ b/Framework/Toolbox/CoordinateSystem3D.cpp	Wed Feb 14 08:58:31 2018 +0100
@@ -21,6 +21,7 @@
 
 #include "CoordinateSystem3D.h"
 
+#include "LinearAlgebra.h"
 #include "GeometryToolbox.h"
 
 #include <Core/Logging.h>
@@ -36,8 +37,8 @@
     // product of each direction cosine vector with itself shall be
     // unity."
     // http://dicom.nema.org/medical/dicom/current/output/chtml/part03/sect_C.7.6.2.html
-    if (!GeometryToolbox::IsNear(boost::numeric::ublas::norm_2(axisX_), 1.0) ||
-        !GeometryToolbox::IsNear(boost::numeric::ublas::norm_2(axisY_), 1.0))
+    if (!LinearAlgebra::IsNear(boost::numeric::ublas::norm_2(axisX_), 1.0) ||
+        !LinearAlgebra::IsNear(boost::numeric::ublas::norm_2(axisY_), 1.0))
     {
       throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
     }
@@ -47,25 +48,25 @@
     // column direction cosine vectors shall be orthogonal, i.e.,
     // their dot product shall be zero."
     // http://dicom.nema.org/medical/dicom/current/output/chtml/part03/sect_C.7.6.2.html
-    if (!GeometryToolbox::IsCloseToZero(boost::numeric::ublas::inner_prod(axisX_, axisY_)))
+    if (!LinearAlgebra::IsCloseToZero(boost::numeric::ublas::inner_prod(axisX_, axisY_)))
     {
       throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
     }
 
-    GeometryToolbox::CrossProduct(normal_, axisX_, axisY_);
+    LinearAlgebra::CrossProduct(normal_, axisX_, axisY_);
 
     d_ = -(normal_[0] * origin_[0] + normal_[1] * origin_[1] + normal_[2] * origin_[2]);
 
     // Just a sanity check, it should be useless by construction
-    assert(GeometryToolbox::IsNear(boost::numeric::ublas::norm_2(normal_), 1.0));
+    assert(LinearAlgebra::IsNear(boost::numeric::ublas::norm_2(normal_), 1.0));
   }
 
 
   void CoordinateSystem3D::SetupCanonical()
   {
-    GeometryToolbox::AssignVector(origin_, 0, 0, 0);
-    GeometryToolbox::AssignVector(axisX_, 1, 0, 0);
-    GeometryToolbox::AssignVector(axisY_, 0, 1, 0);
+    LinearAlgebra::AssignVector(origin_, 0, 0, 0);
+    LinearAlgebra::AssignVector(axisX_, 1, 0, 0);
+    LinearAlgebra::AssignVector(axisY_, 0, 1, 0);
     CheckAndComputeNormal();
   }
 
@@ -88,8 +89,8 @@
     std::string tmpOrientation = Orthanc::Toolbox::StripSpaces(imageOrientationPatient);
 
     Vector orientation;
-    if (!GeometryToolbox::ParseVector(origin_, tmpPosition) ||
-        !GeometryToolbox::ParseVector(orientation, tmpOrientation) ||
+    if (!LinearAlgebra::ParseVector(origin_, tmpPosition) ||
+        !LinearAlgebra::ParseVector(orientation, tmpOrientation) ||
         origin_.size() != 3 ||
         orientation.size() != 6)
     {

--- a/Framework/Toolbox/DicomFrameConverter.cpp	Fri Feb 09 16:00:29 2018 +0100
+++ b/Framework/Toolbox/DicomFrameConverter.cpp	Wed Feb 14 08:58:31 2018 +0100
@@ -21,7 +21,7 @@
 
 #include "DicomFrameConverter.h"
 
-#include "GeometryToolbox.h"
+#include "LinearAlgebra.h"
 
 #include <Core/Images/Image.h>
 #include <Core/Images/ImageProcessing.h>
@@ -48,8 +48,8 @@
     SetDefaultParameters();
 
     Vector c, w;
-    if (GeometryToolbox::ParseVector(c, dicom, Orthanc::DICOM_TAG_WINDOW_CENTER) &&
-        GeometryToolbox::ParseVector(w, dicom, Orthanc::DICOM_TAG_WINDOW_WIDTH) &&
+    if (LinearAlgebra::ParseVector(c, dicom, Orthanc::DICOM_TAG_WINDOW_CENTER) &&
+        LinearAlgebra::ParseVector(w, dicom, Orthanc::DICOM_TAG_WINDOW_WIDTH) &&
         c.size() > 0 && 
         w.size() > 0)
     {

--- a/Framework/Toolbox/DicomStructureSet.cpp	Fri Feb 09 16:00:29 2018 +0100
+++ b/Framework/Toolbox/DicomStructureSet.cpp	Wed Feb 14 08:58:31 2018 +0100
@@ -124,7 +124,7 @@
   {
     std::string value;
     return (dataset.GetStringValue(value, tag) &&
-            GeometryToolbox::ParseVector(target, value));
+            LinearAlgebra::ParseVector(target, value));
   }
 
 
@@ -132,9 +132,9 @@
   {
     if (hasSlice_)
     {
-      if (!GeometryToolbox::IsNear(GeometryToolbox::ProjectAlongNormal(v, geometry_.GetNormal()),
-                                   projectionAlongNormal_, 
-                                   sliceThickness_ / 2.0 /* in mm */))
+      if (!LinearAlgebra::IsNear(GeometryToolbox::ProjectAlongNormal(v, geometry_.GetNormal()),
+                                 projectionAlongNormal_, 
+                                 sliceThickness_ / 2.0 /* in mm */))
       {
         LOG(ERROR) << "This RT-STRUCT contains a point that is off the slice of its instance";
         throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);          
@@ -203,8 +203,8 @@
     
     double d = GeometryToolbox::ProjectAlongNormal(slice.GetOrigin(), geometry_.GetNormal());
 
-    return (GeometryToolbox::IsNear(d, projectionAlongNormal_,
-                                    sliceThickness_ / 2.0));
+    return (LinearAlgebra::IsNear(d, projectionAlongNormal_,
+                                  sliceThickness_ / 2.0));
   }
 
   
@@ -468,7 +468,7 @@
                      DICOM_TAG_CONTOUR_DATA));
 
         Vector points;
-        if (!GeometryToolbox::ParseVector(points, slicesData) ||
+        if (!LinearAlgebra::ParseVector(points, slicesData) ||
             points.size() != 3 * countPoints)
         {
           throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);          
@@ -497,7 +497,7 @@
     const Structure& structure = GetStructure(index);
 
     Vector center;
-    GeometryToolbox::AssignVector(center, 0, 0, 0);
+    LinearAlgebra::AssignVector(center, 0, 0, 0);
     if (structure.polygons_.empty())
     {
       return center;
@@ -614,7 +614,7 @@
     std::string s;
     Vector v;
     if (dataset.CopyToString(s, Orthanc::DICOM_TAG_SLICE_THICKNESS, false) &&
-        GeometryToolbox::ParseVector(v, s) &&
+        LinearAlgebra::ParseVector(v, s) &&
         v.size() > 0)
     {
       thickness = v[0];
@@ -657,7 +657,7 @@
     if (referencedSlices_.empty())
     {
       Vector v;
-      GeometryToolbox::AssignVector(v, 0, 0, 1);
+      LinearAlgebra::AssignVector(v, 0, 0, 1);
       return v;
     }
     else

--- a/Framework/Toolbox/GeometryToolbox.cpp	Fri Feb 09 16:00:29 2018 +0100
+++ b/Framework/Toolbox/GeometryToolbox.cpp	Wed Feb 14 08:58:31 2018 +0100
@@ -23,142 +23,20 @@
 
 #include <Core/Logging.h>
 #include <Core/OrthancException.h>
-#include <Core/Toolbox.h>
 
-#include <stdio.h>
-#include <boost/lexical_cast.hpp>
+#include <cassert>
 
 namespace OrthancStone
 {
   namespace GeometryToolbox
   {
-    void Print(const Vector& v)
-    {
-      for (size_t i = 0; i < v.size(); i++)
-      {
-        printf("%g\n", v[i]);
-        //printf("%8.2f\n", v[i]);
-      }
-      printf("\n");
-    }
-
-
-    void Print(const Matrix& m)
-    {
-      for (size_t i = 0; i < m.size1(); i++)
-      {
-        for (size_t j = 0; j < m.size2(); j++)
-        {
-          printf("%g  ", m(i,j));
-          //printf("%8.2f  ", m(i,j));
-        }
-        printf("\n");        
-      }
-      printf("\n");        
-    }
-
-
-    bool ParseVector(Vector& target,
-                     const std::string& value)
-    {
-      std::vector<std::string> items;
-      Orthanc::Toolbox::TokenizeString(items, value, '\\');
-
-      target.resize(items.size());
-
-      for (size_t i = 0; i < items.size(); i++)
-      {
-        try
-        {
-          target[i] = boost::lexical_cast<double>(Orthanc::Toolbox::StripSpaces(items[i]));
-        }
-        catch (boost::bad_lexical_cast&)
-        {
-          target.clear();
-          return false;
-        }
-      }
-
-      return true;
-    }
-
-
-    bool ParseVector(Vector& target,
-                     const Orthanc::DicomMap& dataset,
-                     const Orthanc::DicomTag& tag)
-    {
-      std::string value;
-      return (dataset.CopyToString(value, tag, false) &&
-              ParseVector(target, value));
-    }
-
-
-    void AssignVector(Vector& v,
-                      double v1,
-                      double v2)
-    {
-      v.resize(2);
-      v[0] = v1;
-      v[1] = v2;
-    }
-
-
-    void AssignVector(Vector& v,
-                      double v1,
-                      double v2,
-                      double v3)
-    {
-      v.resize(3);
-      v[0] = v1;
-      v[1] = v2;
-      v[2] = v3;
-    }
-
-
-    bool IsNear(double x,
-                double y)
-    {
-      // As most input is read as single-precision numbers, we take the
-      // epsilon machine for float32 into consideration to compare numbers
-      return IsNear(x, y, 10.0 * std::numeric_limits<float>::epsilon());
-    }
-
-
-    void NormalizeVector(Vector& u)
-    {
-      double norm = boost::numeric::ublas::norm_2(u);
-      if (!IsCloseToZero(norm))
-      {
-        u = u / norm;
-      }
-    }
-
-
-    void CrossProduct(Vector& result,
-                      const Vector& u,
-                      const Vector& v)
-    {
-      if (u.size() != 3 ||
-          v.size() != 3)
-      {
-        throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
-      }
-
-      result.resize(3);
-
-      result[0] = u[1] * v[2] - u[2] * v[1];
-      result[1] = u[2] * v[0] - u[0] * v[2];
-      result[2] = u[0] * v[1] - u[1] * v[0];
-    }
-
-
     void ProjectPointOntoPlane(Vector& result,
                                const Vector& point,
                                const Vector& planeNormal,
                                const Vector& planeOrigin)
     {
       double norm =  boost::numeric::ublas::norm_2(planeNormal);
-      if (IsCloseToZero(norm))
+      if (LinearAlgebra::IsCloseToZero(norm))
       {
         // Division by zero
         throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
@@ -187,8 +65,8 @@
       double normU = boost::numeric::ublas::norm_2(u);
       double normV = boost::numeric::ublas::norm_2(v);
 
-      if (IsCloseToZero(normU) ||
-          IsCloseToZero(normV))
+      if (LinearAlgebra::IsCloseToZero(normU) ||
+          LinearAlgebra::IsCloseToZero(normV))
       {
         return false;
       }
@@ -198,12 +76,12 @@
       // The angle must be zero, so the cosine must be almost equal to
       // cos(0) == 1 (or cos(180) == -1 if allowOppositeDirection == true)
 
-      if (IsCloseToZero(cosAngle - 1.0))
+      if (LinearAlgebra::IsCloseToZero(cosAngle - 1.0))
       {
         isOpposite = false;
         return true;
       }
-      else if (IsCloseToZero(fabs(cosAngle) - 1.0))
+      else if (LinearAlgebra::IsCloseToZero(fabs(cosAngle) - 1.0))
       {
         isOpposite = true;
         return true;
@@ -237,10 +115,10 @@
 
       // The direction of the line of intersection is orthogonal to the
       // normal of both planes
-      CrossProduct(direction, normal1, normal2);
+      LinearAlgebra::CrossProduct(direction, normal1, normal2);
 
       double norm = boost::numeric::ublas::norm_2(direction);
-      if (IsCloseToZero(norm))
+      if (LinearAlgebra::IsCloseToZero(norm))
       {
         // The two planes are parallel or coincident
         return false;
@@ -250,7 +128,7 @@
       double d2 = -boost::numeric::ublas::inner_prod(normal2, origin2);
       Vector tmp = d2 * normal1 - d1 * normal2;
 
-      CrossProduct(p, tmp, direction);
+      LinearAlgebra::CrossProduct(p, tmp, direction);
       p /= norm;
 
       return true;
@@ -312,19 +190,19 @@
 
       // Create the coordinates of the rectangle
       Vector x[4];
-      AssignVector(x[0], xmin, ymin, 1.0);
-      AssignVector(x[1], xmax, ymin, 1.0);
-      AssignVector(x[2], xmax, ymax, 1.0);
-      AssignVector(x[3], xmin, ymax, 1.0);
+      LinearAlgebra::AssignVector(x[0], xmin, ymin, 1.0);
+      LinearAlgebra::AssignVector(x[1], xmax, ymin, 1.0);
+      LinearAlgebra::AssignVector(x[2], xmax, ymax, 1.0);
+      LinearAlgebra::AssignVector(x[3], xmin, ymax, 1.0);
 
       // Move to homogoneous coordinates in 2D
       Vector p;
 
       {
         Vector a, b;
-        AssignVector(a, ax, ay, 1.0);
-        AssignVector(b, bx, by, 1.0);
-        CrossProduct(p, a, b);
+        LinearAlgebra::AssignVector(a, ax, ay, 1.0);
+        LinearAlgebra::AssignVector(b, bx, by, 1.0);
+        LinearAlgebra::CrossProduct(p, a, b);
       }
 
       uint8_t c = 0;
@@ -349,10 +227,10 @@
       else
       {
         Vector a, b, e;
-        CrossProduct(e, x[i], x[(i + 1) % 4]);
-        CrossProduct(a, p, e);
-        CrossProduct(e, x[j], x[(j + 1) % 4]);
-        CrossProduct(b, p, e);
+        LinearAlgebra::CrossProduct(e, x[i], x[(i + 1) % 4]);
+        LinearAlgebra::CrossProduct(a, p, e);
+        LinearAlgebra::CrossProduct(e, x[j], x[(j + 1) % 4]);
+        LinearAlgebra::CrossProduct(b, p, e);
 
         // Go back to non-homogeneous coordinates
         x1 = a[0] / a[2];
@@ -371,7 +249,7 @@
     {
       Vector v;
 
-      if (ParseVector(v, dicom, Orthanc::DICOM_TAG_PIXEL_SPACING))
+      if (LinearAlgebra::ParseVector(v, dicom, Orthanc::DICOM_TAG_PIXEL_SPACING))
       {
         if (v.size() != 2 ||
             v[0] <= 0 ||
@@ -513,51 +391,5 @@
         return true;
       }
     }
-
-
-    void FillMatrix(Matrix& target,
-                    size_t rows,
-                    size_t columns,
-                    const double values[])
-    {
-      target.resize(rows, columns);
-
-      size_t index = 0;
-
-      for (size_t y = 0; y < rows; y++)
-      {
-        for (size_t x = 0; x < columns; x++, index++)
-        {
-          target(y, x) = values[index];
-        }
-      }
-    }
-
-
-    void FillVector(Vector& target,
-                    size_t size,
-                    const double values[])
-    {
-      target.resize(size);
-
-      for (size_t i = 0; i < size; i++)
-      {
-        target[i] = values[i];
-      }
-    }
-
-
-    void Convert(Matrix& target,
-                 const Vector& source)
-    {
-      const size_t n = source.size();
-
-      target.resize(n, 1);
-
-      for (size_t i = 0; i < n; i++)
-      {
-        target(i, 0) = source[i];
-      }      
-    }
   }
 }

--- a/Framework/Toolbox/GeometryToolbox.h	Fri Feb 09 16:00:29 2018 +0100
+++ b/Framework/Toolbox/GeometryToolbox.h	Wed Feb 14 08:58:31 2018 +0100
@@ -21,67 +21,12 @@
 
 #pragma once
 
-// Patch for ublas in Boost 1.64.0
-// https://github.com/dealii/dealii/issues/4302
-#include <boost/version.hpp>
-#if BOOST_VERSION >= 106300  // or 64, need to check
-#  include <boost/serialization/array_wrapper.hpp>
-#endif
-
-#include <Core/DicomFormat/DicomMap.h>
-
-#include <boost/numeric/ublas/matrix.hpp>
-#include <boost/numeric/ublas/vector.hpp>
-#include <cassert>
+#include "LinearAlgebra.h"
 
 namespace OrthancStone
 {
-  typedef boost::numeric::ublas::matrix<double>   Matrix;
-  typedef boost::numeric::ublas::vector<double>   Vector;
-
   namespace GeometryToolbox
   {
-    void Print(const Vector& v);
-
-    void Print(const Matrix& m);
-
-    bool ParseVector(Vector& target,
-                     const std::string& s);
-
-    bool ParseVector(Vector& target,
-                     const Orthanc::DicomMap& dataset,
-                     const Orthanc::DicomTag& tag);
-
-    void AssignVector(Vector& v,
-                      double v1,
-                      double v2);
-
-    void AssignVector(Vector& v,
-                      double v1,
-                      double v2,
-                      double v3);
-
-    inline bool IsNear(double x,
-                       double y,
-                       double threshold)
-    {
-      return fabs(x - y) < threshold;
-    }
-
-    bool IsNear(double x,
-                double y);
-
-    inline bool IsCloseToZero(double x)
-    {
-      return IsNear(x, 0.0);
-    }
-
-    void NormalizeVector(Vector& u);
-
-    void CrossProduct(Vector& result,
-                      const Vector& u,
-                      const Vector& v);
-
     void ProjectPointOntoPlane(Vector& result,
                                const Vector& point,
                                const Vector& planeNormal,
@@ -142,18 +87,6 @@
                                const Vector& origin,
                                const Vector& direction);
 
-    void FillMatrix(Matrix& target,
-                    size_t rows,
-                    size_t columns,
-                    const double values[]);
-
-    void FillVector(Vector& target,
-                    size_t size,
-                    const double values[]);
-
-    void Convert(Matrix& target,
-                 const Vector& source);
-
     inline float ComputeBilinearInterpolationInternal(float x,
                                                       float y,
                                                       float f00,    // source(x, y)

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Framework/Toolbox/LinearAlgebra.cpp	Wed Feb 14 08:58:31 2018 +0100
@@ -0,0 +1,200 @@
+/**
+ * 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 "LinearAlgebra.h"
+
+#include <Core/Logging.h>
+#include <Core/OrthancException.h>
+#include <Core/Toolbox.h>
+
+#include <stdio.h>
+#include <boost/lexical_cast.hpp>
+
+namespace OrthancStone
+{
+  namespace LinearAlgebra
+  {
+    void Print(const Vector& v)
+    {
+      for (size_t i = 0; i < v.size(); i++)
+      {
+        printf("%g\n", v[i]);
+        //printf("%8.2f\n", v[i]);
+      }
+      printf("\n");
+    }
+
+
+    void Print(const Matrix& m)
+    {
+      for (size_t i = 0; i < m.size1(); i++)
+      {
+        for (size_t j = 0; j < m.size2(); j++)
+        {
+          printf("%g  ", m(i,j));
+          //printf("%8.2f  ", m(i,j));
+        }
+        printf("\n");        
+      }
+      printf("\n");        
+    }
+
+
+    bool ParseVector(Vector& target,
+                     const std::string& value)
+    {
+      std::vector<std::string> items;
+      Orthanc::Toolbox::TokenizeString(items, value, '\\');
+
+      target.resize(items.size());
+
+      for (size_t i = 0; i < items.size(); i++)
+      {
+        try
+        {
+          target[i] = boost::lexical_cast<double>(Orthanc::Toolbox::StripSpaces(items[i]));
+        }
+        catch (boost::bad_lexical_cast&)
+        {
+          target.clear();
+          return false;
+        }
+      }
+
+      return true;
+    }
+
+
+    bool ParseVector(Vector& target,
+                     const Orthanc::DicomMap& dataset,
+                     const Orthanc::DicomTag& tag)
+    {
+      std::string value;
+      return (dataset.CopyToString(value, tag, false) &&
+              ParseVector(target, value));
+    }
+
+
+    void AssignVector(Vector& v,
+                      double v1,
+                      double v2)
+    {
+      v.resize(2);
+      v[0] = v1;
+      v[1] = v2;
+    }
+
+
+    void AssignVector(Vector& v,
+                      double v1,
+                      double v2,
+                      double v3)
+    {
+      v.resize(3);
+      v[0] = v1;
+      v[1] = v2;
+      v[2] = v3;
+    }
+
+
+    bool IsNear(double x,
+                double y)
+    {
+      // As most input is read as single-precision numbers, we take the
+      // epsilon machine for float32 into consideration to compare numbers
+      return IsNear(x, y, 10.0 * std::numeric_limits<float>::epsilon());
+    }
+
+
+    void NormalizeVector(Vector& u)
+    {
+      double norm = boost::numeric::ublas::norm_2(u);
+      if (!IsCloseToZero(norm))
+      {
+        u = u / norm;
+      }
+    }
+
+
+    void CrossProduct(Vector& result,
+                      const Vector& u,
+                      const Vector& v)
+    {
+      if (u.size() != 3 ||
+          v.size() != 3)
+      {
+        throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
+      }
+
+      result.resize(3);
+
+      result[0] = u[1] * v[2] - u[2] * v[1];
+      result[1] = u[2] * v[0] - u[0] * v[2];
+      result[2] = u[0] * v[1] - u[1] * v[0];
+    }
+
+
+    void FillMatrix(Matrix& target,
+                    size_t rows,
+                    size_t columns,
+                    const double values[])
+    {
+      target.resize(rows, columns);
+
+      size_t index = 0;
+
+      for (size_t y = 0; y < rows; y++)
+      {
+        for (size_t x = 0; x < columns; x++, index++)
+        {
+          target(y, x) = values[index];
+        }
+      }
+    }
+
+
+    void FillVector(Vector& target,
+                    size_t size,
+                    const double values[])
+    {
+      target.resize(size);
+
+      for (size_t i = 0; i < size; i++)
+      {
+        target[i] = values[i];
+      }
+    }
+
+
+    void Convert(Matrix& target,
+                 const Vector& source)
+    {
+      const size_t n = source.size();
+
+      target.resize(n, 1);
+
+      for (size_t i = 0; i < n; i++)
+      {
+        target(i, 0) = source[i];
+      }      
+    }
+  }
+}

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Framework/Toolbox/LinearAlgebra.h	Wed Feb 14 08:58:31 2018 +0100
@@ -0,0 +1,96 @@
+/**
+ * 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
+
+// Patch for ublas in Boost 1.64.0
+// https://github.com/dealii/dealii/issues/4302
+#include <boost/version.hpp>
+#if BOOST_VERSION >= 106300  // or 64, need to check
+#  include <boost/serialization/array_wrapper.hpp>
+#endif
+
+#include <Core/DicomFormat/DicomMap.h>
+
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/vector.hpp>
+
+namespace OrthancStone
+{
+  typedef boost::numeric::ublas::matrix<double>   Matrix;
+  typedef boost::numeric::ublas::vector<double>   Vector;
+
+  namespace LinearAlgebra
+  {
+    void Print(const Vector& v);
+
+    void Print(const Matrix& m);
+
+    bool ParseVector(Vector& target,
+                     const std::string& s);
+
+    bool ParseVector(Vector& target,
+                     const Orthanc::DicomMap& dataset,
+                     const Orthanc::DicomTag& tag);
+
+    void AssignVector(Vector& v,
+                      double v1,
+                      double v2);
+
+    void AssignVector(Vector& v,
+                      double v1,
+                      double v2,
+                      double v3);
+
+    inline bool IsNear(double x,
+                       double y,
+                       double threshold)
+    {
+      return fabs(x - y) < threshold;
+    }
+
+    bool IsNear(double x,
+                double y);
+
+    inline bool IsCloseToZero(double x)
+    {
+      return IsNear(x, 0.0);
+    }
+
+    void NormalizeVector(Vector& u);
+
+    void CrossProduct(Vector& result,
+                      const Vector& u,
+                      const Vector& v);
+    
+    void FillMatrix(Matrix& target,
+                    size_t rows,
+                    size_t columns,
+                    const double values[]);
+
+    void FillVector(Vector& target,
+                    size_t size,
+                    const double values[]);
+
+    void Convert(Matrix& target,
+                 const Vector& source);
+  };
+}

--- a/Framework/Toolbox/OrientedBoundingBox.cpp	Fri Feb 09 16:00:29 2018 +0100
+++ b/Framework/Toolbox/OrientedBoundingBox.cpp	Wed Feb 14 08:58:31 2018 +0100
@@ -237,7 +237,7 @@
     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);
+    LinearAlgebra::AssignVector(target, x, y, z);
   }
 
 

--- a/Framework/Toolbox/ParallelSlices.cpp	Fri Feb 09 16:00:29 2018 +0100
+++ b/Framework/Toolbox/ParallelSlices.cpp	Wed Feb 14 08:58:31 2018 +0100
@@ -28,7 +28,7 @@
 {
   ParallelSlices::ParallelSlices()
   {
-    GeometryToolbox::AssignVector(normal_, 0, 0, 1);
+    LinearAlgebra::AssignVector(normal_, 0, 0, 1);
   }
 
 

--- a/Framework/Toolbox/Slice.cpp	Fri Feb 09 16:00:29 2018 +0100
+++ b/Framework/Toolbox/Slice.cpp	Wed Feb 14 08:58:31 2018 +0100
@@ -87,7 +87,7 @@
       return false;
     }
 
-    if (!GeometryToolbox::IsCloseToZero(offset0))
+    if (!LinearAlgebra::IsCloseToZero(offset0))
     {
       LOG(ERROR) << "Invalid syntax";
       return false;
@@ -304,9 +304,9 @@
     return (GeometryToolbox::IsParallelOrOpposite(opposite,
                                                   GetGeometry().GetNormal(),
                                                   plane.GetNormal()) &&
-            GeometryToolbox::IsNear(GetGeometry().ProjectAlongNormal(GetGeometry().GetOrigin()),
-                                    GetGeometry().ProjectAlongNormal(plane.GetOrigin()),
-                                    thickness_ / 2.0));
+            LinearAlgebra::IsNear(GetGeometry().ProjectAlongNormal(GetGeometry().GetOrigin()),
+                                  GetGeometry().ProjectAlongNormal(plane.GetOrigin()),
+                                  thickness_ / 2.0));
   }
 
   

--- a/Framework/Volumes/ImageBuffer3D.cpp	Fri Feb 09 16:00:29 2018 +0100
+++ b/Framework/Volumes/ImageBuffer3D.cpp	Wed Feb 14 08:58:31 2018 +0100
@@ -122,7 +122,7 @@
     computeRange_(computeRange),
     hasRange_(false)
   {
-    GeometryToolbox::AssignVector(voxelDimensions_, 1, 1, 1);
+    LinearAlgebra::AssignVector(voxelDimensions_, 1, 1, 1);
 
     LOG(INFO) << "Created an image of "
               << (GetEstimatedMemorySize() / (1024ll * 1024ll)) << "MB";
@@ -153,7 +153,7 @@
     }
 
     {
-      GeometryToolbox::AssignVector(voxelDimensions_, x, y, z);
+      LinearAlgebra::AssignVector(voxelDimensions_, x, y, z);
     }
   }
 
@@ -168,11 +168,11 @@
         break;
 
       case VolumeProjection_Coronal:
-        GeometryToolbox::AssignVector(result, voxelDimensions_[0], voxelDimensions_[2], voxelDimensions_[1]);
+        LinearAlgebra::AssignVector(result, voxelDimensions_[0], voxelDimensions_[2], voxelDimensions_[1]);
         break;
 
       case VolumeProjection_Sagittal:
-        GeometryToolbox::AssignVector(result, voxelDimensions_[1], voxelDimensions_[2], voxelDimensions_[0]);
+        LinearAlgebra::AssignVector(result, voxelDimensions_[1], voxelDimensions_[2], voxelDimensions_[0]);
         break;
 
       default:

--- a/Framework/Volumes/VolumeReslicer.cpp	Fri Feb 09 16:00:29 2018 +0100
+++ b/Framework/Volumes/VolumeReslicer.cpp	Wed Feb 14 08:58:31 2018 +0100
@@ -849,7 +849,7 @@
         slow.GetNearestCoordinates(qx, qy, qz);
 
         Vector d;
-        GeometryToolbox::AssignVector(d, px - qx, py - qy, pz - qz);
+        LinearAlgebra::AssignVector(d, px - qx, py - qy, pz - qz);
         double norm = boost::numeric::ublas::norm_2(d);
         if (norm > 0.0001)
         {

--- a/Framework/dev.h	Fri Feb 09 16:00:29 2018 +0100
+++ b/Framework/dev.h	Wed Feb 14 08:58:31 2018 +0100
@@ -86,8 +86,8 @@
         return false;
       }
 
-      if (!GeometryToolbox::IsNear(a.GetPixelSpacingX(), b.GetPixelSpacingX()) ||
-          !GeometryToolbox::IsNear(a.GetPixelSpacingY(), b.GetPixelSpacingY()))
+      if (!LinearAlgebra::IsNear(a.GetPixelSpacingX(), b.GetPixelSpacingX()) ||
+          !LinearAlgebra::IsNear(a.GetPixelSpacingY(), b.GetPixelSpacingY()))
       {
         LOG(ERROR) << "The pixel spacing of the slices change across the volume image";
         return false;
@@ -138,8 +138,8 @@
 
       for (size_t i = 1; i < loader.GetSliceCount(); i++)
       {
-        if (!GeometryToolbox::IsNear(spacingZ, GetDistance(loader.GetSlice(i - 1), loader.GetSlice(i)),
-                                     0.001 /* this is expressed in mm */))
+        if (!LinearAlgebra::IsNear(spacingZ, GetDistance(loader.GetSlice(i - 1), loader.GetSlice(i)),
+                                   0.001 /* this is expressed in mm */))
         {
           LOG(ERROR) << "The distance between successive slices is not constant in a volume image";
           SlicedVolumeBase::NotifyGeometryError();

--- a/Resources/CMake/OrthancStoneConfiguration.cmake	Fri Feb 09 16:00:29 2018 +0100
+++ b/Resources/CMake/OrthancStoneConfiguration.cmake	Wed Feb 14 08:58:31 2018 +0100
@@ -181,6 +181,7 @@
   ${ORTHANC_STONE_DIR}/Framework/Toolbox/DownloadStack.cpp
   ${ORTHANC_STONE_DIR}/Framework/Toolbox/Extent2D.cpp
   ${ORTHANC_STONE_DIR}/Framework/Toolbox/GeometryToolbox.cpp
+  ${ORTHANC_STONE_DIR}/Framework/Toolbox/LinearAlgebra.cpp
   ${ORTHANC_STONE_DIR}/Framework/Toolbox/MessagingToolbox.cpp
   ${ORTHANC_STONE_DIR}/Framework/Toolbox/OrientedBoundingBox.cpp
   ${ORTHANC_STONE_DIR}/Framework/Toolbox/OrthancSlicesLoader.cpp
@@ -195,9 +196,9 @@
   ${ORTHANC_STONE_DIR}/Framework/Viewport/WidgetViewport.cpp
   ${ORTHANC_STONE_DIR}/Framework/Volumes/ImageBuffer3D.cpp
   ${ORTHANC_STONE_DIR}/Framework/Volumes/SlicedVolumeBase.cpp
+  ${ORTHANC_STONE_DIR}/Framework/Volumes/StructureSetLoader.cpp
   ${ORTHANC_STONE_DIR}/Framework/Volumes/VolumeLoaderBase.cpp
   ${ORTHANC_STONE_DIR}/Framework/Volumes/VolumeReslicer.cpp
-  ${ORTHANC_STONE_DIR}/Framework/Volumes/StructureSetLoader.cpp
   ${ORTHANC_STONE_DIR}/Framework/Widgets/CairoWidget.cpp
   ${ORTHANC_STONE_DIR}/Framework/Widgets/EmptyWidget.cpp
   ${ORTHANC_STONE_DIR}/Framework/Widgets/LayerWidget.cpp