Robertson update

modules/photo/include/opencv2/photo.hpp

@@ -80,6 +80,8 @@ CV_EXPORTS_W void fastNlMeansDenoisingColoredMulti( InputArrayOfArrays srcImgs,
                                                     float h = 3, float hColor = 3,
                                                     int templateWindowSize = 7, int searchWindowSize = 21);
 
+enum { LDR_SIZE = 256 };
+
 class CV_EXPORTS_W Tonemap : public Algorithm
 {
 public:
@@ -227,9 +229,11 @@ public:
     
     CV_WRAP virtual float getThreshold() const = 0;
     CV_WRAP virtual void setThreshold(float threshold) = 0;
+    
+    CV_WRAP virtual Mat getRadiance() const = 0;
 };
 
-CV_EXPORTS_W Ptr<CalibrateRobertson> createCalibrateRobertson(int samples = 50, float lambda = 10.0f);
+CV_EXPORTS_W Ptr<CalibrateRobertson> createCalibrateRobertson(int max_iter = 30, float threshold = 0.01f);
 
 class CV_EXPORTS_W ExposureMerge : public Algorithm
 {

modules/photo/src/align.cpp

@@ -243,14 +243,14 @@ protected:
     {
         int channels = 0;
         Mat hist; 
-        int hist_size = 256;
+        int hist_size = LDR_SIZE;
-        float range[] = {0, 256} ;
+        float range[] = {0, LDR_SIZE} ;
         const float* ranges[] = {range};
         calcHist(&img, 1, &channels, Mat(), hist, 1, &hist_size, ranges);
         float *ptr = hist.ptr<float>();
         int median = 0, sum = 0;
         int thresh = img.total() / 2;
-        while(sum < thresh && median < 256) {
+        while(sum < thresh && median < LDR_SIZE) {
             sum += static_cast<int>(ptr[median]);
             median++;
         }
@@ -309,7 +309,7 @@ public:
 
         std::vector<Mat> splitted(channels);
         split(images[0], splitted);
-        for(int i = 0; i < images.size() - 1; i++) {
+        for(size_t i = 0; i < images.size() - 1; i++) {
             
             std::vector<Mat> next_splitted(channels);
             split(images[i + 1], next_splitted);
@@ -399,7 +399,7 @@ public:
         split(radiance, splitted);
         std::vector<Mat> resp_split(channels);
         split(response, resp_split);
-        for(int i = 0; i < images.size() - 1; i++) {
+        for(size_t i = 0; i < images.size() - 1; i++) {
             
             std::vector<Mat> next_splitted(channels);
             LUT(images[i + 1], response, radiance);
@@ -430,7 +430,9 @@ public:
 
     virtual void process(InputArrayOfArrays src, OutputArray dst, std::vector<float>& times)
     {
-        process(src, dst, times, linearResponse(3));
+        Mat response = linearResponse(3);
+        response.at<Vec3f>(0) = response.at<Vec3f>(1);
+        process(src, dst, times, response);
     }
 
     CV_WRAP virtual int getThreshold() {return thresh;}

modules/photo/src/calibrate.cpp

@@ -45,7 +45,6 @@
 #include "opencv2/imgproc.hpp"
 //#include "opencv2/highgui.hpp"
 #include "hdr_common.hpp"
-#include <iostream>
 
 namespace cv
 {
@@ -74,7 +73,7 @@ public:
         int channels = images[0].channels();
         int CV_32FCC = CV_MAKETYPE(CV_32F, channels);
 
-        dst.create(256, 1, CV_32FCC);
+        dst.create(LDR_SIZE, 1, CV_32FCC);
         Mat result = dst.getMat();
         
         std::vector<Point> sample_points;
@@ -97,7 +96,7 @@ public:
 
         std::vector<Mat> result_split(channels);
         for(int channel = 0; channel < channels; channel++) {
-            Mat A = Mat::zeros(sample_points.size() * images.size() + 257, 256 + sample_points.size(), CV_32F);
+            Mat A = Mat::zeros(sample_points.size() * images.size() + LDR_SIZE + 1, LDR_SIZE + sample_points.size(), CV_32F);
             Mat B = Mat::zeros(A.rows, 1, CV_32F);
 
             int eq = 0;
@@ -107,12 +106,12 @@ public:
 
                     int val = images[j].ptr()[3*(sample_points[i].y * images[j].cols + sample_points[j].x) + channel];
                     A.at<float>(eq, val) = w.at<float>(val);
-                    A.at<float>(eq, 256 + i) = -w.at<float>(val);
+                    A.at<float>(eq, LDR_SIZE + i) = -w.at<float>(val);
                     B.at<float>(eq, 0) = w.at<float>(val) * log(times[j]);        
                     eq++;
                 }
             }
-            A.at<float>(eq, 128) = 1;
+            A.at<float>(eq, LDR_SIZE / 2) = 1;
             eq++;
 
             for(int i = 0; i < 254; i++) {
@@ -123,7 +122,7 @@ public:
             }
             Mat solution;
             solve(A, B, solution, DECOMP_SVD);
-            solution.rowRange(0, 256).copyTo(result_split[channel]);
+            solution.rowRange(0, LDR_SIZE).copyTo(result_split[channel]);
         }
         merge(result_split, result);
         exp(result, result);
@@ -192,20 +191,14 @@ public:
         int channels = images[0].channels();
         int CV_32FCC = CV_MAKETYPE(CV_32F, channels);
 
-        dst.create(256, 1, CV_32FCC);
+        dst.create(LDR_SIZE, 1, CV_32FCC);
         Mat response = dst.getMat();
-        
+        response = linearResponse(3) / (LDR_SIZE / 2.0f);
-        response = Mat::zeros(256, 1, CV_32FCC);
-        for(int i = 0; i < 256; i++) {
-            for(int c = 0; c < channels; c++) {
-                response.at<Vec3f>(i)[c] = i / 128.0;
-            }
-        }
 
-        Mat card = Mat::zeros(256, 1, CV_32FCC);
+        Mat card = Mat::zeros(LDR_SIZE, 1, CV_32FCC);
-        for(int i = 0; i < images.size(); i++) {
+        for(size_t i = 0; i < images.size(); i++) {
            uchar *ptr = images[i].ptr();
-           for(int pos = 0; pos < images[i].total(); pos++) {
+           for(size_t pos = 0; pos < images[i].total(); pos++) {
                for(int c = 0; c < channels; c++, ptr++) {
                    card.at<Vec3f>(*ptr)[c] += 1;
                }
@@ -213,43 +206,34 @@ public:
         }
         card = 1.0 / card;
 
+        Ptr<MergeRobertson> merge = createMergeRobertson();        
         for(int iter = 0; iter < max_iter; iter++) {
 
-            Scalar channel_err(0, 0, 0);
+            radiance = Mat::zeros(images[0].size(), CV_32FCC);
-            Mat radiance = Mat::zeros(images[0].size(), CV_32FCC);
+            merge->process(images, radiance, times, response);
-            Mat wsum = Mat::zeros(images[0].size(), CV_32FCC);
-            for(int i = 0; i < images.size(); i++) {
-                Mat im, w;
-                LUT(images[i], weight, w);
-                LUT(images[i], response, im);
 
-                Mat err_mat;
+            Mat new_response = Mat::zeros(LDR_SIZE, 1, CV_32FC3);
-                pow(im - times[i] * radiance, 2.0f, err_mat);
+            for(size_t i = 0; i < images.size(); i++) {
-                err_mat = w.mul(err_mat);
-                channel_err += sum(err_mat);
-
-                radiance += times[i] * w.mul(im);
-                wsum += pow(times[i], 2) * w;
-            }
-            float err = (channel_err[0] + channel_err[1] + channel_err[2]) / (channels * radiance.total());
-            radiance = radiance.mul(1 / wsum);
-
-            float* rad_ptr = radiance.ptr<float>();
-            response = Mat::zeros(256, 1, CV_32FC3);
-            for(int i = 0; i < images.size(); i++) {
                 uchar *ptr = images[i].ptr();
-                for(int pos = 0; pos < images[i].total(); pos++) {
+                float* rad_ptr = radiance.ptr<float>();
+                for(size_t pos = 0; pos < images[i].total(); pos++) {
                     for(int c = 0; c < channels; c++, ptr++, rad_ptr++) {
-                        response.at<Vec3f>(*ptr)[c] += times[i] * *rad_ptr;
+                        new_response.at<Vec3f>(*ptr)[c] += times[i] * *rad_ptr;
                     }
                 }
             }
-            response = response.mul(card);
+            new_response = new_response.mul(card);
             for(int c = 0; c < 3; c++) {
-                for(int i = 0; i < 256; i++) {
+                float middle = new_response.at<Vec3f>(LDR_SIZE / 2)[c];
-                    response.at<Vec3f>(i)[c] /= response.at<Vec3f>(128)[c];
+                for(int i = 0; i < LDR_SIZE; i++) {
+                    new_response.at<Vec3f>(i)[c] /= middle;
                 }
             }
+            float diff = sum(sum(abs(new_response - response)))[0] / channels;
+            new_response.copyTo(response);
+            if(diff < threshold) {
+                break;
+            }
         }
     }
 
@@ -259,6 +243,8 @@ public:
     float getThreshold() const { return threshold; }
     void setThreshold(float val) { threshold = val; }
 
+    Mat getRadiance() const { return radiance; }
+
     void write(FileStorage& fs) const
     {
         fs << "name" << name
@@ -278,7 +264,7 @@ protected:
     String name;
     int max_iter;
     float threshold;
-    Mat weight;
+    Mat weight, radiance;
 };
 
 Ptr<CalibrateRobertson> createCalibrateRobertson(int max_iter, float threshold)

modules/photo/src/hdr_common.cpp

@@ -61,21 +61,22 @@ void checkImageDimensions(const std::vector<Mat>& images)
 
 Mat tringleWeights()
 {
-    Mat w(256, 1, CV_32F);
+    Mat w(LDR_SIZE, 1, CV_32F);
-    for(int i = 0; i < 256; i++) {
+    int half = LDR_SIZE / 2;
-        w.at<float>(i) = i < 128 ? i + 1.0f : 256.0f - i;
+    for(int i = 0; i < LDR_SIZE; i++) {
+        w.at<float>(i) = i < half ? i + 1.0f : LDR_SIZE - i;
     }
     return w;
 }
 
 Mat RobertsonWeights()
 {
-    Mat weight(256, 1, CV_32FC3);
+    Mat weight(LDR_SIZE, 1, CV_32FC3);
-    for(int i = 0; i < 256; i++) {
+    float q = (LDR_SIZE - 1) / 4.0f;
-        float value = exp(-4.0f * pow(i - 127.5f, 2.0f) / pow(127.5f, 2.0f));
+    for(int i = 0; i < LDR_SIZE; i++) {
-        for(int c = 0; c < 3; c++) {
+        float value = i / q - 2.0f;
-            weight.at<Vec3f>(i)[c] = value;
+        value = exp(-value * value);
-        }
+        weight.at<Vec3f>(i) = Vec3f::all(value);
     }
     return weight;
 }
@@ -94,19 +95,11 @@ void mapLuminance(Mat src, Mat dst, Mat lum, Mat new_lum, float saturation)
 
 Mat linearResponse(int channels)
 {
-    Mat single_response = Mat(256, 1, CV_32F);
+    Mat response = Mat(LDR_SIZE, 1, CV_MAKETYPE(CV_32F, channels));
-    for(int i = 1; i < 256; i++) {
+    for(int i = 0; i < LDR_SIZE; i++) {
-        single_response.at<float>(i) = static_cast<float>(i);
+        response.at<Vec3f>(i) = Vec3f::all(i);
     }
-    single_response.at<float>(0) = static_cast<float>(1);
+    return response;
-
-    std::vector<Mat> splitted(channels);
-    for(int c = 0; c < channels; c++) {
-        splitted[c] = single_response;
-    }
-    Mat result;
-    merge(splitted, result);
-    return result;
 }
 
 };

modules/photo/src/merge.cpp

@@ -43,7 +43,6 @@
 #include "opencv2/photo.hpp"
 #include "opencv2/imgproc.hpp"
 #include "hdr_common.hpp"
-#include <iostream>
 
 namespace cv
 {
@@ -77,9 +76,10 @@ public:
 
         if(response.empty()) {
             response = linearResponse(channels);
+            response.at<Vec3f>(0) = response.at<Vec3f>(1);
         }
         log(response, response);
-        CV_Assert(response.rows == 256 && response.cols == 1 && 
+        CV_Assert(response.rows == LDR_SIZE && response.cols == 1 && 
                   response.channels() == channels);
 
         Mat exp_values(times);
@@ -312,9 +312,9 @@ public:
 
         Mat response = input_response.getMat();
         if(response.empty()) {
-            response = linearResponse(channels) / 128.0f;
+            response = linearResponse(channels) / (LDR_SIZE / 2.0f);
         }
-        CV_Assert(response.rows == 256 && response.cols == 1 && 
+        CV_Assert(response.rows == LDR_SIZE && response.cols == 1 && 
                   response.channels() == channels);
     
         result = Mat::zeros(images[0].size(), CV_32FCC);

modules/photo/test/test_hdr.cpp

@@ -187,7 +187,22 @@ TEST(Photo_MergeDebevec, regression)
 	Mat result, expected;
 	loadImage(test_path + "merge/debevec.exr", expected);
 	merge->process(images, result, times, response);
-	imwrite("test.exr", result);
+	checkEqual(expected, result, 1e-2f);
+}
+
+TEST(Photo_MergeRobertson, regression)
+{
+	string test_path = string(cvtest::TS::ptr()->get_data_path()) + "hdr/";
+
+	vector<Mat> images;
+	vector<float> times;
+	loadExposureSeq(test_path + "exposures/", images, times);
+
+	Ptr<MergeRobertson> merge = createMergeRobertson();
+
+	Mat result, expected;
+	loadImage(test_path + "merge/robertson.exr", expected);
+	merge->process(images, result, times);
 	checkEqual(expected, result, 1e-2f);
 }
 
@@ -208,3 +223,18 @@ TEST(Photo_CalibrateDebevec, regression)
     minMaxLoc(diff, NULL, &max);
     ASSERT_FALSE(max > 0.1);
 }
+
+TEST(Photo_CalibrateRobertson, regression)
+{
+	string test_path = string(cvtest::TS::ptr()->get_data_path()) + "hdr/";
+
+	vector<Mat> images;
+	vector<float> times;
+	Mat response, expected;
+	loadExposureSeq(test_path + "exposures/", images, times);
+    loadResponseCSV(test_path + "calibrate/robertson.csv", expected);
+
+	Ptr<CalibrateRobertson> calibrate = createCalibrateRobertson();
+	calibrate->process(images, response, times);
+    checkEqual(expected, response, 1e-3f);
+}