--- 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();
+  }
 }