1.1 diff -r 5edfbd3e7a46 -r 1204ebf9340d PTdecode/CImg-1.3.0/examples/pde_TschumperleDeriche2d.cpp 1.2 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 1.3 +++ b/PTdecode/CImg-1.3.0/examples/pde_TschumperleDeriche2d.cpp Mon Aug 03 14:09:20 2009 +0100 1.4 @@ -0,0 +1,235 @@ 1.5 +/* 1.6 + # 1.7 + # File : pde_TschumperleDeriche2D.cpp 1.8 + # ( C++ source file ) 1.9 + # 1.10 + # Description : Implementation of the Tschumperle-Deriche's Regularization 1.11 + # PDE, for 2D multivalued images, as described in the articles below. 1.12 + # This file is a part of the CImg Library project. 1.13 + # ( http://cimg.sourceforge.net ) 1.14 + # 1.15 + # (1) PDE-Based Regularization of Multivalued Images and Applications. 1.16 + # (D. Tschumperle). PhD Thesis. University of Nice-Sophia Antipolis, December 2002. 1.17 + # (2) Diffusion PDE's on Vector-valued Images : Local Approach and Geometric Viewpoint. 1.18 + # (D. Tschumperle and R. Deriche). IEEE Signal Processing Magazine, October 2002. 1.19 + # (3) Vector-Valued Image Regularization with PDE's : A Common Framework for Different Applications. 1.20 + # (D. Tschumperle and R. Deriche). CVPR'2003, Computer Vision and Pattern Recognition, Madison, United States, June 2003. 1.21 + # 1.22 + # This code can be used to perform image restoration, inpainting, magnification or flow visualization. 1.23 + # 1.24 + # NOTE : THIS SOURCE IS DISTRIBUTED FOR EDUCATIONAL PURPOSES ONLY. A BETTER ANISOTROPIC SMOOTHING ALGORITHM CAN BE FOUND 1.25 + # IN THE FILE 'greycstoration.cpp' WHICH IS THE RESULT OF MORE RECENT WORK. 1.26 + # 1.27 + # Copyright : David Tschumperle 1.28 + # ( http://www.greyc.ensicaen.fr/~dtschump/ ) 1.29 + # 1.30 + # License : CeCILL v2.0 1.31 + # ( http://www.cecill.info/licences/Licence_CeCILL_V2-en.html ) 1.32 + # 1.33 + # This software is governed by the CeCILL license under French law and 1.34 + # abiding by the rules of distribution of free software. You can use, 1.35 + # modify and/ or redistribute the software under the terms of the CeCILL 1.36 + # license as circulated by CEA, CNRS and INRIA at the following URL 1.37 + # "http://www.cecill.info". 1.38 + # 1.39 + # As a counterpart to the access to the source code and rights to copy, 1.40 + # modify and redistribute granted by the license, users are provided only 1.41 + # with a limited warranty and the software's author, the holder of the 1.42 + # economic rights, and the successive licensors have only limited 1.43 + # liability. 1.44 + # 1.45 + # In this respect, the user's attention is drawn to the risks associated 1.46 + # with loading, using, modifying and/or developing or reproducing the 1.47 + # software by the user in light of its specific status of free software, 1.48 + # that may mean that it is complicated to manipulate, and that also 1.49 + # therefore means that it is reserved for developers and experienced 1.50 + # professionals having in-depth computer knowledge. Users are therefore 1.51 + # encouraged to load and test the software's suitability as regards their 1.52 + # requirements in conditions enabling the security of their systems and/or 1.53 + # data to be ensured and, more generally, to use and operate it in the 1.54 + # same conditions as regards security. 1.55 + # 1.56 + # The fact that you are presently reading this means that you have had 1.57 + # knowledge of the CeCILL license and that you accept its terms. 1.58 + # 1.59 +*/ 1.60 + 1.61 +#include "CImg.h" 1.62 +using namespace cimg_library; 1.63 + 1.64 +// The lines below are necessary when using a non-standard compiler as visualcpp6. 1.65 +#ifdef cimg_use_visualcpp6 1.66 +#define std 1.67 +#endif 1.68 +#ifdef min 1.69 +#undef min 1.70 +#undef max 1.71 +#endif 1.72 + 1.73 +#ifndef cimg_imagepath 1.74 +#define cimg_imagepath "img/" 1.75 +#endif 1.76 + 1.77 +int main(int argc,char **argv) { 1.78 + 1.79 + // Read command line arguments 1.80 + //----------------------------- 1.81 + cimg_usage("Tschumperle-Deriche's flow for 2D Image Restoration, Inpainting, Magnification or Flow visualization"); 1.82 + const char *file_i = cimg_option("-i",cimg_imagepath "milla.bmp","Input image"); 1.83 + const char *file_m = cimg_option("-m",(char*)NULL,"Mask image (if Inpainting)"); 1.84 + const char *file_f = cimg_option("-f",(char*)NULL,"Flow image (if Flow visualization)"); 1.85 + const char *file_o = cimg_option("-o",(char*)NULL,"Output file"); 1.86 + const double zoom = cimg_option("-zoom",1.0,"Image magnification"); 1.87 + 1.88 + const unsigned int nb_iter = cimg_option("-iter",100000,"Number of iterations"); 1.89 + const double dt = cimg_option("-dt",20.0,"Adapting time step"); 1.90 + const double alpha = cimg_option("-alpha",0.0,"Gradient smoothing"); 1.91 + const double sigma = cimg_option("-sigma",0.5,"Structure tensor smoothing"); 1.92 + const float a1 = cimg_option("-a1",0.5f,"Diffusion limiter along minimal variations"); 1.93 + const float a2 = cimg_option("-a2",0.9f,"Diffusion limiter along maximal variations"); 1.94 + const double noiseg = cimg_option("-ng",0.0,"Add gauss noise before aplying the algorithm"); 1.95 + const double noiseu = cimg_option("-nu",0.0,"Add uniform noise before applying the algorithm"); 1.96 + const double noises = cimg_option("-ns",0.0,"Add salt&pepper noise before applying the algorithm"); 1.97 + const bool stflag = cimg_option("-stats",false,"Display image statistics at each iteration"); 1.98 + const unsigned int save = cimg_option("-save",0,"Iteration saving step"); 1.99 + const unsigned int visu = cimg_option("-visu",10,"Visualization step (0=no visualization)"); 1.100 + const unsigned int init = cimg_option("-init",3,"Inpainting initialization (0=black, 1=white, 2=noise, 3=unchanged)"); 1.101 + const unsigned int skip = cimg_option("-skip",1,"Step of image geometry computation"); 1.102 + bool view_t = cimg_option("-d",false,"View tensor directions (useful for debug)"); 1.103 + double xdt = 0; 1.104 + 1.105 + // Variable initialization 1.106 + //------------------------- 1.107 + CImg<> img, flow; 1.108 + CImg<int> mask; 1.109 + 1.110 + if (file_i) { 1.111 + img = CImg<>(file_i).resize(-100,-100,1,-100); 1.112 + if (file_m) mask = CImg<unsigned char>(file_m).resize(img.dimx(),img.dimy(),1,1); 1.113 + else if (zoom>1) { 1.114 + mask = CImg<int>(img.dimx(),img.dimy(),1,1,-1).resize((int)(img.dimx()*zoom),(int)(img.dimy()*zoom),1,1,4)+1; 1.115 + img.resize((int)(img.dimx()*zoom),(int)(img.dimy()*zoom),1,-100,3); 1.116 + } 1.117 + } else { 1.118 + if (file_f) { 1.119 + flow = CImg<>(file_f); 1.120 + img = CImg<>((int)(flow.dimx()*zoom),(int)(flow.dimy()*zoom),1,1,0).noise(100,2); 1.121 + flow.resize(img.dimx(),img.dimy(),1,2,3); 1.122 + } else throw CImgException("You need to specify at least one input image (option -i), or one flow image (option -f)"); 1.123 + } 1.124 + img.noise(noiseg,0).noise(noiseu,1).noise(noises,2); 1.125 + float initial_min, initial_max = img.maxmin(initial_min); 1.126 + if (mask.data && init!=3) 1.127 + cimg_forXYV(img,x,y,k) if (mask(x,y)) 1.128 + img(x,y,k)=(float)((init? 1.129 + (init==1?initial_max:((initial_max-initial_min)*cimg::rand())): 1.130 + initial_min)); 1.131 + 1.132 + CImgDisplay disp; 1.133 + if (visu) disp.assign(img,"Iterated Image"); 1.134 + CImg<> G(img.dimx(),img.dimy(),1,3,0), T(G), veloc(img), val(2), vec(2,2); 1.135 + 1.136 + // PDE main iteration loop 1.137 + //------------------------- 1.138 + for (unsigned int iter=0; iter<nb_iter && (!disp || (!disp.is_closed && !disp.is_keyQ && !disp.is_keyESC)); iter++) { 1.139 + std::printf("\riter %u , xdt = %g ",iter,xdt); std::fflush(stdout); 1.140 + if (stflag) img.print(); 1.141 + if (disp && disp.key==cimg::keySPACE) { view_t = !view_t; disp.key=0; } 1.142 + 1.143 + if (!(iter%skip)) { 1.144 + // Compute the tensor field T, used to drive the diffusion 1.145 + //--------------------------------------------------------- 1.146 + 1.147 + // When using PDE for flow visualization 1.148 + if (flow.data) cimg_forXY(flow,x,y) { 1.149 + const float 1.150 + u = flow(x,y,0,0), 1.151 + v = flow(x,y,0,1), 1.152 + n = (float)std::sqrt((double)(u*u+v*v)), 1.153 + nn = (n!=0)?n:1; 1.154 + T(x,y,0) = u*u/nn; 1.155 + T(x,y,1) = u*v/nn; 1.156 + T(x,y,2) = v*v/nn; 1.157 + } else { 1.158 + 1.159 + // Compute structure tensor field G 1.160 + CImgList<> grad = img.get_gradient(); 1.161 + if (alpha!=0) cimglist_for(grad,l) grad[l].blur((float)alpha); 1.162 + G.fill(0); 1.163 + cimg_forXYV(img,x,y,k) { 1.164 + const float ix = grad[0](x,y,k), iy = grad[1](x,y,k); 1.165 + G(x,y,0) += ix*ix; 1.166 + G(x,y,1) += ix*iy; 1.167 + G(x,y,2) += iy*iy; 1.168 + } 1.169 + if (sigma!=0) G.blur((float)sigma); 1.170 + 1.171 + // When using PDE for image restoration, inpainting or zooming 1.172 + T.fill(0); 1.173 + if (!mask.data) cimg_forXY(G,x,y) { 1.174 + G.get_tensor_at(x,y).symmetric_eigen(val,vec); 1.175 + const float 1.176 + l1 = (float)std::pow(1.0f+val[0]+val[1],-a1), 1.177 + l2 = (float)std::pow(1.0f+val[0]+val[1],-a2), 1.178 + ux = vec(1,0), 1.179 + uy = vec(1,1); 1.180 + T(x,y,0) = l1*ux*ux + l2*uy*uy; 1.181 + T(x,y,1) = l1*ux*uy - l2*ux*uy; 1.182 + T(x,y,2) = l1*uy*uy + l2*ux*ux; 1.183 + } 1.184 + else cimg_forXY(G,x,y) if (mask(x,y)) { 1.185 + G.get_tensor_at(x,y).symmetric_eigen(val,vec); 1.186 + const float 1.187 + ux = vec(1,0), 1.188 + uy = vec(1,1); 1.189 + T(x,y,0) = ux*ux; 1.190 + T(x,y,1) = ux*uy; 1.191 + T(x,y,2) = uy*uy; 1.192 + } 1.193 + } 1.194 + } 1.195 + 1.196 + // Compute the PDE velocity and update the iterated image 1.197 + //-------------------------------------------------------- 1.198 + CImg_3x3(I,float); 1.199 + veloc.fill(0); 1.200 + cimg_forV(img,k) cimg_for3x3(img,x,y,0,k,I) { 1.201 + const float 1.202 + a = T(x,y,0), 1.203 + b = T(x,y,1), 1.204 + c = T(x,y,2), 1.205 + ixx = Inc+Ipc-2*Icc, 1.206 + iyy = Icn+Icp-2*Icc, 1.207 + ixy = 0.25f*(Ipp+Inn-Ipn-Inp); 1.208 + veloc(x,y,k) = a*ixx + 2*b*ixy + c*iyy; 1.209 + } 1.210 + if (dt>0) { 1.211 + float m, M = veloc.maxmin(m); 1.212 + xdt = dt/cimg::max(cimg::abs(m),cimg::abs(M)); 1.213 + } else xdt=-dt; 1.214 + img+=veloc*xdt; 1.215 + img.cut((float)initial_min,(float)initial_max); 1.216 + 1.217 + // Display and save iterations 1.218 + if (disp && !(iter%visu)) { 1.219 + if (!view_t) img.display(disp); 1.220 + else { 1.221 + const unsigned char white[3] = {255,255,255}; 1.222 + CImg<unsigned char> visu = img.get_resize(disp.dimx(),disp.dimy()).normalize(0,255); 1.223 + CImg<> isophotes(img.dimx(),img.dimy(),1,2,0); 1.224 + cimg_forXY(img,x,y) if (!mask.data || mask(x,y)) { 1.225 + T.get_tensor_at(x,y).symmetric_eigen(val,vec); 1.226 + isophotes(x,y,0) = vec(0,0); 1.227 + isophotes(x,y,1) = vec(0,1); 1.228 + } 1.229 + visu.draw_quiver(isophotes,white,0.5f,10,9,0).display(disp); 1.230 + } 1.231 + } 1.232 + if (save && file_o && !(iter%save)) img.save(file_o,iter); 1.233 + if (disp) disp.resize().display(img); 1.234 + } 1.235 + 1.236 + // Save result and exit. 1.237 + if (file_o) img.save(file_o); 1.238 + return 0; 1.239 +}