New LDR2HDR software

.gitignore

@@ -37,7 +37,7 @@ x64
 CMakeFiles  # intermediate folders
 
 # user files
-*.user
+*.user*
 
 # netbeans projects
 /nbproject

src/CMakeLists.txt

@@ -50,6 +50,7 @@ set(CMAKE_CXX_STANDARD_REQUIRED ON)
 option(ALICEVISION_BUILD_SHARED "Build AliceVision shared libs" ON)
 option(ALICEVISION_BUILD_SFM "Build AliceVision SfM part" ON)
 option(ALICEVISION_BUILD_MVS "Build AliceVision MVS part" ON)
+option(ALICEVISION_BUILD_HDR "Build AliceVision HDR part" ON)
 option(ALICEVISION_BUILD_EXAMPLES "Build AliceVision samples applications." OFF)
 option(ALICEVISION_BUILD_COVERAGE "Enable code coverage generation (gcc only)" OFF)
 trilean_option(ALICEVISION_BUILD_DOC "Build AliceVision documentation" AUTO)

src/aliceVision/CMakeLists.txt

@@ -1,7 +1,10 @@
 # Common modules
 
+add_subdirectory(image)
+add_subdirectory(numeric)
 add_subdirectory(system)
 
+
 # SfM modules
 
 if(ALICEVISION_BUILD_SFM)
@@ -12,14 +15,12 @@ if(ALICEVISION_BUILD_SFM)
   add_subdirectory(geometry)
   add_subdirectory(graph)
   add_subdirectory(gpu)
-  add_subdirectory(image)
   add_subdirectory(keyframe)
   add_subdirectory(linearProgramming)
   add_subdirectory(localization)
   add_subdirectory(matching)
   add_subdirectory(matchingImageCollection)
   add_subdirectory(multiview)
-  add_subdirectory(numeric)
   add_subdirectory(rig)
   add_subdirectory(robustEstimation)
   add_subdirectory(sensorDB)
@@ -35,6 +36,10 @@ if(ALICEVISION_BUILD_SFM)
   endif()
 endif()
 
+if(ALICEVISION_BUILD_HDR)
+add_subdirectory(hdr)
+endif()
+
 # MVS modules
 
 if(ALICEVISION_BUILD_MVS)

src/aliceVision/hdr/CMakeLists.txt

@@ -0,0 +1,31 @@
+
+# Headers
+set(hdr_files_headers
+ rgbCurve.hpp
+ RobertsonCalibrate.hpp
+ hdrMerge.hpp
+ DebevecCalibrate.hpp
+ GrossbergCalibrate.hpp
+ emorCurve.hpp
+)
+
+# Sources
+set(hdr_files_sources
+ rgbCurve.cpp
+ RobertsonCalibrate.cpp
+ hdrMerge.cpp
+ DebevecCalibrate.cpp
+ GrossbergCalibrate.cpp
+ emorCurve.cpp
+)
+
+alicevision_add_library(aliceVision_hdr
+  SOURCES ${hdr_files_headers} ${hdr_files_sources}
+  PRIVATE_LINKS
+    aliceVision_system
+    aliceVision_image
+    ${Boost_FILESYSTEM_LIBRARY}
+)
+
+# Unit tests
+

src/aliceVision/hdr/DebevecCalibrate.cpp

