diff --git a/modules/intensity_transform/CMakeLists.txt b/modules/intensity_transform/CMakeLists.txt index d68ca5e3b..9bd05532d 100644 --- a/modules/intensity_transform/CMakeLists.txt +++ b/modules/intensity_transform/CMakeLists.txt @@ -1,2 +1,2 @@ set(the_description "Intensity transformations") -ocv_define_module(intensity_transform opencv_core WRAP python) +ocv_define_module(intensity_transform opencv_core opencv_imgproc WRAP python) diff --git a/modules/intensity_transform/doc/intensity_transform.bib b/modules/intensity_transform/doc/intensity_transform.bib index a2a0f72d0..96c96dc83 100644 --- a/modules/intensity_transform/doc/intensity_transform.bib +++ b/modules/intensity_transform/doc/intensity_transform.bib @@ -1,16 +1,28 @@ @book{Gonzalez2018, - title={Digital Image Processing 4th Edition}, - author={Rafael C. Gonzalez, Richard E. Woods}, - year={2018}, - publisher={Pearson} + title = {Digital Image Processing 4th Edition}, + author = {Rafael C. Gonzalez, Richard E. Woods}, + year = {2018}, + publisher = {Pearson} } - @misc{lcs435lab, - title={CS425 Lab: Intensity Transformations and Spatial Filtering}, - url={http://www.cs.uregina.ca/Links/class-info/425/Lab3/} + title = {CS425 Lab: Intensity Transformations and Spatial Filtering}, + url = {http://www.cs.uregina.ca/Links/class-info/425/Lab3/} } - @misc{theailearner, - title={Contrast Stretching}, - url={https://theailearner.com/2019/01/30/contrast-stretching/} -} \ No newline at end of file + title = {Contrast Stretching}, + url = {https://theailearner.com/2019/01/30/contrast-stretching/} +} +@article{ying2017bio, + title = {A Bio-Inspired Multi-Exposure Fusion Framework for Low-light Image Enhancement}, + author = {Ying, Zhenqiang and Li, Ge and Gao, Wen}, + journal = {arXiv preprint arXiv:1711.00591}, + year = {2017} +} +@inproceedings{ying2017new, + title = {A New Image Contrast Enhancement Algorithm Using Exposure Fusion Framework}, + author = {Ying, Zhenqiang and Li, Ge and Ren, Yurui and Wang, Ronggang and Wang, Wenmin}, + booktitle = {International Conference on Computer Analysis of Images and Patterns}, + pages = {36--46}, + year = {2017}, + organization = {Springer} +} diff --git a/modules/intensity_transform/include/opencv2/intensity_transform.hpp b/modules/intensity_transform/include/opencv2/intensity_transform.hpp index 12447f5c9..442f4b95a 100644 --- a/modules/intensity_transform/include/opencv2/intensity_transform.hpp +++ b/modules/intensity_transform/include/opencv2/intensity_transform.hpp @@ -17,6 +17,7 @@ * - Log Transformations * - Power-Law (Gamma) Transformations * - Contrast Stretching + * - BIMEF, A Bio-Inspired Multi-Exposure Fusion Framework for Low-light Image Enhancement @cite ying2017bio @cite ying2017new * * Reference from following book and websites: * - Digital Image Processing 4th Edition Chapter 3 [Rafael C. Gonzalez, Richard E. Woods] @cite Gonzalez2018 @@ -71,8 +72,32 @@ CV_EXPORTS_W void autoscaling(const Mat input, Mat& output); */ CV_EXPORTS_W void contrastStretching(const Mat input, Mat& output, const int r1, const int s1, const int r2, const int s2); +/** + * @brief Given an input color image, enhance low-light images using the BIMEF method (@cite ying2017bio @cite ying2017new). + * + * @param input input color image. + * @param output resulting image. + * @param mu enhancement ratio. + * @param a a-parameter in the Camera Response Function (CRF). + * @param b b-parameter in the Camera Response Function (CRF). +*/ +CV_EXPORTS_W void BIMEF(InputArray input, OutputArray output, float mu=0.5f, float a=-0.3293f, float b=1.1258f); + +/** + * @brief Given an input color image, enhance low-light images using the BIMEF method (@cite ying2017bio @cite ying2017new). + * This is an overload function with the exposure ratio given as parameter. + * + * @param input input color image. + * @param output resulting image. + * @param k exposure ratio. + * @param mu enhancement ratio. + * @param a a-parameter in the Camera Response Function (CRF). + * @param b b-parameter in the Camera Response Function (CRF). +*/ +CV_EXPORTS_AS(BIMEF2) void BIMEF(InputArray input, OutputArray output, float k, float mu, float a, float b); + //! @} }} // cv::intensity_transform:: -#endif \ No newline at end of file +#endif diff --git a/modules/intensity_transform/misc/python/test/test_intensity_transform_bimef.py b/modules/intensity_transform/misc/python/test/test_intensity_transform_bimef.py new file mode 100644 index 000000000..0e711eb26 --- /dev/null +++ b/modules/intensity_transform/misc/python/test/test_intensity_transform_bimef.py @@ -0,0 +1,44 @@ +#!/usr/bin/env python + +# Python 2/3 compatibility +from __future__ import print_function + +import os +import cv2 as cv +import numpy as np + +from tests_common import NewOpenCVTests, unittest + +class intensity_transform_test(NewOpenCVTests): + def setUp(self): + super(intensity_transform_test, self).setUp() + try: + result_ = cv.intensity_transform.BIMEF(None) + except cv.error as e: + if e.code == cv.Error.StsNotImplemented: + self.skipTest('BIMEF is not implemented (missing Eigen dependency)') + + @unittest.skipIf('OPENCV_TEST_DATA_PATH' not in os.environ, + "OPENCV_TEST_DATA_PATH is not defined") + def test_BIMEF(self): + filenames = ['P1000205_resize', 'P1010676_resize', 'P1010815_resize'] + + for f in filenames: + img = self.get_sample('cv/intensity_transform/BIMEF/{}.png'.format(f)) + self.assertTrue(img.size > 0) + + img_ref = self.get_sample('cv/intensity_transform/BIMEF/{}_ref.png'.format(f)) + self.assertTrue(img_ref.size > 0) + + img_BIMEF = cv.intensity_transform.BIMEF(img) + self.assertTrue(img_BIMEF.size > 0) + self.assertTrue(img_BIMEF.shape == img_ref.shape) + self.assertTrue(img_BIMEF.dtype == img_ref.dtype) + + RMSE = np.sqrt(cv.norm(img_BIMEF, img_ref, cv.NORM_L2SQR) / (img_ref.shape[0]*img_ref.shape[1]*img_ref.shape[2])) + max_RMSE_threshold = 9.0 + self.assertLessEqual(RMSE, max_RMSE_threshold) + print('BIMEF RMSE:', RMSE) + +if __name__ == '__main__': + NewOpenCVTests.bootstrap() diff --git a/modules/intensity_transform/samples/intensity_transform.cpp b/modules/intensity_transform/samples/intensity_transform.cpp index 6274aaf2a..2b547d402 100644 --- a/modules/intensity_transform/samples/intensity_transform.cpp +++ b/modules/intensity_transform/samples/intensity_transform.cpp @@ -9,31 +9,115 @@ using namespace std; using namespace cv; using namespace cv::intensity_transform; +namespace +{ +static std::string keys = + "{ help h | | Print help message. }" + "{ input i | | Path to the input image. }"; + +// global variables +Mat g_image; + +int g_gamma = 40; +const int g_gammaMax = 500; +Mat g_imgGamma; +const std::string g_gammaWinName = "Gamma Correction"; + +Mat g_contrastStretch; +int g_r1 = 70; +int g_s1 = 15; +int g_r2 = 120; +int g_s2 = 240; +const std::string g_contrastWinName = "Contrast Stretching"; + +Mat g_imgBIMEF; +int g_mu = 50; +const int g_muMax = 100; +const std::string g_BIMEFWinName = "BIMEF"; + +static void onTrackbarGamma(int, void*) +{ + float gamma = g_gamma / 100.0f; + gammaCorrection(g_image, g_imgGamma, gamma); + imshow(g_gammaWinName, g_imgGamma); +} + +static void onTrackbarContrastR1(int, void*) +{ + contrastStretching(g_image, g_contrastStretch, g_r1, g_s1, g_r2, g_s2); + imshow("Contrast Stretching", g_contrastStretch); +} + +static void onTrackbarContrastS1(int, void*) +{ + contrastStretching(g_image, g_contrastStretch, g_r1, g_s1, g_r2, g_s2); + imshow("Contrast Stretching", g_contrastStretch); +} + +static void onTrackbarContrastR2(int, void*) +{ + contrastStretching(g_image, g_contrastStretch, g_r1, g_s1, g_r2, g_s2); + imshow("Contrast Stretching", g_contrastStretch); +} + +static void onTrackbarContrastS2(int, void*) +{ + contrastStretching(g_image, g_contrastStretch, g_r1, g_s1, g_r2, g_s2); + imshow("Contrast Stretching", g_contrastStretch); +} + +static void onTrackbarBIMEF(int, void*) +{ + float mu = g_mu / 100.0f; + BIMEF(g_image, g_imgBIMEF, mu); + imshow(g_BIMEFWinName, g_imgBIMEF); +} +} + int main(int argc, char **argv) { - if (argc != 2) + CommandLineParser parser(argc, argv, keys); + + const std::string inputFilename = parser.get("input"); + parser.about("Use this script to apply intensity transformation on an input image."); + if (parser.has("help") || inputFilename.empty()) { - cerr << "Must input the path of the input image. Ex: intensity_transform image.jpg" << endl; - return -1; + parser.printMessage(); + return 0; } // Read input image - Mat image = imread(argv[1]); + g_image = imread(inputFilename); + + // Create trackbars + namedWindow(g_gammaWinName); + createTrackbar("Gamma value", g_gammaWinName, &g_gamma, g_gammaMax, onTrackbarGamma); + + namedWindow(g_contrastWinName); + createTrackbar("Contrast R1", g_contrastWinName, &g_r1, 256, onTrackbarContrastR1); + createTrackbar("Contrast S1", g_contrastWinName, &g_s1, 256, onTrackbarContrastS1); + createTrackbar("Contrast R2", g_contrastWinName, &g_r2, 256, onTrackbarContrastR2); + createTrackbar("Contrast S2", g_contrastWinName, &g_s2, 256, onTrackbarContrastS2); + + namedWindow(g_BIMEFWinName); + createTrackbar("Enhancement ratio mu", g_BIMEFWinName, &g_mu, g_muMax, onTrackbarBIMEF); // Apply intensity transformations - Mat imgGamma, imgAutoscaled, imgLog, contrastStretch; - gammaCorrection(image, imgGamma, (float)(0.4)); - autoscaling(image, imgAutoscaled); - logTransform(image, imgLog); - contrastStretching(image, contrastStretch, 70, 15, 120, 240); + Mat imgAutoscaled, imgLog; + autoscaling(g_image, imgAutoscaled); + gammaCorrection(g_image, g_imgGamma, g_gamma/100.0f); + logTransform(g_image, imgLog); + contrastStretching(g_image, g_contrastStretch, g_r1, g_s1, g_r2, g_s2); + BIMEF(g_image, g_imgBIMEF, g_mu / 100.0f); // Display intensity transformation results - imshow("Original Image", image); + imshow("Original Image", g_image); imshow("Autoscale", imgAutoscaled); - imshow("Gamma Correction", imgGamma); + imshow(g_gammaWinName, g_imgGamma); imshow("Log Transformation", imgLog); - imshow("Contrast Stretching", contrastStretch); - waitKey(0); + imshow(g_contrastWinName, g_contrastStretch); + imshow(g_BIMEFWinName, g_imgBIMEF); + waitKey(0); return 0; -} \ No newline at end of file +} diff --git a/modules/intensity_transform/src/bimef.cpp b/modules/intensity_transform/src/bimef.cpp new file mode 100644 index 000000000..64a5f55ae --- /dev/null +++ b/modules/intensity_transform/src/bimef.cpp @@ -0,0 +1,572 @@ +// This file is part of OpenCV project. +// It is subject to the license terms in the LICENSE file found in the top-level directory +// of this distribution and at http://opencv.org/license.html. + +/* + * MIT License + * + * Copyright (c) 2017 Zhenqiang.Ying + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ + +#include "precomp.hpp" + +#ifdef HAVE_EIGEN +#include +#include +#include +#endif + +namespace cv { +namespace intensity_transform { + +#ifdef HAVE_EIGEN +static void diff(const Mat_& src, Mat_& srcVDiff, Mat_& srcHDiff) +{ + srcVDiff = Mat_(src.size()); + for (int i = 0; i < src.rows; i++) + { + if (i < src.rows-1) + { + for (int j = 0; j < src.cols; j++) + { + srcVDiff(i,j) = src(i+1,j) - src(i,j); + } + } + else + { + for (int j = 0; j < src.cols; j++) + { + srcVDiff(i,j) = src(0,j) - src(i,j); + } + } + } + + srcHDiff = Mat_(src.size()); + for (int j = 0; j < src.cols-1; j++) + { + for (int i = 0; i < src.rows; i++) + { + srcHDiff(i,j) = src(i,j+1) - src(i,j); + } + } + for (int i = 0; i < src.rows; i++) + { + srcHDiff(i,src.cols-1) = src(i,0) - src(i,src.cols-1); + } +} + +static void computeTextureWeights(const Mat_& x, float sigma, float sharpness, Mat_& W_h, Mat_& W_v) +{ + Mat_ dt0_v, dt0_h; + diff(x, dt0_v, dt0_h); + + Mat_ gauker_h; + Mat_ kernel_h = Mat_::ones(1, static_cast(sigma)); + filter2D(dt0_h, gauker_h, -1, kernel_h, Point(-1,-1), 0, BORDER_CONSTANT); + + Mat_ gauker_v; + Mat_ kernel_v = Mat_::ones(static_cast(sigma), 1); + filter2D(dt0_v, gauker_v, -1, kernel_v, Point(-1,-1), 0, BORDER_CONSTANT); + + W_h = Mat_(gauker_h.size()); + W_v = Mat_(gauker_v.size()); + + for (int i = 0; i < gauker_h.rows; i++) + { + for (int j = 0; j < gauker_h.cols; j++) + { + W_h(i,j) = 1 / (std::abs(gauker_h(i,j)) * std::abs(dt0_h(i,j)) + sharpness); + W_v(i,j) = 1 / (std::abs(gauker_v(i,j)) * std::abs(dt0_v(i,j)) + sharpness); + } + } +} + +template +static Eigen::SparseMatrix spdiags(const Eigen::Matrix &B, + const Eigen::VectorXi &d, int m, int n) { + typedef Eigen::Triplet triplet_t; + std::vector triplets; + triplets.reserve(static_cast(std::min(m,n)*d.size())); + + for (int k = 0; k < d.size(); ++k) { + int diag = d(k); // get diagonal + int i_start = std::max(-diag, 0); // get row of 1st element + int i_end = std::min(m, m-diag-(m-n)); // get row of last element + int j = -std::min(0, -diag); // get col of 1st element + int B_i; // start index i in matrix B + if (m < n) { + B_i = std::max(-diag,0); // m < n + } else { + B_i = std::max(0,diag); // m >= n + } + for (int i = i_start; i < i_end; ++i, ++j, ++B_i) { + triplets.push_back( {i, j, B(B_i,k)} ); + } + } + Eigen::SparseMatrix A(m,n); + A.setFromTriplets(triplets.begin(), triplets.end()); + return A; +} + + +static Mat solveLinearEquation(const Mat_& img, Mat_& W_h_, Mat_& W_v_, float lambda) +{ + Eigen::MatrixXf W_h; + cv2eigen(W_h_, W_h); + Eigen::MatrixXf tempx(W_h.rows(), W_h.cols()); + tempx.block(0, 1, tempx.rows(), tempx.cols()-1) = W_h.block(0, 0, W_h.rows(), W_h.cols()-1); + for (Eigen::Index i = 0; i < tempx.rows(); i++) + { + tempx(i,0) = W_h(i, W_h.cols()-1); + } + + Eigen::MatrixXf W_v; + cv2eigen(W_v_, W_v); + Eigen::MatrixXf tempy(W_v.rows(), W_v.cols()); + tempy.block(1, 0, tempx.rows()-1, tempx.cols()) = W_v.block(0, 0, W_v.rows()-1, W_v.cols()); + for (Eigen::Index j = 0; j < tempy.cols(); j++) + { + tempy(0,j) = W_v(W_v.rows()-1, j); + } + + + Eigen::VectorXf dx(W_h.rows()*W_h.cols()); + Eigen::VectorXf dy(W_v.rows()*W_v.cols()); + + Eigen::VectorXf dxa(tempx.rows()*tempx.cols()); + Eigen::VectorXf dya(tempy.rows()*tempy.cols()); + + //Flatten in a col-major order + for (Eigen::Index j = 0; j < W_h.cols(); j++) + { + for (Eigen::Index i = 0; i < W_h.rows(); i++) + { + dx(j*W_h.rows() + i) = -lambda*W_h(i,j); + dy(j*W_h.rows() + i) = -lambda*W_v(i,j); + + dxa(j*W_h.rows() + i) = -lambda*tempx(i,j); + dya(j*W_h.rows() + i) = -lambda*tempy(i,j); + } + } + + tempx.setZero(); + tempx.col(0) = W_h.col(W_h.cols()-1); + + tempy.setZero(); + tempy.row(0) = W_v.row(W_v.rows()-1); + + W_h.col(W_h.cols()-1).setZero(); + W_v.row(W_v.rows()-1).setZero(); + + Eigen::VectorXf dxd1(tempx.rows()*tempx.cols()); + Eigen::VectorXf dyd1(tempy.rows()*tempy.cols()); + Eigen::VectorXf dxd2(W_h.rows()*W_h.cols()); + Eigen::VectorXf dyd2(W_v.rows()*W_v.cols()); + + //Flatten in a col-major order + for (Eigen::Index j = 0; j < tempx.cols(); j++) + { + for (Eigen::Index i = 0; i < tempx.rows(); i++) + { + dxd1(j*tempx.rows() + i) = -lambda*tempx(i,j); + dyd1(j*tempx.rows() + i) = -lambda*tempy(i,j); + + dxd2(j*tempx.rows() + i) = -lambda*W_h(i,j); + dyd2(j*tempx.rows() + i) = -lambda*W_v(i,j); + } + } + + Eigen::MatrixXf dxd(dxd1.rows(), dxd1.cols()+dxd2.cols()); + dxd << dxd1, dxd2; + + Eigen::MatrixXf dyd(dyd1.rows(), dyd1.cols()+dyd2.cols()); + dyd << dyd1, dyd2; + + const int k = img.rows*img.cols; + const int r = img.rows; + Eigen::Matrix diagx_idx; + diagx_idx << -k+r, -r; + Eigen::SparseMatrix Ax = spdiags(dxd, diagx_idx, k, k); + + Eigen::Matrix diagy_idx; + diagy_idx << -r+1, -1; + Eigen::SparseMatrix Ay = spdiags(dyd, diagy_idx, k, k); + + Eigen::MatrixXf D = (dx + dy + dxa + dya); + D = Eigen::MatrixXf::Ones(D.rows(), D.cols()) - D; + + Eigen::Matrix diag_idx_zero; + diag_idx_zero << 0; + Eigen::SparseMatrix A = (Ax + Ay) + Eigen::SparseMatrix((Ax + Ay).transpose()) + spdiags(D, diag_idx_zero, k, k); + + //CG solver of Eigen + Eigen::ConjugateGradient, Eigen::Lower|Eigen::Upper, Eigen::IncompleteCholesky > cg; + cg.setTolerance(0.1f); + cg.setMaxIterations(50); + cg.compute(A); + Mat_ img_t = img.t(); + Eigen::Map tin(img_t.ptr(), img_t.rows*img_t.cols); + Eigen::VectorXf x = cg.solve(tin); + + Mat_ tout(img.rows, img.cols); + tout.forEach( + [&](float &pixel, const int * position) -> void + { + pixel = x(position[1]*img.rows + position[0]); + } + ); + + return tout; +} + +static Mat_ tsmooth(const Mat_& src, float lambda=0.01f, float sigma=3.0f, float sharpness=0.001f) +{ + Mat_ W_h, W_v; + computeTextureWeights(src, sigma, sharpness, W_h, W_v); + + Mat_ S = solveLinearEquation(src, W_h, W_v, lambda); + + return S; +} + +static Mat_ rgb2gm(const Mat_& I) +{ + Mat_ gm(I.rows, I.cols); + gm.forEach( + [&](float &pixel, const int * position) -> void + { + pixel = std::pow(I(position[0], position[1])[0]*I(position[0], position[1])[1]*I(position[0], position[1])[2], 1/3.0f); + } + ); + + return gm; +} + +static Mat_ applyK(const Mat_& I, float k, float a=-0.3293f, float b=1.1258f) { + float beta = std::exp((1 - std::pow(k, a)) * b); + float gamma = std::pow(k, a); + + Mat_ J(I.size()); + pow(I, gamma, J); + J = J*beta; + + return J; +} + +static Mat_ applyK(const Mat_& I, float k, float a=-0.3293f, float b=1.1258f, float offset=0) { + float beta = std::exp((1 - std::pow(k, a)) * b); + float gamma = std::pow(k, a); + + Mat_ J(I.size()); + pow(I, gamma, J); + + return J * beta + Scalar::all(offset); +} + +static float entropy(const Mat_& I) +{ + Mat_ I_uchar; + I.convertTo(I_uchar, CV_8U, 255); + + std::vector planes; + planes.push_back(I_uchar); + Mat_ hist; + const int histSize = 256; + float range[] = { 0, 256 }; + const float* histRange = { range }; + calcHist(&I_uchar, 1, NULL, Mat(), hist, 1, &histSize, &histRange); + + Mat_ hist_norm = hist / cv::sum(hist)[0]; + + float E = 0; + for (int i = 0; i < hist_norm.rows; i++) + { + if (hist_norm(i,0) > 0) + { + E += hist_norm(i,0) * std::log2(hist_norm(i,0)); + } + } + + return -E; +} + +template static int sgn(T val) +{ + return (T(0) < val) - (val < T(0)); +} + +static double minimize_scalar_bounded(const Mat_& I, double begin, double end, + double xatol=1e-4, int maxiter=500) +{ +// From scipy: https://github.com/scipy/scipy/blob/v1.4.1/scipy/optimize/optimize.py#L1753-L1894 +// """ +// Options +// ------- +// maxiter : int +// Maximum number of iterations to perform. +// disp: int, optional +// If non-zero, print messages. +// 0 : no message printing. +// 1 : non-convergence notification messages only. +// 2 : print a message on convergence too. +// 3 : print iteration results. +// xatol : float +// Absolute error in solution `xopt` acceptable for convergence. +// """ + double x1 = begin, x2 = end; + + if (x1 > x2) { + std::runtime_error("The lower bound exceeds the upper bound."); + } + + double sqrt_eps = std::sqrt(2.2e-16); + double golden_mean = 0.5 * (3.0 - std::sqrt(5.0)); + double a = x1, b = x2; + double fulc = a + golden_mean * (b - a); + double nfc = fulc, xf = fulc; + double rat = 0.0, e = 0.0; + double x = xf; + double fx = -entropy(applyK(I, static_cast(x))); + int num = 1; + double fu = std::numeric_limits::infinity(); + + double ffulc = fx, fnfc = fx; + double xm = 0.5 * (a + b); + double tol1 = sqrt_eps * std::abs(xf) + xatol / 3.0; + double tol2 = 2.0 * tol1; + + for (int iter = 0; iter < maxiter && std::abs(xf - xm) > (tol2 - 0.5 * (b - a)); iter++) + { + int golden = 1; + // Check for parabolic fit + if (std::abs(e) > tol1) { + golden = 0; + double r = (xf - nfc) * (-entropy(applyK(I, static_cast(x))) - ffulc); + double q = (xf - fulc) * (-entropy(applyK(I, static_cast(x))) - fnfc); + double p = (xf - fulc) * q - (xf - nfc) * r; + q = 2.0 * (q - r); + + if (q > 0.0) { + p = -p; + } + q = std::abs(q); + r = e; + e = rat; + + // Check for acceptability of parabola + if (((std::abs(p) < std::abs(0.5*q*r)) && (p > q*(a - xf)) & + (p < q * (b - xf)))) { + rat = (p + 0.0) / q; + x = xf + rat; + + if (((x - a) < tol2) || ((b - x) < tol2)) { + double si = sgn(xm - xf) + ((xm - xf) == 0); + rat = tol1 * si; + } + } else { // do a golden-section step + golden = 1; + } + } + + if (golden) { // do a golden-section step + if (xf >= xm) { + e = a - xf; + } else { + e = b - xf; + } + rat = golden_mean*e; + } + + double si = sgn(rat) + (rat == 0); + x = xf + si * std::max(std::abs(rat), tol1); + fu = -entropy(applyK(I, static_cast(x))); + num += 1; + + if (fu <= fx) { + if (x >= xf) { + a = xf; + } else { + b = xf; + } + + fulc = nfc; + ffulc = fnfc; + nfc = xf; + fnfc = fx; + xf = x; + fx = fu; + } else { + if (x < xf) { + a = x; + } else { + b = x; + } + + if ((fu <= fnfc) || (nfc == xf)) { + fulc = nfc; + ffulc = fnfc; + nfc = x; + fnfc = fu; + } else if ((fu <= ffulc) || (fulc == xf) || (fulc == nfc)) { + fulc = x; + ffulc = fu; + } + } + + xm = 0.5 * (a + b); + tol1 = sqrt_eps * std::abs(xf) + xatol / 3.0; + tol2 = 2.0 * tol1; + } + + return xf; +} + +static Mat_ maxEntropyEnhance(const Mat_& I, const Mat_& isBad, float a, float b) +{ + Mat_ input; + resize(I, input, Size(50,50)); + + Mat_ Y = rgb2gm(input); + + Mat_ isBad_resize; + resize(isBad, isBad_resize, Size(50,50)); + + std::vector Y_vec; + for (int i = 0; i < isBad_resize.rows; i++) + { + for (int j = 0; j < isBad_resize.cols; j++) + { + if (isBad_resize(i,j) >= 0.5) + { + Y_vec.push_back(Y(i,j)); + } + } + } + + if (Y_vec.empty()) + { + return I; + } + + Mat_ Y_mat(static_cast(Y_vec.size()), 1, Y_vec.data()); + float opt_k = static_cast(minimize_scalar_bounded(Y_mat, 1, 7)); + + return applyK(I, opt_k, a, b, -0.01f); +} + +static void BIMEF_impl(InputArray input_, OutputArray output_, float mu, float *k, float a, float b) +{ + CV_INSTRUMENT_REGION() + + Mat input = input_.getMat(); + if (input.empty()) + { + return; + } + CV_CheckTypeEQ(input.type(), CV_8UC3, "Input image must be 8-bits color image (CV_8UC3)."); + + Mat_ imgDouble; + input.convertTo(imgDouble, CV_32F, 1/255.0); + + // t: scene illumination map + Mat_ t_b(imgDouble.size()); + t_b.forEach( + [&](float &pixel, const int * position) -> void + { + pixel = std::max(std::max(imgDouble(position[0], position[1])[0], + imgDouble(position[0], position[1])[1]), + imgDouble(position[0], position[1])[2]); + } + ); + + const float lambda = 0.5; + const float sigma = 5; + + Mat_ t_b_resize; + resize(t_b, t_b_resize, Size(), 0.5, 0.5); + + Mat_ t_our = tsmooth(t_b_resize, lambda, sigma); + resize(t_our, t_our, t_b.size()); + + // k: exposure ratio + Mat_ J; + if (k == NULL) + { + Mat_ isBad(t_our.size()); + isBad.forEach( + [&](uchar &pixel, const int * position) -> void + { + pixel = t_our(position[0], position[1]) < 0.5 ? 1 : 0; + } + ); + + J = maxEntropyEnhance(imgDouble, isBad, a, b); + } + else + { + J = applyK(imgDouble, *k, a, b); + + // fix overflow + J.forEach( + [](Vec3f &pixel, const int * /*position*/) -> void + { + pixel(0) = std::min(1.0f, pixel(0)); + pixel(1) = std::min(1.0f, pixel(1)); + pixel(2) = std::min(1.0f, pixel(2)); + } + ); + } + + // W: Weight Matrix + Mat_ W(t_our.size()); + pow(t_our, mu, W); + + + output_.create(input.size(), CV_8UC3); + Mat output = output_.getMat(); + output.forEach( + [&](Vec3b &pixel, const int * position) -> void + { + float w = W(position[0], position[1]); + pixel(0) = saturate_cast((imgDouble(position[0], position[1])[0] * w + J(position[0], position[1])[0] * (1 - w)) * 255); + pixel(1) = saturate_cast((imgDouble(position[0], position[1])[1] * w + J(position[0], position[1])[1] * (1 - w)) * 255); + pixel(2) = saturate_cast((imgDouble(position[0], position[1])[2] * w + J(position[0], position[1])[2] * (1 - w)) * 255); + } + ); +} +#else +static void BIMEF_impl(InputArray, OutputArray, float, float *, float, float) +{ + CV_Error(Error::StsNotImplemented, "This algorithm requires OpenCV built with the Eigen library."); +} +#endif + +void BIMEF(InputArray input, OutputArray output, float mu, float a, float b) +{ + BIMEF_impl(input, output, mu, NULL, a, b); +} + +void BIMEF(InputArray input, OutputArray output, float k, float mu, float a, float b) +{ + BIMEF_impl(input, output, mu, &k, a, b); +} + +}} // cv::intensity_transform:: diff --git a/modules/intensity_transform/src/precomp.hpp b/modules/intensity_transform/src/precomp.hpp index 6708d5bba..733057145 100644 --- a/modules/intensity_transform/src/precomp.hpp +++ b/modules/intensity_transform/src/precomp.hpp @@ -7,6 +7,7 @@ #include "opencv2/core.hpp" #include "opencv2/intensity_transform.hpp" +#include "opencv2/core/private.hpp" #include #endif diff --git a/modules/intensity_transform/test/test_intensity_transform.cpp b/modules/intensity_transform/test/test_intensity_transform.cpp index 5780904d0..bc0419761 100644 --- a/modules/intensity_transform/test/test_intensity_transform.cpp +++ b/modules/intensity_transform/test/test_intensity_transform.cpp @@ -198,4 +198,42 @@ TEST(intensity_transform_contrastStretching, accuracy) EXPECT_LE(cvtest::norm(res, expectedRes, NORM_INF), 1); } +typedef testing::TestWithParam intensity_transform_BIMEF; + +TEST_P(intensity_transform_BIMEF, accuracy) +{ +#ifdef HAVE_EIGEN + const std::string directory = "intensity_transform/BIMEF/"; + std::string filename = GetParam(); + + const std::string inputFilename = cvtest::findDataFile(directory + filename + ".png"); + Mat img = imread(inputFilename); + EXPECT_TRUE(!img.empty()); + Mat imgBIMEF; + BIMEF(img, imgBIMEF); + + const std::string referenceFilename = cvtest::findDataFile(directory + filename + "_ref.png"); + Mat imgRef = imread(referenceFilename); + EXPECT_TRUE(!imgRef.empty()); + + EXPECT_EQ(imgBIMEF.rows, imgRef.rows); + EXPECT_EQ(imgBIMEF.cols, imgRef.cols); + EXPECT_EQ(imgBIMEF.type(), imgRef.type()); + double rmse = sqrt(cv::norm(imgBIMEF, imgRef, NORM_L2SQR) / (imgRef.total()*imgRef.channels())); + std::cout << "BIMEF, RMSE for " << filename << ": " << rmse << std::endl; + const float max_rmse = 9; + EXPECT_LE(rmse, max_rmse); +#endif +} + +const string BIMEF_accuracy_cases[] = { + "P1000205_resize", + "P1010676_resize", + "P1010815_resize" +}; + +INSTANTIATE_TEST_CASE_P(/*nothing*/, intensity_transform_BIMEF, + testing::ValuesIn(BIMEF_accuracy_cases) +); + }} // namespace