orthanc-stone: Framework/Volumes/VolumeReslicer.cpp comparison

comparison Framework/Volumes/VolumeReslicer.cpp @ 177:83200c4d07ca wasm

fix interpolation

author	Sebastien Jodogne <s.jodogne@gmail.com>
date	Fri, 09 Mar 2018 12:32:01 +0100
parents	316324f42848
children	db21c1810c89

comparison

equal deleted inserted replaced

-:ab9c799f5de1
+:83200c4d07ca
 }
 };
-class VoxelReaderBase
+class VoxelReaderBase : public boost::noncopyable
 {
-protected:
+private:
-const ImageBuffer3D&  source_;
+const ImageBuffer3D&  image_;
-unsigned int          sourceWidth_;
+unsigned int          width_;
-unsigned int          sourceHeight_;
+unsigned int          height_;
-unsigned int          sourceDepth_;
+unsigned int          depth_;
+float                 widthFloat_;
-public:
+float                 heightFloat_;
-VoxelReaderBase(const ImageBuffer3D& source) :
+float                 depthFloat_;
-source_(source),
-sourceWidth_(source.GetWidth()),
+public:
-sourceHeight_(source.GetHeight()),
+VoxelReaderBase(const ImageBuffer3D& image) :
-sourceDepth_(source.GetDepth())
+image_(image),
-{
+width_(image.GetWidth()),
-}
+height_(image.GetHeight()),
+depth_(image.GetDepth()),
-unsigned int GetSourceWidth() const
+widthFloat_(static_cast<float>(image.GetWidth())),
-{
+heightFloat_(static_cast<float>(image.GetHeight())),
-return sourceWidth_;
+depthFloat_(static_cast<float>(image.GetDepth()))
-}
+{
+}
-unsigned int GetSourceHeight() const
-{
+const ImageBuffer3D& GetImage() const
-return sourceHeight_;
+{
-}
+return image_;
+}
-unsigned int GetSourceDepth() const
-{
+const unsigned int GetImageWidth() const
-return sourceDepth_;
+{
-}
+return width_;
+}
-bool GetNearestCoordinates(unsigned int& sourceX,
-unsigned int& sourceY,
+const unsigned int GetImageHeight() const
-unsigned int& sourceZ,
+{
-float worldX,
+return height_;
-float worldY,
+}
-float worldZ) const
-{
+const unsigned int GetImageDepth() const
-if (worldX >= 0 &&
+{
-worldY >= 0 &&
+return depth_;
-worldZ >= 0)
+}
-{
-sourceX = static_cast<unsigned int>(worldX * static_cast<float>(sourceWidth_));
+bool GetNearestCoordinates(unsigned int& imageX,
-sourceY = static_cast<unsigned int>(worldY * static_cast<float>(sourceHeight_));
+unsigned int& imageY,
-sourceZ = static_cast<unsigned int>(worldZ * static_cast<float>(sourceDepth_));
+unsigned int& imageZ,
+float& fractionalX,
-return (sourceX < sourceWidth_ &&
+float& fractionalY,
-sourceY < sourceHeight_ &&
+float& fractionalZ,
-sourceZ < sourceDepth_);
+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;
 }
 public VoxelReaderBase
 {
 public:
 typedef typename PixelTraits<InputFormat>::PixelType   InputPixelType;
-VoxelReader(const ImageBuffer3D& source) :
+VoxelReader(const ImageBuffer3D& image) :
-VoxelReaderBase(source)
+VoxelReaderBase(image)
 {
 }
 ORTHANC_STONE_FORCE_INLINE
-float GetFloatValue(float worldX,
+float GetFloatValue(float volumeX,
-float worldY,
+float volumeY,
-float worldZ) const
+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 volumeX,
-float worldY,
+float volumeY,
-float worldZ) const
+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,
-PixelTraits<InputFormat>::GetVoxel(target, source_, sourceX, sourceY, sourceZ);
+fractionalX, fractionalY, fractionalZ,
+volumeX, volumeY, volumeZ))
+{
+PixelTraits<InputFormat>::GetVoxel(target, GetImage(), imageX, imageY, imageZ);
 }
 else
 {
 PixelTraits<InputFormat>::SetOutOfVolume(target);
 }
 {
 private:
 float outOfVolume_;
 public:
-VoxelReader(const ImageBuffer3D& source) :
+VoxelReader(const ImageBuffer3D& image) :
-VoxelReaderBase(source)
+VoxelReaderBase(image)
 {
 typename PixelTraits<InputFormat>::PixelType value;
 PixelTraits<InputFormat>::SetOutOfVolume(value);
 outOfVolume_ = static_cast<float>(value);
 }
 void SampleVoxels(float& f00,
+float& f01,
 float& f10,
-float& f01,
 float& f11,
-unsigned int sourceX,
+unsigned int imageX,
-unsigned int sourceY,
+unsigned int imageY,
-unsigned int sourceZ) const
+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_ &&
+if (imageX + 1 < GetImageWidth() &&
-sourceY + 1 < sourceHeight_)
+imageY + 1 < GetImageHeight())
 {
-f11 = PixelTraits<InputFormat>::GetFloatVoxel(source_, sourceX + 1, sourceY + 1, sourceZ);
+f11 = PixelTraits<InputFormat>::GetFloatVoxel(GetImage(), imageX + 1, imageY + 1, imageZ);
 }
 else
 {
 f11 = f00;
 }
 float GetOutOfVolume() const
 {
 return outOfVolume_;
 }
-float GetFloatValue(float worldX,
+float GetFloatValue(float volumeX,
-float worldY,
+float volumeY,
-float worldZ) const
+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,
-float f00, f10, f01, f11;
+fractionalX, fractionalY, fractionalZ,
-SampleVoxels(f00, f10, f01, f11, sourceX, sourceY, sourceZ);
+volumeX, volumeY, volumeZ))
-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
 {
 return outOfVolume_;
 }
 }
 };
 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) :
