PhotometricOptimizer.cpp

Go to the documentation of this file.
00001 // -*- c-basic-offset: 4 -*-
00024 #include "PhotometricOptimizer.h"
00025 
00026 #include <fstream>
00027 #include <foreign/levmar/levmar.h>
00028 #include <photometric/ResponseTransform.h>
00029 #include <algorithms/basic/LayerStacks.h>
00030 
00031 #ifdef DEBUG
00032 #define DEBUG_LOG_VIG 1
00033 #endif
00034 
00035 
00036 namespace HuginBase {
00037     
00038 using namespace std;
00039 
00040 
00042 inline double weightHuber(double x, double sigma)
00043 {
00044     if (x > sigma) {
00045         x = sqrt(sigma* (2*x - sigma));
00046     }
00047     return x;
00048 }
00049 
00050 
00051 
00052 PhotometricOptimizer::OptimData::OptimData(const PanoramaData & pano, const OptimizeVector & optvars,
00053                                            const std::vector<vigra_ext::PointPairRGB> & data,
00054                                            double mEstimatorSigma, bool symmetric,
00055                                            int maxIter, AppBase::ProgressDisplay* progress)
00056   : m_pano(pano), m_data(data), huberSigma(mEstimatorSigma), symmetricError(symmetric),
00057     m_maxIter(maxIter), m_progress(progress)
00058 {
00059     assert(pano.getNrOfImages() == optvars.size());
00060         
00061     for (unsigned i=0; i < pano.getNrOfImages(); i++) {
00062         m_imgs.push_back(pano.getSrcImage(i));
00063     }
00064 
00065     std::vector<std::set<std::string> > usedVars(pano.getNrOfImages());
00066 
00067     // create variable map with param <-> var assignments
00068     for (unsigned i=0; i < optvars.size(); i++) 
00069     {
00070         const std::set<std::string> vars = optvars[i];
00071         const SrcPanoImage & img_i = pano.getImage(i);
00072         for (std::set<std::string>::const_iterator it = vars.begin();
00073              it != vars.end(); ++it)
00074         {
00075             VarMapping var;
00076             var.type = *it;
00077             //check if variable is yet included
00078             if(set_contains(usedVars[i],var.type))
00079                 continue;
00080             var.imgs.insert(i);
00081             usedVars[i].insert(var.type);
00082             //now check all linked images and add image nr
00083 #define CheckLinked(name)\
00084             if(img_i.name##isLinked())\
00085             {\
00086                 for(unsigned j=i+1;j<pano.getNrOfImages();j++)\
00087                     if(img_i.name##isLinkedWith(pano.getImage(j)))\
00088                     {\
00089                         var.imgs.insert(j);\
00090                         usedVars[j].insert(var.type);\
00091                     };\
00092             }
00093 
00094             if(var.type=="Eev")
00095             {
00096                 CheckLinked(ExposureValue)
00097             };
00098             if(var.type=="Er")
00099             {
00100                 CheckLinked(WhiteBalanceRed)
00101             };
00102             if(var.type=="Eb")
00103             {
00104                 CheckLinked(WhiteBalanceBlue)
00105             };
00106             if(var.type[0]=='R')
00107             {
00108                 CheckLinked(EMoRParams)
00109             };
00110             if(var.type=="Va" || var.type=="Vb" || var.type=="Vc" || var.type=="Vd")
00111             {
00112                 CheckLinked(RadialVigCorrCoeff)
00113             }
00114             if(var.type=="Vx" || var.type=="Vy")
00115             {
00116                 CheckLinked(RadialVigCorrCenterShift)
00117             };
00118 #undef CheckLinked
00119             m_vars.push_back(var);
00120         }
00121     }
00122 }
00123 
00124 void PhotometricOptimizer::OptimData::ToX(double * x)
00125 {
00126     for (size_t i=0; i < m_vars.size(); i++)
00127     {
00128         assert(m_vars[i].imgs.size() > 0);
00129             // get corresponding image number
00130         unsigned j = *(m_vars[i].imgs.begin());
00131             // get value
00132         x[i] = m_imgs[j].getVar(m_vars[i].type);
00133         // TODO: transform some variables, such as the vignetting center!
00134     }
00135 }
00136 
00137 
00138 void PhotometricOptimizer::OptimData::FromX(double * x)
00139 {
00140     for (size_t i=0; i < m_vars.size(); i++)
00141     {
00142         // TODO: transform some variables, such as the vignetting center!
00143         assert(m_vars[i].imgs.size() > 0);
00144             // copy value int all images
00145         for (std::set<unsigned>::const_iterator it = m_vars[i].imgs.begin();
00146              it != m_vars[i].imgs.end(); ++it)
00147         {
00148             m_imgs[*it].setVar(m_vars[i].type, x[i]);
00149         }
00150     }
00151 }
00152 
00153 
00154 
00155 void PhotometricOptimizer::photometricError(double *p, double *x, int m, int n, void * data)
00156 {
00157 #ifdef DEBUG_LOG_VIG
00158     static int iter = 0;
00159 #endif
00160     typedef Photometric::ResponseTransform<vigra::RGBValue<double> > RespFunc;
00161     typedef Photometric::InvResponseTransform<vigra::RGBValue<double>, vigra::RGBValue<double> > InvRespFunc;
00162 
00163     int xi = 0 ;
00164 
00165     OptimData * dat = static_cast<OptimData*>(data);
00166     dat->FromX(p);
00167 #ifdef DEBUG_LOG_VIG
00168     ostringstream oss;
00169     oss << "vig_log_" << iter;
00170     iter++;
00171     ofstream log(oss.str().c_str());
00172     log << "VIGparams = [";
00173     for (int i = 0; i < m; i++) {
00174         log << p[i] << " ";
00175     }
00176     log << " ]; " << std::endl;
00177     // TODO: print parameters of images.
00178     std::ofstream script("vig_test.pto");
00179     OptimizeVector optvars(dat->m_pano.getNrOfImages());
00180     UIntSet imgs = dat->m_pano.getActiveImages();
00181     dat->m_pano.printPanoramaScript(script, optvars, dat->m_pano.getOptions(), imgs, false, "");
00182 #endif
00183 
00184     size_t nImg = dat->m_imgs.size();
00185     std::vector<RespFunc> resp(nImg);
00186     std::vector<InvRespFunc> invResp(nImg);
00187     for (size_t i=0; i < nImg; i++) {
00188         resp[i] = RespFunc(dat->m_imgs[i]);
00189         invResp[i] = InvRespFunc(dat->m_imgs[i]);
00190         // calculate the monotonicity error
00191         double monErr = 0;
00192         if (dat->m_imgs[i].getResponseType() == SrcPanoImage::RESPONSE_EMOR) {
00193             // calculate monotonicity error
00194             int lutsize = resp[i].m_lutR.size();
00195             for (int j=0; j < lutsize-1; j++)
00196             {
00197                 double d = resp[i].m_lutR[j] - resp[i].m_lutR[j+1];
00198                 if (d > 0) {
00199                     monErr += d*d*lutsize;
00200                 }
00201             }
00202         }
00203         x[xi++] = monErr;
00204                 // enforce a montonous response curves
00205                 resp[i].enforceMonotonicity();
00206                 invResp[i].enforceMonotonicity();
00207     }
00208 
00209 #ifdef DEBUG
00210     double sqerror=0;
00211 #endif
00212     // loop over all points to calculate the error
00213 #ifdef DEBUG_LOG_VIG
00214     log << "VIGval = [ ";
00215 #endif
00216 
00217     for (std::vector<vigra_ext::PointPairRGB>::const_iterator it = dat->m_data.begin();
00218          it != dat->m_data.end(); ++it)
00219     {
00220         vigra::RGBValue<double> l2 = invResp[it->imgNr2](it->i2, it->p2);
00221         vigra::RGBValue<double> i2ini1 = resp[it->imgNr1](l2, it->p1);
00222         vigra::RGBValue<double> error = it->i1 - i2ini1;
00223 
00224 
00225         // if requested, calcuate the error in image 2 as well.
00226         //TODO: weighting dependent on the pixel value? check if outside of i2 range?
00227         vigra::RGBValue<double> l1 = invResp[it->imgNr1](it->i1, it->p1);
00228         vigra::RGBValue<double> i1ini2 = resp[it->imgNr2](l1, it->p2);
00229         vigra::RGBValue<double> error2 = it->i2 - i1ini2;
00230 
00231 #ifdef DEBUG
00232         for (int i=0; i < 3; i++) {
00233             sqerror += error[i]*error[i];
00234             sqerror += error2[i]*error2[i];
00235         }
00236 #endif
00237 
00238         // use huber robust estimator
00239         if (dat->huberSigma > 0) {
00240             for (int i=0; i < 3; i++) {
00241                 x[xi++] = weightHuber(fabs(error[i]), dat->huberSigma);
00242                 x[xi++] = weightHuber(fabs(error2[i]), dat->huberSigma);
00243             }
00244         } else {
00245             x[xi++] = error[0];
00246             x[xi++] = error[1];
00247             x[xi++] = error[2];
00248             x[xi++] = error2[0];
00249             x[xi++] = error2[1];
00250             x[xi++] = error2[2];
00251         }
00252 
00253 #ifdef DEBUG_LOG_VIG
00254         log << it->i1.green()  << " "<< l1.green()  << " " << i1ini2.green() << "   " 
00255              << it->i2.green()  << " "<< l2.green()  << " " << i2ini1.green() << ";  " << std::endl;
00256 #endif
00257 
00258     }
00259 #ifdef DEBUG_LOG_VIG
00260     log << std::endl << "VIGerr = [";
00261     for (int i = 0; i < n; i++) {
00262         log << x[i] << std::endl;
00263     }
00264     log << " ]; " << std::endl;
00265 #endif
00266 #ifdef DEBUG
00267     DEBUG_DEBUG("squared error: " << sqerror);
00268 #endif
00269 }
00270 
00271 int PhotometricOptimizer::photometricVis(double *p, double *x, int m, int n, int iter, double sqerror, void * data)
00272 {
00273     OptimData * dat = static_cast<OptimData*>(data);
00274     char tmp[200];
00275     tmp[199] = 0;
00276     double error = sqrt(sqerror/n)*255;
00277     snprintf(tmp,199, "Iteration: %d, error: %f", iter, error);
00278     return dat->m_progress->updateDisplay(std::string(tmp)) ? 1 : 0 ;
00279 }
00280 
00281 void PhotometricOptimizer::optimizePhotometric(PanoramaData & pano, const OptimizeVector & vars,
00282                                                const std::vector<vigra_ext::PointPairRGB> & correspondences,
00283                                                AppBase::ProgressDisplay* progress,
00284                                                double & error)
00285 {
00286 
00287     OptimizeVector photometricVars;
00288     // keep only the photometric variables
00289     unsigned int optCount=0;
00290     for (OptimizeVector::const_iterator it=vars.begin(); it != vars.end(); ++it)
00291     {
00292         std::set<std::string> cvars;
00293         for (std::set<std::string>::const_iterator itv = (*it).begin();
00294              itv != (*it).end(); ++itv)
00295         {
00296             if ((*itv)[0] == 'E' || (*itv)[0] == 'R' || (*itv)[0] == 'V') {
00297                 cvars.insert(*itv);
00298             }
00299         }
00300         photometricVars.push_back(cvars);
00301         optCount+=cvars.size();
00302     }
00303     //if no variables to optimize return
00304     if(optCount==0)
00305     {
00306         return;
00307     };
00308 
00309     int nMaxIter = 250;
00310     OptimData data(pano, photometricVars, correspondences, 5/255.0, false, nMaxIter, progress);
00311 
00312     double info[LM_INFO_SZ];
00313 
00314     // parameters
00315     int m=data.m_vars.size();
00316     vigra::ArrayVector<double> p(m, 0.0);
00317 
00318     // vector for errors
00319     int n=2*3*correspondences.size()+pano.getNrOfImages();
00320     vigra::ArrayVector<double> x(n, 0.0);
00321 
00322     data.ToX(p.begin());
00323 #ifdef DEBUG
00324     printf("Parameters before optimisation: ");
00325     for(int i=0; i<m; ++i)
00326         printf("%.7g ", p[i]);
00327     printf("\n");
00328 #endif
00329 
00330     // TODO: setup optimisation options with some good defaults.
00331     double optimOpts[5];
00332     
00333     optimOpts[0] = 1E-03;  // init mu
00334     // stop thresholds
00335     optimOpts[1] = 1e-5;   // ||J^T e||_inf
00336     optimOpts[2] = 1e-5;   // ||Dp||_2
00337     optimOpts[3] = 1e-1;   // ||e||_2
00338     // difference mode
00339     optimOpts[4] = LM_DIFF_DELTA;
00340     
00341     dlevmar_dif(&photometricError, &photometricVis, &(p[0]), &(x[0]), m, n, nMaxIter, optimOpts, info, NULL,NULL, &data);  // no jacobian
00342     // copy to source images (data.m_imgs)
00343     data.FromX(p.begin());
00344     // calculate error at solution
00345     data.huberSigma = 0;
00346     photometricError(&(p[0]), &(x[0]), m, n, &data);
00347     error = 0;
00348     for (int i=0; i<n; i++) {
00349         error += x[i]*x[i];
00350     }
00351     error = sqrt(error/n);
00352 
00353 #ifdef DEBUG
00354     printf("Levenberg-Marquardt returned %d in %g iter, reason %g\nSolution: ", ret, info[5], info[6]);
00355     for(int i=0; i<m; ++i)
00356         printf("%.7g ", p[i]);
00357     printf("\n\nMinimization info:\n");
00358     for(int i=0; i<LM_INFO_SZ; ++i)
00359         printf("%g ", info[i]);
00360     printf("\n");
00361 #endif
00362 
00363     // copy settings to panorama
00364     for (unsigned i=0; i<pano.getNrOfImages(); i++) {
00365         pano.setSrcImage(i, data.m_imgs[i]);
00366     }
00367 }
00368 
00369 bool IsHighVignetting(std::vector<double> vigCorr)
00370 {
00371     SrcPanoImage srcImage;
00372     srcImage.setRadialVigCorrCoeff(vigCorr);
00373     srcImage.setSize(vigra::Size2D(500, 500));
00374     Photometric::ResponseTransform<double> transform(srcImage);
00375     transform.enforceMonotonicity();
00376     return transform(1.0, hugin_utils::FDiff2D(0,0))<0.7;
00377 };
00378 
00379 bool CheckStrangeWB(PanoramaData & pano)
00380 {
00381     for(size_t i=0; i<pano.getNrOfImages(); i++)
00382     {
00383         if(pano.getImage(i).getWhiteBalanceBlue()>3)
00384         {
00385             return true;
00386         };
00387         if(pano.getImage(i).getWhiteBalanceRed()>3)
00388         {
00389             return true;
00390         };
00391     };
00392     return false;
00393 };
00394 
00395 void SmartPhotometricOptimizer::smartOptimizePhotometric(PanoramaData & pano, PhotometricOptimizeMode mode,
00396                                                           const std::vector<vigra_ext::PointPairRGB> & correspondences,
00397                                                           AppBase::ProgressDisplay* progress,
00398                                                           double & error)
00399 {
00400     PanoramaOptions opts = pano.getOptions();
00401     UIntSet images;
00402     fill_set(images, 0, pano.getNrOfImages()-1);
00403     std::vector<UIntSet> stacks = getHDRStacks(pano, images, pano.getOptions());
00404     const bool singleStack = (stacks.size() == 1);
00405 
00406     int vars = 0;
00407     if (mode == OPT_PHOTOMETRIC_LDR || mode == OPT_PHOTOMETRIC_LDR_WB)
00408     {
00409         // optimize exposure
00410         vars = OPT_EXP;
00411         optimizePhotometric(pano, 
00412                             createOptVars(pano, vars, opts.colorReferenceImage),
00413                             correspondences, progress, error);
00414     };
00415 
00416     if(!singleStack)
00417     {
00418         //optimize vignetting only if there are more than 1 image stack
00419         // for a single stack vignetting can't be calculated by this optimization
00420         vars |= OPT_VIG;
00421         VariableMapVector oldVars = pano.getVariables();
00422         optimizePhotometric(pano, 
00423                             createOptVars(pano, vars, opts.colorReferenceImage),
00424                             correspondences, progress, error);
00425         // check if vignetting is plausible
00426         if(IsHighVignetting(pano.getImage(0).getRadialVigCorrCoeff()))
00427         {
00428             vars &= ~OPT_VIG;
00429             pano.updateVariables(oldVars);
00430         };
00431     };
00432 
00433     // now take response curve into account
00434     vars |= OPT_RESP;
00435     // also WB if desired
00436     if (mode == OPT_PHOTOMETRIC_LDR_WB || mode == OPT_PHOTOMETRIC_HDR_WB)
00437     {
00438         vars |= OPT_WB;
00439     };
00440     VariableMapVector oldVars = pano.getVariables();
00441     optimizePhotometric(pano, 
00442                         createOptVars(pano, vars, opts.colorReferenceImage),
00443                         correspondences, progress, error);
00444     // now check the results
00445     const bool hasHighVignetting = IsHighVignetting(pano.getImage(0).getRadialVigCorrCoeff());
00446     // @TODO check also response curve, what parameters are considered invalid?
00447     const bool hasHighWB = CheckStrangeWB(pano);
00448     if(hasHighVignetting)
00449     {
00450         vars &= ~OPT_VIG;
00451     };
00452     if(hasHighWB)
00453     {
00454         vars &= ~OPT_WB;
00455     };
00456     if(hasHighVignetting || hasHighWB)
00457     {
00458         // we got strange results, optimize again with less parameters
00459         pano.updateVariables(oldVars);
00460         if(vars>0)
00461         {
00462             optimizePhotometric(pano, 
00463                                 createOptVars(pano, vars, opts.colorReferenceImage),
00464                                 correspondences, progress, error);
00465         };
00466     };
00467 }
00468 
00469 
00470 bool PhotometricOptimizer::runAlgorithm()
00471 {
00472     optimizePhotometric(o_panorama, 
00473                         o_vars, o_correspondences,
00474                         getProgressDisplay(), o_resultError);
00475     
00476     // optimizePhotometric does not tell us if it's cancelled
00477     if (getProgressDisplay()->wasCancelled())
00478     {
00479         cancelAlgorithm();
00480     }
00481     
00482     return wasCancelled(); // let's hope so.
00483 }
00484 
00485 bool SmartPhotometricOptimizer::runAlgorithm()
00486 {
00487     smartOptimizePhotometric(o_panorama,
00488                              o_optMode,
00489                              o_correspondences,
00490                              getProgressDisplay(), o_resultError);
00491     
00492     // smartOptimizePhotometric does not tell us if it's cancelled
00493     if (getProgressDisplay()->wasCancelled())
00494     {
00495         cancelAlgorithm();
00496     };
00497     
00498     return !wasCancelled(); // let's hope so.
00499 }
00500 
00501 } //namespace

Generated on 30 Jul 2015 for Hugintrunk by  doxygen 1.4.7