Mercurial > hg > orthanc-stone
diff Framework/Toolbox/ShearWarpProjectiveTransform.cpp @ 193:4abddd083374 wasm
ShearWarpProjectiveTransform::ApplyAxial()
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Tue, 20 Mar 2018 14:05:39 +0100 |
parents | 46cb2eedc2e0 |
children | dabe9982fca3 |
line wrap: on
line diff
--- a/Framework/Toolbox/ShearWarpProjectiveTransform.cpp Fri Mar 16 17:11:11 2018 +0100 +++ b/Framework/Toolbox/ShearWarpProjectiveTransform.cpp Tue Mar 20 14:05:39 2018 +0100 @@ -26,10 +26,13 @@ #include "FiniteProjectiveCamera.h" #include "GeometryToolbox.h" +#include <Core/Images/PixelTraits.h> +#include <Core/Images/ImageProcessing.h> #include <Core/OrthancException.h> #include <Core/Logging.h> #include <boost/numeric/ublas/matrix_proxy.hpp> +#include <boost/math/special_functions/round.hpp> #include <cassert> @@ -326,4 +329,319 @@ return M_view; } + + + template <Orthanc::PixelFormat SourceFormat, + Orthanc::PixelFormat TargetFormat, + bool MIP> + static void ApplyAxialInternal(Orthanc::ImageAccessor& target, + float& maxValue, + const Matrix& M_view, + const ImageBuffer3D& source, + double pixelSpacing, + unsigned int countSlices, + ImageInterpolation shearInterpolation, + ImageInterpolation warpInterpolation) + { + typedef Orthanc::PixelTraits<SourceFormat> SourceTraits; + typedef Orthanc::PixelTraits<TargetFormat> TargetTraits; + + /** + * Step 1: Precompute some information. + **/ + + if (target.GetFormat() != TargetFormat || + source.GetFormat() != SourceFormat || + !std::numeric_limits<float>::is_iec559 || + sizeof(float) != 4) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); + } + + if (countSlices > source.GetDepth()) + { + countSlices = source.GetDepth(); + } + + if (countSlices == 0) + { + maxValue = 0; + Orthanc::ImageProcessing::Set(target, 0); + return; + } + + LOG(INFO) << "Number of rendered slices: " << countSlices; + + + /** + * Step 2: Extract the shear-warp transform corresponding to + * M_view. + **/ + + // Compute the "world" matrix that maps the source volume to the + // (0,0,0)->(1,1,1) unit cube + Vector origin = source.GetCoordinates(0, 0, 0); + Vector ps = source.GetVoxelDimensions(VolumeProjection_Axial); + Matrix world = LinearAlgebra::Product( + GeometryToolbox::CreateScalingMatrix(1.0 / ps[0], 1.0 / ps[1], 1.0 / ps[2]), + GeometryToolbox::CreateTranslationMatrix(-origin[0], -origin[1], -origin[2])); + + Matrix worldInv; + LinearAlgebra::InvertMatrix(worldInv, world); + + ShearWarpProjectiveTransform shearWarp(LinearAlgebra::Product(M_view, worldInv), + /*LinearAlgebra::IdentityMatrix(4),*/ + source.GetWidth(), + source.GetHeight(), + source.GetDepth(), + pixelSpacing, pixelSpacing, + target.GetWidth(), target.GetHeight()); + + const unsigned int intermediateWidth = shearWarp.GetIntermediateWidth(); + const unsigned int intermediateHeight = shearWarp.GetIntermediateHeight(); + + + /** + * Step 3: Apply the "shear" part of the transform to form the + * intermediate image. The sheared images are accumulated into the + * Float32 image "accumulator". The number of samples available + * for each pixel is stored in the "counter" image. + **/ + + std::auto_ptr<Orthanc::ImageAccessor> accumulator, counter, intermediate; + + accumulator.reset(new Orthanc::Image(Orthanc::PixelFormat_Float32, + intermediateWidth, intermediateHeight, false)); + counter.reset(new Orthanc::Image(Orthanc::PixelFormat_Grayscale16, + intermediateWidth, intermediateHeight, false)); + intermediate.reset(new Orthanc::Image(SourceFormat, intermediateWidth, intermediateHeight, false)); + + Orthanc::ImageProcessing::Set(*accumulator, 0); + Orthanc::ImageProcessing::Set(*counter, 0); + + // Loop around the slices of the volume + for (unsigned int i = 0; i <= countSlices; i++) + { + // (3.a) Compute the shear for this specific slice + unsigned int z = static_cast<unsigned int>( + boost::math::iround(static_cast<double>(i) / + static_cast<double>(countSlices) * + static_cast<double>(source.GetDepth() - 1))); + + double a11, b1, a22, b2, vz; + shearWarp.ComputeShearOnSlice(a11, b1, a22, b2, vz, static_cast<double>(z) + 0.5); + + + { + // (3.b) Detect the "useful" portion of the intermediate image + // for this slice (i.e. the bounding box where the source + // slice is mapped to by the shear), so as to update "counter" + Matrix a = LinearAlgebra::ZeroMatrix(3, 3); + a(0,0) = a11; + a(0,2) = b1; + a(1,1) = a22; + a(1,2) = b2; + a(2,2) = 1; + + unsigned int x1, y1, x2, y2; + if (GetProjectiveTransformExtent(x1, y1, x2, y2, a, + source.GetWidth(), source.GetHeight(), + intermediateWidth, intermediateHeight)) + { + for (unsigned int y = y1; y <= y2; y++) + { + uint16_t* p = reinterpret_cast<uint16_t*>(counter->GetRow(y)) + x1; + for (unsigned int x = x1; x <= x2; x++, p++) + { + if (MIP) + { + // TODO - In the case of MIP, "counter" could be + // reduced to "PixelFormat_Grayscale8" to reduce + // memory usage + *p = 1; + } + else + { + *p += 1; + } + } + } + } + } + + + { + // (3.c) Shear the source slice into a temporary image + ImageBuffer3D::SliceReader reader(source, VolumeProjection_Axial, z); + ApplyAffineTransform(*intermediate, reader.GetAccessor(), + a11, 0, b1, + 0, a22, b2, + shearInterpolation); + } + + + for (unsigned int y = 0; y < intermediateHeight; y++) + { + // (3.d) Accumulate the pixels of the sheared image into "accumulator" + const typename SourceTraits::PixelType* p = + reinterpret_cast<const typename SourceTraits::PixelType*>(intermediate->GetConstRow(y)); + + float* q = reinterpret_cast<float*>(accumulator->GetRow(y)); + + for (unsigned int x = 0; x < intermediateWidth; x++) + { + float pixel = SourceTraits::PixelToFloat(*p); + + if (MIP) + { + // Get maximum for MIP + if (*q < pixel) + { + *q = pixel; + } + } + else + { + *q += pixel; + } + + p++; + q++; + } + } + } + + + /** + * Step 4: The intermediate image (that will be transformed by the + * "warp") is now available as an accumulator image together with + * a counter image. "Flatten" these two images into one. + **/ + + intermediate.reset(new Orthanc::Image + (TargetFormat, intermediateWidth, intermediateHeight, false)); + + maxValue = 0; + + for (unsigned int y = 0; y < intermediateHeight; y++) + { + const float *qacc = reinterpret_cast<const float*>(accumulator->GetConstRow(y)); + const uint16_t *qcount = reinterpret_cast<const uint16_t*>(counter->GetConstRow(y)); + typename TargetTraits::PixelType *p = + reinterpret_cast<typename TargetTraits::PixelType*>(intermediate->GetRow(y)); + + for (unsigned int x = 0; x < intermediateWidth; x++) + { + if (*qcount == 0) + { + TargetTraits::SetZero(*p); + } + else + { + *p = *qacc / static_cast<float>(*qcount); + + if (*p > maxValue) + { + maxValue = *p; + } + } + + p++; + qacc++; + qcount++; + } + } + + // We don't need the accumulator images anymore + accumulator.reset(NULL); + counter.reset(NULL); + + + /** + * Step 6: Apply the "warp" part of the transform to map the + * intermediate image to the final image. + **/ + + Matrix warp; + + { + // (5.a) Compute the "warp" matrix by removing the 3rd row and + // 3rd column from the GetWarp() matrix + // Check out: ../../Resources/Computations/ComputeWarp.py + + Matrix fullWarp = LinearAlgebra::Product + (shearWarp.GetIntrinsicParameters(), shearWarp.GetWarp()); + + const double v[] = { + fullWarp(0,0), fullWarp(0,1), fullWarp(0,3), + fullWarp(1,0), fullWarp(1,1), fullWarp(1,3), + fullWarp(2,0), fullWarp(2,1), fullWarp(2,3) + }; + + LinearAlgebra::FillMatrix(warp, 3, 3, v); + } + + // (5.b) Apply the projective transform to the image + ApplyProjectiveTransform(target, *intermediate, warp, warpInterpolation); + } + + + template <Orthanc::PixelFormat SourceFormat, + Orthanc::PixelFormat TargetFormat> + static void ApplyAxialInternal2(Orthanc::ImageAccessor& target, + float& maxValue, + const Matrix& M_view, + const ImageBuffer3D& source, + bool mip, + double pixelSpacing, + unsigned int countSlices, + ImageInterpolation shearInterpolation, + ImageInterpolation warpInterpolation) + { + if (mip) + { + ApplyAxialInternal<SourceFormat, TargetFormat, true> + (target, maxValue, M_view, source, pixelSpacing, + countSlices, shearInterpolation, warpInterpolation); + } + else + { + ApplyAxialInternal<SourceFormat, TargetFormat, false> + (target, maxValue, M_view, source, pixelSpacing, + countSlices, shearInterpolation, warpInterpolation); + } + } + + + Orthanc::ImageAccessor* + ShearWarpProjectiveTransform::ApplyAxial(float& maxValue, + const Matrix& M_view, + const ImageBuffer3D& source, + Orthanc::PixelFormat targetFormat, + unsigned int targetWidth, + unsigned int targetHeight, + bool mip, + double pixelSpacing, + unsigned int countSlices, + ImageInterpolation shearInterpolation, + ImageInterpolation warpInterpolation) + { + std::auto_ptr<Orthanc::ImageAccessor> target + (new Orthanc::Image(targetFormat, targetWidth, targetHeight, false)); + + if (source.GetFormat() == Orthanc::PixelFormat_Grayscale16 && + targetFormat == Orthanc::PixelFormat_Grayscale16) + { + ApplyAxialInternal2<Orthanc::PixelFormat_Grayscale16, + Orthanc::PixelFormat_Grayscale16> + (*target, maxValue, M_view, source, mip, pixelSpacing, + countSlices, shearInterpolation, warpInterpolation); + } + else + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); + } + + return target.release(); + } }