@@ -0,0 +1,118 @@
+// This file is part of the AliceVision project.
+// Copyright (c) 2019 AliceVision contributors.
+// This Source Code Form is subject to the terms of the Mozilla Public License,
+// v. 2.0. If a copy of the MPL was not distributed with this file,
+// You can obtain one at https://mozilla.org/MPL/2.0/.
+
+#include "DebevecCalibrate.hpp"
+#include <iostream>
+#include <fstream>
+#include <cassert>
+#include <aliceVision/alicevision_omp.hpp>
+#include <aliceVision/system/Logger.hpp>
+
+namespace aliceVision {
+namespace hdr {
+
+using T = Eigen::Triplet<double>;
+
+void DebevecCalibrate::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)
+{
+  //checks
+  for (int g = 0; g < ldrImageGroups.size(); ++g)
+  {
+    assert(ldrImageGroups[g].size() == times[g].size());
+  }
+
+  //set channels count always RGB
+  static const std::size_t channels = 3;
+
+  //initialize response
+  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 nbImages = ldrImagesGroup.size();
+    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);
+      int count=0;
+
+      std::vector<T> tripletList;
+      tripletList.reserve(2 * nbPoints*nbImages + 1 + 3 * channelQuantization);
+
+//      ALICEVISION_LOG_TRACE("filling A and b matrices");
+
+      // 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<nbPoints; ++i)
+        {
+          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));
+
+          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));
+      count += 1;
+
+      // include the smoothness equations
+      for(std::size_t k=0; k<channelQuantization-2; ++k)
+      {
+        float w = weight.getValue(k+1, channel);
+
+        tripletList.push_back(T(count, k, lambda * w));
+        tripletList.push_back(T(count, k+1, -2.f * lambda * w));
+        tripletList.push_back(T(count, k+2, lambda * w));
+
+        count += 1;
+      }
+
+      A.setFromTriplets(tripletList.begin(), tripletList.end());
+
+
+
+      // solve the system using SVD decomposition
+      A.makeCompressed();
+      Eigen::SparseQR<sMat, Eigen::COLAMDOrdering<int>> solver;
+      solver.compute(A);
+      if(solver.info() != Eigen::Success)  return; // decomposition failed
+      Vec x = solver.solve(b);
+      if(solver.info() != Eigen::Success)  return; // solving failed
+
+//      ALICEVISION_LOG_TRACE("system solved");
+
+//      double relative_error = (A*x - b).norm() / b.norm();
+//      ALICEVISION_LOG_TRACE("relative error is : " << relative_error);
+
+      for(std::size_t k=0; k<channelQuantization; ++k)
+        response.setValue(k, channel, x(k));
+
+    }
+  }
+}
+
+} // namespace hdr
+} // namespace aliceVision

src/aliceVision/hdr/DebevecCalibrate.hpp

@@ -0,0 +1,41 @@
+// This file is part of the AliceVision project.
+// Copyright (c) 2019 AliceVision contributors.
+// This Source Code Form is subject to the terms of the Mozilla Public License,
+// v. 2.0. If a copy of the MPL was not distributed with this file,
+// You can obtain one at https://mozilla.org/MPL/2.0/.
+
+#pragma once
+#include <aliceVision/image/all.hpp>
+#include "rgbCurve.hpp"
+#include <aliceVision/numeric/numeric.hpp>
+#include <Eigen/SparseQR>
+
+namespace aliceVision {
+namespace hdr {
+
+class DebevecCalibrate
+{
+public:  
+
+  /**
+   * @brief
+   * @param[in] LDR images groups
+   * @param[in] channel quantization
+   * @param[in] exposure times
+   * @param[in] number of samples
+   * @param[in] calibration weight function
+   * @param[in] lambda (parameter of smoothness)
+   * @param[out] camera response function
+   */
+  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);
+
+};
+
+} // namespace hdr
+} // namespace aliceVision

src/aliceVision/hdr/GrossbergCalibrate.cpp

