added upright mode to SURF (#825)

modules/features2d/include/opencv2/features2d/features2d.hpp

@@ -81,6 +81,7 @@ CV_INLINE CvSURFPoint cvSURFPoint( CvPoint2D32f pt, int laplacian,
 typedef struct CvSURFParams
 {
     int    extended;
+    int    upright;
     double hessianThreshold;
 
     int    nOctaves;
@@ -395,7 +396,7 @@ public:
     CV_WRAP SURF();
     //! the full constructor taking all the necessary parameters
     CV_WRAP SURF(double _hessianThreshold, int _nOctaves=4,
-         int _nOctaveLayers=2, bool _extended=false);
+         int _nOctaveLayers=2, bool _extended=false, bool _upright=false);
 
     //! returns the descriptor size in float's (64 or 128)
     CV_WRAP int descriptorSize() const;
@@ -1519,7 +1520,7 @@ protected:
 class CV_EXPORTS SurfFeatureDetector : public FeatureDetector
 {
 public:
-    SurfFeatureDetector( double hessianThreshold=400., int octaves=3, int octaveLayers=4 );
+    SurfFeatureDetector( double hessianThreshold=400., int octaves=3, int octaveLayers=4, bool upright=false );
     virtual void read( const FileNode& fn );
     virtual void write( FileStorage& fs ) const;
 
@@ -1897,7 +1898,7 @@ protected:
 class CV_EXPORTS SurfDescriptorExtractor : public DescriptorExtractor
 {
 public:
-    SurfDescriptorExtractor( int nOctaves=4, int nOctaveLayers=2, bool extended=false );
+    SurfDescriptorExtractor( int nOctaves=4, int nOctaveLayers=2, bool extended=false, bool upright=false );
 
     virtual void read( const FileNode &fn );
     virtual void write( FileStorage &fs ) const;
@@ -1906,7 +1907,7 @@ public:
     virtual int descriptorType() const;
 
 protected:
-	virtual void computeImpl( const Mat& image, vector<KeyPoint>& keypoints, Mat& descriptors ) const;
+    virtual void computeImpl( const Mat& image, vector<KeyPoint>& keypoints, Mat& descriptors ) const;
 
     SURF surf;
 };

modules/features2d/src/descriptors.cpp

@@ -191,8 +191,8 @@ int SiftDescriptorExtractor::descriptorType() const
 *                                SurfDescriptorExtractor                                 *
 \****************************************************************************************/
 SurfDescriptorExtractor::SurfDescriptorExtractor( int nOctaves,
-                                                  int nOctaveLayers, bool extended )
+                                                  int nOctaveLayers, bool extended, bool upright )
-    : surf( 0.0, nOctaves, nOctaveLayers, extended )
+    : surf( 0.0, nOctaves, nOctaveLayers, extended, upright )
 {}
 
 void SurfDescriptorExtractor::computeImpl( const Mat& image,
@@ -218,8 +218,9 @@ void SurfDescriptorExtractor::read( const FileNode &fn )
     int nOctaves = fn["nOctaves"];
     int nOctaveLayers = fn["nOctaveLayers"];
     bool extended = (int)fn["extended"] != 0;
+    bool upright = (int)fn["upright"] != 0;
 
-    surf = SURF( 0.0, nOctaves, nOctaveLayers, extended );
+    surf = SURF( 0.0, nOctaves, nOctaveLayers, extended, upright );
 }
 
 void SurfDescriptorExtractor::write( FileStorage &fs ) const
@@ -229,6 +230,7 @@ void SurfDescriptorExtractor::write( FileStorage &fs ) const
     fs << "nOctaves" << surf.nOctaves;
     fs << "nOctaveLayers" << surf.nOctaveLayers;
     fs << "extended" << surf.extended;
+    fs << "upright" << surf.upright;
 }
 
 int SurfDescriptorExtractor::descriptorSize() const

modules/features2d/src/detectors.cpp

@@ -407,8 +407,8 @@ void SiftFeatureDetector::detectImpl( const Mat& image, vector<KeyPoint>& keypoi
 /*
  *  SurfFeatureDetector
  */
-SurfFeatureDetector::SurfFeatureDetector( double hessianThreshold, int octaves, int octaveLayers)
+SurfFeatureDetector::SurfFeatureDetector( double hessianThreshold, int octaves, int octaveLayers, bool upright )
-    : surf(hessianThreshold, octaves, octaveLayers)
+    : surf(hessianThreshold, octaves, octaveLayers, false, upright)
 {}
 
 void SurfFeatureDetector::read (const FileNode& fn)
@@ -416,8 +416,9 @@ void SurfFeatureDetector::read (const FileNode& fn)
     double hessianThreshold = fn["hessianThreshold"];
     int octaves = fn["octaves"];
     int octaveLayers = fn["octaveLayers"];
+    bool upright = (int)fn["upright"] != 0;
 
-    surf = SURF( hessianThreshold, octaves, octaveLayers );
+    surf = SURF( hessianThreshold, octaves, octaveLayers, false, upright );
 }
 
 void SurfFeatureDetector::write (FileStorage& fs) const
@@ -427,6 +428,7 @@ void SurfFeatureDetector::write (FileStorage& fs) const
     fs << "hessianThreshold" << surf.hessianThreshold;
     fs << "octaves" << surf.nOctaves;
     fs << "octaveLayers" << surf.nOctaveLayers;
+    fs << "upright" << surf.upright;
 }
 
 void SurfFeatureDetector::detectImpl( const Mat& image, vector<KeyPoint>& keypoints, const Mat& mask ) const

