Mantiuk's tonemapping

modules/highgui/src/grfmt_hdr.cpp

@@ -123,10 +123,9 @@ HdrEncoder::~HdrEncoder()
 
 bool HdrEncoder::write( const Mat& _img, const std::vector<int>& params )
 {
+	CV_Assert(_img.channels() == 3);
 	Mat img;
-	if(_img.depth() == CV_32F) {
-		_img.convertTo(img, CV_32FC3);
-	} else {
+	if(_img.depth() != CV_32F) {
 		_img.convertTo(img, CV_32FC3, 1/255.0f);
 	}
 	CV_Assert(params.empty() || params[0] == HDR_NONE || params[0] == HDR_RLE);

modules/photo/include/opencv2/photo.hpp

@@ -102,17 +102,25 @@ CV_EXPORTS_W Ptr<TonemapLinear> createTonemapLinear(float gamma = 1.0f);
 class CV_EXPORTS_W TonemapDrago : public Tonemap
 {
 public:
+
+    CV_WRAP virtual float getSaturation() const = 0;
+    CV_WRAP virtual void setSaturation(float saturation) = 0;
+
     CV_WRAP virtual float getBias() const = 0;
     CV_WRAP virtual void setBias(float bias) = 0;
 };
 
-CV_EXPORTS_W Ptr<TonemapDrago> createTonemapDrago(float gamma = 1.0f, float bias = 0.85f);
+CV_EXPORTS_W Ptr<TonemapDrago> createTonemapDrago(float gamma = 1.0f, float saturation = 1.0f, float bias = 0.85f);
 
 // "Fast Bilateral Filtering for the Display of High-Dynamic-Range Images", Durand, Dorsey, 2002
 
 class CV_EXPORTS_W TonemapDurand : public Tonemap
 {
 public:
+
+    CV_WRAP virtual float getSaturation() const = 0;
+    CV_WRAP virtual void setSaturation(float saturation) = 0;
+
     CV_WRAP virtual float getContrast() const = 0;
     CV_WRAP virtual void setContrast(float contrast) = 0;
 
@@ -124,7 +132,7 @@ public:
 };
 
 CV_EXPORTS_W Ptr<TonemapDurand> 
-createTonemapDurand(float gamma = 1.0f, float contrast = 4.0f, float sigma_space = 2.0f, float sigma_color = 2.0f);
+createTonemapDurand(float gamma = 1.0f, float saturation = 1.0f, float contrast = 4.0f, float sigma_space = 2.0f, float sigma_color = 2.0f);
 
 // "Dynamic Range Reduction Inspired by Photoreceptor Physiology", Reinhard, Devlin, 2005
 
@@ -144,6 +152,19 @@ public:
 CV_EXPORTS_W Ptr<TonemapReinhardDevlin> 
 createTonemapReinhardDevlin(float gamma = 1.0f, float intensity = 0.0f, float light_adapt = 1.0f, float color_adapt = 0.0f);
 
+class CV_EXPORTS_W TonemapMantiuk : public Tonemap
+{
+public:
+    CV_WRAP virtual float getScale() const = 0;
+    CV_WRAP virtual void setScale(float scale) = 0;
+
+    CV_WRAP virtual float getSaturation() const = 0;
+    CV_WRAP virtual void setSaturation(float saturation) = 0;
+};
+
+CV_EXPORTS_W Ptr<TonemapMantiuk> 
+createTonemapMantiuk(float gamma = 1.0f, float scale = 0.7f, float saturation = 1.0f);
+
 class CV_EXPORTS_W ExposureAlign : public Algorithm
 {
 public:

modules/photo/src/hdr_common.cpp

@@ -71,4 +71,16 @@ Mat tringleWeights()
     return w;
 }
 
+void mapLuminance(Mat src, Mat dst, Mat lum, Mat new_lum, float saturation)
+{
+    std::vector<Mat> channels(3);
+    split(src, channels);
+    for(int i = 0; i < 3; i++) {
+        channels[i] = channels[i].mul(1.0f / lum);
+        pow(channels[i], saturation, channels[i]);
+        channels[i] = channels[i].mul(new_lum);
+    }
+    merge(channels, dst);
+}
+
 };
\ No newline at end of file

modules/photo/src/hdr_common.hpp

@@ -53,6 +53,8 @@ void checkImageDimensions(const std::vector<Mat>& images);
 
 Mat tringleWeights();
 
+void mapLuminance(Mat src, Mat dst, Mat lum, Mat new_lum, float saturation);
+
 };
 
 #endif
\ No newline at end of file

modules/photo/src/merge.cpp

