Improve initialization performance of Brisk

reformatting

Improve initialization performance of Brisk

fix formatting

Improve initialization performance of Brisk

formatting

Improve initialization performance of Brisk

make a lookup table for ring

use cosine/sine lookup table for theta in brisk and utilize trig identity

fix ring lookup table

use cosine/sine lookup table for theta in brisk and utilize trig identity

formatting

use cosine/sine lookup table for theta in brisk and utilize trig identity

move scale radius product to ring loop to ensure it's not recomputed for each rot

revert change

move scale radius product to ring loop to ensure it's not recomputed for each rot

remove rings lookup table

move scale radius product to ring loop to ensure it's not recomputed for each rot

fix formatting of for loop

move scale radius product to ring loop to ensure it's not recomputed for each rot

use sine/cosine approximations for brisk lookup table.

add documentation for sine/cosine lookup tables

Improve initialization performance of BRISK

Author: danielenricocahall

modules/features2d/src/brisk.cpp

@@ -353,13 +353,30 @@ BRISK_Impl::generateKernel(const std::vector<float> &radiusList,
   const int rings = (int)radiusList.size();
   CV_Assert(radiusList.size() != 0 && radiusList.size() == numberList.size());
   points_ = 0; // remember the total number of points
+  double sineThetaLookupTable[n_rot_];
+  double cosThetaLookupTable[n_rot_];
   for (int ring = 0; ring < rings; ring++)
   {
     points_ += numberList[ring];
   }
+
+  // using a sine/cosine approximation for the lookup table
+  // utilizes the trig identities:
+  // sin(a + b) = sin(a)cos(b) + cos(a)sin(b)
+  // cos(a + b) = cos(a)cos(b) - sin(a)sin(b)
+  // and the fact that sin(0) = 0, cos(0) = 1
+  double cosval = 1., sinval = 0.;
+  double dcos = cos(2*CV_PI/double(n_rot_)), dsin = sin(2*CV_PI/double(n_rot_));
+  for( size_t rot = 0; rot < n_rot_; ++rot)
+  {
+    sineThetaLookupTable[rot] = sinval;
+    cosThetaLookupTable[rot] = cosval;
+    double t = sinval*dcos + cosval*dsin;
+    cosval = cosval*dcos - sinval*dsin;
+    sinval = t;
+  }
   // set up the patterns
   patternPoints_ = new BriskPatternPoint[points_ * scales_ * n_rot_];
-  BriskPatternPoint* patternIterator = patternPoints_;
 
   // define the scale discretization:
   static const float lb_scale = (float)(std::log(scalerange_) / std::log(2.0));
@@ -370,46 +387,51 @@ BRISK_Impl::generateKernel(const std::vector<float> &radiusList,
 
   const float sigma_scale = 1.3f;
 
-  for (unsigned int scale = 0; scale < scales_; ++scale)
-  {
-    scaleList_[scale] = (float)std::pow((double) 2.0, (double) (scale * lb_scale_step));
-    sizeList_[scale] = 0;
-
-    // generate the pattern points look-up
-    double alpha, theta;
-    for (size_t rot = 0; rot < n_rot_; ++rot)
-    {
-      theta = double(rot) * 2 * CV_PI / double(n_rot_); // this is the rotation of the feature
-      for (int ring = 0; ring < rings; ++ring)
-      {
-        for (int num = 0; num < numberList[ring]; ++num)
-        {
-          // the actual coordinates on the circle
-          alpha = (double(num)) * 2 * CV_PI / double(numberList[ring]);
-          patternIterator->x = (float)(scaleList_[scale] * radiusList[ring] * cos(alpha + theta)); // feature rotation plus angle of the point
-          patternIterator->y = (float)(scaleList_[scale] * radiusList[ring] * sin(alpha + theta));
-          // and the gaussian kernel sigma
-          if (ring == 0)
-          {
-            patternIterator->sigma = sigma_scale * scaleList_[scale] * 0.5f;
-          }
-          else
-          {
-            patternIterator->sigma = (float)(sigma_scale * scaleList_[scale] * (double(radiusList[ring]))
-                                     * sin(CV_PI / numberList[ring]));
+  for (unsigned int scale = 0; scale < scales_; ++scale) {
+      scaleList_[scale] = (float) std::pow((double) 2.0, (double) (scale * lb_scale_step));
+      sizeList_[scale] = 0;
+      BriskPatternPoint *patternIteratorOuter = patternPoints_ + (scale * n_rot_ * points_);
+      // generate the pattern points look-up
+      for (int ring = 0; ring < rings; ++ring) {
+          double scaleRadiusProduct = scaleList_[scale] * radiusList[ring];
+          float patternSigma = 0.0f;
+          if (ring == 0) {
+              patternSigma = sigma_scale * scaleList_[scale] * 0.5f;
+          } else {
+              patternSigma = (float) (sigma_scale * scaleList_[scale] * (double(radiusList[ring]))
+                                      * sin(CV_PI / numberList[ring]));
           }
           // adapt the sizeList if necessary
-          const unsigned int size = cvCeil(((scaleList_[scale] * radiusList[ring]) + patternIterator->sigma)) + 1;
-          if (sizeList_[scale] < size)
-          {
-            sizeList_[scale] = size;
+          const unsigned int size = cvCeil(((scaleList_[scale] * radiusList[ring]) + patternSigma)) + 1;
+          if (sizeList_[scale] < size) {
+              sizeList_[scale] = size;
+          }
+          for (int num = 0; num < numberList[ring]; ++num) {
+              BriskPatternPoint *patternIterator = patternIteratorOuter;
+              double alpha = (double(num)) * 2 * CV_PI / double(numberList[ring]);
+              double sine_alpha = sin(alpha);
+              double cosine_alpha = cos(alpha);
+
+              for (size_t rot = 0; rot < n_rot_; ++rot) {
+                  double cosine_theta = cosThetaLookupTable[rot];
+                  double sine_theta = sineThetaLookupTable[rot];
+
+                  // the actual coordinates on the circle
+                  // sin(a + b) = sin(a) cos(b) + cos(a) sin(b)
+                  // cos(a + b) = cos(a) cos(b) - sin(a) sin(b)
+                  patternIterator->x = (float) (scaleRadiusProduct *
+                                                (cosine_theta * cosine_alpha -
+                                                 sine_theta * sine_alpha)); // feature rotation plus angle of the point
+                  patternIterator->y = (float) (scaleRadiusProduct *
+                                                (sine_theta * cosine_alpha + cosine_theta * sine_alpha));
+                  patternIterator->sigma = patternSigma;
+                  // and the gaussian kernel sigma
+                  // increment the iterator
+                  patternIterator += points_;
+              }
+              ++patternIteratorOuter;
           }
-
-          // increment the iterator
-          ++patternIterator;
-        }
       }
-    }
   }
 
   // now also generate pairings