Main MRPT website > C++ reference for MRPT 1.5.6
CFeatureExtraction_SIFT.cpp
Go to the documentation of this file.
1 /* +---------------------------------------------------------------------------+
2  | Mobile Robot Programming Toolkit (MRPT) |
3  | http://www.mrpt.org/ |
4  | |
5  | Copyright (c) 2005-2017, Individual contributors, see AUTHORS file |
6  | See: http://www.mrpt.org/Authors - All rights reserved. |
7  | Released under BSD License. See details in http://www.mrpt.org/License |
8  +---------------------------------------------------------------------------+ */
9 
10 #include "vision-precomp.h" // Precompiled headers
11 
12 #include <mrpt/system/threads.h>
13 #include <mrpt/system/os.h>
15 
16 #if MRPT_HAS_OPENCV && MRPT_HAS_SIFT_HESS
17  #include "sift-hess/sift.h"
18  #include "sift-hess/imgfeatures.h"
19  #include "sift-hess/utils.h"
20 #endif
21 
22 
23 // Universal include for all versions of OpenCV
24 #include <mrpt/otherlibs/do_opencv_includes.h>
25 #ifdef HAVE_OPENCV_NONFREE //MRPT_HAS_OPENCV_NONFREE
26 # include <opencv2/nonfree/nonfree.hpp>
27 #endif
28 #ifdef HAVE_OPENCV_FEATURES2D
29 # include <opencv2/features2d/features2d.hpp>
30 #endif
31 #ifdef HAVE_OPENCV_XFEATURES2D
32 # include <opencv2/xfeatures2d.hpp>
33 #endif
34 
35 
36 // TODO: Remove, it's just for GetTempPathA
37 #ifdef MRPT_OS_WINDOWS
38  #include <windows.h>
39 #endif
40 
41 using namespace mrpt;
42 using namespace mrpt::vision;
43 using namespace mrpt::system;
44 using namespace mrpt::utils;
45 using namespace mrpt::math;
46 using namespace std;
47 
48 #if MRPT_HAS_OPENCV
49 # if MRPT_OPENCV_VERSION_NUM>=0x211
50  using namespace cv;
51 # endif
52 #endif
53 
54 /************************* Local Function Prototypes for Hess' SIFT *************************/
55 #if MRPT_HAS_OPENCV && MRPT_HAS_SIFT_HESS
56 extern "C" // This is mandatory since the implementations are in ".c" files
57 {
58 IplImage* create_init_img(const IplImage*, int, double );
59 IplImage*** build_gauss_pyr( IplImage*, int, int, double );
60 IplImage*** build_dog_pyr( IplImage***, int, int );
61 CvSeq* scale_space_extrema( IplImage***, int, int, double, int, CvMemStorage*);
62 int is_extremum( IplImage***, int, int, int, int );
63 struct feature* new_feature( void );
64 void calc_feature_scales( CvSeq*, double, int );
65 void adjust_for_img_dbl( CvSeq* );
66 void calc_feature_oris( CvSeq*, IplImage*** );
67 void compute_descriptors( CvSeq*, IplImage***, int, int );
68 int feature_cmp( void*, void*, void* );
69 void release_pyr( IplImage****, int, int );
70 }
71 #endif // MRPT_HAS_OPENCV
72 
73 
74 /************************************************************************************************
75 * extractFeaturesSIFT *
76 ************************************************************************************************/
78  const CImage &img,
79  CFeatureList &feats,
80  unsigned int init_ID,
81  unsigned int nDesiredFeatures,
82  const TImageROI &ROI) const
83 {
84  bool usingROI = false;
85  if( ROI.xMin != 0 || ROI.xMax != 0 || ROI.yMin != 0 || ROI.yMax != 0 )
86  usingROI = true; // A ROI has been defined
87 
88  // ROI can not be managed properly (yet) with these method, so we extract a subimage
89 
90  // use a smart pointer so we just copy the pointer if the image is grayscale, or we'll create a new one if it was RGB:
91  CImage img_grayscale(img, FAST_REF_OR_CONVERT_TO_GRAY); // Was: auxImgPtr;
92  if( usingROI )
93  {
94  ASSERT_( ROI.xMin >= 0 && ROI.xMin < ROI.xMax && ROI.xMax < img.getWidth() && ROI.yMin >= 0 && ROI.yMax < img.getHeight() && ROI.yMin < ROI.yMax );
95  CImage auximg;
96  img_grayscale.extract_patch( auximg, ROI.xMin, ROI.yMin, ROI.xMax-ROI.xMin+1, ROI.yMax-ROI.yMin+1 ); // Subimage in "auxImg"
97  img_grayscale.swap(auximg);
98  }
99 
100  switch( options.SIFTOptions.implementation )
101  {
102 // --------------------------------------------------------------------------------------
103 // Binary in C# -> OPTIONAL: Feature position already computed
104 // --------------------------------------------------------------------------------------
105 
106  case CSBinary:
107  {
108 #ifdef MRPT_OS_WINDOWS
109 
110  char filImg[2000],filOut[2000],filFeat[2000];
111  char paramImg[2000];
112 
113  GetTempPathA(1000,filOut); os::strcat(filOut,1000,"temp_out.txt"); // OUTPUT FILE
114  GetTempPathA(1000,filImg); os::strcat(filImg,1000,"temp_img.bmp"); // INPUT IMAGE (BMP) FOR BINARY IN (C#)
115 
116  bool onlyDesc = feats.size() > 0 ? true : false;
117 
118  if( onlyDesc )
119  {
120  GetTempPathA(1000,filFeat); os::strcat(filFeat,1000,"temp_feats.txt"); // KEYPOINTS INPUT FILE
121  CMatrix listPoints(feats.size(),2);
122  for (size_t i= 0;i<feats.size();i++)
123  {
124  listPoints(i,0) = feats[i]->x;
125  listPoints(i,1) = feats[i]->y;
126  }
127  listPoints.saveToTextFile( filFeat, MATRIX_FORMAT_FIXED /*Float format*/ );
128  } // end if
129  // -------------------------------------------
130  // CALL TO "extractSIFT.exe"
131  // -------------------------------------------
132  img_grayscale.saveToFile( filImg );
133 
134  // ------------------------------------
135  // Version with "CreateProcess":
136  // ------------------------------------
137  os::strcpy(paramImg,1000,"extractSIFT.exe -i"); os::strcat(paramImg,1000,filImg);
138  os::strcat(paramImg,1000," -f"); os::strcat(paramImg,1000,filOut);
139  os::strcat(paramImg,1000," -l"); os::strcat(paramImg,1000,filFeat);
140 
141  // ------------------------------------
142  // Launch process
143  // ------------------------------------
144  bool ret = mrpt::system::launchProcess( paramImg );
145 
146  if( !ret )
147  THROW_EXCEPTION( "[extractFeaturesSIFT] Could not launch external process... (extractSIFT.exe)" )
148 
149  // Process Results
150  CFeatureList::iterator itFeat = feats.begin();
151  size_t nFeats;
152 
153  CMatrix aux;
154  aux.loadFromTextFile( filOut );
155  std::cout << "[computeSiftFeatures] " << aux.getRowCount() << " features." << std::endl;
156 
157  if( onlyDesc )
158  nFeats = feats.size();
159  else
160  {
161  nFeats = aux.getRowCount();
162  feats.resize( nFeats );
163  }
164 
165  for( size_t i = 0;
166  itFeat != feats.end();
167  i++, itFeat++)
168  {
169  (*itFeat)->type = featSIFT;
170  (*itFeat)->x = usingROI ? aux(i,0) + ROI.xMin : aux(i,0);
171  (*itFeat)->y = usingROI ? aux(i,1) + ROI.yMin : aux(i,1);
172  (*itFeat)->orientation = aux(i,2);
173  (*itFeat)->scale = aux(i,3);
174  (*itFeat)->ID = init_ID + i;
175 
176  // The descriptor:
177  aux.extractRow(i, (*itFeat)->descriptors.SIFT, 4);
178  }
179  remove(filImg);
180  remove(filOut);
181 #else
182  THROW_EXCEPTION("Unfortunately, this SIFT Implementation only runs in Windows OS, try Hess implementation");
183 #endif
184  break;
185  } // end case Binary in C#
186  case VedaldiBinary:
187  {
188 
189  // --------------------------------------------------------------------------------------
190  // Binary by Vedaldi: NOT IMPLEMENTED YET. Input in PGM format
191  // --------------------------------------------------------------------------------------
192 #ifdef MRPT_OS_WINDOWS
193  THROW_EXCEPTION("Usage of Vedaldi Binary not implemented yet, please, try another one");
194 #else
195  THROW_EXCEPTION("Unfortunately, this SIFT Implementation only runs in Windows OS, try Hess implementation");
196 #endif
197  break;
198  } // end case Binary by Vedaldi
199 // --------------------------------------------------------------------------------------
200 // Binary by David Lowe
201 // --------------------------------------------------------------------------------------
202  case LoweBinary: // Binary by Lowe
203  {
204 
205 #ifdef MRPT_OS_WINDOWS
206  char filImg[2000],filOut[2000];
207  char paramImg[2000];
208 
209  feats.clear();
210 
211  GetTempPathA(1000,filOut); os::strcat(filOut,1000,"temp_out.txt"); // OUTPUT FILE
212  GetTempPathA(1000,filImg); os::strcat(filImg,1000,"temp_img.pgm"); // INPUT IMAGE (PGM) FOR ORIGINAL BINARY BY LOWE
213 
214  bool valid = img_grayscale.saveToFile( filImg );
215  if(!valid)
216  THROW_EXCEPTION( "An error occurred when saving input image into a .pgm file");
217 
218  // CONVERT TO UNCOMPRESSED RAW PGM (TODO: Solve in a better way this issue)
219  os::strcpy( paramImg,1000, format( "cmd /C gmic.exe %s -o %s -quiet", filImg, filImg ).c_str() );
220 
221  bool ret = mrpt::system::launchProcess( paramImg );
222 
223  if(!ret)
224  THROW_EXCEPTION("[extractFeaturesSIFT] Could not launch external process... (gmic.exe)");
225 
226  // ------------------------------------
227  // Version with "CreateProcess":
228  // ------------------------------------
229  os::strcpy(paramImg,1000,"cmd /C siftWin32.exe <"); os::strcat(paramImg,1000,filImg);
230  os::strcat(paramImg,1000," >"); os::strcat(paramImg,1000,filOut);
231 
232  ret = mrpt::system::launchProcess( paramImg );
233 
234  if(!ret)
235  THROW_EXCEPTION("[extractFeaturesSIFT] Could not launch external process... (siftWin32.exe)");
236 
237  // ------------------------------------
238  // Process Results
239  // ------------------------------------
240  unsigned int dLen, nFeats;
241  FILE *f = os::fopen( filOut, "rt");
242  if(!f)
243  THROW_EXCEPTION( "Error in extract SIFT with Lowe binary, output file not found!" );
244  fscanf( f,"%u %u", &nFeats, &dLen); // Number of feats and length of the descriptor
245 
246  for( size_t i = 0; i < nFeats; i++ )
247  {
248  CFeaturePtr feat = CFeature::Create();
249 
250  feat->type = featSIFT; // Type
251  feat->ID = init_ID + i; // Identifier
252 
253  // Position, orientation and scale
254  // IMPORTANTE NOTE: Lowe format stores first the 'y' coordinate and then the 'x' one
255  float fx,fy,fo,fs;
256  fscanf( f, "%f %f %f %f", &fy, &fx, &fo, &fs );
257 
258  feat->x = usingROI ? fx + ROI.xMin : fx;
259  feat->y = usingROI ? fy + ROI.yMin : fy;
260  feat->orientation = fo;
261  feat->scale = fs;
262 
263  // The descriptor
264  feat->descriptors.SIFT.resize( dLen );
265  unsigned int c;
266  for(unsigned int k = 0; k < dLen; k++)
267  {
268  fscanf( f, "%u", &c );
269  feat->descriptors.SIFT[k] = (unsigned char)c;
270  }
271 
272  feats.push_back( feat );
273  } // end for
274  os::fclose( f );
275  remove(filImg);
276  remove(filOut);
277 #else
278  THROW_EXCEPTION("Unfortunately, this SIFT Implementation only runs in Windows OS, try Hess implementation");
279 #endif
280  break;
281  } // end case Binary by Lowe
282 // --------------------------------------------------------------------------------------
283 // Hess implementation
284 // --------------------------------------------------------------------------------------
285  case Hess: // Implementation by Robert Hess
286  {
287 
288 #if !MRPT_HAS_SIFT_HESS
289  THROW_EXCEPTION("Method not available since MRPT has been compiled without Hess' SIFT library")
290 #elif MRPT_HAS_OPENCV // OK, we have Hess' sift:
291  IplImage* init_img;
292  IplImage*** gauss_pyr, *** dog_pyr;
293  CvMemStorage* storage;
294  CvSeq* features;
295  int octvs;
296  /* check arguments */
297  ASSERT_(img_grayscale.getWidth() != 0 && img_grayscale.getHeight() != 0);
298  /* build scale space pyramid; smallest dimension of top level is ~4 pixels */
299  const IplImage* ipl_im = img_grayscale.getAs<IplImage>();
300  init_img = create_init_img( ipl_im, SIFT_IMG_DBL, SIFT_SIGMA );
301  octvs = log( (float)(MIN( init_img->width, init_img->height )) ) / log((float)2) - 2;
302  gauss_pyr = build_gauss_pyr( init_img, octvs, SIFT_INTVLS, SIFT_SIGMA );
303  dog_pyr = build_dog_pyr( gauss_pyr, octvs, SIFT_INTVLS );
304  storage = cvCreateMemStorage( 0 );
305  features = scale_space_extrema( dog_pyr, octvs, SIFT_INTVLS,
306  options.SIFTOptions.threshold, // SIFT_CONTR_THR,
307  options.SIFTOptions.edgeThreshold, // SIFT_CURV_THR
308  storage );
309  calc_feature_scales( features, SIFT_SIGMA, SIFT_INTVLS );
310  if( SIFT_IMG_DBL )
311  adjust_for_img_dbl( features );
312  calc_feature_oris( features, gauss_pyr );
313  compute_descriptors( features, gauss_pyr, SIFT_DESCR_WIDTH, SIFT_DESCR_HIST_BINS );
314 
315  /* sort features by decreasing scale and move from CvSeq to array */
316  cvSeqSort( features, (CvCmpFunc)feature_cmp, NULL );
317 
318  /* get only the desired features */
319  if( nDesiredFeatures > 0 )
320  {
321  if( nDesiredFeatures < (unsigned int)features->total )
322  cvSeqPopMulti( features, NULL, features->total - nDesiredFeatures );
323  else
324  cout << "[Warning] Detected less features than the requested " << features->total << " vs " << nDesiredFeatures << endl;
325  } // end if
326 
327  /* convert CvSeq into a FeatureList */
328  convertCvSeqInCFeatureList( features, feats, init_ID, ROI );
329 
330  // clear Hess-features
331  cvClearSeq( features );
332  cvReleaseMemStorage( &storage );
333  cvReleaseImage( &init_img );
334  release_pyr( &gauss_pyr, octvs, SIFT_INTVLS + 3 );
335  release_pyr( &dog_pyr, octvs, SIFT_INTVLS + 2 );
336 #else
337  THROW_EXCEPTION("Method not available since MRPT has been compiled without OpenCV")
338 #endif //MRPT_HAS_OPENCV
339  break;
340  } // end case Hess
341 //***********************************************************************************************
342 // USING OPENCV
343 //***********************************************************************************************
344  case OpenCV:
345  {
346 #if defined(HAVE_OPENCV_NONFREE) || defined(HAVE_OPENCV_XFEATURES2D) //MRPT_HAS_OPENCV_NONFREE
347 
348  #if MRPT_OPENCV_VERSION_NUM >= 0x211 && MRPT_OPENCV_VERSION_NUM < 0x300
349 
350  SiftFeatureDetector SIFTDetector(
351  options.SIFTOptions.threshold, //SIFT::DetectorParams::GET_DEFAULT_THRESHOLD(),
352  options.SIFTOptions.edgeThreshold //SIFT::DetectorParams::GET_DEFAULT_EDGE_THRESHOLD() );
353  );
354 
355  SiftDescriptorExtractor SIFTDescriptor;
356 
357  vector<KeyPoint> cv_feats; // The OpenCV output feature list
358 
359 
360  const IplImage* cGrey = img_grayscale.getAs<IplImage>();
361 
362  Mat theImg = cvarrToMat( cGrey );
363  SIFTDetector.detect( theImg, cv_feats );
364 
365  Mat desc;
366  SIFTDescriptor.compute( theImg, cv_feats, desc );
367 
368  //fromOpenCVToMRPT( theImg, cv_feats, desc, nDesiredFeatures, outList );
369  const size_t N = cv_feats.size();
370  unsigned int nMax = nDesiredFeatures != 0 && N > nDesiredFeatures ? nDesiredFeatures : N;
371  const int offset = (int)this->options.patchSize/2 + 1;
372  const size_t size_2 = options.patchSize/2;
373  const size_t imgH = img.getHeight();
374  const size_t imgW = img.getWidth();
375  unsigned int i = 0;
376  unsigned int cont = 0;
377  TFeatureID nextID = init_ID;
378  feats.clear();
379  while( cont != nMax && i != N )
380  {
381  const int xBorderInf = (int)floor( cv_feats[i].pt.x - size_2 );
382  const int xBorderSup = (int)floor( cv_feats[i].pt.x + size_2 );
383  const int yBorderInf = (int)floor( cv_feats[i].pt.y - size_2 );
384  const int yBorderSup = (int)floor( cv_feats[i].pt.y + size_2 );
385 
386  if( options.patchSize==0 || ( (xBorderSup < (int)imgW) && (xBorderInf > 0) && (yBorderSup < (int)imgH) && (yBorderInf > 0) ) )
387  {
388  CFeaturePtr ft = CFeature::Create();
389  ft->type = featSIFT;
390  ft->ID = nextID++;
391  ft->x = cv_feats[i].pt.x;
392  ft->y = cv_feats[i].pt.y;
393  ft->response = cv_feats[i].response;
394  ft->orientation = cv_feats[i].angle;
395  ft->scale = cv_feats[i].size;
396  ft->patchSize = options.patchSize; // The size of the feature patch
397  ft->descriptors.SIFT.resize( 128 );
398  memcpy( &(ft->descriptors.SIFT[0]), &desc.data[128*i], 128*sizeof(ft->descriptors.SIFT[0]) ); // The descriptor
399 
400  if( options.patchSize > 0 )
401  {
402  img.extract_patch(
403  ft->patch,
404  round( ft->x ) - offset,
405  round( ft->y ) - offset,
406  options.patchSize,
407  options.patchSize ); // Image patch surronding the feature
408  }
409  feats.push_back( ft );
410  ++cont;
411  }
412  ++i;
413  }
414  feats.resize( cont );
415  #else
416  // MRPT_OPENCV_VERSION_NUM >= 0x300
417  using namespace cv;
418  vector<KeyPoint> cv_feats;
419 
420  cv::Ptr<cv::xfeatures2d::SIFT>sift = cv::xfeatures2d::SIFT::create(nDesiredFeatures,3, options.SIFTOptions.threshold, options.SIFTOptions.edgeThreshold,1.6 ); //gb
421  const IplImage* cGrey = img_grayscale.getAs<IplImage>();
422  Mat theImg = cvarrToMat(cGrey);
423  //SIFTDetector.detect(theImg, cv_feats);
424  sift->detect(theImg, cv_feats); //gb
425  Mat desc;
426  //SIFTDescriptor.compute(theImg, cv_feats, desc);
427  sift->compute(theImg, cv_feats, desc);
428 
429  //fromOpenCVToMRPT( theImg, cv_feats, desc, nDesiredFeatures, outList );
430  const size_t N = cv_feats.size();
431  unsigned int nMax = nDesiredFeatures != 0 && N > nDesiredFeatures ? nDesiredFeatures : N;
432  const int offset = (int)this->options.patchSize / 2 + 1;
433  const size_t size_2 = options.patchSize / 2;
434  const size_t imgH = img.getHeight();
435  const size_t imgW = img.getWidth();
436  unsigned int i = 0;
437  unsigned int cont = 0;
438  TFeatureID nextID = init_ID;
439  feats.clear();
440 
441 
442  while (cont != nMax && i != N)
443  {
444  const int xBorderInf = (int)floor(cv_feats[i].pt.x - size_2);
445  const int xBorderSup = (int)floor(cv_feats[i].pt.x + size_2);
446  const int yBorderInf = (int)floor(cv_feats[i].pt.y - size_2);
447  const int yBorderSup = (int)floor(cv_feats[i].pt.y + size_2);
448 
449  if (options.patchSize == 0 || ((xBorderSup < (int)imgW) && (xBorderInf > 0) && (yBorderSup < (int)imgH) && (yBorderInf > 0)))
450  {
451  CFeaturePtr ft = CFeature::Create();
452  ft->type = featSIFT;
453  ft->ID = nextID++;
454  ft->x = cv_feats[i].pt.x;
455  ft->y = cv_feats[i].pt.y;
456  ft->response = cv_feats[i].response;
457  ft->orientation = cv_feats[i].angle;
458  ft->scale = cv_feats[i].size;
459  ft->patchSize = options.patchSize; // The size of the feature patch
460  ft->descriptors.SIFT.resize(128);
461  memcpy(&(ft->descriptors.SIFT[0]), &desc.data[128 * i], 128 * sizeof(ft->descriptors.SIFT[0])); // The descriptor
462 
463  if (options.patchSize > 0)
464  {
465  img.extract_patch(
466  ft->patch,
467  round(ft->x) - offset,
468  round(ft->y) - offset,
469  options.patchSize,
470  options.patchSize); // Image patch surronding the feature
471  }
472  feats.push_back(ft);
473  ++cont;
474  }
475  ++i;
476  }
477  feats.resize(cont);
478  #endif
479 #else
480  THROW_EXCEPTION("This method requires OpenCV >= 2.1.1 with nonfree module")
481 #endif
482  break;
483  } // end case OpenCV
484  return;
485  default:{break;} // end default
486  } // end switch
487 } // end extractFeaturesSIFT
488 
489 
490 /************************************************************************************************
491 * computeSiftDescriptors *
492 ************************************************************************************************/
493 // Compute SIFT descriptors on a set of already localized points
495  CFeatureList &in_features) const
496 {
497 
498  ASSERT_( in_features.size() > 0 );
499  switch( options.SIFTOptions.implementation )
500  {
501  case CSBinary:
502  {
503 #ifdef MRPT_OS_WINDOWS
504  char filImg[2000],filOut[2000],filFeat[2000];
505  char paramImg[2000];
506 
508 
509  // Save to temporary file
510  GetTempPathA(1000,filImg); os::strcat( filImg, 1000, "temp_img.bmp" );
511  GetTempPathA(1000,filOut); os::strcat( filOut,1000, "temp_feats.txt" );
512  GetTempPathA(1000,filFeat); os::strcat( filFeat, 1000, "temp_KLT_feats.txt" );
513 
514  // Fill the input file
515  FILE *fout = os::fopen( filFeat,"wt");
516  for(feat = in_features.begin(); feat != in_features.end(); feat++)
517  os::fprintf( fout, "%.6f %.6f\n", (*feat)->x, (*feat)->y );
518 
519  os::fclose(fout);
520 
521  // -------------------------------------------
522  // CALL TO "extractSIFT.exe"
523  // -------------------------------------------
524  in_img.saveToFile( filImg );
525 
526  // Version with "CreateProcess":
527  // ------------------------------------
528  os::strcpy(paramImg,1000,"extractSIFT.exe -i"); os::strcat(paramImg,1000,filImg);
529  os::strcat(paramImg,1000," -f"); os::strcat(paramImg,1000,filOut);
530  os::strcat(paramImg,1000," -l"); os::strcat(paramImg,1000,filFeat);
531 
532  bool ret = mrpt::system::launchProcess( paramImg );
533 
534  if(!ret)
535  THROW_EXCEPTION("[extractFeaturesSIFT] Could not launch external process... (extractSIFT.exe)");
536 
537  MRPT_START
538  // Load the results:
539  CMatrix aux;
540  aux.loadFromTextFile( filOut );
541  size_t nRows = aux.getRowCount();
542 
543  std::cout << "[computeSiftFeatures1] " << nRows << " features.\n";
544 
545  unsigned int i;
546  float lx,ly;
547  lx = 0.0;
548  ly = 0.0;
549  feat = in_features.begin();
550 
551  for(i = 0; i < nRows; i++)
552  {
553  if( aux(i,0) != lx || aux(i,1) != ly ) // Only one descriptor for feature
554  {
555  (*feat)->type = featSIFT;
556  (*feat)->orientation = aux(i,2);
557  (*feat)->scale = aux(i,3);
558 
559  // The descriptor:
560  aux.extractRow(i, (*feat)->descriptors.SIFT, 4);
561 
562  lx = aux(i,0);
563  ly = aux(i,1);
564 
565  feat++;
566  } // end if
567  } // end for
568  MRPT_END
569 
570  break;
571 #else
572  THROW_EXCEPTION("TO DO @ linux OS!");
573 #endif
574  } // end case CSBinary
575  case Hess: // Implementation by Hess
576  {
577 #if !MRPT_HAS_SIFT_HESS
578  THROW_EXCEPTION("Method not available since MRPT has been compiled without Hess' SIFT library")
579 #elif MRPT_HAS_OPENCV // OK, we have Hess' sift:
580  IplImage* init_img;
581  IplImage*** gauss_pyr, *** dog_pyr;
582  CvMemStorage* storage;
583  CvSeq* features;
584  int octvs; //, n = 0;
585 
586  /* check arguments */
587  ASSERT_(in_img.getWidth() != 0 && in_img.getHeight() != 0);
588 
589  /* build scale space pyramid; smallest dimension of top level is ~4 pixels */
590  const CImage img_grayscale(in_img, FAST_REF_OR_CONVERT_TO_GRAY);
591  const IplImage* ipl_im = img_grayscale.getAs<IplImage>();
592 
593  init_img = create_init_img( ipl_im, SIFT_IMG_DBL, SIFT_SIGMA );
594  octvs = log( (float)(MIN( init_img->width, init_img->height )) ) / log((float)2) - 2;
595  gauss_pyr = build_gauss_pyr( init_img, octvs, SIFT_INTVLS, SIFT_SIGMA );
596  dog_pyr = build_dog_pyr( gauss_pyr, octvs, SIFT_INTVLS );
597 
598  storage = cvCreateMemStorage( 0 );
599  features = static_cast<CvSeq*>(my_scale_space_extrema( in_features, dog_pyr, octvs, SIFT_INTVLS,
600  options.SIFTOptions.threshold, // SIFT_CONTR_THR,
601  options.SIFTOptions.edgeThreshold, // SIFT_CURV_THR
602  storage ));
603  calc_feature_scales( features, SIFT_SIGMA, SIFT_INTVLS );
604  if( SIFT_IMG_DBL )
605  my_adjust_for_img_dbl( features );
606  calc_feature_oris( features, gauss_pyr );
607  compute_descriptors( features, gauss_pyr, SIFT_DESCR_WIDTH, SIFT_DESCR_HIST_BINS );
608 
609  // merge Hess-features and MRPT-features
610  insertCvSeqInCFeatureList( features, in_features );
611 
612  // clear Hess-features
613  cvClearSeq( features );
614 
615  // Free memory
616  cvReleaseMemStorage( &storage );
617  cvReleaseImage( &init_img );
618  release_pyr( &gauss_pyr, octvs, SIFT_INTVLS + 3 );
619  release_pyr( &dog_pyr, octvs, SIFT_INTVLS + 2 );
620 
621  /* sort features by decreasing scale and move from CvSeq to array */
622  //cvSeqSort( features, (CvCmpFunc)feature_cmp, NULL );
623  //n = features->total;
624 
625  //struct feature** feat;
626  //*feat = (feature*)calloc( n, sizeof(struct feature) );
627  //*feat = (feature*)cvCvtSeqToArray( features, *feat, CV_WHOLE_SEQ );
628  //for( i = 0; i < n; i++ )
629  //{
630  // free( (*feat)[i].feature_data );
631  // (*feat)[i].feature_data = NULL;
632  //}
633 
634  //// The number of features must be the same
635  //cout << "NUMBER OF INPUT FEATURES: " << in_features.size() << endl;
636  //cout << "NUMBER OF OUTPUT FEATURES: " << n << endl;
637 
638  //// Transform to CFeatureList:
639  //CFeatureList::iterator itFeat;
640  //int m;
641  //for( itFeat = in_features.begin(), m = 0; itFeat != in_features.end(); itFeat++, m++)
642  //{
643  // (*itFeat)->type = featSIFT;
644  // (*itFeat)->x = (*feat)[m].x;
645  // (*itFeat)->y = (*feat)[m].y;
646  // (*itFeat)->orientation = (*feat)[m].ori;
647  // (*itFeat)->scale = (*feat)[m].scl;
648  // for( unsigned int k = 0; k < FEATURE_MAX_D; k++)
649  // (*itFeat)->descriptors.SIFT[k] = (*feat)[m].descr[k];
650  // (*itFeat)->hasDescriptor = true;
651  //} // end for features
652 #else
653  THROW_EXCEPTION("Method not available since MRPT has been compiled without OpenCV")
654 #endif //MRPT_HAS_OPENCV
655  break;
656  } // end case Hess
657  default:
658  {
659  cout << "SIFT Extraction method not supported for features with already known image coordinates" << endl;
660  break;
661  }
662  } // end switch
663 
664 } // end computeSiftDescriptors
665 
666 
667 // ------------------------------------------------------------------------------------
668 // convertCvSeqInCFeatureList
669 // ------------------------------------------------------------------------------------
670 void CFeatureExtraction::convertCvSeqInCFeatureList( void* features_, CFeatureList &list, unsigned int init_ID, const TImageROI &ROI ) const
671 {
672 #if MRPT_HAS_OPENCV && MRPT_HAS_SIFT_HESS
673  CvSeq* features = reinterpret_cast<CvSeq*>( features_ );
674 
675  // Is there a defined ROI?
676  bool usingROI = false;
677  if( ROI.xMin != 0 || ROI.xMax != 0 || ROI.yMin != 0 || ROI.yMax != 0 )
678  usingROI = true;
679 
680  int n = features->total;
681  ASSERT_(n > 0);
682 
683  list.clear();
684  struct feature* thisFeat;
685  for( int k = 0; k < n; k++ )
686  {
687  thisFeat = (feature*)cvGetSeqElem( features, k );
688  CFeaturePtr feat = CFeature::Create();
689  feat->ID = (TFeatureID)(k + init_ID);
690  feat->x = usingROI ? thisFeat->x + ROI.xMin : thisFeat->x;
691  feat->y = usingROI ? thisFeat->y + ROI.yMin : thisFeat->y;
692  feat->type = featSIFT;
693  feat->orientation = thisFeat->ori;
694  feat->scale = thisFeat->scl;
695  feat->descriptors.SIFT.resize( thisFeat->d );
696  for( int i = 0; i < thisFeat->d; i++ )
697  feat->descriptors.SIFT[i] = (unsigned char)thisFeat->descr[i];
698 
699  list.push_back(feat);
700  } // end for
701 #else
702  THROW_EXCEPTION("Method not available since MRPT has been compiled without OpenCV")
703 #endif //MRPT_HAS_OPENCV
704 }
705 // ------------------------------------------------------------------------------------
706 // insertCvSeqInCFeatureList
707 // ------------------------------------------------------------------------------------
708 void CFeatureExtraction::insertCvSeqInCFeatureList( void* features_, CFeatureList &list, unsigned int init_ID ) const
709 {
710 #if MRPT_HAS_OPENCV && MRPT_HAS_SIFT_HESS
711  CvSeq* features = reinterpret_cast<CvSeq*>( features_ );
712 
713  int n = features->total;
714  ASSERT_(n > 0);
715 
716  CFeatureList::iterator itFeat;
717  struct feature* thisFeat;
718  int k;
719  for( itFeat = list.begin(), k = 0; itFeat != list.end() && k < n; k++ )
720  {
721  thisFeat = (feature*)cvGetSeqElem( features, k );
722  if( (*itFeat)->x == thisFeat->x && (*itFeat)->y == thisFeat->y )
723  {
724  (*itFeat)->ID = (TFeatureID)(k + init_ID);
725  (*itFeat)->orientation = thisFeat->ori;
726  (*itFeat)->scale = thisFeat->scl;
727  (*itFeat)->descriptors.SIFT.resize( thisFeat->d );
728  for( int i = 0; i < thisFeat->d; i++ )
729  (*itFeat)->descriptors.SIFT[i] = (unsigned char)thisFeat->descr[i];
730 
731  itFeat++;
732  } // end if
733  } // end for
734 #else
735  THROW_EXCEPTION("Method not available since MRPT has been compiled without OpenCV")
736 #endif //MRPT_HAS_OPENCV
737 }
738 
739 /*
740 Halves feature coordinates and scale in case the input image was doubled
741 prior to scale space construction.
742 
743 @param features array of features
744 */
745 void CFeatureExtraction::my_adjust_for_img_dbl( void* features_ ) const
746 {
747 #if MRPT_HAS_OPENCV && MRPT_HAS_SIFT_HESS
748  CvSeq* features = reinterpret_cast<CvSeq*>( features_ );
749  struct feature* feat;
750  int i, n;
751 
752  n = features->total;
753  for( i = 0; i < n; i++ )
754  {
755  feat = CV_GET_SEQ_ELEM( struct feature, features, i );
756  feat->scl /= 2.0;
757  }
758 #else
759  THROW_EXCEPTION("Method not available since MRPT has been compiled without OpenCV")
760 #endif //MRPT_HAS_OPENCV
761 }
762 
763 // ------------------------------------------------------------------------------------
764 // my_scale_space_extrema
765 // ------------------------------------------------------------------------------------
767  CFeatureList &featList, void* dog_pyr_,
768  int octvs, int intvls, double contr_thr, int curv_thr,
769  void* storage_ ) const
770 {
771  MRPT_UNUSED_PARAM(contr_thr); MRPT_UNUSED_PARAM(curv_thr);
772 #if MRPT_HAS_OPENCV && MRPT_HAS_SIFT_HESS
773  CvMemStorage* storage = reinterpret_cast<CvMemStorage*>( storage_ );
774  IplImage*** dog_pyr = reinterpret_cast<IplImage***>( dog_pyr_ );
775 
776  CvSeq* features;
777 // double prelim_contr_thr = 0.5 * contr_thr / intvls;
778  struct feature* feat;
779  struct detection_data* ddata;
780  int o, i;
781 
782  CFeatureList::iterator itFeats;
783 
784  features = cvCreateSeq( 0, sizeof(CvSeq), sizeof(struct feature), storage );
785  for(itFeats = featList.begin(); itFeats != featList.end(); itFeats++)
786  {
787  //unsigned int nMax, nMin;
788  double gExt = -1e5;
789  int go = 0, gi = 1;
790  float gr = 0.0, gc = 0.0;
791 
792  if( (*itFeats)->y < 0 || (*itFeats)->x < 0 )
793  continue; // Do not process this feature
794 
795  // Find the best octave and interval where the feature is the 'best' extrema
796  for( o = 0; o < octvs; o++ )
797  {
798  float r = (*itFeats)->y/pow(2.0,o); // Dog pyramid coordinates
799  float c = (*itFeats)->x/pow(2.0,o);
800 
801  for( i = 1; i <= intvls; i++ )
802  {
803  double s = SIFT_SIGMA * pow(2.0, o+i/intvls);
804  double val = s*getLaplacianValue( dog_pyr, o, i, r, c );
805  if( val > gExt )
806  {
807  gExt = val;
808  go = o;
809  gi = i;
810  gr = r;
811  gc = c;
812  }
813  //getTimesExtrema( dog_pyr, o, i, r, c, nMax, nMin ); // Get number of times the point is bigger or lower than the surroundings
814  //if( nMax > gExt || nMin > gExt )
815  //{
816  // gExt = max(nMax,nMin);
817  // go = o;
818  // gi = i;
819  // gr = r;
820  // gc = c;
821  //} // end if
822  } // end for intervals
823  } // end for octaves
824 
825  feat = new_feature(); // Create the new feature
826  ddata = feat_detection_data( feat ); // Feature data
827  feat->img_pt.x = feat->x = (*itFeats)->x; // Feature Coordinates in the original image
828  feat->img_pt.y = feat->y = (*itFeats)->y;
829  ddata->r = (int)gr; // Feature Coordinates in the proper octave
830  ddata->c = (int)gc;
831  ddata->octv = go; // Octave
832  ddata->intvl = gi; // Interval
833  ddata->subintvl = 0.0; // Subinterval (not used)
834 
835  cvSeqPush( features, feat ); // Insert into the feature sequence
836  } // end for features in the list
837 
838  return features;
839 #else
840  THROW_EXCEPTION("Method not available since MRPT has been compiled without OpenCV")
841 #endif //MRPT_HAS_OPENCV
842 }
843 
844 /************************************************************************************************
845 * getLaplaceValue *
846 ************************************************************************************************/
847 double CFeatureExtraction::getLaplacianValue( void* dog_pyr_, int octv, int intvl, float row, float col ) const
848 {
849 #if MRPT_HAS_OPENCV && MRPT_HAS_SIFT_HESS
850  IplImage*** dog_pyr = reinterpret_cast<IplImage***>( dog_pyr_ );
851  double value = 0.0;
852  int j, k;
853 
854  for( j = -1; j <= 1; j++ )
855  for( k = -1; k <= 1; k++ )
856  {
857  // Apply convolution mask for computing the Laplacian:
858  // -1 -1 -1
859  // -1 +8 -1
860  // -1 -1 -1
861  if( k == 0 && j == 0 )
862  value += 8*pixval32f( dog_pyr[octv][intvl], (int)(row + j), (int)(col + k) );
863  else
864  value -= pixval32f( dog_pyr[octv][intvl], (int)(row + j), (int)(col + k) );
865  }
866  return value;
867 #else
868  THROW_EXCEPTION("Method not available since MRPT has been compiled without OpenCV")
869 #endif //MRPT_HAS_OPENCV
870 } // end getLaplacianValue
871 
872 /************************************************************************************************
873 * getTimesExtrema *
874 ************************************************************************************************/
875 //void CFeatureExtraction::getTimesExtrema( IplImage*** dog_pyr, int octv, int intvl, float row, float col, unsigned int &nMin, unsigned int &nMax )
876 void CFeatureExtraction::getTimesExtrema( void* dog_pyr_, int octv, int intvl, float row, float col, unsigned int &nMin, unsigned int &nMax ) const
877 {
878 #if MRPT_HAS_OPENCV && MRPT_HAS_SIFT_HESS
879  IplImage*** dog_pyr = reinterpret_cast<IplImage***>( dog_pyr_ );
880 
881  float val = pixval32f( dog_pyr[octv][intvl], (int)row, (int)col );
882  int i, j, k;
883 
884  nMax = 0;
885  nMin = 0;
886  /* check for extrema */
887  //cout << "VAL: " << pixval32f( dog_pyr[octv][intvl], row, col );
888  for( i = -1; i <= 1; i++ )
889  for( j = -1; j <= 1; j++ )
890  for( k = -1; k <= 1; k++ )
891  {
892  //cout << pixval32f( dog_pyr[octv][intvl+i], row + j, col + k ) << endl;
893  if( val > pixval32f( dog_pyr[octv][intvl+i], (int)(row + j), (int)(col + k) ) )
894  nMax++;
895  if( val < pixval32f( dog_pyr[octv][intvl+i], (int)(row + j), (int)(col + k) ) )
896  nMin++;
897  }
898 #else
899  THROW_EXCEPTION("Method not available since MRPT has been compiled without OpenCV")
900 #endif //MRPT_HAS_OPENCV
901 } // end getTimesExtrema
902 
void BASE_IMPEXP memcpy(void *dest, size_t destSize, const void *src, size_t copyCount) MRPT_NO_THROWS
An OS and compiler independent version of "memcpy".
Definition: os.cpp:358
FILE BASE_IMPEXP * fopen(const char *fileName, const char *mode) MRPT_NO_THROWS
An OS-independent version of fopen.
Definition: os.cpp:255
double scl
scale of a Lowe-style feature
Definition: imgfeatures.h:59
Classes for serialization, sockets, ini-file manipulation, streams, list of properties-values, timewatch, extensions to STL.
Definition: zip.h:16
void internal_computeSiftDescriptors(const mrpt::utils::CImage &in_img, CFeatureList &in_features) const
Compute the SIFT descriptor of the provided features into the input image.
int octv
Definition: sift.h:28
double subintvl
Definition: sift.h:30
void getTimesExtrema(void *dog_pyr, int octvs, int intvls, float row, float col, unsigned int &nMin, unsigned int &nMax) const
Gets the number of times that a point in the image is higher or lower than the surroundings in the im...
This namespace provides a OS-independent interface to many useful functions: filenames manipulation...
Definition: math_frwds.h:29
int BASE_IMPEXP void BASE_IMPEXP fclose(FILE *f)
An OS-independent version of fclose.
Definition: os.cpp:272
A class for storing images as grayscale or RGB bitmaps.
Definition: CImage.h:101
static __inline float pixval32f(IplImage *img, int r, int c)
A function to get a pixel value from a 32-bit floating-point image.
double x
x coord
Definition: imgfeatures.h:54
size_t size() const
Definition: CFeature.h:280
#define THROW_EXCEPTION(msg)
#define SIFT_DESCR_HIST_BINS
default number of bins per histogram in descriptor array
Definition: sift.h:58
GLenum GLsizei n
Definition: glext.h:4618
GLintptr offset
Definition: glext.h:3780
CvPoint2D64f img_pt
location in image
Definition: imgfeatures.h:68
double ori
orientation of a Lowe-style feature
Definition: imgfeatures.h:60
#define MIN(a, b)
Definition: jpegint.h:266
int BASE_IMPEXP fprintf(FILE *fil, const char *format,...) MRPT_NO_THROWS MRPT_printf_format_check(2
An OS-independent version of fprintf.
Definition: os.cpp:412
STL namespace.
char BASE_IMPEXP * strcpy(char *dest, size_t destSize, const char *source) MRPT_NO_THROWS
An OS-independent version of strcpy.
Definition: os.cpp:296
GLdouble s
Definition: glext.h:3602
bool BASE_IMPEXP launchProcess(const std::string &command)
Executes the given command (which may contain a program + arguments), and waits until it finishes...
Definition: threads.cpp:454
const T * getAs() const
Returns a pointer to a const T* containing the image - the idea is to call like "img.getAs<IplImage>()" so we can avoid here including OpenCV&#39;s headers.
Definition: CImage.h:517
float yMin
Y coordinate limits [0,imageHeight)
int intvl
Definition: sift.h:29
A structure for defining a ROI within an image.
void extractFeaturesSIFT(const mrpt::utils::CImage &img, CFeatureList &feats, unsigned int init_ID=0, unsigned int nDesiredFeatures=0, const TImageROI &ROI=TImageROI()) const
Extract features from the image based on the SIFT method.
int d
descriptor length
Definition: imgfeatures.h:61
This base provides a set of functions for maths stuff.
Definition: CArrayNumeric.h:19
void resize(size_t N)
Definition: CFeature.h:283
#define MRPT_END
#define MRPT_UNUSED_PARAM(a)
Can be used to avoid "not used parameters" warnings from the compiler.
double getLaplacianValue(void *dog_pyr, int octvs, int intvls, float row, float col) const
Computes the Laplacian value of the feature in the corresponing image in the pyramid.
const GLubyte * c
Definition: glext.h:5590
GLint GLvoid * img
Definition: glext.h:3645
Scale Invariant Feature Transform [LOWE&#39;04].
int val
Definition: mrpt_jpeglib.h:953
std::string BASE_IMPEXP format(const char *fmt,...) MRPT_printf_format_check(1
A std::string version of C sprintf.
Definition: format.cpp:21
Classes for computer vision, detectors, features, etc.
#define SIFT_DESCR_WIDTH
default width of descriptor histogram array
Definition: sift.h:55
void insertCvSeqInCFeatureList(void *features, CFeatureList &list, unsigned int init_ID=0) const
Append a sequence of openCV features into an MRPT feature list.
void push_back(const CFeaturePtr &f)
Definition: CFeature.h:285
uint64_t TFeatureID
Definition of a feature ID.
holds feature data relevant to detection
Definition: sift.h:24
void swap(CImage &o)
Very efficient swap of two images (just swap the internal pointers)
Definition: CImage.cpp:135
A list of visual features, to be used as output by detectors, as input/output by trackers, etc.
Definition: CFeature.h:211
#define feat_detection_data(f)
Definition: sift.h:94
double descr[FEATURE_MAX_D]
descriptor
Definition: imgfeatures.h:62
float xMin
X coordinate limits [0,imageWidth)
void my_adjust_for_img_dbl(void *features) const
Adjust scale if the image was initially doubled.
void saveToTextFile(const std::string &fileName, bool APPEND=false)
Save feature list to a text file.
Definition: CFeature.cpp:865
#define MRPT_START
void * my_scale_space_extrema(CFeatureList &featList, void *dog_pyr, int octvs, int intvls, double contr_thr, int curv_thr, void *storage) const
Computes extrema in the scale space.
This is the global namespace for all Mobile Robot Programming Toolkit (MRPT) libraries.
GLdouble GLdouble GLdouble r
Definition: glext.h:3618
Structure to represent an affine invariant image feature.
Definition: imgfeatures.h:52
#define SIFT_INTVLS
default number of sampled intervals per octave
Definition: sift.h:40
GLenum GLenum GLvoid * row
Definition: glext.h:3533
void extract_patch(CImage &patch, const unsigned int col=0, const unsigned int row=0, const unsigned int width=1, const unsigned int height=1) const
Extract a patch from this image, saveing it into "patch" (its previous contents will be overwritten)...
Definition: CImage.cpp:1368
#define ASSERT_(f)
char BASE_IMPEXP * strcat(char *dest, size_t destSize, const char *source) MRPT_NO_THROWS
An OS-independent version of strcat.
Definition: os.cpp:281
int round(const T value)
Returns the closer integer (int) to x.
Definition: round.h:26
bool saveToFile(const std::string &fileName, int jpeg_quality=95) const
Save the image to a file, whose format is determined from the extension (internally uses OpenCV)...
Definition: CImage.cpp:299
TInternalFeatList::iterator iterator
Definition: CFeature.h:261
static CFeaturePtr Create()
GLsizei const GLfloat * value
Definition: glext.h:3929
#define SIFT_SIGMA
default sigma for initial gaussian smoothing
Definition: sift.h:43
#define SIFT_IMG_DBL
default threshold on keypoint contrast |D(x)|
Definition: sift.h:52
double y
y coord
Definition: imgfeatures.h:55
This class is a "CSerializable" wrapper for "CMatrixFloat".
Definition: CMatrix.h:30
size_t getWidth() const MRPT_OVERRIDE
Returns the width of the image in pixels.
Definition: CImage.cpp:855
void convertCvSeqInCFeatureList(void *features, CFeatureList &list, unsigned int init_ID=0, const TImageROI &ROI=TImageROI()) const
Converts a sequence of openCV features into an MRPT feature list.
fixed floating point &#39;f&#39;
Definition: math_frwds.h:66
size_t getHeight() const MRPT_OVERRIDE
Returns the height of the image in pixels.
Definition: CImage.cpp:884



Page generated by Doxygen 1.8.14 for MRPT 1.5.6 Git: 4c65e8431 Tue Apr 24 08:18:17 2018 +0200 at lun oct 28 01:35:26 CET 2019