Wed, 05 Aug 2009 15:00:54 +0100
small changes to hexdump code to stop a gcc warning on platforms where sizeof(int) != sizeof(int*) e.g. x86_64
1 /*
2 #
3 # File : image_registration.cpp
4 # ( C++ source file )
5 #
6 # Description : Compute a motion field between two images,
7 # with a multiscale and variational algorithm.
8 # This file is a part of the CImg Library project.
9 # ( http://cimg.sourceforge.net )
10 #
11 # Copyright : David Tschumperle
12 # ( http://www.greyc.ensicaen.fr/~dtschump/ )
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 */
45 #include "CImg.h"
46 using namespace cimg_library;
48 // The lines below are necessary when using a non-standard compiler as visualcpp6.
49 #ifdef cimg_use_visualcpp6
50 #define std
51 #endif
52 #ifdef min
53 #undef min
54 #undef max
55 #endif
57 #ifndef cimg_imagepath
58 #define cimg_imagepath "img/"
59 #endif
61 // animate_warp() : Create warping animation from two images and a motion field
62 //----------------
63 void animate_warp(const CImg<unsigned char>& src, const CImg<unsigned char>& dest, const CImg<>& u,
64 const bool morph, const bool imode, const char *filename,int nb, CImgDisplay& disp) {
65 CImg<unsigned char> visu = CImgList<unsigned char>(src,dest,src).get_append('x'), warp(src);
66 float t=0;
67 for (unsigned int iter=0; !disp || (!disp.is_closed && !disp.is_keyQ); iter++) {
68 if (morph) cimg_forXYV(warp,x,y,k) {
69 const float dx = u(x,y,0), dy = u(x,y,1),
70 I1 = (float)src.linear_atXY(x-t*dx, y-t*dy, k),
71 I2 = (float)dest.linear_atXY(x+(1-t)*dx,y+(1-t)*dy,k);
72 warp(x,y,k) = (unsigned char)((1-t)*I1 + t*I2);
73 } else cimg_forXYV(warp,x,y,k) {
74 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);
75 warp(x,y,k) = (unsigned char)I1;
76 }
77 if (disp) visu.draw_image(2*src.dimx(),warp).display(disp.resize().wait(30));
78 if (filename && *filename && (imode || (int)iter<nb)) {
79 std::fprintf(stderr,"\r > frame %d ",iter);
80 warp.save(filename,iter);
81 }
82 t+=1.0f/nb;
83 if (t<0) { t=0; nb=-nb; }
84 if (t>1) { t=1; nb=-nb; if (filename && *filename) std::exit(0); }
85 }
86 }
88 // get_warp() : Return the image src warped by the motion field u.
89 //------------
90 template<typename T> CImg<T> getwarp(const CImg<T>& src, const CImg<>& u) {
91 CImg<T> warp(src);
92 cimg_forXY(warp,x,y) warp(x,y) = (T)src.linear_atXY(x - u(x,y,0), y - u(x,y,1));
93 return warp;
94 }
96 // optmonoflow() : Register images for one scale ( semi-implicite PDE scheme ) between I2->I1
97 //---------------
98 CImg<> optmonoflow(const CImg<>& I1, const CImg<>& I2, const CImg<>& u0,
99 const float smooth, const float precision, CImgDisplay& disp) {
101 CImg<> u = u0.get_resize(I1.dimx(),I1.dimy(),1,2,3),dI(u);
102 CImg_3x3(I,float);
103 float dt=2,E=1e20f;
105 // compute first derivatives of I2
106 cimg_for3x3(I2,x,y,0,0,I) {
107 dI(x,y,0) = 0.5f*(Inc-Ipc);
108 dI(x,y,1) = 0.5f*(Icn-Icp);
109 }
111 // Main PDE iteration
112 for (unsigned int iter=0; iter<100000; iter++) {
113 std::fprintf(stderr,"\r- Iteration %d - E = %g",iter,E); std::fflush(stderr);
114 const float Eold = E;
115 E = 0;
116 cimg_for3XY(u,x,y) {
117 const float
118 X = x + u(x,y,0),
119 Y = y + u(x,y,1),
120 deltaI = (float)(I2.linear_atXY(X,Y) - I1(x,y));
121 float tmpf = 0;
122 cimg_forV(u,k) {
123 const float
124 ux = 0.5f*(u(_n1x,y,k)-u(_p1x,y,k)),
125 uy = 0.5f*(u(x,_n1y,k)-u(x,_p1y,k));
126 u(x,y,k) = (float)( u(x,y,k) +
127 dt*(
128 -deltaI*dI.linear_atXY(X,Y,k) +
129 smooth* ( u(_n1x,y,k) + u(_p1x,y,k) + u(x,_n1y,k) + u(x,_p1y,k) )
130 )
131 )/(1+4*smooth*dt);
132 tmpf += ux*ux + uy*uy;
133 }
134 E += deltaI*deltaI + smooth * tmpf;
135 }
136 if (cimg::abs(Eold-E)<precision) break;
137 if (Eold<E) dt*=0.5;
138 if (disp) disp.resize();
139 if (disp && disp.is_closed) std::exit(0);
140 if (disp && !(iter%300)) {
141 const unsigned char white = 255;
142 CImg<unsigned char> tmp = getwarp(I1,u).normalize(0,200);
143 tmp.resize(disp.dimx(),disp.dimy()).draw_quiver(u,&white,0.7f,15,-14,0).display(disp);
144 }
145 }
146 return u;
147 }
149 // optflow() : multiscale version of the image registration algorithm
150 //-----------
151 CImg<> optflow(const CImg<>& xsrc, const CImg<>& xdest,
152 const float smooth, const float precision, const unsigned int pnb_scale, CImgDisplay& disp) {
153 const CImg<>
154 src = xsrc.get_pointwise_norm(1).resize(xdest.dimx(),xdest.dimy(),1,1,3).normalize(0,1),
155 dest = xdest.get_pointwise_norm(1).resize(xdest.dimx(),xdest.dimy(),1,1,3).normalize(0,1);
156 CImg<> u = CImg<>(src.dimx(),src.dimy(),1,2).fill(0);
158 const unsigned int nb_scale = pnb_scale>0?pnb_scale:(unsigned int)(2*std::log((double)(cimg::max(src.dimx(),src.dimy()))));
159 for (int scale=nb_scale-1; scale>=0; scale--) {
160 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);
161 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);
162 std::fprintf(stderr," * Scale %d\n",scale);
163 u*=1.5;
164 u = optmonoflow(I1,I2,u,smooth,(float)(precision/std::pow(2.25,1+scale)),disp);
165 std::fprintf(stderr,"\n");
166 }
167 return u;
168 }
170 /*------------------------
172 Main function
174 ------------------------*/
176 int main(int argc,char **argv) {
178 // Read command line parameters
179 cimg_usage("Compute an optical flow between two 2D images, and create a warped animation");
180 const char
181 *name_i1 = cimg_option("-i",cimg_imagepath "sh0r.pgm","Input Image 1 (Destination)"),
182 *name_i2 = cimg_option("-i2",cimg_imagepath "sh1r.pgm","Input Image 2 (Source)"),
183 *name_o = cimg_option("-o",(const char*)NULL,"Output 2D flow (inrimage)"),
184 *name_seq = cimg_option("-o2",(const char*)NULL,"Output Warping Sequence");
185 const float
186 smooth = cimg_option("-s",0.1f,"Flow Smoothness"),
187 precision = cimg_option("-p",0.9f,"Convergence precision");
188 const unsigned int
189 nb = cimg_option("-n",40,"Number of warped frames"),
190 nbscale = cimg_option("-scale",0,"Number of scales (0=auto)");
191 const bool
192 normalize = cimg_option("-equalize",true,"Histogram normalization of the images"),
193 morph = cimg_option("-m",true,"Morphing mode"),
194 imode = cimg_option("-c",true,"Complete interpolation (or last frame is missing)"),
195 dispflag = !cimg_option("-novisu",false,"Visualization");
197 // Init images and display
198 std::fprintf(stderr," - Init images.\n");
199 const CImg<>
200 src(name_i1),
201 dest(CImg<>(name_i2).resize(src,3)),
202 src_blur = normalize?src.get_blur(0.5f).equalize(256):src.get_blur(0.5f),
203 dest_blur = normalize?dest.get_blur(0.5f).equalize(256):dest.get_blur(0.5f);
205 CImgDisplay disp;
206 if (dispflag) {
207 unsigned int w = src.dimx(), h = src.dimy();
208 const unsigned int dmin = cimg::min(w,h), minsiz = 512;
209 if (dmin<minsiz) { w=w*minsiz/dmin; h=h*minsiz/dmin; }
210 const unsigned int dmax = cimg::max(w,h), maxsiz = 1024;
211 if (dmax>maxsiz) { w=w*maxsiz/dmax; h=h*maxsiz/dmax; }
212 disp.assign(w,h,"Estimated Motion",0);
213 }
215 // Run Motion estimation algorithm
216 std::fprintf(stderr," - Compute optical flow.\n");
217 const CImg<> u = optflow(src_blur,dest_blur,smooth,precision,nbscale,disp);
218 if (name_o) u.save(name_o);
219 u.print("Computed flow");
221 // Do morphing animation
222 std::fprintf(stderr," - Create warped animation.\n");
223 CImgDisplay disp2;
224 if (dispflag) {
225 unsigned int w = src.dimx(), h = src.dimy();
226 const unsigned int dmin = cimg::min(w,h), minsiz = 100;
227 if (dmin<minsiz) { w=w*minsiz/dmin; h=h*minsiz/dmin; }
228 const unsigned int dmax = cimg::max(w,h), maxsiz = 1024/3;
229 if (dmax>maxsiz) { w=w*maxsiz/dmax; h=h*maxsiz/dmax; }
230 disp2.assign(3*w,h,"Source/Destination images and Motion animation",0);
231 }
233 animate_warp(src.get_normalize(0,255),dest.get_normalize(0,255),u,morph,imode,name_seq,nb,disp2);
235 std::exit(0);
236 return 0;
237 }