# HG changeset patch
# User Sebastien Jodogne <s.jodogne@gmail.com>
# Date 1521551139 -3600
# Node ID 4abddd08337486cb431ad66c2dc75e11f4040146
# Parent  371da7fe2c0e947a3d07673b01047df5b898d9da
ShearWarpProjectiveTransform::ApplyAxial()

diff -r 371da7fe2c0e -r 4abddd083374 Framework/Toolbox/ShearWarpProjectiveTransform.cpp
--- 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();
+  }
 }
diff -r 371da7fe2c0e -r 4abddd083374 Framework/Toolbox/ShearWarpProjectiveTransform.h
--- a/Framework/Toolbox/ShearWarpProjectiveTransform.h	Fri Mar 16 17:11:11 2018 +0100
+++ b/Framework/Toolbox/ShearWarpProjectiveTransform.h	Tue Mar 20 14:05:39 2018 +0100
@@ -88,5 +88,17 @@
     static Matrix CalibrateView(const Vector& camera,
                                 const Vector& principalPoint,
                                 double angle);
+
+    static Orthanc::ImageAccessor* ApplyAxial(float& maxValue,
+                                              const Matrix& M_view,  // cf. "CalibrateView()"
+                                              const ImageBuffer3D& source,
+                                              Orthanc::PixelFormat targetFormat,
+                                              unsigned int targetWidth,
+                                              unsigned int targetHeight,
+                                              bool mip,
+                                              double pixelSpacing,
+                                              unsigned int countSlices,
+                                              ImageInterpolation shearInterpolation,
+                                              ImageInterpolation warpInterpolation);
   };
 }
diff -r 371da7fe2c0e -r 4abddd083374 Resources/Computations/ComputeWarp.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Resources/Computations/ComputeWarp.py	Tue Mar 20 14:05:39 2018 +0100
@@ -0,0 +1,112 @@
+#!/usr/bin/python
+
+from sympy import *
+from sympy.solvers import solve
+import pprint
+import sys
+
+init_printing(use_unicode=True)
+
+
+# Create a test 3D vector using homogeneous coordinates
+x, y, z, w = symbols('x y z w')
+p = Matrix([ x, y, z, w ])
+
+
+# Create a shear matrix, and a scale/shift "T * S" transform as in
+# Lacroute's thesis (Equation A.16, page 209)
+ex, ey, ew = symbols('ex ey ew')
+sx, sy, tx, ty = symbols('sx sy tx ty')
+
+TS = Matrix([[ sx, 0,  0, tx ],
+             [ 0,  sy, 0, ty ],
+             [ 0,  0,  1, 0  ],
+             [ 0,  0,  0, 1  ]])
+
+pureShear = Matrix([[ 1, 0, ex, 0 ],
+                    [ 0, 1, ey, 0 ],
+                    [ 0, 0, 1,  0 ],
+                    [ 0, 0, ew, 1 ]])
+
+
+# Create a general warp matrix, that corresponds to "M_warp" in
+# Equation (A.17) of Lacroute's thesis:
+ww11, ww12, ww13, ww14, ww21, ww22, ww23, ww24, ww31, ww32, ww33, ww34, ww41, ww42, ww43, ww44 = symbols('ww11 ww12 ww13 ww14 ww21 ww22 ww23 ww24 ww31 ww32 ww33 ww34 ww41 ww42 ww43 ww44')
+
+WW = Matrix([[ ww11, ww12, ww13, ww14 ],
+             [ ww21, ww22, ww23, ww24 ],
+             [ ww31, ww32, ww33, ww34 ],
+             [ ww41, ww43, ww43, ww44 ]])
+
+
+# Create the matrix of intrinsic parameters of the camera
+k11, k22, k14, k24 = symbols('k11 k22 k14 k24')
+K = Matrix([[ k11, 0,   0, k14 ],
+            [ 0,   k22, 0, k24 ],
+            [ 0,   0,   0, 1   ]])
+
+
+# The full decomposition is:
+M_shear = TS * pureShear
+M_warp = K * WW * TS.inv()
+AA = M_warp * M_shear
+
+# Check that the central component "M_warp == K * WW * TS.inv()" that
+# is the left part of "A" is another general warp matrix (i.e. no
+# exception is thrown about incompatible matrix sizes):
+M_warp * p
+
+if (M_warp.cols != 4 or
+    M_warp.rows != 3):
+    raise Exception('Invalid matrix size')
+
+
+# We've just shown that "M_warp" is a general 3x4 matrix. Let's call
+# it W:
+w11, w12, w13, w14, w21, w22, w23, w24, w41, w42, w43, w44 = symbols('w11 w12 w13 w14 w21 w22 w23 w24 w41 w42 w43 w44')
+
+W = Matrix([[ w11, w12, w13, w14 ],
+            [ w21, w22, w23, w24 ],
+            [ w41, w43, w43, w44 ]])
+
+# This shows that it is sufficient to study a decomposition of the
+# following form:
+A = W * M_shear
+print('\nA = W * M_shear =')
+pprint.pprint(A)
+
+sys.stdout.write('\nW = ')
+pprint.pprint(W)
+
+sys.stdout.write('\nM_shear = ')
+pprint.pprint(M_shear)
+
+
+
+# Let's consider one fixed 2D point (i,j) in the intermediate
+# image. The 3D points (x,y,z,1) that are mapped to (i,j) must satisfy
+# the equation "(i,j) == M_shear * (x,y,z,w)". As "M_shear" is
+# invertible, we solve "(x,y,z,w) == inv(M_shear) * (i,j,k,1)".
+
+i, j, k = symbols('i j k')
+l = M_shear.inv() * Matrix([ i, j, k, 1 ])
+
+print('\nLocus for points imaged to some fixed (i,j,k,l) point in the intermediate image:')
+print('x = %s' % l[0])
+print('y = %s' % l[1])
+print('z = %s' % l[2])
+print('w = %s' % l[3])
+
+
+# By inspecting the 4 equations above, we see that the locus entirely
+# depends upon the "k" value that encodes the Z-axis
+
+print('\nGlobal effect of the shear-warp transform on this locus:')
+q = expand(A * l)
+pprint.pprint(q)
+
+print("\nWe can arbitrarily fix the value of 'k', so let's choose 'k=0':")
+pprint.pprint(q.subs(k, 0))
+
+print("\nThis gives the warp transform.")
+print("QED: line after Equation (A.17) on page 209.\n")