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                                                const float imageStepSize,
00281                                                AppBase::ProgressDisplay* progress,
00282                                                double & error)
00283 {
00284 
00285     OptimizeVector photometricVars;
00286     // keep only the photometric variables
00287     unsigned int optCount=0;
00288     for (OptimizeVector::const_iterator it=vars.begin(); it != vars.end(); ++it)
00289     {
00290         std::set<std::string> cvars;
00291         for (std::set<std::string>::const_iterator itv = (*it).begin();
00292              itv != (*it).end(); ++itv)
00293         {
00294             if ((*itv)[0] == 'E' || (*itv)[0] == 'R' || (*itv)[0] == 'V') {
00295                 cvars.insert(*itv);
00296             }
00297         }
00298         photometricVars.push_back(cvars);
00299         optCount+=cvars.size();
00300     }
00301     //if no variables to optimize return
00302     if(optCount==0)
00303     {
00304         return;
00305     };
00306 
00307     int nMaxIter = 250;
00308     OptimData data(pano, photometricVars, correspondences, 5 * imageStepSize, false, nMaxIter, progress);
00309 
00310     double info[LM_INFO_SZ];
00311 
00312     // parameters
00313     int m=data.m_vars.size();
00314     vigra::ArrayVector<double> p(m, 0.0);
00315 
00316     // vector for errors
00317     int n=2*3*correspondences.size()+pano.getNrOfImages();
00318     vigra::ArrayVector<double> x(n, 0.0);
00319 
00320     data.ToX(p.begin());
00321 #ifdef DEBUG
00322     printf("Parameters before optimisation: ");
00323     for(int i=0; i<m; ++i)
00324         printf("%.7g ", p[i]);
00325     printf("\n");
00326 #endif
00327 
00328     // TODO: setup optimisation options with some good defaults.
00329     double optimOpts[5];
00330     
00331     optimOpts[0] = LM_INIT_MU;  // init mu
00332     // stop thresholds
00333     optimOpts[1] = LM_STOP_THRESH;   // ||J^T e||_inf
00334     optimOpts[2] = LM_STOP_THRESH;   // ||Dp||_2
00335     optimOpts[3] = std::pow(imageStepSize*0.1f, 2);   // ||e||_2
00336     // difference mode
00337     optimOpts[4] = LM_DIFF_DELTA;
00338     
00339     dlevmar_dif(&photometricError, &photometricVis, &(p[0]), &(x[0]), m, n, nMaxIter, optimOpts, info, NULL,NULL, &data);  // no jacobian
00340 
00341     // copy to source images (data.m_imgs)
00342     data.FromX(p.begin());
00343     // calculate error at solution
00344     data.huberSigma = 0;
00345     photometricError(&(p[0]), &(x[0]), m, n, &data);
00346     error = 0;
00347     for (int i=0; i<n; i++) {
00348         error += x[i]*x[i];
00349     }
00350     error = sqrt(error/n);
00351 
00352 #ifdef DEBUG
00353     printf("Levenberg-Marquardt returned in %g iter, reason %g\nSolution: ", info[5], info[6]);
00354     for(int i=0; i<m; ++i)
00355         printf("%.7g ", p[i]);
00356     printf("\n\nMinimization info:\n");
00357     for(int i=0; i<LM_INFO_SZ; ++i)
00358         printf("%g ", info[i]);
00359     printf("\n");
00360 #endif
00361 
00362     // copy settings to panorama
00363     for (unsigned i=0; i<pano.getNrOfImages(); i++) {
00364         pano.setSrcImage(i, data.m_imgs[i]);
00365     }
00366 }
00367 
00368 bool IsHighVignetting(std::vector<double> vigCorr)
00369 {
00370     SrcPanoImage srcImage;
00371     srcImage.setRadialVigCorrCoeff(vigCorr);
00372     srcImage.setSize(vigra::Size2D(500, 500));
00373     Photometric::ResponseTransform<double> transform(srcImage);
00374     transform.enforceMonotonicity();
00375     return transform(1.0, hugin_utils::FDiff2D(0,0))<0.7;
00376 };
00377 
00378 bool CheckStrangeWB(PanoramaData & pano)
00379 {
00380     for(size_t i=0; i<pano.getNrOfImages(); i++)
00381     {
00382         if(pano.getImage(i).getWhiteBalanceBlue()>3)
00383         {
00384             return true;
00385         };
00386         if(pano.getImage(i).getWhiteBalanceRed()>3)
00387         {
00388             return true;
00389         };
00390     };
00391     return false;
00392 };
00393 
00394 void SmartPhotometricOptimizer::smartOptimizePhotometric(PanoramaData & pano, PhotometricOptimizeMode mode,
00395                                                           const std::vector<vigra_ext::PointPairRGB> & correspondences,
00396                                                           const float imageStepSize, 
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, imageStepSize, 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, imageStepSize, 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, imageStepSize, 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, imageStepSize, progress, error);
00465         };
00466     };
00467 }
00468 
00469 
00470 bool PhotometricOptimizer::runAlgorithm()
00471 {
00472     optimizePhotometric(o_panorama, 
00473                         o_vars, o_correspondences, o_imageStepSize,
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, o_imageStepSize,
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 7 Dec 2016 for Hugintrunk by  doxygen 1.4.7