VigQuotientEstimator.h

Go to the documentation of this file.
00001 #ifndef _VIG_QUOTIENT_ESTIMATOR_H_
00002 #define _VIG_QUOTIENT_ESTIMATOR_H_
00003 
00004 #include <hugin_math/hugin_math.h>
00005 
00006 #include "RansacParameterEstimator.h"
00007 
00008 #include <iostream>
00009 
00010 #include <boost/random.hpp>
00011 
00012 #include <vigra/rgbvalue.hxx>
00013 #include <vigra/stdimage.hxx>
00014 
00015 #include <vigra_ext/ROIImage.h>
00016 #include <vigra_ext/utils.h>
00017 
00018 #include "ransac.h"
00019 
00020 // uncomment this to use our version of the linear solver
00021 //#define HUGIN_VIG_USE_UBLAS
00022 
00023 #ifdef HUGIN_VIG_USE_UBLAS
00024 
00025 #include <boost/numeric/ublas/vector.hpp>
00026 #include <boost/numeric/ublas/operation.hpp>
00027 #include <boost/numeric/ublas/vector_proxy.hpp>
00028 #include <boost/numeric/ublas/triangular.hpp>
00029 #include <boost/numeric/ublas/matrix.hpp>
00030 #include <boost/numeric/ublas/lu.hpp>
00031 #include <boost/numeric/ublas/io.hpp>
00032 
00033 #else
00034 
00035 #include <hugin_math/lu.h>
00036 
00037 #endif
00038 
00039 namespace vigra_ext
00040 {
00041 
00043 template <class PITER>
00044 double calcVigCorrPoly(PITER p, double r)
00045 {
00046     const unsigned n=3;
00047     double rsq = r*r;
00048     double res = p[n-1] * rsq; 
00049     for (int j=n-2; j >=0; j--) 
00050     {
00051        res = (res + p[j]) * rsq;
00052     }
00053     // add constant offset
00054     return res + 1;
00055 }
00056 
00057 
00058 // structure with various information about the fitting process
00059 struct VigQuotientEstimateResult
00060 {
00061     int nUsedPoints;
00062     double brightnessRMSE;
00063 };
00064 
00075 class VigQuotientEstimator : public RansacParameterEstimator
00076 {
00077 
00078 protected:
00079         double m_delta;
00080 
00081 public:
00082 
00083 #ifdef HUGIN_VIG_USE_UBLAS
00084     typedef boost::numeric::ublas::vector<double> Param;
00085 #else
00086     typedef std::vector<double> Param;
00087 //    typedef double[3] Param;
00088 #endif
00089 
00090     VigQuotientEstimator (double delta)
00091     : RansacParameterEstimator(3), m_delta(delta)
00092     {
00093     }
00094 
00109     virtual bool estimate(const std::vector<const PointPair *> &data,
00110                           Param &p) const
00111     {
00112         const size_t n=3;
00113         if(data.size() < n) {
00114             return false;
00115         }
00116         // estimate point
00117 
00118 
00119 #ifdef HUGIN_VIG_USE_UBLAS
00120         p.resize(n);
00121         // create linear equation 
00122         boost::numeric::ublas::permutation_matrix<size_t> permut(n);
00123         boost::numeric::ublas::matrix<double> A(n,n);
00124 
00125         for (unsigned i = 0; i < n; i++) {
00126             double r1sq = data[i]->r1 * data[i]->r1;
00127             double r2sq = data[i]->r2 * data[i]->r2;
00128             double r1f = r1sq;
00129             double r2f = r2sq;
00130             for (unsigned j=0;j < n; j++) {
00131                 A(i,j) = r2f - r1f;
00132                 r1f*=r1sq;
00133                 r2f*=r2sq;
00134             }
00135             double c = (data[i]->i1 / data[i]->i2);
00136             p[i] =  1 - c;
00137          }
00138 
00139          try {
00140 //             std::cout << "before factorize A:" << A << std::endl << "p: " << p << std::endl;
00141              lu_factorize(A,permut);
00142 //             std::cout << "after factorize A:" << A << std::endl << "perm: " << permut << std::endl;
00143              lu_substitute(A,permut,p); 
00144 //             std::cout << "after substitute A:" << A << std::endl << "p: " << p << std::endl;
00145          } catch (boost::numeric::ublas::singular){
00146              return false;
00147          }
00148          return true;
00149 #else
00150          p.resize(n);
00151          double matrix[n*(n+1)];
00152          for (unsigned i = 0; i < n; i++) {
00153              double r1sq = data[i]->r1 * data[i]->r1;
00154              double r2sq = data[i]->r2 * data[i]->r2;
00155              double r1f = r1sq;
00156              double r2f = r2sq;
00157              for (unsigned j=0;j < n; j++) {
00158                  matrix[i + j*n] = r2f - r1f;
00159                  r1f*=r1sq;
00160                  r2f*=r2sq;
00161              }
00162              double c = (data[i]->i1 / data[i]->i2);
00163              matrix[i + n*n] =  1 - c;
00164          }
00165 
00166          double sol[3];
00167          bool ok =math_lu_solve(matrix, sol, n) != 0;
00168          p[0] = sol[0];
00169          p[1] = sol[1];
00170          p[2] = sol[2];
00171          return ok;
00172 #endif
00173     }
00174 
00175 
00179     bool leastSquaresEstimate(const std::vector<const PointPair *> &data,
00180                               Param & p) const
00181     {
00182         /*
00183         {
00184             std::ofstream of("hugin_vigquot_least_sq_points");
00185             for (std::vector<const PointPair *>::const_iterator it = data.begin(); it != data.end(); ++it) {
00186                 of << (*it)->r1 << " " << (*it)->i1 << " " << (*it)->r2 << " " << (*it)->i2 << std::endl;
00187             }
00188         }
00189         */
00190 
00191         const size_t n=3;
00192         if(data.size() < n) {
00193             return false;
00194         }
00195         
00196 #ifdef HUGIN_VIG_USE_UBLAS
00197         // calculate least squares, based on pseudoinverse.
00198         // A'A x = A' b.
00199         boost::numeric::ublas::permutation_matrix<size_t> permut(n);
00200         boost::numeric::ublas::matrix<double> AtA(n,n);
00201         for (unsigned i=0; i< n; i++){
00202             for (unsigned j=0; j< n; j++){
00203                 AtA(i,j) = 0;
00204             }
00205         }
00206         // is this the right way to clear a boost matrix?
00207 //        AtA = boost::numeric::ublas::zero_matrix<double>(n,n);
00208         p.resize(n);
00209         // is this the right way to clear a boost vector?
00210         for (unsigned i=0; i< n; i++){
00211             p[i] = 0;
00212         }
00213 
00214         for(unsigned k=0; k < data.size(); k++) {
00215             // calculate one row of A
00216             double Arow[n];
00217             double r1sq = data[k]->r1 * data[k]->r1;
00218             double r2sq = data[k]->r2 * data[k]->r2;
00219             double r1f = r1sq;
00220             double r2f = r2sq;
00221             for (unsigned j=0;j < n; j++) {
00222                 Arow[j] = r2f - r1f;
00223                 r1f*=r1sq;
00224                 r2f*=r2sq;
00225             }
00226             double c = (data[k]->i1 / data[k]->i2);
00227             double bRow =  1 - c;
00228 
00229             // add to pseudoinverse
00230             for( unsigned i=0; i<n; ++i)
00231             {
00232                 // calculate  Atb
00233                 p[i]+=Arow[i]*bRow;
00234                 for( unsigned j=0; j<n; ++j)
00235                 {
00236                     AtA(i,j)+=Arow[i]*Arow[j];
00237                 }
00238             }
00239         }
00240 
00241         // solve system
00242         try {
00243             lu_factorize(AtA,permut);
00244             lu_substitute(AtA,permut,p); 
00245 //            std::cout << "lms solution: AtA (after solve):" << AtA << std::endl << "p: " << p << std::endl;
00246         } catch (boost::numeric::ublas::singular){
00247              return false;
00248         }
00249         return true;
00250 #else
00251         hugin_utils::LMS_Solver solver(n);
00252         for(unsigned k=0; k < data.size(); k++) {
00253             // calculate one row of A
00254             double Arow[n];
00255             double r1sq = data[k]->r1 * data[k]->r1;
00256             double r2sq = data[k]->r2 * data[k]->r2;
00257             double r1f = r1sq;
00258             double r2f = r2sq;
00259             for (unsigned j=0;j < n; j++) {
00260                 Arow[j] = r2f - r1f;
00261                 r1f*=r1sq;
00262                 r2f*=r2sq;
00263             }
00264             double c = (data[k]->i1 / data[k]->i2);
00265             double bRow =  1 - c;
00266             solver.addRow(&Arow[0], bRow);
00267         }
00268         return solver.solve(p);
00269 #endif
00270     }
00271 
00278     bool agree(Param &p, const PointPair &data) const
00279     {
00280         const int n = 3;
00281 
00282         double r1sq = data.r1 * data.r1;
00283         double r2sq = data.r2 * data.r2;
00284         double poly1 = p[n-1] * r1sq; 
00285         double poly2 = p[n-1] * r2sq; 
00286         for (int j=n-2; j >=0; j--) 
00287         {
00288             poly1 = (poly1 + p[j]) * r1sq;
00289             poly2 = (poly2 + p[j]) * r2sq;
00290         }
00291         poly1 += 1;
00292         poly2 += 1;
00293         double e;
00294         if (data.r1 < data.r2) {
00295             // I1 should be bigger than I2
00296             e = data.i2 / data.i1 - poly2/poly1;
00297         } else {
00298             e = data.i1 / data.i2 - poly1/poly2;
00299         }
00300         return fabs(e) < m_delta;
00301     }
00302 };
00303 
00304 
00305 
00306 template<class ImageType, class CoordType>
00307 void extractRandomPoints(std::vector<vigra_ext::ROIImage<ImageType, vigra::BImage> *>  &remapped,
00308                          std::vector<CoordType> & imgXCoord,
00309                          std::vector<CoordType> & imgYCoord,
00310                          const std::vector<vigra::Size2D> & imgSize,
00311                          const std::vector<hugin_utils::FDiff2D> & imgCenter,
00312                          unsigned nPointsPerOverlap,
00313                          std::vector<PointPair> & points,
00314                          unsigned & nBadPoints)
00315 {
00316     assert(remapped.size() == imgXCoord.size());
00317     assert(remapped.size() == imgYCoord.size());
00318 
00319     // init random number generator
00320     boost::mt19937 rng;
00321     unsigned nImg = remapped.size();
00322     nBadPoints = 0;
00323     // extract points from all overlaps.
00324     for (unsigned i=0; i < nImg; i++) {
00325         for (unsigned j=i+1; j < nImg; j++) {
00326             // get unifiying region.
00327             vigra::Rect2D roi1 = remapped[i]->boundingBox();
00328             vigra::Rect2D roi2 = remapped[j]->boundingBox();
00329             vigra::Rect2D roi = roi1 & roi2;
00330             if (roi.isEmpty()) {
00331                 continue;
00332             }
00333 
00334             double maxr1 = sqrt(((double)imgSize[i].x)*imgSize[i].x + ((double)imgSize[i].y)*imgSize[i].y) / 2.0;
00335             double maxr2 = sqrt(((double)imgSize[i].x)*imgSize[i].x + ((double)imgSize[i].y)*imgSize[i].y) / 2.0;
00336             // randomly sample points.
00337             boost::uniform_int<> distribx(roi.left(), roi.right()-1);
00338             boost::uniform_int<> distriby(roi.top(), roi.bottom()-1);
00339 
00340             boost::variate_generator<boost::mt19937&, boost::uniform_int<> >
00341                     randX(rng, distribx);             // glues randomness with mapping
00342             boost::variate_generator<boost::mt19937&, boost::uniform_int<> >
00343                     randY(rng, distriby);             // glues randomness with mapping
00344 
00345             for (unsigned maxTry = nPointsPerOverlap*2; nPointsPerOverlap > 0 && maxTry > 0; maxTry--) {
00346                 unsigned x = randX();
00347                 unsigned y = randY();
00348                 if (remapped[i]->getMask(x,y) && remapped[j]->getMask(x,y)) {
00349                     // extract gray value pair..
00350                     hugin_utils::FDiff2D p1(imgXCoord[i](x-roi1.left(), y-roi1.top()),
00351                                imgYCoord[i](x-roi1.left(), y-roi1.top()));
00352                     hugin_utils::FDiff2D p1c = p1 - imgCenter[i];
00353                     p1 = p1/maxr1;
00354                     p1c = p1c/maxr1;
00355                     hugin_utils::FDiff2D p2(imgXCoord[j](x-roi2.left(), y-roi2.top()),
00356                                imgYCoord[j](x-roi2.left(), y-roi2.top()));
00357                     hugin_utils::FDiff2D p2c = p2 - imgCenter[j];
00358                     p2 = p2/maxr2;
00359                     p2c = p2c/maxr2;
00360                     double r1 = sqrt(p1c.x*p1c.x + p1c.y*p1c.y);
00361                     double r2 = sqrt(p2c.x*p2c.x + p2c.y*p2c.y);
00362                     double i1 = remapped[i]->operator()(x,y);
00363                     double i2 = remapped[j]->operator()(x,y);
00364                     if (i1 >= i2 && r1 <= r2) {
00365                         // ok, point is good. i1 is closer to centre, swap point
00366                         // so that i1 < i2
00367                         points.push_back(PointPair(j, i2, p2, r2,   i, i1, p1, r1) );
00368                     } else if(i1 < i2 && r1 > r2) {
00369                         points.push_back(PointPair(i, i1, p1, r1,   j, i2, p2, r2) );
00370                     } else {
00371                         nBadPoints++;
00372                     }
00373 
00374                     nPointsPerOverlap--;
00375                 }
00376             }
00377         }
00378     }
00379     DEBUG_DEBUG("random point extraction: consistent: " << points.size()
00380                 << " inconsistent: " << nBadPoints << " points");
00381 }
00382 
00383 inline VigQuotientEstimateResult
00384 optimizeVignettingQuotient(const std::vector<PointPair> & points,
00385                            double ransacDelta,
00386                            std::vector<double> & vigCoeff)
00387 {
00388     VigQuotientEstimateResult res;
00389 
00390     // do optimisation..
00391     // the quotient error should be < 0.1
00392     VigQuotientEstimator vqEst(ransacDelta);
00393 
00394 #ifdef HUGIN_VIG_USE_UBLAS
00395     boost::numeric::ublas::vector<double> vigCoeff2(3);
00396     for (int i=0; i < 3; i++) vigCoeff2[i] = vigCoeff[i];
00397 
00398     res.nUsedPoints = Ransac::compute(vigCoeff2, vqEst, points, 0.99, 0.3);
00399 
00400     for (int i=0; i < 3; i++) vigCoeff[i] = vigCoeff2[i];
00401 #else
00402     std::vector<int> inliers_idx;
00403     std::vector<const PointPair *> inliers = Ransac::compute(vigCoeff, inliers_idx, vqEst, points, 0.99, 0.3);
00404     res.nUsedPoints = inliers.size();
00405 #endif
00406 
00407     double sqerror=0;
00408     // calculate brightness RMSE
00409     // not a real residual, since the brightness error was not minimized. instead
00410     // the quotient error was minimized. However, a brightness difference is observeable
00411     // in the final panorama and can be intuitively understood
00412     // (well, actually its the RMSE, but thats a detail ;)
00413     for (std::vector<PointPair>::const_iterator pnt=points.begin(); pnt != points.end();
00414          pnt++)
00415     {
00416 #ifdef HUGIN_VIG_USE_UBLAS
00417         sqerror += hugin_utils::sqr(calcVigCorrPoly(vigCoeff2, pnt->r1)*pnt->i1 - calcVigCorrPoly(vigCoeff, pnt->r2)*pnt->i2);
00418 #else
00419         sqerror += hugin_utils::sqr(calcVigCorrPoly(vigCoeff, pnt->r1)*pnt->i1 - calcVigCorrPoly(vigCoeff, pnt->r2)*pnt->i2);
00420 #endif
00421     }
00422     res.brightnessRMSE = sqrt(sqerror/points.size());
00423     return res;
00424 }
00425 
00426 #if 0
00427 template <class Image, class ImageSource>
00428 void runPointExtractor(Panorama & pano, ImageSource & imgLoader)
00429 {
00430     std::vector<FImage *> grayImgs;
00431     std::vector<BImage *> maskImgs;
00432     std::vector<FImage *> lapImgs;
00433     std::vector<InterpolImg> srcImgs;
00434     std::vector<SrcPanoImage> srcDescr;
00435 
00436     for (unsigned i=0; i < m_pano.getNrOfImages(); i++) {
00437         SrcPanoImage src = m_pano.getSrcImage(i);
00438         RGBFImage * img = imageLoader(src);
00439         ImageCache::Entry * cacheEntry =  ImageCache::getInstance().getSmallImage(src.getFilename().c_str());
00440             if (!cacheEntry) {
00441                 throw std::runtime_error("could not retrieve small source image for vignetting optimisation");
00442             }
00443             wxImage * srcImgWX = cacheEntry->image;
00444             // image view
00445             BasicImageView<RGBValue<unsigned char> > srcImgInCache((RGBValue<unsigned char> *)srcImgWX->GetData(),
00446                                      srcImgWX->GetWidth(),
00447                                      srcImgWX->GetHeight());
00448 
00449             src.resize(srcImgInCache.size());
00450             BImage * grayImg = new BImage(srcImgInCache.size());
00451             // change to grayscale
00452             // just use green channel
00453             vigra::copyImage(vigra::srcImageRange(srcImgInCache, vigra::GreenAccessor<RGBValue<unsigned char> >()),
00454                              vigra::destImage(*grayImg));
00455             BImage *maskImg = new BImage(srcImgInCache.size().x,srcImgInCache.size().y , 255);
00456 
00457             // TODO need to do gamma correction here!
00458             if (src.getGamma() != 1.0) {
00459                 vigra_ext::applyGammaCorrection(srcImageRange(*grayImg), destImage(*grayImg),
00460                                                 src.getGamma(), 255);
00461             }
00462             // TODO: cut out the masked parts.
00463 //            applyMask(src, vigra::destImageRange(*maskImg));
00464 
00465             vigra_ext::interp_spline36 interp;
00466             // create interpolating accessor to images.
00467             srcImgs.push_back(InterpolImg(vigra::srcImageRange(*grayImg),
00468                               vigra::srcImage(*maskImg), interp,
00469                               false) );
00470             grayImgs.push_back(grayImg);
00471             maskImgs.push_back(maskImg);
00472             srcDescr.push_back(src);
00473             if (uniformRadiusSample) {
00474                 vigra::FImage * lap = new vigra::FImage(grayImg->size());
00475                 laplacianOfGaussian(srcImageRange(*grayImg), destImage(*lap), 1);
00476                 lapImgs.push_back(lap);
00477             }
00478 
00479         }
00480 
00481 #endif
00482 
00483 }
00484 
00485 
00486 
00487 #endif //_VIG_QUOTIENT_PARAM_ESTIMATOR_H_

Generated on 1 Nov 2014 for Hugintrunk by  doxygen 1.4.7