Merge pull request #23305 from bhavitp:fix/calib3d/undistortion_grid

Resolves https://github.com/opencv/opencv/issues/23304

Fixes the incorrect pixel grid
Switches type to double to avoid precision loss as all callers use doubles

### Pull Request Readiness Checklist

See details at https://github.com/opencv/opencv/wiki/How_to_contribute#making-a-good-pull-request

- [X] I agree to contribute to the project under Apache 2 License.
- [X] To the best of my knowledge, the proposed patch is not based on a code under GPL or another license that is incompatible with OpenCV
- [X] The PR is proposed to the proper branch
- [X] There is a reference to the original bug report and related work
- [X] There is accuracy test, performance test and test data in opencv_extra repository, if applicable
      Patch to opencv_extra has the same branch name.
- [X] The feature is well documented and sample code can be built with the project CMake

modules/calib3d/src/calibration.cpp

@@ -2258,28 +2258,28 @@ double cvStereoCalibrate( const CvMat* _objectPoints, const CvMat* _imagePoints1
 static void
 icvGetRectangles( const CvMat* cameraMatrix, const CvMat* distCoeffs,
                  const CvMat* R, const CvMat* newCameraMatrix, CvSize imgSize,
-                 cv::Rect_<float>& inner, cv::Rect_<float>& outer )
+                 cv::Rect_<double>& inner, cv::Rect_<double>& outer )
 {
     const int N = 9;
     int x, y, k;
-    cv::Ptr<CvMat> _pts(cvCreateMat(1, N*N, CV_32FC2));
+    cv::Ptr<CvMat> _pts(cvCreateMat(1, N*N, CV_64FC2));
-    CvPoint2D32f* pts = (CvPoint2D32f*)(_pts->data.ptr);
+    CvPoint2D64f* pts = (CvPoint2D64f*)(_pts->data.ptr);
 
     for( y = k = 0; y < N; y++ )
         for( x = 0; x < N; x++ )
-            pts[k++] = cvPoint2D32f((float)x*imgSize.width/(N-1),
+            pts[k++] = cvPoint2D64f((double)x*(imgSize.width-1)/(N-1),
-                                    (float)y*imgSize.height/(N-1));
+                                    (double)y*(imgSize.height-1)/(N-1));
 
     cvUndistortPoints(_pts, _pts, cameraMatrix, distCoeffs, R, newCameraMatrix);
 
-    float iX0=-FLT_MAX, iX1=FLT_MAX, iY0=-FLT_MAX, iY1=FLT_MAX;
+    double iX0=-FLT_MAX, iX1=FLT_MAX, iY0=-FLT_MAX, iY1=FLT_MAX;
-    float oX0=FLT_MAX, oX1=-FLT_MAX, oY0=FLT_MAX, oY1=-FLT_MAX;
+    double oX0=FLT_MAX, oX1=-FLT_MAX, oY0=FLT_MAX, oY1=-FLT_MAX;
     // find the inscribed rectangle.
     // the code will likely not work with extreme rotation matrices (R) (>45%)
     for( y = k = 0; y < N; y++ )
         for( x = 0; x < N; x++ )
         {
-            CvPoint2D32f p = pts[k++];
+            CvPoint2D64f p = pts[k++];
             oX0 = MIN(oX0, p.x);
             oX1 = MAX(oX1, p.x);
             oY0 = MIN(oY0, p.y);
@@ -2294,8 +2294,8 @@ icvGetRectangles( const CvMat* cameraMatrix, const CvMat* distCoeffs,
             if( y == N-1 )
                 iY1 = MIN(iY1, p.y);
         }
-    inner = cv::Rect_<float>(iX0, iY0, iX1-iX0, iY1-iY0);
+    inner = cv::Rect_<double>(iX0, iY0, iX1-iX0, iY1-iY0);
-    outer = cv::Rect_<float>(oX0, oY0, oX1-oX0, oY1-oY0);
+    outer = cv::Rect_<double>(oX0, oY0, oX1-oX0, oY1-oY0);
 }
 
 
@@ -2308,7 +2308,7 @@ void cvStereoRectify( const CvMat* _cameraMatrix1, const CvMat* _cameraMatrix2,
 {
     double _om[3], _t[3] = {0}, _uu[3]={0,0,0}, _r_r[3][3], _pp[3][4];
     double _ww[3], _wr[3][3], _z[3] = {0,0,0}, _ri[3][3];
-    cv::Rect_<float> inner1, inner2, outer1, outer2;
+    cv::Rect_<double> inner1, inner2, outer1, outer2;
 
     CvMat om  = cvMat(3, 1, CV_64F, _om);
     CvMat t   = cvMat(3, 1, CV_64F, _t);
@@ -2515,7 +2515,7 @@ void cvGetOptimalNewCameraMatrix( const CvMat* cameraMatrix, const CvMat* distCo
                                   CvMat* newCameraMatrix, CvSize newImgSize,
                                   CvRect* validPixROI, int centerPrincipalPoint )
 {
-    cv::Rect_<float> inner, outer;
+    cv::Rect_<double> inner, outer;
     newImgSize = newImgSize.width*newImgSize.height != 0 ? newImgSize : imgSize;
 
     double M[3][3];
@@ -2545,10 +2545,10 @@ void cvGetOptimalNewCameraMatrix( const CvMat* cameraMatrix, const CvMat* distCo
 
         if( validPixROI )
         {
-            inner = cv::Rect_<float>((float)((inner.x - cx0)*s + cx),
+            inner = cv::Rect_<double>((double)((inner.x - cx0)*s + cx),
-                                     (float)((inner.y - cy0)*s + cy),
+                                      (double)((inner.y - cy0)*s + cy),
-                                     (float)(inner.width*s),
+                                      (double)(inner.width*s),
-                                     (float)(inner.height*s));
+                                      (double)(inner.height*s));
             cv::Rect r(cvCeil(inner.x), cvCeil(inner.y), cvFloor(inner.width), cvFloor(inner.height));
             r &= cv::Rect(0, 0, newImgSize.width, newImgSize.height);
             *validPixROI = cvRect(r);

modules/calib3d/test/test_undistort.cpp

@@ -157,6 +157,104 @@ void CV_DefaultNewCameraMatrixTest::prepare_to_validation( int /*test_case_idx*/
 
 //---------
 
+class CV_GetOptimalNewCameraMatrixNoDistortionTest : public cvtest::ArrayTest
+{
+public:
+    CV_GetOptimalNewCameraMatrixNoDistortionTest();
+protected:
+    int prepare_test_case (int test_case_idx);
+    void prepare_to_validation(int test_case_idx);
+    void get_test_array_types_and_sizes(int test_case_idx, vector<vector<Size> >& sizes, vector<vector<int> >& types);
+    void run_func();
+
+private:
+    cv::Mat camera_mat;
+    cv::Mat distortion_coeffs;
+    cv::Mat new_camera_mat;
+
+    cv::Size img_size;
+    double alpha;
+    bool center_principal_point;
+
+    int matrix_type;
+
+    static const int MAX_X = 2048;
+    static const int MAX_Y = 2048;
+};
+
+CV_GetOptimalNewCameraMatrixNoDistortionTest::CV_GetOptimalNewCameraMatrixNoDistortionTest()
+{
+    test_array[INPUT].push_back(NULL); // camera_mat
+    test_array[INPUT].push_back(NULL); // distortion_coeffs
+    test_array[OUTPUT].push_back(NULL); // new_camera_mat
+    test_array[REF_OUTPUT].push_back(NULL);
+
+    alpha = 0.0;
+    center_principal_point = false;
+    matrix_type = 0;
+}
+
+void CV_GetOptimalNewCameraMatrixNoDistortionTest::get_test_array_types_and_sizes(int test_case_idx, vector<vector<Size> >& sizes, vector<vector<int> >& types)
+{
+    cvtest::ArrayTest::get_test_array_types_and_sizes(test_case_idx, sizes, types);
+    RNG& rng = ts->get_rng();
+    matrix_type = types[INPUT][0] = types[INPUT][1] = types[OUTPUT][0] = types[REF_OUTPUT][0] = cvtest::randInt(rng)%2 ? CV_64F : CV_32F;
+    sizes[INPUT][0] = sizes[OUTPUT][0] = sizes[REF_OUTPUT][0] = cvSize(3,3);
+    sizes[INPUT][1] = cvSize(1,4);
+}
+
+int CV_GetOptimalNewCameraMatrixNoDistortionTest::prepare_test_case(int test_case_idx)
+{
+    int code = cvtest::ArrayTest::prepare_test_case( test_case_idx );
+
+    if (code <= 0)
+        return code;
+
+    RNG& rng = ts->get_rng();
+
+    alpha = cvtest::randReal(rng);
+    center_principal_point = ((cvtest::randInt(rng) % 2)!=0);
+
+    // Generate random camera matrix. Use floating point precision for source to avoid precision loss
+    img_size.width = cvtest::randInt(rng) % MAX_X + 1;
+    img_size.height = cvtest::randInt(rng) % MAX_Y + 1;
+    const float aspect_ratio = static_cast<float>(img_size.width) / img_size.height;
+    float cam_array[9] = {0,0,0,0,0,0,0,0,1};
+    cam_array[2] = static_cast<float>((img_size.width - 1)*0.5);  // center
+    cam_array[5] = static_cast<float>((img_size.height - 1)*0.5); // center
+    cam_array[0] = static_cast<float>(MAX(img_size.width, img_size.height)/(0.9 - cvtest::randReal(rng)*0.6));
+    cam_array[4] = aspect_ratio*cam_array[0];
+
+    Mat& input_camera_mat = test_mat[INPUT][0];
+    cvtest::convert(Mat(3, 3, CV_32F, cam_array), input_camera_mat, input_camera_mat.type());
+    camera_mat = input_camera_mat;
+
+    // Generate zero distortion matrix
+    const Mat zero_dist_coeffs = Mat::zeros(1, 4, CV_32F);
+    Mat& input_dist_coeffs = test_mat[INPUT][1];
+    cvtest::convert(zero_dist_coeffs, input_dist_coeffs, input_dist_coeffs.type());
+    distortion_coeffs = input_dist_coeffs;
+
+    return code;
+}
+
+void CV_GetOptimalNewCameraMatrixNoDistortionTest::run_func()
+{
+    new_camera_mat = cv::getOptimalNewCameraMatrix(camera_mat, distortion_coeffs, img_size, alpha, img_size, NULL, center_principal_point);
+}
+
+void CV_GetOptimalNewCameraMatrixNoDistortionTest::prepare_to_validation(int /*test_case_idx*/)
+{
+    const Mat& src = test_mat[INPUT][0];
+    Mat& dst = test_mat[REF_OUTPUT][0];
+    cvtest::copy(src, dst);
+
+    Mat& output = test_mat[OUTPUT][0];
+    cvtest::convert(new_camera_mat, output, output.type());
+}
+
+//---------
+
 class CV_UndistortPointsTest : public cvtest::ArrayTest
 {
 public:
@@ -935,6 +1033,7 @@ double CV_InitUndistortRectifyMapTest::get_success_error_level( int /*test_case_
 //////////////////////////////////////////////////////////////////////////////////////////////////////
 
 TEST(Calib3d_DefaultNewCameraMatrix, accuracy) { CV_DefaultNewCameraMatrixTest test; test.safe_run(); }
+TEST(Calib3d_GetOptimalNewCameraMatrixNoDistortion, accuracy) { CV_GetOptimalNewCameraMatrixNoDistortionTest test; test.safe_run(); }
 TEST(Calib3d_UndistortPoints, accuracy) { CV_UndistortPointsTest test; test.safe_run(); }
 TEST(Calib3d_InitUndistortRectifyMap, accuracy) { CV_InitUndistortRectifyMapTest test; test.safe_run(); }