# HG changeset patch
# User Sebastien Jodogne <s.jodogne@gmail.com>
# Date 1521112321 -3600
# Node ID 2cbfb08f3a957687f382f8f9ff7893b15a035c60
# Parent  ff85568745579e060f11ccb4b260bb41248f3c52
ImageGeometry.cpp

diff -r ff8556874557 -r 2cbfb08f3a95 Framework/Toolbox/ImageGeometry.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Framework/Toolbox/ImageGeometry.cpp	Thu Mar 15 12:12:01 2018 +0100
@@ -0,0 +1,556 @@
+/**
+ * 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 "ImageGeometry.h"
+
+#include "Extent2D.h"
+#include "SubpixelReader.h"
+
+#include <Core/Images/ImageProcessing.h>
+#include <Core/Logging.h>
+
+
+namespace OrthancStone
+{
+  static void AddTransformedPoint(Extent2D& extent,
+                                  const Matrix& a,
+                                  double x,
+                                  double y)
+  {
+    assert(a.size1() == 3 &&
+           a.size2() == 3);
+
+    Vector p = LinearAlgebra::Product(a, LinearAlgebra::CreateVector(x, y, 1));
+
+    if (!LinearAlgebra::IsCloseToZero(p[2]))
+    {
+      extent.AddPoint(p[0] / p[2], p[1] / p[2]);
+    }
+  }
+  
+  
+  bool GetPerpectiveTransformExtent(unsigned int& x1,
+                                    unsigned int& y1,
+                                    unsigned int& x2,
+                                    unsigned int& y2,
+                                    const Matrix& a,
+                                    unsigned int sourceWidth,
+                                    unsigned int sourceHeight,
+                                    unsigned int targetWidth,
+                                    unsigned int targetHeight)
+  {
+    if (targetWidth == 0 ||
+        targetHeight == 0)
+    {
+      return false;
+    }
+    
+    Extent2D extent;
+    AddTransformedPoint(extent, a, 0, 0);
+    AddTransformedPoint(extent, a, sourceWidth, 0);
+    AddTransformedPoint(extent, a, 0, sourceHeight);
+    AddTransformedPoint(extent, a, sourceWidth, sourceHeight);
+
+    if (extent.IsEmpty())
+    {
+      return false;
+    }
+    else
+    {
+      int tmp;
+
+      tmp = std::floor(extent.GetX1());
+      if (tmp < 0)
+      {
+        x1 = 0;
+      }
+      else
+      {
+        x1 = static_cast<unsigned int>(tmp);
+      }
+
+      tmp = std::floor(extent.GetY1());
+      if (tmp < 0)
+      {
+        y1 = 0;
+      }
+      else
+      {
+        y1 = static_cast<unsigned int>(tmp);
+      }
+
+      tmp = std::ceil(extent.GetX2());
+      if (tmp < 0)
+      {
+        return false;
+      }
+      else if (static_cast<unsigned int>(tmp) >= targetWidth)
+      {
+        x2 = targetWidth - 1;
+      }
+      else
+      {
+        x2 = static_cast<unsigned int>(tmp);
+      }
+
+      tmp = std::ceil(extent.GetY2());
+      if (tmp < 0)
+      {
+        return false;
+      }
+      else if (static_cast<unsigned int>(tmp) >= targetHeight)
+      {
+        y2 = targetHeight - 1;
+      }
+      else
+      {
+        y2 = static_cast<unsigned int>(tmp);
+      }
+
+      return (x1 <= x2 &&
+              y1 <= y2);
+    }
+  }
+
+
+  template <typename Reader,
+            bool HasOffsetX,
+            bool HasOffsetY>
+  static void ApplyAffineTransformToRow(typename Reader::PixelType* p,
+                                        Reader& reader,
+                                        unsigned int x1,
+                                        unsigned int x2,
+                                        float positionX,
+                                        float positionY,
+                                        float offsetX,
+                                        float offsetY)
+  {
+    typename Reader::PixelType value;
+
+    for (unsigned int x = x1; x <= x2; x++, p++)
+    {     
+      if (reader.GetValue(value, positionX, positionY))
+      {
+        *p = value;
+      }
+      else
+      {
+        Reader::Traits::SetZero(*p);
+      }
+
+      if (HasOffsetX)
+      {
+        positionX += offsetX;
+      }
+
+      if (HasOffsetY)
+      {
+        positionY += offsetY;
+      }
+    }
+  }
+  
+
+  template <Orthanc::PixelFormat Format,
+            ImageInterpolation Interpolation>
+  static void ApplyAffineInternal(Orthanc::ImageAccessor& target,
+                                  const Orthanc::ImageAccessor& source,
+                                  const Matrix& a)
+  {
+    assert(target.GetFormat() == Format &&
+           source.GetFormat() == Format);
+    
+    typedef SubpixelReader<Format, Interpolation>  Reader;
+    typedef typename Reader::PixelType             PixelType;
+
+    if (Format == Orthanc::PixelFormat_RGB24)
+    {
+      Orthanc::ImageProcessing::Set(target, 0, 0, 0, 255);
+    }
+    else
+    {
+      Orthanc::ImageProcessing::Set(target, 0);
+    }
+
+    Matrix inva;
+    if (!LinearAlgebra::InvertMatrixUnsafe(inva, a))
+    {
+      // Singular matrix
+      return;
+    }
+
+    Reader reader(source);
+
+    unsigned int x1, y1, x2, y2;
+
+    if (GetPerpectiveTransformExtent(x1, y1, x2, y2, a,
+                                     source.GetWidth(), source.GetHeight(),
+                                     target.GetWidth(), target.GetHeight()))
+    {
+      const size_t targetPitch = target.GetPitch();
+      uint8_t *targetRow = reinterpret_cast<uint8_t*>(reinterpret_cast<PixelType*>(target.GetRow(y1)) + x1);
+
+      for (unsigned int y = y1; y <= y2; y++)
+      {
+        Vector start;
+        LinearAlgebra::AssignVector(start, static_cast<double>(x1) + 0.5,
+                                    static_cast<double>(y) + 0.5, 1);
+        start = boost::numeric::ublas::prod(inva, start);
+        assert(LinearAlgebra::IsNear(1.0, start(2)));
+
+        Vector offset;
+        LinearAlgebra::AssignVector(offset, static_cast<double>(x1) + 1.5,
+                                    static_cast<double>(y) + 0.5, 1);
+        offset = boost::numeric::ublas::prod(inva, offset) - start;
+        assert(LinearAlgebra::IsNear(0.0, offset(2)));
+
+        float startX = static_cast<float>(start[0]);
+        float startY = static_cast<float>(start[1]);
+        float offsetX = static_cast<float>(offset[0]);
+        float offsetY = static_cast<float>(offset[1]);
+
+        PixelType* pixel = reinterpret_cast<PixelType*>(targetRow);
+        if (LinearAlgebra::IsCloseToZero(offsetX))
+        {
+          ApplyAffineTransformToRow<Reader, false, true>
+            (pixel, reader, x1, x2, startX, startY, offsetX, offsetY);
+        }
+        else if (LinearAlgebra::IsCloseToZero(offsetY))
+        {
+          ApplyAffineTransformToRow<Reader, true, false>
+            (pixel, reader, x1, x2, startX, startY, offsetX, offsetY);
+        }
+        else
+        {
+          ApplyAffineTransformToRow<Reader, true, true>
+            (pixel, reader, x1, x2, startX, startY, offsetX, offsetY);
+        }
+
+        targetRow += targetPitch;
+      }
+    }    
+  }
+
+
+  void ApplyAffineTransform(Orthanc::ImageAccessor& target,
+                            const Orthanc::ImageAccessor& source,
+                            double a11,
+                            double a12,
+                            double b1,
+                            double a21,
+                            double a22,
+                            double b2,
+                            ImageInterpolation interpolation)
+  {
+    if (source.GetFormat() != target.GetFormat())
+    {
+      throw Orthanc::OrthancException(Orthanc::ErrorCode_IncompatibleImageFormat);
+    }
+
+    if (interpolation != ImageInterpolation_Nearest &&
+        interpolation != ImageInterpolation_Bilinear)
+    {
+      throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
+    }
+
+    Matrix a;
+    a.resize(3, 3);
+    a(0, 0) = a11;
+    a(0, 1) = a12;
+    a(0, 2) = b1;
+    a(1, 0) = a21;
+    a(1, 1) = a22;
+    a(1, 2) = b2;
+    a(2, 0) = 0;
+    a(2, 1) = 0;
+    a(2, 2) = 1;
+
+    switch (source.GetFormat())
+    {
+      case Orthanc::PixelFormat_Grayscale8:
+        switch (interpolation)
+        {
+          case ImageInterpolation_Nearest:
+            ApplyAffineInternal<Orthanc::PixelFormat_Grayscale8, 
+                                ImageInterpolation_Nearest>(target, source, a);
+            break;
+
+          case ImageInterpolation_Bilinear:
+            ApplyAffineInternal<Orthanc::PixelFormat_Grayscale8, 
+                                ImageInterpolation_Bilinear>(target, source, a);
+            break;
+
+          default:
+            throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
+        }
+        break;
+
+      case Orthanc::PixelFormat_Grayscale16:
+        switch (interpolation)
+        {
+          case ImageInterpolation_Nearest:
+            ApplyAffineInternal<Orthanc::PixelFormat_Grayscale16, 
+                                ImageInterpolation_Nearest>(target, source, a);
+            break;
+
+          case ImageInterpolation_Bilinear:
+            ApplyAffineInternal<Orthanc::PixelFormat_Grayscale16, 
+                                ImageInterpolation_Bilinear>(target, source, a);
+            break;
+
+          default:
+            throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
+        }
+        break;
+
+      case Orthanc::PixelFormat_SignedGrayscale16:
+        switch (interpolation)
+        {
+          case ImageInterpolation_Nearest:
+            ApplyAffineInternal<Orthanc::PixelFormat_SignedGrayscale16, 
+                                ImageInterpolation_Nearest>(target, source, a);
+            break;
+
+          case ImageInterpolation_Bilinear:
+            ApplyAffineInternal<Orthanc::PixelFormat_SignedGrayscale16, 
+                                ImageInterpolation_Bilinear>(target, source, a);
+            break;
+
+          default:
+            throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
+        }
+        break;
+
+      case Orthanc::PixelFormat_RGB24:
+        switch (interpolation)
+        {
+          case ImageInterpolation_Nearest:
+            ApplyAffineInternal<Orthanc::PixelFormat_RGB24, 
+                                ImageInterpolation_Nearest>(target, source, a);
+            break;
+
+          default:
+            throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
+        }
+        break;
+
+      default:
+        throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
+    }
+  }
+
+
+  template <Orthanc::PixelFormat Format,
+            ImageInterpolation Interpolation>
+  static void ApplyPerspectiveInternal(Orthanc::ImageAccessor& target,
+                                       const Orthanc::ImageAccessor& source,
+                                       const Matrix& a,
+                                       const Matrix& inva)
+  {
+    assert(target.GetFormat() == Format &&
+           source.GetFormat() == Format);
+    
+    typedef SubpixelReader<Format, Interpolation> Reader;
+    typedef typename Reader::PixelType            PixelType;
+
+    Reader reader(source);
+    unsigned int x1, y1, x2, y2;
+
+    const float floatWidth = source.GetWidth();
+    const float floatHeight = source.GetHeight();
+
+    if (GetPerpectiveTransformExtent(x1, y1, x2, y2, a,
+                                     source.GetWidth(), source.GetHeight(),
+                                     target.GetWidth(), target.GetHeight()))
+    {
+      const size_t targetPitch = target.GetPitch();
+      uint8_t *targetRow = reinterpret_cast<uint8_t*>(reinterpret_cast<PixelType*>(target.GetRow(y1)) + x1);
+
+      for (unsigned int y = y1; y <= y2; y++)
+      {
+        PixelType *p = reinterpret_cast<PixelType*>(targetRow);
+
+        for (unsigned int x = x1; x <= x2; x++)
+        {
+          Vector v;
+          LinearAlgebra::AssignVector(v, static_cast<double>(x) + 0.5, 
+                                      static_cast<double>(y) + 0.5, 1);
+
+          Vector vv = LinearAlgebra::Product(inva, v);
+
+          assert(!LinearAlgebra::IsCloseToZero(vv[2]));
+          const double w = 1.0 / vv[2];
+          const float sourceX = static_cast<float>(vv[0] * w);
+          const float sourceY = static_cast<float>(vv[1] * w);
+          
+          // Make sure no integer overflow will occur after truncation
+          // (the static_cast<unsigned int> could otherwise throw an
+          // exception in WebAssembly if strong perspective effects)
+          if (sourceX < floatWidth &&
+              sourceY < floatHeight)
+          { 
+            reader.GetValue(*p, sourceX, sourceY);
+          }
+          else
+          {
+            Reader::Traits::SetZero(*p);
+          }
+
+          p++;
+        }
+
+        targetRow += targetPitch;
+      }
+    }
+  }
+
+    
+  void ApplyPerspectiveTransform(Orthanc::ImageAccessor& target,
+                                 const Orthanc::ImageAccessor& source,
+                                 const Matrix& a,
+                                 ImageInterpolation interpolation)
+  {
+    if (source.GetFormat() != target.GetFormat())
+    {
+      throw Orthanc::OrthancException(Orthanc::ErrorCode_IncompatibleImageFormat);
+    }
+
+    if (a.size1() != 3 ||
+        a.size2() != 3)
+    {
+      throw Orthanc::OrthancException(Orthanc::ErrorCode_IncompatibleImageSize);
+    }
+
+    if (interpolation != ImageInterpolation_Nearest &&
+        interpolation != ImageInterpolation_Bilinear)
+    {
+      throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
+    }
+
+    // Check whether we are dealing with an affine transform
+    if (LinearAlgebra::IsCloseToZero(a(2, 0)) &&
+        LinearAlgebra::IsCloseToZero(a(2, 1)))
+    {
+      double w = a(2, 2);
+      if (LinearAlgebra::IsCloseToZero(w))
+      {
+        LOG(ERROR) << "Singular perspective matrix";
+        throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
+      }
+      else
+      {
+        ApplyAffineTransform(target, source, 
+                             a(0, 0) / w, a(0, 1) / w, a(0, 2) / w,
+                             a(1, 0) / w, a(1, 1) / w, a(1, 2) / w,
+                             interpolation);
+        return;
+      }
+    }
+
+    if (target.GetFormat() == Orthanc::PixelFormat_RGB24)
+    {
+      Orthanc::ImageProcessing::Set(target, 0, 0, 0, 255);
+    }
+    else
+    {
+      Orthanc::ImageProcessing::Set(target, 0);
+    }
+
+    Matrix inva;
+    if (!LinearAlgebra::InvertMatrixUnsafe(inva, a))
+    {
+      return;
+    }
+
+    switch (source.GetFormat())
+    {
+      case Orthanc::PixelFormat_Grayscale8:
+        switch (interpolation)
+        {
+          case ImageInterpolation_Nearest:
+            ApplyPerspectiveInternal<Orthanc::PixelFormat_Grayscale8, 
+                                     ImageInterpolation_Nearest>(target, source, a, inva);
+            break;
+
+          case ImageInterpolation_Bilinear:
+            ApplyPerspectiveInternal<Orthanc::PixelFormat_Grayscale8, 
+                                     ImageInterpolation_Bilinear>(target, source, a, inva);
+            break;
+
+          default:
+            throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
+        }
+        break;
+
+      case Orthanc::PixelFormat_Grayscale16:
+        switch (interpolation)
+        {
+          case ImageInterpolation_Nearest:
+            ApplyPerspectiveInternal<Orthanc::PixelFormat_Grayscale16, 
+                                     ImageInterpolation_Nearest>(target, source, a, inva);
+            break;
+
+          case ImageInterpolation_Bilinear:
+            ApplyPerspectiveInternal<Orthanc::PixelFormat_Grayscale16, 
+                                     ImageInterpolation_Bilinear>(target, source, a, inva);
+            break;
+
+          default:
+            throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
+        }
+        break;
+
+      case Orthanc::PixelFormat_SignedGrayscale16:
+        switch (interpolation)
+        {
+          case ImageInterpolation_Nearest:
+            ApplyPerspectiveInternal<Orthanc::PixelFormat_SignedGrayscale16, 
+                                     ImageInterpolation_Nearest>(target, source, a, inva);
+            break;
+
+          case ImageInterpolation_Bilinear:
+            ApplyPerspectiveInternal<Orthanc::PixelFormat_SignedGrayscale16, 
+                                     ImageInterpolation_Bilinear>(target, source, a, inva);
+            break;
+
+          default:
+            throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
+        }
+        break;
+
+      case Orthanc::PixelFormat_RGB24:
+        switch (interpolation)
+        {
+          case ImageInterpolation_Nearest:
+            ApplyPerspectiveInternal<Orthanc::PixelFormat_RGB24, 
+                                     ImageInterpolation_Nearest>(target, source, a, inva);
+            break;
+
+          default:
+            throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
+        }
+        break;
+
+      default:
+        throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
+    }
+  }
+}
diff -r ff8556874557 -r 2cbfb08f3a95 Framework/Toolbox/ImageGeometry.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Framework/Toolbox/ImageGeometry.h	Thu Mar 15 12:12:01 2018 +0100
@@ -0,0 +1,59 @@
+/**
+ * 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 "../Enumerations.h"
+#include "LinearAlgebra.h"
+
+#include <Core/Images/ImageAccessor.h>
+
+
+namespace OrthancStone
+{
+  // Returns the "useful" portion of the target image when applying a
+  // 3x3 perspective transform "a" (i.e. the bounding box where points
+  // of the source image are mapped to)
+  bool GetPerpectiveTransformExtent(unsigned int& x1,
+                                    unsigned int& y1,
+                                    unsigned int& x2,
+                                    unsigned int& y2,
+                                    const Matrix& a,
+                                    unsigned int sourceWidth,
+                                    unsigned int sourceHeight,
+                                    unsigned int targetWidth,
+                                    unsigned int targetHeight);
+
+  void ApplyAffineTransform(Orthanc::ImageAccessor& target,
+                            const Orthanc::ImageAccessor& source,
+                            double a11,
+                            double a12,
+                            double b1,
+                            double a21,
+                            double a22,
+                            double b2,
+                            ImageInterpolation interpolation);
+
+  void ApplyPerspectiveTransform(Orthanc::ImageAccessor& target,
+                                 const Orthanc::ImageAccessor& source,
+                                 const Matrix& a,
+                                 ImageInterpolation interpolation);
+}
diff -r ff8556874557 -r 2cbfb08f3a95 Framework/Toolbox/SubpixelReader.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Framework/Toolbox/SubpixelReader.h	Thu Mar 15 12:12:01 2018 +0100
@@ -0,0 +1,212 @@
+/**
+ * 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 "../Enumerations.h"
+#include "GeometryToolbox.h"
+
+#include <Core/Images/ImageTraits.h>
+
+#include <boost/noncopyable.hpp>
+#include <cmath>
+
+namespace OrthancStone
+{
+  namespace Internals
+  {
+    class SubpixelReaderBase : public boost::noncopyable
+    {
+    private:
+      const Orthanc::ImageAccessor& source_;
+      unsigned int                  width_;
+      unsigned int                  height_;
+
+    public:
+      SubpixelReaderBase(const Orthanc::ImageAccessor& source) :
+        source_(source),
+        width_(source.GetWidth()),
+        height_(source.GetHeight())
+      {
+      }
+
+      ORTHANC_FORCE_INLINE
+      const Orthanc::ImageAccessor& GetSource() const
+      {
+        return source_;
+      }
+      
+      ORTHANC_FORCE_INLINE
+      unsigned int GetWidth() const
+      {
+        return width_;
+      }
+      
+      ORTHANC_FORCE_INLINE
+      unsigned int GetHeight() const
+      {
+        return height_;
+      }
+    };
+  }
+
+    
+  template <Orthanc::PixelFormat Format,
+            ImageInterpolation Interpolation>
+  class SubpixelReader;
+
+    
+  template <Orthanc::PixelFormat Format>
+  class SubpixelReader<Format, ImageInterpolation_Nearest> : 
+    public Internals::SubpixelReaderBase
+  {
+  public:
+    typedef Orthanc::PixelTraits<Format>  Traits;
+    typedef typename Traits::PixelType    PixelType;
+
+    SubpixelReader(const Orthanc::ImageAccessor& source) :
+      SubpixelReaderBase(source)
+    {
+    }
+
+    inline bool GetValue(PixelType& target,
+                         float x,
+                         float y) const;
+  };
+    
+    
+  template <Orthanc::PixelFormat Format>
+  class SubpixelReader<Format, ImageInterpolation_Bilinear> : 
+    public Internals::SubpixelReaderBase
+  {
+  public:
+    typedef Orthanc::PixelTraits<Format>  Traits;
+    typedef typename Traits::PixelType    PixelType;
+
+    SubpixelReader(const Orthanc::ImageAccessor& source) :
+      SubpixelReaderBase(source)
+    {
+    }
+
+    inline bool GetValue(PixelType& target,
+                         float x,
+                         float y);
+  };
+
+
+
+  template <Orthanc::PixelFormat Format>
+  bool SubpixelReader<Format, ImageInterpolation_Nearest>::GetValue(PixelType& target,
+                                                                    float x,
+                                                                    float y) const
+  {
+    if (x < 0 ||
+        y < 0)
+    {
+      return false;
+    }
+    else
+    {
+      unsigned int ux = static_cast<unsigned int>(std::floor(x));
+      unsigned int uy = static_cast<unsigned int>(std::floor(y));
+
+      if (ux < GetWidth() &&
+          uy < GetHeight())
+      {
+        Orthanc::ImageTraits<Format>::GetPixel(target, GetSource(), ux, uy);
+        return true;
+      }
+      else
+      {
+        return false;
+      }
+    }
+  }
+
+
+
+  template <Orthanc::PixelFormat Format>
+  bool SubpixelReader<Format, ImageInterpolation_Bilinear>::GetValue(PixelType& target,
+                                                                     float x,
+                                                                     float y)
+  {
+    x -= 0.5f;
+    y -= 0.5f;
+        
+    if (x < 0 ||
+        y < 0)
+    {
+      return false;
+    }
+    else
+    {
+      unsigned int ux = static_cast<unsigned int>(std::floor(x));
+      unsigned int uy = static_cast<unsigned int>(std::floor(y));
+
+      float f00, f01, f10, f11;
+
+      if (ux < GetWidth() &&
+          uy < GetHeight())
+      {
+        f00 = Orthanc::ImageTraits<Format>::GetFloatPixel(GetSource(), ux, uy);
+      }
+      else
+      {
+        return false;
+      }
+
+      if (ux + 1 < GetWidth())
+      {
+        f01 = Orthanc::ImageTraits<Format>::GetFloatPixel(GetSource(), ux + 1, uy);
+      }
+      else
+      {
+        f01 = f00;
+      }
+
+      if (uy + 1 < GetHeight())
+      {
+        f10 = Orthanc::ImageTraits<Format>::GetFloatPixel(GetSource(), ux, uy + 1);
+      }
+      else
+      {
+        f10 = f00;
+      }
+
+      if (ux + 1 < GetWidth() &&
+          uy + 1 < GetHeight())
+      {
+        f11 = Orthanc::ImageTraits<Format>::GetFloatPixel(GetSource(), ux + 1, uy + 1);
+      }
+      else
+      {
+        f11 = f00;
+      }
+
+      float ax = x - static_cast<float>(ux);
+      float ay = y - static_cast<float>(uy);
+      float value = GeometryToolbox::ComputeBilinearInterpolationUnitSquare(ax, ay, f00, f01, f10, f11);
+          
+      Traits::FloatToPixel(target, value);
+      return true;
+    }
+  }
+}
diff -r ff8556874557 -r 2cbfb08f3a95 Resources/CMake/OrthancStoneConfiguration.cmake
--- a/Resources/CMake/OrthancStoneConfiguration.cmake	Wed Mar 14 17:48:02 2018 +0100
+++ b/Resources/CMake/OrthancStoneConfiguration.cmake	Thu Mar 15 12:12:01 2018 +0100
@@ -178,6 +178,7 @@
   ${ORTHANC_STONE_DIR}/Framework/Toolbox/Extent2D.cpp
   ${ORTHANC_STONE_DIR}/Framework/Toolbox/FiniteProjectiveCamera.cpp
   ${ORTHANC_STONE_DIR}/Framework/Toolbox/GeometryToolbox.cpp
+  ${ORTHANC_STONE_DIR}/Framework/Toolbox/ImageGeometry.cpp
   ${ORTHANC_STONE_DIR}/Framework/Toolbox/LinearAlgebra.cpp
   ${ORTHANC_STONE_DIR}/Framework/Toolbox/MessagingToolbox.cpp
   ${ORTHANC_STONE_DIR}/Framework/Toolbox/OrientedBoundingBox.cpp