@@ -109,7 +109,7 @@ public:
         Mat response(256, 3, CV_32F);
         for(int i = 0; i < 256; i++) {
             for(int j = 0; j < 3; j++) {
-				response.at<float>(i, j) = max(i, 1);
+                response.at<float>(i, j) = static_cast<float>(max(i, 1));
             }
         }
         process(src, dst, times, response);

modules/photo/src/tonemap.cpp

@@ -43,6 +43,7 @@
 #include "precomp.hpp"
 #include "opencv2/photo.hpp"
 #include "opencv2/imgproc.hpp"
+#include "hdr_common.hpp"
 
 namespace cv
 {
@@ -101,10 +102,11 @@ Ptr<TonemapLinear> createTonemapLinear(float gamma)
 class TonemapDragoImpl : public TonemapDrago
 {
 public:
-	TonemapDragoImpl(float gamma, float bias) : 
+    TonemapDragoImpl(float gamma, float saturation, float bias) : 
         gamma(gamma), 
+        saturation(saturation),
         bias(bias),
-		name("TonemapLinear")
+        name("TonemapDrago")
     {
     }
 
@@ -135,17 +137,9 @@ public:
         pow(gray_img / static_cast<float>(max), logf(bias) / logf(0.5f), div);
         log(2.0f + 8.0f * div, div);
         map = map.mul(1.0f / div);
-		map = map.mul(1.0f / gray_img);
         div.release();
-		gray_img.release();
 
-		std::vector<Mat> channels(3);
-		split(img, channels);
-		for(int i = 0; i < 3; i++) {
-			channels[i] = channels[i].mul(map);
-		}
-		map.release();
-		merge(channels, img);
+        mapLuminance(img, img, gray_img, map, saturation);
         
         linear->setGamma(gamma);
         linear->process(img, img);
@@ -154,6 +148,9 @@ public:
     float getGamma() const { return gamma; }
     void setGamma(float val) { gamma = val; }
 
+    float getSaturation() const { return saturation; }
+    void setSaturation(float val) { saturation = val; }
+
     float getBias() const { return bias; }
     void setBias(float val) { bias = val; }
 
@@ -161,7 +158,8 @@ public:
     {
         fs << "name" << name
            << "gamma" << gamma
-		   << "bias" << bias;
+           << "bias" << bias
+           << "saturation" << saturation;
     }
 
     void read(const FileNode& fn)
@@ -170,23 +168,25 @@ public:
         CV_Assert(n.isString() && String(n) == name);
         gamma = fn["gamma"];
         bias = fn["bias"];
+        saturation = fn["saturation"];
     }
 
 protected:
     String name;
-	float gamma, bias;
+    float gamma, saturation, bias;
 };
 
-Ptr<TonemapDrago> createTonemapDrago(float gamma, float bias)
+Ptr<TonemapDrago> createTonemapDrago(float gamma, float saturation, float bias)
 {
-	return new TonemapDragoImpl(gamma, bias);
+    return new TonemapDragoImpl(gamma, saturation, bias);
 }
  
 class TonemapDurandImpl : public TonemapDurand
 {
 public:
-	TonemapDurandImpl(float gamma, float contrast, float sigma_color, float sigma_space) : 
+    TonemapDurandImpl(float gamma, float saturation, float contrast, float sigma_color, float sigma_space) : 
         gamma(gamma), 
+        saturation(saturation),
         contrast(contrast),
         sigma_color(sigma_color),
         sigma_space(sigma_space),
@@ -199,10 +199,12 @@ public:
         Mat src = _src.getMat();
         CV_Assert(!src.empty());
         _dst.create(src.size(), CV_32FC3);
-		Mat dst = _dst.getMat();
+        Mat img = _dst.getMat();
+        Ptr<TonemapLinear> linear = createTonemapLinear(1.0f);
+        linear->process(src, img);
 
         Mat gray_img;
-		cvtColor(src, gray_img, COLOR_RGB2GRAY);
+        cvtColor(img, gray_img, COLOR_RGB2GRAY);
         Mat log_img;
         log(gray_img, log_img);
         Mat map_img;
@@ -211,24 +213,19 @@ public:
         double min, max;
         minMaxLoc(map_img, &min, &max);
         float scale = contrast / static_cast<float>(max - min);
-
         exp(map_img * (scale - 1.0f) + log_img, map_img);
         log_img.release();
-		map_img = map_img.mul(1.0f / gray_img);
-		gray_img.release();
 
-		std::vector<Mat> channels(3);
-		split(src, channels);
-		for(int i = 0; i < 3; i++) {
-			channels[i] = channels[i].mul(map_img);
-		}
-		merge(channels, dst);
-		pow(dst, 1.0f / gamma, dst);
+        mapLuminance(img, img, gray_img, map_img, saturation);
+        pow(img, 1.0f / gamma, img);
     }
 
     float getGamma() const { return gamma; }
     void setGamma(float val) { gamma = val; }
 
