PanoDetectorLogic.cpp

Go to the documentation of this file.
00001 // -*- c-basic-offset: 4 ; tab-width: 4 -*-
00002 /*
00003 * Copyright (C) 2007-2008 Anael Orlinski
00004 *
00005 * This file is part of Panomatic.
00006 *
00007 * Panomatic is free software; you can redistribute it and/or modify
00008 * it under the terms of the GNU General Public License as published by
00009 * the Free Software Foundation; either version 2 of the License, or
00010 * (at your option) any later version.
00011 *
00012 * Panomatic is distributed in the hope that it will be useful,
00013 * but WITHOUT ANY WARRANTY; without even the implied warranty of
00014 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015 * GNU General Public License for more details.
00016 *
00017 * You should have received a copy of the GNU General Public License
00018 * along with Panomatic; if not, write to the Free Software
00019 * <http://www.gnu.org/licenses/>.
00020 */
00021 
00022 #include "ImageImport.h"
00023 
00024 #include "PanoDetector.h"
00025 #include <iostream>
00026 #include <fstream>
00027 #include <vigra/distancetransform.hxx>
00028 #include "vigra_ext/impexalpha.hxx"
00029 
00030 #include <localfeatures/Sieve.h>
00031 #include <localfeatures/PointMatch.h>
00032 #include <localfeatures/RansacFiltering.h>
00033 #include <localfeatures/KeyPointIO.h>
00034 #include <localfeatures/CircularKeyPointDescriptor.h>
00035 
00036 /*
00037 #include "KDTree.h"
00038 #include "KDTreeImpl.h"
00039 */
00040 #include "Utils.h"
00041 #include "Tracer.h"
00042 
00043 #include <algorithms/nona/ComputeImageROI.h>
00044 #include <algorithms/optimizer/PTOptimizer.h>
00045 #include <nona/RemappedPanoImage.h>
00046 #include <nona/ImageRemapper.h>
00047 
00048 #include <time.h>
00049 
00050 #define TRACE_IMG(X) {if (iPanoDetector.getVerbose() > 1) { TRACE_INFO("i" << ioImgInfo._number << " : " << X << endl);} }
00051 #define TRACE_PAIR(X) {if (iPanoDetector.getVerbose() > 1){ TRACE_INFO("i" << ioMatchData._i1->_number << " <> " \
00052                 "i" << ioMatchData._i2->_number << " : " << X << endl)}}
00053 
00054 using namespace std;
00055 using namespace lfeat;
00056 using namespace HuginBase;
00057 using namespace AppBase;
00058 using namespace HuginBase::Nona;
00059 using namespace hugin_utils;
00060 
00061 
00062 // define a Keypoint insertor
00063 class KeyPointVectInsertor : public lfeat::KeyPointInsertor
00064 {
00065 public:
00066     explicit KeyPointVectInsertor(KeyPointVect_t& iVect) : _v(iVect) {};
00067     inline virtual void operator()(const lfeat::KeyPoint& k)
00068     {
00069         _v.push_back(KeyPointPtr(new lfeat::KeyPoint(k)));
00070     }
00071 
00072 private:
00073     KeyPointVect_t& _v;
00074 
00075 };
00076 
00077 
00078 // define a sieve extractor
00079 class SieveExtractorKP : public lfeat::SieveExtractor<KeyPointPtr>
00080 {
00081 public:
00082     explicit SieveExtractorKP(KeyPointVect_t& iV) : _v(iV) {};
00083     inline virtual void operator()(const KeyPointPtr& k)
00084     {
00085         _v.push_back(k);
00086     }
00087 private:
00088     KeyPointVect_t& _v;
00089 };
00090 
00091 class SieveExtractorMatch : public lfeat::SieveExtractor<lfeat::PointMatchPtr>
00092 {
00093 public:
00094     explicit SieveExtractorMatch(lfeat::PointMatchVector_t& iM) : _m(iM) {};
00095     inline virtual void operator()(const lfeat::PointMatchPtr& m)
00096     {
00097         _m.push_back(m);
00098     }
00099 private:
00100     lfeat::PointMatchVector_t& _m;
00101 };
00102 
00103 bool PanoDetector::LoadKeypoints(ImgData& ioImgInfo, const PanoDetector& iPanoDetector)
00104 {
00105     TRACE_IMG("Loading keypoints...");
00106 
00107     ImageInfo info = lfeat::loadKeypoints(ioImgInfo._keyfilename, ioImgInfo._kp);
00108     ioImgInfo._loadFail = (info.filename.size() == 0);
00109 
00110     // update ImgData
00111     if(ioImgInfo._needsremap)
00112     {
00113         ioImgInfo._detectWidth = max(info.width,info.height);
00114         ioImgInfo._detectHeight = max(info.width,info.height);
00115         ioImgInfo._projOpts.setWidth(ioImgInfo._detectWidth);
00116         ioImgInfo._projOpts.setHeight(ioImgInfo._detectHeight);
00117     }
00118     else
00119     {
00120         ioImgInfo._detectWidth = info.width;
00121         ioImgInfo._detectHeight = info.height;
00122     };
00123     ioImgInfo._descLength = info.dimensions;
00124 
00125     return true;
00126 }
00127 
00128 vigra::RGBValue<float> gray2RGB(float const& v)
00129 {
00130     return vigra::RGBValue<float>(v,v,v);
00131 }
00132 
00133 template <class SrcImageIterator, class SrcAccessor>
00134 void applyMaskAndCrop(vigra::triple<SrcImageIterator, SrcImageIterator, SrcAccessor> img, const HuginBase::SrcPanoImage& SrcImg)
00135 {
00136     vigra::Diff2D imgSize = img.second - img.first;
00137 
00138     // create dest y iterator
00139     SrcImageIterator yd(img.first);
00140     // loop over the image and transform
00141     for(int y=0; y < imgSize.y; ++y, ++yd.y)
00142     {
00143         // create x iterators
00144         SrcImageIterator xd(yd);
00145         for(int x=0; x < imgSize.x; ++x, ++xd.x)
00146         {
00147             if(!SrcImg.isInside(vigra::Point2D(x,y)))
00148             {
00149                 *xd=0;
00150             };
00151         }
00152     }
00153 }
00154 
00155 // save some intermediate images to disc if defined
00156 //#define DEBUG_LOADING_REMAPPING
00157 bool PanoDetector::AnalyzeImage(ImgData& ioImgInfo, const PanoDetector& iPanoDetector)
00158 {
00159     vigra::DImage final_img;
00160     vigra::BImage final_mask;
00161 
00162     try
00163     {
00164         ioImgInfo._loadFail=false;
00165 
00166         TRACE_IMG("Load image...");
00167         vigra::ImageImportInfo aImageInfo(ioImgInfo._name.c_str());
00168         vigra::FRGBImage RGBimg(aImageInfo.width(), aImageInfo.height());
00169         vigra::BImage mask;
00170 
00171         if (aImageInfo.numBands() == 1)
00172         {
00173             // grayscale image, convert to RGB. This is kind of stupid, but celeste wants RGB images...
00174             vigra::FImage tmpImg(aImageInfo.width(), aImageInfo.height());
00175             if (aImageInfo.numExtraBands() == 0)
00176             {
00177                 vigra::importImage(aImageInfo, destImage(tmpImg));
00178             }
00179             else if (aImageInfo.numExtraBands() == 1)
00180             {
00181                 mask.resize(aImageInfo.size());
00182                 importImageAlpha(aImageInfo, destImage(tmpImg), destImage(mask));
00183             }
00184             else
00185             {
00186                 TRACE_INFO("Image with multiple alpha channels are not supported");
00187                 ioImgInfo._loadFail = true;
00188                 return false;
00189             }
00190             //vigra::GrayToRGBAccessor<vigra::RGBValue<float> > ga;
00191             RGBimg.resize(aImageInfo.size());
00192             vigra::transformImage(srcImageRange(tmpImg), destImage(RGBimg), &gray2RGB);
00193         }
00194         else
00195         {
00196             if(aImageInfo.numExtraBands() == 1)
00197             {
00198                 mask.resize(aImageInfo.size());
00199                 importImageAlpha(aImageInfo, destImage(RGBimg), destImage(mask));
00200             }
00201             else
00202             {
00203                 if (aImageInfo.numExtraBands() == 0)
00204                 {
00205                     vigra::importImage(aImageInfo, destImage(RGBimg));
00206                 }
00207                 else
00208                 {
00209                     TRACE_INFO("Image with multiple alpha channels are not supported");
00210                     ioImgInfo._loadFail = true;
00211                     return false;
00212                 };
00213             };
00214         }
00215 
00216         //convert image, so that all (rgb) values are between 0 and 1
00217         if(aImageInfo.getPixelType() == std::string("FLOAT") || aImageInfo.getPixelType() == std::string("DOUBLE"))
00218         {
00219             vigra::RGBToGrayAccessor<vigra::RGBValue<float> > ga;
00220             vigra::FindMinMax<float> minmax;   // init functor
00221             vigra::inspectImage(srcImageRange(RGBimg, ga),minmax);
00222             double minVal = minmax.min;
00223             double maxVal = minmax.max;
00224             vigra_ext::applyMapping(srcImageRange(RGBimg), destImage(RGBimg), minVal, maxVal, 0);
00225         }
00226         else
00227         {
00228             vigra::transformImage(srcImageRange(RGBimg), destImage(RGBimg),
00229                                   vigra::functor::Arg1()/vigra::functor::Param(vigra_ext::getMaxValForPixelType(aImageInfo.getPixelType())));
00230         };
00231 
00232         if(ioImgInfo._needsremap)
00233         {
00234             TRACE_IMG("Remap image...");
00235             RemappedPanoImage<vigra::FRGBImage,vigra::BImage>* remapped=new RemappedPanoImage<vigra::FRGBImage,vigra::BImage>;
00236             vigra::FImage ffImg;
00237             ProgressDisplay* progress=new DummyProgressDisplay();
00238             remapped->setPanoImage(iPanoDetector._panoramaInfoCopy.getImage(ioImgInfo._number),
00239                                    ioImgInfo._projOpts, ioImgInfo._projOpts.getROI());
00240             if(mask.width()>0)
00241             {
00242                 remapped->remapImage(vigra::srcImageRange(RGBimg), vigra::srcImage(mask), vigra_ext::INTERP_CUBIC, progress, true);
00243             }
00244             else
00245             {
00246                 remapped->remapImage(vigra::srcImageRange(RGBimg), vigra_ext::INTERP_CUBIC, progress, true);
00247             };
00248             RGBimg.resize(0,0);
00249             mask.resize(0,0);
00250             RGBimg=remapped->m_image;
00251             mask=remapped->m_mask;
00252             delete remapped;
00253             delete progress;
00254         }
00255         else
00256         {
00257             const SrcPanoImage& SrcImg=iPanoDetector._panoramaInfoCopy.getImage(ioImgInfo._number);
00258             if(SrcImg.hasActiveMasks() || (SrcImg.getCropMode()!=SrcPanoImage::NO_CROP && !SrcImg.getCropRect().isEmpty()))
00259             {
00260                 if(mask.width()!=aImageInfo.width() || mask.height()!=aImageInfo.height())
00261                 {
00262                     mask.resize(aImageInfo.size().width(),aImageInfo.size().height(),255);
00263                 };
00264                 //copy mask and crop from pto file into alpha layer
00265                 applyMaskAndCrop(vigra::destImageRange(mask), SrcImg);
00266             };
00267         };
00268 
00269         if(iPanoDetector.getCeleste())
00270         {
00271             vigra::FRGBImage scaled(ioImgInfo._detectWidth,ioImgInfo._detectHeight);
00272             if(iPanoDetector._downscale && (!ioImgInfo._needsremap))
00273             {
00274                 TRACE_IMG("Resize image...");
00275                 vigra::resizeImageNoInterpolation(srcImageRange(RGBimg),destImageRange(scaled));
00276                 if(mask.width()>0 && mask.height()>0)
00277                 {
00278                     final_mask.resize(ioImgInfo._detectWidth, ioImgInfo._detectHeight);
00279                     vigra::resizeImageNoInterpolation(srcImageRange(mask),destImageRange(final_mask));
00280                 };
00281             }
00282             else
00283             {
00284                 vigra::copyImage(srcImageRange(RGBimg),destImage(scaled));
00285                 if(mask.width()>0 && mask.height()>0)
00286                 {
00287                     final_mask.resize(ioImgInfo._detectWidth, ioImgInfo._detectHeight);
00288                     vigra::copyImage(srcImageRange(mask),destImage(final_mask));
00289                 };
00290             };
00291             RGBimg.resize(0,0);
00292             mask.resize(0,0);
00293             TRACE_IMG("Mask areas with clouds...");
00294             vigra::UInt16RGBImage image16(scaled.size());
00295             vigra::transformImage(srcImageRange(scaled),destImage(image16),vigra::functor::Arg1()*vigra::functor::Param(65535));
00296             int radius=iPanoDetector.getCelesteRadius();
00297             if(iPanoDetector._downscale)
00298             {
00299                 radius>>= 1;
00300             };
00301             if(radius<2)
00302             {
00303                 radius=2;
00304             };
00305             vigra::BImage celeste_mask=celeste::getCelesteMask(iPanoDetector.svmModel,image16,radius,iPanoDetector.getCelesteThreshold(),800,true,false);
00306 #ifdef DEBUG_LOADING_REMAPPING
00307             // DEBUG: export celeste mask
00308             std::ostringstream maskfilename;
00309             maskfilename << ioImgInfo._name << "_celeste_mask.JPG";
00310             vigra::ImageExportInfo maskexinfo(maskfilename.str().c_str());
00311             vigra::exportImage(srcImageRange(celeste_mask), maskexinfo);
00312 #endif
00313             image16.resize(0,0);
00314             if(final_mask.width()>0)
00315             {
00316                 vigra::copyImageIf(srcImageRange(celeste_mask),srcImage(final_mask),destImage(final_mask));
00317             }
00318             else
00319             {
00320                 final_mask.resize(ioImgInfo._detectWidth,ioImgInfo._detectHeight);
00321                 vigra::copyImage(srcImageRange(celeste_mask),destImage(final_mask));
00322             };
00323             celeste_mask.resize(0,0);
00324             TRACE_IMG("Convert to greyscale double...");
00325             final_img.resize(ioImgInfo._detectWidth,ioImgInfo._detectHeight);
00326             vigra::copyImage(scaled.upperLeft(), scaled.lowerRight(), vigra::RGBToGrayAccessor<vigra::RGBValue<double> >(),
00327                              final_img.upperLeft(), vigra::DImage::Accessor());
00328             scaled.resize(0,0);
00329         }
00330         else
00331         {
00332             //without celeste
00333             final_img.resize(ioImgInfo._detectWidth,ioImgInfo._detectHeight);
00334             if (iPanoDetector._downscale && !ioImgInfo._needsremap)
00335             {
00336                 // Downscale and convert to grayscale double format
00337                 TRACE_IMG("Resize to greyscale double...");
00338                 vigra::resizeImageNoInterpolation(
00339                     RGBimg.upperLeft(),
00340                     RGBimg.upperLeft() + vigra::Diff2D(aImageInfo.width(), aImageInfo.height()),
00341                     vigra::RGBToGrayAccessor<vigra::RGBValue<double> >(),
00342                     final_img.upperLeft(),
00343                     final_img.lowerRight(),
00344                     vigra::DImage::Accessor());
00345                 RGBimg.resize(0,0);
00346                 //downscale mask
00347                 if(mask.width()>0 && mask.height()>0)
00348                 {
00349                     final_mask.resize(ioImgInfo._detectWidth, ioImgInfo._detectHeight);
00350                     vigra::resizeImageNoInterpolation(srcImageRange(mask),destImageRange(final_mask));
00351                     mask.resize(0,0);
00352                 };
00353             }
00354             else
00355             {
00356                 // convert to grayscale
00357                 TRACE_IMG("Convert to greyscale double...");
00358                 vigra::copyImage(RGBimg.upperLeft(), RGBimg.lowerRight(), vigra::RGBToGrayAccessor<vigra::RGBValue<double> >(),
00359                                  final_img.upperLeft(), vigra::DImage::Accessor());
00360                 RGBimg.resize(0,0);
00361                 if(mask.width()>0 && mask.height()>0)
00362                 {
00363                     final_mask.resize(ioImgInfo._detectWidth, ioImgInfo._detectHeight);
00364                     vigra::copyImage(srcImageRange(mask),destImage(final_mask));
00365                     mask.resize(0,0);
00366                 };
00367             };
00368         };
00369 
00370         //now scale image from 0..1 to 0..255
00371         vigra::transformImage(srcImageRange(final_img),destImage(final_img),vigra::functor::Arg1()*vigra::functor::Param(255));
00372 
00373 #ifdef DEBUG_LOADING_REMAPPING
00374         // DEBUG: export remapped
00375         std::ostringstream filename;
00376         filename << ioImgInfo._name << "_grey.JPG";
00377         vigra::ImageExportInfo exinfo(filename.str().c_str());
00378         vigra::exportImage(srcImageRange(final_img), exinfo);
00379 #endif
00380 
00381         // Build integral image
00382         TRACE_IMG("Build integral image...");
00383         ioImgInfo._ii.init(final_img);
00384         final_img.resize(0,0);
00385 
00386         // compute distance map
00387         if(final_mask.width()>0 && final_mask.height()>0)
00388         {
00389             TRACE_IMG("Build distance map...");
00390             //apply threshold, in case loaded mask contains other values than 0 and 255
00391             vigra::transformImage(srcImageRange(final_mask), destImage(final_mask),
00392                                   vigra::Threshold<vigra::BImage::PixelType, vigra::BImage::PixelType>(1, 255, 0, 255));
00393             ioImgInfo._distancemap.resize(final_mask.width(),final_mask.height(),0);
00394             vigra::distanceTransform(srcImageRange(final_mask), destImage(ioImgInfo._distancemap), 255, 2);
00395 #ifdef DEBUG_LOADING_REMAPPING
00396             std::ostringstream maskfilename;
00397             maskfilename << ioImgInfo._name << "_mask.JPG";
00398             vigra::ImageExportInfo maskexinfo(maskfilename.str().c_str());
00399             vigra::exportImage(srcImageRange(final_mask), maskexinfo);
00400             std::ostringstream distfilename;
00401             distfilename << ioImgInfo._name << "_distancemap.JPG";
00402             vigra::ImageExportInfo distexinfo(distfilename.str().c_str());
00403             vigra::exportImage(srcImageRange(ioImgInfo._distancemap), distexinfo);
00404 #endif
00405             final_mask.resize(0,0);
00406         };
00407     }
00408     catch (std::exception& e)
00409     {
00410         TRACE_INFO("An error happened while loading image : caught exception: " << e.what() << endl);
00411         ioImgInfo._loadFail=true;
00412         return false;
00413     }
00414 
00415     return true;
00416 }
00417 
00418 
00419 bool PanoDetector::FindKeyPointsInImage(ImgData& ioImgInfo, const PanoDetector& iPanoDetector)
00420 {
00421     TRACE_IMG("Find keypoints...");
00422 
00423     // setup the detector
00424     KeyPointDetector aKP;
00425 
00426     // detect the keypoints
00427     KeyPointVectInsertor aInsertor(ioImgInfo._kp);
00428     aKP.detectKeypoints(ioImgInfo._ii, aInsertor);
00429 
00430     TRACE_IMG("Found "<< ioImgInfo._kp.size() << " interest points.");
00431 
00432     return true;
00433 }
00434 
00435 bool PanoDetector::FilterKeyPointsInImage(ImgData& ioImgInfo, const PanoDetector& iPanoDetector)
00436 {
00437     TRACE_IMG("Filtering keypoints...");
00438 
00439     lfeat::Sieve<lfeat::KeyPointPtr, lfeat::KeyPointPtrSort > aSieve(iPanoDetector.getSieve1Width(),
00440             iPanoDetector.getSieve1Height(),
00441             iPanoDetector.getSieve1Size());
00442 
00443     // insert the points in the Sieve if they are not masked
00444     double aXF = (double)iPanoDetector.getSieve1Width() / (double)ioImgInfo._detectWidth;
00445     double aYF = (double)iPanoDetector.getSieve1Height() / (double)ioImgInfo._detectHeight;
00446 
00447     const bool distmap_valid=(ioImgInfo._distancemap.width()>0 && ioImgInfo._distancemap.height()>0);
00448     for (size_t i = 0; i < ioImgInfo._kp.size(); ++i)
00449     {
00450         KeyPointPtr& aK = ioImgInfo._kp[i];
00451         if(distmap_valid)
00452         {
00453             if(aK->_x > 0 && aK->_x < ioImgInfo._distancemap.width() && aK->_y > 0 && aK->_y < ioImgInfo._distancemap.height()
00454                     && ioImgInfo._distancemap((int)(aK->_x),(int)(aK->_y)) >aK->_scale*8)
00455             {
00456                 //cout << " dist from border:" << ioImgInfo._distancemap((int)(aK->_x),(int)(aK->_y)) << " required dist: " << aK->_scale*12 << std::endl;
00457                 aSieve.insert(aK, (int)(aK->_x * aXF), (int)(aK->_y * aYF));
00458             }
00459         }
00460         else
00461         {
00462             aSieve.insert(aK, (int)(aK->_x * aXF), (int)(aK->_y * aYF));
00463         };
00464     }
00465 
00466     // pull remaining values from the sieve
00467     ioImgInfo._kp.clear();
00468 
00469     // make an extractor and pull the points
00470     SieveExtractorKP aSieveExt(ioImgInfo._kp);
00471     aSieve.extract(aSieveExt);
00472 
00473     TRACE_IMG("Kept " << ioImgInfo._kp.size() << " interest points.");
00474 
00475     return true;
00476 
00477 }
00478 
00479 bool PanoDetector::MakeKeyPointDescriptorsInImage(ImgData& ioImgInfo, const PanoDetector& iPanoDetector)
00480 {
00481     TRACE_IMG("Make keypoint descriptors...");
00482 
00483     // build a keypoint descriptor
00484     CircularKeyPointDescriptor aKPD(ioImgInfo._ii);
00485 
00486     // vector for keypoints with more than one orientation
00487     KeyPointVect_t kp_new_ori;
00488     for (size_t j = 0; j < ioImgInfo._kp.size(); ++j)
00489     {
00490         KeyPointPtr& aK = ioImgInfo._kp[j];
00491         double angles[4];
00492         int nAngles = aKPD.assignOrientation(*aK, angles);
00493         for (int i=0; i < nAngles; i++)
00494         {
00495             // duplicate Keypoint with additional angles
00496             KeyPointPtr aKn = KeyPointPtr ( new lfeat::KeyPoint ( *aK ) );
00497             aKn->_ori = angles[i];
00498             kp_new_ori.push_back(aKn);
00499         }
00500     }
00501     ioImgInfo._kp.insert(ioImgInfo._kp.end(), kp_new_ori.begin(), kp_new_ori.end());
00502 
00503     for (size_t i = 0; i < ioImgInfo._kp.size(); ++i)
00504     {
00505         aKPD.makeDescriptor(*(ioImgInfo._kp[i]));
00506     }
00507     // store the descriptor length
00508     ioImgInfo._descLength = aKPD.getDescriptorLength();
00509     return true;
00510 }
00511 
00512 bool PanoDetector::RemapBackKeypoints(ImgData& ioImgInfo, const PanoDetector& iPanoDetector)
00513 {
00514 
00515     double scale=iPanoDetector._downscale ? 2.0:1.0;
00516 
00517     if (!ioImgInfo._needsremap)
00518     {
00519         if(scale != 1.0)
00520         {
00521             for (size_t i = 0; i < ioImgInfo._kp.size(); ++i)
00522             {
00523                 KeyPointPtr& aK = ioImgInfo._kp[i];
00524                 aK->_x *= scale;
00525                 aK->_y *= scale;
00526                 aK->_scale *= scale;
00527             }
00528         };
00529     }
00530     else
00531     {
00532         TRACE_IMG("Remapping back keypoints...");
00533         HuginBase::PTools::Transform trafo1;
00534         trafo1.createTransform(iPanoDetector._panoramaInfoCopy.getSrcImage(ioImgInfo._number),
00535                                ioImgInfo._projOpts);
00536 
00537         int dx1 = ioImgInfo._projOpts.getROI().left();
00538         int dy1 = ioImgInfo._projOpts.getROI().top();
00539 
00540         for (size_t i = 0; i < ioImgInfo._kp.size(); ++i)
00541         {
00542             KeyPointPtr& aK = ioImgInfo._kp[i];
00543             double xout, yout;
00544             if(trafo1.transformImgCoord(xout, yout, aK->_x + dx1, aK->_y+ dy1))
00545             {
00546                 // downscaling is take care of by the remapping transform
00547                 // no need for multiplying the scale factor...
00548                 aK->_x=xout;
00549                 aK->_y=yout;
00550                 aK->_scale *= scale;
00551             }
00552         }
00553     }
00554     return true;
00555 }
00556 
00557 bool PanoDetector::BuildKDTreesInImage(ImgData& ioImgInfo, const PanoDetector& iPanoDetector)
00558 {
00559     TRACE_IMG("Build KDTree...");
00560 
00561     if(ioImgInfo._kp.size()==0)
00562     {
00563         return false;
00564     };
00565     // build a vector of KDElemKeyPointPtr
00566 
00567     // create feature vector matrix for flann
00568     ioImgInfo._flann_descriptors = flann::Matrix<double>(new double[ioImgInfo._kp.size()*ioImgInfo._descLength],
00569                                    ioImgInfo._kp.size(), ioImgInfo._descLength);
00570     for (size_t i = 0; i < ioImgInfo._kp.size(); ++i)
00571     {
00572         memcpy(ioImgInfo._flann_descriptors[i], ioImgInfo._kp[i]->_vec, sizeof(double)*ioImgInfo._descLength);
00573     }
00574 
00575     // build query structure
00576     ioImgInfo._flann_index = new flann::Index<flann::L2<double> > (ioImgInfo._flann_descriptors, flann::KDTreeIndexParams(4));
00577     ioImgInfo._flann_index->buildIndex();
00578 
00579     return true;
00580 }
00581 
00582 bool PanoDetector::FreeMemoryInImage(ImgData& ioImgInfo, const PanoDetector& iPanoDetector)
00583 {
00584     TRACE_IMG("Freeing memory...");
00585 
00586     ioImgInfo._ii.clean();
00587     ioImgInfo._distancemap.resize(0,0);
00588 
00589     return true;
00590 }
00591 
00592 
00593 bool PanoDetector::FindMatchesInPair(MatchData& ioMatchData, const PanoDetector& iPanoDetector)
00594 {
00595     TRACE_PAIR("Find Matches...");
00596 
00597     // retrieve the KDTree of image 2
00598     flann::Index<flann::L2<double> > * index2 = ioMatchData._i2->_flann_index;
00599 
00600     // retrieve query points from image 1
00601     flann::Matrix<double> & query = ioMatchData._i1->_flann_descriptors;
00602 
00603     // storage for sorted 2 best matches
00604     int nn = 2;
00605     flann::Matrix<int> indices(new int[query.rows*nn], query.rows, nn);
00606     flann::Matrix<double> dists(new double[query.rows*nn], query.rows, nn);
00607 
00608     // perform matching using flann
00609     index2->knnSearch(query, indices, dists, nn, flann::SearchParams(iPanoDetector.getKDTreeSearchSteps()));
00610 
00611     //typedef KDTreeSpace::BestMatch<KDElemKeyPoint>            BM_t;
00612     //std::set<BM_t, std::greater<BM_t> >       aBestMatches;
00613 
00614     // store the matches already found to avoid 2 points in image1
00615     // match the same point in image2
00616     // both matches will be removed.
00617     set<int> aAlreadyMatched;
00618     set<int> aBadMatch;
00619 
00620     // unfiltered vector of matches
00621     typedef std::pair<KeyPointPtr, int> TmpPair_t;
00622     std::vector<TmpPair_t>      aUnfilteredMatches;
00623 
00624     //PointMatchVector_t aMatches;
00625 
00626     // go through all the keypoints of image 1
00627     for (unsigned aKIt = 0; aKIt < query.rows; ++aKIt)
00628     {
00629         // accept the match if the second match is far enough
00630         // put a lower value for stronger matching default 0.15
00631         if (dists[aKIt][0] > iPanoDetector.getKDTreeSecondDistance()  * dists[aKIt][1])
00632         {
00633             continue;
00634         }
00635 
00636         // check if the kdtree match number is already in the already matched set
00637         if (aAlreadyMatched.find(indices[aKIt][0]) != aAlreadyMatched.end())
00638         {
00639             // add to delete list and continue
00640             aBadMatch.insert(indices[aKIt][0]);
00641             continue;
00642         }
00643 
00644         // TODO: add check for duplicate matches (can happen if a keypoint gets multiple orientations)
00645 
00646         // add the match number in already matched set
00647         aAlreadyMatched.insert(indices[aKIt][0]);
00648 
00649         // add the match to the unfiltered list
00650         aUnfilteredMatches.push_back(TmpPair_t(ioMatchData._i1->_kp[aKIt], indices[aKIt][0]));
00651     }
00652 
00653     // now filter and fill the vector of matches
00654     for (size_t i = 0; i < aUnfilteredMatches.size(); ++i)
00655     {
00656         TmpPair_t& aP = aUnfilteredMatches[i];
00657         // if the image2 match number is in the badmatch set, skip it.
00658         if (aBadMatch.find(aP.second) != aBadMatch.end())
00659         {
00660             continue;
00661         }
00662 
00663         // add the match in the output vector
00664         ioMatchData._matches.push_back(lfeat::PointMatchPtr( new lfeat::PointMatch(aP.first, ioMatchData._i2->_kp[aP.second])));
00665     }
00666 
00667     TRACE_PAIR("Found " << ioMatchData._matches.size() << " matches.");
00668     return true;
00669 }
00670 
00671 bool PanoDetector::RansacMatchesInPair(MatchData& ioMatchData, const PanoDetector& iPanoDetector)
00672 {
00673     // Use panotools model for wide angle lenses
00674     RANSACOptimizer::Mode rmode = iPanoDetector._ransacMode;
00675     if (rmode == RANSACOptimizer::HOMOGRAPHY ||
00676             (rmode == RANSACOptimizer::AUTO && iPanoDetector._panoramaInfo->getImage(ioMatchData._i1->_number).getHFOV() < 65 &&
00677              iPanoDetector._panoramaInfo->getImage(ioMatchData._i2->_number).getHFOV() < 65))
00678     {
00679         return RansacMatchesInPairHomography(ioMatchData, iPanoDetector);
00680     }
00681     else
00682     {
00683         return RansacMatchesInPairCam(ioMatchData, iPanoDetector);
00684     }
00685 }
00686 
00687 // new code with fisheye aware ransac
00688 bool PanoDetector::RansacMatchesInPairCam(MatchData& ioMatchData, const PanoDetector& iPanoDetector)
00689 {
00690     TRACE_PAIR("RANSAC Filtering with Panorama model...");
00691 
00692     if (ioMatchData._matches.size() < (unsigned int)iPanoDetector.getMinimumMatches())
00693     {
00694         TRACE_PAIR("Too few matches ... removing all of them.");
00695         ioMatchData._matches.clear();
00696         return true;
00697     }
00698 
00699     if (ioMatchData._matches.size() < 6)
00700     {
00701         TRACE_PAIR("Not enough matches for RANSAC filtering.");
00702         return true;
00703     }
00704 
00705     // setup a panorama project with the two images.
00706     // is this threadsafe (is this read only access?)
00707     UIntSet imgs;
00708     int pano_i1 = ioMatchData._i1->_number;
00709     int pano_i2 = ioMatchData._i2->_number;
00710     imgs.insert(pano_i1);
00711     imgs.insert(pano_i2);
00712     int pano_local_i1 = 0;
00713     int pano_local_i2 = 1;
00714     if (pano_i1 > pano_i2)
00715     {
00716         pano_local_i1 = 1;
00717         pano_local_i2 = 0;
00718     }
00719 
00720     // perform ransac matching.
00721     // ARGH the panotools optimizer uses global variables is not reentrant
00722     std::vector<int> inliers;
00723 #pragma omp critical
00724     {
00725         PanoramaData* panoSubset = iPanoDetector._panoramaInfo->getNewSubset(imgs);
00726 
00727         // create control point vector
00728         CPVector controlPoints(ioMatchData._matches.size());
00729         for (size_t i = 0; i < ioMatchData._matches.size(); ++i)
00730         {
00731             PointMatchPtr& aM=ioMatchData._matches[i];
00732             controlPoints[i] = ControlPoint(pano_local_i1, aM->_img1_x, aM->_img1_y,
00733                                             pano_local_i2, aM->_img2_x, aM->_img2_y);
00734         }
00735         panoSubset->setCtrlPoints(controlPoints);
00736 
00737 
00738         PT_setProgressFcn(ptProgress);
00739         PT_setInfoDlgFcn(ptinfoDlg);
00740 
00741         RANSACOptimizer::Mode rmode = iPanoDetector._ransacMode;
00742         if (rmode == RANSACOptimizer::AUTO)
00743         {
00744             rmode = RANSACOptimizer::RPY;
00745         }
00746         inliers = HuginBase::RANSACOptimizer::findInliers(*panoSubset, pano_local_i1, pano_local_i2,
00747                   iPanoDetector.getRansacDistanceThreshold(), rmode);
00748         PT_setProgressFcn(NULL);
00749         PT_setInfoDlgFcn(NULL);
00750         delete panoSubset;
00751     }
00752 
00753     TRACE_PAIR("Removed " << ioMatchData._matches.size() - inliers.size() << " matches. " << inliers.size() << " remaining.");
00754     if (inliers.size() < 0.5 * ioMatchData._matches.size())
00755     {
00756         // more than 50% of matches were removed, ignore complete pair...
00757         TRACE_PAIR("RANSAC found more than 50% outliers, removing all matches");
00758         ioMatchData._matches.clear();
00759         return true;
00760     }
00761 
00762 
00763     if (inliers.size() < (unsigned int)iPanoDetector.getMinimumMatches())
00764     {
00765         TRACE_PAIR("Too few matches ... removing all of them.");
00766         ioMatchData._matches.clear();
00767         return true;
00768     }
00769 
00770     // keep only inlier matches
00771     PointMatchVector_t aInlierMatches;
00772     aInlierMatches.reserve(inliers.size());
00773 
00774     for (size_t i = 0; i < inliers.size(); ++i)
00775     {
00776         aInlierMatches.push_back(ioMatchData._matches[inliers[i]]);
00777     }
00778     ioMatchData._matches = aInlierMatches;
00779 
00780     /*
00781     if (iPanoDetector.getTest())
00782         TestCode::drawRansacMatches(ioMatchData._i1->_name, ioMatchData._i2->_name, ioMatchData._matches,
00783                                                                 aRemovedMatches, aRansacFilter, iPanoDetector.getDownscale());
00784     */
00785 
00786     return true;
00787 }
00788 
00789 // homography based ransac matching
00790 bool PanoDetector::RansacMatchesInPairHomography(MatchData& ioMatchData, const PanoDetector& iPanoDetector)
00791 {
00792     TRACE_PAIR("RANSAC Filtering...");
00793 
00794     if (ioMatchData._matches.size() < (unsigned int)iPanoDetector.getMinimumMatches())
00795     {
00796         TRACE_PAIR("Too few matches ... removing all of them.");
00797         ioMatchData._matches.clear();
00798         return true;
00799     }
00800 
00801     if (ioMatchData._matches.size() < 6)
00802     {
00803         TRACE_PAIR("Not enough matches for RANSAC filtering.");
00804         return true;
00805     }
00806 
00807     PointMatchVector_t aRemovedMatches;
00808 
00809     Ransac aRansacFilter;
00810     aRansacFilter.setIterations(iPanoDetector.getRansacIterations());
00811     int thresholdDistance=iPanoDetector.getRansacDistanceThreshold();
00812     //increase RANSAC distance if the image were remapped to not exclude
00813     //too much points in this case
00814     if(ioMatchData._i1->_needsremap || ioMatchData._i2->_needsremap)
00815     {
00816         thresholdDistance*=5;
00817     }
00818     aRansacFilter.setDistanceThreshold(thresholdDistance);
00819     aRansacFilter.filter(ioMatchData._matches, aRemovedMatches);
00820 
00821 
00822     TRACE_PAIR("Removed " << aRemovedMatches.size() << " matches. " << ioMatchData._matches.size() << " remaining.");
00823 
00824     if (aRemovedMatches.size() > ioMatchData._matches.size())
00825     {
00826         // more than 50% of matches were removed, ignore complete pair...
00827         TRACE_PAIR("More than 50% outliers, removing all matches");
00828         ioMatchData._matches.clear();
00829         return true;
00830     }
00831 
00832     if (iPanoDetector.getTest())
00833         TestCode::drawRansacMatches(ioMatchData._i1->_name, ioMatchData._i2->_name, ioMatchData._matches,
00834                                     aRemovedMatches, aRansacFilter, iPanoDetector.getDownscale());
00835 
00836     return true;
00837 
00838 }
00839 
00840 
00841 bool PanoDetector::FilterMatchesInPair(MatchData& ioMatchData, const PanoDetector& iPanoDetector)
00842 {
00843     TRACE_PAIR("Clustering matches...");
00844 
00845     if (ioMatchData._matches.size() < 2)
00846     {
00847         return true;
00848     }
00849 
00850     // compute min,max of x,y for image1
00851 
00852     double aMinX = numeric_limits<double>::max();
00853     double aMinY = numeric_limits<double>::max();
00854     double aMaxX = -numeric_limits<double>::max();
00855     double aMaxY = -numeric_limits<double>::max();
00856 
00857     for (size_t i = 0; i < ioMatchData._matches.size(); ++i)
00858     {
00859         PointMatchPtr& aM = ioMatchData._matches[i];
00860         if (aM->_img1_x < aMinX)
00861         {
00862             aMinX = aM->_img1_x;
00863         }
00864         if (aM->_img1_x > aMaxX)
00865         {
00866             aMaxX = aM->_img1_x;
00867         }
00868 
00869         if (aM->_img1_y < aMinY)
00870         {
00871             aMinY = aM->_img1_y;
00872         }
00873         if (aM->_img1_y > aMaxY)
00874         {
00875             aMaxY = aM->_img1_y;
00876         }
00877     }
00878 
00879     double aSizeX = aMaxX - aMinX + 2; // add 2 so max/aSize is strict < 1
00880     double aSizeY = aMaxY - aMinY + 2;
00881 
00882     //
00883 
00884     Sieve<PointMatchPtr, PointMatchPtrSort> aSieve(iPanoDetector.getSieve2Width(),
00885             iPanoDetector.getSieve2Height(),
00886             iPanoDetector.getSieve2Size());
00887 
00888     // insert the points in the Sieve
00889     double aXF = (double)iPanoDetector.getSieve2Width() / aSizeX;
00890     double aYF = (double)iPanoDetector.getSieve2Height() / aSizeY;
00891     for (size_t i = 0; i < ioMatchData._matches.size(); ++i)
00892     {
00893         PointMatchPtr& aM = ioMatchData._matches[i];
00894         aSieve.insert(aM, (int)((aM->_img1_x - aMinX) * aXF), (int)((aM->_img1_y - aMinY) * aYF));
00895     }
00896 
00897     // pull remaining values from the sieve
00898     ioMatchData._matches.clear();
00899 
00900     // make an extractor and pull the points
00901     SieveExtractorMatch aSieveExt(ioMatchData._matches);
00902     aSieve.extract(aSieveExt);
00903 
00904     TRACE_PAIR("Kept " << ioMatchData._matches.size() << " matches.");
00905     return true;
00906 }
00907 
00908 void PanoDetector::writeOutput()
00909 {
00910     // Write output pto file
00911 
00912     ofstream aOut(_outputFile.c_str(), ios_base::trunc);
00913     if( !aOut )
00914     {
00915         cerr << "ERROR : "
00916              << "Couldn't open file '" << _outputFile << "'!" << endl; //STS
00917         return;
00918     }
00919 
00920     aOut << "# pto project file generated by Hugin's cpfind" << endl << endl;
00921 
00922     _panoramaInfo->removeDuplicateCtrlPoints();
00923     AppBase::DocumentData::ReadWriteError err = _panoramaInfo->writeData(aOut);
00924     if (err != AppBase::DocumentData::SUCCESSFUL)
00925     {
00926         cerr << "ERROR couldn't write to output file '" << _outputFile << "'!" << endl;
00927         return;
00928     }
00929 }
00930 
00931 void PanoDetector::writeKeyfile(ImgData& imgInfo)
00932 {
00933     // Write output keyfile
00934 
00935     ofstream aOut(imgInfo._keyfilename.c_str(), ios_base::trunc);
00936 
00937     SIFTFormatWriter writer(aOut);
00938 
00939     int origImgWidth =  _panoramaInfo->getImage(imgInfo._number).getSize().width();
00940     int origImgHeight =  _panoramaInfo->getImage(imgInfo._number).getSize().height();
00941 
00942     ImageInfo img_info(imgInfo._name, origImgWidth, origImgHeight);
00943 
00944     writer.writeHeader ( img_info, imgInfo._kp.size(), imgInfo._descLength );
00945 
00946     for(size_t i=0; i<imgInfo._kp.size(); ++i)
00947     {
00948         KeyPointPtr& aK=imgInfo._kp[i];
00949         writer.writeKeypoint ( aK->_x, aK->_y, aK->_scale, aK->_ori, aK->_score,
00950                                imgInfo._descLength, aK->_vec );
00951     }
00952     writer.writeFooter();
00953 }
00954 

Generated on 3 Aug 2015 for Hugintrunk by  doxygen 1.4.7