changeset 177:83200c4d07ca wasm

fix interpolation
author Sebastien Jodogne <s.jodogne@gmail.com>
date Fri, 09 Mar 2018 12:32:01 +0100
parents ab9c799f5de1
children 6dafcdec4b87
files Framework/Toolbox/GeometryToolbox.h Framework/Volumes/VolumeReslicer.cpp UnitTestsSources/UnitTestsMain.cpp
diffstat 3 files changed, 208 insertions(+), 190 deletions(-) [+]
line wrap: on
line diff
--- a/Framework/Toolbox/GeometryToolbox.h	Fri Mar 09 10:06:59 2018 +0100
+++ b/Framework/Toolbox/GeometryToolbox.h	Fri Mar 09 12:32:01 2018 +0100
@@ -95,44 +95,37 @@
                                const Vector& origin,
                                const Vector& direction);
 
-    inline float ComputeBilinearInterpolationInternal(float x,
-                                                      float y,
-                                                      float f00,    // source(x, y)
-                                                      float f01,    // source(x + 1, y)
-                                                      float f10,    // source(x, y + 1)
-                                                      float f11);   // source(x + 1, y + 1)
+    inline float ComputeBilinearInterpolationUnitSquare(float x,
+                                                        float y,
+                                                        float f00,    // source(0, 0)
+                                                        float f01,    // source(1, 0)
+                                                        float f10,    // source(0, 1)
+                                                        float f11);   // source(1, 1)
 
-    inline float ComputeBilinearInterpolation(float x,
-                                              float y,
-                                              float f00,    // source(x, y)
-                                              float f01,    // source(x + 1, y)
-                                              float f10,    // source(x, y + 1)
-                                              float f11);   // source(x + 1, y + 1)
-
-    inline float ComputeTrilinearInterpolation(float x,
-                                               float y,
-                                               float z,
-                                               float f000,   // source(x, y, z)
-                                               float f001,   // source(x + 1, y, z)
-                                               float f010,   // source(x, y + 1, z)
-                                               float f011,   // source(x + 1, y + 1, z)
-                                               float f100,   // source(x, y, z + 1)
-                                               float f101,   // source(x + 1, y, z + 1)
-                                               float f110,   // source(x, y + 1, z + 1)
-                                               float f111);  // source(x + 1, y + 1, z + 1)
+    inline float ComputeTrilinearInterpolationUnitSquare(float x,
+                                                         float y,
+                                                         float z,
+                                                         float f000,   // source(0, 0, 0)
+                                                         float f001,   // source(1, 0, 0)
+                                                         float f010,   // source(0, 1, 0)
+                                                         float f011,   // source(1, 1, 0)
+                                                         float f100,   // source(0, 0, 1)
+                                                         float f101,   // source(1, 0, 1)
+                                                         float f110,   // source(0, 1, 1)
+                                                         float f111);  // source(1, 1, 1)
   };
 }
 
 