+    float getSaturation() const { return saturation; }
+    void setSaturation(float val) { saturation = val; }
+
     float getContrast() const { return contrast; }
     void setContrast(float val) { contrast = val; }
 
@@ -244,7 +241,8 @@ public:
            << "gamma" << gamma
            << "contrast" << contrast 
            << "sigma_color" << sigma_color 
-		   << "sigma_space" << sigma_space;
+           << "sigma_space" << sigma_space
+           << "saturation" << saturation;
     }
 
     void read(const FileNode& fn)
@@ -255,16 +253,17 @@ public:
         contrast = fn["contrast"];
         sigma_color = fn["sigma_color"];
         sigma_space = fn["sigma_space"];
+        saturation = fn["saturation"];
     }
 
 protected:
     String name;
-	float gamma, contrast, sigma_color, sigma_space;
+    float gamma, saturation, contrast, sigma_color, sigma_space;
 };
 
-Ptr<TonemapDurand> createTonemapDurand(float gamma, float contrast, float sigma_color, float sigma_space)
+Ptr<TonemapDurand> createTonemapDurand(float gamma, float saturation, float contrast, float sigma_color, float sigma_space)
 {
-	return new TonemapDurandImpl(gamma, contrast, sigma_color, sigma_space);
+    return new TonemapDurandImpl(gamma, saturation, contrast, sigma_color, sigma_space);
 }
 
 class TonemapReinhardDevlinImpl : public TonemapReinhardDevlin
@@ -285,7 +284,6 @@ public:
         CV_Assert(!src.empty());
         _dst.create(src.size(), CV_32FC3);
         Mat img = _dst.getMat();
-		
         Ptr<TonemapLinear> linear = createTonemapLinear(1.0f);
         linear->process(src, img);
         
@@ -363,4 +361,172 @@ Ptr<TonemapReinhardDevlin> createTonemapReinhardDevlin(float gamma, float contra
     return new TonemapReinhardDevlinImpl(gamma, contrast, sigma_color, sigma_space);
 }
 
