modules/features2d/src/freak.cpp

   1 //  freak.cpp
   2 //
   3 //  Copyright (C) 2011-2012  Signal processing laboratory 2, EPFL,
   4 //  Kirell Benzi (kirell.benzi@epfl.ch),
   5 //  Raphael Ortiz (raphael.ortiz@a3.epfl.ch)
   6 //  Alexandre Alahi (alexandre.alahi@epfl.ch)
   7 //  and Pierre Vandergheynst (pierre.vandergheynst@epfl.ch)
   8 //
   9 //  Redistribution and use in source and binary forms, with or without modification,
  10 //  are permitted provided that the following conditions are met:
  11 //
  12 //   * Redistribution's of source code must retain the above copyright notice,
  13 //     this list of conditions and the following disclaimer.
  14 //
  15 //   * Redistribution's in binary form must reproduce the above copyright notice,
  16 //     this list of conditions and the following disclaimer in the documentation
  17 //     and/or other materials provided with the distribution.
  18 //
  19 //   * The name of the copyright holders may not be used to endorse or promote products
  20 //     derived from this software without specific prior written permission.
  21 //
  22 //  This software is provided by the copyright holders and contributors "as is" and
  23 //  any express or implied warranties, including, but not limited to, the implied
  24 //  warranties of merchantability and fitness for a particular purpose are disclaimed.
  25 //  In no event shall the Intel Corporation or contributors be liable for any direct,
  26 //  indirect, incidental, special, exemplary, or consequential damages
  27 //  (including, but not limited to, procurement of substitute goods or services;
  28 //  loss of use, data, or profits; or business interruption) however caused
  29 //  and on any theory of liability, whether in contract, strict liability,
  30 //  or tort (including negligence or otherwise) arising in any way out of
  31 //  the use of this software, even if advised of the possibility of such damage.
  32
  33 #include "precomp.hpp"
  34 #include <fstream>
  35 #include <stdlib.h>
  36 #include <algorithm>
  37 #include <iostream>
  38 #include <bitset>
  39 #include <sstream>
  40 #include <algorithm>
  41 #include <iomanip>
  42 #include <string.h>
  43
  44 namespace cv
  45 {
  46
  47 static const double FREAK_SQRT2 = 1.4142135623731;
  48 static const double FREAK_INV_SQRT2 = 1.0 / FREAK_SQRT2;
  49 static const double FREAK_LOG2 = 0.693147180559945;
  50 static const int FREAK_NB_ORIENTATION = 256;
  51 static const int FREAK_NB_POINTS = 43;
  52 static const int FREAK_SMALLEST_KP_SIZE = 7; // smallest size of keypoints
  53 static const int FREAK_NB_SCALES = FREAK::NB_SCALES;
  54 static const int FREAK_NB_PAIRS = FREAK::NB_PAIRS;
  55 static const int FREAK_NB_ORIENPAIRS = FREAK::NB_ORIENPAIRS;
  56
  57 // default pairs
  58 static const int FREAK_DEF_PAIRS[FREAK::NB_PAIRS] =
  59 {
  60      404,431,818,511,181,52,311,874,774,543,719,230,417,205,11,
  61      560,149,265,39,306,165,857,250,8,61,15,55,717,44,412,
  62      592,134,761,695,660,782,625,487,549,516,271,665,762,392,178,
  63      796,773,31,672,845,548,794,677,654,241,831,225,238,849,83,
  64      691,484,826,707,122,517,583,731,328,339,571,475,394,472,580,
  65      381,137,93,380,327,619,729,808,218,213,459,141,806,341,95,
  66      382,568,124,750,193,749,706,843,79,199,317,329,768,198,100,
  67      466,613,78,562,783,689,136,838,94,142,164,679,219,419,366,
  68      418,423,77,89,523,259,683,312,555,20,470,684,123,458,453,833,
  69      72,113,253,108,313,25,153,648,411,607,618,128,305,232,301,84,
  70      56,264,371,46,407,360,38,99,176,710,114,578,66,372,653,
  71      129,359,424,159,821,10,323,393,5,340,891,9,790,47,0,175,346,
  72      236,26,172,147,574,561,32,294,429,724,755,398,787,288,299,
  73      769,565,767,722,757,224,465,723,498,467,235,127,802,446,233,
  74      544,482,800,318,16,532,801,441,554,173,60,530,713,469,30,
  75      212,630,899,170,266,799,88,49,512,399,23,500,107,524,90,
  76      194,143,135,192,206,345,148,71,119,101,563,870,158,254,214,
  77      276,464,332,725,188,385,24,476,40,231,620,171,258,67,109,
  78      844,244,187,388,701,690,50,7,850,479,48,522,22,154,12,659,
  79      736,655,577,737,830,811,174,21,237,335,353,234,53,270,62,
  80      182,45,177,245,812,673,355,556,612,166,204,54,248,365,226,
  81      242,452,700,685,573,14,842,481,468,781,564,416,179,405,35,
  82      819,608,624,367,98,643,448,2,460,676,440,240,130,146,184,
  83      185,430,65,807,377,82,121,708,239,310,138,596,730,575,477,
  84      851,797,247,27,85,586,307,779,326,494,856,324,827,96,748,
  85      13,397,125,688,702,92,293,716,277,140,112,4,80,855,839,1,
  86      413,347,584,493,289,696,19,751,379,76,73,115,6,590,183,734,
  87      197,483,217,344,330,400,186,243,587,220,780,200,793,246,824,
  88      41,735,579,81,703,322,760,720,139,480,490,91,814,813,163,
  89      152,488,763,263,425,410,576,120,319,668,150,160,302,491,515,
  90      260,145,428,97,251,395,272,252,18,106,358,854,485,144,550,
  91      131,133,378,68,102,104,58,361,275,209,697,582,338,742,589,
  92      325,408,229,28,304,191,189,110,126,486,211,547,533,70,215,
  93      670,249,36,581,389,605,331,518,442,822
  94 };
  95
  96 // used to sort pairs during pairs selection
  97 struct PairStat
  98 {
  99     double mean;
 100     int idx;
 101 };
 102
 103 struct sortMean
 104 {
 105     bool operator()( const PairStat& a, const PairStat& b ) const
 106     {
 107         return a.mean < b.mean;
 108     }
 109 };
 110
 111 void FREAK::buildPattern()
 112 {
 113     if( patternScale == patternScale0 && nOctaves == nOctaves0 && !patternLookup.empty() )
 114         return;
 115
 116     nOctaves0 = nOctaves;
 117     patternScale0 = patternScale;
 118
 119     patternLookup.resize(FREAK_NB_SCALES*FREAK_NB_ORIENTATION*FREAK_NB_POINTS);
 120     double scaleStep = std::pow(2.0, (double)(nOctaves)/FREAK_NB_SCALES ); // 2 ^ ( (nOctaves-1) /nbScales)
 121     double scalingFactor, alpha, beta, theta = 0;
 122
 123     // pattern definition, radius normalized to 1.0 (outer point position+sigma=1.0)
 124     const int n[8] = {6,6,6,6,6,6,6,1}; // number of points on each concentric circle (from outer to inner)
 125     const double bigR(2.0/3.0); // bigger radius
 126     const double smallR(2.0/24.0); // smaller radius
 127     const double unitSpace( (bigR-smallR)/21.0 ); // define spaces between concentric circles (from center to outer: 1,2,3,4,5,6)
 128     // radii of the concentric cirles (from outer to inner)
 129     const double radius[8] = {bigR, bigR-6*unitSpace, bigR-11*unitSpace, bigR-15*unitSpace, bigR-18*unitSpace, bigR-20*unitSpace, smallR, 0.0};
 130     // sigma of pattern points (each group of 6 points on a concentric cirle has the same sigma)
 131     const double sigma[8] = {radius[0]/2.0, radius[1]/2.0, radius[2]/2.0,
 132                              radius[3]/2.0, radius[4]/2.0, radius[5]/2.0,
 133                              radius[6]/2.0, radius[6]/2.0
 134                             };
 135     // fill the lookup table
 136     for( int scaleIdx=0; scaleIdx < FREAK_NB_SCALES; ++scaleIdx )
 137     {
 138         patternSizes[scaleIdx] = 0; // proper initialization
 139         scalingFactor = std::pow(scaleStep,scaleIdx); //scale of the pattern, scaleStep ^ scaleIdx
 140
 141         for( int orientationIdx = 0; orientationIdx < FREAK_NB_ORIENTATION; ++orientationIdx )
 142         {
 143             theta = double(orientationIdx)* 2*CV_PI/double(FREAK_NB_ORIENTATION); // orientation of the pattern
 144             int pointIdx = 0;
 145
 146             PatternPoint* patternLookupPtr = &patternLookup[0];
 147             for( size_t i = 0; i < 8; ++i )
 148             {
 149                 for( int k = 0 ; k < n[i]; ++k )
 150                 {
 151                     beta = CV_PI/n[i] * (i%2); // orientation offset so that groups of points on each circles are staggered
 152                     alpha = double(k)* 2*CV_PI/double(n[i])+beta+theta;
 153
 154                     // add the point to the look-up table
 155                     PatternPoint& point = patternLookupPtr[ scaleIdx*FREAK_NB_ORIENTATION*FREAK_NB_POINTS+orientationIdx*FREAK_NB_POINTS+pointIdx ];
 156                     point.x = static_cast<float>(radius[i] * cos(alpha) * scalingFactor * patternScale);
 157                     point.y = static_cast<float>(radius[i] * sin(alpha) * scalingFactor * patternScale);
 158                     point.sigma = static_cast<float>(sigma[i] * scalingFactor * patternScale);
 159
 160                     // adapt the sizeList if necessary
 161                     const int sizeMax = static_cast<int>(ceil((radius[i]+sigma[i])*scalingFactor*patternScale)) + 1;
 162                     if( patternSizes[scaleIdx] < sizeMax )
 163                         patternSizes[scaleIdx] = sizeMax;
 164
 165                     ++pointIdx;
 166                 }
 167             }
 168         }
 169     }
 170
 171     // build the list of orientation pairs
 172     orientationPairs[0].i=0; orientationPairs[0].j=3; orientationPairs[1].i=1; orientationPairs[1].j=4; orientationPairs[2].i=2; orientationPairs[2].j=5;
 173     orientationPairs[3].i=0; orientationPairs[3].j=2; orientationPairs[4].i=1; orientationPairs[4].j=3; orientationPairs[5].i=2; orientationPairs[5].j=4;
 174     orientationPairs[6].i=3; orientationPairs[6].j=5; orientationPairs[7].i=4; orientationPairs[7].j=0; orientationPairs[8].i=5; orientationPairs[8].j=1;
 175
 176     orientationPairs[9].i=6; orientationPairs[9].j=9; orientationPairs[10].i=7; orientationPairs[10].j=10; orientationPairs[11].i=8; orientationPairs[11].j=11;
 177     orientationPairs[12].i=6; orientationPairs[12].j=8; orientationPairs[13].i=7; orientationPairs[13].j=9; orientationPairs[14].i=8; orientationPairs[14].j=10;
 178     orientationPairs[15].i=9; orientationPairs[15].j=11; orientationPairs[16].i=10; orientationPairs[16].j=6; orientationPairs[17].i=11; orientationPairs[17].j=7;
 179
 180     orientationPairs[18].i=12; orientationPairs[18].j=15; orientationPairs[19].i=13; orientationPairs[19].j=16; orientationPairs[20].i=14; orientationPairs[20].j=17;
 181     orientationPairs[21].i=12; orientationPairs[21].j=14; orientationPairs[22].i=13; orientationPairs[22].j=15; orientationPairs[23].i=14; orientationPairs[23].j=16;
 182     orientationPairs[24].i=15; orientationPairs[24].j=17; orientationPairs[25].i=16; orientationPairs[25].j=12; orientationPairs[26].i=17; orientationPairs[26].j=13;
 183
 184     orientationPairs[27].i=18; orientationPairs[27].j=21; orientationPairs[28].i=19; orientationPairs[28].j=22; orientationPairs[29].i=20; orientationPairs[29].j=23;
 185     orientationPairs[30].i=18; orientationPairs[30].j=20; orientationPairs[31].i=19; orientationPairs[31].j=21; orientationPairs[32].i=20; orientationPairs[32].j=22;
 186     orientationPairs[33].i=21; orientationPairs[33].j=23; orientationPairs[34].i=22; orientationPairs[34].j=18; orientationPairs[35].i=23; orientationPairs[35].j=19;
 187
 188     orientationPairs[36].i=24; orientationPairs[36].j=27; orientationPairs[37].i=25; orientationPairs[37].j=28; orientationPairs[38].i=26; orientationPairs[38].j=29;
 189     orientationPairs[39].i=30; orientationPairs[39].j=33; orientationPairs[40].i=31; orientationPairs[40].j=34; orientationPairs[41].i=32; orientationPairs[41].j=35;
 190     orientationPairs[42].i=36; orientationPairs[42].j=39; orientationPairs[43].i=37; orientationPairs[43].j=40; orientationPairs[44].i=38; orientationPairs[44].j=41;
 191
 192     for( unsigned m = FREAK_NB_ORIENPAIRS; m--; )
 193     {
 194         const float dx = patternLookup[orientationPairs[m].i].x-patternLookup[orientationPairs[m].j].x;
 195         const float dy = patternLookup[orientationPairs[m].i].y-patternLookup[orientationPairs[m].j].y;
 196         const float norm_sq = (dx*dx+dy*dy);
 197         orientationPairs[m].weight_dx = int((dx/(norm_sq))*4096.0+0.5);
 198         orientationPairs[m].weight_dy = int((dy/(norm_sq))*4096.0+0.5);
 199     }
 200
 201     // build the list of description pairs
 202     std::vector<DescriptionPair> allPairs;
 203     for( unsigned int i = 1; i < (unsigned int)FREAK_NB_POINTS; ++i )
 204     {
 205         // (generate all the pairs)
 206         for( unsigned int j = 0; (unsigned int)j < i; ++j )
 207         {
 208             DescriptionPair pair = {(uchar)i,(uchar)j};
 209             allPairs.push_back(pair);
 210         }
 211     }
 212     // Input vector provided
 213     if( !selectedPairs0.empty() )
 214     {
 215         if( (int)selectedPairs0.size() == FREAK_NB_PAIRS )
 216         {
 217             for( int i = 0; i < FREAK_NB_PAIRS; ++i )
 218                  descriptionPairs[i] = allPairs[selectedPairs0.at(i)];
 219         }
 220         else
 221         {
 222             CV_Error(Error::StsVecLengthErr, "Input vector does not match the required size");
 223         }
 224     }
 225     else // default selected pairs
 226     {
 227         for( int i = 0; i < FREAK_NB_PAIRS; ++i )
 228              descriptionPairs[i] = allPairs[FREAK_DEF_PAIRS[i]];
 229     }
 230 }
 231
 232 void FREAK::computeImpl( InputArray _image, std::vector<KeyPoint>& keypoints, OutputArray _descriptors ) const
 233 {
 234     Mat image = _image.getMat();
 235     if( image.empty() )
 236         return;
 237     if( keypoints.empty() )
 238         return;
 239
 240     ((FREAK*)this)->buildPattern();
 241
 242     Mat imgIntegral;
 243     integral(image, imgIntegral);
 244     std::vector<int> kpScaleIdx(keypoints.size()); // used to save pattern scale index corresponding to each keypoints
 245     const std::vector<int>::iterator ScaleIdxBegin = kpScaleIdx.begin(); // used in std::vector erase function
 246     const std::vector<cv::KeyPoint>::iterator kpBegin = keypoints.begin(); // used in std::vector erase function
 247     const float sizeCst = static_cast<float>(FREAK_NB_SCALES/(FREAK_LOG2* nOctaves));
 248     uchar pointsValue[FREAK_NB_POINTS];
 249     int thetaIdx = 0;
 250     int direction0;
 251     int direction1;
 252
 253     // compute the scale index corresponding to the keypoint size and remove keypoints close to the border
 254     if( scaleNormalized )
 255     {
 256         for( size_t k = keypoints.size(); k--; )
 257         {
 258             //Is k non-zero? If so, decrement it and continue"
 259             kpScaleIdx[k] = std::max( (int)(std::log(keypoints[k].size/FREAK_SMALLEST_KP_SIZE)*sizeCst+0.5) ,0);
 260             if( kpScaleIdx[k] >= FREAK_NB_SCALES )
 261                 kpScaleIdx[k] = FREAK_NB_SCALES-1;
 262
 263             if( keypoints[k].pt.x <= patternSizes[kpScaleIdx[k]] || //check if the description at this specific position and scale fits inside the image
 264                  keypoints[k].pt.y <= patternSizes[kpScaleIdx[k]] ||
 265                  keypoints[k].pt.x >= image.cols-patternSizes[kpScaleIdx[k]] ||
 266                  keypoints[k].pt.y >= image.rows-patternSizes[kpScaleIdx[k]]
 267                )
 268             {
 269                 keypoints.erase(kpBegin+k);
 270                 kpScaleIdx.erase(ScaleIdxBegin+k);
 271             }
 272         }
 273     }
 274     else
 275     {
 276         const int scIdx = std::max( (int)(1.0986122886681*sizeCst+0.5) ,0);
 277         for( size_t k = keypoints.size(); k--; )
 278         {
 279             kpScaleIdx[k] = scIdx; // equivalent to the formule when the scale is normalized with a constant size of keypoints[k].size=3*SMALLEST_KP_SIZE
 280             if( kpScaleIdx[k] >= FREAK_NB_SCALES )
 281             {
 282                 kpScaleIdx[k] = FREAK_NB_SCALES-1;
 283             }
 284             if( keypoints[k].pt.x <= patternSizes[kpScaleIdx[k]] ||
 285                 keypoints[k].pt.y <= patternSizes[kpScaleIdx[k]] ||
 286                 keypoints[k].pt.x >= image.cols-patternSizes[kpScaleIdx[k]] ||
 287                 keypoints[k].pt.y >= image.rows-patternSizes[kpScaleIdx[k]]
 288                )
 289             {
 290                 keypoints.erase(kpBegin+k);
 291                 kpScaleIdx.erase(ScaleIdxBegin+k);
 292             }
 293         }
 294     }
 295
 296     // allocate descriptor memory, estimate orientations, extract descriptors
 297     if( !extAll )
 298     {
 299         // extract the best comparisons only
 300         _descriptors.create((int)keypoints.size(), FREAK_NB_PAIRS/8, CV_8U);
 301         _descriptors.setTo(Scalar::all(0));
 302         Mat descriptors = _descriptors.getMat();
 303 #if CV_SSE2
 304         __m128i* ptr= (__m128i*) (descriptors.data+(keypoints.size()-1)*descriptors.step[0]);
 305 #else
 306         std::bitset<FREAK_NB_PAIRS>* ptr = (std::bitset<FREAK_NB_PAIRS>*) (descriptors.data+(keypoints.size()-1)*descriptors.step[0]);
 307 #endif
 308         for( size_t k = keypoints.size(); k--; )
 309         {
 310             // estimate orientation (gradient)
 311             if( !orientationNormalized )
 312             {
 313                 thetaIdx = 0; // assign 0° to all keypoints
 314                 keypoints[k].angle = 0.0;
 315             }
 316             else
 317             {
 318                 // get the points intensity value in the un-rotated pattern
 319                 for( int i = FREAK_NB_POINTS; i--; )
 320                 {
 321                     pointsValue[i] = meanIntensity(image, imgIntegral, keypoints[k].pt.x,keypoints[k].pt.y, kpScaleIdx[k], 0, i);
 322                 }
 323                 direction0 = 0;
 324                 direction1 = 0;
 325                 for( int m = 45; m--; )
 326                 {
 327                     //iterate through the orientation pairs
 328                     const int delta = (pointsValue[ orientationPairs[m].i ]-pointsValue[ orientationPairs[m].j ]);
 329                     direction0 += delta*(orientationPairs[m].weight_dx)/2048;
 330                     direction1 += delta*(orientationPairs[m].weight_dy)/2048;
 331                 }
 332
 333                 keypoints[k].angle = static_cast<float>(atan2((float)direction1,(float)direction0)*(180.0/CV_PI));//estimate orientation
 334                 thetaIdx = int(FREAK_NB_ORIENTATION*keypoints[k].angle*(1/360.0)+0.5);
 335                 if( thetaIdx < 0 )
 336                     thetaIdx += FREAK_NB_ORIENTATION;
 337
 338                 if( thetaIdx >= FREAK_NB_ORIENTATION )
 339                     thetaIdx -= FREAK_NB_ORIENTATION;
 340             }
 341             // extract descriptor at the computed orientation
 342             for( int i = FREAK_NB_POINTS; i--; )
 343             {
 344                 pointsValue[i] = meanIntensity(image, imgIntegral, keypoints[k].pt.x,keypoints[k].pt.y, kpScaleIdx[k], thetaIdx, i);
 345             }
 346 #if CV_SSE2
 347             // note that comparisons order is modified in each block (but first 128 comparisons remain globally the same-->does not affect the 128,384 bits segmanted matching strategy)
 348             int cnt = 0;
 349             for( int n = FREAK_NB_PAIRS/128; n-- ; )
 350             {
 351                 __m128i result128 = _mm_setzero_si128();
 352                 for( int m = 128/16; m--; cnt += 16 )
 353                 {
 354                     __m128i operand1 = _mm_set_epi8(
 355                         pointsValue[descriptionPairs[cnt+0].i],
 356                         pointsValue[descriptionPairs[cnt+1].i],
 357                         pointsValue[descriptionPairs[cnt+2].i],
 358                         pointsValue[descriptionPairs[cnt+3].i],
 359                         pointsValue[descriptionPairs[cnt+4].i],
 360                         pointsValue[descriptionPairs[cnt+5].i],
 361                         pointsValue[descriptionPairs[cnt+6].i],
 362                         pointsValue[descriptionPairs[cnt+7].i],
 363                         pointsValue[descriptionPairs[cnt+8].i],
 364                         pointsValue[descriptionPairs[cnt+9].i],
 365                         pointsValue[descriptionPairs[cnt+10].i],
 366                         pointsValue[descriptionPairs[cnt+11].i],
 367                         pointsValue[descriptionPairs[cnt+12].i],
 368                         pointsValue[descriptionPairs[cnt+13].i],
 369                         pointsValue[descriptionPairs[cnt+14].i],
 370                         pointsValue[descriptionPairs[cnt+15].i]);
 371
 372                     __m128i operand2 = _mm_set_epi8(
 373                         pointsValue[descriptionPairs[cnt+0].j],
 374                         pointsValue[descriptionPairs[cnt+1].j],
 375                         pointsValue[descriptionPairs[cnt+2].j],
 376                         pointsValue[descriptionPairs[cnt+3].j],
 377                         pointsValue[descriptionPairs[cnt+4].j],
 378                         pointsValue[descriptionPairs[cnt+5].j],
 379                         pointsValue[descriptionPairs[cnt+6].j],
 380                         pointsValue[descriptionPairs[cnt+7].j],
 381                         pointsValue[descriptionPairs[cnt+8].j],
 382                         pointsValue[descriptionPairs[cnt+9].j],
 383                         pointsValue[descriptionPairs[cnt+10].j],
 384                         pointsValue[descriptionPairs[cnt+11].j],
 385                         pointsValue[descriptionPairs[cnt+12].j],
 386                         pointsValue[descriptionPairs[cnt+13].j],
 387                         pointsValue[descriptionPairs[cnt+14].j],
 388                         pointsValue[descriptionPairs[cnt+15].j]);
 389
 390                     __m128i workReg = _mm_min_epu8(operand1, operand2); // emulated "not less than" for 8-bit UNSIGNED integers
 391                     workReg = _mm_cmpeq_epi8(workReg, operand2);        // emulated "not less than" for 8-bit UNSIGNED integers
 392
 393                     workReg = _mm_and_si128(_mm_set1_epi16(short(0x8080 >> m)), workReg); // merge the last 16 bits with the 128bits std::vector until full
 394                     result128 = _mm_or_si128(result128, workReg);
 395                 }
 396                 (*ptr) = result128;
 397                 ++ptr;
 398             }
 399             ptr -= 8;
 400 #else
 401             // extracting descriptor preserving the order of SSE version
 402             int cnt = 0;
 403             for( int n = 7; n < FREAK_NB_PAIRS; n += 128)
 404             {
 405                 for( int m = 8; m--; )
 406                 {
 407                     int nm = n-m;
 408                     for(int kk = nm+15*8; kk >= nm; kk-=8, ++cnt)
 409                     {
 410                         ptr->set(kk, pointsValue[descriptionPairs[cnt].i] >= pointsValue[descriptionPairs[cnt].j]);
 411                     }
 412                 }
 413             }
 414             --ptr;
 415 #endif
 416         }
 417     }
 418     else // extract all possible comparisons for selection
 419     {
 420         _descriptors.create((int)keypoints.size(), 128, CV_8U);
 421         _descriptors.setTo(Scalar::all(0));
 422         Mat descriptors = _descriptors.getMat();
 423         std::bitset<1024>* ptr = (std::bitset<1024>*) (descriptors.data+(keypoints.size()-1)*descriptors.step[0]);
 424
 425         for( size_t k = keypoints.size(); k--; )
 426         {
 427             //estimate orientation (gradient)
 428             if( !orientationNormalized )
 429             {
 430                 thetaIdx = 0;//assign 0° to all keypoints
 431                 keypoints[k].angle = 0.0;
 432             }
 433             else
 434             {
 435                 //get the points intensity value in the un-rotated pattern
 436                 for( int i = FREAK_NB_POINTS;i--; )
 437                     pointsValue[i] = meanIntensity(image, imgIntegral, keypoints[k].pt.x,keypoints[k].pt.y, kpScaleIdx[k], 0, i);
 438
 439                 direction0 = 0;
 440                 direction1 = 0;
 441                 for( int m = 45; m--; )
 442                 {
 443                     //iterate through the orientation pairs
 444                     const int delta = (pointsValue[ orientationPairs[m].i ]-pointsValue[ orientationPairs[m].j ]);
 445                     direction0 += delta*(orientationPairs[m].weight_dx)/2048;
 446                     direction1 += delta*(orientationPairs[m].weight_dy)/2048;
 447                 }
 448
 449                 keypoints[k].angle = static_cast<float>(atan2((float)direction1,(float)direction0)*(180.0/CV_PI)); //estimate orientation
 450                 thetaIdx = int(FREAK_NB_ORIENTATION*keypoints[k].angle*(1/360.0)+0.5);
 451
 452                 if( thetaIdx < 0 )
 453                     thetaIdx += FREAK_NB_ORIENTATION;
 454
 455                 if( thetaIdx >= FREAK_NB_ORIENTATION )
 456                     thetaIdx -= FREAK_NB_ORIENTATION;
 457             }
 458             // get the points intensity value in the rotated pattern
 459             for( int i = FREAK_NB_POINTS; i--; )
 460             {
 461                 pointsValue[i] = meanIntensity(image, imgIntegral, keypoints[k].pt.x,
 462                                              keypoints[k].pt.y, kpScaleIdx[k], thetaIdx, i);
 463             }
 464
 465             int cnt(0);
 466             for( int i = 1; i < FREAK_NB_POINTS; ++i )
 467             {
 468                 //(generate all the pairs)
 469                 for( int j = 0; j < i; ++j )
 470                 {
 471                     ptr->set(cnt, pointsValue[i] >= pointsValue[j] );
 472                     ++cnt;
 473                 }
 474             }
 475             --ptr;
 476         }
 477     }
 478 }
 479
 480 // simply take average on a square patch, not even gaussian approx
 481 uchar FREAK::meanIntensity( InputArray _image, InputArray _integral,
 482                             const float kp_x,
 483                             const float kp_y,
 484                             const unsigned int scale,
 485                             const unsigned int rot,
 486                             const unsigned int point) const
 487 {
 488     Mat image = _image.getMat(), integral = _integral.getMat();
 489     // get point position in image
 490     const PatternPoint& FreakPoint = patternLookup[scale*FREAK_NB_ORIENTATION*FREAK_NB_POINTS + rot*FREAK_NB_POINTS + point];
 491     const float xf = FreakPoint.x+kp_x;
 492     const float yf = FreakPoint.y+kp_y;
 493     const int x = int(xf);
 494     const int y = int(yf);
 495     const int& imagecols = image.cols;
 496
 497     // get the sigma:
 498     const float radius = FreakPoint.sigma;
 499
 500     // calculate output:
 501     if( radius < 0.5 )
 502     {
 503         // interpolation multipliers:
 504         const int r_x = static_cast<int>((xf-x)*1024);
 505         const int r_y = static_cast<int>((yf-y)*1024);
 506         const int r_x_1 = (1024-r_x);
 507         const int r_y_1 = (1024-r_y);
 508         uchar* ptr = image.data+x+y*imagecols;
 509         unsigned int ret_val;
 510         // linear interpolation:
 511         ret_val = (r_x_1*r_y_1*int(*ptr));
 512         ptr++;
 513         ret_val += (r_x*r_y_1*int(*ptr));
 514         ptr += imagecols;
 515         ret_val += (r_x*r_y*int(*ptr));
 516         ptr--;
 517         ret_val += (r_x_1*r_y*int(*ptr));
 518         //return the rounded mean
 519         ret_val += 2 * 1024 * 1024;
 520         return static_cast<uchar>(ret_val / (4 * 1024 * 1024));
 521     }
 522
 523     // expected case:
 524
 525     // calculate borders
 526     const int x_left = int(xf-radius+0.5);
 527     const int y_top = int(yf-radius+0.5);
 528     const int x_right = int(xf+radius+1.5);//integral image is 1px wider
 529     const int y_bottom = int(yf+radius+1.5);//integral image is 1px higher
 530     int ret_val;
 531
 532     ret_val = integral.at<int>(y_bottom,x_right);//bottom right corner
 533     ret_val -= integral.at<int>(y_bottom,x_left);
 534     ret_val += integral.at<int>(y_top,x_left);
 535     ret_val -= integral.at<int>(y_top,x_right);
 536     ret_val = ret_val/( (x_right-x_left)* (y_bottom-y_top) );
 537     //~ std::cout<<integral.step[1]<<std::endl;
 538     return static_cast<uchar>(ret_val);
 539 }
 540
 541 // pair selection algorithm from a set of training images and corresponding keypoints
 542 std::vector<int> FREAK::selectPairs(const std::vector<Mat>& images
 543                                         , std::vector<std::vector<KeyPoint> >& keypoints
 544                                         , const double corrTresh
 545                                         , bool verbose )
 546 {
 547     extAll = true;
 548     // compute descriptors with all pairs
 549     Mat descriptors;
 550
 551     if( verbose )
 552         std::cout << "Number of images: " << images.size() << std::endl;
 553
 554     for( size_t i = 0;i < images.size(); ++i )
 555     {
 556         Mat descriptorsTmp;
 557         computeImpl(images[i],keypoints[i],descriptorsTmp);
 558         descriptors.push_back(descriptorsTmp);
 559     }
 560
 561     if( verbose )
 562         std::cout << "number of keypoints: " << descriptors.rows << std::endl;
 563
 564     //descriptor in floating point format (each bit is a float)
 565     Mat descriptorsFloat = Mat::zeros(descriptors.rows, 903, CV_32F);
 566
 567     std::bitset<1024>* ptr = (std::bitset<1024>*) (descriptors.data+(descriptors.rows-1)*descriptors.step[0]);
 568     for( int m = descriptors.rows; m--; )
 569     {
 570         for( int n = 903; n--; )
 571         {
 572             if( ptr->test(n) == true )
 573                 descriptorsFloat.at<float>(m,n)=1.0f;
 574         }
 575         --ptr;
 576     }
 577
 578     std::vector<PairStat> pairStat;
 579     for( int n = 903; n--; )
 580     {
 581         // the higher the variance, the better --> mean = 0.5
 582         PairStat tmp = { fabs( mean(descriptorsFloat.col(n))[0]-0.5 ) ,n};
 583         pairStat.push_back(tmp);
 584     }
 585
 586     std::sort( pairStat.begin(),pairStat.end(), sortMean() );
 587
 588     std::vector<PairStat> bestPairs;
 589     for( int m = 0; m < 903; ++m )
 590     {
 591         if( verbose )
 592             std::cout << m << ":" << bestPairs.size() << " " << std::flush;
 593         double corrMax(0);
 594
 595         for( size_t n = 0; n < bestPairs.size(); ++n )
 596         {
 597             int idxA = bestPairs[n].idx;
 598             int idxB = pairStat[m].idx;
 599             double corr(0);
 600             // compute correlation between 2 pairs
 601             corr = fabs(compareHist(descriptorsFloat.col(idxA), descriptorsFloat.col(idxB), HISTCMP_CORREL));
 602
 603             if( corr > corrMax )
 604             {
 605                 corrMax = corr;
 606                 if( corrMax >= corrTresh )
 607                     break;
 608             }
 609         }
 610
 611         if( corrMax < corrTresh/*0.7*/ )
 612             bestPairs.push_back(pairStat[m]);
 613
 614         if( bestPairs.size() >= 512 )
 615         {
 616             if( verbose )
 617                 std::cout << m << std::endl;
 618             break;
 619         }
 620     }
 621
 622     std::vector<int> idxBestPairs;
 623     if( (int)bestPairs.size() >= FREAK_NB_PAIRS )
 624     {
 625         for( int i = 0; i < FREAK_NB_PAIRS; ++i )
 626             idxBestPairs.push_back(bestPairs[i].idx);
 627     }
 628     else
 629     {
 630         if( verbose )
 631             std::cout << "correlation threshold too small (restrictive)" << std::endl;
 632         CV_Error(Error::StsError, "correlation threshold too small (restrictive)");
 633     }
 634     extAll = false;
 635     return idxBestPairs;
 636 }
 637
 638
 639 /*
 640 // create an image showing the brisk pattern
 641 void FREAKImpl::drawPattern()
 642 {
 643     Mat pattern = Mat::zeros(1000, 1000, CV_8UC3) + Scalar(255,255,255);
 644     int sFac = 500 / patternScale;
 645     for( int n = 0; n < kNB_POINTS; ++n )
 646     {
 647         PatternPoint& pt = patternLookup[n];
 648         circle(pattern, Point( pt.x*sFac,pt.y*sFac)+Point(500,500), pt.sigma*sFac, Scalar(0,0,255),2);
 649         // rectangle(pattern, Point( (pt.x-pt.sigma)*sFac,(pt.y-pt.sigma)*sFac)+Point(500,500), Point( (pt.x+pt.sigma)*sFac,(pt.y+pt.sigma)*sFac)+Point(500,500), Scalar(0,0,255),2);
 650
 651         circle(pattern, Point( pt.x*sFac,pt.y*sFac)+Point(500,500), 1, Scalar(0,0,0),3);
 652         std::ostringstream oss;
 653         oss << n;
 654         putText( pattern, oss.str(), Point( pt.x*sFac,pt.y*sFac)+Point(500,500), FONT_HERSHEY_SIMPLEX,0.5, Scalar(0,0,0), 1);
 655     }
 656     imshow( "FreakDescriptorExtractor pattern", pattern );
 657     waitKey(0);
 658 }
 659 */
 660
 661 // -------------------------------------------------
 662 /* FREAK interface implementation */
 663 FREAK::FREAK( bool _orientationNormalized, bool _scaleNormalized
 664             , float _patternScale, int _nOctaves, const std::vector<int>& _selectedPairs )
 665     : orientationNormalized(_orientationNormalized), scaleNormalized(_scaleNormalized),
 666     patternScale(_patternScale), nOctaves(_nOctaves), extAll(false), nOctaves0(0), selectedPairs0(_selectedPairs)
 667 {
 668 }
 669
 670 FREAK::~FREAK()
 671 {
 672 }
 673
 674 int FREAK::descriptorSize() const
 675 {
 676     return FREAK_NB_PAIRS / 8; // descriptor length in bytes
 677 }
 678
 679 int FREAK::descriptorType() const
 680 {
 681     return CV_8U;
 682 }
 683
 684 int FREAK::defaultNorm() const
 685 {
 686     return NORM_HAMMING;
 687 }
 688
 689 } // END NAMESPACE CV