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