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

Generated on 29 Jun 2016 for Hugintrunk by  doxygen 1.4.7