Updating STAR detector and FREAK descriptor to work with large and/or 16-bit images

modules/features2d/src/freak.cpp

@@ -224,13 +224,128 @@ void FREAK::computeImpl( const Mat& image, std::vector<KeyPoint>& keypoints, Mat
 
     ((FREAK*)this)->buildPattern();
 
+    // Convert to gray if not already
+    Mat grayImage = image;
+    if( image.channels() > 1 )
+        cvtColor( image, grayImage, COLOR_BGR2GRAY );
+
+    // Use 32-bit integers if we won't overflow in the integral image
+    if ((image.depth() == CV_8U || image.depth() == CV_8S) &&
+        (image.rows * image.cols) < 8388608 ) // 8388608 = 2 ^ (32 - 8(bit depth) - 1(sign bit))
+    {
+        // Create the integral image appropriate for our type & usage
+        if (image.depth() == CV_8U)
+            computeDescriptors<uchar, int>(grayImage, keypoints, descriptors);
+        else if (image.depth() == CV_8S)
+            computeDescriptors<char, int>(grayImage, keypoints, descriptors);
+        else
+            CV_Error( Error::StsUnsupportedFormat, "" );
+    } else {
+        // Create the integral image appropriate for our type & usage
+        if ( image.depth() == CV_8U )
+            computeDescriptors<uchar, double>(grayImage, keypoints, descriptors);
+        else if ( image.depth() == CV_8S )
+            computeDescriptors<char, double>(grayImage, keypoints, descriptors);
+        else if ( image.depth() == CV_16U )
+            computeDescriptors<ushort, double>(grayImage, keypoints, descriptors);
+        else if ( image.depth() == CV_16S )
+            computeDescriptors<short, double>(grayImage, keypoints, descriptors);
+        else
+            CV_Error( Error::StsUnsupportedFormat, "" );
+    }
+}
+
+template <typename srcMatType>
+void FREAK::extractDescriptor(srcMatType *pointsValue, void ** ptr) const
+{
+    std::bitset<FREAK_NB_PAIRS>** ptrScalar = (std::bitset<FREAK_NB_PAIRS>**) ptr;
+
+    // extracting descriptor preserving the order of SSE version
+    int cnt = 0;
+    for( int n = 7; n < FREAK_NB_PAIRS; n += 128)
+    {
+        for( int m = 8; m--; )
+        {
+            int nm = n-m;
+            for(int kk = nm+15*8; kk >= nm; kk-=8, ++cnt)
+            {
+                (*ptrScalar)->set(kk, pointsValue[descriptionPairs[cnt].i] >= pointsValue[descriptionPairs[cnt].j]);
+            }
+        }
+    }
+    --(*ptrScalar);
+}
+
+#if CV_SSE2
+template <>
+void FREAK::extractDescriptor(uchar *pointsValue, void ** ptr) const
+{
+    __m128i** ptrSSE = (__m128i**) ptr;
+
+    // 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)
+    int cnt = 0;
+    for( int n = FREAK_NB_PAIRS/128; n-- ; )
+    {
+        __m128i result128 = _mm_setzero_si128();
+        for( int m = 128/16; m--; cnt += 16 )
+        {
+            __m128i operand1 = _mm_set_epi8(pointsValue[descriptionPairs[cnt+0].i],
+                                            pointsValue[descriptionPairs[cnt+1].i],
+                                            pointsValue[descriptionPairs[cnt+2].i],
+                                            pointsValue[descriptionPairs[cnt+3].i],
+                                            pointsValue[descriptionPairs[cnt+4].i],
+                                            pointsValue[descriptionPairs[cnt+5].i],
+                                            pointsValue[descriptionPairs[cnt+6].i],
+                                            pointsValue[descriptionPairs[cnt+7].i],
+                                            pointsValue[descriptionPairs[cnt+8].i],
+                                            pointsValue[descriptionPairs[cnt+9].i],
+                                            pointsValue[descriptionPairs[cnt+10].i],
+                                            pointsValue[descriptionPairs[cnt+11].i],
+                                            pointsValue[descriptionPairs[cnt+12].i],
+                                            pointsValue[descriptionPairs[cnt+13].i],
+                                            pointsValue[descriptionPairs[cnt+14].i],
+                                            pointsValue[descriptionPairs[cnt+15].i]);
+
+            __m128i operand2 = _mm_set_epi8(pointsValue[descriptionPairs[cnt+0].j],
+                                            pointsValue[descriptionPairs[cnt+1].j],
+                                            pointsValue[descriptionPairs[cnt+2].j],
+                                            pointsValue[descriptionPairs[cnt+3].j],
+                                            pointsValue[descriptionPairs[cnt+4].j],
+                                            pointsValue[descriptionPairs[cnt+5].j],
+                                            pointsValue[descriptionPairs[cnt+6].j],
+                                            pointsValue[descriptionPairs[cnt+7].j],
+                                            pointsValue[descriptionPairs[cnt+8].j],
+                                            pointsValue[descriptionPairs[cnt+9].j],
+                                            pointsValue[descriptionPairs[cnt+10].j],
+                                            pointsValue[descriptionPairs[cnt+11].j],
+                                            pointsValue[descriptionPairs[cnt+12].j],
+                                            pointsValue[descriptionPairs[cnt+13].j],
+                                            pointsValue[descriptionPairs[cnt+14].j],
+                                            pointsValue[descriptionPairs[cnt+15].j]);
+
+            __m128i workReg = _mm_min_epu8(operand1, operand2); // emulated "not less than" for 8-bit UNSIGNED integers
+            workReg = _mm_cmpeq_epi8(workReg, operand2);        // emulated "not less than" for 8-bit UNSIGNED integers
+
+            workReg = _mm_and_si128(_mm_set1_epi16(short(0x8080 >> m)), workReg); // merge the last 16 bits with the 128bits std::vector until full
+            result128 = _mm_or_si128(result128, workReg);
+        }
+        (**ptrSSE) = result128;
+        ++(*ptrSSE);
+    }
+    (*ptrSSE) -= 8;
+}
+#endif
+
+template <typename srcMatType, typename iiMatType>
+void FREAK::computeDescriptors( const Mat& image, std::vector<KeyPoint>& keypoints, Mat& descriptors ) const {
+
     Mat imgIntegral;
-    integral(image, imgIntegral);
+    integral(image, imgIntegral, DataType<iiMatType>::type);
     std::vector<int> kpScaleIdx(keypoints.size()); // used to save pattern scale index corresponding to each keypoints
     const std::vector<int>::iterator ScaleIdxBegin = kpScaleIdx.begin(); // used in std::vector erase function
     const std::vector<cv::KeyPoint>::iterator kpBegin = keypoints.begin(); // used in std::vector erase function
     const float sizeCst = static_cast<float>(FREAK_NB_SCALES/(FREAK_LOG2* nOctaves));
-    uchar pointsValue[FREAK_NB_POINTS];
+    srcMatType pointsValue[FREAK_NB_POINTS];
     int thetaIdx = 0;
     int direction0;
     int direction1;
@@ -275,11 +390,8 @@ void FREAK::computeImpl( const Mat& image, std::vector<KeyPoint>& keypoints, Mat
     if( !extAll ) {
         // extract the best comparisons only
         descriptors = cv::Mat::zeros((int)keypoints.size(), FREAK_NB_PAIRS/8, CV_8U);
-#if CV_SSE2
-        __m128i* ptr= (__m128i*) (descriptors.data+(keypoints.size()-1)*descriptors.step[0]);
-#else
-        std::bitset<FREAK_NB_PAIRS>* ptr = (std::bitset<FREAK_NB_PAIRS>*) (descriptors.data+(keypoints.size()-1)*descriptors.step[0]);
-#endif
+        void *ptr = descriptors.data+(keypoints.size()-1)*descriptors.step[0];
+
         for( size_t k = keypoints.size(); k--; ) {
             // estimate orientation (gradient)
             if( !orientationNormalized ) {
@@ -289,7 +401,9 @@ void FREAK::computeImpl( const Mat& image, std::vector<KeyPoint>& keypoints, Mat
             else {
                 // get the points intensity value in the un-rotated pattern
                 for( int i = FREAK_NB_POINTS; i--; ) {
-                    pointsValue[i] = meanIntensity(image, imgIntegral, keypoints[k].pt.x,keypoints[k].pt.y, kpScaleIdx[k], 0, i);
+                    pointsValue[i] = meanIntensity<srcMatType, iiMatType>(image, imgIntegral,
+                                                                          keypoints[k].pt.x, keypoints[k].pt.y,
+                                                                          kpScaleIdx[k], 0, i);
                 }
                 direction0 = 0;
                 direction1 = 0;
@@ -310,78 +424,13 @@ void FREAK::computeImpl( const Mat& image, std::vector<KeyPoint>& keypoints, Mat
             }
             // extract descriptor at the computed orientation
             for( int i = FREAK_NB_POINTS; i--; ) {
-                pointsValue[i] = meanIntensity(image, imgIntegral, keypoints[k].pt.x,keypoints[k].pt.y, kpScaleIdx[k], thetaIdx, i);
+                pointsValue[i] = meanIntensity<srcMatType, iiMatType>(image, imgIntegral,
+                                                                      keypoints[k].pt.x, keypoints[k].pt.y,
+                                                                      kpScaleIdx[k], thetaIdx, i);
             }
-#if CV_SSE2
-            // 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)
-            int cnt = 0;
-            for( int n = FREAK_NB_PAIRS/128; n-- ; )
-            {
-                __m128i result128 = _mm_setzero_si128();
-                for( int m = 128/16; m--; cnt += 16 )
-                {
-                    __m128i operand1 = _mm_set_epi8(
-                        pointsValue[descriptionPairs[cnt+0].i],
-                        pointsValue[descriptionPairs[cnt+1].i],
-                        pointsValue[descriptionPairs[cnt+2].i],
-                        pointsValue[descriptionPairs[cnt+3].i],
-                        pointsValue[descriptionPairs[cnt+4].i],
-                        pointsValue[descriptionPairs[cnt+5].i],
-                        pointsValue[descriptionPairs[cnt+6].i],
-                        pointsValue[descriptionPairs[cnt+7].i],
-                        pointsValue[descriptionPairs[cnt+8].i],
-                        pointsValue[descriptionPairs[cnt+9].i],
-                        pointsValue[descriptionPairs[cnt+10].i],
-                        pointsValue[descriptionPairs[cnt+11].i],
-                        pointsValue[descriptionPairs[cnt+12].i],
-                        pointsValue[descriptionPairs[cnt+13].i],
-                        pointsValue[descriptionPairs[cnt+14].i],
-                        pointsValue[descriptionPairs[cnt+15].i]);
 
-                    __m128i operand2 = _mm_set_epi8(
-                        pointsValue[descriptionPairs[cnt+0].j],
-                        pointsValue[descriptionPairs[cnt+1].j],
-                        pointsValue[descriptionPairs[cnt+2].j],
-                        pointsValue[descriptionPairs[cnt+3].j],
-                        pointsValue[descriptionPairs[cnt+4].j],
-                        pointsValue[descriptionPairs[cnt+5].j],
-                        pointsValue[descriptionPairs[cnt+6].j],
-                        pointsValue[descriptionPairs[cnt+7].j],
-                        pointsValue[descriptionPairs[cnt+8].j],
-                        pointsValue[descriptionPairs[cnt+9].j],
-                        pointsValue[descriptionPairs[cnt+10].j],
-                        pointsValue[descriptionPairs[cnt+11].j],
-                        pointsValue[descriptionPairs[cnt+12].j],
-                        pointsValue[descriptionPairs[cnt+13].j],
-                        pointsValue[descriptionPairs[cnt+14].j],
-                        pointsValue[descriptionPairs[cnt+15].j]);
-
-                    __m128i workReg = _mm_min_epu8(operand1, operand2); // emulated "not less than" for 8-bit UNSIGNED integers
-                    workReg = _mm_cmpeq_epi8(workReg, operand2);        // emulated "not less than" for 8-bit UNSIGNED integers
-
-                    workReg = _mm_and_si128(_mm_set1_epi16(short(0x8080 >> m)), workReg); // merge the last 16 bits with the 128bits std::vector until full
-                    result128 = _mm_or_si128(result128, workReg);
-                }
-                (*ptr) = result128;
-                ++ptr;
-            }
-            ptr -= 8;
-#else
-            // extracting descriptor preserving the order of SSE version
-            int cnt = 0;
-            for( int n = 7; n < FREAK_NB_PAIRS; n += 128)
-            {
-                for( int m = 8; m--; )
-                {
-                    int nm = n-m;
-                    for(int kk = nm+15*8; kk >= nm; kk-=8, ++cnt)
-                    {
-                        ptr->set(kk, pointsValue[descriptionPairs[cnt].i] >= pointsValue[descriptionPairs[cnt].j]);
-                    }
-                }
-            }
-            --ptr;
-#endif
+            // Extract descriptor
+            extractDescriptor<srcMatType>(pointsValue, &ptr);
         }
     }
     else { // extract all possible comparisons for selection
@@ -397,7 +446,9 @@ void FREAK::computeImpl( const Mat& image, std::vector<KeyPoint>& keypoints, Mat
             else {
                 //get the points intensity value in the un-rotated pattern
                 for( int i = FREAK_NB_POINTS;i--; )
-                    pointsValue[i] = meanIntensity(image, imgIntegral, keypoints[k].pt.x,keypoints[k].pt.y, kpScaleIdx[k], 0, i);
+                    pointsValue[i] = meanIntensity<srcMatType, iiMatType>(image, imgIntegral,
+                                                                          keypoints[k].pt.x,keypoints[k].pt.y,
+                                                                          kpScaleIdx[k], 0, i);
 
                 direction0 = 0;
                 direction1 = 0;
@@ -419,8 +470,9 @@ void FREAK::computeImpl( const Mat& image, std::vector<KeyPoint>& keypoints, Mat
             }
             // get the points intensity value in the rotated pattern
             for( int i = FREAK_NB_POINTS; i--; ) {
-                pointsValue[i] = meanIntensity(image, imgIntegral, keypoints[k].pt.x,
-                                             keypoints[k].pt.y, kpScaleIdx[k], thetaIdx, i);
+                pointsValue[i] = meanIntensity<srcMatType, iiMatType>(image, imgIntegral,
+                                                                      keypoints[k].pt.x, keypoints[k].pt.y,
+                                                                      kpScaleIdx[k], thetaIdx, i);
             }
 
             int cnt(0);
@@ -437,19 +489,19 @@ void FREAK::computeImpl( const Mat& image, std::vector<KeyPoint>& keypoints, Mat
 }
 
 // simply take average on a square patch, not even gaussian approx
-uchar FREAK::meanIntensity( const cv::Mat& image, const cv::Mat& integral,
-                            const float kp_x,
-                            const float kp_y,
-                            const unsigned int scale,
-                            const unsigned int rot,
-                            const unsigned int point) const {
+template <typename imgType, typename iiType>
+imgType FREAK::meanIntensity( const cv::Mat& image, const cv::Mat& integral,
+                              const float kp_x,
+                              const float kp_y,
+                              const unsigned int scale,
+                              const unsigned int rot,
+                              const unsigned int point) const {
     // get point position in image
     const PatternPoint& FreakPoint = patternLookup[scale*FREAK_NB_ORIENTATION*FREAK_NB_POINTS + rot*FREAK_NB_POINTS + point];
     const float xf = FreakPoint.x+kp_x;
     const float yf = FreakPoint.y+kp_y;
     const int x = int(xf);
     const int y = int(yf);
-    const int& imagecols = image.cols;
 
     // get the sigma:
     const float radius = FreakPoint.sigma;
@@ -461,19 +513,15 @@ uchar FREAK::meanIntensity( const cv::Mat& image, const cv::Mat& integral,
         const int r_y = static_cast<int>((yf-y)*1024);
         const int r_x_1 = (1024-r_x);
         const int r_y_1 = (1024-r_y);
-        uchar* ptr = image.data+x+y*imagecols;
         unsigned int ret_val;
         // linear interpolation:
-        ret_val = (r_x_1*r_y_1*int(*ptr));
-        ptr++;
-        ret_val += (r_x*r_y_1*int(*ptr));
-        ptr += imagecols;
-        ret_val += (r_x*r_y*int(*ptr));
-        ptr--;
-        ret_val += (r_x_1*r_y*int(*ptr));
+        ret_val = r_x_1*r_y_1*int(image.at<imgType>(y  , x  ))
+                + r_x  *r_y_1*int(image.at<imgType>(y  , x+1))
+                + r_x_1*r_y  *int(image.at<imgType>(y+1, x  ))
+                + r_x  *r_y  *int(image.at<imgType>(y+1, x+1));
         //return the rounded mean
         ret_val += 2 * 1024 * 1024;
-        return static_cast<uchar>(ret_val / (4 * 1024 * 1024));
+        return static_cast<imgType>(ret_val / (4 * 1024 * 1024));
     }
 
     // expected case:
@@ -483,15 +531,15 @@ uchar FREAK::meanIntensity( const cv::Mat& image, const cv::Mat& integral,
     const int y_top = int(yf-radius+0.5);
     const int x_right = int(xf+radius+1.5);//integral image is 1px wider
     const int y_bottom = int(yf+radius+1.5);//integral image is 1px higher
-    int ret_val;
+    iiType ret_val;
 
-    ret_val = integral.at<int>(y_bottom,x_right);//bottom right corner
-    ret_val -= integral.at<int>(y_bottom,x_left);
-    ret_val += integral.at<int>(y_top,x_left);
-    ret_val -= integral.at<int>(y_top,x_right);
+    ret_val = integral.at<iiType>(y_bottom,x_right);//bottom right corner
+    ret_val -= integral.at<iiType>(y_bottom,x_left);
+    ret_val += integral.at<iiType>(y_top,x_left);
+    ret_val -= integral.at<iiType>(y_top,x_right);
     ret_val = ret_val/( (x_right-x_left)* (y_bottom-y_top) );
     //~ std::cout<<integral.step[1]<<std::endl;
-    return static_cast<uchar>(ret_val);
+    return static_cast<imgType>(ret_val);
 }
 
 // pair selection algorithm from a set of training images and corresponding keypoints