different interpolation by double image (#23124)

* different interpolation by double image

* fixing scaling mapping

* fixing a test

* added an option to enable previous interpolation

* added doxygen entries for the new parameter

* ASSERT_TRUE -> ASSERT_EQ

* changed log message when using old upscale mode

modules/features2d/include/opencv2/features2d.hpp

@@ -286,10 +286,14 @@ public:
 
     @param sigma The sigma of the Gaussian applied to the input image at the octave \#0. If your image
     is captured with a weak camera with soft lenses, you might want to reduce the number.
+
+    @param enable_precise_upscale Whether to enable precise upscaling in the scale pyramid, which maps
+    index \f$\texttt{x}\f$ to \f$\texttt{2x}\f$. This prevents localization bias. The option
+    to disable it (which is deprecated and issues a warning) is provided to keep the original behavior.
     */
     CV_WRAP static Ptr<SIFT> create(int nfeatures = 0, int nOctaveLayers = 3,
         double contrastThreshold = 0.04, double edgeThreshold = 10,
-        double sigma = 1.6);
+        double sigma = 1.6, bool enable_precise_upscale = true);
 
     /** @brief Create SIFT with specified descriptorType.
     @param nfeatures The number of best features to retain. The features are ranked by their scores
@@ -313,10 +317,14 @@ public:
     is captured with a weak camera with soft lenses, you might want to reduce the number.
 
     @param descriptorType The type of descriptors. Only CV_32F and CV_8U are supported.
+
+    @param enable_precise_upscale Whether to enable precise upscaling in the scale pyramid, which maps
+    index \f$\texttt{x}\f$ to \f$\texttt{2x}\f$. This prevents localization bias. The option
+    to disable it (which is deprecated and issues a warning) is provided to keep the original behavior.
     */
     CV_WRAP static Ptr<SIFT> create(int nfeatures, int nOctaveLayers,
         double contrastThreshold, double edgeThreshold,
-        double sigma, int descriptorType);
+        double sigma, int descriptorType, bool enable_precise_upscale = true);
 
     CV_WRAP virtual String getDefaultName() const CV_OVERRIDE;
 

modules/features2d/src/sift.dispatch.cpp

@@ -72,6 +72,7 @@
 #include "precomp.hpp"
 #include <opencv2/core/hal/hal.hpp>
 #include <opencv2/core/utils/tls.hpp>
+#include <opencv2/core/utils/logger.hpp>
 
 #include "sift.simd.hpp"
 #include "sift.simd_declarations.hpp" // defines CV_CPU_DISPATCH_MODES_ALL=AVX2,...,BASELINE based on CMakeLists.txt content
@@ -88,7 +89,8 @@ class SIFT_Impl : public SIFT
 public:
     explicit SIFT_Impl( int nfeatures = 0, int nOctaveLayers = 3,
                           double contrastThreshold = 0.04, double edgeThreshold = 10,
-                          double sigma = 1.6, int descriptorType = CV_32F );
+                          double sigma = 1.6, int descriptorType = CV_32F,
+                          bool enable_precise_upscale = true );
 
     //! returns the descriptor size in floats (128)
     int descriptorSize() const CV_OVERRIDE;
@@ -136,24 +138,25 @@ protected:
     CV_PROP_RW double edgeThreshold;
     CV_PROP_RW double sigma;
     CV_PROP_RW int descriptor_type;
+    CV_PROP_RW bool enable_precise_upscale;
 };
 
 Ptr<SIFT> SIFT::create( int _nfeatures, int _nOctaveLayers,
-                     double _contrastThreshold, double _edgeThreshold, double _sigma )
+                     double _contrastThreshold, double _edgeThreshold, double _sigma, bool enable_precise_upscale )
 {
     CV_TRACE_FUNCTION();
 
-    return makePtr<SIFT_Impl>(_nfeatures, _nOctaveLayers, _contrastThreshold, _edgeThreshold, _sigma, CV_32F);
+    return makePtr<SIFT_Impl>(_nfeatures, _nOctaveLayers, _contrastThreshold, _edgeThreshold, _sigma, CV_32F, enable_precise_upscale);
 }
 
 Ptr<SIFT> SIFT::create( int _nfeatures, int _nOctaveLayers,
-                     double _contrastThreshold, double _edgeThreshold, double _sigma, int _descriptorType )
+                     double _contrastThreshold, double _edgeThreshold, double _sigma, int _descriptorType, bool enable_precise_upscale )
 {
     CV_TRACE_FUNCTION();
 
     // SIFT descriptor supports 32bit floating point and 8bit unsigned int.
     CV_Assert(_descriptorType == CV_32F || _descriptorType == CV_8U);
-    return makePtr<SIFT_Impl>(_nfeatures, _nOctaveLayers, _contrastThreshold, _edgeThreshold, _sigma, _descriptorType);
+    return makePtr<SIFT_Impl>(_nfeatures, _nOctaveLayers, _contrastThreshold, _edgeThreshold, _sigma, _descriptorType, enable_precise_upscale);
 }
 
 String SIFT::getDefaultName() const
