--- a/Framework/Toolbox/SubpixelReader.h	Thu Mar 15 12:12:01 2018 +0100
+++ b/Framework/Toolbox/SubpixelReader.h	Thu Mar 15 15:19:24 2018 +0100
@@ -90,7 +90,12 @@
     inline bool GetValue(PixelType& target,
                          float x,
                          float y) const;
+
+    inline bool GetFloatValue(float& target,
+                              float x,
+                              float y) const;
   };
+
     
     
   template <Orthanc::PixelFormat Format>
@@ -106,9 +111,13 @@
     {
     }
 
+    inline bool GetFloatValue(float& target,
+                              float x,
+                              float y) const;
+
     inline bool GetValue(PixelType& target,
                          float x,
-                         float y);
+                         float y) const;
   };
 
 
@@ -144,9 +153,49 @@
 
 
   template <Orthanc::PixelFormat Format>
+  bool SubpixelReader<Format, ImageInterpolation_Nearest>::GetFloatValue(float& target,
+                                                                         float x,
+                                                                         float y) const
+  {
+    PixelType value;
+    
+    if (GetValue(value, x, y))
+    {
+      target = Traits::PixelToFloat(value);
+      return true;
+    }
+    else
+    {
+      return false;
+    }
+  }
+
+
+
+  template <Orthanc::PixelFormat Format>
   bool SubpixelReader<Format, ImageInterpolation_Bilinear>::GetValue(PixelType& target,
                                                                      float x,
-                                                                     float y)
+                                                                     float y) const
+  {
+    float value;
+
+    if (GetFloatValue(value, x, y))
+    {
+      Traits::FloatToPixel(target, value);
+      return true;
+    }
+    else
+    {
+      return false;
+    }
+  }
+
+
+
+  template <Orthanc::PixelFormat Format>
+  bool SubpixelReader<Format, ImageInterpolation_Bilinear>::GetFloatValue(float& target,
+                                                                          float x,
+                                                                          float y) const
   {
     x -= 0.5f;
     y -= 0.5f;
@@ -203,9 +252,8 @@
 
       float ax = x - static_cast<float>(ux);
       float ay = y - static_cast<float>(uy);
-      float value = GeometryToolbox::ComputeBilinearInterpolationUnitSquare(ax, ay, f00, f01, f10, f11);
+      target = GeometryToolbox::ComputeBilinearInterpolationUnitSquare(ax, ay, f00, f01, f10, f11);
           
-      Traits::FloatToPixel(target, value);
       return true;
     }
   }

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Framework/Toolbox/SubvoxelReader.h	Thu Mar 15 15:19:24 2018 +0100
@@ -0,0 +1,419 @@
+/**
+ * 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 "../Volumes/ImageBuffer3D.h"
+#include "GeometryToolbox.h"
+
+#include <Core/Images/ImageTraits.h>
+
+#include <boost/noncopyable.hpp>
+#include <cmath>
+
+namespace OrthancStone
+{
+  namespace Internals
+  {
+    class SubvoxelReaderBase : public boost::noncopyable
+    {
+    private:
+      const ImageBuffer3D& source_;
+      unsigned int         width_;
+      unsigned int         height_;
+      unsigned int         depth_;
+
+    public:
+      SubvoxelReaderBase(const ImageBuffer3D& source) :
+        source_(source),
+        width_(source.GetWidth()),
+        height_(source.GetHeight()),
+        depth_(source.GetDepth())
+      {
+      }
+
+      ORTHANC_FORCE_INLINE
+      const Orthanc::ImageAccessor& GetSource() const
+      {
+        return source_.GetInternalImage();
+      }
+      
+      ORTHANC_FORCE_INLINE
+      unsigned int GetWidth() const
+      {
+        return width_;
+      }
+      
+      ORTHANC_FORCE_INLINE
+      unsigned int GetHeight() const
+      {
+        return height_;
+      }
+      
+      ORTHANC_FORCE_INLINE
+      unsigned int GetDepth() const
+      {
+        return depth_;
+      }
+      
+      ORTHANC_FORCE_INLINE
+      unsigned int ComputeRow(unsigned int y,
+                              unsigned int z) const
+      {
+        return z * height_ + y;
+      }
+    };
+  }
+
+    
+  template <Orthanc::PixelFormat Format,
+            ImageInterpolation Interpolation>
+  class SubvoxelReader;
+
+    
+  template <Orthanc::PixelFormat Format>
+  class SubvoxelReader<Format, ImageInterpolation_Nearest> : 
+    public Internals::SubvoxelReaderBase
+  {
+  public:
+    typedef Orthanc::PixelTraits<Format>  Traits;
+    typedef typename Traits::PixelType    PixelType;
+
+    SubvoxelReader(const ImageBuffer3D& source) :
+      SubvoxelReaderBase(source)
+    {
+    }
+
+    inline bool GetValue(PixelType& target,
+                         float x,
+                         float y,
+                         float z) const;
+    
+    inline bool GetFloatValue(float& target,
+                              float x,
+                              float y,
+                              float z) const;
+  };
+    
+    
+  template <Orthanc::PixelFormat Format>
+  class SubvoxelReader<Format, ImageInterpolation_Bilinear> : 
+    public Internals::SubvoxelReaderBase
+  {
+  public:
+    typedef Orthanc::PixelTraits<Format>  Traits;
+    typedef typename Traits::PixelType    PixelType;
+
+    SubvoxelReader(const ImageBuffer3D& source) :
+      SubvoxelReaderBase(source)
+    {
+    }
+
+    inline bool Sample(float& f00,
+                       float& f01,
+                       float& f10,
+                       float& f11,
+                       unsigned int ux,
+                       unsigned int uy,
+                       unsigned int uz) const;
+
+    inline bool GetValue(PixelType& target,
+                         float x,
+                         float y,
+                         float z) const;
+
+    inline bool GetFloatValue(float& target,
+                              float x,
+                              float y,
+                              float z) const;
+  };
+    
+
+  template <Orthanc::PixelFormat Format>
+  class SubvoxelReader<Format, ImageInterpolation_Trilinear> : 
+    public Internals::SubvoxelReaderBase
+  {
+  private:
+    SubvoxelReader<Format, ImageInterpolation_Bilinear>   bilinear_;
+
+  public:
+    typedef Orthanc::PixelTraits<Format>  Traits;
+    typedef typename Traits::PixelType    PixelType;
+
+    SubvoxelReader(const ImageBuffer3D& source) :
+      SubvoxelReaderBase(source),
+      bilinear_(source)
+    {
+    }
+
+    inline bool GetValue(PixelType& target,
+                         float x,
+                         float y,
+                         float z) const;
+
+    inline bool GetFloatValue(float& target,
+                              float x,
+                              float y,
+                              float z) const;
+  };
+
+
+
+  template <Orthanc::PixelFormat Format>
+  bool SubvoxelReader<Format, ImageInterpolation_Nearest>::GetValue(PixelType& target,
+                                                                    float x,
+                                                                    float y,
+                                                                    float z) const
+  {
+    if (x < 0 ||
+        y < 0 ||
+        z < 0)
+    {
+      return false;
+    }
+    else
+    {
+      unsigned int ux = static_cast<unsigned int>(std::floor(x));
+      unsigned int uy = static_cast<unsigned int>(std::floor(y));
+      unsigned int uz = static_cast<unsigned int>(std::floor(z));
+
+      if (ux < GetWidth() &&
+          uy < GetHeight() &&
+          uz < GetDepth())
+      {
+        Orthanc::ImageTraits<Format>::GetPixel(target, GetSource(), ux, ComputeRow(uy, uz));
+        return true;
+      }
+      else
+      {
+        return false;
+      }
+    }
+  }
+
+
+  template <Orthanc::PixelFormat Format>
+  bool SubvoxelReader<Format, ImageInterpolation_Nearest>::GetFloatValue(float& target,
+                                                                         float x,
+                                                                         float y,
+                                                                         float z) const
+  {
+    PixelType value;
+    
+    if (GetValue(value, x, y, z))
+    {
+      target = Traits::PixelToFloat(value);
+      return true;
+    }
+    else
+    {
+      return false;
+    }
+  }
+
+
+
+  template <Orthanc::PixelFormat Format>
+  bool SubvoxelReader<Format, ImageInterpolation_Bilinear>::Sample(float& f00,
+                                                                   float& f01,
+                                                                   float& f10,
+                                                                   float& f11,
+                                                                   unsigned int ux,
+                                                                   unsigned int uy,
+                                                                   unsigned int uz) const
+  {
+    if (ux < GetWidth() &&
+        uy < GetHeight() &&
+        uz < GetDepth())
+    {
+      f00 = Orthanc::ImageTraits<Format>::GetFloatPixel(GetSource(), ux, ComputeRow(uy, uz));
+    }
+    else
+    {
+      // Pixel is out of the volume
+      return false;
+    }
+
+    if (ux + 1 < GetWidth())
+    {
+      f01 = Orthanc::ImageTraits<Format>::GetFloatPixel(GetSource(), ux + 1, ComputeRow(uy, uz));
+    }
+    else
+    {
+      f01 = f00;
+    }
+
+    if (uy + 1 < GetHeight())
+    {
+      f10 = Orthanc::ImageTraits<Format>::GetFloatPixel(GetSource(), ux, ComputeRow(uy + 1, uz));
+    }
+    else
+    {
+      f10 = f00;
+    }
+
+    if (ux + 1 < GetWidth() &&
+        uy + 1 < GetHeight())
+    {
+      f11 = Orthanc::ImageTraits<Format>::GetFloatPixel(GetSource(), ux + 1, ComputeRow(uy + 1, uz));
+    }
+    else
+    {
+      f11 = f00;
+    }
+
+    return true;
+  }
+
+
+
+  template <Orthanc::PixelFormat Format>
+  bool SubvoxelReader<Format, ImageInterpolation_Bilinear>::GetFloatValue(float& target,
+                                                                          float x,
+                                                                          float y,
+                                                                          float z) const
+  {
+    x -= 0.5f;
+    y -= 0.5f;
+        
+    if (x < 0 ||
+        y < 0 ||
+        z < 0)
+    {
+      return false;
+    }
+    else
+    {
+      unsigned int ux = static_cast<unsigned int>(std::floor(x));
+      unsigned int uy = static_cast<unsigned int>(std::floor(y));
+      unsigned int uz = static_cast<unsigned int>(std::floor(z));
+
+      float f00, f01, f10, f11;
+      if (Sample(f00, f01, f10, f11, ux, uy, uz))
+      {
+        float ax = x - static_cast<float>(ux);
+        float ay = y - static_cast<float>(uy);
+
+        target = GeometryToolbox::ComputeBilinearInterpolationUnitSquare(ax, ay, f00, f01, f10, f11);
+        return true;
+      }
+      else
+      {
+        return false;
+      }
+    }
+  }
+
+
+
+  template <Orthanc::PixelFormat Format>
+  bool SubvoxelReader<Format, ImageInterpolation_Bilinear>::GetValue(PixelType& target,
+                                                                     float x,
+                                                                     float y,
+                                                                     float z) const
+  {
+    float value;
+
+    if (GetFloatValue(value, x, y, z))
+    {
+      Traits::FloatToPixel(target, value);
+      return true;
+    }
+    else
+    {
+      return false;
+    }
+  }
+
+
+
+  template <Orthanc::PixelFormat Format>
+  bool SubvoxelReader<Format, ImageInterpolation_Trilinear>::GetFloatValue(float& target,
+                                                                           float x,
+                                                                           float y,
+                                                                           float z) const
+  {
+    x -= 0.5f;
+    y -= 0.5f;
+    z -= 0.5f;
+        
+    if (x < 0 ||
+        y < 0 ||
+        z < 0)
+    {
+      return false;
+    }
+    else
+    {
+      unsigned int ux = static_cast<unsigned int>(std::floor(x));
+      unsigned int uy = static_cast<unsigned int>(std::floor(y));
+      unsigned int uz = static_cast<unsigned int>(std::floor(z));
+
+      float f000, f001, f010, f011;
+      if (bilinear_.Sample(f000, f001, f010, f011, ux, uy, uz))
+      {
+        const float ax = x - static_cast<float>(ux);
+        const float ay = y - static_cast<float>(uy);
+
+        float f100, f101, f110, f111;
+
+        if (bilinear_.Sample(f100, f101, f110, f111, ux, uy, uz + 1))
+        {
+          const float az = z - static_cast<float>(uz);
+          target = GeometryToolbox::ComputeTrilinearInterpolationUnitSquare
+            (ax, ay, az, f000, f001, f010, f011, f100, f101, f110, f111);
+        }
+        else
+        {
+          target = GeometryToolbox::ComputeBilinearInterpolationUnitSquare
+            (ax, ay, f000, f001, f010, f011);
+        }
+
+        return true;
+      }
+      else
+      {
+        return false;
+      }
+    }
+  }
+
+
+
+  template <Orthanc::PixelFormat Format>
+  bool SubvoxelReader<Format, ImageInterpolation_Trilinear>::GetValue(PixelType& target,
+                                                                      float x,
+                                                                      float y,
+                                                                      float z) const
+  {
+    float value;
+
+    if (GetFloatValue(value, x, y, z))
+    {
+      Traits::FloatToPixel(target, value);
+      return true;
+    }
+    else
+    {
+      return false;
+    }
+  }
+}

--- a/Framework/Volumes/VolumeReslicer.cpp	Thu Mar 15 12:12:01 2018 +0100
+++ b/Framework/Volumes/VolumeReslicer.cpp	Thu Mar 15 15:19:24 2018 +0100
@@ -1,6 +1,7 @@
 #include "VolumeReslicer.h"
 
 #include "../Toolbox/GeometryToolbox.h"
+#include "../Toolbox/SubvoxelReader.h"
 
 #include <Core/Images/ImageTraits.h>
 #include <Core/Logging.h>
@@ -116,286 +117,6 @@
     };
 
 
-
-    class VoxelReaderBase : public boost::noncopyable
-    {
-    private:
-      const Orthanc::ImageAccessor&  image_;
-      unsigned int                   width_;
-      unsigned int                   height_;
-      unsigned int                   depth_;
-      float                          widthFloat_;
-      float                          heightFloat_;
-      float                          depthFloat_;
-
-    public:
-      VoxelReaderBase(const ImageBuffer3D& image) :
-        image_(image.GetInternalImage()),
-        width_(image.GetWidth()),
-        height_(image.GetHeight()),
-        depth_(image.GetDepth()),
-        widthFloat_(static_cast<float>(image.GetWidth())),
-        heightFloat_(static_cast<float>(image.GetHeight())),
-        depthFloat_(static_cast<float>(image.GetDepth()))
-      {
-      }
-
-      const Orthanc::ImageAccessor& GetImage() const
-      {
-        return image_;
-      }
-
-      const unsigned int GetImageWidth() const
-      {
-        return width_;
-      }
-
-      const unsigned int GetImageHeight() const
-      {
-        return height_;
-      }
-
-      const unsigned int GetImageDepth() const
-      {
-        return depth_;
-      }
-
-      bool GetNearestCoordinates(unsigned int& imageX,
-                                 unsigned int& imageY,
-                                 unsigned int& imageZ,
-                                 float& fractionalX,
-                                 float& fractionalY,
-                                 float& fractionalZ,
-                                 float volumeX,
-                                 float volumeY,
-                                 float volumeZ) const
-      {
-        if (volumeX >= 0 &&
-            volumeY >= 0 &&
-            volumeZ >= 0)
-        {
-          const float x = volumeX * widthFloat_;
-          const float y = volumeY * heightFloat_;
-          const float z = volumeZ * depthFloat_;
-          
-          imageX = static_cast<unsigned int>(std::floor(x));
-          imageY = static_cast<unsigned int>(std::floor(y));
-          imageZ = static_cast<unsigned int>(std::floor(z));
-
-          if (imageX < width_ &&
-              imageY < height_ &&
-              imageZ < depth_)
-          {
-            fractionalX = x - static_cast<float>(imageX);
-            fractionalY = y - static_cast<float>(imageY);
-            fractionalZ = z - static_cast<float>(imageZ);
-            return true;
-          }
-          else
-          {
-            return false;
-          }
-        }
-        else
-        {
-          return false;
-        }
-      }
-    };
-
-
-    template <Orthanc::PixelFormat InputFormat,
-              ImageInterpolation Interpolation>
-    class VoxelReader;
-
-    
-    template <Orthanc::PixelFormat InputFormat>
-    class VoxelReader<InputFormat, ImageInterpolation_Nearest> :
-      public VoxelReaderBase
-    {
-    public:
-      typedef typename Orthanc::PixelTraits<InputFormat>::PixelType   InputPixelType;
-
-      VoxelReader(const ImageBuffer3D& image) :
-        VoxelReaderBase(image)
-      {
-      }
-
-      ORTHANC_FORCE_INLINE
-      float GetFloatValue(float volumeX,
-                          float volumeY,
-                          float volumeZ) const
-      {
-        InputPixelType value;
-        GetValue(value, volumeX, volumeY, volumeZ);
-        return static_cast<float>(value);
-      }
-
-      ORTHANC_FORCE_INLINE
-      void GetValue(InputPixelType& target,
-                    float volumeX,
-                    float volumeY,
-                    float volumeZ) const
-      {
-        unsigned int imageX, imageY, imageZ;
-        float fractionalX, fractionalY, fractionalZ;  // unused
-
-        if (GetNearestCoordinates(imageX, imageY, imageZ,
-                                  fractionalX, fractionalY, fractionalZ,
-                                  volumeX, volumeY, volumeZ))
-        {
-          Orthanc::ImageTraits<InputFormat>::GetPixel(target, GetImage(), imageX, 
-                                                      imageY + imageZ * GetImageHeight());
-        }
-        else
-        {
-          target = std::numeric_limits<InputPixelType>::min();
-        }
-      }
-    };
-    
-    
-    template <Orthanc::PixelFormat InputFormat>
-    class VoxelReader<InputFormat, ImageInterpolation_Bilinear> :
-      public VoxelReaderBase
-    {
-    private:
-      float outOfVolume_;
-      
-    public:
-      VoxelReader(const ImageBuffer3D& image) :
-        VoxelReaderBase(image)
-      {
-        typedef typename Orthanc::PixelTraits<InputFormat>::PixelType Pixel;
-        outOfVolume_ = static_cast<float>(std::numeric_limits<Pixel>::min());
-      }
-
-      void SampleVoxels(float& f00,
-                        float& f01,
-                        float& f10,
-                        float& f11,
-                        unsigned int imageX,
-                        unsigned int imageY,
-                        unsigned int imageZ) const
-      {
-        f00 = Orthanc::ImageTraits<InputFormat>::GetFloatPixel
-          (GetImage(), imageX, imageY + imageZ * GetImageHeight());
-
-        if (imageX + 1 < GetImageWidth())
-        {
-          f01 = Orthanc::ImageTraits<InputFormat>::GetFloatPixel
-            (GetImage(), imageX + 1, imageY + imageZ * GetImageHeight());
-        }
-        else
-        {
-          f01 = f00;
-        }
-
-        if (imageY + 1 < GetImageWidth())
-        {
-          f10 = Orthanc::ImageTraits<InputFormat>::GetFloatPixel
-            (GetImage(), imageX, imageY + 1 + imageZ * GetImageHeight());
-        }
-        else
-        {
-          f10 = f00;
-        }
-
-        if (imageX + 1 < GetImageWidth() &&
-            imageY + 1 < GetImageHeight())
-        {
-          f11 = Orthanc::ImageTraits<InputFormat>::GetFloatPixel
-            (GetImage(), imageX + 1, imageY + 1 + imageZ * GetImageHeight());
-        }
-        else
-        {
-          f11 = f00;
-        }
-      }
-
-      float GetOutOfVolume() const
-      {
-        return outOfVolume_;
-      }
-      
-      float GetFloatValue(float volumeX,
-                          float volumeY,
-                          float volumeZ) const
-      {
-        unsigned int imageX, imageY, imageZ;
-        float fractionalX, fractionalY, fractionalZ;
-
-        if (GetNearestCoordinates(imageX, imageY, imageZ,
-                                  fractionalX, fractionalY, fractionalZ,
-                                  volumeX, volumeY, volumeZ))
-        {
-          float f00, f01, f10, f11;
-          SampleVoxels(f00, f01, f10, f11, imageX, imageY, imageZ);
-          return GeometryToolbox::ComputeBilinearInterpolationUnitSquare
-            (fractionalX, fractionalY, f00, f01, f10, f11);
-        }
-        else
-        {
-          return outOfVolume_;
-        }
-      }
-    };
-
-
-    template <Orthanc::PixelFormat InputFormat>
-    class VoxelReader<InputFormat, ImageInterpolation_Trilinear> :
-      public VoxelReaderBase
-    {
-    private:
-      typedef VoxelReader<InputFormat, ImageInterpolation_Bilinear>  Bilinear;
-
-      Bilinear      bilinear_;
-      unsigned int  imageDepth_;
-
-    public:
-      VoxelReader(const ImageBuffer3D& image) :
-        VoxelReaderBase(image),
-        bilinear_(image),
-        imageDepth_(image.GetDepth())
-      {
-      }
-
-      float GetFloatValue(float volumeX,
-                          float volumeY,
-                          float volumeZ) const
-      {
-        unsigned int imageX, imageY, imageZ;
-        float fractionalX, fractionalY, fractionalZ;
-
-        if (GetNearestCoordinates(imageX, imageY, imageZ,
-                                  fractionalX, fractionalY, fractionalZ,
-                                  volumeX, volumeY, volumeZ))
-        {
-          float f000, f001, f010, f011;
-          bilinear_.SampleVoxels(f000, f001, f010, f011, imageX, imageY, imageZ);
-
-          if (imageZ + 1 < imageDepth_)
-          {
-            float f100, f101, f110, f111;
-            bilinear_.SampleVoxels(f100, f101, f110, f111, imageX, imageY, imageZ + 1);
-            return GeometryToolbox::ComputeTrilinearInterpolationUnitSquare
-              (fractionalX, fractionalY, fractionalZ,
-               f000, f001, f010, f011, f100, f101, f110, f111);
-          }
-          else
-          {
-            return GeometryToolbox::ComputeBilinearInterpolationUnitSquare
-              (fractionalX, fractionalY, f000, f001, f010, f011);
-          }
-        }
-        else
-        {
-          return bilinear_.GetOutOfVolume();
-        }
-      }
-    };
-
-
     template <typename VoxelReader,
               typename PixelWriter,
               TransferFunction Function>
@@ -424,9 +145,14 @@
                  float volumeY,
                  float volumeZ)
       {
-        typename VoxelReader::InputPixelType image;
-        reader_.GetValue(image, volumeX, volumeY, volumeZ);
-        writer_.SetValue(pixel, image);
+        typename VoxelReader::PixelType value;
+
+        if (!reader_.GetValue(value, volumeX, volumeY, volumeZ))
+        {
+          VoxelReader::Traits::SetMinValue(value);
+        }
+
+        writer_.SetValue(pixel, value);
       }        
     };
 
@@ -438,12 +164,14 @@
     private:
       VoxelReader  reader_;
       PixelWriter  writer_;
+      float        outOfVolume_;
       
     public:
       PixelShader(const ImageBuffer3D& image,
                   float /* scaling */,
                   float /* offset */) :
