Mercurial > hg > orthanc-stone
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)); }