PointSampler.h

Go to the documentation of this file.
00001 // -*- c-basic-offset: 4 -*-
00024 #ifndef _POINTSAMPLER_H
00025 #define _POINTSAMPLER_H
00026 
00027 #include <ctime>
00028 
00029 #include <hugin_shared.h>
00030 #include <hugin_config.h>
00031 #include <vigra/stdimage.hxx>
00032 #include <algorithms/PanoramaAlgorithm.h>
00033 
00034 #include <random>
00035 #include <functional>
00036 #include <vigra_ext/utils.h>
00037 #include <appbase/ProgressDisplay.h>
00038 #include <panodata/PanoramaData.h>
00039 
00040 namespace HuginBase
00041 {
00044     class IMPEX LimitIntensity
00045     {
00046         public:
00048             enum LimitType{LIMIT_UINT8, LIMIT_UINT16, LIMIT_FLOAT};
00050             LimitIntensity();
00052             LimitIntensity(LimitType limit);
00054             const float GetMinI() const {
00055                 return m_minI;
00056             };
00058             const float GetMaxI() const {
00059                 return m_maxI;
00060             }
00061         private:
00063             float m_minI, m_maxI;
00064     };
00065     typedef std::vector<LimitIntensity> LimitIntensityVector;
00066 
00067     class IMPEX PointSampler : public TimeConsumingPanoramaAlgorithm
00068     {
00069         protected:
00071             PointSampler(PanoramaData& panorama, AppBase::ProgressDisplay* progressDisplay,
00072                          std::vector<vigra::FRGBImage*> images, LimitIntensityVector limits,
00073                          int nPoints)
00074                 : TimeConsumingPanoramaAlgorithm(panorama, progressDisplay),
00075                   o_images(images), o_numPoints(nPoints), m_limits(limits)
00076             {};
00077         
00078         public:        
00080             virtual ~PointSampler() {};
00081 
00082         protected:
00084             typedef vigra_ext::ImageInterpolator<vigra::FRGBImage::const_traverser,
00085                                                  vigra::FRGBImage::ConstAccessor,
00086                                                  vigra_ext::interp_cubic           > InterpolImg;
00087             
00089             void sampleAndExtractPoints(AppBase::ProgressDisplay* progress);
00090             
00092             virtual void samplePoints(const std::vector<InterpolImg>& imgs,
00093                                       const std::vector<vigra::FImage*>& voteImgs,
00094                                       const PanoramaData& pano,
00095                                       const LimitIntensityVector limitI,
00096                                       std::vector<std::multimap<double,vigra_ext::PointPairRGB> >& radiusHist,
00097                                       unsigned& nGoodPoints,
00098                                       unsigned& nBadPoints,
00099                                       AppBase::ProgressDisplay*) = 0;
00100             
00105             template<class PointPairClass>
00106             static void sampleRadiusUniform(const std::vector<std::multimap<double,PointPairClass> >& radiusHist,
00107                                             unsigned nPoints,
00108                                             std::vector<PointPairClass>& selectedPoints,
00109                                             AppBase::ProgressDisplay*);
00110             
00111         public:
00113             virtual bool modifiesPanoramaData() const
00114                 { return false; }
00115             
00117             virtual bool runAlgorithm();
00118 
00119             PointSampler & execute()
00120             {
00121                                 run();
00122                 return *this;
00123             }
00124             
00125         public:
00127             typedef std::vector<vigra_ext::PointPairRGB> PointPairs;
00128             
00130             PointPairs getResultPoints()
00131             {
00132                 //[TODO] debug unsuccessful
00133                 
00134                 return o_resultPoints;
00135             }
00136             
00137             
00138         protected:
00139             std::vector<vigra::FRGBImage*> o_images;
00140             int o_numPoints;
00141             PointPairs o_resultPoints;
00142             LimitIntensityVector m_limits;
00143     };
00144 
00145 
00149     class AllPointSampler : public  PointSampler
00150     {
00151         public:
00153             AllPointSampler(PanoramaData& panorama, AppBase::ProgressDisplay* progressDisplay,
00154                                std::vector<vigra::FRGBImage*> images, LimitIntensityVector limits,
00155                                int nPoints)
00156              : PointSampler(panorama, progressDisplay, images, limits, nPoints)
00157             {};
00158 
00160             virtual ~AllPointSampler() {};
00161             
00162             
00163         public:
00167             template <class Img, class VoteImg, class PP>
00168             static void sampleAllPanoPoints(const std::vector<Img> &imgs,
00169                                      const std::vector<VoteImg *> &voteImgs,
00170                                      const PanoramaData& pano,
00171                                      int nPoints,
00172                                      const LimitIntensityVector limitI,
00173                                      std::vector<std::multimap<double, PP > > & radiusHist,
00174                                      unsigned & nGoodPoints,
00175                                      unsigned & nBadPoints,
00176                                      AppBase::ProgressDisplay* progress);
00177             
00178         protected:
00180             virtual void samplePoints(const std::vector<InterpolImg>& imgs,
00181                                       const std::vector<vigra::FImage*>& voteImgs,
00182                                       const PanoramaData& pano,
00183                                       const LimitIntensityVector limitI,
00184                                       std::vector<std::multimap<double,vigra_ext::PointPairRGB> >& radiusHist,
00185                                       unsigned& nGoodPoints,
00186                                       unsigned& nBadPoints,
00187                                       AppBase::ProgressDisplay* progress)
00188             {
00189                 sampleAllPanoPoints(imgs,
00190                                     voteImgs,
00191                                     pano,
00192                                     o_numPoints,
00193                                     limitI,
00194                                     radiusHist,
00195                                     nGoodPoints,
00196                                     nBadPoints,
00197                                     progress);
00198             }
00199             
00200     };
00201 
00202 
00206     class RandomPointSampler : public  PointSampler
00207     {
00208         public:
00210             RandomPointSampler(PanoramaData& panorama, AppBase::ProgressDisplay* progressDisplay,
00211                                std::vector<vigra::FRGBImage*> images, LimitIntensityVector limits,
00212                                int nPoints)
00213              : PointSampler(panorama, progressDisplay, images, limits, nPoints)
00214             {};
00215 
00217             virtual ~RandomPointSampler() {};
00218             
00219             
00220         public:
00221             template <class Img, class VoteImg, class PP>
00222             static void sampleRandomPanoPoints(const std::vector<Img>& imgs,
00223                                                const std::vector<VoteImg *> &voteImgs,
00224                                                const PanoramaData& pano,
00225                                                int nPoints,
00226                                                const LimitIntensityVector limitI,
00227                                                std::vector<std::multimap<double, PP > > & radiusHist,
00228                                                unsigned & nBadPoints,
00229                                                AppBase::ProgressDisplay* progress);
00230             
00231         protected:
00233             virtual void samplePoints(const std::vector<InterpolImg>& imgs,
00234                                       const std::vector<vigra::FImage*>& voteImgs,
00235                                       const PanoramaData& pano,
00236                                       const LimitIntensityVector limitI,
00237                                       std::vector<std::multimap<double,vigra_ext::PointPairRGB> >& radiusHist,
00238                                       unsigned& nGoodPoints,
00239                                       unsigned& nBadPoints,
00240                                       AppBase::ProgressDisplay* progress)
00241             {
00242                  sampleRandomPanoPoints(imgs,
00243                                         voteImgs,
00244                                         pano,
00245                                         5*o_numPoints,
00246                                         limitI,
00247                                         radiusHist,
00248                                         nBadPoints,
00249                                         progress);
00250             }
00251             
00252     };
00253 
00254     
00255     
00256 } // namespace
00257     
00258     
00259 
00260 
00261     
00262 //==============================================================================
00263 //  templated methods
00264 
00265 #include <panotools/PanoToolsInterface.h>
00266 
00267 namespace HuginBase {
00268 
00269 template<class PointPairClass>
00270 void PointSampler::sampleRadiusUniform(const std::vector<std::multimap<double, PointPairClass> >& radiusHist,
00271                                        unsigned nPoints,
00272                                        std::vector<PointPairClass> &selectedPoints,
00273                                        AppBase::ProgressDisplay* progress)
00274 {
00275     typedef std::multimap<double,PointPairClass> PointPairMap;
00276     
00277     // reserve selected points..
00278     int drawsPerBin = nPoints / radiusHist.size();
00279     // reallocate output vector.
00280     selectedPoints.reserve(drawsPerBin*radiusHist.size());
00281     for (typename std::vector<PointPairMap>::const_iterator bin = radiusHist.begin();
00282          bin != radiusHist.end(); ++bin) 
00283     {
00284         unsigned i=drawsPerBin;
00285         // copy into vector
00286         for (typename PointPairMap::const_iterator it= (*bin).begin();
00287              it != (*bin).end(); ++it) 
00288         {
00289             selectedPoints.push_back(it->second);
00290             // do not extract more than drawsPerBin Point pairs.
00291             --i;
00292             if (i == 0)
00293                 break;
00294         }
00295     }
00296 }
00297     
00298     
00299     
00300 template <class Img, class VoteImg, class PP>
00301 void AllPointSampler::sampleAllPanoPoints(const std::vector<Img> &imgs,
00302                                           const std::vector<VoteImg *> &voteImgs,
00303                                           const PanoramaData& pano,
00304                                           int nPoints,
00305                                           const LimitIntensityVector limitI,
00306                                           //std::vector<vigra_ext::PointPair> &points,
00307                                           std::vector<std::multimap<double, PP > > & radiusHist,
00308                                           unsigned & nGoodPoints,
00309                                           unsigned & nBadPoints,
00310                                           AppBase::ProgressDisplay* progress)
00311 {
00312     typedef typename Img::PixelType PixelType;
00313 
00314     // use 10 bins
00315     radiusHist.resize(10);
00316     const unsigned pairsPerBin = nPoints / radiusHist.size();
00317 
00318     nGoodPoints = 0;
00319     nBadPoints = 0;
00320     vigra_precondition(imgs.size() > 1, "sampleAllPanoPoints: At least two images required");
00321     
00322     const unsigned nImg = imgs.size();
00323 
00324     vigra::Size2D srcSize = pano.getImage(0).getSize();
00325     double maxr = sqrt(((double)srcSize.x)*srcSize.x + ((double)srcSize.y)*srcSize.y) / 2.0;
00326 
00327     // create an array of transforms.
00328     //std::vector<SpaceTransform> transf(imgs.size());
00329     std::vector<PTools::Transform*> transf(imgs.size());
00330 
00331     // initialize transforms, and interpolating accessors
00332     for(unsigned i=0; i < nImg; i++) {
00333         vigra_precondition(pano.getImage(i).getSize() == srcSize, "images need to have the same size");
00334         transf[i] = new PTools::Transform;
00335         transf[i]->createTransform(pano.getImage(i), pano.getOptions());
00336     }
00337 
00338     const vigra::Rect2D roi = pano.getOptions().getROI();
00339     for (int y=roi.top(); y < roi.bottom(); ++y) {
00340         for (int x=roi.left(); x < roi.right(); ++x) {
00341             hugin_utils::FDiff2D panoPnt(x,y);
00342             for (unsigned i=0; i< nImg-1; i++) {
00343             // transform pixel
00344                 hugin_utils::FDiff2D p1;
00345                 if(!transf[i]->transformImgCoord(p1, panoPnt))
00346                     continue;
00347                 vigra::Point2D p1Int(p1.toDiff2D());
00348                 // is inside:
00349                 if (!pano.getImage(i).isInside(p1Int)) {
00350                     // point is outside image
00351                     continue;
00352                 }
00353                 PixelType i1;
00354                 vigra::UInt8 maskI;
00355                 if (imgs[i](p1.x,p1.y, i1, maskI)){
00356                     float im1 = vigra_ext::getMaxComponent(i1);
00357                     if (limitI[i].GetMinI() > im1 || limitI[i].GetMaxI() < im1 || maskI == 0) {
00358                         // ignore pixels that are too dark or bright
00359                         continue;
00360                     }
00361                     double r1 = hugin_utils::norm((p1 - pano.getImage(i).getRadialVigCorrCenter()) / maxr);
00362 
00363                     // check inner image
00364                     for (unsigned j=i+1; j < nImg; j++) {
00365                         hugin_utils::FDiff2D p2;
00366                         if(!transf[j]->transformImgCoord(p2, panoPnt))
00367                             continue;
00368                         vigra::Point2D p2Int(p2.toDiff2D());
00369                         if (!pano.getImage(j).isInside(p2Int)) {
00370                             // point is outside image
00371                             continue;
00372                         }
00373                         PixelType i2;
00374                         vigra::UInt8 maskI2;
00375                         if (imgs[j](p2.x, p2.y, i2, maskI2)){
00376                             float im2 = vigra_ext::getMaxComponent(i2);
00377                             if (limitI[j].GetMinI() > im2 || limitI[j].GetMaxI() < im2 || maskI2 == 0) {
00378                                 // ignore pixels that are too dark or bright
00379                                 continue;
00380                             }
00381                             double r2 = hugin_utils::norm((p2 - pano.getImage(j).getRadialVigCorrCenter()) / maxr);
00382                             // add pixel
00383                             const VoteImg & vimg1 =  *voteImgs[i];
00384                             const VoteImg & vimg2 =  *voteImgs[j];
00385                             double laplace = hugin_utils::sqr(vimg1[p1Int]) + hugin_utils::sqr(vimg2[p2Int]);
00386                             int bin1 = (int)(r1*10);
00387                             int bin2 = (int)(r2*10);
00388                             // a center shift might lead to radi > 1.
00389                             if (bin1 > 9) bin1 = 9;
00390                             if (bin2 > 9) bin2 = 9;
00391 
00392                             PP pp;
00393                             if (im1 <= im2) {
00394                                 // choose i1 to be smaller than i2
00395                                 pp = PP(i, i1, p1, r1,   j, i2, p2, r2);
00396                             } else {
00397                                 pp = PP(j, i2, p2, r2,   i, i1, p1, r1);
00398                             }
00399 
00400                             // decide which bin should be used.
00401                             std::multimap<double, PP> * map1 = &radiusHist[bin1];
00402                             std::multimap<double, PP> * map2 = &radiusHist[bin2];
00403                             std::multimap<double, PP> * destMap;
00404                             if (map1->empty()) {
00405                                 destMap = map1;
00406                             } else if (map2->empty()) {
00407                                 destMap = map2;
00408                             } else if (map1->size() < map2->size()) {
00409                                 destMap = map1;
00410                             } else if (map1->size() > map2->size()) {
00411                                 destMap = map2;
00412                             } else if (map1->rbegin()->first > map2->rbegin()->first) {
00413                                 // heuristic: insert into bin with higher maximum laplacian filter response
00414                                 // (higher probablity of misregistration).
00415                                 destMap = map1;
00416                             } else {
00417                                 destMap = map2;
00418                             }
00419                             // insert
00420                             destMap->insert(std::make_pair(laplace,pp));
00421                             // remove last element if too many elements have been gathered
00422                             if (destMap->size() > pairsPerBin) {
00423                                 destMap->erase((--(destMap->end())));
00424                             }
00425                             nGoodPoints++;
00426                         }
00427                     }
00428                 }
00429             }
00430         }
00431     }
00432 
00433     for(unsigned i=0; i < nImg; i++) {
00434         delete transf[i];
00435     }
00436 }
00437 
00438 
00439 
00440 template <class Img, class VoteImg, class PP>
00441 void RandomPointSampler::sampleRandomPanoPoints(const std::vector<Img>& imgs,
00442                                                 const std::vector<VoteImg *> &voteImgs,
00443                                                 const PanoramaData& pano,
00444                                                 int nPoints,
00445                                                 const LimitIntensityVector limitI,
00446                                                 //std::vector<PP> &points,
00447                                                 std::vector<std::multimap<double, PP > > & radiusHist,
00448                                                 unsigned & nBadPoints,
00449                                                 AppBase::ProgressDisplay* progress)
00450 {
00451     typedef typename Img::PixelType PixelType;
00452 
00453     vigra_precondition(imgs.size() > 1, "sampleRandomPanoPoints: At least two images required");
00454     
00455     const unsigned nImg = imgs.size();
00456 
00457     const unsigned nBins = radiusHist.size();
00458     const unsigned pairsPerBin = nPoints / nBins;
00459 
00460     // create an array of transforms.
00461     //std::vector<SpaceTransform> transf(imgs.size());
00462     std::vector<PTools::Transform *> transf(imgs.size());
00463     std::vector<double> maxr(imgs.size());
00464 
00465     // initialize transforms, and interpolating accessors
00466     for(unsigned i=0; i < nImg; i++) {
00467         // same size is not needed?
00468 //        vigra_precondition(src[i].getSize() == srcSize, "images need to have the same size");
00469         transf[i] = new PTools::Transform;
00470         transf[i]->createTransform(pano.getImage(i), pano.getOptions());
00471         vigra::Size2D srcSize = pano.getImage(i).getSize();
00472         maxr[i] = sqrt(((double)srcSize.x)*srcSize.x + ((double)srcSize.y)*srcSize.y) / 2.0;
00473     }
00474     // init random number generator
00475     const vigra::Rect2D roi = pano.getOptions().getROI();
00476     std::mt19937 rng(static_cast<unsigned int>(std::time(0)));
00477     std::uniform_int_distribution<unsigned int> distribx(roi.left(), roi.right()-1);
00478     std::uniform_int_distribution<unsigned int> distriby(roi.top(), roi.bottom()-1);
00479     auto randX = std::bind(distribx, std::ref(rng));
00480     auto randY = std::bind(distriby, std::ref(rng));
00481 
00482     for (unsigned maxTry = nPoints*5; nPoints > 0 && maxTry > 0; maxTry--) {
00483         unsigned x = randX();
00484         unsigned y = randY();
00485         hugin_utils::FDiff2D panoPnt(x,y);
00486         for (unsigned i=0; i< nImg-1; i++) {
00487             // transform pixel
00488             PixelType i1;
00489             hugin_utils::FDiff2D p1;
00490             if(!transf[i]->transformImgCoord(p1, panoPnt))
00491                 continue;
00492             vigra::Point2D p1Int(p1.toDiff2D());
00493             // check if pixel is valid
00494             if (!pano.getImage(i).isInside(p1Int)) {
00495                 // point is outside image
00496                 continue;
00497             }
00498             vigra::UInt8 maskI;
00499             if ( imgs[i](p1.x,p1.y, i1, maskI)){
00500                 float im1 = vigra_ext::getMaxComponent(i1);
00501                 if (limitI[i].GetMinI() > im1 || limitI[i].GetMaxI() < im1) {
00502                     // ignore pixels that are too dark or bright
00503                     continue;
00504                 }
00505                 double r1 = hugin_utils::norm((p1 - pano.getImage(i).getRadialVigCorrCenter()) / maxr[i]);
00506                 for (unsigned j=i+1; j < nImg; j++) {
00507                     PixelType i2;
00508                     hugin_utils::FDiff2D p2;
00509                     if(!transf[j]->transformImgCoord(p2, panoPnt))
00510                         continue;
00511                     // check if a pixel is inside the source image
00512                     vigra::Point2D p2Int(p2.toDiff2D());
00513                     if (!pano.getImage(j).isInside(p2Int)) {
00514                         // point is outside image
00515                         continue;
00516                     }
00517                     vigra::UInt8 maskI2;
00518                     if (imgs[j](p2.x, p2.y, i2, maskI2)){
00519                         float im2 = vigra_ext::getMaxComponent(i2);
00520                         if (limitI[j].GetMinI() > im2 || limitI[j].GetMaxI() < im2) {
00521                             // ignore pixels that are too dark or bright
00522                             continue;
00523                         }
00524                         // TODO: add check for gradient radius.
00525                         double r2 = hugin_utils::norm((p2 - pano.getImage(j).getRadialVigCorrCenter()) / maxr[j]);
00526 #if 0
00527                         // add pixel
00528                         if (im1 <= im2) {
00529                             points.push_back(PP(i, i1, p1, r1,   j, i2, p2, r2) );
00530                         } else {
00531                             points.push_back(PP(j, i2, p2, r2,   i, i1, p1, r1) );
00532                         }
00533 #else
00534                             // add pixel
00535                             const VoteImg & vimg1 =  *voteImgs[i];
00536                             const VoteImg & vimg2 =  *voteImgs[j];
00537                             double laplace = hugin_utils::sqr(vimg1[p1Int]) + hugin_utils::sqr(vimg2[p2Int]);
00538                             size_t bin1 = (size_t)(r1*nBins);
00539                             size_t bin2 = (size_t)(r2*nBins);
00540                             // a center shift might lead to radi > 1.
00541                             if (bin1+1 > nBins) bin1 = nBins-1;
00542                             if (bin2+1 > nBins) bin2 = nBins-1;
00543 
00544                             PP pp;
00545                             if (im1 <= im2) {
00546                                 // choose i1 to be smaller than i2
00547                                 pp = PP(i, i1, p1, r1,   j, i2, p2, r2);
00548                             } else {
00549                                 pp = PP(j, i2, p2, r2,   i, i1, p1, r1);
00550                             }
00551 
00552                             // decide which bin should be used.
00553                             std::multimap<double, PP> * map1 = &radiusHist[bin1];
00554                             std::multimap<double, PP> * map2 = &radiusHist[bin2];
00555                             std::multimap<double, PP> * destMap;
00556                             if (map1->empty()) {
00557                                 destMap = map1;
00558                             } else if (map2->empty()) {
00559                                 destMap = map2;
00560                             } else if (map1->size() < map2->size()) {
00561                                 destMap = map1;
00562                             } else if (map1->size() > map2->size()) {
00563                                 destMap = map2;
00564                                                         } else if (map1->rbegin()->first > map2->rbegin()->first) {
00565                                 // heuristic: insert into bin with higher maximum laplacian filter response
00566                                 // (higher probablity of misregistration).
00567                                 destMap = map1;
00568                             } else {
00569                                 destMap = map2;
00570                             }
00571                             // insert
00572                             destMap->insert(std::make_pair(laplace,pp));
00573                             // remove last element if too many elements have been gathered
00574                             if (destMap->size() > pairsPerBin) {
00575                                 destMap->erase((--(destMap->end())));
00576                             }
00577 //                            nGoodPoints++;
00578 #endif
00579                         nPoints--;
00580                     }
00581                 }
00582             }
00583         }
00584     }
00585     for(unsigned i=0; i < imgs.size(); i++) {
00586         delete transf[i];
00587     }
00588    
00589 }
00590 
00591 
00592 } // namespace
00593 #endif //_H

Generated on 4 Dec 2016 for Hugintrunk by  doxygen 1.4.7