modules/features2d/src/surf.cpp

@@ -114,6 +114,7 @@ CvSURFParams cvSURFParams(double threshold, int extended)
     CvSURFParams params;
     params.hessianThreshold = threshold;
     params.extended = extended;
+    params.upright = 0;
     params.nOctaves = 4;
     params.nOctaveLayers = 2;
     return params;
@@ -630,94 +631,130 @@ struct SURFInvoker
                 kp->size = -1;
                 continue;
             }
-            icvResizeHaarPattern( dx_s, dx_t, NX, 4, grad_wav_size, sum->cols );
-            icvResizeHaarPattern( dy_s, dy_t, NY, 4, grad_wav_size, sum->cols );
-            for( kk = 0, nangle = 0; kk < nOriSamples; kk++ )
-            {
-                const int* ptr;
-                float vx, vy;
-                x = cvRound( center.x + apt[kk].x*s - (float)(grad_wav_size-1)/2 );
-                y = cvRound( center.y + apt[kk].y*s - (float)(grad_wav_size-1)/2 );
-                if( (unsigned)y >= (unsigned)(sum->rows - grad_wav_size) ||
-                   (unsigned)x >= (unsigned)(sum->cols - grad_wav_size) )
-                    continue;
-                ptr = sum_ptr + x + y*sum_cols;
-                vx = icvCalcHaarPattern( ptr, dx_t, 2 );
-                vy = icvCalcHaarPattern( ptr, dy_t, 2 );
-                X[nangle] = vx*aptw[kk]; Y[nangle] = vy*aptw[kk];
-                nangle++;
-            }
-            if ( nangle == 0 )
-            {
-                /* No gradient could be sampled because the keypoint is too
-                 * near too one or more of the sides of the image. As we
-                 * therefore cannot find a dominant direction, we skip this
-                 * keypoint and mark it for later deletion from the sequence. */
-                kp->size = -1;
-                continue;
-            }
-            matX.cols = matY.cols = _angle.cols = nangle;
-            cvCartToPolar( &matX, &matY, 0, &_angle, 1 );
 