-float OrthancStone::GeometryToolbox::ComputeBilinearInterpolationInternal(float x,
-                                                                          float y,
-                                                                          float f00,
-                                                                          float f01,
-                                                                          float f10,
-                                                                          float f11)
+float OrthancStone::GeometryToolbox::ComputeBilinearInterpolationUnitSquare(float x,
+                                                                            float y,
+                                                                            float f00,
+                                                                            float f01,
+                                                                            float f10,
+                                                                            float f11)
 {
-  // This function only works on fractional parts
-  assert(x >= 0 && y >= 0 && x < 1 && y < 1);
+  // This function only works within the unit square
+  assert(x >= 0 && y >= 0 && x <= 1 && y <= 1);
 
   // https://en.wikipedia.org/wiki/Bilinear_interpolation#Unit_square
   return (f00 * (1.0f - x) * (1.0f - y) +
@@ -142,42 +135,23 @@
 }
 
 
-float OrthancStone::GeometryToolbox::ComputeBilinearInterpolation(float x,
-                                                                  float y,
-                                                                  float f00,
-                                                                  float f01,
-                                                                  float f10,
-                                                                  float f11)
+float OrthancStone::GeometryToolbox::ComputeTrilinearInterpolationUnitSquare(float x,
+                                                                             float y,
+                                                                             float z,
+                                                                             float f000,
+                                                                             float f001,
+                                                                             float f010,
+                                                                             float f011,
+                                                                             float f100,
+                                                                             float f101,
+                                                                             float f110,
+                                                                             float f111)
 {
-  // Compute the fractional part of (x,y)
-  float xx = x - std::floor(x);
-  float yy = y - std::floor(y);
-
-  return ComputeBilinearInterpolationInternal(xx, yy, f00, f01, f10, f11);
-}
-
-
-float OrthancStone::GeometryToolbox::ComputeTrilinearInterpolation(float x,
-                                                                   float y,
-                                                                   float z,
-                                                                   float f000,
-                                                                   float f001,
-                                                                   float f010,
-                                                                   float f011,
-                                                                   float f100,
-                                                                   float f101,
-                                                                   float f110,
-                                                                   float f111)
-{
-  float xx = x - std::floor(x);
-  float yy = y - std::floor(y);
-  float zz = z - std::floor(z);
-
   // "In practice, a trilinear interpolation is identical to two
   // bilinear interpolation combined with a linear interpolation"
   // https://en.wikipedia.org/wiki/Trilinear_interpolation#Method
-  float a = ComputeBilinearInterpolationInternal(xx, yy, f000, f001, f010, f011);
-  float b = ComputeBilinearInterpolationInternal(xx, yy, f100, f101, f110, f111);
+  float a = ComputeBilinearInterpolationUnitSquare(x, y, f000, f001, f010, f011);
+  float b = ComputeBilinearInterpolationUnitSquare(x, y, f100, f101, f110, f111);
 
-  return (1.0f - zz) * a + zz * b;
+  return (1.0f - z) * a + z * b;
 }
--- a/Framework/Volumes/VolumeReslicer.cpp	Fri Mar 09 10:06:59 2018 +0100
+++ b/Framework/Volumes/VolumeReslicer.cpp	Fri Mar 09 12:32:01 2018 +0100
@@ -234,56 +234,84 @@
 
 
 
-    class VoxelReaderBase
+    class VoxelReaderBase : public boost::noncopyable
     {
-    protected:
-      const ImageBuffer3D&  source_;
-      unsigned int          sourceWidth_;
-      unsigned int          sourceHeight_;
-      unsigned int          sourceDepth_;
+    private:
+      const ImageBuffer3D&  image_;
+      unsigned int          width_;
+      unsigned int          height_;
+      unsigned int          depth_;
+      float                 widthFloat_;
+      float                 heightFloat_;
+      float                 depthFloat_;
 
     public:
-      VoxelReaderBase(const ImageBuffer3D& source) :
-        source_(source),
-        sourceWidth_(source.GetWidth()),
-        sourceHeight_(source.GetHeight()),
-        sourceDepth_(source.GetDepth())
+      VoxelReaderBase(const ImageBuffer3D& image) :
+        image_(image),
+        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()))
       {
       }
 
-      unsigned int GetSourceWidth() const
+      const ImageBuffer3D& GetImage() const
       {
-        return sourceWidth_;
+        return image_;
       }
 
-      unsigned int GetSourceHeight() const
+      const unsigned int GetImageWidth() const
       {
-        return sourceHeight_;
+        return width_;
       }
 
-      unsigned int GetSourceDepth() const
+      const unsigned int GetImageHeight() const
       {
-        return sourceDepth_;
+        return height_;
+      }
+
+      const unsigned int GetImageDepth() const
+      {
+        return depth_;
       }
 
