generalize number of channels

plus minor edits and fixes

modules/photo/include/opencv2/photo.hpp

@@ -591,7 +591,7 @@ public:
 @param samples number of pixel locations to use
 @param lambda smoothness term weight. Greater values produce smoother results, but can alter the
 response.
-@param random if true sample pixel locations are chosen at random, otherwise the form a
+@param random if true sample pixel locations are chosen at random, otherwise they form a
 rectangular grid.
  */
 CV_EXPORTS_W Ptr<CalibrateDebevec> createCalibrateDebevec(int samples = 70, float lambda = 10.0f, bool random = false);

modules/photo/src/calibrate.cpp

@@ -43,7 +43,6 @@
 #include "precomp.hpp"
 #include "opencv2/photo.hpp"
 #include "opencv2/imgproc.hpp"
-//#include "opencv2/highgui.hpp"
 #include "hdr_common.hpp"
 
 namespace cv
@@ -57,7 +56,7 @@ public:
         samples(_samples),
         lambda(_lambda),
         random(_random),
-        w(tringleWeights())
+        w(triangleWeights())
     {
     }
 
@@ -65,6 +64,7 @@ public:
     {
         CV_INSTRUMENT_REGION()
 
+        // check inputs
         std::vector<Mat> images;
         src.getMatVector(images);
         Mat times = _times.getMat();
@@ -72,62 +72,88 @@ public:
         CV_Assert(images.size() == times.total());
         checkImageDimensions(images);
         CV_Assert(images[0].depth() == CV_8U);
+        CV_Assert(times.type() == CV_32FC1);
 
+        // create output
         int channels = images[0].channels();
         int CV_32FCC = CV_MAKETYPE(CV_32F, channels);
+        int rows = images[0].rows;
+        int cols = images[0].cols;
 
         dst.create(LDR_SIZE, 1, CV_32FCC);
         Mat result = dst.getMat();
 
-        std::vector<Point> sample_points;
+        // pick pixel locations (either random or in a rectangular grid)
+        std::vector<Point> points;
+        points.reserve(samples);
         if(random) {
             for(int i = 0; i < samples; i++) {
-                sample_points.push_back(Point(rand() % images[0].cols, rand() % images[0].rows));
+                points.push_back(Point(rand() % cols, rand() % rows));
             }
         } else {
-            int x_points = static_cast<int>(sqrt(static_cast<double>(samples) * images[0].cols / images[0].rows));
+            int x_points = static_cast<int>(sqrt(static_cast<double>(samples) * cols / rows));
+            CV_Assert(0 < x_points && x_points <= cols);
             int y_points = samples / x_points;
-            int step_x = images[0].cols / x_points;
-            int step_y = images[0].rows / y_points;
-
+            CV_Assert(0 < y_points && y_points <= rows);
+            int step_x = cols / x_points;
+            int step_y = rows / y_points;
             for(int i = 0, x = step_x / 2; i < x_points; i++, x += step_x) {
                 for(int j = 0, y = step_y / 2; j < y_points; j++, y += step_y) {
-                    if( 0 <= x && x < images[0].cols &&
-                        0 <= y && y < images[0].rows )
-                        sample_points.push_back(Point(x, y));
+                    if( 0 <= x && x < cols && 0 <= y && y < rows ) {
+                        points.push_back(Point(x, y));
+                    }
                 }
             }
+            // we can have slightly less grid points than specified
+            //samples = static_cast<int>(points.size());
         }
 
+        // we need enough equations to ensure a sufficiently overdetermined system
+        // (maybe only as a warning)
+        //CV_Assert(points.size() * (images.size() - 1) >= LDR_SIZE);
+
+        // solve for imaging system response function, over each channel separately
         std::vector<Mat> result_split(channels);
-        for(int channel = 0; channel < channels; channel++) {
-            Mat A = Mat::zeros((int)sample_points.size() * (int)images.size() + LDR_SIZE + 1, LDR_SIZE + (int)sample_points.size(), CV_32F);
+        for(int ch = 0; ch < channels; ch++) {
+            // initialize system of linear equations
+            Mat A = Mat::zeros((int)points.size() * (int)images.size() + LDR_SIZE + 1,
+                LDR_SIZE + (int)points.size(), CV_32F);
             Mat B = Mat::zeros(A.rows, 1, CV_32F);
 
-            int eq = 0;
-            for(size_t i = 0; i < sample_points.size(); i++) {
+            // include the data−fitting equations
+            int k = 0;
+            for(size_t i = 0; i < points.size(); i++) {
                 for(size_t j = 0; j < images.size(); j++) {
-
-                    int val = images[j].ptr()[3*(sample_points[i].y * images[j].cols + sample_points[i].x) + channel];
-                    A.at<float>(eq, val) = w.at<float>(val);
-                    A.at<float>(eq, LDR_SIZE + (int)i) = -w.at<float>(val);
-                    B.at<float>(eq, 0) = w.at<float>(val) * log(times.at<float>((int)j));
-                    eq++;
+                    // val = images[j].at<Vec3b>(points[i].y, points[i].x)[ch]
+                    int val = images[j].ptr()[channels*(points[i].y * cols + points[i].x) + ch];
+                    float wij = w.at<float>(val);
+                    A.at<float>(k, val) = wij;
+                    A.at<float>(k, LDR_SIZE + (int)i) = -wij;
+                    B.at<float>(k, 0) = wij * log(times.at<float>((int)j));
+                    k++;
                 }
             }
-            A.at<float>(eq, LDR_SIZE / 2) = 1;
-            eq++;
 
-            for(int i = 0; i < 254; i++) {
-                A.at<float>(eq, i) = lambda * w.at<float>(i + 1);
-                A.at<float>(eq, i + 1) = -2 * lambda * w.at<float>(i + 1);
-                A.at<float>(eq, i + 2) = lambda * w.at<float>(i + 1);
-                eq++;
+            // fix the curve by setting its middle value to 0
+            A.at<float>(k, LDR_SIZE / 2) = 1;
+            k++;
+
+            // include the smoothness equations
+            for(int i = 0; i < (LDR_SIZE - 2); i++) {
+                float wi = w.at<float>(i + 1);
+                A.at<float>(k, i) = lambda * wi;
+                A.at<float>(k, i + 1) = -2 * lambda * wi;
+                A.at<float>(k, i + 2) = lambda * wi;
+                k++;
             }
+
+            // solve the overdetermined system using SVD (least-squares problem)
             Mat solution;
             solve(A, B, solution, DECOMP_SVD);
-            solution.rowRange(0, LDR_SIZE).copyTo(result_split[channel]);
+            solution.rowRange(0, LDR_SIZE).copyTo(result_split[ch]);
         }
+
+        // combine log-exposures and take its exponent
         merge(result_split, result);
         exp(result, result);
     }
@@ -161,11 +187,11 @@ public:
     }
 
 protected:
-    String name;
-    int samples;
-    float lambda;
-    bool random;
-    Mat w;
+    String name;  // calibration algorithm identifier
+    int samples;  // number of pixel locations to sample
+    float lambda; // constant that determines the amount of smoothness
+    bool random;  // whether to sample locations randomly or in a grid shape
+    Mat w;        // weighting function for corresponding pixel values
 };
 
 Ptr<CalibrateDebevec> createCalibrateDebevec(int samples, float lambda, bool random)

modules/photo/src/hdr_common.cpp

@@ -59,8 +59,9 @@ void checkImageDimensions(const std::vector<Mat>& images)
     }
 }
 
-Mat tringleWeights()
+Mat triangleWeights()
 {
+    // hat function
     Mat w(LDR_SIZE, 1, CV_32F);
     int half = LDR_SIZE / 2;
     for(int i = 0; i < LDR_SIZE; i++) {

modules/photo/src/hdr_common.hpp

@@ -50,7 +50,7 @@ namespace cv
 
 void checkImageDimensions(const std::vector<Mat>& images);
 
-Mat tringleWeights();
+Mat triangleWeights();
 
 void mapLuminance(Mat src, Mat dst, Mat lum, Mat new_lum, float saturation);
 

modules/photo/src/merge.cpp

@@ -52,7 +52,7 @@ class MergeDebevecImpl : public MergeDebevec
 public:
     MergeDebevecImpl() :
         name("MergeDebevec"),
-        weights(tringleWeights())
+        weights(triangleWeights())
     {
     }