00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #include <gvars3/instances.h>
00021 #include <cvd/image_io.h>
00022 #include <cvd/convolution.h>
00023 #include <TooN/wls.h>
00024 #include <TooN/SVD.h>
00025 #include <tr1/tuple>
00026
00027 #include "storm_imagery.h"
00028 #include "debug.h"
00029 #include "utility.h"
00030
00031 using namespace CVD;
00032 using namespace TooN;
00033 using namespace GVars3;
00034 using namespace std;
00035 using namespace std::tr1;
00036
00037
00038
00039
00040
00041
00042
00043 vector<Image<float> > load_and_preprocess_images2(const vector<string>& names)
00044 {
00045 vector<Image<float> > ims;
00046
00047 for(unsigned int i=0; i < names.size(); i++)
00048 {
00049 Image<float> im = img_load(names[i]);
00050 ims.push_back(im);
00051
00052 if(ims.back().size() != ims[0].size())
00053 {
00054 cerr << "Error with image " << names[i] << ": all images must be the same size!\n";
00055 exit(1);
00056 }
00057 }
00058 double mean, variance;
00059 tie(mean, variance) = mean_and_variance(ims);
00060 {
00061 for(unsigned int i=0; i < ims.size(); i++)
00062 transform(ims[i].begin(), ims[i].end(), ims[i].begin(), bind2nd(minus<double>(), mean));
00063 for(unsigned int i=0; i < ims.size(); i++)
00064 transform(ims[i].begin(), ims[i].end(), ims[i].begin(), bind1st(multiplies<double>(), 1/ sqrt(variance)));
00065 }
00066
00067 tie(mean, variance) = mean_and_variance(ims);
00068
00069 cerr << "Rescaled:\n";
00070 cerr << "mean = " << mean << endl;
00071 cerr << "std = " << sqrt(variance) << endl;
00072
00073
00074
00075
00076
00077
00078 ImageRef size = ims[0].size();
00079 Vector<10> p = Zeros;
00080 p[6]=-3;
00081 p[9]=-4;
00082
00083 Image<Vector<6> > monomials(size);
00084 Image<double> polynomial(size);
00085 for(int yy=0; yy < size.y; yy++)
00086 for(int xx=0; xx < size.x; xx++)
00087 {
00088 double x = xx *2./ size.x -1 ;
00089 double x2 = x*x;
00090 double y = yy *2./size.y - 1;
00091 double y2 = yy;
00092 monomials[yy][xx] = makeVector(1, x, y, x2, x*y, y2);
00093 }
00094
00095
00096 for(int i=0;i < 100;i++)
00097 {
00098 for(int yy=0; yy < size.y; yy++)
00099 for(int xx=0; xx < size.x; xx++)
00100 polynomial[yy][xx] = monomials[yy][xx] * p.slice<0,6>();
00101
00102 WLS<10, double, SQSVD> wls;
00103 for(unsigned int i=0; i < ims.size(); i++)
00104 for(int yy=0; yy < size.y; yy++)
00105 for(int xx=0; xx < size.x; xx++)
00106 {
00107 double t = i *1. / ims.size();
00108 double func = polynomial[yy][xx] * (exp(p[6]*t) + p[8]*exp(p[9]*t)) + p[7];
00109
00110 Vector<10> mJ;
00111
00112 mJ.slice<0,6>() = exp(p[6]*t)* monomials[yy][xx];
00113
00114 mJ[6] = polynomial[yy][xx] * exp(p[6]*t) * t;
00115
00116 mJ[7] = 1;
00117
00118 mJ[8] = polynomial[yy][xx] * exp(p[9]*t);
00119 mJ[9] = polynomial[yy][xx] * exp(p[9]*t) * t * p[8];
00120
00121 double err = ims[i][yy][xx] - func;
00122
00123 double w;
00124
00125
00126 if(err > 0)
00127 w = .01 / (abs(err) + .01);
00128 else
00129 w = 1;
00130
00131 wls.add_mJ(func - ims[i][yy][xx], -mJ, w);
00132 }
00133
00134 wls.add_prior(10);
00135 wls.compute();
00136
00137 p += wls.get_mu();
00138
00139 cout << p << endl << endl;
00140 }
00141
00142 for(unsigned int i=0; i < ims.size(); i++)
00143 for(int yy=0; yy < size.y; yy++)
00144 for(int xx=0; xx < size.x; xx++)
00145 {
00146 double x = xx *2./ size.x -1 ;
00147 double x2 = x*x;
00148 double y = yy *2./size.y - 1;
00149 double y2 = yy;
00150 double t = i *1. / ims.size();
00151 Vector<6> f = makeVector(1, x, y, x2, x*y, y2);
00152
00153 double func = f * p.slice<0,6>() * (exp(p[6]*t) + p[8]*exp(p[9]*t)) + p[7];
00154 ims[i][yy][xx] -= func;
00155 }
00156
00157 tie(mean, variance) = mean_and_variance(ims);
00158
00159
00160 cerr << "The mean should be small compared to std:\n";
00161 cerr << "mean = " << mean << endl;
00162 cerr << "std = " << sqrt(variance) << endl;
00163
00164
00165 {
00166 for(unsigned int i=0; i < ims.size(); i++)
00167 transform(ims[i].begin(), ims[i].end(), ims[i].begin(), bind1st(multiplies<double>(), 1/ sqrt(variance)));
00168 }
00169 tie(mean, variance) = mean_and_variance(ims);
00170
00171 cerr << "Rescaled:\n";
00172 cerr << "mean = " << mean << endl;
00173 cerr << "std = " << sqrt(variance) << endl;
00174
00175 return ims;
00176 }
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191 Image<float> preprocess_image(const Image<float>& im)
00192 {
00193 float wide = GV3::get<float>("preprocess.lpf", 0., -1);
00194 bool p = GV3::get<bool>("preprocess.skip", 0, -1);
00195
00196
00197 if(!p)
00198 {
00199 Image<float> f(im.size(), 0), fwide(im.size(), 0);
00200 convolveGaussian_fir(im, fwide, wide);
00201 for(int r=1; r < im.size().y-1; r++)
00202 for(int c=1; c < im.size().x-1; c++)
00203 f[r][c] = im[r][c] - fwide[r][c];
00204
00205 return f;
00206 }
00207 else
00208 return im.copy_from_me();
00209 }
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223 vector<Image<float> > load_and_preprocess_images(const vector<string>& names)
00224 {
00225 vector<Image<float> > ims;
00226
00227
00228
00229
00230 for(unsigned int i=0; i < names.size(); i++)
00231 {
00232 Image<float> im = img_load(names[i]);
00233
00234 ims.push_back(preprocess_image(im));
00235
00236 if(ims.back().size() != ims[0].size())
00237 {
00238 cerr << "Error with image " << names[i] << ": all images must be the same size!\n";
00239 exit(1);
00240 }
00241 }
00242 return ims;
00243 }
00244
00245
00246
00247
00248
00249
00250 pair<float, float> mean_and_variance(const vector<Image<float> >& images)
00251 {
00252 assert_same_size(images);
00253
00254 double sum=0, sum2=0, area=0;
00255
00256 for(unsigned int i=0; i < images.size(); i++)
00257 {
00258 area += images[i].size().area();
00259 for(int r=0; r < images[i].size().y; r++)
00260 for(int c=0; c < images[i].size().x; c++)
00261 {
00262 sum += images[i][r][c];
00263 sum2 += sq(images[i][r][c]);
00264 }
00265 }
00266
00267 sum /= area;
00268 sum2 /= area;
00269 return make_pair(sum, sum2 - sq(sum));
00270 }
00271
00272
00273
00274
00275
00276 pair<double, double> auto_fixed_scaling(const vector<Image<float> >& ims, double frac)
00277 {
00278 assert_same_size(ims);
00279
00280
00281 Image<double> ave(ims[0].size());
00282 ave.fill(0);
00283 for(unsigned int i=0; i < ims.size(); i++)
00284 for(int y=0; y < ave.size().y; y++)
00285 for(int x=0; x < ave.size().x; x++)
00286 ave[y][x] += ims[i][y][x];
00287
00288
00289 vector<pair<double, ImageRef> > pixels;
00290 for(int y=0; y < ave.size().y; y++)
00291 for(int x=0; x < ave.size().x; x++)
00292 pixels.push_back(make_pair(ave[y][x], ImageRef(x,y)));
00293
00294 int npix = (int) floor(frac *pixels.size() + 0.5);
00295 npix = max(0, min(npix, (int) pixels.size()));
00296
00297 nth_element(pixels.begin(), pixels.begin() + npix, pixels.end());
00298
00299 pixels.resize(npix);
00300
00301
00302 double sum=0, sum2=0;
00303
00304 for(unsigned int i=0; i < ims.size(); i++)
00305 {
00306 for(unsigned int j=0; j < pixels.size(); j++)
00307 {
00308 sum += ims[i][pixels[j].second];
00309 sum2 += sq(ims[i][pixels[j].second]);
00310 }
00311 }
00312
00313 double num = 1.0 * pixels.size() * ims.size();
00314 double mean = sum / num;
00315 double std = sqrt(((sum2/num) - mean*mean) * num / (num-1));
00316
00317 cout << "Automatic determination of fixed scaling:" << endl
00318 << "mean = " << mean << endl
00319 << "std = " << std << endl
00320 << "sqrt(mean) = " << sqrt(mean*255)/255 << endl;
00321
00322 return make_pair(mean, std);
00323 }
00324
00325
00326
00327
00328
00329
00330
00331 vector<Image<float> > load_and_normalize_images(const vector<string>& files)
00332 {
00333
00334 vector<Image<float> > ims = load_and_preprocess_images(files);
00335 double mean, variance;
00336 tie(mean, variance) = mean_and_variance(ims);
00337
00338 if(GV3::get<bool>("preprocess.fixed_scaling", 0, FATAL_IF_NOT_DEFINED))
00339 {
00340 bool skip = GV3::get<bool>("preprocess.skip");
00341 if(!skip)
00342 {
00343 cerr << "WARNING WARNING WARNING WARNING!!!!!!!!!!!!!!!\n";
00344 cerr << "preprocessing and fixed scaling selected!!!\n";
00345 exit(1);
00346 }
00347
00348 double sub, div;
00349 if(GV3::get<bool>("preprocess.fixed_scaling.auto", 0, FATAL_IF_NOT_DEFINED))
00350 {
00351 double frac = GV3::get<double>("preprocess.fixed_scaling.auto.proportion", 0, FATAL_IF_NOT_DEFINED);
00352 tie(sub, div) = auto_fixed_scaling(ims, frac);
00353 }
00354 else
00355 {
00356 sub = GV3::get<double>("preprocess.fixed_scaling.subtract", 0, FATAL_IF_NOT_DEFINED);
00357 div = GV3::get<double>("preprocess.fixed_scaling.divide", 0, FATAL_IF_NOT_DEFINED);
00358 }
00359
00360 for(unsigned int i=0; i < ims.size(); i++)
00361 for(Image<float>::iterator j=ims[i].begin(); j != ims[i].end(); j++)
00362 *j = (*j - sub)/div;
00363 }
00364 else
00365 {
00366
00367 cerr << "The mean should be small compared to std:\n";
00368 cerr << "mean = " << mean << endl;
00369 cerr << "std = " << sqrt(variance) << endl;
00370
00371
00372 {
00373 for(unsigned int i=0; i < ims.size(); i++)
00374 transform(ims[i].begin(), ims[i].end(), ims[i].begin(), bind1st(multiplies<double>(), 1/ sqrt(variance)));
00375 }
00376 }
00377
00378 tie(mean, variance) = mean_and_variance(ims);
00379
00380
00381 cerr << "Rescaled:\n";
00382 cerr << "mean = " << mean << endl;
00383 cerr << "std = " << sqrt(variance) << endl;
00384
00385 return ims;
00386 }
00387
00388
00389