[hdr] add Grossberg algorithm

alicevision · May 24, 2019 · 1f9daa4 · 1f9daa4
1 parent 87a630b
commit 1f9daa4
Show file tree

Hide file tree

Showing 25 changed files with 18,034 additions and 292 deletions.
diff --git a/CMakeLists.txt.user.2fd304d b/CMakeLists.txt.user.2fd304d
diff --git a/CMakeLists.txt.user.4b48bbd b/CMakeLists.txt.user.4b48bbd
diff --git a/CMakeLists.txt.user.65857b6 b/CMakeLists.txt.user.65857b6
diff --git a/CMakeLists.txt.user.g21835 b/CMakeLists.txt.user.g21835
diff --git a/src/aliceVision/hdr/CMakeLists.txt b/src/aliceVision/hdr/CMakeLists.txt
@@ -3,18 +3,21 @@
 set(hdr_files_headers
  rgbCurve.hpp
  RobertsonCalibrate.hpp
- RobertsonMerge.hpp
+ hdrMerge.hpp
  DebevecCalibrate.hpp
- DebevecMerge.hpp
+ GrossbergCalibrate.hpp
+ emorCurve.hpp
+ eiquadprog.hpp
 )
 
 # Sources
 set(hdr_files_sources
  rgbCurve.cpp
  RobertsonCalibrate.cpp
- RobertsonMerge.cpp
+ hdrMerge.cpp
  DebevecCalibrate.cpp
- DebevecMerge.cpp
+ GrossbergCalibrate.cpp
+ emorCurve.cpp
 )
 
 alicevision_add_library(aliceVision_hdr

diff --git a/src/aliceVision/hdr/DebevecCalibrate.cpp b/src/aliceVision/hdr/DebevecCalibrate.cpp
@@ -13,14 +13,15 @@
 namespace aliceVision {
 namespace hdr {
 
-//using T = Eigen::Triplet<float>;
 using T = Eigen::Triplet<double>;
 
 void DebevecCalibrate::process(const std::vector< std::vector< image::Image<image::RGBfColor> > > &ldrImageGroups,
-             const std::vector< std::vector<float> > &times, const int nbPoints,
-             const rgbCurve &weight,
-             const float lambda,
-             rgbCurve &response)
+                               const std::size_t channelQuantization,
+                               const std::vector< std::vector<float> > &times,
+                               const int nbPoints,
+                               const rgbCurve &weight,
+                               const float lambda,
+                               rgbCurve &response)
 {
   //checks
   for (int g = 0; g < ldrImageGroups.size(); ++g)
@@ -32,66 +33,54 @@ void DebevecCalibrate::process(const std::vector< std::vector< image::Image<imag
   static const std::size_t channels = 3;
 
   //initialize response
-  response = rgbCurve(_channelQuantization);
+  response = rgbCurve(channelQuantization);
   response.setLinear();
   response.normalize();
 
-
   for(unsigned int g=0; g<ldrImageGroups.size(); ++g)
   {
     const std::vector< image::Image<image::RGBfColor> > &ldrImagesGroup = ldrImageGroups.at(g);
     const std::vector<float> &ldrTimes = times.at(g);
-    const int nbCols = ldrImagesGroup.at(0).Width();
-    const int nbRows = ldrImagesGroup.at(0).Height();
     const int nbImages = ldrImagesGroup.size();
-    const int nbPixels = nbCols*nbRows;
-    const int step = std::floor(nbPixels/nbPoints);
+    const int step = std::floor(ldrImagesGroup.at(0).Width() * ldrImagesGroup.at(0).Height() / nbPoints);
 
     for(unsigned int channel=0; channel<channels; ++channel)
     {
-      sMat A(nbPoints*nbImages + _channelQuantization + 1, _channelQuantization + nbPoints);
-      Vec b = Vec::Zero(nbPoints*nbImages + _channelQuantization + 1);
+      sMat A(nbPoints*nbImages + channelQuantization + 1, channelQuantization + nbPoints);
+      Vec b = Vec::Zero(nbPoints*nbImages + channelQuantization + 1);
       int count=0;
 
       std::vector<T> tripletList;
-      tripletList.reserve(2 * nbPoints*nbImages + 1 + 3 * _channelQuantization);
+      tripletList.reserve(2 * nbPoints*nbImages + 1 + 3 * channelQuantization);
 
       std::cout << "filling A and b matrices" << std::endl;
 
       // include the data-fitting equations
       for(unsigned int j=0; j<nbImages; ++j)
       {
         const image::Image<image::RGBfColor> &image = ldrImagesGroup.at(j);
-//        for(unsigned int i=0; i<nbCols*nbRows; ++i)
         for(unsigned int i=0; i<nbPoints; ++i)
         {
-//          int index = static_cast<int>( std::round((_channelQuantization-1) * std::max(0.f, std::min(1.f, image(step*i)(channel)))) );
-          std::size_t index = static_cast<std::size_t>( std::round((_channelQuantization-1) * std::max(0.f, std::min(1.f, image(step*i)(channel)))) );
-          float w_ij = weight.getValue(index, channel);
-//          double w_ij = static_cast<double>(weight.getValue(index, channel));
-
-//          float w_ij = weight(image(step*i)(channel), channel);
-
+          float sample = std::max(0.f, std::min(1.f, image(step*i)(channel)));
+          float w_ij = weight(sample, channel);
+          std::size_t index = std::round(sample * (channelQuantization - 1));
 
           tripletList.push_back(T(count, index, w_ij));
-          tripletList.push_back(T(count, _channelQuantization+i, -w_ij));
+          tripletList.push_back(T(count, channelQuantization+i, -w_ij));
 
           b(count) = w_ij * std::log(ldrTimes.at(j));
           count += 1;
         }
       }
 
       // fix the curve by setting its middle value to zero
-      tripletList.push_back(T(count, std::floor(_channelQuantization/2), 1.f));
+      tripletList.push_back(T(count, std::floor(channelQuantization/2), 1.f));
       count += 1;
 
       // include the smoothness equations
-//      for(float k=0; k<_channelQuantization-2; ++k)
-      for(std::size_t k=0; k<_channelQuantization-2; ++k)
+      for(std::size_t k=0; k<channelQuantization-2; ++k)
       {
-//        float w = weight(static_cast<float>((k+1)/_channelQuantization), channel);
         float w = weight.getValue(k+1, channel);
-//        double w = static_cast<double>(weight.getValue(k+1, channel));
 
         tripletList.push_back(T(count, k, lambda * w));
         tripletList.push_back(T(count, k+1, -2.f * lambda * w));
@@ -102,7 +91,7 @@ void DebevecCalibrate::process(const std::vector< std::vector< image::Image<imag
 
       A.setFromTriplets(tripletList.begin(), tripletList.end());
 
-      std::cout << "solving Ax=b system" << std::endl;
+
 
       // solve the system using SVD decomposition
       A.makeCompressed();
@@ -114,9 +103,10 @@ void DebevecCalibrate::process(const std::vector< std::vector< image::Image<imag
 
       std::cout << "system solved" << std::endl;
 
-      for(std::size_t k=0; k<_channelQuantization; ++k)
-//      for(float k=0; k<_channelQuantization; ++k) {
-//        response(static_cast<float>(k/_channelQuantization), channel) = x(k);
+      double relative_error = (A*x - b).norm() / b.norm();
+      std::cout << "relative error is : " << relative_error << std::endl;
+
+      for(std::size_t k=0; k<channelQuantization; ++k)
         response.setValue(k, channel, x(k));
 
     }

diff --git a/src/aliceVision/hdr/DebevecCalibrate.hpp b/src/aliceVision/hdr/DebevecCalibrate.hpp
@@ -17,33 +17,23 @@ class DebevecCalibrate
 {
 public:  
 
-  /**
-   * @brief
-   * @param[in] maxIter
-   * @param[in] threshold
-   */
-  DebevecCalibrate(std::size_t channelQuantization = std::pow(2, 12)) :   //RAW 12 bit precision, 2^12 values between black and white point
-    _channelQuantization(channelQuantization)
-  {}
-
   /**
    * @brief
    * @param[in] groups
    * @param[in] times
+   * @param[in] nbPoints (number of samples for calibration)
    * @param[in] weight
    * @param[in] lambda (parameter of smoothness)
-   * @param[in] targetTime (target exposure time)
    * @param[out] response
    */
   void process(const std::vector< std::vector< image::Image<image::RGBfColor> > > &ldrImageGroups,
+               const std::size_t channelQuantization,
                const std::vector< std::vector<float> > &times,
                const int nbPoints,
                const rgbCurve &weight,
                const float lambda,
                rgbCurve &response);
 
-private:
-  std::size_t _channelQuantization;
 };
 
 } // namespace hdr

diff --git a/src/aliceVision/hdr/DebevecMerge.cpp b/src/aliceVision/hdr/DebevecMerge.cpp
diff --git a/src/aliceVision/hdr/DebevecMerge.hpp b/src/aliceVision/hdr/DebevecMerge.hpp