orthanc-stone: Framework/Toolbox/ShearWarpProjectiveTransform.cpp comparison

comparison 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

comparison

equal deleted inserted replaced

-:371da7fe2c0e
+:4abddd083374
 #include "ImageGeometry.h"
 #include "Extent2D.h"
 #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>
 namespace OrthancStone
 {
 }
 }
 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();
+}
 }

Mercurial > hg > orthanc-stone

comparison Framework/Toolbox/ShearWarpProjectiveTransform.cpp @ 193:4abddd083374 wasm