# HG changeset patch
# User Sebastien Jodogne <s.jodogne@gmail.com>
# Date 1521216671 -3600
# Node ID 371da7fe2c0e947a3d07673b01047df5b898d9da
# Parent  46cb2eedc2e0bee31eeefc3cde11fab16d6736e7
FiniteProjectiveCamera::ApplyRaytracer

diff -r 46cb2eedc2e0 -r 371da7fe2c0e Framework/Toolbox/FiniteProjectiveCamera.cpp
--- a/Framework/Toolbox/FiniteProjectiveCamera.cpp	Fri Mar 16 15:01:52 2018 +0100
+++ b/Framework/Toolbox/FiniteProjectiveCamera.cpp	Fri Mar 16 17:11:11 2018 +0100
@@ -22,9 +22,12 @@
 #include "FiniteProjectiveCamera.h"
 
 #include "GeometryToolbox.h"
+#include "SubpixelReader.h"
 
 #include <Core/Logging.h>
 #include <Core/OrthancException.h>
+#include <Core/Images/Image.h>
+#include <Core/Images/ImageProcessing.h>
 
 namespace OrthancStone
 {
@@ -288,4 +291,168 @@
       }
     }
   }
+
+
+  template <Orthanc::PixelFormat TargetFormat,
+            Orthanc::PixelFormat SourceFormat,
+            bool MIP>
+  static void ApplyRaytracerInternal(Orthanc::ImageAccessor& target,
+                                     const FiniteProjectiveCamera& camera,
+                                     const ImageBuffer3D& source,
+                                     VolumeProjection projection)
+  {
+    if (source.GetFormat() != SourceFormat ||
+        target.GetFormat() != TargetFormat ||
+        !std::numeric_limits<float>::is_iec559 ||
+        sizeof(float) != 4)
+    {
+      throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
+    }
+
+    LOG(WARNING) << "Input volume size: " << source.GetWidth() << "x"
+                 << source.GetHeight() << "x" << source.GetDepth();
+    LOG(WARNING) << "Input pixel format: " << Orthanc::EnumerationToString(source.GetFormat());
+    LOG(WARNING) << "Output image size: " << target.GetWidth() << "x" << target.GetHeight();
+    LOG(WARNING) << "Output pixel format: " << Orthanc::EnumerationToString(target.GetFormat());
+
+    std::auto_ptr<OrthancStone::ParallelSlices> slices(source.GetGeometry(projection));
+    const OrthancStone::Vector pixelSpacing = source.GetVoxelDimensions(projection);
+    const unsigned int targetWidth = target.GetWidth();
+    const unsigned int targetHeight = target.GetHeight();
+
+    Orthanc::Image accumulator(Orthanc::PixelFormat_Float32, targetWidth, targetHeight, false);
+    Orthanc::Image counter(Orthanc::PixelFormat_Grayscale16, targetWidth, targetHeight, false);
+    Orthanc::ImageProcessing::Set(accumulator, 0);
+    Orthanc::ImageProcessing::Set(counter, 0);
+
+    typedef SubpixelReader<SourceFormat, ImageInterpolation_Nearest>  SourceReader;
+    
+    for (size_t z = 0; z < slices->GetSliceCount(); z++)
+    {
+      LOG(INFO) << "Applying raytracer on slice: " << z << "/" << slices->GetSliceCount();
+      
+      const OrthancStone::CoordinateSystem3D& slice = slices->GetSlice(z);
+      OrthancStone::ImageBuffer3D::SliceReader sliceReader(source, projection, z);
+
+      SourceReader pixelReader(sliceReader.GetAccessor());
+      
+      for (unsigned int y = 0; y < targetHeight; y++)
+      {
+        float *qacc = reinterpret_cast<float*>(accumulator.GetRow(y));
+        uint16_t *qcount = reinterpret_cast<uint16_t*>(counter.GetRow(y));
+
+        for (unsigned int x = 0; x < targetWidth; x++)
+        {
+          // Backproject the ray originating from the center of the target pixel
+          OrthancStone::Vector direction = camera.GetRayDirection(static_cast<double>(x + 0.5),
+                                                                  static_cast<double>(y + 0.5));
+
+          // Compute the 3D intersection of the ray with the slice plane
+          OrthancStone::Vector p;
+          if (slice.IntersectLine(p, camera.GetCenter(), direction))
+          {
+            // Compute the 2D coordinates of the intersections, in slice coordinates
+            double ix, iy;
+            slice.ProjectPoint(ix, iy, p);
+
+            ix /= pixelSpacing[0];
+            iy /= pixelSpacing[1];
+
+            // Read and accumulate the value of the pixel
+            float pixel;
+            if (pixelReader.GetFloatValue(pixel, ix, iy))
+            {
+              if (MIP)
+              {
+                // MIP rendering
+                if (*qcount == 0)
+                {
+                  (*qacc) = pixel;
+                  (*qcount) = 1;
+                }
+                else if (pixel > *qacc)
+                {
+                  (*qacc) = pixel;
+                }
+              }
+              else
+              {
+                // Mean intensity
+                (*qacc) += pixel;
+                (*qcount) ++;
+              }
+            }
+          }
+
+          qacc++;
+          qcount++;
+        }
+      }
+    }
+
+
+    typedef Orthanc::PixelTraits<TargetFormat>  TargetTraits;
+
+    // "Flatten" the accumulator image to create the target image
+    for (unsigned int y = 0; y < targetHeight; 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*>(target.GetRow(y));
+
+      for (unsigned int x = 0; x < targetWidth; x++)
+      {
+        if (*qcount == 0)
+        {
+          TargetTraits::SetZero(*p);
+        }
+        else
+        {
+          TargetTraits::FloatToPixel(*p, *qacc / static_cast<float>(*qcount));
+        }
+        
+        p++;
+        qacc++;
+        qcount++;
+      }
+    }
+  }
+
+
+  Orthanc::ImageAccessor*
+  FiniteProjectiveCamera::ApplyRaytracer(const ImageBuffer3D& source,
+                                         Orthanc::PixelFormat targetFormat,
+                                         unsigned int targetWidth,
+                                         unsigned int targetHeight,
+                                         bool mip) const
+  {
+    // TODO - We consider the axial projection of the volume, but we
+    // should choose the projection that is the "most perpendicular"
+    // to the line joining the camera center and the principal point
+    const VolumeProjection projection = VolumeProjection_Axial;
+
+    std::auto_ptr<Orthanc::ImageAccessor> target
+      (new Orthanc::Image(targetFormat, targetWidth, targetHeight, false));
+    
+    if (targetFormat == Orthanc::PixelFormat_Grayscale16 &&
+        source.GetFormat() == Orthanc::PixelFormat_Grayscale16 && mip)
+    {
+      ApplyRaytracerInternal<Orthanc::PixelFormat_Grayscale16,
+                             Orthanc::PixelFormat_Grayscale16, true>
+        (*target, *this, source, projection);
+    }
+    else if (targetFormat == Orthanc::PixelFormat_Grayscale16 &&
+             source.GetFormat() == Orthanc::PixelFormat_Grayscale16 && !mip)
+    {
+      ApplyRaytracerInternal<Orthanc::PixelFormat_Grayscale16,
+                             Orthanc::PixelFormat_Grayscale16, false>
+        (*target, *this, source, projection);
+    }
+    else
+    {
+      throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
+    }
+
+    return target.release();
+  }
 }