-            float bestx = 0, besty = 0, descriptor_mod = 0;
+            float descriptor_dir = 90.f;
-            for( i = 0; i < 360; i += ORI_SEARCH_INC )
+            if (params->upright == 0)
             {
-                float sumx = 0, sumy = 0, temp_mod;
+                icvResizeHaarPattern( dx_s, dx_t, NX, 4, grad_wav_size, sum->cols );
-                for( j = 0; j < nangle; j++ )
+                icvResizeHaarPattern( dy_s, dy_t, NY, 4, grad_wav_size, sum->cols );
+                for( kk = 0, nangle = 0; kk < nOriSamples; kk++ )
                 {
-                    int d = std::abs(cvRound(angle[j]) - i);
+                    const int* ptr;
-                    if( d < ORI_WIN/2 || d > 360-ORI_WIN/2 )
+                    float vx, vy;
+                    x = cvRound( center.x + apt[kk].x*s - (float)(grad_wav_size-1)/2 );
+                    y = cvRound( center.y + apt[kk].y*s - (float)(grad_wav_size-1)/2 );
+                    if( (unsigned)y >= (unsigned)(sum->rows - grad_wav_size) ||
+                        (unsigned)x >= (unsigned)(sum->cols - grad_wav_size) )
+                        continue;
+                    ptr = sum_ptr + x + y*sum_cols;
+                    vx = icvCalcHaarPattern( ptr, dx_t, 2 );
+                    vy = icvCalcHaarPattern( ptr, dy_t, 2 );
+                    X[nangle] = vx*aptw[kk]; Y[nangle] = vy*aptw[kk];
+                    nangle++;
+                }
+                if ( nangle == 0 )
+                {
+                    /* No gradient could be sampled because the keypoint is too
+                    * near too one or more of the sides of the image. As we
+                    * therefore cannot find a dominant direction, we skip this
+                    * keypoint and mark it for later deletion from the sequence. */
+                    kp->size = -1;
+                    continue;
+                }
+                matX.cols = matY.cols = _angle.cols = nangle;
+                cvCartToPolar( &matX, &matY, 0, &_angle, 1 );
+
+                float bestx = 0, besty = 0, descriptor_mod = 0;
+                for( i = 0; i < 360; i += ORI_SEARCH_INC )
+                {
+                    float sumx = 0, sumy = 0, temp_mod;
+                    for( j = 0; j < nangle; j++ )
                     {
-                        sumx += X[j];
+                        int d = std::abs(cvRound(angle[j]) - i);
-                        sumy += Y[j];
+                        if( d < ORI_WIN/2 || d > 360-ORI_WIN/2 )
+                        {
+                            sumx += X[j];
+                            sumy += Y[j];
+                        }
+                    }
+                    temp_mod = sumx*sumx + sumy*sumy;
+                    if( temp_mod > descriptor_mod )
+                    {
+                        descriptor_mod = temp_mod;
+                        bestx = sumx;
+                        besty = sumy;
                     }
                 }
-                temp_mod = sumx*sumx + sumy*sumy;
+                descriptor_dir = cvFastArctan( besty, bestx );
-                if( temp_mod > descriptor_mod )
-                {
-                    descriptor_mod = temp_mod;
-                    bestx = sumx;
-                    besty = sumy;
-                }
             }
-            float descriptor_dir = cvFastArctan( besty, bestx );
             kp->dir = descriptor_dir;
             if( !descriptors )
                 continue;
-            descriptor_dir *= (float)(CV_PI/180);
+            
             /* Extract a window of pixels around the keypoint of size 20s */
             int win_size = (int)((PATCH_SZ+1)*s);
             CV_Assert( winbuf->cols >= win_size*win_size );
             CvMat win = cvMat(win_size, win_size, CV_8U, winbuf->data.ptr);
-            float sin_dir = sin(descriptor_dir);
-            float cos_dir = cos(descriptor_dir) ;
 
