1 /*
2 #
3 # File : chlpca.cpp
4 # ( C++ source file )
5 #
6 # Description : Example of use for the CImg plugin 'plugins/chlpca.h'.
7 # This file is a part of the CImg Library project.
8 # ( http://cimg.eu )
9 #
10 # Copyright : Jerome Boulanger
11 # ( http://www.irisa.fr/vista/Equipe/People/Jerome.Boulanger.html )
12 #
13 #
14 # License : CeCILL v2.0
15 # ( http://www.cecill.info/licences/Licence_CeCILL_V2-en.html )
16 #
17 # This software is governed by the CeCILL license under French law and
18 # abiding by the rules of distribution of free software. You can use,
19 # modify and/ or redistribute the software under the terms of the CeCILL
20 # license as circulated by CEA, CNRS and INRIA at the following URL
21 # "http://www.cecill.info".
22 #
23 # As a counterpart to the access to the source code and rights to copy,
24 # modify and redistribute granted by the license, users are provided only
25 # with a limited warranty and the software's author, the holder of the
26 # economic rights, and the successive licensors have only limited
27 # liability.
28 #
29 # In this respect, the user's attention is drawn to the risks associated
30 # with loading, using, modifying and/or developing or reproducing the
31 # software by the user in light of its specific status of free software,
32 # that may mean that it is complicated to manipulate, and that also
33 # therefore means that it is reserved for developers and experienced
34 # professionals having in-depth computer knowledge. Users are therefore
35 # encouraged to load and test the software's suitability as regards their
36 # requirements in conditions enabling the security of their systems and/or
37 # data to be ensured and, more generally, to use and operate it in the
38 # same conditions as regards security.
39 #
40 # The fact that you are presently reading this means that you have had
41 # knowledge of the CeCILL license and that you accept its terms.
42 #
43 */
44
45 #ifndef cimg_plugin_chlpca
46 #define cimg_plugin_chlpca
47
48 // Define some useful macros.
49
50 //! Some loops
51 #define cimg_for_step1(bound,i,step) for (int i = 0; i<(int)(bound); i+=step)
52 #define cimg_for_stepX(img,x,step) cimg_for_step1((img)._width,x,step)
53 #define cimg_for_stepY(img,y,step) cimg_for_step1((img)._height,y,step)
54 #define cimg_for_stepZ(img,z,step) cimg_for_step1((img)._depth,z,step)
55 #define cimg_for_stepXY(img,x,y,step) cimg_for_stepY(img,y,step) cimg_for_stepX(img,x,step)
56 #define cimg_for_stepXYZ(img,x,y,step) cimg_for_stepZ(img,z,step) cimg_for_stepY(img,y,step) cimg_for_stepX(img,x,step)
57
58 //! Loop for point J(xj,yj) in the neighborhood of a point I(xi,yi) of size (2*rx+1,2*ry+1)
59 /**
60 Point J is kept inside the boundaries of the image img.
61 example of summing the pixels values in a neighborhood 11x11
62 cimg_forXY(img,xi,yi) cimg_for_windowXY(img,xi,yi,xj,yj,5,5) dest(yi,yi) += src(xj,yj);
63 **/
64 #define cimg_forXY_window(img,xi,yi,xj,yj,rx,ry) \
65 for (int yi0 = std::max(0,yi-ry), yi1=std::min(yi + ry,(int)img.height() - 1), yj=yi0;yj<=yi1;++yj) \
66 for (int xi0 = std::max(0,xi-rx), xi1=std::min(xi + rx,(int)img.width() - 1), xj=xi0;xj<=xi1;++xj)
67
68 #define cimg_forXYZ_window(img,xi,yi,zi,xj,yj,zj,rx,ry,rz) \
69 for (int zi0 = std::max(0,zi-rz), zi1=std::min(zi + rz,(int)img.depth() - 1) , zj=zi0;zj<=zi1;++zj) \
70 for (int yi0 = std::max(0,yi-ry), yi1=std::min(yi + ry,(int)img.height() - 1), yj=yi0;yj<=yi1;++yj) \
71 for (int xi0 = std::max(0,xi-rx), xi1=std::min(xi + rx,(int)img.width() - 1) , xj=xi0;xj<=xi1;++xj)
72
73 //! Crop a patch in the image around position x,y,z and return a column vector
74 /**
75 \param x x-coordinate of the center of the patch
76 \param y y-coordinate of the center of the patch
77 \param z z-coordinate of the center of the patch
78 \param px the patch half width
79 \param px the patch half height
80 \param px the patch half depth
81 \return img.get_crop(x0,y0,z0,x1,y1,z1).unroll('y');
82 **/
get_patch(int x,int y,int z,int px,int py,int pz)83 CImg<T> get_patch(int x, int y, int z,
84 int px, int py, int pz) const {
85 if (depth() == 1){
86 const int x0 = x - px, y0 = y - py, x1 = x + px, y1 = y + py;
87 return get_crop(x0, y0, x1, y1).unroll('y');
88 } else {
89 const int
90 x0 = x - px, y0 = y - py, z0 = z - pz,
91 x1 = x + px, y1 = y + py, z1 = z + pz;
92 return get_crop(x0, y0, z0, x1, y1, z1).unroll('y');
93 }
94 }
95
96 //! Extract a local patch dictionnary around point xi,yi,zi
get_patch_dictionnary(const int xi,const int yi,const int zi,const int px,const int py,const int pz,const int wx,const int wy,const int wz,int & idc)97 CImg<T> get_patch_dictionnary(const int xi, const int yi, const int zi,
98 const int px, const int py, const int pz,
99 const int wx, const int wy, const int wz,
100 int & idc) const {
101 const int
102 n = (2*wx + 1) * (2*wy + 1) * (2 * (depth()==1?0:wz) + 1),
103 d = (2*px + 1) * (2*py + 1) * (2 * (depth()==1?0:px) + 1) * spectrum();
104 CImg<> S(n, d);
105 int idx = 0;
106 if (depth() == 1) {
107 cimg_forXY_window((*this), xi, yi, xj, yj, wx, wy){
108 CImg<T> patch = get_patch(xj, yj, 0, px, py, 1);
109 cimg_forY(S,y) S(idx,y) = patch(y);
110 if (xj==xi && yj==yi) idc = idx;
111 idx++;
112 }
113 } else {
114 cimg_forXYZ_window((*this), xi,yi,zi,xj,yj,zj,wx,wy,wz){
115 CImg<T> patch = get_patch(xj, yj, zj, px, py, pz);
116 cimg_forY(S,y) S(idx,y) = patch(y);
117 if (xj==xi && yj==yi && zj==zi) idc = idx;
118 idx++;
119 }
120 }
121 S.columns(0, idx - 1);
122 return S;
123 }
124
125 //! Add a patch to the image
126 /**
127 \param x x-coordinate of the center of the patch
128 \param y y-coordinate of the center of the patch
129 \param z z-coordinate of the center of the patch
130 \param img the patch as a 1D column vector
131 \param px the patch half width
132 \param px the patch half height
133 \param px the patch half depth
134 **/
add_patch(const int xi,const int yi,const int zi,const CImg<T> & patch,const int px,const int py,const int pz)135 CImg<T> & add_patch(const int xi, const int yi, const int zi,
136 const CImg<T> & patch,
137 const int px, const int py, const int pz) {
138 const int
139 x0 = xi - px, y0 = yi - py, z0 = (depth() == 1 ? 0 : zi - pz),
140 sx = 2 * px + 1, sy = 2 * py + 1, sz = (depth() == 1 ? 1 : 2 * pz +1);
141 draw_image(x0, y0, z0, 0, patch.get_resize(sx, sy, sz, spectrum(), -1), -1);
142 return (*this);
143 }
144
145 //! Add a constant patch to the image
146 /**
147 \param x x-coordinate of the center of the patch
148 \param y y-coordinate of the center of the patch
149 \param z z-coordinate of the center of the patch
150 \param value in the patch
151 \param px the patch half width
152 \param px the patch half height
153 \param px the patch half depth
154 **/
add_patch(const int xi,const int yi,const int zi,const T value,const int px,const int py,const int pz)155 CImg<T> & add_patch(const int xi, const int yi, const int zi, const T value,
156 const int px, const int py, const int pz) {
157 const int
158 x0 = xi - px, y0 = yi - py, z0 = (depth() == 1 ? 0 : zi - pz),
159 x1 = xi + px, y1 = yi + py, z1 = (depth() == 1 ? 0 : zi + pz);
160 draw_rectangle(x0, y0, z0, 0, x1, y1, z1, spectrum()-1, value, -1);
161 return (*this);
162 }
163
164 //! CHLPCA denoising from the PhD thesis of Hu Haijuan
165 /**
166 \param px the patch half width
167 \param py the patch half height
168 \param pz the patch half depth
169 \param wx the training region half width
170 \param wy the training region half height
171 \param wz the training region half depth
172 \param nstep the subsampling of the image domain
173 \param nsim the number of patches used for training as a factor of the patch size
174 \param lambda_min the threshold on the eigen values of the PCA for dimension reduction
175 \param threshold the threshold on the value of the coefficients
176 \param pca_use_svd if true use the svd approach to perform the pca otherwise use the covariance method
177 \note please cite the PhD thesis of Hu Haijuan http://www.univ-ubs.fr/soutenance-de-these-hu-haijuan-337653.kjsp?RH=1318498222799
178 **/
get_chlpca(const int px,const int py,const int pz,const int wx,const int wy,const int wz,const int nstep,const float nsim,const float lambda_min,const float threshold,const float noise_std,const bool pca_use_svd)179 CImg<T> get_chlpca(const int px, const int py, const int pz,
180 const int wx, const int wy, const int wz,
181 const int nstep, const float nsim,
182 const float lambda_min, const float threshold,
183 const float noise_std, const bool pca_use_svd) const {
184 const int
185 nd = (2*px + 1) * (2*py + 1) * (depth()==1?1:2*pz + 1) * spectrum(),
186 K = (int)(nsim * nd);
187 #ifdef DEBUG
188 fprintf(stderr,"chlpca: p:%dx%dx%d,w:%dx%dx%d,nd:%d,K:%d\n",
189 2*px + 1,2*py + 1,2*pz + 1,2*wx + 1,2*wy + 1,2*wz + 1,nd,K);
190 #endif
191 float sigma;
192 if (noise_std<0) sigma = (float)std::sqrt(variance_noise());
193 else sigma = noise_std;
194 CImg<T> dest(*this), count(*this);
195 dest.fill(0);
196 count.fill(0);
197 cimg_for_stepZ(*this,zi,(depth()==1||pz==0)?1:nstep){
198 #ifdef cimg_use_openmp
199
200 cimg_pragma_openmp(parallel for)
201 #endif
202 cimg_for_stepXY((*this),xi,yi,nstep){
203 // extract the training region X
204 int idc = 0;
205 CImg<T> S = get_patch_dictionnary(xi,yi,zi,px,py,pz,wx,wy,wz,idc);
206 // select the K most similar patches within the training set
207 CImg<T> Sk(S);
208 CImg<unsigned int> a_index(S.width());
209 if (K < Sk.width() - 1){
210 CImg<T> mse(S.width());
211 CImg<unsigned int> perms;
212 cimg_forX(S,x) { mse(x) = (T)S.get_column(idc).MSE(S.get_column(x)); }
213 mse.sort(perms,true);
214 cimg_foroff(perms,i) {
215 cimg_forY(S,j) Sk(i,j) = S(perms(i),j);
216 a_index(perms(i)) = i;
217 }
218 Sk.columns(0, K);
219 perms.threshold(K);
220 } else {
221 cimg_foroff(a_index,i) a_index(i)=i;
222 }
223 // centering the patches
224 CImg<T> M(1, Sk.height(), 1, 1, 0);
225 cimg_forXY(Sk,x,y) { M(y) += Sk(x,y); }
226 M /= (T)Sk.width();
227 cimg_forXY(Sk,x,y) { Sk(x,y) -= M(y); }
228 // compute the principal component of the training set S
229 CImg<T> P, lambda;
230 if (pca_use_svd) {
231 CImg<T> V;
232 Sk.get_transpose().SVD(V,lambda,P,true,100);
233 } else {
234 (Sk * Sk.get_transpose()).symmetric_eigen(lambda, P);
235 lambda.sqrt();
236 }
237 // dimension reduction
238 int s = 0;
239 const T tx = (T)(std::sqrt((double)Sk.width()-1.0) * lambda_min * sigma);
240 while((lambda(s) > tx) && (s < ((int)lambda.size() - 1))) { s++; }
241 P.columns(0,s);
242 // project all the patches on the basis (compute scalar product)
243 Sk = P.get_transpose() * Sk;
244 // threshold the coefficients
245 if (threshold > 0) { Sk.threshold(threshold, 1); }
246 // project back to pixel space
247 Sk = P * Sk;
248 // recenter the patches
249 cimg_forXY(Sk,x,y) { Sk(x,y) += M(y); }
250 int j = 0;
251 cimg_forXYZ_window((*this),xi,yi,zi,xj,yj,zj,wx,wy,wz){
252 const int id = a_index(j);
253 if (id < Sk.width()) {
254 dest.add_patch(xj, yj, zj, Sk.get_column(id), px, py, pz);
255 count.add_patch(xj, yj, zj, (T)1, px, py, pz);
256 }
257 j++;
258 }
259 }
260 }
261 cimg_foroff(dest, i) {
262 if(count(i) != 0) { dest(i) /= count(i); }
263 else { dest(i) = (*this)(i); }
264 }
265 return dest;
266 }
267
268 //! CHLPCA denoising from the PhD thesis of Hu Haijuan
269 /**
270 \param px the patch half width
271 \param px the patch half height
272 \param px the patch half depth
273 \param wx the training region half width
274 \param wy the training region half height
275 \param wz the training region half depth
276 \param nstep the subsampling of the image domain
277 \param nsim the number of patches used for training as a factor of the patch size
278 \param lambda_min the threshold on the eigen values of the PCA for dimension reduction
279 \param threshold the threshold on the value of the coefficients
280 \param pca_use_svd if true use the svd approach to perform the pca otherwise use the covariance method
281 \note please cite the PhD thesis of Hu Haijuan http://www.univ-ubs.fr/soutenance-de-these-hu-haijuan-337653.kjsp?RH=1318498222799
282 **/
chlpca(const int px,const int py,const int pz,const int wx,const int wy,const int wz,const int nstep,const float nsim,const float lambda_min,const float threshold,const float noise_std,const bool pca_use_svd)283 CImg<T> & chlpca(const int px, const int py, const int pz,
284 const int wx, const int wy, const int wz,
285 const int nstep, const float nsim,
286 const float lambda_min, const float threshold,
287 const float noise_std, const bool pca_use_svd) {
288 (*this) = get_chlpca(px, py, pz, wx, wy, wz, nstep, nsim, lambda_min,
289 threshold, noise_std, pca_use_svd);
290 return (*this);
291 }
292
293 //! CHLPCA denoising from the PhD thesis of Hu Haijuan
294 /**
295 \param p the patch half size
296 \param w the training region half size
297 \param nstep the subsampling of the image domain
298 \param nsim the number of patches used for training as a factor of the patch size
299 \param lambda_min the threshold on the eigen values of the PCA for dimension reduction
300 \param threshold the threshold on the value of the coefficients
301 \param pca_use_svd if true use the svd approach to perform the pca otherwise use the covariance method
302 \note please cite the PhD thesis of Hu Haijuan http://www.univ-ubs.fr/soutenance-de-these-hu-haijuan-337653.kjsp?RH=1318498222799
303 **/
304 CImg<T> get_chlpca(const int p=3, const int w=10,
305 const int nstep=5, const float nsim=10,
306 const float lambda_min=2, const float threshold = -1,
307 const float noise_std=-1, const bool pca_use_svd=true) const {
308 if (depth()==1) return get_chlpca(p, p, 0, w, w, 0, nstep, nsim, lambda_min,
309 threshold, noise_std, pca_use_svd);
310 else return get_chlpca(p, p, p, w, w, w, nstep, nsim, lambda_min,
311 threshold, noise_std, pca_use_svd);
312 }
313
314 CImg<T> chlpca(const int p=3, const int w=10,
315 const int nstep=5, const float nsim=10,
316 const float lambda_min=2, const float threshold = -1,
317 const float noise_std=-1, const bool pca_use_svd=true) {
318 (*this) = get_chlpca(p, w, nstep, nsim, lambda_min,
319 threshold, noise_std, pca_use_svd);
320 return (*this);
321 }
322
323 #endif /* cimg_plugin_chlpca */
324