From ccc71ac190b3a23f0b8837001c1f142beb3975ae Mon Sep 17 00:00:00 2001 From: Alex Leontiev Date: Mon, 19 Aug 2013 17:19:52 +0800 Subject: [PATCH 1/4] Primal-dual algorithm This is an implementation of primal-dual algorithm, based on the C++ source code by Vadim Pisarevsky. It was extended to handle the denoising based on multiple observations. It also contains documentation and tests. --- modules/optim/doc/optim.rst | 1 + modules/optim/doc/primal_dual_algorithm.rst | 51 ++++++ modules/optim/include/opencv2/optim.hpp | 1 + modules/optim/src/denoise_tvl1.cpp | 183 ++++++++++++++++++++ modules/optim/test/test_denoise_tvl1.cpp | 81 +++++++++ 5 files changed, 317 insertions(+) create mode 100644 modules/optim/doc/primal_dual_algorithm.rst create mode 100644 modules/optim/src/denoise_tvl1.cpp create mode 100644 modules/optim/test/test_denoise_tvl1.cpp diff --git a/modules/optim/doc/optim.rst b/modules/optim/doc/optim.rst index cdbaeac25d..b3c7a740b1 100644 --- a/modules/optim/doc/optim.rst +++ b/modules/optim/doc/optim.rst @@ -9,3 +9,4 @@ optim. Generic numerical optimization linear_programming downhill_simplex_method + primal_dual_algorithm diff --git a/modules/optim/doc/primal_dual_algorithm.rst b/modules/optim/doc/primal_dual_algorithm.rst new file mode 100644 index 0000000000..1eda00f596 --- /dev/null +++ b/modules/optim/doc/primal_dual_algorithm.rst @@ -0,0 +1,51 @@ +Primal-Dual Algorithm +======================= + +.. highlight:: cpp + +optim::denoise_TVL1 +--------------------------------- + +Primal-dual algorithm is an algorithm for solving special types of variational +problems (that is, finding a function to minimize some functional) +. As the image denoising, in particular, may be seen as the variational +problem, primal-dual algorithm then can be used to perform denoising and this +is exactly what is implemented. + +It should be noted, that this implementation was taken from the July 2013 blog entry [Mordvintsev]_, which also contained +(slightly more general) ready-to-use +source code on Python. Subsequently, that code was rewritten on C++ with the usage of openCV by Vadim Pisarevsky +at the end of July 2013 and finally it was slightly adapted by later authors. + +Although the thorough discussion and justification +of the algorithm involved may be found in [ChambolleEtAl]_, it might make sense to skim over it here, following [Mordvintsev]_. To +begin with, we consider the 1-byte gray-level images as the functions from the rectangular domain of pixels +(it may be seen as set :math:`\left\{(x,y)\in\mathbb{N}\times\mathbb{N}\mid 1\leq x\leq n,\;1\leq y\leq m\right\}` +for some :math:`m,\;n\in\mathbb{N}`) into :math:`\{0,1,\dots,255\}`. We shall denote the noised images as :math:`f_i` and with this +view, given some image :math:`x` of the same size, we may measure how bad it is by the formula + +.. math:: + \left\|\left\|\nabla x\right\|\right\| + \lambda\sum_i\left\|\left\|x-f_i\right\|\right\| + +:math:`\|\|\cdot\|\|` here denotes :math:`L_2`-norm and as you see, the first addend states that we want our image to be smooth +(ideally, having zero gradient, thus being constant) and the second states that we want our result to be close to the observations we've got. +If we treat :math:`x` as a function, this is exactly the functional what we seek to minimize and here the Primal-Dual algorithm comes +into play. + +.. ocv:function:: void optim::denoise_TVL1(const std::vector& observations,Mat& result, double lambda, int niters) + + :param observations: This array should contain one or more noised versions of the image that is to be restored. + + :param result: Here the denoised image will be stored. There is no need to do pre-allocation of storage space, as it will be automatically allocated, if necessary. + + :param lambda: Corresponds to :math:`\lambda` in the formulas above. As it is enlarged, the smooth (blurred) images are treated + more favorably than detailed (but maybe more noised) ones. Roughly speaking, as it becomes smaller, the result will be + more blur but more sever outliers will be removed. + + + :param niters: Number of iterations that the algorithm will run. Of course, as more iterations as better, but it is hard to quantitatively refine this statement, so just use the default and increase it if the results are poor. + + +.. [ChambolleEtAl] A. Chambolle, V. Caselles, M. Novaga, D. Cremers and T. Pock, An Introduction to Total Variation for Image Analysis, http://hal.archives-ouvertes.fr/docs/00/43/75/81/PDF/preprint.pdf (pdf) + +.. [Mordvintsev] Alexander Mordvintsev, ROF and TV-L1 denoising with Primal-Dual algorithm, http://znah.net/rof-and-tv-l1-denoising-with-primal-dual-algorithm.html (blog entry) diff --git a/modules/optim/include/opencv2/optim.hpp b/modules/optim/include/opencv2/optim.hpp index 9acd95d02e..715372b69a 100644 --- a/modules/optim/include/opencv2/optim.hpp +++ b/modules/optim/include/opencv2/optim.hpp @@ -96,6 +96,7 @@ enum }; CV_EXPORTS_W int solveLP(const Mat& Func, const Mat& Constr, Mat& z); +CV_EXPORTS_W void denoise_TVL1(const std::vector& observations,Mat& result, double lambda=1.0, int niters=30); }}// cv #endif diff --git a/modules/optim/src/denoise_tvl1.cpp b/modules/optim/src/denoise_tvl1.cpp new file mode 100644 index 0000000000..47f759affe --- /dev/null +++ b/modules/optim/src/denoise_tvl1.cpp @@ -0,0 +1,183 @@ +#include "precomp.hpp" +#define ALEX_DEBUG +#include "debug.hpp" +#include +#include + +#define ABSCLIP(val,threshold) MIN(MAX((val),-(threshold)),(threshold)) + +namespace cv{namespace optim{ + + class AddFloatToCharScaled{ + public: + AddFloatToCharScaled(float scale):_scale(scale){} + inline float operator()(float a,uchar b){ + return a+_scale*((float)b); + } + private: + float _scale; + }; + + void solve_TVL1(const Mat& img, Mat& res, double _clambda, int niters) + { + const float L2 = 8.0f, tau = 0.02f, sigma = 1./(L2*tau), theta = 1.f, img_scale = 1.f/255; + float clambda = (float)_clambda, threshold = clambda*tau; + const int workdepth = CV_32F; + + int i, x, y, rows=img.rows, cols=img.cols; + + Mat X, P = Mat::zeros(rows, cols, CV_MAKETYPE(workdepth, 2)); + img.convertTo(X, workdepth, 1./255); + + for( i = 0; i < niters; i++ ) + { + float currsigma = i == 0 ? 1 + sigma : sigma; + + // P_ = P + sigma*nabla(X) + // P(x,y) = P_(x,y)/max(||P(x,y)||,1) + for( y = 0; y < rows; y++ ) + { + const float* x_curr = X.ptr(y); + const float* x_next = X.ptr(std::min(y+1, rows-1)); + Point2f* p_curr = P.ptr(y); + float dx, dy, m; + for( x = 0; x < cols-1; x++ ) + { + dx = (x_curr[x+1] - x_curr[x])*currsigma + p_curr[x].x; + dy = (x_next[x] - x_curr[x])*currsigma + p_curr[x].y; + m = 1.f/std::max(std::sqrt(dx*dx + dy*dy), 1.f); + p_curr[x].x = dx*m; + p_curr[x].y = dy*m; + } + dy = (x_next[x] - x_curr[x])*currsigma + p_curr[x].y; + m = 1.f/std::max(std::abs(dy), 1.f); + p_curr[x].x = 0.f; + p_curr[x].y = dy*m; + } + + // X1 = X + tau*(-nablaT(P)) + // X2 = X1 + clip(img - X1, -clambda*tau, clambda*tau) + // X = X2 + theta*(X2 - X) + for( y = 0; y < rows; y++ ) + { + const uchar* img_curr = img.ptr(y); + float* x_curr = X.ptr(y); + const Point2f* p_curr = P.ptr(y); + const Point2f* p_prev = P.ptr(std::max(y - 1, 0)); + + x = 0; + float x_new = x_curr[x] + tau*(p_curr[x].y - p_prev[x].y); + x_new += std::min(std::max(img_curr[x]*img_scale - x_new, -threshold), threshold); + x_curr[x] = x_new + theta*(x_new - x_curr[x]); + + for( x = 1; x < cols; x++ ) + { + x_new = x_curr[x] + tau*(p_curr[x].x - p_curr[x-1].x + p_curr[x].y - p_prev[x].y); + x_new += std::min(std::max(img_curr[x]*img_scale - x_new, -threshold), threshold); + x_curr[x] = x_new + theta*(x_new - x_curr[x]); + } + } + } + + res.create(X.rows,X.cols,CV_8U); + X.convertTo(res, CV_8U, 255); + } + + void denoise_TVL1(const std::vector& observations,Mat& result, double lambda, int niters){ + + CV_Assert(observations.size()>0 && niters>0 && lambda>0); +#if 0 + solve_TVL1(observations[0],result,lambda,niters); + return; +#endif + + const float L2 = 8.0f, tau = 0.02f, sigma = 1./(L2*tau), theta = 1.f, img_scale = 1.f/255; + float clambda = (float)lambda, threshold = clambda*tau; + float s=0; + const int workdepth = CV_32F; + + int i, x, y, rows=observations[0].rows, cols=observations[0].cols,count; + for(i=1;i > Rs(observations.size()); + for(count=0;count(y); + const float* x_next = X.ptr(std::min(y+1, rows-1)); + Point2f* p_curr = P.ptr(y); + float dx, dy, m; + for( x = 0; x < cols-1; x++ ) + { + dx = (x_curr[x+1] - x_curr[x])*currsigma + p_curr[x].x; + dy = (x_next[x] - x_curr[x])*currsigma + p_curr[x].y; + m = 1.f/std::max(std::sqrt(dx*dx + dy*dy), 1.f); + p_curr[x].x = dx*m; + p_curr[x].y = dy*m; + } + dy = (x_next[x] - x_curr[x])*currsigma + p_curr[x].y; + m = 1.f/std::max(std::abs(dy), 1.f); + p_curr[x].x = 0.f; + p_curr[x].y = dy*m; + } + + + //Rs = clip(Rs + sigma*(X-imgs), -clambda, clambda) + for(count=0;count,MatConstIterator_,MatIterator_,AddFloatToCharScaled>( + Rs[count].begin(),Rs[count].end(),observations[count].begin(), + Rs[count].begin(),AddFloatToCharScaled(-sigma/255.0)); + Rs[count]+=sigma*X; + min(Rs[count],clambda,Rs[count]); + max(Rs[count],-clambda,Rs[count]); + } + + for( y = 0; y < rows; y++ ) + { + const uchar* img_curr = observations[0].ptr(y); + float* x_curr = X.ptr(y); + const Point2f* p_curr = P.ptr(y); + const Point2f* p_prev = P.ptr(std::max(y - 1, 0)); + + // X1 = X + tau*(-nablaT(P)) + x = 0; + s=0.0; + for(count=0;count(x,y),val,(val==image.at(x,y))?"true":"false"); + return (image.at(x,y)==val); +} + +TEST(Optim_denoise_tvl1, regression_basic){ + cv::RNG rng(42); + cv::Mat img = cv::imread("lena.jpg", 0), noisy,res; + if(img.rows!=512 || img.cols!=512){ + printf("\tplease, put lena.jpg from samples/c in the current folder\n"); + printf("\tnow, the test will fail...\n"); + ASSERT_TRUE(false); + } + + const int obs_num=5; + std::vector images(obs_num,cv::Mat()); + for(int i=0;i Date: Sun, 1 Sep 2013 01:02:06 +0800 Subject: [PATCH 2/4] Minor fixes As the opencv's build-bot did not want to compile this revision, I had to do some changes. In particular, 1) Removed unsigned int vs int comparisons, that were treated as errors 2) Removed unused variables and functions 3) Removed functions without previous declaration 4) Fixed whitespaces --- modules/optim/doc/primal_dual_algorithm.rst | 7 +- modules/optim/src/denoise_tvl1.cpp | 86 ++------------------- modules/optim/test/test_denoise_tvl1.cpp | 1 - 3 files changed, 10 insertions(+), 84 deletions(-) diff --git a/modules/optim/doc/primal_dual_algorithm.rst b/modules/optim/doc/primal_dual_algorithm.rst index 1eda00f596..09d736f242 100644 --- a/modules/optim/doc/primal_dual_algorithm.rst +++ b/modules/optim/doc/primal_dual_algorithm.rst @@ -12,7 +12,7 @@ problems (that is, finding a function to minimize some functional) problem, primal-dual algorithm then can be used to perform denoising and this is exactly what is implemented. -It should be noted, that this implementation was taken from the July 2013 blog entry [Mordvintsev]_, which also contained +It should be noted, that this implementation was taken from the July 2013 blog entry [Mordvintsev]_, which also contained (slightly more general) ready-to-use source code on Python. Subsequently, that code was rewritten on C++ with the usage of openCV by Vadim Pisarevsky at the end of July 2013 and finally it was slightly adapted by later authors. @@ -38,10 +38,7 @@ into play. :param result: Here the denoised image will be stored. There is no need to do pre-allocation of storage space, as it will be automatically allocated, if necessary. - :param lambda: Corresponds to :math:`\lambda` in the formulas above. As it is enlarged, the smooth (blurred) images are treated - more favorably than detailed (but maybe more noised) ones. Roughly speaking, as it becomes smaller, the result will be - more blur but more sever outliers will be removed. - + :param lambda: Corresponds to :math:`\lambda` in the formulas above. As it is enlarged, the smooth (blurred) images are treated more favorably than detailed (but maybe more noised) ones. Roughly speaking, as it becomes smaller, the result will be more blur but more sever outliers will be removed. :param niters: Number of iterations that the algorithm will run. Of course, as more iterations as better, but it is hard to quantitatively refine this statement, so just use the default and increase it if the results are poor. diff --git a/modules/optim/src/denoise_tvl1.cpp b/modules/optim/src/denoise_tvl1.cpp index 47f759affe..c8ba64bb66 100644 --- a/modules/optim/src/denoise_tvl1.cpp +++ b/modules/optim/src/denoise_tvl1.cpp @@ -1,5 +1,5 @@ #include "precomp.hpp" -#define ALEX_DEBUG +#undef ALEX_DEBUG #include "debug.hpp" #include #include @@ -18,93 +18,24 @@ namespace cv{namespace optim{ float _scale; }; - void solve_TVL1(const Mat& img, Mat& res, double _clambda, int niters) - { - const float L2 = 8.0f, tau = 0.02f, sigma = 1./(L2*tau), theta = 1.f, img_scale = 1.f/255; - float clambda = (float)_clambda, threshold = clambda*tau; - const int workdepth = CV_32F; - - int i, x, y, rows=img.rows, cols=img.cols; - - Mat X, P = Mat::zeros(rows, cols, CV_MAKETYPE(workdepth, 2)); - img.convertTo(X, workdepth, 1./255); - - for( i = 0; i < niters; i++ ) - { - float currsigma = i == 0 ? 1 + sigma : sigma; - - // P_ = P + sigma*nabla(X) - // P(x,y) = P_(x,y)/max(||P(x,y)||,1) - for( y = 0; y < rows; y++ ) - { - const float* x_curr = X.ptr(y); - const float* x_next = X.ptr(std::min(y+1, rows-1)); - Point2f* p_curr = P.ptr(y); - float dx, dy, m; - for( x = 0; x < cols-1; x++ ) - { - dx = (x_curr[x+1] - x_curr[x])*currsigma + p_curr[x].x; - dy = (x_next[x] - x_curr[x])*currsigma + p_curr[x].y; - m = 1.f/std::max(std::sqrt(dx*dx + dy*dy), 1.f); - p_curr[x].x = dx*m; - p_curr[x].y = dy*m; - } - dy = (x_next[x] - x_curr[x])*currsigma + p_curr[x].y; - m = 1.f/std::max(std::abs(dy), 1.f); - p_curr[x].x = 0.f; - p_curr[x].y = dy*m; - } - - // X1 = X + tau*(-nablaT(P)) - // X2 = X1 + clip(img - X1, -clambda*tau, clambda*tau) - // X = X2 + theta*(X2 - X) - for( y = 0; y < rows; y++ ) - { - const uchar* img_curr = img.ptr(y); - float* x_curr = X.ptr(y); - const Point2f* p_curr = P.ptr(y); - const Point2f* p_prev = P.ptr(std::max(y - 1, 0)); - - x = 0; - float x_new = x_curr[x] + tau*(p_curr[x].y - p_prev[x].y); - x_new += std::min(std::max(img_curr[x]*img_scale - x_new, -threshold), threshold); - x_curr[x] = x_new + theta*(x_new - x_curr[x]); - - for( x = 1; x < cols; x++ ) - { - x_new = x_curr[x] + tau*(p_curr[x].x - p_curr[x-1].x + p_curr[x].y - p_prev[x].y); - x_new += std::min(std::max(img_curr[x]*img_scale - x_new, -threshold), threshold); - x_curr[x] = x_new + theta*(x_new - x_curr[x]); - } - } - } - - res.create(X.rows,X.cols,CV_8U); - X.convertTo(res, CV_8U, 255); - } - void denoise_TVL1(const std::vector& observations,Mat& result, double lambda, int niters){ CV_Assert(observations.size()>0 && niters>0 && lambda>0); -#if 0 - solve_TVL1(observations[0],result,lambda,niters); - return; -#endif - const float L2 = 8.0f, tau = 0.02f, sigma = 1./(L2*tau), theta = 1.f, img_scale = 1.f/255; - float clambda = (float)lambda, threshold = clambda*tau; + const float L2 = 8.0f, tau = 0.02f, sigma = 1./(L2*tau), theta = 1.f; + float clambda = (float)lambda; float s=0; const int workdepth = CV_32F; int i, x, y, rows=observations[0].rows, cols=observations[0].cols,count; - for(i=1;i > Rs(observations.size()); - for(count=0;count,MatConstIterator_,MatIterator_,AddFloatToCharScaled>( Rs[count].begin(),Rs[count].end(),observations[count].begin(), Rs[count].begin(),AddFloatToCharScaled(-sigma/255.0)); @@ -147,7 +78,6 @@ namespace cv{namespace optim{ for( y = 0; y < rows; y++ ) { - const uchar* img_curr = observations[0].ptr(y); float* x_curr = X.ptr(y); const Point2f* p_curr = P.ptr(y); const Point2f* p_prev = P.ptr(std::max(y - 1, 0)); @@ -155,7 +85,7 @@ namespace cv{namespace optim{ // X1 = X + tau*(-nablaT(P)) x = 0; s=0.0; - for(count=0;count Date: Sun, 1 Sep 2013 07:11:31 +0800 Subject: [PATCH 3/4] Minor fixes Fixed integer vs unsigned integer comparison in .cpp test source and trailing whitespaces in source code --- modules/optim/src/denoise_tvl1.cpp | 4 ++-- modules/optim/test/test_denoise_tvl1.cpp | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/modules/optim/src/denoise_tvl1.cpp b/modules/optim/src/denoise_tvl1.cpp index c8ba64bb66..97435a231c 100644 --- a/modules/optim/src/denoise_tvl1.cpp +++ b/modules/optim/src/denoise_tvl1.cpp @@ -7,7 +7,7 @@ #define ABSCLIP(val,threshold) MIN(MAX((val),-(threshold)),(threshold)) namespace cv{namespace optim{ - + class AddFloatToCharScaled{ public: AddFloatToCharScaled(float scale):_scale(scale){} @@ -91,7 +91,7 @@ namespace cv{namespace optim{ float x_new = x_curr[x] + tau*(p_curr[x].y - p_prev[x].y)-tau*s; // X = X2 + theta*(X2 - X) x_curr[x] = x_new + theta*(x_new - x_curr[x]); - + for(x = 1; x < cols; x++ ) { diff --git a/modules/optim/test/test_denoise_tvl1.cpp b/modules/optim/test/test_denoise_tvl1.cpp index dfb1f7e69e..7ce541b7d6 100644 --- a/modules/optim/test/test_denoise_tvl1.cpp +++ b/modules/optim/test/test_denoise_tvl1.cpp @@ -41,7 +41,7 @@ TEST(Optim_denoise_tvl1, regression_basic){ const int obs_num=5; std::vector images(obs_num,cv::Mat()); - for(int i=0;i Date: Sun, 1 Sep 2013 13:59:15 +0800 Subject: [PATCH 4/4] Eliminate use of 32-bit floating pt type Replace all "float" by "double" (64-bit) to avoid "lose precision" warnings. --- modules/optim/src/denoise_tvl1.cpp | 44 ++++++++++++------------ modules/optim/test/test_denoise_tvl1.cpp | 2 +- 2 files changed, 23 insertions(+), 23 deletions(-) diff --git a/modules/optim/src/denoise_tvl1.cpp b/modules/optim/src/denoise_tvl1.cpp index 97435a231c..b11ebc0587 100644 --- a/modules/optim/src/denoise_tvl1.cpp +++ b/modules/optim/src/denoise_tvl1.cpp @@ -10,22 +10,22 @@ namespace cv{namespace optim{ class AddFloatToCharScaled{ public: - AddFloatToCharScaled(float scale):_scale(scale){} - inline float operator()(float a,uchar b){ - return a+_scale*((float)b); + AddFloatToCharScaled(double scale):_scale(scale){} + inline double operator()(double a,uchar b){ + return a+_scale*((double)b); } private: - float _scale; + double _scale; }; void denoise_TVL1(const std::vector& observations,Mat& result, double lambda, int niters){ CV_Assert(observations.size()>0 && niters>0 && lambda>0); - const float L2 = 8.0f, tau = 0.02f, sigma = 1./(L2*tau), theta = 1.f; - float clambda = (float)lambda; - float s=0; - const int workdepth = CV_32F; + const double L2 = 8.0, tau = 0.02, sigma = 1./(L2*tau), theta = 1.0; + double clambda = (double)lambda; + double s=0; + const int workdepth = CV_64F; int i, x, y, rows=observations[0].rows, cols=observations[0].cols,count; for(i=1;i<(int)observations.size();i++){ @@ -34,41 +34,41 @@ namespace cv{namespace optim{ Mat X, P = Mat::zeros(rows, cols, CV_MAKETYPE(workdepth, 2)); observations[0].convertTo(X, workdepth, 1./255); - std::vector< Mat_ > Rs(observations.size()); + std::vector< Mat_ > Rs(observations.size()); for(count=0;count<(int)Rs.size();count++){ Rs[count]=Mat::zeros(rows,cols,workdepth); } for( i = 0; i < niters; i++ ) { - float currsigma = i == 0 ? 1 + sigma : sigma; + double currsigma = i == 0 ? 1 + sigma : sigma; // P_ = P + sigma*nabla(X) // P(x,y) = P_(x,y)/max(||P(x,y)||,1) for( y = 0; y < rows; y++ ) { - const float* x_curr = X.ptr(y); - const float* x_next = X.ptr(std::min(y+1, rows-1)); - Point2f* p_curr = P.ptr(y); - float dx, dy, m; + const double* x_curr = X.ptr(y); + const double* x_next = X.ptr(std::min(y+1, rows-1)); + Point2d* p_curr = P.ptr(y); + double dx, dy, m; for( x = 0; x < cols-1; x++ ) { dx = (x_curr[x+1] - x_curr[x])*currsigma + p_curr[x].x; dy = (x_next[x] - x_curr[x])*currsigma + p_curr[x].y; - m = 1.f/std::max(std::sqrt(dx*dx + dy*dy), 1.f); + m = 1.0/std::max(std::sqrt(dx*dx + dy*dy), 1.0); p_curr[x].x = dx*m; p_curr[x].y = dy*m; } dy = (x_next[x] - x_curr[x])*currsigma + p_curr[x].y; - m = 1.f/std::max(std::abs(dy), 1.f); - p_curr[x].x = 0.f; + m = 1.0/std::max(std::abs(dy), 1.0); + p_curr[x].x = 0.0; p_curr[x].y = dy*m; } //Rs = clip(Rs + sigma*(X-imgs), -clambda, clambda) for(count=0;count<(int)Rs.size();count++){ - std::transform,MatConstIterator_,MatIterator_,AddFloatToCharScaled>( + std::transform,MatConstIterator_,MatIterator_,AddFloatToCharScaled>( Rs[count].begin(),Rs[count].end(),observations[count].begin(), Rs[count].begin(),AddFloatToCharScaled(-sigma/255.0)); Rs[count]+=sigma*X; @@ -78,9 +78,9 @@ namespace cv{namespace optim{ for( y = 0; y < rows; y++ ) { - float* x_curr = X.ptr(y); - const Point2f* p_curr = P.ptr(y); - const Point2f* p_prev = P.ptr(std::max(y - 1, 0)); + double* x_curr = X.ptr(y); + const Point2d* p_curr = P.ptr(y); + const Point2d* p_prev = P.ptr(std::max(y - 1, 0)); // X1 = X + tau*(-nablaT(P)) x = 0; @@ -88,7 +88,7 @@ namespace cv{namespace optim{ for(count=0;count<(int)Rs.size();count++){ s=s+Rs[count](y,x); } - float x_new = x_curr[x] + tau*(p_curr[x].y - p_prev[x].y)-tau*s; + double x_new = x_curr[x] + tau*(p_curr[x].y - p_prev[x].y)-tau*s; // X = X2 + theta*(X2 - X) x_curr[x] = x_new + theta*(x_new - x_curr[x]); diff --git a/modules/optim/test/test_denoise_tvl1.cpp b/modules/optim/test/test_denoise_tvl1.cpp index 7ce541b7d6..2721a7666d 100644 --- a/modules/optim/test/test_denoise_tvl1.cpp +++ b/modules/optim/test/test_denoise_tvl1.cpp @@ -8,7 +8,7 @@ void make_noisy(const cv::Mat& img, cv::Mat& noisy, double sigma, double pepper_ cv::addWeighted(img, 1, noise, 1, -128, noisy); cv::randn(noise, cv::Scalar::all(0), cv::Scalar::all(2)); noise *= 255; - cv::randu(mask, 0, round(1./pepper_salt_ratio)); + cv::randu(mask, 0, cvRound(1./pepper_salt_ratio)); cv::Mat half = mask.colRange(0, img.cols/2); half = cv::Scalar::all(1); noise.setTo(128, mask);