-      bool GetNearestCoordinates(unsigned int& sourceX,
-                                 unsigned int& sourceY,
-                                 unsigned int& sourceZ,
-                                 float worldX,
-                                 float worldY,
-                                 float worldZ) const
+      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 (worldX >= 0 &&
-            worldY >= 0 &&
-            worldZ >= 0)
+        if (volumeX >= 0 &&
+            volumeY >= 0 &&
+            volumeZ >= 0)
         {
-          sourceX = static_cast<unsigned int>(worldX * static_cast<float>(sourceWidth_));
-          sourceY = static_cast<unsigned int>(worldY * static_cast<float>(sourceHeight_));
-          sourceZ = static_cast<unsigned int>(worldZ * static_cast<float>(sourceDepth_));
+          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));
 
-          return (sourceX < sourceWidth_ &&
-                  sourceY < sourceHeight_ &&
-                  sourceZ < sourceDepth_);
+          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
         {
@@ -305,32 +333,35 @@
     public:
       typedef typename PixelTraits<InputFormat>::PixelType   InputPixelType;
 
-      VoxelReader(const ImageBuffer3D& source) :
-        VoxelReaderBase(source)
+      VoxelReader(const ImageBuffer3D& image) :
+        VoxelReaderBase(image)
       {
       }
 
       ORTHANC_STONE_FORCE_INLINE
-      float GetFloatValue(float worldX,
-                          float worldY,
-                          float worldZ) const
+      float GetFloatValue(float volumeX,
+                          float volumeY,
+                          float volumeZ) const
       {
         InputPixelType value;
-        GetValue(value, worldX, worldY, worldZ);
+        GetValue(value, volumeX, volumeY, volumeZ);
         return static_cast<float>(value);
       }
 
       ORTHANC_STONE_FORCE_INLINE
       void GetValue(InputPixelType& target,
-                    float worldX,
-                    float worldY,
-                    float worldZ) const
+                    float volumeX,
+                    float volumeY,
+                    float volumeZ) const
       {
-        unsigned int sourceX, sourceY, sourceZ;
+        unsigned int imageX, imageY, imageZ;
+        float fractionalX, fractionalY, fractionalZ;  // unused
 
-        if (GetNearestCoordinates(sourceX, sourceY, sourceZ, worldX, worldY, worldZ))
+        if (GetNearestCoordinates(imageX, imageY, imageZ,
+                                  fractionalX, fractionalY, fractionalZ,
+                                  volumeX, volumeY, volumeZ))
         {
-          PixelTraits<InputFormat>::GetVoxel(target, source_, sourceX, sourceY, sourceZ);
+          PixelTraits<InputFormat>::GetVoxel(target, GetImage(), imageX, imageY, imageZ);
         }
         else
         {
@@ -348,8 +379,8 @@
       float outOfVolume_;
       
     public:
-      VoxelReader(const ImageBuffer3D& source) :
-        VoxelReaderBase(source)
+      VoxelReader(const ImageBuffer3D& image) :
+        VoxelReaderBase(image)
       {
         typename PixelTraits<InputFormat>::PixelType value;
         PixelTraits<InputFormat>::SetOutOfVolume(value);
@@ -357,37 +388,37 @@
       }
 
       void SampleVoxels(float& f00,
+                        float& f01,
                         float& f10,
-                        float& f01,
                         float& f11,
-                        unsigned int sourceX,
-                        unsigned int sourceY,
-                        unsigned int sourceZ) const
+                        unsigned int imageX,
+                        unsigned int imageY,
+                        unsigned int imageZ) const
       {
-        f00 = PixelTraits<InputFormat>::GetFloatVoxel(source_, sourceX, sourceY, sourceZ);
+        f00 = PixelTraits<InputFormat>::GetFloatVoxel(GetImage(), imageX, imageY, imageZ);
 
-        if (sourceX + 1 < sourceWidth_)
+        if (imageX + 1 < GetImageWidth())
         {
-          f01 = PixelTraits<InputFormat>::GetFloatVoxel(source_, sourceX + 1, sourceY, sourceZ);
+          f01 = PixelTraits<InputFormat>::GetFloatVoxel(GetImage(), imageX + 1, imageY, imageZ);
         }
         else
         {
           f01 = f00;
         }
 
-        if (sourceY + 1 < sourceHeight_)
+        if (imageY + 1 < GetImageWidth())
         {
-          f10 = PixelTraits<InputFormat>::GetFloatVoxel(source_, sourceX, sourceY + 1, sourceZ);
+          f10 = PixelTraits<InputFormat>::GetFloatVoxel(GetImage(), imageX, imageY + 1, imageZ);
         }
         else
         {
           f10 = f00;
         }
 
-        if (sourceX + 1 < sourceWidth_ &&
-            sourceY + 1 < sourceHeight_)
+        if (imageX + 1 < GetImageWidth() &&
+            imageY + 1 < GetImageHeight())
         {
-          f11 = PixelTraits<InputFormat>::GetFloatVoxel(source_, sourceX + 1, sourceY + 1, sourceZ);
+          f11 = PixelTraits<InputFormat>::GetFloatVoxel(GetImage(), imageX + 1, imageY + 1, imageZ);
         }
         else
         {
@@ -400,17 +431,21 @@
         return outOfVolume_;
       }
       
-      float GetFloatValue(float worldX,
-                          float worldY,
-                          float worldZ) const
+      float GetFloatValue(float volumeX,
+                          float volumeY,
+                          float volumeZ) const
       {
-        unsigned int sourceX, sourceY, sourceZ;
+        unsigned int imageX, imageY, imageZ;
+        float fractionalX, fractionalY, fractionalZ;
 
-        if (GetNearestCoordinates(sourceX, sourceY, sourceZ, worldX, worldY, worldZ))
+        if (GetNearestCoordinates(imageX, imageY, imageZ,
+                                  fractionalX, fractionalY, fractionalZ,
+                                  volumeX, volumeY, volumeZ))
         {
-          float f00, f10, f01, f11;
-          SampleVoxels(f00, f10, f01, f11, sourceX, sourceY, sourceZ);
-          return GeometryToolbox::ComputeBilinearInterpolation(worldX, worldY, f00, f01, f10, f11);
+          float f00, f01, f10, f11;
+          SampleVoxels(f00, f01, f10, f11, imageX, imageY, imageZ);
+          return GeometryToolbox::ComputeBilinearInterpolationUnitSquare
+            (fractionalX, fractionalY, f00, f01, f10, f11);
         }
         else
         {
@@ -421,42 +456,49 @@
 
 
     template <Orthanc::PixelFormat InputFormat>
-    class VoxelReader<InputFormat, ImageInterpolation_Trilinear>
+    class VoxelReader<InputFormat, ImageInterpolation_Trilinear> :
+      public VoxelReaderBase
     {
     private:
       typedef VoxelReader<InputFormat, ImageInterpolation_Bilinear>  Bilinear;
 
       Bilinear      bilinear_;
-      unsigned int  sourceDepth_;
+      unsigned int  imageDepth_;
 
     public:
-      VoxelReader(const ImageBuffer3D& source) :
-        bilinear_(source),
-        sourceDepth_(source.GetDepth())
+      VoxelReader(const ImageBuffer3D& image) :
+        VoxelReaderBase(image),
+        bilinear_(image),
+        imageDepth_(image.GetDepth())
       {
       }
 
-      float GetFloatValue(float worldX,
-                          float worldY,
-                          float worldZ) const
+      float GetFloatValue(float volumeX,
+                          float volumeY,
+                          float volumeZ) const
       {
-        unsigned int sourceX, sourceY, sourceZ;
+        unsigned int imageX, imageY, imageZ;
+        float fractionalX, fractionalY, fractionalZ;
 
-        if (bilinear_.GetNearestCoordinates(sourceX, sourceY, sourceZ, worldX, worldY, worldZ))
+        if (GetNearestCoordinates(imageX, imageY, imageZ,
+                                  fractionalX, fractionalY, fractionalZ,
+                                  volumeX, volumeY, volumeZ))
         {
-          float f000, f010, f001, f011, f100, f110, f101, f111;
-
-          bilinear_.SampleVoxels(f000, f010, f001, f011, sourceX, sourceY, sourceZ);
+          float f000, f001, f010, f011;
+          bilinear_.SampleVoxels(f000, f001, f010, f011, imageX, imageY, imageZ);
 
-          if (sourceZ + 1 < sourceDepth_)
+          if (imageZ + 1 < imageDepth_)
           {
-            bilinear_.SampleVoxels(f100, f110, f101, f111, sourceX, sourceY, sourceZ + 1);
-            return GeometryToolbox::ComputeTrilinearInterpolation
-              (worldX, worldY, worldZ, f000, f001, f010, f011, f100, f110, f101, f111);
+            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::ComputeBilinearInterpolation(worldX, worldY, f000, f001, f010, f011);
+            return GeometryToolbox::ComputeBilinearInterpolationUnitSquare
+              (fractionalX, fractionalY, f000, f001, f010, f011);
           }
         }
         else
@@ -491,13 +533,13 @@
       
       ORTHANC_STONE_FORCE_INLINE
       void Apply(typename PixelWriter::OutputPixelType* pixel,
-                 float worldX,
-                 float worldY,
-                 float worldZ)
+                 float volumeX,
+                 float volumeY,
+                 float volumeZ)
       {
-        typename VoxelReader::InputPixelType source;
-        reader_.GetValue(source, worldX, worldY, worldZ);
-        writer_.SetValue(pixel, source);
+        typename VoxelReader::InputPixelType image;
+        reader_.GetValue(image, volumeX, volumeY, volumeZ);
+        writer_.SetValue(pixel, image);
       }        
     };
 
@@ -520,11 +562,11 @@
       
       ORTHANC_STONE_FORCE_INLINE
       void Apply(typename PixelWriter::OutputPixelType* pixel,
-                 float worldX,
-                 float worldY,
-                 float worldZ)
+                 float volumeX,
+                 float volumeY,
+                 float volumeZ)
       {
-        writer_.SetFloatValue(pixel, reader_.GetFloatValue(worldX, worldY, worldZ));
+        writer_.SetFloatValue(pixel, reader_.GetFloatValue(volumeX, volumeY, volumeZ));
       }        
     };
 
@@ -551,11 +593,11 @@
       
       ORTHANC_STONE_FORCE_INLINE
       void Apply(typename PixelWriter::OutputPixelType* pixel,
-                 float worldX,
-                 float worldY,
-                 float worldZ)
+                 float volumeX,
+                 float volumeY,
+                 float volumeZ)
       {
-        writer_.SetFloatValue(pixel, scaling_ * reader_.GetFloatValue(worldX, worldY, worldZ) + offset_);
+        writer_.SetFloatValue(pixel, scaling_ * reader_.GetFloatValue(volumeX, volumeY, volumeZ) + offset_);
       }        
     };
 
