1 /* 2 # 3 # File : nlmeans.h 4 # ( C++ header file - CImg plug-in ) 5 # 6 # Description : CImg plugin that implements the non-local mean filter. 7 # This file is a part of the CImg Library project. 8 # ( http://cimg.eu ) 9 # 10 # [1] Buades, A.; Coll, B.; Morel, J.-M.: A non-local algorithm for image denoising 11 # IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2005. CVPR 2005. 12 # Volume 2, 20-25 June 2005 Page(s):60 - 65 vol. 2 13 # 14 # [2] Buades, A. Coll, B. and Morel, J.: A review of image denoising algorithms, with a new one. 15 # Multiscale Modeling and Simulation: A SIAM Interdisciplinary Journal 4 (2004) 490-530 16 # 17 # [3] Gasser, T. Sroka,L. Jennen Steinmetz,C. Residual variance and residual pattern nonlinear regression. 18 # Biometrika 73 (1986) 625-659 19 # 20 # Copyright : Jerome Boulanger 21 # ( http://www.irisa.fr/vista/Equipe/People/Jerome.Boulanger.html ) 22 # 23 # License : CeCILL v2.0 24 # ( http://www.cecill.info/licences/Licence_CeCILL_V2-en.html ) 25 # 26 # This software is governed by the CeCILL license under French law and 27 # abiding by the rules of distribution of free software. You can use, 28 # modify and/ or redistribute the software under the terms of the CeCILL 29 # license as circulated by CEA, CNRS and INRIA at the following URL 30 # "http://www.cecill.info". 31 # 32 # As a counterpart to the access to the source code and rights to copy, 33 # modify and redistribute granted by the license, users are provided only 34 # with a limited warranty and the software's author, the holder of the 35 # economic rights, and the successive licensors have only limited 36 # liability. 37 # 38 # In this respect, the user's attention is drawn to the risks associated 39 # with loading, using, modifying and/or developing or reproducing the 40 # software by the user in light of its specific status of free software, 41 # that may mean that it is complicated to manipulate, and that also 42 # therefore means that it is reserved for developers and experienced 43 # professionals having in-depth computer knowledge. Users are therefore 44 # encouraged to load and test the software's suitability as regards their 45 # requirements in conditions enabling the security of their systems and/or 46 # data to be ensured and, more generally, to use and operate it in the 47 # same conditions as regards security. 48 # 49 # The fact that you are presently reading this means that you have had 50 # knowledge of the CeCILL license and that you accept its terms. 51 # 52 */ 53 54 #ifndef cimg_plugin_nlmeans 55 #define cimg_plugin_nlmeans 56 57 //! NL-Means denoising algorithm. 58 /** 59 This is the in-place version of get_nlmean(). 60 **/ 61 CImg<T>& nlmeans(int patch_size=1, double lambda=-1, double alpha=3, double sigma=-1, int sampling=1){ 62 if (!is_empty()){ 63 if (sigma<0) sigma = std::sqrt(variance_noise()); // noise variance estimation 64 const double np = (2*patch_size + 1)*(2*patch_size + 1)*spectrum()/(double)sampling; 65 if (lambda<0) {// Bandwidth estimation 66 if (np<100) 67 lambda = ((((((1.1785e-12*np - 5.1827e-10)*np + 9.5946e-08)*np - 68 9.7798e-06)*np + 6.0756e-04)*np - 0.0248)*np + 1.9203)*np + 7.9599; 69 else 70 lambda = (-7.2611e-04*np + 1.3213)*np + 15.2726; 71 } 72 #if cimg_debug>=1 73 std::fprintf(stderr,"Size of the patch : %dx%d \n", 74 2*patch_size + 1,2*patch_size + 1); 75 std::fprintf(stderr,"Size of window where similar patch are looked for : %dx%d \n", 76 (int)(alpha*(2*patch_size + 1)),(int)(alpha*(2*patch_size + 1))); 77 std::fprintf(stderr,"Bandwidth of the kernel : %fx%f^2 \n", 78 lambda,sigma); 79 std::fprintf(stderr,"Noise standard deviation estimated to : %f \n", 80 sigma); 81 #endif 82 83 CImg<T> dest(width(),height(),depth(),spectrum(),0); 84 double *uhat = new double[spectrum()]; 85 const double h2 = -.5/(lambda*sigma*sigma); // [Kervrann] notations 86 if (depth()!=1){ // 3D case 87 const CImg<> P = (*this).get_blur(1); // inspired from Mahmoudi&Sapiro SPletter dec 05 88 const int n_simu = 64; 89 CImg<> tmp(n_simu,n_simu,n_simu); 90 const double sig = std::sqrt(tmp.fill(0.f).noise(sigma).blur(1).pow(2.).sum()/(n_simu*n_simu*n_simu)); 91 const int 92 patch_size_z = 0, 93 pxi = (int)(alpha*patch_size), 94 pyi = (int)(alpha*patch_size), 95 pzi = 2; //Define the size of the neighborhood in z 96 for (int zi = 0; zi<depth(); ++zi) { 97 #if cimg_debug>=1 98 std::fprintf(stderr,"\rProcessing : %3d %%",(int)((float)zi/(float)depth()*100.));fflush(stdout); 99 #endif 100 for (int yi = 0; yi<height(); ++yi) 101 for (int xi = 0; xi<width(); ++xi) { 102 cimg_forC(*this,v) uhat[v] = 0; 103 float sw = 0, wmax = -1; 104 for (int zj = std::max(0,zi - pzi); zj<std::min(depth(),zi + pzi + 1); ++zj) 105 for (int yj = std::max(0,yi - pyi); yj<std::min(height(),yi + pyi + 1); ++yj) 106 for (int xj = std::max(0,xi - pxi); xj<std::min(width(),xi + pxi + 1); ++xj) 107 if (cimg::abs(P(xi,yi,zi) - P(xj,yj,zj))/sig<3) { 108 double d = 0; 109 int n = 0; 110 if (xi!=xj && yi!=yj && zi!=zj){ 111 for (int kz = -patch_size_z; kz<patch_size_z + 1; kz+=sampling) { 112 int 113 zj_ = zj + kz, 114 zi_ = zi + kz; 115 if (zj_>=0 && zj_<depth() && zi_>=0 && zi_<depth()) 116 for (int ky = -patch_size; ky<=patch_size; ky+=sampling) { 117 int 118 yj_ = yj + ky, 119 yi_ = yi + ky; 120 if (yj_>=0 && yj_<height() && yi_>=0 && yi_<height()) 121 for (int kx = -patch_size; kx<=patch_size; kx+=sampling) { 122 int 123 xj_ = xj + kx, 124 xi_ = xi + kx; 125 if (xj_>=0 && xj_<width() && xi_>=0 && xi_<width()) 126 cimg_forC(*this,v) { 127 double d1 = (*this)(xj_,yj_,zj_,v) - (*this)(xi_,yi_,zi_,v); 128 d+=d1*d1; 129 ++n; 130 } 131 } 132 } 133 } 134 float w = (float)std::exp(d*h2); 135 wmax = w>wmax?w:wmax; 136 cimg_forC(*this,v) uhat[v]+=w*(*this)(xj,yj,zj,v); 137 sw+=w; 138 } 139 } 140 // add the central pixel 141 cimg_forC(*this,v) uhat[v]+=wmax*(*this)(xi,yi,zi,v); 142 sw+=wmax; 143 if (sw) cimg_forC(*this,v) dest(xi,yi,zi,v) = (T)(uhat[v]/=sw); 144 else cimg_forC(*this,v) dest(xi,yi,zi,v) = (*this)(xi,yi,zi,v); 145 } 146 } 147 } 148 else { // 2D case 149 const CImg<> P = (*this).get_blur(1); // inspired from Mahmoudi&Sapiro SPletter dec 05 150 const int n_simu = 512; 151 CImg<> tmp(n_simu,n_simu); 152 const double sig = std::sqrt(tmp.fill(0.f).noise(sigma).blur(1).pow(2.).sum()/(n_simu*n_simu)); 153 const int 154 pxi = (int)(alpha*patch_size), 155 pyi = (int)(alpha*patch_size); //Define the size of the neighborhood 156 for (int yi = 0; yi<height(); ++yi) { 157 #if cimg_debug>=1 158 std::fprintf(stderr,"\rProcessing : %3d %%",(int)((float)yi/(float)height()*100.));fflush(stdout); 159 #endif 160 for (int xi = 0; xi<width(); ++xi) { 161 cimg_forC(*this,v) uhat[v] = 0; 162 float sw = 0, wmax = -1; 163 for (int yj = std::max(0,yi - pyi); yj<std::min(height(),yi + pyi + 1); ++yj) 164 for (int xj = std::max(0,xi - pxi); xj<std::min(width(),xi + pxi + 1); ++xj) 165 if (cimg::abs(P(xi,yi) - P(xj,yj))/sig<3.) { 166 double d = 0; 167 int n = 0; 168 if (!(xi==xj && yi==yj)) //{ 169 for (int ky = -patch_size; ky<patch_size + 1; ky+=sampling) { 170 int 171 yj_ = yj + ky, 172 yi_ = yi + ky; 173 if (yj_>=0 && yj_<height() && yi_>=0 && yi_<height()) 174 for (int kx = -patch_size; kx<patch_size + 1; kx+=sampling) { 175 int 176 xj_ = xj + kx, 177 xi_ = xi + kx; 178 if (xj_>=0 && xj_<width() && xi_>=0 && xi_<width()) 179 cimg_forC(*this,v) { 180 double d1 = (*this)(xj_,yj_,v) - (*this)(xi_,yi_,v); 181 d+=d1*d1; 182 n++; 183 } 184 } 185 //} 186 float w = (float)std::exp(d*h2); 187 cimg_forC(*this,v) uhat[v]+=w*(*this)(xj,yj,v); 188 wmax = w>wmax?w:wmax; // Store the maximum of the weights 189 sw+=w; // Compute the sum of the weights 190 } 191 } 192 // add the central pixel with the maximum weight 193 cimg_forC(*this,v) uhat[v]+=wmax*(*this)(xi,yi,v); 194 sw+=wmax; 195 196 // Compute the estimate for the current pixel 197 if (sw) cimg_forC(*this,v) dest(xi,yi,v) = (T)(uhat[v]/=sw); 198 else cimg_forC(*this,v) dest(xi,yi,v) = (*this)(xi,yi,v); 199 } 200 } // main loop 201 } // 2d 202 delete [] uhat; 203 dest.move_to(*this); 204 #if cimg_debug>=1 205 std::fprintf(stderr,"\n"); // make a new line 206 #endif 207 } // is empty 208 return *this; 209 } 210 211 //! Get the result of the NL-Means denoising algorithm. 212 /** 213 \param patch_size = radius of the patch (1=3x3 by default) 214 \param lambda = bandwidth ( -1 by default : automatic selection) 215 \param alpha = size of the region where similar patch are searched (3 x patch_size = 9x9 by default) 216 \param sigma = noise standard deviation (-1 = estimation) 217 \param sampling = sampling of the patch (1 = uses all point, 2 = uses one point on 4, etc) 218 If the image has three dimensions then the patch is only in 2D and the neighborhood extent in time is only 5. 219 If the image has several channel (color images), the distance between the two patch is computed using 220 all the channels. 221 The greater the patch is the best is the result. 222 Lambda parameter is function of the size of the patch size. The automatic Lambda parameter is taken 223 in the Chi2 table at a significiance level of 0.01. This diffear from the original paper [1]. 224 The weighted average becomes then: 225 \f$$ \hat{f}(x,y) = \sum_{x',y'} \frac{1}{Z} exp(\frac{P(x,y)-P(x',y')}{2 \lambda \sigma^2}) f(x',y') $$\f 226 where \f$ P(x,y) $\f denotes the patch in (x,y) location. 227 228 An a priori is also used to increase the speed of the algorithm in the spirit of Sapiro et al. SPletter dec 05 229 230 This very basic version of the Non-Local Means algorithm provides an output image which contains 231 some residual noise with a relatively small variance (\f$\sigma<5$\f). 232 233 [1] A non-local algorithm for image denoising 234 Buades, A.; Coll, B.; Morel, J.-M.; 235 Computer Vision and Pattern Recognition, 2005. CVPR 2005. IEEE Computer Society Conference on 236 Volume 2, 20-25 June 2005 Page(s):60 - 65 vol. 2 237 **/ 238 CImg<T> get_nlmeans( int patch_size=1, double lambda=-1, double alpha=3 ,double sigma=-1, int sampling=1) const { 239 return CImg<T>(*this).nlmeans(patch_size,lambda,alpha,sigma,sampling); 240 } 241 242 #endif /* cimg_plugin_nlmeans */ 243