1 /*M/////////////////////////////////////////////////////////////////////////////////////// 2 // 3 // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. 4 // 5 // By downloading, copying, installing or using the software you agree to this license. 6 // If you do not agree to this license, do not download, install, 7 // copy or use the software. 8 // 9 // 10 // License Agreement 11 // For Open Source Computer Vision Library 12 // 13 // Copyright (C) 2013, OpenCV Foundation, all rights reserved. 14 // Third party copyrights are property of their respective owners. 15 // 16 // Redistribution and use in source and binary forms, with or without modification, 17 // are permitted provided that the following conditions are met: 18 // 19 // * Redistribution's of source code must retain the above copyright notice, 20 // this list of conditions and the following disclaimer. 21 // 22 // * Redistribution's in binary form must reproduce the above copyright notice, 23 // this list of conditions and the following disclaimer in the documentation 24 // and/or other materials provided with the distribution. 25 // 26 // * The name of the copyright holders may not be used to endorse or promote products 27 // derived from this software without specific prior written permission. 28 // 29 // This software is provided by the copyright holders and contributors "as is" and 30 // any express or implied warranties, including, but not limited to, the implied 31 // warranties of merchantability and fitness for a particular purpose are disclaimed. 32 // In no event shall the OpenCV Foundation or contributors be liable for any direct, 33 // indirect, incidental, special, exemplary, or consequential damages 34 // (including, but not limited to, procurement of substitute goods or services; 35 // loss of use, data, or profits; or business interruption) however caused 36 // and on any theory of liability, whether in contract, strict liability, 37 // or tort (including negligence or otherwise) arising in any way out of 38 // the use of this software, even if advised of the possibility of such damage. 39 // 40 //M*/ 41 #include "precomp.hpp" 42 #include <vector> 43 #include <algorithm> 44 45 #define ABSCLIP(val,threshold) MIN(MAX((val),-(threshold)),(threshold)) 46 47 namespace cv{ 48 49 class AddFloatToCharScaled{ 50 public: 51 AddFloatToCharScaled(double scale):_scale(scale){} 52 inline double operator()(double a,uchar b){ 53 return a+_scale*((double)b); 54 } 55 private: 56 double _scale; 57 }; 58 59 #ifndef OPENCV_NOSTL 60 using std::transform; 61 #else 62 template <class InputIterator, class InputIterator2, class OutputIterator, class BinaryOperator> 63 static OutputIterator transform (InputIterator first1, InputIterator last1, InputIterator2 first2, 64 OutputIterator result, BinaryOperator binary_op) 65 { 66 while (first1 != last1) 67 { 68 *result = binary_op(*first1, *first2); 69 ++result; ++first1; ++first2; 70 } 71 return result; 72 } 73 #endif 74 void denoise_TVL1(const std::vector<Mat>& observations,Mat& result, double lambda, int niters){ 75 76 CV_Assert(observations.size()>0 && niters>0 && lambda>0); 77 78 const double L2 = 8.0, tau = 0.02, sigma = 1./(L2*tau), theta = 1.0; 79 double clambda = (double)lambda; 80 double s=0; 81 const int workdepth = CV_64F; 82 83 int i, x, y, rows=observations[0].rows, cols=observations[0].cols,count; 84 for(i=1;i<(int)observations.size();i++){ 85 CV_Assert(observations[i].rows==rows && observations[i].cols==cols); 86 } 87 88 Mat X, P = Mat::zeros(rows, cols, CV_MAKETYPE(workdepth, 2)); 89 observations[0].convertTo(X, workdepth, 1./255); 90 std::vector< Mat_<double> > Rs(observations.size()); 91 for(count=0;count<(int)Rs.size();count++){ 92 Rs[count]=Mat::zeros(rows,cols,workdepth); 93 } 94 95 for( i = 0; i < niters; i++ ) 96 { 97 double currsigma = i == 0 ? 1 + sigma : sigma; 98 99 // P_ = P + sigma*nabla(X) 100 // P(x,y) = P_(x,y)/max(||P(x,y)||,1) 101 for( y = 0; y < rows; y++ ) 102 { 103 const double* x_curr = X.ptr<double>(y); 104 const double* x_next = X.ptr<double>(std::min(y+1, rows-1)); 105 Point2d* p_curr = P.ptr<Point2d>(y); 106 double dx, dy, m; 107 for( x = 0; x < cols-1; x++ ) 108 { 109 dx = (x_curr[x+1] - x_curr[x])*currsigma + p_curr[x].x; 110 dy = (x_next[x] - x_curr[x])*currsigma + p_curr[x].y; 111 m = 1.0/std::max(std::sqrt(dx*dx + dy*dy), 1.0); 112 p_curr[x].x = dx*m; 113 p_curr[x].y = dy*m; 114 } 115 dy = (x_next[x] - x_curr[x])*currsigma + p_curr[x].y; 116 m = 1.0/std::max(std::abs(dy), 1.0); 117 p_curr[x].x = 0.0; 118 p_curr[x].y = dy*m; 119 } 120 121 122 //Rs = clip(Rs + sigma*(X-imgs), -clambda, clambda) 123 for(count=0;count<(int)Rs.size();count++){ 124 transform<MatIterator_<double>,MatConstIterator_<uchar>,MatIterator_<double>,AddFloatToCharScaled>( 125 Rs[count].begin(),Rs[count].end(),observations[count].begin<uchar>(), 126 Rs[count].begin(),AddFloatToCharScaled(-sigma/255.0)); 127 Rs[count]+=sigma*X; 128 min(Rs[count],clambda,Rs[count]); 129 max(Rs[count],-clambda,Rs[count]); 130 } 131 132 for( y = 0; y < rows; y++ ) 133 { 134 double* x_curr = X.ptr<double>(y); 135 const Point2d* p_curr = P.ptr<Point2d>(y); 136 const Point2d* p_prev = P.ptr<Point2d>(std::max(y - 1, 0)); 137 138 // X1 = X + tau*(-nablaT(P)) 139 x = 0; 140 s=0.0; 141 for(count=0;count<(int)Rs.size();count++){ 142 s=s+Rs[count](y,x); 143 } 144 double x_new = x_curr[x] + tau*(p_curr[x].y - p_prev[x].y)-tau*s; 145 // X = X2 + theta*(X2 - X) 146 x_curr[x] = x_new + theta*(x_new - x_curr[x]); 147 148 149 for(x = 1; x < cols; x++ ) 150 { 151 s=0.0; 152 for(count=0;count<(int)Rs.size();count++){ 153 s+=Rs[count](y,x); 154 } 155 // X1 = X + tau*(-nablaT(P)) 156 x_new = x_curr[x] + tau*(p_curr[x].x - p_curr[x-1].x + p_curr[x].y - p_prev[x].y)-tau*s; 157 // X = X2 + theta*(X2 - X) 158 x_curr[x] = x_new + theta*(x_new - x_curr[x]); 159 } 160 } 161 } 162 163 result.create(X.rows,X.cols,CV_8U); 164 X.convertTo(result, CV_8U, 255); 165 } 166 } 167
