Mon, 03 Aug 2009 23:41:04 +0100
added dep/*.d and obj/*.o to hgignore
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.sourceforge.net )
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 */
54 #ifndef cimg_plugin_nlmeans
55 #define cimg_plugin_nlmeans
57 #include "noise_analysis.h"
59 //! NL-Means denoising algorithm.
60 /**
61 This is the in-place version of get_nlmean().
62 **/
63 CImg<T>& nlmeans(int patch_size=1, double lambda=-1, double alpha=3, double sigma=-1, int sampling=1){
64 if (!is_empty()){
65 if (sigma<0) sigma=std::sqrt(noise_variance()); // noise variance estimation
66 double np=(2*patch_size+1)*(2*patch_size+1)*dimv()/(double)sampling;
67 if (lambda<0) {// Bandwidth estimation
68 if (np<100)
69 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;
70 else
71 lambda = (-7.2611e-04*np+ 1.3213)*np+ 15.2726;
72 }
73 #if cimg_debug>=1
74 std::fprintf(stderr,"Size of the patch : %dx%d \n",
75 2*patch_size+1,2*patch_size+1);
76 std::fprintf(stderr,"Size of window where similar patch are looked for : %dx%d \n",
77 (int)(alpha*(2*patch_size+1)),(int)(alpha*(2*patch_size+1)));
78 std::fprintf(stderr,"Bandwidth of the kernel : %fx%f^2 \n",
79 lambda,sigma);
80 std::fprintf(stderr,"Noise standard deviation estimated to : %f \n",
81 sigma);
82 #endif
84 CImg<T> dest(dimx(),dimy(),dimz(),dimv(),0);
85 double * uhat = new double[dimv()];
86 double h2=-.5/(lambda*sigma*sigma); // [Kervrann] notations
87 if (dimz()!=1){// 3D case
88 CImg<> P=(*this).get_blur(1); // inspired from Mahmoudi&Sapiro SPletter dec 05
89 int n_simu=64;
90 CImg<> tmp(n_simu,n_simu,n_simu);
91 double sig = std::sqrt(tmp.fill(0.f).noise(sigma).blur(1).pow(2.).sum()/(n_simu*n_simu*n_simu));
92 int patch_size_z=0;
93 int 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<dimz();zi++){
97 #if cimg_debug>=1
98 std::fprintf(stderr,"\rProcessing : %3d %%",(int)((float)zi/(float)dimz()*100.));fflush(stdout);
99 #endif
100 for (int yi=0;yi<dimy();yi++)
101 for (int xi=0;xi<dimx();xi++){
102 for (int v=0;v<dimv();v++) uhat[v] = 0;
103 float sw=0,wmax=-1;
104 for (int zj=cimg::max(0,zi-pzi);zj<cimg::min(dimz(),zi+pzi+1);zj++)
105 for (int yj=cimg::max(0,yi-pyi);yj<cimg::min(dimy(),yi+pyi+1);yj++)
106 for (int xj=cimg::max(0,xi-pxi);xj<cimg::min(dimx(),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 zj_ = zj+kz;
113 int zi_ = zi+kz;
114 if (zj_>=0 && zj_<dimz() && zi_>=0 && zi_<dimz())
115 for (int ky=-patch_size;ky<=patch_size;ky+=sampling){
116 int yj_ = yj+ky;
117 int yi_ = yi+ky;
118 if (yj_>=0 && yj_<dimy() && yi_>=0 && yi_<dimy())
119 for (int kx=-patch_size;kx<=patch_size;kx+=sampling){
120 int xj_ = xj+kx;
121 int xi_ = xi+kx;
122 if (xj_>=0 && xj_<dimx() && xi_>=0 && xi_<dimx())
123 for (int v=0;v<dimv();v++){
124 double d1 = (*this)(xj_,yj_,zj_,v)-(*this)(xi_,yi_,zi_,v);
125 d += d1*d1;
126 n++;
127 }
128 }
129 }
130 }
131 float w = (float)std::exp(d*h2);
132 wmax = w>wmax?w:wmax;
133 for (int v=0;v<dimv();v++) uhat[v] += w*(*this)(xj,yj,zj,v);
134 sw += w;
135 }
136 }
137 // add the central pixel
138 { for (int v=0;v<dimv();v++) uhat[v] += wmax*(*this)(xi,yi,zi,v); }
139 sw += wmax;
140 { for (int v=0;v<dimv();v++) dest(xi,yi,zi,v)= (T) (uhat[v] /= sw); }
141 }
142 }
143 }
144 else{ // 2D case
145 CImg<> P=(*this).get_blur(1); // inspired from Mahmoudi&Sapiro SPletter dec 05
146 int n_simu=512;
147 CImg<> tmp(n_simu,n_simu);
148 double sig = std::sqrt(tmp.fill(0.f).noise(sigma).blur(1).pow(2.).sum()/(n_simu*n_simu));
149 int pxi=(int)(alpha*patch_size),pyi=(int)(alpha*patch_size);//Define the size of the neighborhood
150 for (int yi=0;yi<dimy();yi++){
151 #if cimg_debug>=1
152 std::fprintf(stderr,"\rProcessing : %3d %%",(int)((float)yi/(float)dimy()*100.));fflush(stdout);
153 #endif
154 for (int xi=0;xi<dimx();xi++){
155 for (int v=0;v<dimv();v++) uhat[v] = 0;
156 float sw=0,wmax=-1;
157 for (int yj=cimg::max(0,yi-pyi);yj<cimg::min(dimy(),yi+pyi+1);yj++)
158 for (int xj=cimg::max(0,xi-pxi);xj<cimg::min(dimx(),xi+pxi+1);xj++)
159 if( cimg::abs(P(xi,yi)-P(xj,yj))/sig < 3.){
160 double d = 0;
161 int n = 0;
162 if (!(xi==xj && yi==yj))
163 for (int ky=-patch_size;ky<patch_size+1;ky+=sampling){
164 int yj_ = yj+ky;
165 int yi_ = yi+ky;
166 if (yj_>=0 && yj_<dimy() && yi_>=0 && yi_<dimy())
167 for (int kx=-patch_size;kx<patch_size+1;kx+=sampling){
168 int xj_ = xj+kx;
169 int xi_ = xi+kx;
170 if (xj_>=0 && xj_<dimx() && xi_>=0 && xi_<dimx())
171 for (int v=0;v<dimv();v++){
172 double d1 = (*this)(xj_,yj_,v)-(*this)(xi_,yi_,v);
173 d += d1*d1;
174 n++;
175 }
176 }
177 }
178 float w=(float)std::exp(d*h2);
179 for (int v=0;v<dimv();v++) uhat[v] += w*(*this)(xj,yj,v);
180 wmax = w>wmax?w:wmax; // Store the maximum of the weights
181 sw += w; // Compute the sum of the weights
182 }
183 // add the central pixel with the maximum weight
184 { for (int v=0;v<dimv();v++) uhat[v] += wmax*(*this)(xi,yi,v); }
185 sw += wmax;
186 // Compute the estimate for the current pixel
187 { for (int v=0;v<dimv();v++) dest(xi,yi,v)= (T) (uhat[v] /= sw); }
188 }
189 }// main loop
190 }// 2d
191 delete [] uhat;
192 *this=dest;
193 #if cimg_debug>=1
194 std::fprintf(stderr,"\n"); // make a new line
195 #endif
196 }// is empty
197 return *this;
198 }
200 //! Get the result of the NL-Means denoising algorithm.
201 /**
202 \param patch_size = radius of the patch (1=3x3 by default)
203 \param lambda = bandwidth ( -1 by default : automatic selection)
204 \param alpha = size of the region where similar patch are searched (3 x patch_size = 9x9 by default)
205 \param sigma = noise standard deviation (-1 = estimation)
206 \param sampling = sampling of the patch (1 = uses all point, 2 = uses one point on 4, etc)
207 If the image has three dimensions then the patch is only in 2D and the neighborhood extent in time is only 5.
208 If the image has several channel (color images), the distance between the two patch is computed using
209 all the channels.
210 The greater the patch is the best is the result.
211 Lambda parameter is function of the size of the patch size. The automatic Lambda parameter is taken
212 in the Chi2 table at a significiance level of 0.01. This diffear from the original paper [1]. The weighted average becomes then:
213 \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
214 where \f$ P(x,y) $\f denotes the patch in (x,y) location.
216 An a priori is also used to increase the speed of the algorithm in the spirit of Sapiro et al. SPletter dec 05
218 This very basic version of the Non-Local Means algorithm provides an output image which contains
219 some residual noise with a relatively small variance (\f$\sigma<5$\f).
221 [1] A non-local algorithm for image denoising
222 Buades, A.; Coll, B.; Morel, J.-M.;
223 Computer Vision and Pattern Recognition, 2005. CVPR 2005. IEEE Computer Society Conference on
224 Volume 2, 20-25 June 2005 Page(s):60 - 65 vol. 2
225 **/
226 CImg<T> get_nlmeans( int patch_size=1, double lambda=-1, double alpha=3 ,double sigma=-1, int sampling=1) {
227 return CImg<T>(*this).nlmeans(patch_size,lambda,alpha,sigma,sampling);
228 }
230 #endif