-        reader_(image)
+        reader_(image),
+        outOfVolume_(static_cast<float>(std::numeric_limits<typename VoxelReader::PixelType>::min()))
       {
       }
       
@@ -453,8 +181,15 @@
                  float volumeY,
                  float volumeZ)
       {
-        writer_.SetFloatValue(pixel, reader_.GetFloatValue(volumeX, volumeY, volumeZ));
-      }        
+        float value;
+
+        if (!reader_.GetFloatValue(value, volumeX, volumeY, volumeZ))
+        {
+          value = outOfVolume_;
+        }
+
+        writer_.SetFloatValue(pixel, value);
+      }
     };
 
     
@@ -467,6 +202,7 @@
       PixelWriter  writer_;
       float        scaling_;
       float        offset_;
+      float        outOfVolume_;
       
     public:
       PixelShader(const ImageBuffer3D& image,
@@ -474,7 +210,8 @@
                   float offset) :
         reader_(image),
         scaling_(scaling),
-        offset_(offset)
+        offset_(offset),
+        outOfVolume_(static_cast<float>(std::numeric_limits<typename VoxelReader::PixelType>::min()))
       {
       }
       
@@ -484,7 +221,18 @@
                  float volumeY,
                  float volumeZ)
       {
-        writer_.SetFloatValue(pixel, scaling_ * reader_.GetFloatValue(volumeX, volumeY, volumeZ) + offset_);
+        float value;
+
+        if (reader_.GetFloatValue(value, volumeX, volumeY, volumeZ))
+        {
+          value = scaling_ * value + offset_;
+        }
+        else
+        {
+          value = outOfVolume_;
+        }
+
+        writer_.SetFloatValue(pixel, value);
       }        
     };
 
