# HG changeset patch
# User Benjamin Golinvaux <bgo@osimis.io>
# Date 1580131742 -3600
# Node ID 5a2d5380148d4068af55ca3e895f5fa18c778bc1
# Parent  69177b10e2b9c5ccf9ed069fac8b683fc9120b69
Added outlier rejection for min max computation in multiframe volume loader

diff -r 69177b10e2b9 -r 5a2d5380148d Framework/Loaders/OrthancMultiframeVolumeLoader.cpp
--- a/Framework/Loaders/OrthancMultiframeVolumeLoader.cpp	Tue Jan 21 16:52:37 2020 +0100
+++ b/Framework/Loaders/OrthancMultiframeVolumeLoader.cpp	Mon Jan 27 14:29:02 2020 +0100
@@ -263,7 +263,8 @@
   }
 
   template <typename T>
-  void OrthancMultiframeVolumeLoader::CopyPixelData(const std::string& pixelData)
+  void OrthancMultiframeVolumeLoader::CopyPixelDataAndComputeDistribution(
+    const std::string& pixelData, std::map<T,uint64_t>& distribution)
   {
     ImageBuffer3D& target = volume_->GetPixelData();
       
@@ -283,43 +284,209 @@
       return;
     }
 
-    const uint8_t* source = reinterpret_cast<const uint8_t*>(pixelData.c_str());
-
-    for (unsigned int z = 0; z < depth; z++)
+    // first pass to initialize map
     {
-      ImageBuffer3D::SliceWriter writer(target, VolumeProjection_Axial, z);
+      const uint8_t* source = reinterpret_cast<const uint8_t*>(pixelData.c_str());
 
-      assert (writer.GetAccessor().GetWidth() == width &&
-              writer.GetAccessor().GetHeight() == height);
+      for (unsigned int z = 0; z < depth; z++)
+      {
+        for (unsigned int y = 0; y < height; y++)
+        {
+          for (unsigned int x = 0; x < width; x++)
+          {
+            T value;
+            CopyPixel(value, source);
+            distribution[value] = 0;
+            source += bpp;
+          }
+        }
+      }
+    }
+
+    {
+      const uint8_t* source = reinterpret_cast<const uint8_t*>(pixelData.c_str());
 
-      for (unsigned int y = 0; y < height; y++)
+      for (unsigned int z = 0; z < depth; z++)
       {
-        assert(sizeof(T) == Orthanc::GetBytesPerPixel(target.GetFormat()));
+        ImageBuffer3D::SliceWriter writer(target, VolumeProjection_Axial, z);
 
-        T* target = reinterpret_cast<T*>(writer.GetAccessor().GetRow(y));
+        assert(writer.GetAccessor().GetWidth() == width &&
+          writer.GetAccessor().GetHeight() == height);
+
+        for (unsigned int y = 0; y < height; y++)
+        {
+          assert(sizeof(T) == Orthanc::GetBytesPerPixel(target.GetFormat()));
 
-        for (unsigned int x = 0; x < width; x++)
-        {
-          CopyPixel(*target, source);
-          target ++;
-          source += bpp;
+          T* target = reinterpret_cast<T*>(writer.GetAccessor().GetRow(y));
+
+          for (unsigned int x = 0; x < width; x++)
+          {
+            CopyPixel(*target, source);
+
+            distribution[*target] += 1;
+
+            target++;
+            source += bpp;
+          }
         }
       }
     }
   }
 