+VoxelReader(const ImageBuffer3D& image) :
-bilinear_(source),
+VoxelReaderBase(image),
-sourceDepth_(source.GetDepth())
+bilinear_(image),
-{
+imageDepth_(image.GetDepth())
-}
+{
+}
-float GetFloatValue(float worldX,
-float worldY,
+float GetFloatValue(float volumeX,
-float worldZ) const
+float volumeY,
-{
+float volumeZ) const
-unsigned int sourceX, sourceY, sourceZ;
+{
+unsigned int imageX, imageY, imageZ;
-if (bilinear_.GetNearestCoordinates(sourceX, sourceY, sourceZ, worldX, worldY, worldZ))
+float fractionalX, fractionalY, fractionalZ;
-{
-float f000, f010, f001, f011, f100, f110, f101, f111;
+if (GetNearestCoordinates(imageX, imageY, imageZ,
+fractionalX, fractionalY, fractionalZ,
-bilinear_.SampleVoxels(f000, f010, f001, f011, sourceX, sourceY, sourceZ);
+volumeX, volumeY, volumeZ))
+{
-if (sourceZ + 1 < sourceDepth_)
+float f000, f001, f010, f011;
+bilinear_.SampleVoxels(f000, f001, f010, f011, imageX, imageY, imageZ);
+if (imageZ + 1 < imageDepth_)
 {
-bilinear_.SampleVoxels(f100, f110, f101, f111, sourceX, sourceY, sourceZ + 1);
+float f100, f101, f110, f111;
-return GeometryToolbox::ComputeTrilinearInterpolation
+bilinear_.SampleVoxels(f100, f101, f110, f111, imageX, imageY, imageZ + 1);
-(worldX, worldY, worldZ, f000, f001, f010, f011, f100, f110, f101, f111);
+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
 {
 return bilinear_.GetOutOfVolume();
 {
 }
 ORTHANC_STONE_FORCE_INLINE
 void Apply(typename PixelWriter::OutputPixelType* pixel,
-float worldX,
+float volumeX,
-float worldY,
+float volumeY,
-float worldZ)
+float volumeZ)
 {
-typename VoxelReader::InputPixelType source;
+typename VoxelReader::InputPixelType image;
-reader_.GetValue(source, worldX, worldY, worldZ);
+reader_.GetValue(image, volumeX, volumeY, volumeZ);
-writer_.SetValue(pixel, source);
+writer_.SetValue(pixel, image);
 }
 };
 template <typename VoxelReader,
 {
 }
 ORTHANC_STONE_FORCE_INLINE
 void Apply(typename PixelWriter::OutputPixelType* pixel,
-float worldX,
+float volumeX,
-float worldY,
+float volumeY,
-float worldZ)
+float volumeZ)
 {
-writer_.SetFloatValue(pixel, reader_.GetFloatValue(worldX, worldY, worldZ));
+writer_.SetFloatValue(pixel, reader_.GetFloatValue(volumeX, volumeY, volumeZ));
 }
 };
 template <typename VoxelReader,
 {
 }
 ORTHANC_STONE_FORCE_INLINE
 void Apply(typename PixelWriter::OutputPixelType* pixel,
-float worldX,
+float volumeX,
-float worldY,
+float volumeY,
-float worldZ)
+float volumeZ)
 {
-writer_.SetFloatValue(pixel, scaling_ * reader_.GetFloatValue(worldX, worldY, worldZ) + offset_);
+writer_.SetFloatValue(pixel, scaling_ * reader_.GetFloatValue(volumeX, volumeY, volumeZ) + offset_);
 }
 };
 position_[1] += offset_[1];
 position_[2] += offset_[2];
 }
 ORTHANC_STONE_FORCE_INLINE
-void GetNearestCoordinates(float& x,
+void GetVolumeCoordinates(float& x,
 float& y,
 float& z) const
 {
 x = position_[0];
 y = position_[1];
 z = position_[2];
 }
 void Next()
 {
 x_++;
 }
-void GetNearestCoordinates(float& x,
+void GetVolumeCoordinates(float& x,
 float& y,
 float& z) const
 {
 assert(x_ < slice_.GetWidth());
 const double width = static_cast<double>(slice_.GetWidth());
 const double height = static_cast<double>(slice_.GetHeight());
 RowIterator it(slice, extent, plane, box, y);
 for (unsigned int x = 0; x < outputWidth; x++, p++)
 {
-float worldX, worldY, worldZ;
+float volumeX, volumeY, volumeZ;
-it.GetNearestCoordinates(worldX, worldY, worldZ);
+it.GetVolumeCoordinates(volumeX, volumeY, volumeZ);
-shader.Apply(p, worldX, worldY, worldZ);
+shader.Apply(p, volumeX, volumeY, volumeZ);
 it.Next();
 }
 }
 }
 SlowRowIterator slow(*slice_, extent_, plane, box, y);
 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);
 double norm = boost::numeric::ublas::norm_2(d);
 if (norm > 0.0001)

Mercurial > hg > orthanc-stone

comparison Framework/Volumes/VolumeReslicer.cpp @ 177:83200c4d07ca wasm