@@ -616,13 +364,17 @@
                              float scaling,
                              float offset)
     {
-      typedef VoxelReader<InputFormat, Interpolation>   Reader;
-      typedef PixelWriter<InputFormat, OutputFormat>    Writer;
-      typedef PixelShader<Reader, Writer, Function>     Shader;
+      typedef SubvoxelReader<InputFormat, Interpolation>   Reader;
+      typedef PixelWriter<InputFormat, OutputFormat>       Writer;
+      typedef PixelShader<Reader, Writer, Function>        Shader;
 
       const unsigned int outputWidth = slice.GetWidth();
       const unsigned int outputHeight = slice.GetHeight();
 
+      const float sourceWidth = static_cast<float>(source.GetWidth());
+      const float sourceHeight = static_cast<float>(source.GetHeight());
+      const float sourceDepth = static_cast<float>(source.GetDepth());
+
       Shader shader(source, scaling, offset);
 
       for (unsigned int y = 0; y < outputHeight; y++)
@@ -636,7 +388,11 @@
         {
           float volumeX, volumeY, volumeZ;
           it.GetVolumeCoordinates(volumeX, volumeY, volumeZ);
-          shader.Apply(p, volumeX, volumeY, volumeZ);
+
+          shader.Apply(p, 
+                       volumeX * sourceWidth, 
+                       volumeY * sourceHeight, 
+                       volumeZ * sourceDepth);
           it.Next();
         }
       }