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.empty());
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.empty());
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     for (size_t x = 0; x < 250; x += 10)
00375     {
00376         const double vigFactor = transform.calcVigFactor(hugin_utils::FDiff2D(x, x));
00377         if (vigFactor<0.7 || vigFactor>1.1)
00378         {
00379             return true;
00380         };
00381     };
00382     return false;
00383 };
00384 
00385 bool CheckStrangeWB(PanoramaData & pano)
00386 {
00387     for(size_t i=0; i<pano.getNrOfImages(); i++)
00388     {
00389         if(pano.getImage(i).getWhiteBalanceBlue()>3)
00390         {
00391             return true;
00392         };
00393         if(pano.getImage(i).getWhiteBalanceRed()>3)
00394         {
00395             return true;
00396         };
00397     };
00398     return false;
00399 };
00400 
00401 void SmartPhotometricOptimizer::smartOptimizePhotometric(PanoramaData & pano, PhotometricOptimizeMode mode,
00402                                                           const std::vector<vigra_ext::PointPairRGB> & correspondences,
00403                                                           const float imageStepSize, 
00404                                                           AppBase::ProgressDisplay* progress,
00405                                                           double & error)
00406 {
00407     PanoramaOptions opts = pano.getOptions();
00408     UIntSet images;
00409     fill_set(images, 0, pano.getNrOfImages()-1);
00410     std::vector<UIntSet> stacks = getHDRStacks(pano, images, pano.getOptions());
00411     const bool singleStack = (stacks.size() == 1);
00412 
00413     int vars = 0;
00414     if (mode == OPT_PHOTOMETRIC_LDR || mode == OPT_PHOTOMETRIC_LDR_WB)
00415     {
00416         // optimize exposure
00417         vars = OPT_EXP;
00418         optimizePhotometric(pano, 
00419                             createOptVars(pano, vars, opts.colorReferenceImage),
00420                             correspondences, imageStepSize, progress, error);
00421     };
00422 
00423     if(!singleStack)
00424     {
00425         //optimize vignetting only if there are more than 1 image stack
00426         // for a single stack vignetting can't be calculated by this optimization
00427         vars |= OPT_VIG;
00428         VariableMapVector oldVars = pano.getVariables();
00429         optimizePhotometric(pano, 
00430                             createOptVars(pano, vars, opts.colorReferenceImage),
00431                             correspondences, imageStepSize, progress, error);
00432         // check if vignetting is plausible
00433         if(IsHighVignetting(pano.getImage(0).getRadialVigCorrCoeff()))
00434         {
00435             vars &= ~OPT_VIG;
00436             pano.updateVariables(oldVars);
00437         };
00438     };
00439 
00440     // now take response curve into account
00441     vars |= OPT_RESP;
00442     // also WB if desired
00443     if (mode == OPT_PHOTOMETRIC_LDR_WB || mode == OPT_PHOTOMETRIC_HDR_WB)
00444     {
00445         vars |= OPT_WB;
00446     };
00447     VariableMapVector oldVars = pano.getVariables();
00448     optimizePhotometric(pano, 
00449                         createOptVars(pano, vars, opts.colorReferenceImage),
00450                         correspondences, imageStepSize, progress, error);
00451     // now check the results
00452     const bool hasHighVignetting = IsHighVignetting(pano.getImage(0).getRadialVigCorrCoeff());
00453     // @TODO check also response curve, what parameters are considered invalid?
00454     const bool hasHighWB = CheckStrangeWB(pano);
00455     if(hasHighVignetting)
00456     {
00457         vars &= ~OPT_VIG;
00458     };
00459     if(hasHighWB)
00460     {
00461         vars &= ~OPT_WB;
00462     };
00463     if(hasHighVignetting || hasHighWB)
00464     {
00465         // we got strange results, optimize again with less parameters
00466         pano.updateVariables(oldVars);
00467         if(vars>0)
00468         {
00469             optimizePhotometric(pano, 
00470                                 createOptVars(pano, vars, opts.colorReferenceImage),
00471                                 correspondences, imageStepSize, progress, error);
00472         };
00473     };
00474 }
00475 
00476 
00477 bool PhotometricOptimizer::runAlgorithm()
00478 {
00479     optimizePhotometric(o_panorama, 
00480                         o_vars, o_correspondences, o_imageStepSize,
00481                         getProgressDisplay(), o_resultError);
00482     
00483     // optimizePhotometric does not tell us if it's cancelled
00484     if (getProgressDisplay()->wasCancelled())
00485     {
00486         cancelAlgorithm();
00487     }
00488     
00489     return wasCancelled(); // let's hope so.
00490 }
00491 
00492 bool SmartPhotometricOptimizer::runAlgorithm()
00493 {
00494     smartOptimizePhotometric(o_panorama,
00495                              o_optMode,
00496                              o_correspondences, o_imageStepSize,
00497                              getProgressDisplay(), o_resultError);
00498     
00499     // smartOptimizePhotometric does not tell us if it's cancelled
00500     if (getProgressDisplay()->wasCancelled())
00501     {
00502         cancelAlgorithm();
00503     };
00504     
00505     return !wasCancelled(); // let's hope so.
00506 }
00507 
00508 } //namespace

Generated on 21 Feb 2018 for Hugintrunk by  doxygen 1.4.7