@@ -0,0 +1,138 @@
+// This file is part of the AliceVision project.
+// Copyright (c) 2019 AliceVision contributors.
+// This Source Code Form is subject to the terms of the Mozilla Public License,
+// v. 2.0. If a copy of the MPL was not distributed with this file,
+// You can obtain one at https://mozilla.org/MPL/2.0/.
+
+#include "GrossbergCalibrate.hpp"
+#include <iostream>
+#include <cassert>
+#include <aliceVision/alicevision_omp.hpp>
+#include <aliceVision/system/Logger.hpp>
+#include <Eigen/Dense>
+
+
+namespace aliceVision {
+namespace hdr {
+
+
+GrossbergCalibrate::GrossbergCalibrate(const unsigned int dimension)
+{
+    _dimension = dimension;
+}
+
+void GrossbergCalibrate::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,
+                                 rgbCurve &response)
+{
+
+    //checks
+    for (int g = 0; g < ldrImageGroups.size(); ++g)
+    {
+      assert(ldrImageGroups[g].size() == times[g].size());
+    }
+
+    //set channels count always RGB
+    static const std::size_t channels = 3;
+
+    //initialize response with g0 from invEmor
+    response = rgbCurve(channelQuantization);
+    const double* ptrf0 = getEmorInvCurve(0);
+    std::vector<double> f0;
+    f0.assign(ptrf0, ptrf0 + channelQuantization);
+
+    Mat H(channelQuantization, _dimension);
+    std::vector<std::vector<double>> hCurves(_dimension);
+
+    for(unsigned int i=0; i<_dimension; ++i)
+    {
+        const double *h = getEmorInvCurve(i+1);
+        hCurves[i].assign(h, h + channelQuantization);
+        H.col(i) = Eigen::Map<Vec>(hCurves[i].data(), channelQuantization);
+    }
+
+    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 nbImages = ldrImagesGroup.size();
+        const int nbPixels = ldrImagesGroup.at(0).Width() * ldrImagesGroup.at(0).Height();
+        const int step = std::floor(nbPixels/nbPoints);
+
+
+        for(unsigned int channel=0; channel<channels; ++channel)
+        {
+            response.getCurve(channel) = std::vector<float>(f0.begin(), f0.end());
+
+            Mat A = Mat::Zero(nbPoints*(nbImages-1), _dimension);
+            Vec b = Vec::Zero(nbPoints*(nbImages-1));
+
+//            ALICEVISION_LOG_TRACE("filling A and b matrices");
+
+            for(unsigned int j=0; j<nbImages-1; ++j)
+            {
+                const image::Image<image::RGBfColor> &image1 = ldrImagesGroup.at(j);
+                const image::Image<image::RGBfColor> &image2 = ldrImagesGroup.at(j+1);
+                const double k = ldrTimes.at(j+1)/ldrTimes.at(j);
+
+                // fill A and b matrices with the equations
+                for(unsigned int l=0; l<nbPoints; ++l)
+                {
+                    double sample1 = std::max(0.f, std::min(1.f, image1(step*l)(channel)));
+                    double sample2 = std::max(0.f, std::min(1.f, image2(step*l)(channel)));
+
+
+                    std::size_t index1 = std::round((channelQuantization-1) * sample1);
+                    std::size_t index2 = std::round((channelQuantization-1) * sample2);
+
+//                    double w = std::min(weight(sample1, channel), weight(sample2, channel));
+
+                    double w = (weight(sample1, channel) + weight(sample2, channel)) / 2.0;
+
+//                    double w_R = (weight(sample1, 0) + weight(sample2, 0)) / 2.0;
+//                    double w_G = (weight(sample1, 1) + weight(sample2, 1)) / 2.0;
+//                    double w_B = (weight(sample1, 2) + weight(sample2, 2)) / 2.0;
+
+//                    double w = (w_R + w_G + w_B) / 3.0;
+//                    double w = (w_R * w_G * w_B);
+//                    double w = std::min(std::min(w_R, w_G), w_B);
+
+                    b(j*nbPoints + l) = w * (f0.at(index2) - k * f0.at(index1));
+                    for(unsigned int i=0; i<_dimension; ++i)
+                      A(j*nbPoints + l, i) = w * (k * H(index1, i) - H(index2, i));
+
+                }
+            }
+
+//            ALICEVISION_LOG_TRACE("solving Ax=b system");
+
+            // solve the system using QR decomposition
+            Eigen::HouseholderQR<Mat> solver(A);
+            Vec c = solver.solve(b);
+
+//            ALICEVISION_LOG_TRACE("system solved");
+
+//            double relative_error = (A*c - b).norm() / b.norm();
+//            ALICEVISION_LOG_TRACE("relative error is : " << relative_error);
+
+            for(unsigned int i=0; i<_dimension; ++i)
+            {
+              std::vector<double> temp_hCurve = hCurves[i];
+              for(auto &value : temp_hCurve)
+                value *= c(i);
+
+//              ALICEVISION_LOG_TRACE(c(i));
+
+              std::transform(response.getCurve(channel).begin(), response.getCurve(channel).end(), temp_hCurve.begin(), response.getCurve(channel).begin(), std::plus<float>());
+            }
+        }
+    }
+}
+
+
+} // namespace hdr
+} // namespace aliceVision