1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 1.2 +++ b/PTdecode/CImg-1.3.0/plugins/nlmeans.h Mon Aug 03 14:09:20 2009 +0100 1.3 @@ -0,0 +1,230 @@ 1.4 +/* 1.5 + # 1.6 + # File : nlmeans.h 1.7 + # ( C++ header file - CImg plug-in ) 1.8 + # 1.9 + # Description : CImg plugin that implements the non-local mean filter. 1.10 + # This file is a part of the CImg Library project. 1.11 + # ( http://cimg.sourceforge.net ) 1.12 + # 1.13 + # [1] Buades, A.; Coll, B.; Morel, J.-M.: A non-local algorithm for image denoising 1.14 + # IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2005. CVPR 2005. 1.15 + # Volume 2, 20-25 June 2005 Page(s):60 - 65 vol. 2 1.16 + # 1.17 + # [2] Buades, A. Coll, B. and Morel, J.: A review of image denoising algorithms, with a new one. 1.18 + # Multiscale Modeling and Simulation: A SIAM Interdisciplinary Journal 4 (2004) 490-530 1.19 + # 1.20 + # [3] Gasser, T. Sroka,L. Jennen Steinmetz,C. Residual variance and residual pattern nonlinear regression. 1.21 + # Biometrika 73 (1986) 625-659 1.22 + # 1.23 + # Copyright : Jerome Boulanger 1.24 + # ( http://www.irisa.fr/vista/Equipe/People/Jerome.Boulanger.html ) 1.25 + # 1.26 + # License : CeCILL v2.0 1.27 + # ( http://www.cecill.info/licences/Licence_CeCILL_V2-en.html ) 1.28 + # 1.29 + # This software is governed by the CeCILL license under French law and 1.30 + # abiding by the rules of distribution of free software. You can use, 1.31 + # modify and/ or redistribute the software under the terms of the CeCILL 1.32 + # license as circulated by CEA, CNRS and INRIA at the following URL 1.33 + # "http://www.cecill.info". 1.34 + # 1.35 + # As a counterpart to the access to the source code and rights to copy, 1.36 + # modify and redistribute granted by the license, users are provided only 1.37 + # with a limited warranty and the software's author, the holder of the 1.38 + # economic rights, and the successive licensors have only limited 1.39 + # liability. 1.40 + # 1.41 + # In this respect, the user's attention is drawn to the risks associated 1.42 + # with loading, using, modifying and/or developing or reproducing the 1.43 + # software by the user in light of its specific status of free software, 1.44 + # that may mean that it is complicated to manipulate, and that also 1.45 + # therefore means that it is reserved for developers and experienced 1.46 + # professionals having in-depth computer knowledge. Users are therefore 1.47 + # encouraged to load and test the software's suitability as regards their 1.48 + # requirements in conditions enabling the security of their systems and/or 1.49 + # data to be ensured and, more generally, to use and operate it in the 1.50 + # same conditions as regards security. 1.51 + # 1.52 + # The fact that you are presently reading this means that you have had 1.53 + # knowledge of the CeCILL license and that you accept its terms. 1.54 + # 1.55 +*/ 1.56 + 1.57 +#ifndef cimg_plugin_nlmeans 1.58 +#define cimg_plugin_nlmeans 1.59 + 1.60 +#include "noise_analysis.h" 1.61 + 1.62 +//! NL-Means denoising algorithm. 1.63 +/** 1.64 + This is the in-place version of get_nlmean(). 1.65 +**/ 1.66 +CImg<T>& nlmeans(int patch_size=1, double lambda=-1, double alpha=3, double sigma=-1, int sampling=1){ 1.67 + if (!is_empty()){ 1.68 + if (sigma<0) sigma=std::sqrt(noise_variance()); // noise variance estimation 1.69 + double np=(2*patch_size+1)*(2*patch_size+1)*dimv()/(double)sampling; 1.70 + if (lambda<0) {// Bandwidth estimation 1.71 + if (np<100) 1.72 + lambda =(((((( 1.1785e-12*np -5.1827e-10)*np+ 9.5946e-08)*np -9.7798e-06)*np+ 6.0756e-04)*np -0.0248)*np+ 1.9203)*np +7.9599; 1.73 + else 1.74 + lambda = (-7.2611e-04*np+ 1.3213)*np+ 15.2726; 1.75 + } 1.76 +#if cimg_debug>=1 1.77 + std::fprintf(stderr,"Size of the patch : %dx%d \n", 1.78 + 2*patch_size+1,2*patch_size+1); 1.79 + std::fprintf(stderr,"Size of window where similar patch are looked for : %dx%d \n", 1.80 + (int)(alpha*(2*patch_size+1)),(int)(alpha*(2*patch_size+1))); 1.81 + std::fprintf(stderr,"Bandwidth of the kernel : %fx%f^2 \n", 1.82 + lambda,sigma); 1.83 + std::fprintf(stderr,"Noise standard deviation estimated to : %f \n", 1.84 + sigma); 1.85 +#endif 1.86 + 1.87 + CImg<T> dest(dimx(),dimy(),dimz(),dimv(),0); 1.88 + double * uhat = new double[dimv()]; 1.89 + double h2=-.5/(lambda*sigma*sigma); // [Kervrann] notations 1.90 + if (dimz()!=1){// 3D case 1.91 + CImg<> P=(*this).get_blur(1); // inspired from Mahmoudi&Sapiro SPletter dec 05 1.92 + int n_simu=64; 1.93 + CImg<> tmp(n_simu,n_simu,n_simu); 1.94 + double sig = std::sqrt(tmp.fill(0.f).noise(sigma).blur(1).pow(2.).sum()/(n_simu*n_simu*n_simu)); 1.95 + int patch_size_z=0; 1.96 + int pxi=(int)(alpha*patch_size), 1.97 + pyi=(int)(alpha*patch_size), 1.98 + pzi=2;//Define the size of the neighborhood in z 1.99 + for (int zi=0;zi<dimz();zi++){ 1.100 +#if cimg_debug>=1 1.101 + std::fprintf(stderr,"\rProcessing : %3d %%",(int)((float)zi/(float)dimz()*100.));fflush(stdout); 1.102 +#endif 1.103 + for (int yi=0;yi<dimy();yi++) 1.104 + for (int xi=0;xi<dimx();xi++){ 1.105 + for (int v=0;v<dimv();v++) uhat[v] = 0; 1.106 + float sw=0,wmax=-1; 1.107 + for (int zj=cimg::max(0,zi-pzi);zj<cimg::min(dimz(),zi+pzi+1);zj++) 1.108 + for (int yj=cimg::max(0,yi-pyi);yj<cimg::min(dimy(),yi+pyi+1);yj++) 1.109 + for (int xj=cimg::max(0,xi-pxi);xj<cimg::min(dimx(),xi+pxi+1);xj++) 1.110 + if( cimg::abs(P(xi,yi,zi)-P(xj,yj,zj))/sig < 3){ 1.111 + double d = 0; 1.112 + int n = 0; 1.113 + if (xi!=xj && yi!=yj && zi!=zj){ 1.114 + for (int kz=-patch_size_z;kz<patch_size_z+1;kz+=sampling){ 1.115 + int zj_ = zj+kz; 1.116 + int zi_ = zi+kz; 1.117 + if (zj_>=0 && zj_<dimz() && zi_>=0 && zi_<dimz()) 1.118 + for (int ky=-patch_size;ky<=patch_size;ky+=sampling){ 1.119 + int yj_ = yj+ky; 1.120 + int yi_ = yi+ky; 1.121 + if (yj_>=0 && yj_<dimy() && yi_>=0 && yi_<dimy()) 1.122 + for (int kx=-patch_size;kx<=patch_size;kx+=sampling){ 1.123 + int xj_ = xj+kx; 1.124 + int xi_ = xi+kx; 1.125 + if (xj_>=0 && xj_<dimx() && xi_>=0 && xi_<dimx()) 1.126 + for (int v=0;v<dimv();v++){ 1.127 + double d1 = (*this)(xj_,yj_,zj_,v)-(*this)(xi_,yi_,zi_,v); 1.128 + d += d1*d1; 1.129 + n++; 1.130 + } 1.131 + } 1.132 + } 1.133 + } 1.134 + float w = (float)std::exp(d*h2); 1.135 + wmax = w>wmax?w:wmax; 1.136 + for (int v=0;v<dimv();v++) uhat[v] += w*(*this)(xj,yj,zj,v); 1.137 + sw += w; 1.138 + } 1.139 + } 1.140 + // add the central pixel 1.141 + { for (int v=0;v<dimv();v++) uhat[v] += wmax*(*this)(xi,yi,zi,v); } 1.142 + sw += wmax; 1.143 + { for (int v=0;v<dimv();v++) dest(xi,yi,zi,v)= (T) (uhat[v] /= sw); } 1.144 + } 1.145 + } 1.146 + } 1.147 + else{ // 2D case 1.148 + CImg<> P=(*this).get_blur(1); // inspired from Mahmoudi&Sapiro SPletter dec 05 1.149 + int n_simu=512; 1.150 + CImg<> tmp(n_simu,n_simu); 1.151 + double sig = std::sqrt(tmp.fill(0.f).noise(sigma).blur(1).pow(2.).sum()/(n_simu*n_simu)); 1.152 + int pxi=(int)(alpha*patch_size),pyi=(int)(alpha*patch_size);//Define the size of the neighborhood 1.153 + for (int yi=0;yi<dimy();yi++){ 1.154 +#if cimg_debug>=1 1.155 + std::fprintf(stderr,"\rProcessing : %3d %%",(int)((float)yi/(float)dimy()*100.));fflush(stdout); 1.156 +#endif 1.157 + for (int xi=0;xi<dimx();xi++){ 1.158 + for (int v=0;v<dimv();v++) uhat[v] = 0; 1.159 + float sw=0,wmax=-1; 1.160 + for (int yj=cimg::max(0,yi-pyi);yj<cimg::min(dimy(),yi+pyi+1);yj++) 1.161 + for (int xj=cimg::max(0,xi-pxi);xj<cimg::min(dimx(),xi+pxi+1);xj++) 1.162 + if( cimg::abs(P(xi,yi)-P(xj,yj))/sig < 3.){ 1.163 + double d = 0; 1.164 + int n = 0; 1.165 + if (!(xi==xj && yi==yj)) 1.166 + for (int ky=-patch_size;ky<patch_size+1;ky+=sampling){ 1.167 + int yj_ = yj+ky; 1.168 + int yi_ = yi+ky; 1.169 + if (yj_>=0 && yj_<dimy() && yi_>=0 && yi_<dimy()) 1.170 + for (int kx=-patch_size;kx<patch_size+1;kx+=sampling){ 1.171 + int xj_ = xj+kx; 1.172 + int xi_ = xi+kx; 1.173 + if (xj_>=0 && xj_<dimx() && xi_>=0 && xi_<dimx()) 1.174 + for (int v=0;v<dimv();v++){ 1.175 + double d1 = (*this)(xj_,yj_,v)-(*this)(xi_,yi_,v); 1.176 + d += d1*d1; 1.177 + n++; 1.178 + } 1.179 + } 1.180 + } 1.181 + float w=(float)std::exp(d*h2); 1.182 + for (int v=0;v<dimv();v++) uhat[v] += w*(*this)(xj,yj,v); 1.183 + wmax = w>wmax?w:wmax; // Store the maximum of the weights 1.184 + sw += w; // Compute the sum of the weights 1.185 + } 1.186 + // add the central pixel with the maximum weight 1.187 + { for (int v=0;v<dimv();v++) uhat[v] += wmax*(*this)(xi,yi,v); } 1.188 + sw += wmax; 1.189 + // Compute the estimate for the current pixel 1.190 + { for (int v=0;v<dimv();v++) dest(xi,yi,v)= (T) (uhat[v] /= sw); } 1.191 + } 1.192 + }// main loop 1.193 + }// 2d 1.194 + delete [] uhat; 1.195 + *this=dest; 1.196 +#if cimg_debug>=1 1.197 + std::fprintf(stderr,"\n"); // make a new line 1.198 +#endif 1.199 + }// is empty 1.200 + return *this; 1.201 +} 1.202 + 1.203 +//! Get the result of the NL-Means denoising algorithm. 1.204 +/** 1.205 + \param patch_size = radius of the patch (1=3x3 by default) 1.206 + \param lambda = bandwidth ( -1 by default : automatic selection) 1.207 + \param alpha = size of the region where similar patch are searched (3 x patch_size = 9x9 by default) 1.208 + \param sigma = noise standard deviation (-1 = estimation) 1.209 + \param sampling = sampling of the patch (1 = uses all point, 2 = uses one point on 4, etc) 1.210 + If the image has three dimensions then the patch is only in 2D and the neighborhood extent in time is only 5. 1.211 + If the image has several channel (color images), the distance between the two patch is computed using 1.212 + all the channels. 1.213 + The greater the patch is the best is the result. 1.214 + Lambda parameter is function of the size of the patch size. The automatic Lambda parameter is taken 1.215 + in the Chi2 table at a significiance level of 0.01. This diffear from the original paper [1]. The weighted average becomes then: 1.216 + \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 1.217 + where \f$ P(x,y) $\f denotes the patch in (x,y) location. 1.218 + 1.219 + An a priori is also used to increase the speed of the algorithm in the spirit of Sapiro et al. SPletter dec 05 1.220 + 1.221 + This very basic version of the Non-Local Means algorithm provides an output image which contains 1.222 + some residual noise with a relatively small variance (\f$\sigma<5$\f). 1.223 + 1.224 + [1] A non-local algorithm for image denoising 1.225 + Buades, A.; Coll, B.; Morel, J.-M.; 1.226 + Computer Vision and Pattern Recognition, 2005. CVPR 2005. IEEE Computer Society Conference on 1.227 + Volume 2, 20-25 June 2005 Page(s):60 - 65 vol. 2 1.228 +**/ 1.229 +CImg<T> get_nlmeans( int patch_size=1, double lambda=-1, double alpha=3 ,double sigma=-1, int sampling=1) { 1.230 + return CImg<T>(*this).nlmeans(patch_size,lambda,alpha,sigma,sampling); 1.231 +} 1.232 + 1.233 +#endif