sift.dispatch.cpp 23 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622
  1. // This file is part of OpenCV project.
  2. // It is subject to the license terms in the LICENSE file found in the top-level directory
  3. // of this distribution and at http://opencv.org/license.html.
  4. //
  5. // Copyright (c) 2006-2010, Rob Hess <hess@eecs.oregonstate.edu>
  6. // Copyright (C) 2009, Willow Garage Inc., all rights reserved.
  7. // Copyright (C) 2020, Intel Corporation, all rights reserved.
  8. /**********************************************************************************************\
  9. Implementation of SIFT is based on the code from http://blogs.oregonstate.edu/hess/code/sift/
  10. Below is the original copyright.
  11. Patent US6711293 expired in March 2020.
  12. // Copyright (c) 2006-2010, Rob Hess <hess@eecs.oregonstate.edu>
  13. // All rights reserved.
  14. // The following patent has been issued for methods embodied in this
  15. // software: "Method and apparatus for identifying scale invariant features
  16. // in an image and use of same for locating an object in an image," David
  17. // G. Lowe, US Patent 6,711,293 (March 23, 2004). Provisional application
  18. // filed March 8, 1999. Asignee: The University of British Columbia. For
  19. // further details, contact David Lowe (lowe@cs.ubc.ca) or the
  20. // University-Industry Liaison Office of the University of British
  21. // Columbia.
  22. // Note that restrictions imposed by this patent (and possibly others)
  23. // exist independently of and may be in conflict with the freedoms granted
  24. // in this license, which refers to copyright of the program, not patents
  25. // for any methods that it implements. Both copyright and patent law must
  26. // be obeyed to legally use and redistribute this program and it is not the
  27. // purpose of this license to induce you to infringe any patents or other
  28. // property right claims or to contest validity of any such claims. If you
  29. // redistribute or use the program, then this license merely protects you
  30. // from committing copyright infringement. It does not protect you from
  31. // committing patent infringement. So, before you do anything with this
  32. // program, make sure that you have permission to do so not merely in terms
  33. // of copyright, but also in terms of patent law.
  34. // Please note that this license is not to be understood as a guarantee
  35. // either. If you use the program according to this license, but in
  36. // conflict with patent law, it does not mean that the licensor will refund
  37. // you for any losses that you incur if you are sued for your patent
  38. // infringement.
  39. // Redistribution and use in source and binary forms, with or without
  40. // modification, are permitted provided that the following conditions are
  41. // met:
  42. // * Redistributions of source code must retain the above copyright and
  43. // patent notices, this list of conditions and the following
  44. // disclaimer.
  45. // * Redistributions in binary form must reproduce the above copyright
  46. // notice, this list of conditions and the following disclaimer in
  47. // the documentation and/or other materials provided with the
  48. // distribution.
  49. // * Neither the name of Oregon State University nor the names of its
  50. // contributors may be used to endorse or promote products derived
  51. // from this software without specific prior written permission.
  52. // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
  53. // IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
  54. // TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
  55. // PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
  56. // HOLDER BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
  57. // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
  58. // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
  59. // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
  60. // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
  61. // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  62. // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  63. \**********************************************************************************************/
  64. #include "precomp.hpp"
  65. #include <opencv2/core/hal/hal.hpp>
  66. #include <opencv2/core/utils/tls.hpp>
  67. #include <opencv2/core/utils/logger.hpp>
  68. #include "sift.simd.hpp"
  69. #include "sift.simd_declarations.hpp" // defines CV_CPU_DISPATCH_MODES_ALL=AVX2,...,BASELINE based on CMakeLists.txt content
  70. namespace cv {
  71. /*!
  72. SIFT implementation.
  73. The class implements SIFT algorithm by D. Lowe.
  74. */
  75. class SIFT_Impl : public SIFT
  76. {
  77. public:
  78. explicit SIFT_Impl( int nfeatures = 0, int nOctaveLayers = 3,
  79. double contrastThreshold = 0.04, double edgeThreshold = 10,
  80. double sigma = 1.6, int descriptorType = CV_32F,
  81. bool enable_precise_upscale = true );
  82. //! returns the descriptor size in floats (128)
  83. int descriptorSize() const CV_OVERRIDE;
  84. //! returns the descriptor type
  85. int descriptorType() const CV_OVERRIDE;
  86. //! returns the default norm type
  87. int defaultNorm() const CV_OVERRIDE;
  88. //! finds the keypoints and computes descriptors for them using SIFT algorithm.
  89. //! Optionally it can compute descriptors for the user-provided keypoints
  90. void detectAndCompute(InputArray img, InputArray mask,
  91. std::vector<KeyPoint>& keypoints,
  92. OutputArray descriptors,
  93. bool useProvidedKeypoints = false) CV_OVERRIDE;
  94. void buildGaussianPyramid( const Mat& base, std::vector<Mat>& pyr, int nOctaves ) const;
  95. void buildDoGPyramid( const std::vector<Mat>& pyr, std::vector<Mat>& dogpyr ) const;
  96. void findScaleSpaceExtrema( const std::vector<Mat>& gauss_pyr, const std::vector<Mat>& dog_pyr,
  97. std::vector<KeyPoint>& keypoints ) const;
  98. void read( const FileNode& fn) CV_OVERRIDE;
  99. void write( FileStorage& fs) const CV_OVERRIDE;
  100. void setNFeatures(int maxFeatures) CV_OVERRIDE { nfeatures = maxFeatures; }
  101. int getNFeatures() const CV_OVERRIDE { return nfeatures; }
  102. void setNOctaveLayers(int nOctaveLayers_) CV_OVERRIDE { nOctaveLayers = nOctaveLayers_; }
  103. int getNOctaveLayers() const CV_OVERRIDE { return nOctaveLayers; }
  104. void setContrastThreshold(double contrastThreshold_) CV_OVERRIDE { contrastThreshold = contrastThreshold_; }
  105. double getContrastThreshold() const CV_OVERRIDE { return contrastThreshold; }
  106. void setEdgeThreshold(double edgeThreshold_) CV_OVERRIDE { edgeThreshold = edgeThreshold_; }
  107. double getEdgeThreshold() const CV_OVERRIDE { return edgeThreshold; }
  108. void setSigma(double sigma_) CV_OVERRIDE { sigma = sigma_; }
  109. double getSigma() const CV_OVERRIDE { return sigma; }
  110. protected:
  111. CV_PROP_RW int nfeatures;
  112. CV_PROP_RW int nOctaveLayers;
  113. CV_PROP_RW double contrastThreshold;
  114. CV_PROP_RW double edgeThreshold;
  115. CV_PROP_RW double sigma;
  116. CV_PROP_RW int descriptor_type;
  117. CV_PROP_RW bool enable_precise_upscale;
  118. };
  119. Ptr<SIFT> SIFT::create( int _nfeatures, int _nOctaveLayers,
  120. double _contrastThreshold, double _edgeThreshold, double _sigma, bool enable_precise_upscale )
  121. {
  122. CV_TRACE_FUNCTION();
  123. return makePtr<SIFT_Impl>(_nfeatures, _nOctaveLayers, _contrastThreshold, _edgeThreshold, _sigma, CV_32F, enable_precise_upscale);
  124. }
  125. Ptr<SIFT> SIFT::create( int _nfeatures, int _nOctaveLayers,
  126. double _contrastThreshold, double _edgeThreshold, double _sigma, int _descriptorType, bool enable_precise_upscale )
  127. {
  128. CV_TRACE_FUNCTION();
  129. // SIFT descriptor supports 32bit floating point and 8bit unsigned int.
  130. CV_Assert(_descriptorType == CV_32F || _descriptorType == CV_8U);
  131. return makePtr<SIFT_Impl>(_nfeatures, _nOctaveLayers, _contrastThreshold, _edgeThreshold, _sigma, _descriptorType, enable_precise_upscale);
  132. }
  133. String SIFT::getDefaultName() const
  134. {
  135. return (Feature2D::getDefaultName() + ".SIFT");
  136. }
  137. static inline void
  138. unpackOctave(const KeyPoint& kpt, int& octave, int& layer, float& scale)
  139. {
  140. octave = kpt.octave & 255;
  141. layer = (kpt.octave >> 8) & 255;
  142. octave = octave < 128 ? octave : (-128 | octave);
  143. scale = octave >= 0 ? 1.f/(1 << octave) : (float)(1 << -octave);
  144. }
  145. static Mat createInitialImage( const Mat& img, bool doubleImageSize, float sigma, bool enable_precise_upscale )
  146. {
  147. CV_TRACE_FUNCTION();
  148. Mat gray, gray_fpt;
  149. if( img.channels() == 3 || img.channels() == 4 )
  150. {
  151. cvtColor(img, gray, COLOR_BGR2GRAY);
  152. gray.convertTo(gray_fpt, DataType<sift_wt>::type, SIFT_FIXPT_SCALE, 0);
  153. }
  154. else
  155. img.convertTo(gray_fpt, DataType<sift_wt>::type, SIFT_FIXPT_SCALE, 0);
  156. float sig_diff;
  157. if( doubleImageSize )
  158. {
  159. sig_diff = sqrtf( std::max(sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA * 4, 0.01f) );
  160. Mat dbl;
  161. if (enable_precise_upscale) {
  162. dbl.create(Size(gray_fpt.cols*2, gray_fpt.rows*2), gray_fpt.type());
  163. Mat H = Mat::zeros(2, 3, CV_32F);
  164. H.at<float>(0, 0) = 0.5f;
  165. H.at<float>(1, 1) = 0.5f;
  166. cv::warpAffine(gray_fpt, dbl, H, dbl.size(), INTER_LINEAR | WARP_INVERSE_MAP, BORDER_REFLECT);
  167. } else {
  168. #if DoG_TYPE_SHORT
  169. resize(gray_fpt, dbl, Size(gray_fpt.cols*2, gray_fpt.rows*2), 0, 0, INTER_LINEAR_EXACT);
  170. #else
  171. resize(gray_fpt, dbl, Size(gray_fpt.cols*2, gray_fpt.rows*2), 0, 0, INTER_LINEAR);
  172. #endif
  173. }
  174. Mat result;
  175. GaussianBlur(dbl, result, Size(), sig_diff, sig_diff);
  176. return result;
  177. }
  178. else
  179. {
  180. sig_diff = sqrtf( std::max(sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA, 0.01f) );
  181. Mat result;
  182. GaussianBlur(gray_fpt, result, Size(), sig_diff, sig_diff);
  183. return result;
  184. }
  185. }
  186. void SIFT_Impl::buildGaussianPyramid( const Mat& base, std::vector<Mat>& pyr, int nOctaves ) const
  187. {
  188. CV_TRACE_FUNCTION();
  189. std::vector<double> sig(nOctaveLayers + 3);
  190. pyr.resize(nOctaves*(nOctaveLayers + 3));
  191. // precompute Gaussian sigmas using the following formula:
  192. // \sigma_{total}^2 = \sigma_{i}^2 + \sigma_{i-1}^2
  193. sig[0] = sigma;
  194. double k = std::pow( 2., 1. / nOctaveLayers );
  195. for( int i = 1; i < nOctaveLayers + 3; i++ )
  196. {
  197. double sig_prev = std::pow(k, (double)(i-1))*sigma;
  198. double sig_total = sig_prev*k;
  199. sig[i] = std::sqrt(sig_total*sig_total - sig_prev*sig_prev);
  200. }
  201. for( int o = 0; o < nOctaves; o++ )
  202. {
  203. for( int i = 0; i < nOctaveLayers + 3; i++ )
  204. {
  205. Mat& dst = pyr[o*(nOctaveLayers + 3) + i];
  206. if( o == 0 && i == 0 )
  207. dst = base;
  208. // base of new octave is halved image from end of previous octave
  209. else if( i == 0 )
  210. {
  211. const Mat& src = pyr[(o-1)*(nOctaveLayers + 3) + nOctaveLayers];
  212. resize(src, dst, Size(src.cols/2, src.rows/2),
  213. 0, 0, INTER_NEAREST);
  214. }
  215. else
  216. {
  217. const Mat& src = pyr[o*(nOctaveLayers + 3) + i-1];
  218. GaussianBlur(src, dst, Size(), sig[i], sig[i]);
  219. }
  220. }
  221. }
  222. }
  223. class buildDoGPyramidComputer : public ParallelLoopBody
  224. {
  225. public:
  226. buildDoGPyramidComputer(
  227. int _nOctaveLayers,
  228. const std::vector<Mat>& _gpyr,
  229. std::vector<Mat>& _dogpyr)
  230. : nOctaveLayers(_nOctaveLayers),
  231. gpyr(_gpyr),
  232. dogpyr(_dogpyr) { }
  233. void operator()( const cv::Range& range ) const CV_OVERRIDE
  234. {
  235. CV_TRACE_FUNCTION();
  236. const int begin = range.start;
  237. const int end = range.end;
  238. for( int a = begin; a < end; a++ )
  239. {
  240. const int o = a / (nOctaveLayers + 2);
  241. const int i = a % (nOctaveLayers + 2);
  242. const Mat& src1 = gpyr[o*(nOctaveLayers + 3) + i];
  243. const Mat& src2 = gpyr[o*(nOctaveLayers + 3) + i + 1];
  244. Mat& dst = dogpyr[o*(nOctaveLayers + 2) + i];
  245. subtract(src2, src1, dst, noArray(), DataType<sift_wt>::type);
  246. }
  247. }
  248. private:
  249. int nOctaveLayers;
  250. const std::vector<Mat>& gpyr;
  251. std::vector<Mat>& dogpyr;
  252. };
  253. void SIFT_Impl::buildDoGPyramid( const std::vector<Mat>& gpyr, std::vector<Mat>& dogpyr ) const
  254. {
  255. CV_TRACE_FUNCTION();
  256. int nOctaves = (int)gpyr.size()/(nOctaveLayers + 3);
  257. dogpyr.resize( nOctaves*(nOctaveLayers + 2) );
  258. parallel_for_(Range(0, nOctaves * (nOctaveLayers + 2)), buildDoGPyramidComputer(nOctaveLayers, gpyr, dogpyr));
  259. }
  260. class findScaleSpaceExtremaComputer : public ParallelLoopBody
  261. {
  262. public:
  263. findScaleSpaceExtremaComputer(
  264. int _o,
  265. int _i,
  266. int _threshold,
  267. int _idx,
  268. int _step,
  269. int _cols,
  270. int _nOctaveLayers,
  271. double _contrastThreshold,
  272. double _edgeThreshold,
  273. double _sigma,
  274. const std::vector<Mat>& _gauss_pyr,
  275. const std::vector<Mat>& _dog_pyr,
  276. TLSData<std::vector<KeyPoint> > &_tls_kpts_struct)
  277. : o(_o),
  278. i(_i),
  279. threshold(_threshold),
  280. idx(_idx),
  281. step(_step),
  282. cols(_cols),
  283. nOctaveLayers(_nOctaveLayers),
  284. contrastThreshold(_contrastThreshold),
  285. edgeThreshold(_edgeThreshold),
  286. sigma(_sigma),
  287. gauss_pyr(_gauss_pyr),
  288. dog_pyr(_dog_pyr),
  289. tls_kpts_struct(_tls_kpts_struct) { }
  290. void operator()( const cv::Range& range ) const CV_OVERRIDE
  291. {
  292. CV_TRACE_FUNCTION();
  293. std::vector<KeyPoint>& kpts = tls_kpts_struct.getRef();
  294. CV_CPU_DISPATCH(findScaleSpaceExtrema, (o, i, threshold, idx, step, cols, nOctaveLayers, contrastThreshold, edgeThreshold, sigma, gauss_pyr, dog_pyr, kpts, range),
  295. CV_CPU_DISPATCH_MODES_ALL);
  296. }
  297. private:
  298. int o, i;
  299. int threshold;
  300. int idx, step, cols;
  301. int nOctaveLayers;
  302. double contrastThreshold;
  303. double edgeThreshold;
  304. double sigma;
  305. const std::vector<Mat>& gauss_pyr;
  306. const std::vector<Mat>& dog_pyr;
  307. TLSData<std::vector<KeyPoint> > &tls_kpts_struct;
  308. };
  309. //
  310. // Detects features at extrema in DoG scale space. Bad features are discarded
  311. // based on contrast and ratio of principal curvatures.
  312. void SIFT_Impl::findScaleSpaceExtrema( const std::vector<Mat>& gauss_pyr, const std::vector<Mat>& dog_pyr,
  313. std::vector<KeyPoint>& keypoints ) const
  314. {
  315. CV_TRACE_FUNCTION();
  316. const int nOctaves = (int)gauss_pyr.size()/(nOctaveLayers + 3);
  317. const int threshold = cvFloor(0.5 * contrastThreshold / nOctaveLayers * 255 * SIFT_FIXPT_SCALE);
  318. keypoints.clear();
  319. TLSDataAccumulator<std::vector<KeyPoint> > tls_kpts_struct;
  320. for( int o = 0; o < nOctaves; o++ )
  321. for( int i = 1; i <= nOctaveLayers; i++ )
  322. {
  323. const int idx = o*(nOctaveLayers+2)+i;
  324. const Mat& img = dog_pyr[idx];
  325. const int step = (int)img.step1();
  326. const int rows = img.rows, cols = img.cols;
  327. parallel_for_(Range(SIFT_IMG_BORDER, rows-SIFT_IMG_BORDER),
  328. findScaleSpaceExtremaComputer(
  329. o, i, threshold, idx, step, cols,
  330. nOctaveLayers,
  331. contrastThreshold,
  332. edgeThreshold,
  333. sigma,
  334. gauss_pyr, dog_pyr, tls_kpts_struct));
  335. }
  336. std::vector<std::vector<KeyPoint>*> kpt_vecs;
  337. tls_kpts_struct.gather(kpt_vecs);
  338. for (size_t i = 0; i < kpt_vecs.size(); ++i) {
  339. keypoints.insert(keypoints.end(), kpt_vecs[i]->begin(), kpt_vecs[i]->end());
  340. }
  341. }
  342. static
  343. void calcSIFTDescriptor(
  344. const Mat& img, Point2f ptf, float ori, float scl,
  345. int d, int n, Mat& dst, int row
  346. )
  347. {
  348. CV_TRACE_FUNCTION();
  349. CV_CPU_DISPATCH(calcSIFTDescriptor, (img, ptf, ori, scl, d, n, dst, row),
  350. CV_CPU_DISPATCH_MODES_ALL);
  351. }
  352. class calcDescriptorsComputer : public ParallelLoopBody
  353. {
  354. public:
  355. calcDescriptorsComputer(const std::vector<Mat>& _gpyr,
  356. const std::vector<KeyPoint>& _keypoints,
  357. Mat& _descriptors,
  358. int _nOctaveLayers,
  359. int _firstOctave)
  360. : gpyr(_gpyr),
  361. keypoints(_keypoints),
  362. descriptors(_descriptors),
  363. nOctaveLayers(_nOctaveLayers),
  364. firstOctave(_firstOctave) { }
  365. void operator()( const cv::Range& range ) const CV_OVERRIDE
  366. {
  367. CV_TRACE_FUNCTION();
  368. const int begin = range.start;
  369. const int end = range.end;
  370. static const int d = SIFT_DESCR_WIDTH, n = SIFT_DESCR_HIST_BINS;
  371. for ( int i = begin; i<end; i++ )
  372. {
  373. KeyPoint kpt = keypoints[i];
  374. int octave, layer;
  375. float scale;
  376. unpackOctave(kpt, octave, layer, scale);
  377. CV_Assert(octave >= firstOctave && layer <= nOctaveLayers+2);
  378. float size=kpt.size*scale;
  379. Point2f ptf(kpt.pt.x*scale, kpt.pt.y*scale);
  380. const Mat& img = gpyr[(octave - firstOctave)*(nOctaveLayers + 3) + layer];
  381. float angle = 360.f - kpt.angle;
  382. if(std::abs(angle - 360.f) < FLT_EPSILON)
  383. angle = 0.f;
  384. calcSIFTDescriptor(img, ptf, angle, size*0.5f, d, n, descriptors, i);
  385. }
  386. }
  387. private:
  388. const std::vector<Mat>& gpyr;
  389. const std::vector<KeyPoint>& keypoints;
  390. Mat& descriptors;
  391. int nOctaveLayers;
  392. int firstOctave;
  393. };
  394. static void calcDescriptors(const std::vector<Mat>& gpyr, const std::vector<KeyPoint>& keypoints,
  395. Mat& descriptors, int nOctaveLayers, int firstOctave )
  396. {
  397. CV_TRACE_FUNCTION();
  398. parallel_for_(Range(0, static_cast<int>(keypoints.size())), calcDescriptorsComputer(gpyr, keypoints, descriptors, nOctaveLayers, firstOctave));
  399. }
  400. //////////////////////////////////////////////////////////////////////////////////////////
  401. SIFT_Impl::SIFT_Impl( int _nfeatures, int _nOctaveLayers,
  402. double _contrastThreshold, double _edgeThreshold, double _sigma, int _descriptorType, bool _enable_precise_upscale)
  403. : nfeatures(_nfeatures), nOctaveLayers(_nOctaveLayers),
  404. contrastThreshold(_contrastThreshold), edgeThreshold(_edgeThreshold), sigma(_sigma), descriptor_type(_descriptorType),
  405. enable_precise_upscale(_enable_precise_upscale)
  406. {
  407. if (!enable_precise_upscale) {
  408. CV_LOG_ONCE_INFO(NULL, "precise upscale disabled, this is now deprecated as it was found to induce a location bias");
  409. }
  410. }
  411. int SIFT_Impl::descriptorSize() const
  412. {
  413. return SIFT_DESCR_WIDTH*SIFT_DESCR_WIDTH*SIFT_DESCR_HIST_BINS;
  414. }
  415. int SIFT_Impl::descriptorType() const
  416. {
  417. return descriptor_type;
  418. }
  419. int SIFT_Impl::defaultNorm() const
  420. {
  421. return NORM_L2;
  422. }
  423. void SIFT_Impl::detectAndCompute(InputArray _image, InputArray _mask,
  424. std::vector<KeyPoint>& keypoints,
  425. OutputArray _descriptors,
  426. bool useProvidedKeypoints)
  427. {
  428. CV_TRACE_FUNCTION();
  429. int firstOctave = -1, actualNOctaves = 0, actualNLayers = 0;
  430. Mat image = _image.getMat(), mask = _mask.getMat();
  431. if( image.empty() || image.depth() != CV_8U )
  432. CV_Error( Error::StsBadArg, "image is empty or has incorrect depth (!=CV_8U)" );
  433. if( !mask.empty() && mask.type() != CV_8UC1 )
  434. CV_Error( Error::StsBadArg, "mask has incorrect type (!=CV_8UC1)" );
  435. if( useProvidedKeypoints )
  436. {
  437. firstOctave = 0;
  438. int maxOctave = INT_MIN;
  439. for( size_t i = 0; i < keypoints.size(); i++ )
  440. {
  441. int octave, layer;
  442. float scale;
  443. unpackOctave(keypoints[i], octave, layer, scale);
  444. firstOctave = std::min(firstOctave, octave);
  445. maxOctave = std::max(maxOctave, octave);
  446. actualNLayers = std::max(actualNLayers, layer-2);
  447. }
  448. firstOctave = std::min(firstOctave, 0);
  449. CV_Assert( firstOctave >= -1 && actualNLayers <= nOctaveLayers );
  450. actualNOctaves = maxOctave - firstOctave + 1;
  451. }
  452. Mat base = createInitialImage(image, firstOctave < 0, (float)sigma, enable_precise_upscale);
  453. std::vector<Mat> gpyr;
  454. int nOctaves = actualNOctaves > 0 ? actualNOctaves : cvRound(std::log( (double)std::min( base.cols, base.rows ) ) / std::log(2.) - 2) - firstOctave;
  455. //double t, tf = getTickFrequency();
  456. //t = (double)getTickCount();
  457. buildGaussianPyramid(base, gpyr, nOctaves);
  458. //t = (double)getTickCount() - t;
  459. //printf("pyramid construction time: %g\n", t*1000./tf);
  460. if( !useProvidedKeypoints )
  461. {
  462. std::vector<Mat> dogpyr;
  463. buildDoGPyramid(gpyr, dogpyr);
  464. //t = (double)getTickCount();
  465. findScaleSpaceExtrema(gpyr, dogpyr, keypoints);
  466. KeyPointsFilter::removeDuplicatedSorted( keypoints );
  467. if( nfeatures > 0 )
  468. KeyPointsFilter::retainBest(keypoints, nfeatures);
  469. //t = (double)getTickCount() - t;
  470. //printf("keypoint detection time: %g\n", t*1000./tf);
  471. if( firstOctave < 0 )
  472. for( size_t i = 0; i < keypoints.size(); i++ )
  473. {
  474. KeyPoint& kpt = keypoints[i];
  475. float scale = 1.f/(float)(1 << -firstOctave);
  476. kpt.octave = (kpt.octave & ~255) | ((kpt.octave + firstOctave) & 255);
  477. kpt.pt *= scale;
  478. kpt.size *= scale;
  479. }
  480. if( !mask.empty() )
  481. KeyPointsFilter::runByPixelsMask( keypoints, mask );
  482. }
  483. else
  484. {
  485. // filter keypoints by mask
  486. //KeyPointsFilter::runByPixelsMask( keypoints, mask );
  487. }
  488. if( _descriptors.needed() )
  489. {
  490. //t = (double)getTickCount();
  491. int dsize = descriptorSize();
  492. _descriptors.create((int)keypoints.size(), dsize, descriptor_type);
  493. Mat descriptors = _descriptors.getMat();
  494. calcDescriptors(gpyr, keypoints, descriptors, nOctaveLayers, firstOctave);
  495. //t = (double)getTickCount() - t;
  496. //printf("descriptor extraction time: %g\n", t*1000./tf);
  497. }
  498. }
  499. void SIFT_Impl::read( const FileNode& fn)
  500. {
  501. // if node is empty, keep previous value
  502. if (!fn["nfeatures"].empty())
  503. fn["nfeatures"] >> nfeatures;
  504. if (!fn["nOctaveLayers"].empty())
  505. fn["nOctaveLayers"] >> nOctaveLayers;
  506. if (!fn["contrastThreshold"].empty())
  507. fn["contrastThreshold"] >> contrastThreshold;
  508. if (!fn["edgeThreshold"].empty())
  509. fn["edgeThreshold"] >> edgeThreshold;
  510. if (!fn["sigma"].empty())
  511. fn["sigma"] >> sigma;
  512. if (!fn["descriptorType"].empty())
  513. fn["descriptorType"] >> descriptor_type;
  514. }
  515. void SIFT_Impl::write( FileStorage& fs) const
  516. {
  517. if(fs.isOpened())
  518. {
  519. fs << "name" << getDefaultName();
  520. fs << "nfeatures" << nfeatures;
  521. fs << "nOctaveLayers" << nOctaveLayers;
  522. fs << "contrastThreshold" << contrastThreshold;
  523. fs << "edgeThreshold" << edgeThreshold;
  524. fs << "sigma" << sigma;
  525. fs << "descriptorType" << descriptor_type;
  526. }
  527. }
  528. }