diff -r 46cb2eedc2e0 -r 371da7fe2c0e Framework/Toolbox/FiniteProjectiveCamera.h
--- a/Framework/Toolbox/FiniteProjectiveCamera.h	Fri Mar 16 15:01:52 2018 +0100
+++ b/Framework/Toolbox/FiniteProjectiveCamera.h	Fri Mar 16 17:11:11 2018 +0100
@@ -22,6 +22,7 @@
 #pragma once
 
 #include "LinearAlgebra.h"
+#include "../Volumes/ImageBuffer3D.h"
 
 namespace OrthancStone
 {
@@ -106,5 +107,11 @@
     // Apply the camera to a 3D point "v" that is possibly at
     // infinity. The result is a 2D point in homogeneous coordinates.
     Vector ApplyGeneral(const Vector& v) const;
+
+    Orthanc::ImageAccessor* ApplyRaytracer(const ImageBuffer3D& source,
+                                           Orthanc::PixelFormat targetFormat,
+                                           unsigned int targetWidth,
+                                           unsigned int targetHeight,
+                                           bool mip) const;
   };
 }
diff -r 46cb2eedc2e0 -r 371da7fe2c0e Framework/Volumes/ImageBuffer3D.cpp
--- a/Framework/Volumes/ImageBuffer3D.cpp	Fri Mar 16 15:01:52 2018 +0100
+++ b/Framework/Volumes/ImageBuffer3D.cpp	Fri Mar 16 17:11:11 2018 +0100
@@ -30,7 +30,7 @@
 namespace OrthancStone
 {
   Orthanc::ImageAccessor ImageBuffer3D::GetAxialSliceAccessor(unsigned int slice,
-                                                              bool readOnly)
+                                                              bool readOnly) const
   {
     if (slice >= depth_)
     {
@@ -55,7 +55,7 @@
 
 
   Orthanc::ImageAccessor ImageBuffer3D::GetCoronalSliceAccessor(unsigned int slice,
-                                                                bool readOnly)
+                                                                bool readOnly) const
   {
     if (slice >= height_)
     {
@@ -353,7 +353,7 @@
   }
 
 
-  ImageBuffer3D::SliceReader::SliceReader(ImageBuffer3D& that,
+  ImageBuffer3D::SliceReader::SliceReader(const ImageBuffer3D& that,
                                           VolumeProjection projection,
                                           unsigned int slice)
   {
diff -r 46cb2eedc2e0 -r 371da7fe2c0e Framework/Volumes/ImageBuffer3D.h
--- a/Framework/Volumes/ImageBuffer3D.h	Fri Mar 16 15:01:52 2018 +0100
+++ b/Framework/Volumes/ImageBuffer3D.h	Fri Mar 16 17:11:11 2018 +0100
@@ -49,10 +49,10 @@
     void ExtendImageRange(const Orthanc::ImageAccessor& slice);
 
     Orthanc::ImageAccessor GetAxialSliceAccessor(unsigned int slice,
-                                                 bool readOnly);
+                                                 bool readOnly) const;
 
     Orthanc::ImageAccessor GetCoronalSliceAccessor(unsigned int slice,
-                                                   bool readOnly);
+                                                   bool readOnly) const;
 
     Orthanc::Image*  ExtractSagittalSlice(unsigned int slice) const;
 
@@ -172,7 +172,7 @@
       std::auto_ptr<Orthanc::Image>  sagittal_;  // Unused for axial and coronal
 
     public:
-      SliceReader(ImageBuffer3D& that,
+      SliceReader(const ImageBuffer3D& that,
                   VolumeProjection projection,
                   unsigned int slice);