1.1 diff -r 5edfbd3e7a46 -r 1204ebf9340d PTdecode/CImg-1.3.0/examples/image_registration.cpp 1.2 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 1.3 +++ b/PTdecode/CImg-1.3.0/examples/image_registration.cpp Mon Aug 03 14:09:20 2009 +0100 1.4 @@ -0,0 +1,237 @@ 1.5 +/* 1.6 + # 1.7 + # File : image_registration.cpp 1.8 + # ( C++ source file ) 1.9 + # 1.10 + # Description : Compute a motion field between two images, 1.11 + # with a multiscale and variational algorithm. 1.12 + # This file is a part of the CImg Library project. 1.13 + # ( http://cimg.sourceforge.net ) 1.14 + # 1.15 + # Copyright : David Tschumperle 1.16 + # ( http://www.greyc.ensicaen.fr/~dtschump/ ) 1.17 + # 1.18 + # License : CeCILL v2.0 1.19 + # ( http://www.cecill.info/licences/Licence_CeCILL_V2-en.html ) 1.20 + # 1.21 + # This software is governed by the CeCILL license under French law and 1.22 + # abiding by the rules of distribution of free software. You can use, 1.23 + # modify and/ or redistribute the software under the terms of the CeCILL 1.24 + # license as circulated by CEA, CNRS and INRIA at the following URL 1.25 + # "http://www.cecill.info". 1.26 + # 1.27 + # As a counterpart to the access to the source code and rights to copy, 1.28 + # modify and redistribute granted by the license, users are provided only 1.29 + # with a limited warranty and the software's author, the holder of the 1.30 + # economic rights, and the successive licensors have only limited 1.31 + # liability. 1.32 + # 1.33 + # In this respect, the user's attention is drawn to the risks associated 1.34 + # with loading, using, modifying and/or developing or reproducing the 1.35 + # software by the user in light of its specific status of free software, 1.36 + # that may mean that it is complicated to manipulate, and that also 1.37 + # therefore means that it is reserved for developers and experienced 1.38 + # professionals having in-depth computer knowledge. Users are therefore 1.39 + # encouraged to load and test the software's suitability as regards their 1.40 + # requirements in conditions enabling the security of their systems and/or 1.41 + # data to be ensured and, more generally, to use and operate it in the 1.42 + # same conditions as regards security. 1.43 + # 1.44 + # The fact that you are presently reading this means that you have had 1.45 + # knowledge of the CeCILL license and that you accept its terms. 1.46 + # 1.47 +*/ 1.48 + 1.49 +#include "CImg.h" 1.50 +using namespace cimg_library; 1.51 + 1.52 +// The lines below are necessary when using a non-standard compiler as visualcpp6. 1.53 +#ifdef cimg_use_visualcpp6 1.54 +#define std 1.55 +#endif 1.56 +#ifdef min 1.57 +#undef min 1.58 +#undef max 1.59 +#endif 1.60 + 1.61 +#ifndef cimg_imagepath 1.62 +#define cimg_imagepath "img/" 1.63 +#endif 1.64 + 1.65 +// animate_warp() : Create warping animation from two images and a motion field 1.66 +//---------------- 1.67 +void animate_warp(const CImg<unsigned char>& src, const CImg<unsigned char>& dest, const CImg<>& u, 1.68 + const bool morph, const bool imode, const char *filename,int nb, CImgDisplay& disp) { 1.69 + CImg<unsigned char> visu = CImgList<unsigned char>(src,dest,src).get_append('x'), warp(src); 1.70 + float t=0; 1.71 + for (unsigned int iter=0; !disp || (!disp.is_closed && !disp.is_keyQ); iter++) { 1.72 + if (morph) cimg_forXYV(warp,x,y,k) { 1.73 + const float dx = u(x,y,0), dy = u(x,y,1), 1.74 + I1 = (float)src.linear_atXY(x-t*dx, y-t*dy, k), 1.75 + I2 = (float)dest.linear_atXY(x+(1-t)*dx,y+(1-t)*dy,k); 1.76 + warp(x,y,k) = (unsigned char)((1-t)*I1 + t*I2); 1.77 + } else cimg_forXYV(warp,x,y,k) { 1.78 + const float dx = u(x,y,0), dy = u(x,y,1), I1 = (float)src.linear_atXY(x-t*dx, y-t*dy, 0,k); 1.79 + warp(x,y,k) = (unsigned char)I1; 1.80 + } 1.81 + if (disp) visu.draw_image(2*src.dimx(),warp).display(disp.resize().wait(30)); 1.82 + if (filename && *filename && (imode || (int)iter<nb)) { 1.83 + std::fprintf(stderr,"\r > frame %d ",iter); 1.84 + warp.save(filename,iter); 1.85 + } 1.86 + t+=1.0f/nb; 1.87 + if (t<0) { t=0; nb=-nb; } 1.88 + if (t>1) { t=1; nb=-nb; if (filename && *filename) std::exit(0); } 1.89 + } 1.90 +} 1.91 + 1.92 +// get_warp() : Return the image src warped by the motion field u. 1.93 +//------------ 1.94 +template<typename T> CImg<T> getwarp(const CImg<T>& src, const CImg<>& u) { 1.95 + CImg<T> warp(src); 1.96 + cimg_forXY(warp,x,y) warp(x,y) = (T)src.linear_atXY(x - u(x,y,0), y - u(x,y,1)); 1.97 + return warp; 1.98 +} 1.99 + 1.100 +// optmonoflow() : Register images for one scale ( semi-implicite PDE scheme ) between I2->I1 1.101 +//--------------- 1.102 +CImg<> optmonoflow(const CImg<>& I1, const CImg<>& I2, const CImg<>& u0, 1.103 + const float smooth, const float precision, CImgDisplay& disp) { 1.104 + 1.105 + CImg<> u = u0.get_resize(I1.dimx(),I1.dimy(),1,2,3),dI(u); 1.106 + CImg_3x3(I,float); 1.107 + float dt=2,E=1e20f; 1.108 + 1.109 + // compute first derivatives of I2 1.110 + cimg_for3x3(I2,x,y,0,0,I) { 1.111 + dI(x,y,0) = 0.5f*(Inc-Ipc); 1.112 + dI(x,y,1) = 0.5f*(Icn-Icp); 1.113 + } 1.114 + 1.115 + // Main PDE iteration 1.116 + for (unsigned int iter=0; iter<100000; iter++) { 1.117 + std::fprintf(stderr,"\r- Iteration %d - E = %g",iter,E); std::fflush(stderr); 1.118 + const float Eold = E; 1.119 + E = 0; 1.120 + cimg_for3XY(u,x,y) { 1.121 + const float 1.122 + X = x + u(x,y,0), 1.123 + Y = y + u(x,y,1), 1.124 + deltaI = (float)(I2.linear_atXY(X,Y) - I1(x,y)); 1.125 + float tmpf = 0; 1.126 + cimg_forV(u,k) { 1.127 + const float 1.128 + ux = 0.5f*(u(_n1x,y,k)-u(_p1x,y,k)), 1.129 + uy = 0.5f*(u(x,_n1y,k)-u(x,_p1y,k)); 1.130 + u(x,y,k) = (float)( u(x,y,k) + 1.131 + dt*( 1.132 + -deltaI*dI.linear_atXY(X,Y,k) + 1.133 + smooth* ( u(_n1x,y,k) + u(_p1x,y,k) + u(x,_n1y,k) + u(x,_p1y,k) ) 1.134 + ) 1.135 + )/(1+4*smooth*dt); 1.136 + tmpf += ux*ux + uy*uy; 1.137 + } 1.138 + E += deltaI*deltaI + smooth * tmpf; 1.139 + } 1.140 + if (cimg::abs(Eold-E)<precision) break; 1.141 + if (Eold<E) dt*=0.5; 1.142 + if (disp) disp.resize(); 1.143 + if (disp && disp.is_closed) std::exit(0); 1.144 + if (disp && !(iter%300)) { 1.145 + const unsigned char white = 255; 1.146 + CImg<unsigned char> tmp = getwarp(I1,u).normalize(0,200); 1.147 + tmp.resize(disp.dimx(),disp.dimy()).draw_quiver(u,&white,0.7f,15,-14,0).display(disp); 1.148 + } 1.149 + } 1.150 + return u; 1.151 +} 1.152 + 1.153 +// optflow() : multiscale version of the image registration algorithm 1.154 +//----------- 1.155 +CImg<> optflow(const CImg<>& xsrc, const CImg<>& xdest, 1.156 + const float smooth, const float precision, const unsigned int pnb_scale, CImgDisplay& disp) { 1.157 + const CImg<> 1.158 + src = xsrc.get_pointwise_norm(1).resize(xdest.dimx(),xdest.dimy(),1,1,3).normalize(0,1), 1.159 + dest = xdest.get_pointwise_norm(1).resize(xdest.dimx(),xdest.dimy(),1,1,3).normalize(0,1); 1.160 + CImg<> u = CImg<>(src.dimx(),src.dimy(),1,2).fill(0); 1.161 + 1.162 + const unsigned int nb_scale = pnb_scale>0?pnb_scale:(unsigned int)(2*std::log((double)(cimg::max(src.dimx(),src.dimy())))); 1.163 + for (int scale=nb_scale-1; scale>=0; scale--) { 1.164 + const CImg<> I1 = src.get_resize((int)(src.dimx()/std::pow(1.5,scale)), (int)(src.dimy()/std::pow(1.5,scale)) ,1,1,3); 1.165 + const CImg<> I2 = dest.get_resize((int)(src.dimx()/std::pow(1.5,scale)), (int)(src.dimy()/std::pow(1.5,scale)) ,1,1,3); 1.166 + std::fprintf(stderr," * Scale %d\n",scale); 1.167 + u*=1.5; 1.168 + u = optmonoflow(I1,I2,u,smooth,(float)(precision/std::pow(2.25,1+scale)),disp); 1.169 + std::fprintf(stderr,"\n"); 1.170 + } 1.171 + return u; 1.172 +} 1.173 + 1.174 +/*------------------------ 1.175 + 1.176 + Main function 1.177 + 1.178 + ------------------------*/ 1.179 + 1.180 +int main(int argc,char **argv) { 1.181 + 1.182 + // Read command line parameters 1.183 + cimg_usage("Compute an optical flow between two 2D images, and create a warped animation"); 1.184 + const char 1.185 + *name_i1 = cimg_option("-i",cimg_imagepath "sh0r.pgm","Input Image 1 (Destination)"), 1.186 + *name_i2 = cimg_option("-i2",cimg_imagepath "sh1r.pgm","Input Image 2 (Source)"), 1.187 + *name_o = cimg_option("-o",(const char*)NULL,"Output 2D flow (inrimage)"), 1.188 + *name_seq = cimg_option("-o2",(const char*)NULL,"Output Warping Sequence"); 1.189 + const float 1.190 + smooth = cimg_option("-s",0.1f,"Flow Smoothness"), 1.191 + precision = cimg_option("-p",0.9f,"Convergence precision"); 1.192 + const unsigned int 1.193 + nb = cimg_option("-n",40,"Number of warped frames"), 1.194 + nbscale = cimg_option("-scale",0,"Number of scales (0=auto)"); 1.195 + const bool 1.196 + normalize = cimg_option("-equalize",true,"Histogram normalization of the images"), 1.197 + morph = cimg_option("-m",true,"Morphing mode"), 1.198 + imode = cimg_option("-c",true,"Complete interpolation (or last frame is missing)"), 1.199 + dispflag = !cimg_option("-novisu",false,"Visualization"); 1.200 + 1.201 + // Init images and display 1.202 + std::fprintf(stderr," - Init images.\n"); 1.203 + const CImg<> 1.204 + src(name_i1), 1.205 + dest(CImg<>(name_i2).resize(src,3)), 1.206 + src_blur = normalize?src.get_blur(0.5f).equalize(256):src.get_blur(0.5f), 1.207 + dest_blur = normalize?dest.get_blur(0.5f).equalize(256):dest.get_blur(0.5f); 1.208 + 1.209 + CImgDisplay disp; 1.210 + if (dispflag) { 1.211 + unsigned int w = src.dimx(), h = src.dimy(); 1.212 + const unsigned int dmin = cimg::min(w,h), minsiz = 512; 1.213 + if (dmin<minsiz) { w=w*minsiz/dmin; h=h*minsiz/dmin; } 1.214 + const unsigned int dmax = cimg::max(w,h), maxsiz = 1024; 1.215 + if (dmax>maxsiz) { w=w*maxsiz/dmax; h=h*maxsiz/dmax; } 1.216 + disp.assign(w,h,"Estimated Motion",0); 1.217 + } 1.218 + 1.219 + // Run Motion estimation algorithm 1.220 + std::fprintf(stderr," - Compute optical flow.\n"); 1.221 + const CImg<> u = optflow(src_blur,dest_blur,smooth,precision,nbscale,disp); 1.222 + if (name_o) u.save(name_o); 1.223 + u.print("Computed flow"); 1.224 + 1.225 + // Do morphing animation 1.226 + std::fprintf(stderr," - Create warped animation.\n"); 1.227 + CImgDisplay disp2; 1.228 + if (dispflag) { 1.229 + unsigned int w = src.dimx(), h = src.dimy(); 1.230 + const unsigned int dmin = cimg::min(w,h), minsiz = 100; 1.231 + if (dmin<minsiz) { w=w*minsiz/dmin; h=h*minsiz/dmin; } 1.232 + const unsigned int dmax = cimg::max(w,h), maxsiz = 1024/3; 1.233 + if (dmax>maxsiz) { w=w*maxsiz/dmax; h=h*maxsiz/dmax; } 1.234 + disp2.assign(3*w,h,"Source/Destination images and Motion animation",0); 1.235 + } 1.236 + 1.237 + animate_warp(src.get_normalize(0,255),dest.get_normalize(0,255),u,morph,imode,name_seq,nb,disp2); 1.238 + 1.239 + std::exit(0); 1.240 + return 0; 1.241 +}