+  template <typename T>
+  void OrthancMultiframeVolumeLoader::ComputeMinMaxWithOutlierRejection(
+    const std::map<T, uint64_t>& distribution)
+  {
+    if (distribution.size() == 0)
+    {
+      LOG(ERROR) << "ComputeMinMaxWithOutlierRejection -- Volume image empty.";
+    }
+    else
+    {
+      ImageBuffer3D& target = volume_->GetPixelData();
+
+      const uint64_t bpp = target.GetBytesPerPixel();
+      const uint64_t width = target.GetWidth();
+      const uint64_t height = target.GetHeight();
+      const uint64_t depth = target.GetDepth();
+      const uint64_t voxelCount = width * height * depth;
+
+      // now that we have distribution[pixelValue] == numberOfPixelsWithValue
+      // compute number of values and check (assertion) that it is equal to 
+      // width * height * depth 
+      {
+        typename std::map<T, uint64_t>::const_iterator it = distribution.begin();
+        uint64_t totalCount = 0;
+        distributionRawMin_ = static_cast<float>(it->first);
+
+        while (it != distribution.end())
+        {
+          T pixelValue = it->first;
+          uint64_t count = it->second;
+          totalCount += count;
+          it++;
+          if (it == distribution.end())
+            distributionRawMax_ = static_cast<float>(pixelValue);
+        }
+        LOG(INFO) << "Volume image. First distribution value = " 
+          << static_cast<float>(distributionRawMin_) 
+          << " | Last distribution value = " 
+          << static_cast<float>(distributionRawMax_);
+
+        if (totalCount != voxelCount)
+        {
+          LOG(ERROR) << "Internal error in dose distribution computation. TC (" 
+            << totalCount << ") != VoxC (" << voxelCount;
+          throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
+        }
+      }
+
+      // compute the number of voxels to reject at each end of the distribution
+      uint64_t endRejectionCount = static_cast<uint64_t>(
+        outliersHalfRejectionRate_ * voxelCount);
+
+      if (endRejectionCount > voxelCount)
+      {
+        LOG(ERROR) << "Internal error in dose distribution computation."
+          << " endRejectionCount = " << endRejectionCount
+          << " | voxelCount = " << voxelCount;
+        throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
+      }
+
+      // this will contain the actual distribution minimum after outlier 
+      // rejection
+      T resultMin = 0;
+
+      // then start from start and remove pixel values up to 
+      // endRejectionCount voxels rejected
+      {
+        typename std::map<T, uint64_t>::const_iterator it = distribution.begin();
+        
+        uint64_t currentCount = 0;
+
+        while (it != distribution.end())
+        {
+          T pixelValue = it->first;
+          uint64_t count = it->second;
+
+          // if this pixelValue crosses the rejection threshold, let's set it
+          // and exit the loop
+          if ((currentCount <= endRejectionCount) &&
+              (currentCount + count > endRejectionCount))
+          {
+            resultMin = pixelValue;
+            break;
+          }
+          else
+          {
+            currentCount += count;
+          }
+          // and continue walking along the distribution
+          it++;
+        }
+      }
+
+      // this will contain the actual distribution maximum after outlier 
+      // rejection
+      T resultMax = 0;
+      // now start from END and remove pixel values up to 
+      // endRejectionCount voxels rejected
+      {
+        typename std::map<T, uint64_t>::const_reverse_iterator it = distribution.rbegin();
+
+        uint64_t currentCount = 0;
+
+        while (it != distribution.rend())
+        {
+          T pixelValue = it->first;
+          uint64_t count = it->second;
+
+          if ((currentCount <= endRejectionCount) &&
+              (currentCount + count > endRejectionCount))
+          {
+            resultMax = pixelValue;
+            break;
+          }
+          else
+          {
+            currentCount += count;
+          }
+          // and continue walking along the distribution
+          it++;
+        }
+      }
+      if (resultMin > resultMax)
+      {
+        LOG(ERROR) << "Internal error in dose distribution computation! " << 
+          "resultMin (" << resultMin << ") > resultMax (" << resultMax << ")";
+        throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
+      }
+      computedDistributionMin_ = static_cast<float>(resultMin);
+      computedDistributionMax_ = static_cast<float>(resultMax);
+    }
+  }
+
+  template <typename T>
+  void OrthancMultiframeVolumeLoader::CopyPixelDataAndComputeMinMax(
+    const std::string& pixelData)
+  {
+    std::map<T, uint64_t> distribution;
+    CopyPixelDataAndComputeDistribution(pixelData, distribution);
+    ComputeMinMaxWithOutlierRejection(distribution);
+  }
+
   void OrthancMultiframeVolumeLoader::SetUncompressedPixelData(const std::string& pixelData)
   {
     switch (volume_->GetPixelData().GetFormat())
     {
       case Orthanc::PixelFormat_Grayscale32:
-        CopyPixelData<uint32_t>(pixelData);
+        CopyPixelDataAndComputeMinMax<uint32_t>(pixelData);
         break;
       case Orthanc::PixelFormat_Grayscale16:
-        CopyPixelData<uint16_t>(pixelData);
+        CopyPixelDataAndComputeMinMax<uint16_t>(pixelData);
         break;
       case Orthanc::PixelFormat_SignedGrayscale16:
-        CopyPixelData<int16_t>(pixelData);
+        CopyPixelDataAndComputeMinMax<int16_t>(pixelData);
         break;
       default:
         throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
@@ -341,13 +508,20 @@
     return volume_->GetGeometry();
   }
 
-  OrthancMultiframeVolumeLoader::OrthancMultiframeVolumeLoader(boost::shared_ptr<DicomVolumeImage> volume,
-                                                               IOracle& oracle,
-                                                               IObservable& oracleObservable) :
+  OrthancMultiframeVolumeLoader::OrthancMultiframeVolumeLoader(
+    boost::shared_ptr<DicomVolumeImage> volume,
+    IOracle& oracle,
+    IObservable& oracleObservable,
+    float outliersHalfRejectionRate) :
     LoaderStateMachine(oracle, oracleObservable),
     IObservable(oracleObservable.GetBroker()),
     volume_(volume),
-    pixelDataLoaded_(false)
+    pixelDataLoaded_(false),
+    outliersHalfRejectionRate_(outliersHalfRejectionRate),
+    distributionRawMin_(0),
+    distributionRawMax_(0),
+    computedDistributionMin_(0),
+    computedDistributionMax_(0)
   {
     if (volume.get() == NULL)
     {
@@ -360,6 +534,29 @@
     LOG(TRACE) << "OrthancMultiframeVolumeLoader::~OrthancMultiframeVolumeLoader()";
   }
 
+
+  void OrthancMultiframeVolumeLoader::GetDistributionMinMax
+  (float& minValue, float& maxValue) const
+  {
+    if (distributionRawMin_ == 0 && distributionRawMax_ == 0)
+    {
+      LOG(WARNING) << "GetDistributionMinMaxWithOutliersRejection called before computation!";
+    }
+    minValue = distributionRawMin_;
+    maxValue = distributionRawMax_;
+  }
+  
+  void OrthancMultiframeVolumeLoader::GetDistributionMinMaxWithOutliersRejection
+    (float& minValue, float& maxValue) const
+  {
+    if (computedDistributionMin_ == 0 && computedDistributionMax_ == 0)
+    {
+      LOG(WARNING) << "GetDistributionMinMaxWithOutliersRejection called before computation!";
+    }
+    minValue = computedDistributionMin_;
+    maxValue = computedDistributionMax_;
+  }
+
   void OrthancMultiframeVolumeLoader::LoadInstance(const std::string& instanceId)
   {
     Start();
diff -r 69177b10e2b9 -r 5a2d5380148d Framework/Loaders/OrthancMultiframeVolumeLoader.h
--- a/Framework/Loaders/OrthancMultiframeVolumeLoader.h	Tue Jan 21 16:52:37 2020 +0100
+++ b/Framework/Loaders/OrthancMultiframeVolumeLoader.h	Mon Jan 27 14:29:02 2020 +0100
@@ -43,6 +43,11 @@
     std::string                          instanceId_;
     std::string                          transferSyntaxUid_;
     bool                                 pixelDataLoaded_;
+    float                                outliersHalfRejectionRate_;
+    float                                distributionRawMin_;
+    float                                distributionRawMax_;
+    float                                computedDistributionMin_;
+    float                                computedDistributionMax_;
 
     const std::string& GetInstanceId() const;
 
@@ -52,8 +57,35 @@
 
     void SetGeometry(const Orthanc::DicomMap& dicom);
 
+
+    /**
+    This method will :
+    
+    - copy the pixel values from the response to the volume image
+    - compute the maximum and minimum value while discarding the
+      outliersHalfRejectionRate_ fraction of the outliers from both the start 
+      and the end of the distribution.
+
+      In English, this means that, if the volume dataset contains a few extreme
+      values very different from the rest (outliers) that we want to get rid of,
+      this method allows to do so.
+
+      If you supply 0.005, for instance, it means 1% of the extreme values will
+      be rejected (0.5% on each side of the distribution)
+    */
     template <typename T>
-    void CopyPixelData(const std::string& pixelData);
+    void CopyPixelDataAndComputeMinMax(const std::string& pixelData);
+      
+    /** Service method for CopyPixelDataAndComputeMinMax*/
+    template <typename T>
+    void CopyPixelDataAndComputeDistribution(
+      const std::string& pixelData, 
+      std::map<T, uint64_t>& distribution);
+
+    /** Service method for CopyPixelDataAndComputeMinMax*/
+    template <typename T>
+    void ComputeMinMaxWithOutlierRejection(
+      const std::map<T, uint64_t>& distribution);
 
     void SetUncompressedPixelData(const std::string& pixelData);
 
@@ -63,7 +95,8 @@
   public:
     OrthancMultiframeVolumeLoader(boost::shared_ptr<DicomVolumeImage> volume,
                                   IOracle& oracle,
-                                  IObservable& oracleObservable);
+                                  IObservable& oracleObservable,
+                                  float outliersHalfRejectionRate = 0.0005);
     
     virtual ~OrthancMultiframeVolumeLoader();
 
@@ -72,6 +105,12 @@
       return pixelDataLoaded_;
     }
 
+    void GetDistributionMinMax
+      (float& minValue, float& maxValue) const;
+
+    void GetDistributionMinMaxWithOutliersRejection
+      (float& minValue, float& maxValue) const;
+
     void LoadInstance(const std::string& instanceId);
   };
 }