Mercurial > hg > orthanc-stone
diff OrthancStone/Sources/Volumes/VolumeReslicer.cpp @ 1512:244ad1e4e76a
reorganization of folders
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Tue, 07 Jul 2020 16:21:02 +0200 |
parents | Framework/Volumes/VolumeReslicer.cpp@30deba7bc8e2 |
children | 4fb8fdf03314 |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/OrthancStone/Sources/Volumes/VolumeReslicer.cpp Tue Jul 07 16:21:02 2020 +0200 @@ -0,0 +1,839 @@ +/** + * Stone of Orthanc + * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics + * Department, University Hospital of Liege, Belgium + * Copyright (C) 2017-2020 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 "VolumeReslicer.h" + +#include "../Toolbox/GeometryToolbox.h" +#include "../Toolbox/SubvoxelReader.h" + +#include <Images/ImageTraits.h> +#include <Logging.h> +#include <OrthancException.h> + +#include <boost/math/special_functions/round.hpp> + +namespace OrthancStone +{ + // Anonymous namespace to avoid clashes between compilation modules + namespace + { + enum TransferFunction + { + TransferFunction_Copy, + TransferFunction_Float, + TransferFunction_Linear + }; + + + template <Orthanc::PixelFormat InputFormat, + Orthanc::PixelFormat OutputFormat, + ImageInterpolation Interpolation, + TransferFunction Function> + class PixelShader; + + + template <Orthanc::PixelFormat Format> + class PixelShader<Format, + Format, + ImageInterpolation_Nearest, + TransferFunction_Copy> + { + private: + typedef SubvoxelReader<Format, ImageInterpolation_Nearest> VoxelReader; + typedef Orthanc::PixelTraits<Format> PixelWriter; + + VoxelReader reader_; + + public: + PixelShader(const ImageBuffer3D& image, + float /* scaling */, + float /* offset */) : + reader_(image) + { + } + + ORTHANC_FORCE_INLINE + void Apply(typename PixelWriter::PixelType* pixel, + float volumeX, + float volumeY, + float volumeZ) + { + typename VoxelReader::PixelType value; + + if (!reader_.GetValue(value, volumeX, volumeY, volumeZ)) + { + VoxelReader::Traits::SetMinValue(value); + } + + *pixel = value; + } + }; + + + template <Orthanc::PixelFormat InputFormat, + Orthanc::PixelFormat OutputFormat> + class PixelShader<InputFormat, + OutputFormat, + ImageInterpolation_Nearest, + TransferFunction_Copy> + { + private: + typedef SubvoxelReader<InputFormat, ImageInterpolation_Nearest> VoxelReader; + typedef Orthanc::PixelTraits<OutputFormat> PixelWriter; + + VoxelReader reader_; + + public: + PixelShader(const ImageBuffer3D& image, + float /* scaling */, + float /* offset */) : + reader_(image) + { + } + + ORTHANC_FORCE_INLINE + void Apply(typename PixelWriter::PixelType* pixel, + float volumeX, + float volumeY, + float volumeZ) + { + typename VoxelReader::PixelType value; + + if (!reader_.GetValue(value, volumeX, volumeY, volumeZ)) + { + VoxelReader::Traits::SetMinValue(value); + } + + PixelWriter::FloatToPixel(*pixel, VoxelReader::Traits::PixelToFloat(value)); + } + }; + + + template <Orthanc::PixelFormat InputFormat, + Orthanc::PixelFormat OutputFormat, + ImageInterpolation Interpolation> + class PixelShader<InputFormat, + OutputFormat, + Interpolation, + TransferFunction_Float> + { + private: + typedef SubvoxelReader<InputFormat, Interpolation> VoxelReader; + typedef Orthanc::PixelTraits<OutputFormat> PixelWriter; + + VoxelReader reader_; + float outOfVolume_; + + public: + PixelShader(const ImageBuffer3D& image, + float /* scaling */, + float /* offset */) : + reader_(image), + outOfVolume_(static_cast<float>(std::numeric_limits<typename VoxelReader::PixelType>::min())) + { + } + + ORTHANC_FORCE_INLINE + void Apply(typename PixelWriter::PixelType* pixel, + float volumeX, + float volumeY, + float volumeZ) + { + float value; + + if (!reader_.GetFloatValue(value, volumeX, volumeY, volumeZ)) + { + value = outOfVolume_; + } + + PixelWriter::FloatToPixel(*pixel, value); + } + }; + + + template <Orthanc::PixelFormat InputFormat, + Orthanc::PixelFormat OutputFormat, + ImageInterpolation Interpolation> + class PixelShader<InputFormat, + OutputFormat, + Interpolation, + TransferFunction_Linear> + { + private: + typedef SubvoxelReader<InputFormat, Interpolation> VoxelReader; + typedef Orthanc::PixelTraits<OutputFormat> PixelWriter; + + VoxelReader reader_; + float scaling_; + float offset_; + float outOfVolume_; + + public: + PixelShader(const ImageBuffer3D& image, + float scaling, + float offset) : + reader_(image), + scaling_(scaling), + offset_(offset), + outOfVolume_(static_cast<float>(std::numeric_limits<typename VoxelReader::PixelType>::min())) + { + } + + ORTHANC_FORCE_INLINE + void Apply(typename PixelWriter::PixelType* pixel, + float volumeX, + float volumeY, + float volumeZ) + { + float value; + + if (reader_.GetFloatValue(value, volumeX, volumeY, volumeZ)) + { + value = scaling_ * value + offset_; + } + else + { + value = outOfVolume_; + } + + PixelWriter::FloatToPixel(*pixel, value); + } + }; + + + + class FastRowIterator : public boost::noncopyable + { + private: + float position_[3]; + float offset_[3]; + + public: + FastRowIterator(const Orthanc::ImageAccessor& slice, + const Extent2D& extent, + const CoordinateSystem3D& plane, + const OrientedVolumeBoundingBox& box, + unsigned int y) + { + const double width = static_cast<double>(slice.GetWidth()); + const double height = static_cast<double>(slice.GetHeight()); + assert(y < height); + + Vector q1 = plane.MapSliceToWorldCoordinates + (extent.GetX1() + extent.GetWidth() * static_cast<double>(0) / static_cast<double>(width + 1), + extent.GetY1() + extent.GetHeight() * static_cast<double>(y) / static_cast<double>(height + 1)); + + Vector q2 = plane.MapSliceToWorldCoordinates + (extent.GetX1() + extent.GetWidth() * static_cast<double>(width - 1) / static_cast<double>(width + 1), + extent.GetY1() + extent.GetHeight() * static_cast<double>(y) / static_cast<double>(height + 1)); + + Vector r1, r2; + box.ToInternalCoordinates(r1, q1); + box.ToInternalCoordinates(r2, q2); + + position_[0] = static_cast<float>(r1[0]); + position_[1] = static_cast<float>(r1[1]); + position_[2] = static_cast<float>(r1[2]); + + Vector tmp = (r2 - r1) / static_cast<double>(width - 1); + offset_[0] = static_cast<float>(tmp[0]); + offset_[1] = static_cast<float>(tmp[1]); + offset_[2] = static_cast<float>(tmp[2]); + } + + ORTHANC_FORCE_INLINE + void Next() + { + position_[0] += offset_[0]; + position_[1] += offset_[1]; + position_[2] += offset_[2]; + } + + ORTHANC_FORCE_INLINE + void GetVolumeCoordinates(float& x, + float& y, + float& z) const + { + x = position_[0]; + y = position_[1]; + z = position_[2]; + } + }; + + + class SlowRowIterator : public boost::noncopyable + { + private: + const Orthanc::ImageAccessor& slice_; + const Extent2D& extent_; + const CoordinateSystem3D& plane_; + const OrientedVolumeBoundingBox& box_; + unsigned int x_; + unsigned int y_; + + public: + SlowRowIterator(const Orthanc::ImageAccessor& slice, + const Extent2D& extent, + const CoordinateSystem3D& plane, + const OrientedVolumeBoundingBox& box, + unsigned int y) : + slice_(slice), + extent_(extent), + plane_(plane), + box_(box), + x_(0), + y_(y) + { + assert(y_ < slice_.GetHeight()); + } + + void Next() + { + 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()); + + Vector q = plane_.MapSliceToWorldCoordinates + (extent_.GetX1() + extent_.GetWidth() * static_cast<double>(x_) / (width + 1.0), + extent_.GetY1() + extent_.GetHeight() * static_cast<double>(y_) / (height + 1.0)); + + Vector r; + box_.ToInternalCoordinates(r, q); + + x = static_cast<float>(r[0]); + y = static_cast<float>(r[1]); + z = static_cast<float>(r[2]); + } + }; + + + template <typename RowIterator, + Orthanc::PixelFormat InputFormat, + Orthanc::PixelFormat OutputFormat, + ImageInterpolation Interpolation, + TransferFunction Function> + static void ProcessImage(Orthanc::ImageAccessor& slice, + const Extent2D& extent, + const ImageBuffer3D& source, + const CoordinateSystem3D& plane, + const OrientedVolumeBoundingBox& box, + float scaling, + float offset) + { + typedef PixelShader<InputFormat, OutputFormat, Interpolation, 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++) + { + typedef typename Orthanc::ImageTraits<OutputFormat>::PixelType PixelType; + PixelType* p = reinterpret_cast<PixelType*>(slice.GetRow(y)); + + RowIterator it(slice, extent, plane, box, y); + + for (unsigned int x = 0; x < outputWidth; x++, p++) + { + float volumeX, volumeY, volumeZ; + it.GetVolumeCoordinates(volumeX, volumeY, volumeZ); + + shader.Apply(p, + volumeX * sourceWidth, + volumeY * sourceHeight, + volumeZ * sourceDepth); + it.Next(); + } + } + } + + + template <typename RowIterator, + Orthanc::PixelFormat InputFormat, + Orthanc::PixelFormat OutputFormat> + static void ProcessImage(Orthanc::ImageAccessor& slice, + const Extent2D& extent, + const ImageBuffer3D& source, + const CoordinateSystem3D& plane, + const OrientedVolumeBoundingBox& box, + ImageInterpolation interpolation, + bool hasLinearFunction, + float scaling, + float offset) + { + if (hasLinearFunction) + { + switch (interpolation) + { + case ImageInterpolation_Nearest: + ProcessImage<RowIterator, InputFormat, OutputFormat, + ImageInterpolation_Nearest, TransferFunction_Linear> + (slice, extent, source, plane, box, scaling, offset); + break; + + case ImageInterpolation_Bilinear: + ProcessImage<RowIterator, InputFormat, OutputFormat, + ImageInterpolation_Bilinear, TransferFunction_Linear> + (slice, extent, source, plane, box, scaling, offset); + break; + + case ImageInterpolation_Trilinear: + ProcessImage<RowIterator, InputFormat, OutputFormat, + ImageInterpolation_Trilinear, TransferFunction_Linear> + (slice, extent, source, plane, box, scaling, offset); + break; + + default: + throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); + } + } + else + { + switch (interpolation) + { + case ImageInterpolation_Nearest: + ProcessImage<RowIterator, InputFormat, OutputFormat, + ImageInterpolation_Nearest, TransferFunction_Copy> + (slice, extent, source, plane, box, 0, 0); + break; + + case ImageInterpolation_Bilinear: + ProcessImage<RowIterator, InputFormat, OutputFormat, + ImageInterpolation_Bilinear, TransferFunction_Float> + (slice, extent, source, plane, box, 0, 0); + break; + + case ImageInterpolation_Trilinear: + ProcessImage<RowIterator, InputFormat, OutputFormat, + ImageInterpolation_Trilinear, TransferFunction_Float> + (slice, extent, source, plane, box, 0, 0); + break; + + default: + throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); + } + } + } + + + template <typename RowIterator> + static void ProcessImage(Orthanc::ImageAccessor& slice, + const Extent2D& extent, + const ImageBuffer3D& source, + const CoordinateSystem3D& plane, + const OrientedVolumeBoundingBox& box, + ImageInterpolation interpolation, + bool hasLinearFunction, + float scaling, + float offset) + { + if (source.GetFormat() == Orthanc::PixelFormat_Grayscale16 && + slice.GetFormat() == Orthanc::PixelFormat_Grayscale8) + { + ProcessImage<RowIterator, + Orthanc::PixelFormat_Grayscale16, + Orthanc::PixelFormat_Grayscale8> + (slice, extent, source, plane, box, interpolation, hasLinearFunction, scaling, offset); + } + else if (source.GetFormat() == Orthanc::PixelFormat_Grayscale16 && + slice.GetFormat() == Orthanc::PixelFormat_Grayscale16) + { + ProcessImage<RowIterator, + Orthanc::PixelFormat_Grayscale16, + Orthanc::PixelFormat_Grayscale16> + (slice, extent, source, plane, box, interpolation, hasLinearFunction, scaling, offset); + } + else if (source.GetFormat() == Orthanc::PixelFormat_SignedGrayscale16 && + slice.GetFormat() == Orthanc::PixelFormat_BGRA32) + { + ProcessImage<RowIterator, + Orthanc::PixelFormat_SignedGrayscale16, + Orthanc::PixelFormat_BGRA32> + (slice, extent, source, plane, box, interpolation, hasLinearFunction, scaling, offset); + } + else if (source.GetFormat() == Orthanc::PixelFormat_Grayscale16 && + slice.GetFormat() == Orthanc::PixelFormat_BGRA32) + { + ProcessImage<RowIterator, + Orthanc::PixelFormat_Grayscale16, + Orthanc::PixelFormat_BGRA32> + (slice, extent, source, plane, box, interpolation, hasLinearFunction, scaling, offset); + } + else + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); + } + } + } + + + + void VolumeReslicer::CheckIterators(const ImageBuffer3D& source, + const CoordinateSystem3D& plane, + const OrientedVolumeBoundingBox& box) const + { + for (unsigned int y = 0; y < slice_->GetHeight(); y++) + { + FastRowIterator fast(*slice_, extent_, plane, box, y); + SlowRowIterator slow(*slice_, extent_, plane, box, y); + + for (unsigned int x = 0; x < slice_->GetWidth(); x++) + { + float px, py, pz; + fast.GetVolumeCoordinates(px, py, pz); + + float 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) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); + } + + fast.Next(); + slow.Next(); + } + } + } + + + void VolumeReslicer::Reset() + { + success_ = false; + extent_.Reset(); + slice_.reset(NULL); + } + + + float VolumeReslicer::GetMinOutputValue() const + { + switch (outputFormat_) + { + case Orthanc::PixelFormat_Grayscale8: + case Orthanc::PixelFormat_Grayscale16: + case Orthanc::PixelFormat_BGRA32: + return 0.0f; + break; + + default: + throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); + } + } + + + float VolumeReslicer::GetMaxOutputValue() const + { + switch (outputFormat_) + { + case Orthanc::PixelFormat_Grayscale8: + case Orthanc::PixelFormat_BGRA32: + return static_cast<float>(std::numeric_limits<uint8_t>::max()); + break; + + case Orthanc::PixelFormat_Grayscale16: + return static_cast<float>(std::numeric_limits<uint16_t>::max()); + break; + + default: + throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); + } + } + + + VolumeReslicer::VolumeReslicer() : + outputFormat_(Orthanc::PixelFormat_Grayscale8), + interpolation_(ImageInterpolation_Nearest), + fastMode_(true), + success_(false) + { + ResetLinearFunction(); + } + + + void VolumeReslicer::GetLinearFunction(float& scaling, + float& offset) const + { + if (hasLinearFunction_) + { + scaling = scaling_; + offset = offset_; + } + else + { + scaling = 1.0f; + offset = 0.0f; + } + } + + + void VolumeReslicer::ResetLinearFunction() + { + Reset(); + hasLinearFunction_ = false; + scaling_ = 1.0f; + offset_ = 0.0f; + } + + + void VolumeReslicer::SetLinearFunction(float scaling, + float offset) + { + Reset(); + hasLinearFunction_ = true; + scaling_ = scaling; + offset_ = offset; + } + + + void VolumeReslicer::SetWindow(float low, + float high) + { + //printf("Range in pixel values: %f->%f\n", low, high); + float scaling = (GetMaxOutputValue() - GetMinOutputValue()) / (high - low); + float offset = GetMinOutputValue() - scaling * low; + + SetLinearFunction(scaling, offset); + + /*float x = scaling_ * low + offset_; + float y = scaling_ * high + offset_; + printf("%f %f (should be %f->%f)\n", x, y, GetMinOutputValue(), GetMaxOutputValue());*/ + } + + + void VolumeReslicer::FitRange(const ImageBuffer3D& image) + { + float minInputValue, maxInputValue; + + if (!image.GetRange(minInputValue, maxInputValue) || + maxInputValue < 1) + { + ResetLinearFunction(); + } + else + { + SetWindow(minInputValue, maxInputValue); + } + } + + + void VolumeReslicer::SetWindowing(ImageWindowing windowing, + const ImageBuffer3D& image, + float rescaleSlope, + float rescaleIntercept) + { + if (windowing == ImageWindowing_Custom) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); + } + + float center, width; + ComputeWindowing(center, width, windowing, 0, 0); + + float a = (center - width / 2.0f - rescaleIntercept) / rescaleSlope; + float b = (center + width / 2.0f - rescaleIntercept) / rescaleSlope; + SetWindow(a, b); + } + + + void VolumeReslicer::SetOutputFormat(Orthanc::PixelFormat format) + { + if (format != Orthanc::PixelFormat_Grayscale8 && + format != Orthanc::PixelFormat_Grayscale16 && + format != Orthanc::PixelFormat_BGRA32) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); + } + + if (hasLinearFunction_) + { + LOG(WARNING) << "Calls to VolumeReslicer::SetOutputFormat() should be done before VolumeReslicer::FitRange()"; + } + + outputFormat_ = format; + Reset(); + } + + + void VolumeReslicer::SetInterpolation(ImageInterpolation interpolation) + { + if (interpolation != ImageInterpolation_Nearest && + interpolation != ImageInterpolation_Bilinear && + interpolation != ImageInterpolation_Trilinear) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); + } + + interpolation_ = interpolation; + Reset(); + } + + + const Extent2D& VolumeReslicer::GetOutputExtent() const + { + if (success_) + { + return extent_; + } + else + { + LOG(ERROR) << "VolumeReslicer::GetOutputExtent(): (!success_)"; + + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); + } + } + + + const Orthanc::ImageAccessor& VolumeReslicer::GetOutputSlice() const + { + if (success_) + { + assert(slice_.get() != NULL); + return *slice_; + } + else + { + LOG(ERROR) << "VolumeReslicer::GetOutputSlice(): (!success_)"; + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); + } + } + + + Orthanc::ImageAccessor* VolumeReslicer::ReleaseOutputSlice() + { + if (success_) + { + assert(slice_.get() != NULL); + success_ = false; + return slice_.release(); + } + else + { + LOG(ERROR) << "VolumeReslicer::ReleaseOutputSlice(): (!success_)"; + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); + } + } + + + void VolumeReslicer::Apply(const ImageBuffer3D& source, + const VolumeImageGeometry& geometry, + const CoordinateSystem3D& plane) + { + // Choose the default voxel size as the finest voxel dimension + // of the source volumetric image + const OrthancStone::Vector dim = + geometry.GetVoxelDimensions(OrthancStone::VolumeProjection_Axial); + double voxelSize = dim[0]; + + if (dim[1] < voxelSize) + { + voxelSize = dim[1]; + } + + if (dim[2] < voxelSize) + { + voxelSize = dim[2]; + } + + if (voxelSize <= 0) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); + } + + Apply(source, geometry, plane, voxelSize); + } + + + void VolumeReslicer::Apply(const ImageBuffer3D& source, + const VolumeImageGeometry& geometry, + const CoordinateSystem3D& plane, + double voxelSize) + { + Reset(); + pixelSpacing_ = voxelSize; + + // Firstly, compute the intersection of the source volumetric + // image with the reslicing plane. This leads to a polygon with 3 + // to 6 vertices. We compute the extent of the intersection + // polygon, with respect to the coordinate system of the reslicing + // plane. + OrientedVolumeBoundingBox box(geometry); + + if (!box.ComputeExtent(extent_, plane)) + { + // The plane does not intersect with the bounding box of the volume + slice_.reset(new Orthanc::Image(outputFormat_, 0, 0, false)); + success_ = true; + return; + } + + // Secondly, the extent together with the voxel size gives the + // size of the output image + unsigned int width = boost::math::iround(extent_.GetWidth() / voxelSize); + unsigned int height = boost::math::iround(extent_.GetHeight() / voxelSize); + + slice_.reset(new Orthanc::Image(outputFormat_, width, height, false)); + + //CheckIterators(source, plane, box); + + if (fastMode_) + { + ProcessImage<FastRowIterator>(*slice_, extent_, source, plane, box, + interpolation_, hasLinearFunction_, scaling_, offset_); + } + else + { + ProcessImage<SlowRowIterator>(*slice_, extent_, source, plane, box, + interpolation_, hasLinearFunction_, scaling_, offset_); + } + + success_ = true; + } + + + double VolumeReslicer::GetPixelSpacing() const + { + if (success_) + { + return pixelSpacing_; + } + else + { + LOG(ERROR) << "VolumeReslicer::GetPixelSpacing(): (!success_)"; + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); + } + } +}