-            /* Subpixel interpolation version (slower). Subpixel not required since
+            if (params->upright == 0)
-             the pixels will all get averaged when we scale down to 20 pixels */
-            /*
-             float w[] = { cos_dir, sin_dir, center.x,
-             -sin_dir, cos_dir , center.y };
-             CvMat W = cvMat(2, 3, CV_32F, w);
-             cvGetQuadrangleSubPix( img, &win, &W );
-             */
-
-            /* Nearest neighbour version (faster) */
-            float win_offset = -(float)(win_size-1)/2;
-            float start_x = center.x + win_offset*cos_dir + win_offset*sin_dir;
-            float start_y = center.y - win_offset*sin_dir + win_offset*cos_dir;
-            uchar* WIN = win.data.ptr;
-            for( i = 0; i < win_size; i++, start_x += sin_dir, start_y += cos_dir )
             {
-                float pixel_x = start_x;
+            	descriptor_dir *= (float)(CV_PI/180);
-                float pixel_y = start_y;
+                float sin_dir = sin(descriptor_dir);
-                for( j = 0; j < win_size; j++, pixel_x += cos_dir, pixel_y -= sin_dir )
+                float cos_dir = cos(descriptor_dir);
+
+                /* Subpixel interpolation version (slower). Subpixel not required since
+                the pixels will all get averaged when we scale down to 20 pixels */
+                /*  
+                float w[] = { cos_dir, sin_dir, center.x,
+                -sin_dir, cos_dir , center.y };
+                CvMat W = cvMat(2, 3, CV_32F, w);
+                cvGetQuadrangleSubPix( img, &win, &W );
+                */
+
+                /* Nearest neighbour version (faster) */
+                float win_offset = -(float)(win_size-1)/2;
+                float start_x = center.x + win_offset*cos_dir + win_offset*sin_dir;
+                float start_y = center.y - win_offset*sin_dir + win_offset*cos_dir;
+                uchar* WIN = win.data.ptr;
+                for( i = 0; i < win_size; i++, start_x += sin_dir, start_y += cos_dir )
                 {
-                    int x = std::min(std::max(cvRound(pixel_x), 0), img->cols-1);
+                    float pixel_x = start_x;
-                    int y = std::min(std::max(cvRound(pixel_y), 0), img->rows-1);
+                    float pixel_y = start_y;
-                    WIN[i*win_size + j] = img->data.ptr[y*img->step + x];
+                    for( j = 0; j < win_size; j++, pixel_x += cos_dir, pixel_y -= sin_dir )
+                    {
+                        int x = std::min(std::max(cvRound(pixel_x), 0), img->cols-1);
+                        int y = std::min(std::max(cvRound(pixel_y), 0), img->rows-1);
+                        WIN[i*win_size + j] = img->data.ptr[y*img->step + x];
+                    }
                 }
             }
+            else
+            {
+                /* extract rect - slightly optimized version of the code above
+                   TODO: find faster code, as this is simply an extract rect operation, 
+                         e.g. by using cvGetSubRect, problem is the border processing */
+                // descriptor_dir == 90 grad
+                // sin_dir == 1
+                // cos_dir == 0
                 
+                float win_offset = -(float)(win_size-1)/2;
+                int start_x = cvRound(center.x + win_offset);
+                int start_y = cvRound(center.y - win_offset);
+                uchar* WIN = win.data.ptr;
+                for( i = 0; i < win_size; i++, start_x++ )
+                {
+                    int pixel_x = start_x;
+                    int pixel_y = start_y;
+                    for( j=0; j<win_size; j++, pixel_y-- )
+                    {
+                        x = MAX( pixel_x, 0 );
+                        y = MAX( pixel_y, 0 );
+                        x = MIN( x, img->cols-1 );
+                        y = MIN( y, img->rows-1 );
+                        WIN[i*win_size + j] = img->data.ptr[y*img->step+x];
+                    }
+                }               
+            }
             /* Scale the window to size PATCH_SZ so each pixel's size is s. This
              makes calculating the gradients with wavelets of size 2s easy */
             cvResize( &win, &_patch, CV_INTER_AREA );
@@ -886,8 +923,8 @@ cvExtractSURF( const CvArr* _img, const CvArr* _mask,
         cv::parallel_for(cv::BlockedRange(0, N),
                      cv::SURFInvoker(&params, keypoints, descriptors, img, sum) );
 #else
-	cv::SURFInvoker invoker(&params, keypoints, descriptors, img, sum);
+	    cv::SURFInvoker invoker(&params, keypoints, descriptors, img, sum);
-	invoker(cv::BlockedRange(0, N));
+	    invoker(cv::BlockedRange(0, N));
 #endif
     }
 
@@ -924,14 +961,16 @@ SURF::SURF()
 {
     hessianThreshold = 100;
     extended = 1;
+    upright = 0;
     nOctaves = 4;
     nOctaveLayers = 2;
 }
 
-SURF::SURF(double _threshold, int _nOctaves, int _nOctaveLayers, bool _extended)
+SURF::SURF(double _threshold, int _nOctaves, int _nOctaveLayers, bool _extended, bool _upright)
 {
     hessianThreshold = _threshold;
     extended = _extended;
+    upright = _upright;
     nOctaves = _nOctaves;
     nOctaveLayers = _nOctaveLayers;
 }