+class TonemapMantiukImpl : public TonemapMantiuk
+{
+public:
+    TonemapMantiukImpl(float gamma, float scale, float saturation) : 
+        gamma(gamma), 
+        scale(scale),
+        saturation(saturation),
+        name("TonemapMantiuk")
+    {
+    }
+
+    void process(InputArray _src, OutputArray _dst) 
+    {
+        Mat src = _src.getMat();
+        CV_Assert(!src.empty());
+        _dst.create(src.size(), CV_32FC3);
+        Mat img = _dst.getMat();
+        Ptr<TonemapLinear> linear = createTonemapLinear(1.0f);
+        linear->process(src, img);
+
+        Mat gray_img;
+        cvtColor(img, gray_img, COLOR_RGB2GRAY);
+        Mat log_img;
+        log(gray_img, log_img);
+
+        std::vector<Mat> x_contrast, y_contrast;
+        getContrast(log_img, x_contrast, y_contrast);
+
+        for(size_t i = 0; i < x_contrast.size(); i++) {
+            mapContrast(x_contrast[i], scale);
+            mapContrast(y_contrast[i], scale);
+        }
+
+        Mat right(src.size(), CV_32F);
+        calculateSum(x_contrast, y_contrast, right);
+
+        Mat p, r, product, x = log_img;
+        calculateProduct(x, r);
+        r = right - r;
+        r.copyTo(p);
+
+        const float target_error = 1e-3f;
+        float target_norm = static_cast<float>(right.dot(right)) * powf(target_error, 2.0f);
+        int max_iterations = 100;
+        float rr = static_cast<float>(r.dot(r));
+
+        for(int i = 0; i < max_iterations; i++)
+        {
+            calculateProduct(p, product);
+            float alpha = rr / static_cast<float>(p.dot(product));
+
+            r -= alpha * product;
+            x += alpha * p;
+
+            float new_rr = static_cast<float>(r.dot(r));
+            p = r + (new_rr / rr) * p;
+            rr = new_rr;
+
+            if(rr < target_norm) {
+                break;
+            }
+        }
+        exp(x, x);
+        mapLuminance(img, img, gray_img, x, saturation);
+
+        linear = createTonemapLinear(gamma);
+        linear->process(img, img);
+    }
+
+    float getGamma() const { return gamma; }
+    void setGamma(float val) { gamma = val; }
+
+    float getScale() const { return scale; }
+    void setScale(float val) { scale = val; }
+
+    float getSaturation() const { return saturation; }
+    void setSaturation(float val) { saturation = val; }
+
+    void write(FileStorage& fs) const
+    {
+        fs << "name" << name
+           << "gamma" << gamma
+           << "scale" << scale 
+           << "saturation" << saturation;
+    }
+
+    void read(const FileNode& fn)
+    {
+        FileNode n = fn["name"];
+        CV_Assert(n.isString() && String(n) == name);
+        gamma = fn["gamma"];
+        scale = fn["scale"];
+        saturation = fn["saturation"];
+    }
+
+protected:
+    String name;
+    float gamma, scale, saturation;
+
+    void signedPow(Mat src, float power, Mat& dst)
+    {
+        Mat sign = (src > 0);
+        sign.convertTo(sign, CV_32F, 1/255.0f);
+        sign = sign * 2 - 1;
+        pow(abs(src), power, dst);
+        dst = dst.mul(sign);
+    }
+
+    void mapContrast(Mat& contrast, float scale)
+    {
+        const float response_power = 0.4185f;
+        signedPow(contrast, response_power, contrast);
+        contrast *= scale;
+        signedPow(contrast, 1.0f / response_power, contrast);
+    }
+
+    void getGradient(Mat src, Mat& dst, int pos)
+    {
+        dst = Mat::zeros(src.size(), CV_32F);
+        Mat a, b;
+        Mat grad = src.colRange(1, src.cols) - src.colRange(0, src.cols - 1);
+        grad.copyTo(dst.colRange(pos, src.cols + pos - 1));
+        if(pos == 1) {
+            src.col(0).copyTo(dst.col(0));
+        }
+    }
+
+    void getContrast(Mat src, std::vector<Mat>& x_contrast, std::vector<Mat>& y_contrast)
+    {
+        int levels = static_cast<int>(logf(static_cast<float>(min(src.rows, src.cols))) / logf(2.0f));
+        x_contrast.resize(levels);
+        y_contrast.resize(levels);
+
+        Mat layer;
+        src.copyTo(layer);
+        for(int i = 0; i < levels; i++) {
+            getGradient(layer, x_contrast[i], 0);
+            getGradient(layer.t(), y_contrast[i], 0);
+            resize(layer, layer, Size(layer.cols / 2, layer.rows / 2));
+        }
+    }
+
+    void calculateSum(std::vector<Mat>& x_contrast, std::vector<Mat>& y_contrast, Mat& sum)
+    {
+        sum = Mat::zeros(x_contrast[x_contrast.size() - 1].size(), CV_32F);
+        for(int i = x_contrast.size() - 1; i >= 0; i--)
+        {
+            Mat grad_x, grad_y;
+            getGradient(x_contrast[i], grad_x, 1);
+            getGradient(y_contrast[i], grad_y, 1);
+            resize(sum, sum, x_contrast[i].size());
+            sum += grad_x + grad_y.t();
+        }
+    }
+
+    void calculateProduct(Mat src, Mat& dst)
+    {
+        std::vector<Mat> x_contrast, y_contrast;
+        getContrast(src, x_contrast, y_contrast);
+        calculateSum(x_contrast, y_contrast, dst);
+    }
+};
+
+Ptr<TonemapMantiuk> createTonemapMantiuk(float gamma, float scale, float saturation)
+{
+    return new TonemapMantiukImpl(gamma, scale, saturation);
+}
+
 }
\ No newline at end of file

modules/photo/test/test_hdr.cpp

@@ -91,12 +91,11 @@ void loadResponseCSV(String path, Mat& response)
 
 TEST(Photo_Tonemap, regression)
 {
-	string test_path = string(cvtest::TS::ptr()->get_data_path()) + "hdr/";
+	string test_path = string(cvtest::TS::ptr()->get_data_path()) + "hdr/tonemap/";
 	
 	Mat img, expected, result;
-	loadImage(test_path + "rle.hdr", img);
+	loadImage(test_path + "image.hdr", img);
 	float gamma = 2.2f;
-	test_path += "tonemap/";
 	
 	Ptr<TonemapLinear> linear = createTonemapLinear(gamma);
 	linear->process(img, result);
@@ -121,6 +120,12 @@ TEST(Photo_Tonemap, regression)
 	loadImage(test_path + "reinharddevlin.png", expected);
 	result.convertTo(result, CV_8UC3, 255);
 	checkEqual(result, expected, 0);
+
+	Ptr<TonemapMantiuk> mantiuk = createTonemapMantiuk(gamma);
+	mantiuk->process(img, result);
+	loadImage(test_path + "mantiuk.png", expected);
+	result.convertTo(result, CV_8UC3, 255);
+	checkEqual(result, expected, 0);
 }
 
 TEST(Photo_AlignMTB, regression)