Mon, 03 Aug 2009 14:21:23 +0100
added keepme for ptdecode obj directory
1 /*
2 #
3 # File : wavelet_atrous.cpp
4 # ( C++ source file )
5 #
6 # Description : Performs a 2D or 3D 'a trous' wavelet transform
7 # (using a cubic spline) on an image or a video sequence.
8 # This file is a part of the CImg Library project.
9 # ( http://cimg.sourceforge.net )
10 #
11 # Author : Renaud Peteri
12 # ( Renaud.Peteri(at)mines-paris.org )
13 #
14 # Institution : CWI, Amsterdam
15 #
16 # Date : February 2005
17 #
18 # References : Starck, J.-L., Murtagh, F. and Bijaoui, A.,
19 # Image Processing and Data Analysis: The Multiscale Approach,
20 # Cambridge University Press, 1998.
21 # (Hardback and softback, ISBN 0-521-59084-1 and 0-521-59914-8.)
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 #include "CImg.h"
55 using namespace cimg_library;
57 // The lines below are necessary when using a non-standard compiler as visualcpp6.
58 #ifdef cimg_use_visualcpp6
59 #define std
60 #endif
61 #ifdef min
62 #undef min
63 #undef max
64 #endif
66 #ifndef cimg_imagepath
67 #define cimg_imagepath "img/"
68 #endif
70 CImg<float> mask_x(const unsigned char scale) {
71 unsigned char d1 = (unsigned char)std::pow(2.0,(double)(scale-1));
72 unsigned char d2 = (unsigned char)std::pow(2.0,(double)(scale));
73 unsigned char cx = (unsigned char)std::pow(2.0,(double)(scale));
74 unsigned char res = (unsigned char)std::pow(2.0,(double)scale);
75 CImg<float> m(2*res +1,1,1);m.fill(0);
76 m(cx) = 6.0;
77 m(cx-d1) = m(cx+d1) =4.0;
78 m(cx-d2) = m(cx+d2) =1.0;
79 m /= 16.0;
80 return m;
81 }
83 CImg<float> mask_y(const unsigned char scale) {
84 unsigned char d1 = (unsigned char)std::pow(2.0,(double)(scale-1));
85 unsigned char d2 = (unsigned char)std::pow(2.0,(double)(scale));
86 unsigned char cy = (unsigned char)std::pow(2.0,(double)(scale));
87 unsigned char res = (unsigned char)std::pow(2.0,(double)scale);
88 CImg<float> m(1,2*res +1);m.fill(0);
89 m(0,cy) = 6.0;
90 m(0,cy-d1) = m(0,cy+d1) =4.0;
91 m(0,cy-d2) = m(0,cy+d2) =1.0;
92 m /= 16.0;
93 return m;
94 }
96 CImg<float> mask_t(const unsigned char scale) {
97 unsigned char d1 = (unsigned char)std::pow(2.0,(double)(scale-1));
98 unsigned char d2 = (unsigned char)std::pow(2.0,(double)(scale));
99 unsigned char ct = (unsigned char)std::pow(2.0,(double)(scale));
100 unsigned char res = (unsigned char)std::pow(2.0,(double)scale);
101 CImg<float> m(1,1,2*res +1);m.fill(0);
102 m(0,0,ct) = 6.0;
103 m(0,0,ct-d1) = m(0,0,ct+d1) =4.0;
104 m(0,0,ct-d2) = m(0,0,ct+d2) =1.0;
105 m /= 16.0;
106 return m;
107 }
109 /*------------------
110 Main procedure
111 ----------------*/
112 int main(int argc,char **argv) {
114 cimg_usage("Perform an 'a trous' wavelet transform (using a cubic spline) on an image or on a video sequence.\n"
115 "This wavelet transform is undecimated and produces 2 images/videos at each scale. For an example of\n"
116 "decomposition on a video, try -i img/trees.inr (sequence from the MIT).\n"
117 "\t(Type -h for help)");
119 // Read command line parameters
120 const char
121 *name_i = cimg_option("-i",cimg_imagepath "lena.pgm","Input image or video"),
122 *name_o = cimg_option("-o","","Name of the multiscale analysis output"),
123 *axe_dec = cimg_option("-axe",(char*)NULL,"Perform the multiscale decomposition in just one direction ('x', 'y' or 't')");
124 const unsigned int
125 s = cimg_option("-s",3,"Scale of decomposition");
127 const bool help = cimg_option("-h",false,"Display Help");
128 if(help) exit(0);
130 // Initialize Image Data
131 std::fprintf(stderr," - Load image sequence '%s'...\n",cimg::basename(name_i));
132 const CImg<float> texture_in(name_i);
133 CImg<float> mask_conv_x, mask_conv_y, mask_conv_t;
134 CImgList<float> res(s, texture_in.dimx(),texture_in.dimy(),texture_in.dimz());
135 CImgList<float> wav(s,texture_in.dimx(), texture_in.dimy(), texture_in.dimz());
136 cimglist_for(res,l) { res(l).fill(0.0); wav(l).fill(0.0);}
137 unsigned int i;
139 if (!axe_dec){
140 // Perform the multiscale decomposition in all directions
141 for(i=0;i<s;i++){
142 std::fprintf(stderr," - Performing scale %u ...\n",i+1);
143 if(i==0){ res(i) = texture_in;} else { res(i) = res(i-1);}
144 mask_conv_x = mask_x(i+1);
145 res(i) = res(i).get_convolve(mask_conv_x);
146 mask_conv_y = mask_y(i+1);
147 res(i) = res(i).get_convolve(mask_conv_y);
148 mask_conv_t = mask_t(i+1);
149 res(i) = res(i).get_convolve(mask_conv_t);
150 if(i==0){wav(i) = texture_in - res(i);} // res(0) and wav(0) are the 1st scale of decompostion
151 else {wav(i) = res(i-1) - res(i);}
152 } }
154 if (axe_dec) {
155 // Perform the multiscale decomposition in just one direction
156 char c;
157 c = cimg::uncase(axe_dec[0]);
158 fprintf(stderr," - Decompose the image along axe '%c'\n",c); fflush(stdout);
160 switch(c) {
161 case 'x': {
162 for(i=0;i<s;i++) {
163 std::fprintf(stderr," - Performing scale %u ...\n",i+1);
164 if(i==0){ res(i) = texture_in;} else { res(i) = res(i-1);}
165 mask_conv_x = mask_x(i+1);
166 res(i) = res(i).get_convolve(mask_conv_x);
167 if(i==0){wav(i) = texture_in - res(i);}
168 else {wav(i) = res(i-1) - res(i);}}}
169 break;
171 case 'y': {
172 for(i=0;i<s;i++) {
173 std::fprintf(stderr," - Performing scale %u ...\n",i+1);
174 if(i==0){ res(i) = texture_in;} else { res(i) = res(i-1);}
175 mask_conv_y = mask_y(i+1);
176 res(i) = res(i).get_convolve(mask_conv_y);
177 if(i==0){wav(i) = texture_in - res(i);}
178 else {wav(i) = res(i-1) - res(i);}}}
179 break;
181 case 't': {
182 for(i=0;i<s;i++) {
183 std::fprintf(stderr," - Performing scale %u ...\n",i+1);
184 if(i==0){ res(i) = texture_in;} else { res(i) = res(i-1);}
185 mask_conv_t = mask_t(i+1);
186 res(i) = res(i).get_convolve(mask_conv_t);
187 if(i==0){wav(i) = texture_in - res(i);}
188 else {wav(i) = res(i-1) - res(i);}}}
189 break;
191 default: throw CImgException("Error, unknow decompostion axe '%c', try 'x', 'y' or 't'",c);
192 }
193 fputc('\n',stderr);
194 }
196 if (*name_o){
197 // Save the Multi-Scale Analysis
198 std::fprintf(stderr," - Saving of all output sequences : %s in the msa/ directory... \n",cimg::basename(name_o));
199 int count = 1; // res0 = original image
200 char filename[256] = "", filename_wav[256] = "";
201 char STmp[3] = "";
202 system("mkdir msa");
203 for(i=0;i<s;i++){
204 strcpy( filename, "msa/res" );
205 strcpy( filename_wav, "msa/wav" );
206 if( count < 10 )
207 { strcat( filename, "0" );strcat( filename_wav, "0" );}
208 sprintf( STmp, "%d_", count );
209 strcat( filename, STmp ); strcat( filename_wav, STmp );
210 strcat( filename,name_o);strcat( filename_wav,name_o);
211 res(i).save(filename);
212 wav(i).save(filename_wav);
213 count++;
214 }
215 }
217 // Result visualization
218 const float value = 255;
219 for(i=0;i<s;i++) {
220 res[i].normalize(0,255).draw_text(2,2,"Scale %d",&value,0,1,11,i);
221 wav[i].normalize(0,255).draw_text(2,2,"Scale %d",&value,0,1,11,i);
222 }
224 CImgDisplay disp(res,"Approximations levels by increasing scale",0);
225 CImgDisplay disp2(wav,"Wavelet coefficients by increasing scale",0);
226 while ( !disp.is_closed && !disp.is_keyQ && !disp.is_keyESC &&
227 !disp2.is_closed && !disp2.is_keyQ && !disp2.is_keyESC ) {
228 if (disp.is_resized) disp.resize().display(res);
229 if (disp2.is_resized) disp2.resize().display(wav);
230 CImgDisplay::wait(disp,disp2);
231 }
233 return 0;
234 }