@@ -170,7 +173,7 @@ unpackOctave(const KeyPoint& kpt, int& octave, int& layer, float& scale)
     scale = octave >= 0 ? 1.f/(1 << octave) : (float)(1 << -octave);
 }
 
-static Mat createInitialImage( const Mat& img, bool doubleImageSize, float sigma )
+static Mat createInitialImage( const Mat& img, bool doubleImageSize, float sigma, bool enable_precise_upscale )
 {
     CV_TRACE_FUNCTION();
 
@@ -188,12 +191,22 @@ static Mat createInitialImage( const Mat& img, bool doubleImageSize, float sigma
     if( doubleImageSize )
     {
         sig_diff = sqrtf( std::max(sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA * 4, 0.01f) );
+
         Mat dbl;
+        if (enable_precise_upscale) {
+            dbl.create(Size(gray_fpt.cols*2, gray_fpt.rows*2), gray_fpt.type());
+            Mat H = Mat::zeros(2, 3, CV_32F);
+            H.at<float>(0, 0) = 0.5f;
+            H.at<float>(1, 1) = 0.5f;
+
+            cv::warpAffine(gray_fpt, dbl, H, dbl.size(), INTER_LINEAR | WARP_INVERSE_MAP, BORDER_REFLECT);
+        } else {
 #if DoG_TYPE_SHORT
-        resize(gray_fpt, dbl, Size(gray_fpt.cols*2, gray_fpt.rows*2), 0, 0, INTER_LINEAR_EXACT);
+            resize(gray_fpt, dbl, Size(gray_fpt.cols*2, gray_fpt.rows*2), 0, 0, INTER_LINEAR_EXACT);
 #else
-        resize(gray_fpt, dbl, Size(gray_fpt.cols*2, gray_fpt.rows*2), 0, 0, INTER_LINEAR);
+            resize(gray_fpt, dbl, Size(gray_fpt.cols*2, gray_fpt.rows*2), 0, 0, INTER_LINEAR);
 #endif
+        }
         Mat result;
         GaussianBlur(dbl, result, Size(), sig_diff, sig_diff);
         return result;
@@ -459,10 +472,14 @@ static void calcDescriptors(const std::vector<Mat>& gpyr, const std::vector<KeyP
 //////////////////////////////////////////////////////////////////////////////////////////
 
 SIFT_Impl::SIFT_Impl( int _nfeatures, int _nOctaveLayers,
-           double _contrastThreshold, double _edgeThreshold, double _sigma, int _descriptorType )
+           double _contrastThreshold, double _edgeThreshold, double _sigma, int _descriptorType, bool _enable_precise_upscale)
     : nfeatures(_nfeatures), nOctaveLayers(_nOctaveLayers),
-    contrastThreshold(_contrastThreshold), edgeThreshold(_edgeThreshold), sigma(_sigma), descriptor_type(_descriptorType)
+    contrastThreshold(_contrastThreshold), edgeThreshold(_edgeThreshold), sigma(_sigma), descriptor_type(_descriptorType),
+    enable_precise_upscale(_enable_precise_upscale)
 {
+    if (!enable_precise_upscale) {
+        CV_LOG_ONCE_INFO(NULL, "precise upscale disabled, this is now deprecated as it was found to induce a location bias");
+    }
 }
 
 int SIFT_Impl::descriptorSize() const
@@ -516,7 +533,7 @@ void SIFT_Impl::detectAndCompute(InputArray _image, InputArray _mask,
         actualNOctaves = maxOctave - firstOctave + 1;
     }
 
-    Mat base = createInitialImage(image, firstOctave < 0, (float)sigma);
+    Mat base = createInitialImage(image, firstOctave < 0, (float)sigma, enable_precise_upscale);
     std::vector<Mat> gpyr;
     int nOctaves = actualNOctaves > 0 ? actualNOctaves : cvRound(std::log( (double)std::min( base.cols, base.rows ) ) / std::log(2.) - 2) - firstOctave;
 

modules/features2d/test/test_descriptors_regression.impl.hpp

@@ -7,6 +7,34 @@ namespace opencv_test { namespace {
 /****************************************************************************************\
 *                     Regression tests for descriptor extractors.                        *
 \****************************************************************************************/
+static void double_image(Mat& src, Mat& dst) {
+
+    dst.create(Size(src.cols*2, src.rows*2), src.type());
+
+    Mat H = Mat::zeros(2, 3, CV_32F);
+    H.at<float>(0, 0) = 0.5f;
+    H.at<float>(1, 1) = 0.5f;
+    cv::warpAffine(src, dst, H, dst.size(), INTER_LINEAR | WARP_INVERSE_MAP, BORDER_REFLECT);
+
+}
+
+static Mat prepare_img(bool rows_indexed) {
+    int rows = 5;
+    int columns = 5;
+    Mat img(rows, columns, CV_32F);
+
+    for (int i = 0; i < rows; i++) {
+        for (int j = 0; j < columns; j++) {
+            if (rows_indexed) {
+                img.at<float>(i, j) = (float)i;
+            } else {
+                img.at<float>(i, j) = (float)j;
+            }
+        }
+    }
+    return img;
+}
+
 static void writeMatInBin( const Mat& mat, const string& filename )
 {
     FILE* f = fopen( filename.c_str(), "wb");
@@ -145,6 +173,25 @@ protected:
             ts->set_failed_test_info( cvtest::TS::FAIL_INVALID_TEST_DATA );
         }
 
+        image = prepare_img(false);
+        Mat dbl;
+        try
+        {
+            double_image(image, dbl);
+
+            Mat downsized_back(dbl.rows/2, dbl.cols/2, CV_32F);
+            resize(dbl, downsized_back, Size(dbl.cols/2, dbl.rows/2), 0, 0, INTER_NEAREST);
+
+            cv::Mat diff = (image != downsized_back);
+            ASSERT_EQ(0, cv::norm(image, downsized_back, NORM_INF));
+        }
+        catch(...)
+        {
+            ts->printf( cvtest::TS::LOG, "double_image() must not generate exception (1).\n");
+            ts->printf( cvtest::TS::LOG, "double_image() when downsized back by NEAREST must generate the same original image (1).\n");
+            ts->set_failed_test_info( cvtest::TS::FAIL_INVALID_TEST_DATA );
+        }
+
         // Several images.
         vector<Mat> images;
         vector<vector<KeyPoint> > keypointsCollection;

modules/features2d/test/test_detectors_invariance.cpp

@@ -37,7 +37,7 @@ INSTANTIATE_TEST_CASE_P(AKAZE_DESCRIPTOR_KAZE, DetectorRotationInvariance,
  */
 
 INSTANTIATE_TEST_CASE_P(SIFT, DetectorScaleInvariance,
-                        Value(IMAGE_BIKES, SIFT::create(0, 3, 0.09), 0.65f, 0.98f));
+                        Value(IMAGE_BIKES, SIFT::create(0, 3, 0.09), 0.60f, 0.98f));
 
 INSTANTIATE_TEST_CASE_P(BRISK, DetectorScaleInvariance,
                         Value(IMAGE_BIKES, BRISK::create(), 0.08f, 0.49f));

modules/features2d/test/test_detectors_invariance.impl.hpp

@@ -25,7 +25,6 @@ void matchKeyPoints(const vector<KeyPoint>& keypoints0, const Mat& H,
         perspectiveTransform(Mat(points0), points0t, H);
 
     matches.clear();
-    vector<uchar> usedMask(keypoints1.size(), 0);
     for(int i0 = 0; i0 < static_cast<int>(keypoints0.size()); i0++)
     {
         int nearestPointIndex = -1;
@@ -33,8 +32,6 @@ void matchKeyPoints(const vector<KeyPoint>& keypoints0, const Mat& H,
         const float r0 =  0.5f * keypoints0[i0].size;
         for(size_t i1 = 0; i1 < keypoints1.size(); i1++)
         {
-            if(nearestPointIndex >= 0 && usedMask[i1])
-                continue;
 
             float r1 = 0.5f * keypoints1[i1].size;
             float intersectRatio = calcIntersectRatio(points0t.at<Point2f>(i0), r0,
@@ -47,8 +44,6 @@ void matchKeyPoints(const vector<KeyPoint>& keypoints0, const Mat& H,
         }
 
         matches.push_back(DMatch(i0, nearestPointIndex, maxIntersectRatio));
-        if(nearestPointIndex >= 0)
-            usedMask[nearestPointIndex] = 1;
     }
 }
 

modules/features2d/test/test_invariance_utils.hpp

@@ -75,8 +75,8 @@ void scaleKeyPoints(const vector<KeyPoint>& src, vector<KeyPoint>& dst, float sc
     dst.resize(src.size());
     for (size_t i = 0; i < src.size(); i++) {
         dst[i] = src[i];
-        dst[i].pt.x *= scale;
-        dst[i].pt.y *= scale;
+        dst[i].pt.x = dst[i].pt.x * scale + (scale - 1.0f) / 2.0f;
+        dst[i].pt.y = dst[i].pt.y * scale + (scale - 1.0f) / 2.0f;
         dst[i].size *= scale;
     }
 }