@@ -609,9 +651,9 @@
       }
 
       ORTHANC_STONE_FORCE_INLINE
-      void GetNearestCoordinates(float& x,
-                                 float& y,
-                                 float& z) const
+      void GetVolumeCoordinates(float& x,
+                                float& y,
+                                float& z) const
       {
         x = position_[0];
         y = position_[1];
@@ -651,9 +693,9 @@
         x_++;
       }
 
-      void GetNearestCoordinates(float& x,
-                                 float& y,
-                                 float& z) const
+      void GetVolumeCoordinates(float& x,
+                                float& y,
+                                float& z) const
       {
         assert(x_ < slice_.GetWidth());
         
@@ -704,9 +746,9 @@
 
         for (unsigned int x = 0; x < outputWidth; x++, p++)
         {
-          float worldX, worldY, worldZ;
-          it.GetNearestCoordinates(worldX, worldY, worldZ);
-          shader.Apply(p, worldX, worldY, worldZ);
+          float volumeX, volumeY, volumeZ;
+          it.GetVolumeCoordinates(volumeX, volumeY, volumeZ);
+          shader.Apply(p, volumeX, volumeY, volumeZ);
           it.Next();
         }
       }
@@ -845,10 +887,10 @@
       for (unsigned int x = 0; x < slice_->GetWidth(); x++)
       {
         float px, py, pz;
-        fast.GetNearestCoordinates(px, py, pz);
+        fast.GetVolumeCoordinates(px, py, pz);
 
         float qx, qy, qz;
-        slow.GetNearestCoordinates(qx, qy, qz);
+        slow.GetVolumeCoordinates(qx, qy, qz);
 
         Vector d;
         LinearAlgebra::AssignVector(d, px - qx, py - qy, pz - qz);
--- a/UnitTestsSources/UnitTestsMain.cpp	Fri Mar 09 10:06:59 2018 +0100
+++ b/UnitTestsSources/UnitTestsMain.cpp	Fri Mar 09 12:32:01 2018 +0100
@@ -136,28 +136,30 @@
 
 TEST(GeometryToolbox, Interpolation)
 {
+  using namespace OrthancStone::GeometryToolbox;
+  
   // https://en.wikipedia.org/wiki/Bilinear_interpolation#Application_in_image_processing
-  ASSERT_FLOAT_EQ(146.1f, OrthancStone::GeometryToolbox::ComputeBilinearInterpolation
-                  (20.5f, 14.2f, 91, 210, 162, 95));
-
-  ASSERT_NEAR(146.1f, OrthancStone::GeometryToolbox::ComputeBilinearInterpolation
-              (-20.5f, 14.2f, 91, 210, 162, 95), 0.0001f);
+  ASSERT_FLOAT_EQ(146.1f, ComputeBilinearInterpolationUnitSquare(0.5f, 0.2f, 91, 210, 162, 95));
 
-  ASSERT_NEAR(146.1f, OrthancStone::GeometryToolbox::ComputeBilinearInterpolation
-              (-20.5f, -8.8f, 91, 210, 162, 95), 0.0001f);
+  ASSERT_FLOAT_EQ(91,  ComputeBilinearInterpolationUnitSquare(0, 0, 91, 210, 162, 95));
+  ASSERT_FLOAT_EQ(210, ComputeBilinearInterpolationUnitSquare(1, 0, 91, 210, 162, 95));
+  ASSERT_FLOAT_EQ(162, ComputeBilinearInterpolationUnitSquare(0, 1, 91, 210, 162, 95));
+  ASSERT_FLOAT_EQ(95,  ComputeBilinearInterpolationUnitSquare(1, 1, 91, 210, 162, 95));
 
-  ASSERT_NEAR(146.1f, OrthancStone::GeometryToolbox::ComputeBilinearInterpolation
-              (20.5f, -8.8f, 91, 210, 162, 95), 0.0001f);
-
-  ASSERT_FLOAT_EQ(123.35f, OrthancStone::GeometryToolbox::ComputeTrilinearInterpolation
-                  (20.5f, 15.2f, 5.7f,
+  ASSERT_FLOAT_EQ(123.35f, ComputeTrilinearInterpolationUnitSquare
+                  (0.5f, 0.2f, 0.7f,
                    91, 210, 162, 95,
                    51, 190, 80, 92));
 
-  ASSERT_FLOAT_EQ(123.35f, OrthancStone::GeometryToolbox::ComputeTrilinearInterpolation
-                  (20.5f, 15.2f, -6.3f,
-                   91, 210, 162, 95,
-                   51, 190, 80, 92));
+  ASSERT_FLOAT_EQ(ComputeBilinearInterpolationUnitSquare(0.5f, 0.2f, 91, 210, 162, 95),
+                  ComputeTrilinearInterpolationUnitSquare(0.5f, 0.2f, 0,
+                                                          91, 210, 162, 95,
+                                                          51, 190, 80, 92));
+
+  ASSERT_FLOAT_EQ(ComputeBilinearInterpolationUnitSquare(0.5f, 0.2f, 51, 190, 80, 92),
+                  ComputeTrilinearInterpolationUnitSquare(0.5f, 0.2f, 1,
+                                                          91, 210, 162, 95,
+                                                          51, 190, 80, 92));
 }