1 // -*- c-basic-offset: 4 ; tab-width: 4 -*-
2 /*
3 * Copyright (C) 2007-2008 Anael Orlinski
4 *
5 * This file is part of Panomatic.
6 *
7 * Panomatic is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 2 of the License, or
10 * (at your option) any later version.
11 *
12 * Panomatic is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with Panomatic; if not, write to the Free Software
19 * <http://www.gnu.org/licenses/>.
20 */
21 
22 #include "ImageImport.h"
23 
24 #include "PanoDetector.h"
25 #include <iostream>
26 #include <fstream>
27 #include <vigra/distancetransform.hxx>
28 #include "vigra_ext/impexalpha.hxx"
29 #include "vigra_ext/cms.h"
30 
31 #include <localfeatures/Sieve.h>
32 #include <localfeatures/PointMatch.h>
33 #include <localfeatures/RansacFiltering.h>
34 #include <localfeatures/KeyPointIO.h>
35 #include <localfeatures/CircularKeyPointDescriptor.h>
36 
37 /*
38 #include "KDTree.h"
39 #include "KDTreeImpl.h"
40 */
41 #include "Utils.h"
42 #include "Tracer.h"
43 
44 #include <algorithms/nona/ComputeImageROI.h>
45 #include <algorithms/optimizer/PTOptimizer.h>
46 #include <nona/RemappedPanoImage.h>
47 #include <nona/ImageRemapper.h>
48 
49 #include <time.h>
50 
51 #define TRACE_IMG(X) {if (iPanoDetector.getVerbose() > 1) { TRACE_INFO("i" << ioImgInfo._number << " : " << X << std::endl);} }
52 #define TRACE_PAIR(X) {if (iPanoDetector.getVerbose() > 1){ TRACE_INFO("i" << ioMatchData._i1->_number << " <> " \
53                 "i" << ioMatchData._i2->_number << " : " << X << std::endl)}}
54 
55 // define a Keypoint insertor
56 class KeyPointVectInsertor : public lfeat::KeyPointInsertor
57 {
58 public:
KeyPointVectInsertor(lfeat::KeyPointVect_t & iVect)59     explicit KeyPointVectInsertor(lfeat::KeyPointVect_t& iVect) : _v(iVect) {};
operator ()(const lfeat::KeyPoint & k)60     inline virtual void operator()(const lfeat::KeyPoint& k)
61     {
62         _v.push_back(lfeat::KeyPointPtr(new lfeat::KeyPoint(k)));
63     }
64 
65 private:
66     lfeat::KeyPointVect_t& _v;
67 
68 };
69 
70 
71 // define a sieve extractor
72 class SieveExtractorKP : public lfeat::SieveExtractor<lfeat::KeyPointPtr>
73 {
74 public:
SieveExtractorKP(lfeat::KeyPointVect_t & iV)75     explicit SieveExtractorKP(lfeat::KeyPointVect_t& iV) : _v(iV) {};
operator ()(const lfeat::KeyPointPtr & k)76     inline virtual void operator()(const lfeat::KeyPointPtr& k)
77     {
78         _v.push_back(k);
79     }
80 private:
81     lfeat::KeyPointVect_t& _v;
82 };
83 
84 class SieveExtractorMatch : public lfeat::SieveExtractor<lfeat::PointMatchPtr>
85 {
86 public:
SieveExtractorMatch(lfeat::PointMatchVector_t & iM)87     explicit SieveExtractorMatch(lfeat::PointMatchVector_t& iM) : _m(iM) {};
operator ()(const lfeat::PointMatchPtr & m)88     inline virtual void operator()(const lfeat::PointMatchPtr& m)
89     {
90         _m.push_back(m);
91     }
92 private:
93     lfeat::PointMatchVector_t& _m;
94 };
95 
LoadKeypoints(ImgData & ioImgInfo,const PanoDetector & iPanoDetector)96 bool PanoDetector::LoadKeypoints(ImgData& ioImgInfo, const PanoDetector& iPanoDetector)
97 {
98     TRACE_IMG("Loading keypoints...");
99 
100     lfeat::ImageInfo info = lfeat::loadKeypoints(ioImgInfo._keyfilename, ioImgInfo._kp);
101     ioImgInfo._loadFail = (info.filename.empty());
102 
103     // update ImgData
104     if(ioImgInfo.NeedsRemapping())
105     {
106         ioImgInfo._detectWidth = std::max(info.width,info.height);
107         ioImgInfo._detectHeight = std::max(info.width,info.height);
108         ioImgInfo._projOpts.setWidth(ioImgInfo._detectWidth);
109         ioImgInfo._projOpts.setHeight(ioImgInfo._detectHeight);
110     }
111     else
112     {
113         ioImgInfo._detectWidth = info.width;
114         ioImgInfo._detectHeight = info.height;
115     };
116     ioImgInfo._descLength = info.dimensions;
117 
118     return true;
119 }
120 
121 /** apply the mask and the crop of the given SrcImg to given mask image */
122 template <class SrcImageIterator, class SrcAccessor>
applyMaskAndCrop(vigra::triple<SrcImageIterator,SrcImageIterator,SrcAccessor> img,const HuginBase::SrcPanoImage & SrcImg)123 void applyMaskAndCrop(vigra::triple<SrcImageIterator, SrcImageIterator, SrcAccessor> img, const HuginBase::SrcPanoImage& SrcImg)
124 {
125     vigra::Diff2D imgSize = img.second - img.first;
126 
127     // create dest y iterator
128     SrcImageIterator yd(img.first);
129     // loop over the image and transform
130     for(int y=0; y < imgSize.y; ++y, ++yd.y)
131     {
132         // create x iterators
133         SrcImageIterator xd(yd);
134         for(int x=0; x < imgSize.x; ++x, ++xd.x)
135         {
136             if(!SrcImg.isInside(vigra::Point2D(x,y)))
137             {
138                 *xd=0;
139             };
140         }
141     }
142 }
143 
144 /** functor to scale image on the fly during other operations */
145 template <class T>
146 struct ScaleFunctor
147 {
148     typedef T result_type;
ScaleFunctorScaleFunctor149     explicit ScaleFunctor(double scale) { m_scale = scale; };
150 
operator ()ScaleFunctor151     T operator()(const T & a) const
152     {
153         return m_scale*a;
154     }
155 
156     template <class T2>
operator ()ScaleFunctor157     T2 operator()(const T2 & a, const hugin_utils::FDiff2D & p) const
158     {
159         return m_scale*a;
160     }
161 
162     template <class T2, class A>
hdrWeightScaleFunctor163     A hdrWeight(T2 v, A a) const
164     {
165         return a;
166     }
167 
168 private:
169     double m_scale;
170 };
171 
172 /** helper function to remap image to given projection, you can supply a pixelTransform,
173  *  which will be applied during remapping, this is intended for scaling a image during
174  *  remapping, but this means also, that no photometric corrections are applied,
175  *  if this is wanted you need to supply a suitable pixelTransform */
176 template <class ImageType, class PixelTransform>
RemapImage(const HuginBase::SrcPanoImage & srcImage,const HuginBase::PanoramaOptions & options,size_t detectWidth,size_t detectHeight,ImageType * & image,vigra::BImage * & mask,const PixelTransform & pixelTransform,ImageType * & finalImage,vigra::BImage * & finalMask)177 void RemapImage(const HuginBase::SrcPanoImage& srcImage, const HuginBase::PanoramaOptions& options,
178     size_t detectWidth, size_t detectHeight,
179     ImageType*& image, vigra::BImage*& mask,
180     const PixelTransform& pixelTransform,
181     ImageType*& finalImage, vigra::BImage*& finalMask)
182 {
183     AppBase::DummyProgressDisplay dummy;
184     HuginBase::PTools::Transform transform;
185     transform.createTransform(srcImage, options);
186     finalImage = new ImageType(detectWidth, detectHeight);
187     finalMask = new vigra::BImage(detectWidth, detectHeight, vigra::UInt8(0));
188     if (srcImage.hasActiveMasks() || (srcImage.getCropMode() != HuginBase::SrcPanoImage::NO_CROP && !srcImage.getCropRect().isEmpty()))
189     {
190         if (!mask)
191         {
192             // image has no mask, create full mask
193             mask = new vigra::BImage(image->size(), vigra::UInt8(255));
194         };
195         applyMaskAndCrop(vigra::destImageRange(*mask), srcImage);
196     };
197     if (mask)
198     {
199         vigra_ext::transformImageAlpha(vigra::srcImageRange(*image), vigra::srcImage(*mask), vigra::destImageRange(*finalImage), vigra::destImage(*finalMask),
200             options.getROI().upperLeft(), transform, pixelTransform, false, vigra_ext::INTERP_CUBIC, &dummy);
201         delete mask;
202         mask = NULL;
203     }
204     else
205     {
206         vigra_ext::transformImage(vigra::srcImageRange(*image), vigra::destImageRange(*finalImage), vigra::destImage(*finalMask),
207             options.getROI().upperLeft(), transform, pixelTransform, false, vigra_ext::INTERP_CUBIC, &dummy);
208     };
209     delete image;
210     image = NULL;
211 }
212 
213 /** downscale image if requested, optimized code for non-downscale version to prevent unnecessary copying
214  *  the image data */
215 template <class ImageType>
HandleDownscaleImage(const HuginBase::SrcPanoImage & srcImage,ImageType * & image,vigra::BImage * & mask,size_t detectWidth,size_t detectHeight,bool downscale,ImageType * & finalImage,vigra::BImage * & finalMask)216 void HandleDownscaleImage(const HuginBase::SrcPanoImage& srcImage, ImageType*& image, vigra::BImage*& mask,
217     size_t detectWidth, size_t detectHeight, bool downscale,
218     ImageType*& finalImage, vigra::BImage*& finalMask)
219 {
220     if (srcImage.hasActiveMasks() || (srcImage.getCropMode() != HuginBase::SrcPanoImage::NO_CROP && !srcImage.getCropRect().isEmpty()))
221     {
222         if (!mask)
223         {
224             // image has no mask, create full mask
225             mask = new vigra::BImage(image->size(), vigra::UInt8(255));
226         };
227         //copy mask and crop from pto file into alpha layer
228         applyMaskAndCrop(vigra::destImageRange(*mask), srcImage);
229     };
230     if (downscale)
231     {
232         // Downscale image
233         finalImage = new ImageType(detectWidth, detectHeight);
234         vigra::resizeImageNoInterpolation(vigra::srcImageRange(*image), vigra::destImageRange(*finalImage));
235         delete image;
236         image = NULL;
237         //downscale mask
238         if (mask)
239         {
240             finalMask = new vigra::BImage(detectWidth, detectHeight);
241             vigra::resizeImageNoInterpolation(vigra::srcImageRange(*mask), vigra::destImageRange(*finalMask));
242             delete mask;
243             mask = NULL;
244         };
245     }
246     else
247     {
248         // simply copy pointer instead of copying the whole image data
249         finalImage = image;
250         if (mask)
251         {
252             finalMask = mask;
253         };
254     };
255 };
256 
257 // save some intermediate images to disc if defined
258 // #define DEBUG_LOADING_REMAPPING
AnalyzeImage(ImgData & ioImgInfo,const PanoDetector & iPanoDetector)259 bool PanoDetector::AnalyzeImage(ImgData& ioImgInfo, const PanoDetector& iPanoDetector)
260 {
261     vigra::DImage* final_img = NULL;
262     vigra::BImage* final_mask = NULL;
263 
264     try
265     {
266         ioImgInfo._loadFail=false;
267 
268         TRACE_IMG("Load image...");
269         vigra::ImageImportInfo aImageInfo(ioImgInfo._name.c_str());
270         if (aImageInfo.numExtraBands() > 1)
271         {
272             TRACE_INFO("Image with multiple alpha channels are not supported");
273             ioImgInfo._loadFail = true;
274             return false;
275         };
276         // remark: it would be possible to handle all cases with the same code
277         // but this would mean that in some cases there are unnecessary
278         // range conversions and image data copying actions needed
279         // so we use specialed code for several cases to reduce memory usage
280         // and prevent unnecessary range adaptions
281         if (aImageInfo.isGrayscale())
282         {
283             // gray scale image
284             vigra::DImage* image = new vigra::DImage(aImageInfo.size());
285             vigra::BImage* mask = NULL;
286             // load gray scale image
287             if (aImageInfo.numExtraBands() == 1)
288             {
289                 mask=new vigra::BImage(aImageInfo.size());
290                 vigra::importImageAlpha(aImageInfo, vigra::destImage(*image), vigra::destImage(*mask));
291             }
292             else
293             {
294                 vigra::importImage(aImageInfo, vigra::destImage(*image));
295             };
296             // adopt range
297             double minVal = 0;
298             double maxVal;
299             if (aImageInfo.getPixelType() == std::string("FLOAT") || aImageInfo.getPixelType() == std::string("DOUBLE"))
300             {
301                 vigra::FindMinMax<float> minmax;   // init functor
302                 vigra::inspectImage(vigra::srcImageRange(*image), minmax);
303                 minVal = minmax.min;
304                 maxVal = minmax.max;
305             }
306             else
307             {
308                 maxVal = vigra_ext::getMaxValForPixelType(aImageInfo.getPixelType());
309             };
310             bool range255 = (fabs(maxVal - 255) < 0.01 && fabs(minVal) < 0.01);
311             if (aImageInfo.getICCProfile().empty())
312             {
313                 // no icc profile, cpfind expects images in 0 ..255 range
314                 TRACE_IMG("Rescale range...");
315                 if (!range255)
316                 {
317                     vigra::transformImage(vigra::srcImageRange(*image), vigra::destImage(*image),
318                         vigra::linearRangeMapping(minVal, maxVal, 0.0, 255.0));
319                 };
320                 range255 = true;
321             }
322             else
323             {
324                 // apply ICC profile
325                 TRACE_IMG("Applying icc profile...");
326                 // lcms expects for double datatype all values between 0 and 1
327                 vigra::transformImage(vigra::srcImageRange(*image), vigra::destImage(*image),
328                     vigra::linearRangeMapping(minVal, maxVal, 0.0, 1.0));
329                 range255 = false;
330                 HuginBase::Color::ApplyICCProfile(*image, aImageInfo.getICCProfile(), TYPE_GRAY_DBL);
331             };
332             if (ioImgInfo.NeedsRemapping())
333             {
334                 // remap image
335                 TRACE_IMG("Remapping image...");
336                 if (range255)
337                 {
338                     RemapImage(iPanoDetector._panoramaInfoCopy.getImage(ioImgInfo._number), ioImgInfo._projOpts,
339                         ioImgInfo._detectWidth, ioImgInfo._detectHeight, image, mask, vigra_ext::PassThroughFunctor<double>(),
340                         final_img, final_mask);
341                 }
342                 else
343                 {
344                     // images has been scaled to 0..1 range before, scale back to 0..255 range
345                     RemapImage(iPanoDetector._panoramaInfoCopy.getImage(ioImgInfo._number), ioImgInfo._projOpts,
346                         ioImgInfo._detectWidth, ioImgInfo._detectHeight, image, mask, ScaleFunctor<double>(255.0),
347                         final_img, final_mask);
348                 };
349             }
350             else
351             {
352                 if (range255)
353                 {
354                     TRACE_IMG("Downscale and transform to suitable grayscale...");
355                     HandleDownscaleImage(iPanoDetector._panoramaInfoCopy.getImage(ioImgInfo._number), image, mask,
356                         ioImgInfo._detectWidth, ioImgInfo._detectHeight, ioImgInfo.IsDownscale(),
357                         final_img, final_mask);
358                 }
359                 else
360                 {
361                     TRACE_IMG("Transform to suitable grayscale...");
362                     HandleDownscaleImage(iPanoDetector._panoramaInfoCopy.getImage(ioImgInfo._number), image, mask,
363                         ioImgInfo._detectWidth, ioImgInfo._detectHeight, ioImgInfo.IsDownscale(),
364                         final_img, final_mask);
365                     vigra::transformImage(vigra::srcImageRange(*final_img), vigra::destImage(*final_img), vigra::linearRangeMapping(0, 1, 0, 255));
366                 };
367             };
368             if (iPanoDetector.getCeleste())
369             {
370                 TRACE_IMG("Celeste does not work with grayscale images. Skipping...");
371             };
372         }
373         else
374         {
375             if (aImageInfo.isColor())
376             {
377                 // rgb images
378                 // prepare radius parameter for celeste
379                 int radius = 1;
380                 if (iPanoDetector.getCeleste())
381                 {
382                     radius = iPanoDetector.getCelesteRadius();
383                     if (iPanoDetector._downscale)
384                     {
385                         radius >>= 1;
386                     };
387                     if (radius < 2)
388                     {
389                         radius = 2;
390                     };
391                 };
392                 switch (aImageInfo.pixelType())
393                 {
394                     case vigra::ImageImportInfo::UINT8:
395                         // special variant for unsigned 8 bit images
396                         {
397                             vigra::BRGBImage* rgbImage=new vigra::BRGBImage(aImageInfo.size());
398                             vigra::BImage* mask = NULL;
399                             // load image
400                             if (aImageInfo.numExtraBands() == 1)
401                             {
402                                 mask=new vigra::BImage(aImageInfo.size());
403                                 vigra::importImageAlpha(aImageInfo, vigra::destImage(*rgbImage), vigra::destImage(*mask));
404                             }
405                             else
406                             {
407                                 vigra::importImage(aImageInfo, vigra::destImage(*rgbImage));
408                             };
409                             // apply icc profile
410                             if (!aImageInfo.getICCProfile().empty())
411                             {
412                                 TRACE_IMG("Applying icc profile...");
413                                 HuginBase::Color::ApplyICCProfile(*rgbImage, aImageInfo.getICCProfile(), TYPE_RGB_8);
414                             };
415                             vigra::BRGBImage* scaled = NULL;
416                             if (ioImgInfo.NeedsRemapping())
417                             {
418                                 // remap image
419                                 TRACE_IMG("Remapping image...");
420                                 RemapImage(iPanoDetector._panoramaInfoCopy.getImage(ioImgInfo._number), ioImgInfo._projOpts,
421                                     ioImgInfo._detectWidth, ioImgInfo._detectHeight, rgbImage, mask,
422                                     vigra_ext::PassThroughFunctor<vigra::RGBValue<vigra::UInt8> >(),
423                                     scaled, final_mask);
424                             }
425                             else
426                             {
427                                 if (ioImgInfo.IsDownscale())
428                                 {
429                                     TRACE_IMG("Downscale image...");
430                                 };
431                                 HandleDownscaleImage(iPanoDetector._panoramaInfoCopy.getImage(ioImgInfo._number), rgbImage, mask,
432                                     ioImgInfo._detectWidth, ioImgInfo._detectHeight, ioImgInfo.IsDownscale(),
433                                     scaled, final_mask);
434                             };
435                             if (iPanoDetector.getCeleste())
436                             {
437                                 TRACE_IMG("Mask areas with clouds...");
438                                 vigra::UInt16RGBImage* image16=new vigra::UInt16RGBImage(scaled->size());
439                                 vigra::transformImage(vigra::srcImageRange(*scaled), vigra::destImage(*image16),
440                                     vigra::linearIntensityTransform<vigra::RGBValue<vigra::UInt16> >(255));
441                                 vigra::BImage* celeste_mask = celeste::getCelesteMask(iPanoDetector.svmModel, *image16, radius, iPanoDetector.getCelesteThreshold(), 800, true, false);
442 #ifdef DEBUG_LOADING_REMAPPING
443                                 // DEBUG: export celeste mask
444                                 std::ostringstream maskfilename;
445                                 maskfilename << ioImgInfo._name << "_celeste_mask.JPG";
446                                 vigra::ImageExportInfo maskexinfo(maskfilename.str().c_str());
447                                 vigra::exportImage(vigra::srcImageRange(*celeste_mask), maskexinfo);
448 #endif
449                                 delete image16;
450                                 if (final_mask)
451                                 {
452                                     vigra::copyImageIf(vigra::srcImageRange(*celeste_mask), vigra::srcImage(*final_mask), vigra::destImage(*final_mask));
453                                 }
454                                 else
455                                 {
456                                     final_mask = celeste_mask;
457                                 };
458                             };
459                             // scale to greyscale
460                             TRACE_IMG("Convert to greyscale double...");
461                             final_img = new vigra::DImage(scaled->size());
462                             vigra::copyImage(vigra::srcImageRange(*scaled, vigra::RGBToGrayAccessor<vigra::RGBValue<vigra::UInt8> >()),
463                                 vigra::destImage(*final_img));
464                             delete scaled;
465                         };
466                         break;
467                     case vigra::ImageImportInfo::UINT16:
468                         // special variant for unsigned 16 bit images
469                         {
470                             vigra::UInt16RGBImage* rgbImage = new vigra::UInt16RGBImage(aImageInfo.size());
471                             vigra::BImage* mask = NULL;
472                             // load image
473                             if (aImageInfo.numExtraBands() == 1)
474                             {
475                                 mask = new vigra::BImage(aImageInfo.size());
476                                 vigra::importImageAlpha(aImageInfo, vigra::destImage(*rgbImage), vigra::destImage(*mask));
477                             }
478                             else
479                             {
480                                 vigra::importImage(aImageInfo, vigra::destImage(*rgbImage));
481                             };
482                             // apply icc profile
483                             if (!aImageInfo.getICCProfile().empty())
484                             {
485                                 TRACE_IMG("Applying icc profile...");
486                                 HuginBase::Color::ApplyICCProfile(*rgbImage, aImageInfo.getICCProfile(), TYPE_RGB_16);
487                             };
488                             vigra::UInt16RGBImage* scaled = NULL;
489                             if (ioImgInfo.NeedsRemapping())
490                             {
491                                 // remap image
492                                 TRACE_IMG("Remapping image...");
493                                 RemapImage(iPanoDetector._panoramaInfoCopy.getImage(ioImgInfo._number), ioImgInfo._projOpts,
494                                     ioImgInfo._detectWidth, ioImgInfo._detectHeight, rgbImage, mask,
495                                     vigra_ext::PassThroughFunctor<vigra::RGBValue<vigra::UInt16> >(),
496                                     scaled, final_mask);
497                             }
498                             else
499                             {
500                                 if (ioImgInfo.IsDownscale())
501                                 {
502                                     TRACE_IMG("Downscale image...");
503                                 };
504                                 HandleDownscaleImage(iPanoDetector._panoramaInfoCopy.getImage(ioImgInfo._number), rgbImage, mask,
505                                     ioImgInfo._detectWidth, ioImgInfo._detectHeight, ioImgInfo.IsDownscale(),
506                                     scaled, final_mask);
507                             };
508                             if (iPanoDetector.getCeleste())
509                             {
510                                 TRACE_IMG("Mask areas with clouds...");
511                                 vigra::BImage* celeste_mask = celeste::getCelesteMask(iPanoDetector.svmModel, *scaled, radius, iPanoDetector.getCelesteThreshold(), 800, true, false);
512 #ifdef DEBUG_LOADING_REMAPPING
513                                 // DEBUG: export celeste mask
514                                 std::ostringstream maskfilename;
515                                 maskfilename << ioImgInfo._name << "_celeste_mask.JPG";
516                                 vigra::ImageExportInfo maskexinfo(maskfilename.str().c_str());
517                                 vigra::exportImage(vigra::srcImageRange(*celeste_mask), maskexinfo);
518 #endif
519                                 if (final_mask)
520                                 {
521                                     vigra::copyImageIf(vigra::srcImageRange(*celeste_mask), vigra::srcImage(*final_mask), vigra::destImage(*final_mask));
522                                 }
523                                 else
524                                 {
525                                     final_mask = celeste_mask;
526                                 };
527                             };
528                             // scale to greyscale
529                             TRACE_IMG("Convert to greyscale double...");
530                             final_img = new vigra::DImage(scaled->size());
531                             // keypoint finder expext 0..255 range
532                             vigra::transformImage(vigra::srcImageRange(*scaled, vigra::RGBToGrayAccessor<vigra::RGBValue<vigra::UInt16> >()),
533                                 vigra::destImage(*final_img), vigra::functor::Arg1() / vigra::functor::Param(255.0));
534                             delete scaled;
535                         };
536                         break;
537                     default:
538                         // double variant for all other cases
539                         {
540                             vigra::DRGBImage* rgbImage = new vigra::DRGBImage(aImageInfo.size());
541                             vigra::BImage* mask = NULL;
542                             // load image
543                             if (aImageInfo.numExtraBands() == 1)
544                             {
545                                 mask = new vigra::BImage(aImageInfo.size());
546                                 vigra::importImageAlpha(aImageInfo, vigra::destImage(*rgbImage), vigra::destImage(*mask));
547                             }
548                             else
549                             {
550                                 vigra::importImage(aImageInfo, vigra::destImage(*rgbImage));
551                             };
552                             // range adaption
553                             double minVal = 0;
554                             double maxVal;
555                             const bool isDouble = aImageInfo.getPixelType() == std::string("FLOAT") || aImageInfo.getPixelType() == std::string("DOUBLE");
556                             if (isDouble)
557                             {
558                                 vigra::FindMinMax<float> minmax;   // init functor
559                                 vigra::inspectImage(vigra::srcImageRange(*rgbImage, vigra::RGBToGrayAccessor<vigra::RGBValue<double> >()), minmax);
560                                 minVal = minmax.min;
561                                 maxVal = minmax.max;
562                             }
563                             else
564                             {
565                                 maxVal = vigra_ext::getMaxValForPixelType(aImageInfo.getPixelType());
566                             };
567                             bool range255 = (fabs(maxVal - 255) < 0.01 && fabs(minVal) < 0.01);
568                             if (aImageInfo.getICCProfile().empty())
569                             {
570                                 // no icc profile, cpfind expects images in 0 ..255 range
571                                 TRACE_IMG("Rescale range...");
572                                 if (!range255)
573                                 {
574                                     int mapping = 0;
575                                     if (isDouble && iPanoDetector._panoramaInfoCopy.getImage(ioImgInfo._number).getResponseType() == HuginBase::BaseSrcPanoImage::RESPONSE_LINEAR)
576                                     {
577                                         // switch to log mapping for double/float images with linear response type
578                                         mapping = 1;
579                                     };
580                                     vigra_ext::applyMapping(vigra::srcImageRange(*rgbImage), vigra::destImage(*rgbImage), minVal, maxVal, mapping);
581                                 };
582                                 range255 = true;
583                             }
584                             else
585                             {
586                                 // apply ICC profile
587                                 TRACE_IMG("Applying icc profile...");
588                                 // lcms expects for double datatype all values between 0 and 1
589                                 vigra::transformImage(vigra::srcImageRange(*rgbImage), vigra::destImage(*rgbImage),
590                                     vigra_ext::LinearTransform<vigra::RGBValue<double> >(1.0 / maxVal - minVal, -minVal));
591                                 range255 = false;
592                                 HuginBase::Color::ApplyICCProfile(*rgbImage, aImageInfo.getICCProfile(), TYPE_RGB_DBL);
593                             };
594                             vigra::DRGBImage* scaled;
595                             if (ioImgInfo.NeedsRemapping())
596                             {
597                                 // remap image
598                                 TRACE_IMG("Remapping image...");
599                                 RemapImage(iPanoDetector._panoramaInfoCopy.getImage(ioImgInfo._number), ioImgInfo._projOpts,
600                                     ioImgInfo._detectWidth, ioImgInfo._detectHeight, rgbImage, mask, vigra_ext::PassThroughFunctor<double>(),
601                                     scaled, final_mask);
602                             }
603                             else
604                             {
605                                 TRACE_IMG("Transform to suitable grayscale...");
606                                 HandleDownscaleImage(iPanoDetector._panoramaInfoCopy.getImage(ioImgInfo._number), rgbImage, mask,
607                                     ioImgInfo._detectWidth, ioImgInfo._detectHeight, ioImgInfo.IsDownscale(),
608                                     scaled, final_mask);
609                             };
610                             if (iPanoDetector.getCeleste())
611                             {
612                                 TRACE_IMG("Mask areas with clouds...");
613                                 vigra::UInt16RGBImage* image16 = new vigra::UInt16RGBImage(scaled->size());
614                                 if (range255)
615                                 {
616                                     vigra::transformImage(vigra::srcImageRange(*scaled), vigra::destImage(*image16),
617                                         vigra::linearIntensityTransform<vigra::RGBValue<vigra::UInt16> >(255));
618                                 }
619                                 else
620                                 {
621                                     vigra::transformImage(vigra::srcImageRange(*scaled), vigra::destImage(*image16),
622                                         vigra::linearIntensityTransform<vigra::RGBValue<vigra::UInt16> >(65535));
623                                 };
624                                 vigra::BImage* celeste_mask = celeste::getCelesteMask(iPanoDetector.svmModel, *image16, radius, iPanoDetector.getCelesteThreshold(), 800, true, false);
625 #ifdef DEBUG_LOADING_REMAPPING
626                                 // DEBUG: export celeste mask
627                                 std::ostringstream maskfilename;
628                                 maskfilename << ioImgInfo._name << "_celeste_mask.JPG";
629                                 vigra::ImageExportInfo maskexinfo(maskfilename.str().c_str());
630                                 vigra::exportImage(vigra::srcImageRange(*celeste_mask), maskexinfo);
631 #endif
632                                 delete image16;
633                                 if (final_mask)
634                                 {
635                                     vigra::copyImageIf(vigra::srcImageRange(*celeste_mask), vigra::srcImage(*final_mask), vigra::destImage(*final_mask));
636                                 }
637                                 else
638                                 {
639                                     final_mask = celeste_mask;
640                                 };
641                             };
642                             // scale to greyscale
643                             TRACE_IMG("Convert to greyscale double...");
644                             final_img = new vigra::DImage(scaled->size());
645                             // keypoint finder expext 0..255 range
646                             if (range255)
647                             {
648                                 vigra::copyImage(vigra::srcImageRange(*scaled, vigra::RGBToGrayAccessor<vigra::RGBValue<double> >()), vigra::destImage(*final_img));
649                             }
650                             else
651                             {
652                                 vigra::transformImage(vigra::srcImageRange(*scaled, vigra::RGBToGrayAccessor<vigra::RGBValue<double> >()),
653                                     vigra::destImage(*final_img), vigra::functor::Arg1() * vigra::functor::Param(255.0));
654                             };
655                             delete scaled;
656                         };
657                         break;
658                 };
659             }
660             else
661             {
662                 TRACE_INFO("Cpfind works only with grayscale or RGB images");
663                 ioImgInfo._loadFail = true;
664                 return false;
665             };
666         };
667 
668 #ifdef DEBUG_LOADING_REMAPPING
669         // DEBUG: export remapped
670         std::ostringstream filename;
671         filename << ioImgInfo._name << "_grey.JPG";
672         vigra::ImageExportInfo exinfo(filename.str().c_str());
673         vigra::exportImage(vigra::srcImageRange(*final_img), exinfo);
674 #endif
675 
676         // Build integral image
677         TRACE_IMG("Build integral image...");
678         ioImgInfo._ii.init(*final_img);
679         delete final_img;
680 
681         // compute distance map
682         if(final_mask)
683         {
684             TRACE_IMG("Build distance map...");
685             //apply threshold, in case loaded mask contains other values than 0 and 255
686             vigra::transformImage(vigra::srcImageRange(*final_mask), vigra::destImage(*final_mask),
687                                   vigra::Threshold<vigra::BImage::PixelType, vigra::BImage::PixelType>(1, 255, 0, 255));
688             ioImgInfo._distancemap.resize(final_mask->width(), final_mask->height(), 0);
689             vigra::distanceTransform(vigra::srcImageRange(*final_mask), vigra::destImage(ioImgInfo._distancemap), 255, 2);
690 #ifdef DEBUG_LOADING_REMAPPING
691             std::ostringstream maskfilename;
692             maskfilename << ioImgInfo._name << "_mask.JPG";
693             vigra::ImageExportInfo maskexinfo(maskfilename.str().c_str());
694             vigra::exportImage(vigra::srcImageRange(*final_mask), maskexinfo);
695             std::ostringstream distfilename;
696             distfilename << ioImgInfo._name << "_distancemap.JPG";
697             vigra::ImageExportInfo distexinfo(distfilename.str().c_str());
698             vigra::exportImage(vigra::srcImageRange(ioImgInfo._distancemap), distexinfo);
699 #endif
700             delete final_mask;
701         };
702     }
703     catch (std::exception& e)
704     {
705         TRACE_INFO("An error happened while loading image : caught exception: " << e.what() << std::endl);
706         ioImgInfo._loadFail=true;
707         return false;
708     }
709 
710     return true;
711 }
712 
713 
FindKeyPointsInImage(ImgData & ioImgInfo,const PanoDetector & iPanoDetector)714 bool PanoDetector::FindKeyPointsInImage(ImgData& ioImgInfo, const PanoDetector& iPanoDetector)
715 {
716     TRACE_IMG("Find keypoints...");
717 
718     // setup the detector
719     KeyPointDetector aKP;
720 
721     // detect the keypoints
722     KeyPointVectInsertor aInsertor(ioImgInfo._kp);
723     aKP.detectKeypoints(ioImgInfo._ii, aInsertor);
724 
725     TRACE_IMG("Found "<< ioImgInfo._kp.size() << " interest points.");
726 
727     return true;
728 }
729 
FilterKeyPointsInImage(ImgData & ioImgInfo,const PanoDetector & iPanoDetector)730 bool PanoDetector::FilterKeyPointsInImage(ImgData& ioImgInfo, const PanoDetector& iPanoDetector)
731 {
732     TRACE_IMG("Filtering keypoints...");
733 
734     lfeat::Sieve<lfeat::KeyPointPtr, lfeat::KeyPointPtrSort > aSieve(iPanoDetector.getSieve1Width(),
735             iPanoDetector.getSieve1Height(),
736             iPanoDetector.getSieve1Size());
737 
738     // insert the points in the Sieve if they are not masked
739     double aXF = (double)iPanoDetector.getSieve1Width() / (double)ioImgInfo._detectWidth;
740     double aYF = (double)iPanoDetector.getSieve1Height() / (double)ioImgInfo._detectHeight;
741 
742     const bool distmap_valid=(ioImgInfo._distancemap.width()>0 && ioImgInfo._distancemap.height()>0);
743     for (size_t i = 0; i < ioImgInfo._kp.size(); ++i)
744     {
745         lfeat::KeyPointPtr& aK = ioImgInfo._kp[i];
746         if(distmap_valid)
747         {
748             if(aK->_x > 0 && aK->_x < ioImgInfo._distancemap.width() && aK->_y > 0 && aK->_y < ioImgInfo._distancemap.height()
749                     && ioImgInfo._distancemap((int)(aK->_x),(int)(aK->_y)) >aK->_scale*8)
750             {
751                 //cout << " dist from border:" << ioImgInfo._distancemap((int)(aK->_x),(int)(aK->_y)) << " required dist: " << aK->_scale*12 << std::endl;
752                 aSieve.insert(aK, (int)(aK->_x * aXF), (int)(aK->_y * aYF));
753             }
754         }
755         else
756         {
757             aSieve.insert(aK, (int)(aK->_x * aXF), (int)(aK->_y * aYF));
758         };
759     }
760 
761     // pull remaining values from the sieve
762     ioImgInfo._kp.clear();
763 
764     // make an extractor and pull the points
765     SieveExtractorKP aSieveExt(ioImgInfo._kp);
766     aSieve.extract(aSieveExt);
767 
768     TRACE_IMG("Kept " << ioImgInfo._kp.size() << " interest points.");
769 
770     return true;
771 
772 }
773 
MakeKeyPointDescriptorsInImage(ImgData & ioImgInfo,const PanoDetector & iPanoDetector)774 bool PanoDetector::MakeKeyPointDescriptorsInImage(ImgData& ioImgInfo, const PanoDetector& iPanoDetector)
775 {
776     TRACE_IMG("Make keypoint descriptors...");
777 
778     // build a keypoint descriptor
779     lfeat::CircularKeyPointDescriptor aKPD(ioImgInfo._ii);
780 
781     // vector for keypoints with more than one orientation
782     lfeat::KeyPointVect_t kp_new_ori;
783     for (size_t j = 0; j < ioImgInfo._kp.size(); ++j)
784     {
785         lfeat::KeyPointPtr& aK = ioImgInfo._kp[j];
786         double angles[4];
787         int nAngles = aKPD.assignOrientation(*aK, angles);
788         for (int i=0; i < nAngles; i++)
789         {
790             // duplicate Keypoint with additional angles
791             lfeat::KeyPointPtr aKn = lfeat::KeyPointPtr ( new lfeat::KeyPoint ( *aK ) );
792             aKn->_ori = angles[i];
793             kp_new_ori.push_back(aKn);
794         }
795     }
796     ioImgInfo._kp.insert(ioImgInfo._kp.end(), kp_new_ori.begin(), kp_new_ori.end());
797 
798     for (size_t i = 0; i < ioImgInfo._kp.size(); ++i)
799     {
800         aKPD.makeDescriptor(*(ioImgInfo._kp[i]));
801     }
802     // store the descriptor length
803     ioImgInfo._descLength = aKPD.getDescriptorLength();
804     return true;
805 }
806 
RemapBackKeypoints(ImgData & ioImgInfo,const PanoDetector & iPanoDetector)807 bool PanoDetector::RemapBackKeypoints(ImgData& ioImgInfo, const PanoDetector& iPanoDetector)
808 {
809     if (ioImgInfo.IsDownscale())
810     {
811         for (size_t i = 0; i < ioImgInfo._kp.size(); ++i)
812         {
813             lfeat::KeyPointPtr& aK = ioImgInfo._kp[i];
814             aK->_x *= 2.0;
815             aK->_y *= 2.0;
816             aK->_scale *= 2.0;
817         };
818     }
819     else
820     {
821         if (ioImgInfo.NeedsRemapping())
822         {
823             TRACE_IMG("Remapping back keypoints...");
824             HuginBase::PTools::Transform trafo1;
825             trafo1.createTransform(iPanoDetector._panoramaInfoCopy.getSrcImage(ioImgInfo._number),
826                 ioImgInfo._projOpts);
827 
828             int dx1 = ioImgInfo._projOpts.getROI().left();
829             int dy1 = ioImgInfo._projOpts.getROI().top();
830 
831             for (size_t i = 0; i < ioImgInfo._kp.size(); ++i)
832             {
833                 lfeat::KeyPointPtr& aK = ioImgInfo._kp[i];
834                 double xout, yout;
835                 if (trafo1.transformImgCoord(xout, yout, aK->_x + dx1, aK->_y + dy1))
836                 {
837                     // downscaling is take care of by the remapping transform
838                     // no need for multiplying the scale factor...
839                     aK->_x = xout;
840                     aK->_y = yout;
841                 };
842             };
843         };
844     };
845     return true;
846 }
847 
BuildKDTreesInImage(ImgData & ioImgInfo,const PanoDetector & iPanoDetector)848 bool PanoDetector::BuildKDTreesInImage(ImgData& ioImgInfo, const PanoDetector& iPanoDetector)
849 {
850     TRACE_IMG("Build KDTree...");
851 
852     if(ioImgInfo._kp.empty())
853     {
854         return false;
855     };
856     // build a vector of KDElemKeyPointPtr
857 
858     // create feature vector matrix for flann
859     ioImgInfo._flann_descriptors = flann::Matrix<double>(new double[ioImgInfo._kp.size()*ioImgInfo._descLength],
860                                    ioImgInfo._kp.size(), ioImgInfo._descLength);
861     for (size_t i = 0; i < ioImgInfo._kp.size(); ++i)
862     {
863         memcpy(ioImgInfo._flann_descriptors[i], ioImgInfo._kp[i]->_vec, sizeof(double)*ioImgInfo._descLength);
864     }
865 
866     // build query structure
867     ioImgInfo._flann_index = new flann::Index<flann::L2<double> > (ioImgInfo._flann_descriptors, flann::KDTreeIndexParams(4));
868     ioImgInfo._flann_index->buildIndex();
869 
870     return true;
871 }
872 
FreeMemoryInImage(ImgData & ioImgInfo,const PanoDetector & iPanoDetector)873 bool PanoDetector::FreeMemoryInImage(ImgData& ioImgInfo, const PanoDetector& iPanoDetector)
874 {
875     TRACE_IMG("Freeing memory...");
876 
877     ioImgInfo._ii.clean();
878     ioImgInfo._distancemap.resize(0,0);
879 
880     return true;
881 }
882 
883 
FindMatchesInPair(MatchData & ioMatchData,const PanoDetector & iPanoDetector)884 bool PanoDetector::FindMatchesInPair(MatchData& ioMatchData, const PanoDetector& iPanoDetector)
885 {
886     TRACE_PAIR("Find Matches...");
887 
888     // retrieve the KDTree of image 2
889     flann::Index<flann::L2<double> > * index2 = ioMatchData._i2->_flann_index;
890 
891     // retrieve query points from image 1
892     flann::Matrix<double> & query = ioMatchData._i1->_flann_descriptors;
893 
894     // storage for sorted 2 best matches
895     int nn = 2;
896     flann::Matrix<int> indices(new int[query.rows*nn], query.rows, nn);
897     flann::Matrix<double> dists(new double[query.rows*nn], query.rows, nn);
898 
899     // perform matching using flann
900     index2->knnSearch(query, indices, dists, nn, flann::SearchParams(iPanoDetector.getKDTreeSearchSteps()));
901 
902     //typedef KDTreeSpace::BestMatch<KDElemKeyPoint>		BM_t;
903     //std::set<BM_t, std::greater<BM_t> >	aBestMatches;
904 
905     // store the matches already found to avoid 2 points in image1
906     // match the same point in image2
907     // both matches will be removed.
908     std::set<int> aAlreadyMatched;
909     std::set<int> aBadMatch;
910 
911     // unfiltered vector of matches
912     typedef std::pair<lfeat::KeyPointPtr, int> TmpPair_t;
913     std::vector<TmpPair_t>	aUnfilteredMatches;
914 
915     //PointMatchVector_t aMatches;
916 
917     // go through all the keypoints of image 1
918     for (unsigned aKIt = 0; aKIt < query.rows; ++aKIt)
919     {
920         // accept the match if the second match is far enough
921         // put a lower value for stronger matching default 0.15
922         if (dists[aKIt][0] > iPanoDetector.getKDTreeSecondDistance()  * dists[aKIt][1])
923         {
924             continue;
925         }
926 
927         // check if the kdtree match number is already in the already matched set
928         if (aAlreadyMatched.find(indices[aKIt][0]) != aAlreadyMatched.end())
929         {
930             // add to delete list and continue
931             aBadMatch.insert(indices[aKIt][0]);
932             continue;
933         }
934 
935         // TODO: add check for duplicate matches (can happen if a keypoint gets multiple orientations)
936 
937         // add the match number in already matched set
938         aAlreadyMatched.insert(indices[aKIt][0]);
939 
940         // add the match to the unfiltered list
941         aUnfilteredMatches.push_back(TmpPair_t(ioMatchData._i1->_kp[aKIt], indices[aKIt][0]));
942     }
943 
944     // now filter and fill the vector of matches
945     for (size_t i = 0; i < aUnfilteredMatches.size(); ++i)
946     {
947         TmpPair_t& aP = aUnfilteredMatches[i];
948         // if the image2 match number is in the badmatch set, skip it.
949         if (aBadMatch.find(aP.second) != aBadMatch.end())
950         {
951             continue;
952         }
953 
954         // add the match in the output vector
955         ioMatchData._matches.push_back(lfeat::PointMatchPtr( new lfeat::PointMatch(aP.first, ioMatchData._i2->_kp[aP.second])));
956     }
957 
958     delete[] indices.ptr();
959     delete[] dists.ptr();
960     TRACE_PAIR("Found " << ioMatchData._matches.size() << " matches.");
961     return true;
962 }
963 
RansacMatchesInPair(MatchData & ioMatchData,const PanoDetector & iPanoDetector)964 bool PanoDetector::RansacMatchesInPair(MatchData& ioMatchData, const PanoDetector& iPanoDetector)
965 {
966     // Use panotools model for wide angle lenses
967     HuginBase::RANSACOptimizer::Mode rmode = iPanoDetector._ransacMode;
968     if (rmode == HuginBase::RANSACOptimizer::HOMOGRAPHY ||
969         (rmode == HuginBase::RANSACOptimizer::AUTO && iPanoDetector._panoramaInfo->getImage(ioMatchData._i1->_number).getHFOV() < 65 &&
970              iPanoDetector._panoramaInfo->getImage(ioMatchData._i2->_number).getHFOV() < 65))
971     {
972         return RansacMatchesInPairHomography(ioMatchData, iPanoDetector);
973     }
974     else
975     {
976         return RansacMatchesInPairCam(ioMatchData, iPanoDetector);
977     }
978 }
979 
980 // new code with fisheye aware ransac
RansacMatchesInPairCam(MatchData & ioMatchData,const PanoDetector & iPanoDetector)981 bool PanoDetector::RansacMatchesInPairCam(MatchData& ioMatchData, const PanoDetector& iPanoDetector)
982 {
983     TRACE_PAIR("RANSAC Filtering with Panorama model...");
984 
985     if (ioMatchData._matches.size() < (unsigned int)iPanoDetector.getMinimumMatches())
986     {
987         TRACE_PAIR("Too few matches ... removing all of them.");
988         ioMatchData._matches.clear();
989         return true;
990     }
991 
992     if (ioMatchData._matches.size() < 6)
993     {
994         TRACE_PAIR("Not enough matches for RANSAC filtering.");
995         return true;
996     }
997 
998     // setup a panorama project with the two images.
999     // is this threadsafe (is this read only access?)
1000     HuginBase::UIntSet imgs;
1001     int pano_i1 = ioMatchData._i1->_number;
1002     int pano_i2 = ioMatchData._i2->_number;
1003     imgs.insert(pano_i1);
1004     imgs.insert(pano_i2);
1005     int pano_local_i1 = 0;
1006     int pano_local_i2 = 1;
1007     if (pano_i1 > pano_i2)
1008     {
1009         pano_local_i1 = 1;
1010         pano_local_i2 = 0;
1011     }
1012 
1013     // perform ransac matching.
1014     // ARGH the panotools optimizer uses global variables is not reentrant
1015     std::vector<int> inliers;
1016 #pragma omp critical
1017     {
1018         HuginBase::PanoramaData* panoSubset = iPanoDetector._panoramaInfo->getNewSubset(imgs);
1019 
1020         // create control point vector
1021         HuginBase::CPVector controlPoints(ioMatchData._matches.size());
1022         for (size_t i = 0; i < ioMatchData._matches.size(); ++i)
1023         {
1024             lfeat::PointMatchPtr& aM=ioMatchData._matches[i];
1025             controlPoints[i] = HuginBase::ControlPoint(pano_local_i1, aM->_img1_x, aM->_img1_y,
1026                                             pano_local_i2, aM->_img2_x, aM->_img2_y);
1027         }
1028         panoSubset->setCtrlPoints(controlPoints);
1029 
1030 
1031         PT_setProgressFcn(ptProgress);
1032         PT_setInfoDlgFcn(ptinfoDlg);
1033 
1034         HuginBase::RANSACOptimizer::Mode rmode = iPanoDetector._ransacMode;
1035         if (rmode == HuginBase::RANSACOptimizer::AUTO)
1036         {
1037             rmode = HuginBase::RANSACOptimizer::RPY;
1038         }
1039         // the RANSAC uses the distance in the image for determination of valid parameter
1040         // so make the threshold depending on the image size, use the given pixel distance relative to a 12 MPix image with 4000x3000 pixel
1041         const double threshold = iPanoDetector.getRansacDistanceThreshold() / 5000.0 * hypot(panoSubset->getImage(pano_local_i2).getWidth(), panoSubset->getImage(pano_local_i2).getHeight());
1042         inliers = HuginBase::RANSACOptimizer::findInliers(*panoSubset, pano_local_i1, pano_local_i2,
1043                   threshold, rmode);
1044         PT_setProgressFcn(NULL);
1045         PT_setInfoDlgFcn(NULL);
1046         delete panoSubset;
1047     }
1048 
1049     TRACE_PAIR("Removed " << ioMatchData._matches.size() - inliers.size() << " matches. " << inliers.size() << " remaining.");
1050     if (inliers.size() < 0.5 * ioMatchData._matches.size())
1051     {
1052         // more than 50% of matches were removed, ignore complete pair...
1053         TRACE_PAIR("RANSAC found more than 50% outliers, removing all matches");
1054         ioMatchData._matches.clear();
1055         return true;
1056     }
1057 
1058 
1059     if (inliers.size() < (unsigned int)iPanoDetector.getMinimumMatches())
1060     {
1061         TRACE_PAIR("Too few matches ... removing all of them.");
1062         ioMatchData._matches.clear();
1063         return true;
1064     }
1065 
1066     // keep only inlier matches
1067     lfeat::PointMatchVector_t aInlierMatches;
1068     aInlierMatches.reserve(inliers.size());
1069 
1070     for (size_t i = 0; i < inliers.size(); ++i)
1071     {
1072         aInlierMatches.push_back(ioMatchData._matches[inliers[i]]);
1073     }
1074     ioMatchData._matches = aInlierMatches;
1075 
1076     /*
1077     if (iPanoDetector.getTest())
1078     	TestCode::drawRansacMatches(ioMatchData._i1->_name, ioMatchData._i2->_name, ioMatchData._matches,
1079     								aRemovedMatches, aRansacFilter, iPanoDetector.getDownscale());
1080     */
1081 
1082     return true;
1083 }
1084 
1085 // homography based ransac matching
RansacMatchesInPairHomography(MatchData & ioMatchData,const PanoDetector & iPanoDetector)1086 bool PanoDetector::RansacMatchesInPairHomography(MatchData& ioMatchData, const PanoDetector& iPanoDetector)
1087 {
1088     TRACE_PAIR("RANSAC Filtering...");
1089 
1090     if (ioMatchData._matches.size() < (unsigned int)iPanoDetector.getMinimumMatches())
1091     {
1092         TRACE_PAIR("Too few matches ... removing all of them.");
1093         ioMatchData._matches.clear();
1094         return true;
1095     }
1096 
1097     if (ioMatchData._matches.size() < 6)
1098     {
1099         TRACE_PAIR("Not enough matches for RANSAC filtering.");
1100         return true;
1101     }
1102 
1103     lfeat::PointMatchVector_t aRemovedMatches;
1104 
1105     lfeat::Ransac aRansacFilter;
1106     aRansacFilter.setIterations(iPanoDetector.getRansacIterations());
1107     int thresholdDistance=iPanoDetector.getRansacDistanceThreshold();
1108     //increase RANSAC distance if the image were remapped to not exclude
1109     //too much points in this case
1110     if(ioMatchData._i1->NeedsRemapping() || ioMatchData._i2->NeedsRemapping())
1111     {
1112         thresholdDistance*=5;
1113     }
1114     aRansacFilter.setDistanceThreshold(thresholdDistance);
1115     aRansacFilter.filter(ioMatchData._matches, aRemovedMatches);
1116 
1117 
1118     TRACE_PAIR("Removed " << aRemovedMatches.size() << " matches. " << ioMatchData._matches.size() << " remaining.");
1119 
1120     if (aRemovedMatches.size() > ioMatchData._matches.size())
1121     {
1122         // more than 50% of matches were removed, ignore complete pair...
1123         TRACE_PAIR("More than 50% outliers, removing all matches");
1124         ioMatchData._matches.clear();
1125         return true;
1126     }
1127 
1128     if (iPanoDetector.getTest())
1129         TestCode::drawRansacMatches(ioMatchData._i1->_name, ioMatchData._i2->_name, ioMatchData._matches,
1130                                     aRemovedMatches, aRansacFilter, iPanoDetector.getDownscale());
1131 
1132     return true;
1133 
1134 }
1135 
1136 
FilterMatchesInPair(MatchData & ioMatchData,const PanoDetector & iPanoDetector)1137 bool PanoDetector::FilterMatchesInPair(MatchData& ioMatchData, const PanoDetector& iPanoDetector)
1138 {
1139     TRACE_PAIR("Clustering matches...");
1140 
1141     if (ioMatchData._matches.size() < 2)
1142     {
1143         return true;
1144     }
1145 
1146     // compute min,max of x,y for image1
1147 
1148     double aMinX = std::numeric_limits<double>::max();
1149     double aMinY = std::numeric_limits<double>::max();
1150     double aMaxX = -std::numeric_limits<double>::max();
1151     double aMaxY = -std::numeric_limits<double>::max();
1152 
1153     for (size_t i = 0; i < ioMatchData._matches.size(); ++i)
1154     {
1155         lfeat::PointMatchPtr& aM = ioMatchData._matches[i];
1156         if (aM->_img1_x < aMinX)
1157         {
1158             aMinX = aM->_img1_x;
1159         }
1160         if (aM->_img1_x > aMaxX)
1161         {
1162             aMaxX = aM->_img1_x;
1163         }
1164 
1165         if (aM->_img1_y < aMinY)
1166         {
1167             aMinY = aM->_img1_y;
1168         }
1169         if (aM->_img1_y > aMaxY)
1170         {
1171             aMaxY = aM->_img1_y;
1172         }
1173     }
1174 
1175     double aSizeX = aMaxX - aMinX + 2; // add 2 so max/aSize is strict < 1
1176     double aSizeY = aMaxY - aMinY + 2;
1177 
1178     //
1179 
1180     lfeat::Sieve<lfeat::PointMatchPtr, lfeat::PointMatchPtrSort> aSieve(iPanoDetector.getSieve2Width(),
1181             iPanoDetector.getSieve2Height(),
1182             iPanoDetector.getSieve2Size());
1183 
1184     // insert the points in the Sieve
1185     double aXF = (double)iPanoDetector.getSieve2Width() / aSizeX;
1186     double aYF = (double)iPanoDetector.getSieve2Height() / aSizeY;
1187     for (size_t i = 0; i < ioMatchData._matches.size(); ++i)
1188     {
1189         lfeat::PointMatchPtr& aM = ioMatchData._matches[i];
1190         aSieve.insert(aM, (int)((aM->_img1_x - aMinX) * aXF), (int)((aM->_img1_y - aMinY) * aYF));
1191     }
1192 
1193     // pull remaining values from the sieve
1194     ioMatchData._matches.clear();
1195 
1196     // make an extractor and pull the points
1197     SieveExtractorMatch aSieveExt(ioMatchData._matches);
1198     aSieve.extract(aSieveExt);
1199 
1200     TRACE_PAIR("Kept " << ioMatchData._matches.size() << " matches.");
1201     return true;
1202 }
1203 
writeOutput()1204 void PanoDetector::writeOutput()
1205 {
1206     // Write output pto file
1207 
1208     std::ofstream aOut(_outputFile.c_str(), std::ios_base::trunc);
1209     if( !aOut )
1210     {
1211         std::cerr << "ERROR : "
1212              << "Couldn't open file '" << _outputFile << "'!" << std::endl; //STS
1213         return;
1214     }
1215 
1216     aOut << "# pto project file generated by Hugin's cpfind" << std::endl << std::endl;
1217 
1218     _panoramaInfo->removeDuplicateCtrlPoints();
1219     AppBase::DocumentData::ReadWriteError err = _panoramaInfo->writeData(aOut);
1220     if (err != AppBase::DocumentData::SUCCESSFUL)
1221     {
1222         std::cerr << "ERROR couldn't write to output file '" << _outputFile << "'!" << std::endl;
1223         return;
1224     }
1225 }
1226 
writeKeyfile(ImgData & imgInfo)1227 void PanoDetector::writeKeyfile(ImgData& imgInfo)
1228 {
1229     // Write output keyfile
1230 
1231     std::ofstream aOut(imgInfo._keyfilename.c_str(), std::ios_base::trunc);
1232 
1233     lfeat::SIFTFormatWriter writer(aOut);
1234 
1235     int origImgWidth =  _panoramaInfo->getImage(imgInfo._number).getSize().width();
1236     int origImgHeight =  _panoramaInfo->getImage(imgInfo._number).getSize().height();
1237 
1238     lfeat::ImageInfo img_info(imgInfo._name, origImgWidth, origImgHeight);
1239 
1240     writer.writeHeader ( img_info, imgInfo._kp.size(), imgInfo._descLength );
1241 
1242     for(size_t i=0; i<imgInfo._kp.size(); ++i)
1243     {
1244         lfeat::KeyPointPtr& aK=imgInfo._kp[i];
1245         writer.writeKeypoint ( aK->_x, aK->_y, aK->_scale, aK->_ori, aK->_score,
1246                                imgInfo._descLength, aK->_vec );
1247     }
1248     writer.writeFooter();
1249 }
1250 
1251