ThreeB 1.1

- Main Page
- Related Pages
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

storm_imagery.cc
Go to the documentation of this file.

    00001 /*
    00002     This file is part of B-cubed.
    00003
    00004     Copyright (C) 2009, 2010, 2011, Edward Rosten and Susan Cox
    00005
    00006     B-cubed is free software; you can redistribute it and/or
    00007     modify it under the terms of the GNU Lesser General Public
    00008     License as published by the Free Software Foundation; either
    00009     version 3.0 of the License, or (at your option) any later version.
    00010
    00011     This library is distributed in the hope that it will be useful,
    00012     but WITHOUT ANY WARRANTY; without even the implied warranty of
    00013     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    00014     Lesser General Public License for more details.
    00015
    00016     You should have received a copy of the GNU General Public License
    00017     along with this program.  If not, see <http://www.gnu.org/licenses/>
    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 <tr1/tuple>
    00025
    00026 #include "storm_imagery.h"
    00027 #include "debug.h"
    00028 #include "utility.h"
    00029
    00030 using namespace CVD;
    00031 using namespace TooN;
    00032 using namespace GVars3;
    00033 using namespace std;
    00034 using namespace std::tr1;
    00035
    00036 /**Load all images from disk and do the initial preprocessing.
    00037
    00038 @param names List of filenames to load.
    00039 @return preprocessed images.
    00040 @ingroup gStormImages
    00041 **/
    00042 vector<Image<float> > load_and_preprocess_images2(const vector<string>& names)
    00043 {
    00044     vector<Image<float> > ims;
    00045     //Load images
    00046     for(unsigned int i=0; i < names.size(); i++)
    00047     {
    00048         Image<float> im = img_load(names[i]);
    00049         ims.push_back(im);
    00050
    00051         if(ims.back().size() != ims[0].size())
    00052         {
    00053             cerr << "Error with image " << names[i] << ":  all images must be the same size!\n";
    00054             exit(1);
    00055         }
    00056     }
    00057     double mean, variance;
    00058     tie(mean, variance) = mean_and_variance(ims);
    00059     {
    00060         for(unsigned int i=0; i < ims.size(); i++)
    00061             transform(ims[i].begin(), ims[i].end(), ims[i].begin(), bind2nd(minus<double>(), mean));
    00062         for(unsigned int i=0; i < ims.size(); i++)
    00063             transform(ims[i].begin(), ims[i].end(), ims[i].begin(), bind1st(multiplies<double>(), 1/ sqrt(variance)));
    00064     }
    00065
    00066     tie(mean, variance) = mean_and_variance(ims);
    00067
    00068     cerr << "Rescaled:\n";
    00069     cerr << "mean = " << mean << endl;
    00070     cerr << "std  = " << sqrt(variance) << endl;
    00071
    00072
    00073
    00074     //Normalize...
    00075
    00076     //Fit the background model
    00077     ImageRef size = ims[0].size();
    00078     Vector<10> p = Zeros;
    00079     p[6]=-3;
    00080     p[9]=-4;
    00081
    00082     Image<Vector<6> > monomials(size);
    00083     Image<double> polynomial(size);
    00084     for(int yy=0; yy < size.y; yy++)
    00085         for(int xx=0; xx < size.x; xx++)
    00086         {
    00087             double x = xx *2./ size.x -1 ;
    00088             double x2 = x*x;
    00089             double y = yy *2./size.y - 1;
    00090             double y2 = yy;
    00091             monomials[yy][xx] = makeVector(1, x, y, x2, x*y, y2);
    00092         }
    00093
    00094
    00095     for(int i=0;i < 100;i++)
    00096     {
    00097         for(int yy=0; yy < size.y; yy++)
    00098             for(int xx=0; xx < size.x; xx++)
    00099                 polynomial[yy][xx] = monomials[yy][xx] *  p.slice<0,6>();
    00100
    00101         WLS<10> wls;
    00102         for(unsigned int i=0; i < ims.size(); i++)
    00103             for(int yy=0; yy < size.y; yy++)
    00104                 for(int xx=0; xx < size.x; xx++)
    00105                 {
    00106                     double t = i *1. / ims.size();
    00107                     double func = polynomial[yy][xx] * (exp(p[6]*t) + p[8]*exp(p[9]*t)) + p[7];
    00108
    00109                     Vector<10> mJ;
    00110
    00111                     mJ.slice<0,6>() = exp(p[6]*t)* monomials[yy][xx];
    00112                     //mJ.slice<3,3>() = Zeros;
    00113                     mJ[6] = polynomial[yy][xx] * exp(p[6]*t) * t;
    00114                     //mJ[6] = func  * t;
    00115                     mJ[7] = 1;
    00116
    00117                     mJ[8] = polynomial[yy][xx] * exp(p[9]*t);
    00118                     mJ[9] = polynomial[yy][xx] * exp(p[9]*t) * t * p[8];
    00119
    00120                     double err = ims[i][yy][xx] - func;
    00121
    00122                     double w;
    00123
    00124
    00125                     if(err > 0)
    00126                         w = .01 / (abs(err) + .01);
    00127                     else
    00128                         w = 1;
    00129
    00130                     wls.add_mJ(func - ims[i][yy][xx], -mJ, w);
    00131                 }
    00132
    00133         wls.add_prior(10);
    00134         wls.compute();
    00135
    00136         p += wls.get_mu();
    00137
    00138         cout << p << endl << endl;
    00139     }
    00140
    00141     for(unsigned int i=0; i < ims.size(); i++)
    00142         for(int yy=0; yy < size.y; yy++)
    00143             for(int xx=0; xx < size.x; xx++)
    00144             {
    00145                 double x = xx *2./ size.x -1 ;
    00146                 double x2 = x*x;
    00147                 double y = yy *2./size.y - 1;
    00148                 double y2 = yy;
    00149                 double t = i *1. / ims.size();
    00150                 Vector<6> f = makeVector(1, x, y, x2, x*y, y2);
    00151
    00152                 double func = f * p.slice<0,6>() * (exp(p[6]*t) + p[8]*exp(p[9]*t)) + p[7];
    00153                 ims[i][yy][xx] -= func;
    00154             }
    00155
    00156     tie(mean, variance) = mean_and_variance(ims);
    00157
    00158     //A sanity check.
    00159     cerr << "The mean should be small compared to std:\n";
    00160     cerr << "mean = " << mean << endl;
    00161     cerr << "std  = " << sqrt(variance) << endl;
    00162
    00163     //Scale by the variance.
    00164     {
    00165         for(unsigned int i=0; i < ims.size(); i++)
    00166             transform(ims[i].begin(), ims[i].end(), ims[i].begin(), bind1st(multiplies<double>(), 1/ sqrt(variance)));
    00167     }
    00168     tie(mean, variance) = mean_and_variance(ims);
    00169
    00170     cerr << "Rescaled:\n";
    00171     cerr << "mean = " << mean << endl;
    00172     cerr << "std  = " << sqrt(variance) << endl;
    00173
    00174     return ims;
    00175 }
    00176
    00177
    00178
    00179
    00180 /**Load all images from disk and do the initial preprocessing. Currently
    00181 this is a high pass filter to make the resultimg images zero mean.
    00182
    00183 The filter is controlled with the \c preprocess.lpf and \c preprocess.skip Gvars
    00184
    00185 See also load_and_preprocess_image()
    00186
    00187 @param names List of filenames to load.
    00188 @return preprocessed images.
    00189 @ingroup gStormImages
    00190 **/
    00191 vector<Image<float> > load_and_preprocess_images(const vector<string>& names)
    00192 {
    00193     vector<Image<float> > ims;
    00194
    00195     //float wide = GV3::get<float>("preprocess.lpf", 0., -1);
    00196     //bool p = GV3::get<bool>("preprocess.skip", 0, -1);
    00197
    00198     for(unsigned int i=0; i < names.size(); i++)
    00199     {
    00200         Image<float> im = img_load(names[i]);
    00201
    00202         ims.push_back(preprocess_image(im));
    00203
    00204         if(ims.back().size() != ims[0].size())
    00205         {
    00206             cerr << "Error with image " << names[i] << ":  all images must be the same size!\n";
    00207             exit(1);
    00208         }
    00209     }
    00210     return ims;
    00211 }
    00212
    00213 /**Compute the mean and variance of the (on average) darkest pixels, in order
    00214 to find the correct scaling, by examining hte background.
    00215 */
    00216 pair<double, double> auto_fixed_scaling(const vector<Image<float> >& ims, double frac)
    00217 {
    00218     assert_same_size(ims);
    00219
    00220     //Compute the mean image (ish)
    00221     Image<double> ave(ims[0].size());
    00222     ave.fill(0);
    00223     for(unsigned int i=0; i < ims.size(); i++)
    00224         for(int y=0; y < ave.size().y; y++)
    00225             for(int x=0; x < ave.size().x; x++)
    00226                 ave[y][x] += ims[i][y][x];
    00227
    00228     //Find the smallest N% of the pixels...
    00229     vector<pair<double, ImageRef> > pixels;
    00230     for(int y=0; y < ave.size().y; y++)
    00231         for(int x=0; x < ave.size().x; x++)
    00232             pixels.push_back(make_pair(ave[y][x], ImageRef(x,y)));
    00233
    00234     int npix = (int) floor(frac *pixels.size() + 0.5);
    00235     npix = max(0, min(npix, (int) pixels.size()));
    00236
    00237     nth_element(pixels.begin(), pixels.begin() + npix, pixels.end());
    00238
    00239     pixels.resize(npix);
    00240
    00241     //Now compute the mean and variance of those pixels.
    00242     double sum=0, sum2=0;
    00243
    00244     for(unsigned int i=0; i < ims.size(); i++)
    00245     {
    00246         for(unsigned int j=0; j < pixels.size(); j++)
    00247         {
    00248             sum += ims[i][pixels[j].second];
    00249             sum2 += sq(ims[i][pixels[j].second]);
    00250         }
    00251     }
    00252
    00253     double num = 1.0 * pixels.size() * ims.size();
    00254     double mean = sum / num;
    00255     double std  = sqrt(((sum2/num) - mean*mean) * num / (num-1));
    00256
    00257     cout << "Automatic determination of fixed scaling:" << endl
    00258          << "mean       = " << mean << endl
    00259          << "std        = " << std << endl
    00260          << "sqrt(mean) = " << sqrt(mean*255)/255 << endl;
    00261
    00262     return make_pair(mean, std);
    00263 }
    00264
    00265 /**Wrapper for load_and_preprocess_images() to allow more flexible behaviour.
    00266
    00267 @param files List of filenames to load.
    00268 @return preprocessed images.
    00269 @ingroup gStormImages
    00270 **/
    00271 vector<Image<float> > load_and_normalize_images(const vector<string>& files)
    00272 {
    00273     //Load the raw data, and then load the spot parameters.
    00274     vector<Image<float> > ims = load_and_preprocess_images(files);
    00275     double mean, variance;
    00276     tie(mean, variance) = mean_and_variance(ims);
    00277
    00278     if(GV3::get<bool>("preprocess.fixed_scaling", 0, FATAL_IF_NOT_DEFINED))
    00279     {
    00280         bool skip = GV3::get<bool>("preprocess.skip");
    00281         if(!skip)
    00282         {
    00283             cerr << "WARNING WARNING WARNING WARNING!!!!!!!!!!!!!!!\n";
    00284             cerr << "preprocessing and fixed scaling selected!!!\n";
    00285             exit(1);
    00286         }
    00287
    00288         double sub, div;
    00289         if(GV3::get<bool>("preprocess.fixed_scaling.auto", 0, FATAL_IF_NOT_DEFINED))
    00290         {
    00291             double frac = GV3::get<double>("preprocess.fixed_scaling.auto.proportion", 0, FATAL_IF_NOT_DEFINED);
    00292             tie(sub, div) = auto_fixed_scaling(ims, frac);
    00293         }
    00294         else
    00295         {
    00296             sub = GV3::get<double>("preprocess.fixed_scaling.subtract", 0, FATAL_IF_NOT_DEFINED);
    00297             div = GV3::get<double>("preprocess.fixed_scaling.divide", 0, FATAL_IF_NOT_DEFINED);
    00298         }
    00299
    00300         for(unsigned int i=0; i < ims.size(); i++)
    00301             for(Image<float>::iterator j=ims[i].begin(); j != ims[i].end(); j++)
    00302                 *j = (*j - sub)/div;
    00303     }
    00304     else
    00305     {
    00306         //A sanity check.
    00307         cerr << "The mean should be small compared to std:\n";
    00308         cerr << "mean = " << mean << endl;
    00309         cerr << "std  = " << sqrt(variance) << endl;
    00310
    00311         //Scale by the variance.
    00312         {
    00313             for(unsigned int i=0; i < ims.size(); i++)
    00314                 transform(ims[i].begin(), ims[i].end(), ims[i].begin(), bind1st(multiplies<double>(), 1/ sqrt(variance)));
    00315         }
    00316     }
    00317
    00318     tie(mean, variance) = mean_and_variance(ims);
    00319
    00320     //A sanity check.
    00321     cerr << "Rescaled:\n";
    00322     cerr << "mean = " << mean << endl;
    00323     cerr << "std  = " << sqrt(variance) << endl;
    00324
    00325     return ims;
    00326 }
    00327
    00328
    00329