--- freak.old.cpp	2013-10-14 08:28:48.000000000 -0700
+++ freak.cpp	2013-10-15 22:38:17.000000000 -0700
@@ -214,23 +214,56 @@
              descriptionPairs[i] = allPairs[FREAK_DEF_PAIRS[i]];
     }
 }
-
+    
 void FREAK::computeImpl( const Mat& image, std::vector<KeyPoint>& keypoints, Mat& descriptors ) const {
-
+    
     if( image.empty() )
         return;
     if( keypoints.empty() )
         return;
-
+    
     ((FREAK*)this)->buildPattern();
+    
+    // Use 32-bit integers if we won't overflow in the integral image
+    if ((image.type() == CV_8U || image.type() == 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.type() == CV_8U)
+            computeDescriptors<uchar, int>(image, keypoints, descriptors);
+        else if (image.type() == CV_8S)
+            computeDescriptors<char, int>(image, keypoints, descriptors);
+        else
+            CV_Error( CV_StsUnsupportedFormat, "" );
+    } else {
+        // Create the integral image appropriate for our type & usage
+        if ( image.type() == CV_8U )
+            computeDescriptors<uchar, double>(image, keypoints, descriptors);
+        else if ( image.type() == CV_8S )
+            computeDescriptors<char, double>(image, keypoints, descriptors);
+        else if ( image.type() == CV_16U )
+            computeDescriptors<ushort, double>(image, keypoints, descriptors);
+        else if ( image.type() == CV_16S )
+            computeDescriptors<short, double>(image, keypoints, descriptors);
+        else
+            CV_Error( CV_StsUnsupportedFormat, "" );
+    }
+}
 
+template <typename srcMatType, typename iiMatType>
+void FREAK::computeDescriptors( const Mat& image, std::vector<KeyPoint>& keypoints, Mat& descriptors ) const {
+    
     Mat imgIntegral;
-    integral(image, imgIntegral);
+    if (sizeof(iiMatType) == sizeof(double))
+        integral(image, imgIntegral, CV_64F);
+    else
+        integral(image, imgIntegral, CV_32S);
+    
     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;
@@ -276,10 +309,10 @@
         // 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]);
+        __m128i* ptrSSE = (__m128i*) (descriptors.data+(keypoints.size()-1)*descriptors.step[0]);
 #endif
+        std::bitset<FREAK_NB_PAIRS>* ptrScalar = (std::bitset<FREAK_NB_PAIRS>*) (descriptors.data+(keypoints.size()-1)*descriptors.step[0]);
+        
         for( size_t k = keypoints.size(); k--; ) {
             // estimate orientation (gradient)
             if( !orientationNormalized ) {
@@ -289,7 +322,9 @@
             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 +345,83 @@
             }
             // 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 )
+            if (sizeof(srcMatType) == sizeof(char)) {
+                // 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 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 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
+                        __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);
+                        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;
                 }
-                (*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)
+                ptrSSE -= 8;
+            } else
+#endif
             {
-                for( int m = 8; m--; )
+                // extracting descriptor preserving the order of SSE version
+                int cnt = 0;
+                for( int n = 7; n < FREAK_NB_PAIRS; n += 128)
                 {
-                    int nm = n-m;
-                    for(int kk = nm+15*8; kk >= nm; kk-=8, ++cnt)
+                    for( int m = 8; m--; )
                     {
-                        ptr->set(kk, pointsValue[descriptionPairs[cnt].i] >= pointsValue[descriptionPairs[cnt].j]);
+                        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;
             }
-            --ptr;
-#endif
         }
     }
     else { // extract all possible comparisons for selection
@@ -397,7 +437,9 @@
             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 +461,9 @@
             }
             // 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,7 +480,8 @@
 }
 
 // simply take average on a square patch, not even gaussian approx
-uchar FREAK::meanIntensity( const cv::Mat& image, const cv::Mat& integral,
+template <typename imgType, typename iiType>
+int FREAK::meanIntensity( const cv::Mat& image, const cv::Mat& integral,
                             const float kp_x,
                             const float kp_y,
                             const unsigned int scale,
@@ -461,19 +505,15 @@
         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<int>(ret_val / (4 * 1024 * 1024));
     }
 
     // expected case:
@@ -483,15 +523,15 @@
     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<int>(ret_val);
 }
 
 // pair selection algorithm from a set of training images and corresponding keypoints
