Change weighting fuction, fix linear_least_square_solver, refactor va…

…riable names

HDRI/debevec_weight.cpp

@@ -3,31 +3,26 @@
 namespace {
 constexpr int z_min = 0;
 constexpr int z_max = 255;
-constexpr int range = z_max - z_min + 1;
+constexpr int size = z_max - z_min + 1;
+constexpr auto z_mid = static_cast<int>(0.5 * (z_max + z_min));
 }  // namespace
 
 namespace HDRI {
 
-debevec_weight::debevec_weight() {
-  table.resize(range);
-  for (auto i = 0; i < range; ++i) {
-    if (i <= (0.5 * (z_max + z_min))) {
-      table[i] = static_cast<std::uint8_t>(i - z_min);
-    } else {
-      table[i] = static_cast<std::uint8_t>(z_max - i);
-    }
-  }
-}
-
-auto debevec_weight::get_size() const noexcept -> size_t { return range; }
+auto debevec_weight::get_size() noexcept -> std::size_t { return size; }
 
-auto debevec_weight::get_weight(int index) const noexcept -> double {
+auto debevec_weight::get_weight(const int z_value) noexcept -> double {
   // cap ?
 
-  if (index < z_min || index > z_max) {
-    return 0;
+  if (z_value < z_min || z_value > z_max) {
+    return 0.0;
   }
-  return table[index];
+
+  if (z_value <= z_mid) {
+    return z_value - z_min;
+  }
+
+  return z_max - z_value;
 }
 
 }  // namespace HDRI

HDRI/debevec_weight.hpp

@@ -1,17 +1,14 @@
 #pragma once
 
-#include <vector>
+#include <cstddef>
 
 namespace HDRI {
 
 class debevec_weight final {
  public:
-  debevec_weight();
-  [[nodiscard]] double get_weight(const int index) const noexcept;
-  [[nodiscard]] std::size_t get_size() const noexcept;
-
- private:
-  std::vector<std::uint8_t> table;
+  debevec_weight() = delete;
+  [[nodiscard]] static double get_weight(const int index) noexcept;
+  [[nodiscard]] static std::size_t get_size() noexcept;
 };
 
 }  // namespace HDRI

HDRI/hdr_image.cpp

@@ -53,7 +53,7 @@ std::vector<cv::Mat> shrink_images(const std::vector<HDRI::raw_image> &input) no
   std::vector<cv::Mat> out;
   out.reserve(std::size(input));
 
-  constexpr size_t ratio = 100UL;
+  constexpr size_t ratio = 50UL;
 
   for (const auto &img : input) {
     const auto &image_data = img.get_image_data();
@@ -108,13 +108,13 @@ std::array<std::vector<std::vector<int>>, 3UL> convert_to_z(const std::vector<st
     z_colors.resize(image_size);
   }
 
-  for (size_t i = 0; i < image_size; ++i) {  // image pixel
+  for (auto i = 0UL; i < image_size; ++i) {  // number of pixels in an image
 
     z_values[0UL][i].resize(num_images);
     z_values[1UL][i].resize(num_images);
     z_values[2UL][i].resize(num_images);
 
-    for (size_t j = 0; j < num_images; ++j) {  // num of iamge
+    for (auto j = 0UL; j < num_images; ++j) {  // number of images
       z_values[0UL][i][j] = pixel[j][i].b;
       z_values[1UL][i][j] = pixel[j][i].g;
       z_values[2UL][i][j] = pixel[j][i].r;
@@ -151,15 +151,14 @@ void hdr_image::compute_curves(const std::vector<HDRI::raw_image> &raw_images,
   // convert
   const auto z_values = convert_to_z(raw_pixel, shrink_mat[0UL].total(), shrink_mat.size());
 
-  debevec_weight dwf;
   constexpr auto lambda = 10;
 
-  const auto weight_function = [&dwf](const auto color_index) { return dwf.get_weight(color_index); };
+  const auto weight_function = [](const auto color) { return debevec_weight::get_weight(color); };
 
   std::array<std::future<cv::Mat>, 3UL> futures;
   for (auto c = 0U; c < std::size(z_values); ++c) {
     futures[c] = std::async(std::launch::async, linear_least_square_solver<decltype(weight_function)>, weight_function,
-                            z_values[c], exposure, dwf.get_size(), lambda);
+                            z_values[c], exposure, debevec_weight::get_size(), lambda);
   }
 
   for (auto c = 0U; c < std::size(futures); ++c) {
@@ -169,9 +168,7 @@ void hdr_image::compute_curves(const std::vector<HDRI::raw_image> &raw_images,
 
 void hdr_image::compute_radiance(const std::vector<HDRI::raw_image> &raw_images,
                                  const std::vector<double> &exposure) noexcept {
-  debevec_weight dwf;
-  radiance = construct_radiance(
-      raw_images, curves, [&dwf](auto color_index) { return dwf.get_weight(color_index); }, exposure);
+  radiance = construct_radiance(raw_images, curves, debevec_weight::get_weight, exposure);
 }
 
 const cv::Mat &hdr_image::get_radiance() const noexcept { return radiance; }

HDRI/linear_least_square_solver.hpp

@@ -15,45 +15,48 @@ template <typename weight_function>
   // From the paper.
 
   /// NOTE: Ax = b
+  /// z_value[i][j]: the pixel value of location i in image j
 
-  const auto n = range;
+  const auto num_plxels{std::size(z_value)};
+  const auto num_images{std::size(z_value[0UL])};
 
-  cv::Mat A = cv::Mat::zeros(z_value.size() * z_value[0UL].size() + n + 1, n + z_value.size(), CV_64F);
+  cv::Mat A = cv::Mat::zeros(num_plxels * num_images + 1UL + (range - 1UL - 1UL), range + num_plxels, CV_64F);
   cv::Mat b = cv::Mat::zeros(A.size().height, 1, CV_64F);
 
-  int k = 0;
-  for (auto i = 0; i < z_value.size(); ++i) {       // N  (pixel)
-    for (auto j = 0; j < z_value[i].size(); ++j) {  // P  (image)
+  size_t current_row = 0UL;
+  for (auto i = 0; i < num_plxels; ++i) {    // N  (pixel)
+    for (auto j = 0; j < num_images; ++j) {  // P  (image)
 
       const auto w_ij = std::invoke(std::forward<weight_function>(weighting), z_value[i][j]);
-      A.at<double>(k, z_value[i][j]) = w_ij;
-      A.at<double>(k, n + i) = -w_ij;
-      b.at<double>(k, 0) = w_ij * std::log(delta_time[j]);
+      A.at<double>(current_row, z_value[i][j]) = w_ij;
+      A.at<double>(current_row, range + i) = -w_ij;
+      b.at<double>(current_row, 0) = w_ij * std::log(delta_time[j]);
 
-      ++k;
+      ++current_row;
     }
   }
 
-  A.at<double>(k, 128) = 1;  // set the middle to Zero
-  ++k;
+  A.at<double>(current_row, 127) = 1;  // set the middle to Zero: g(127) = 0
+  ++current_row;
 
   // get x !
 
-  // g(z - 1) -2g(z) + g(z + 1)
-  for (auto i = 0; i < n - 1; ++i) {
-    const auto result = std::invoke(std::forward<weight_function>(weighting), i + 1);
+  // w(z)g"(z) = w(z)(g(z - 1) -2g(z) + g(z + 1))
+  for (auto z = 1UL; z <= range - 1UL - 1UL; ++z) {
+    const auto result = std::invoke(std::forward<weight_function>(weighting), z);
 
-    A.at<double>(k, i) = lambda * result;
-    A.at<double>(k, i + 1) = -2 * lambda * result;
-    A.at<double>(k, i + 2) = lambda * result;
+    A.at<double>(current_row, z - 1) = lambda * result;
+    A.at<double>(current_row, z) = -2.0 * lambda * result;
+    A.at<double>(current_row, z + 1) = lambda * result;
 
-    ++k;
+    ++current_row;
   }
 
   cv::Mat x;
   cv::solve(A, b, x, cv::DECOMP_SVD);
 
   // get g and lE
+  // g(z) : the log exposure corresponding to pixel value z
 
   return x;
 }