opencv/modules/calib3d/test/test_decompose_projection.cpp
Maxim Smolskiy e3ff6ce0cc
Merge pull request #25182 from MaximSmolskiy:increase-decomposeProjectionMatrix-precision-for-small-scales
Increase decomposeProjectionMatrix precision for small scales #25182 

### Pull Request Readiness Checklist

Fix #23733

It is checked before that `|s| > DBL_EPSILON` (if not, then `s` will be equal to `0`, but `c` will be equal to `1`, then `z` will be equal to `1` and there will be no any problems with small values), so `sqrt(c^2 + s^2) >= |s| > DBL_EPSILON` and thus small values are already taken into account before and there is no need to add `DBL_EPSILON` to `c^2 + s^2` (and I think adding `DBL_EPSILON^2` instead of `DBL_EPSILON` would be more correct).

I ran `Python` script from issue.

`NumPy` and `SciPy` results
```
*Numpy*
P (case 1): 
[[1. 0. 0. 0.]
 [0. 1. 0. 0.]
 [0. 0. 1. 0.]]
Numpy, R.T @ R - I: 
[[0. 0. 0.]
 [0. 0. 0.]
 [0. 0. 0.]]
Numpy scaled 1e-6, R.T @ R - I: 
[[0. 0. 0.]
 [0. 0. 0.]
 [0. 0. 0.]]
P (case 2): 
[[52. -7.  4. 12.]
 [-6. 49. 12.  8.]
 [ 4. 17.  1.  0.]]
Numpy, R.T @ R - I: 
[[-8.88178420e-16 -3.86302608e-16 -2.52050796e-17]
 [-3.86302608e-16 -5.55111512e-16  7.11675423e-18]
 [-2.52050796e-17  7.11675423e-18 -5.55111512e-16]]
Numpy scaled 1e-6, R.T @ R - I: 
[[ 2.22044605e-16 -2.00683644e-16 -1.90063998e-17]
 [-2.00683644e-16  0.00000000e+00  1.16926308e-17]
 [-1.90063998e-17  1.16926308e-17  2.22044605e-16]]

*Scipy*
P (case 1): 
[[1. 0. 0. 0.]
 [0. 1. 0. 0.]
 [0. 0. 1. 0.]]
Scipy, R.T @ R - I: 
[[0. 0. 0.]
 [0. 0. 0.]
 [0. 0. 0.]]
Scipy scaled 1e-6, R.T @ R - I: 
[[0. 0. 0.]
 [0. 0. 0.]
 [0. 0. 0.]]
P (case 2): 
[[52. -7.  4. 12.]
 [-6. 49. 12.  8.]
 [ 4. 17.  1.  0.]]
Scipy, R.T @ R - I: 
[[-1.11022302e-16 -8.74062812e-18 -1.26178867e-17]
 [-8.74062812e-18 -1.11022302e-16  2.07820373e-17]
 [-1.26178867e-17  2.07820373e-17 -1.11022302e-16]]
Scipy scaled 1e-6, R.T @ R - I: 
[[0.00000000e+00 4.04691435e-17 1.12452918e-16]
 [4.04691435e-17 4.44089210e-16 3.74164141e-16]
 [1.12452918e-16 3.74164141e-16 4.44089210e-16]]

*Numpy*
Numpy, P' - P:
 [[ 1.35525272e-20 -9.31736242e-21  8.47032947e-22  3.38813179e-21]
 [-3.38813179e-21  6.77626358e-21  1.69406589e-21  0.00000000e+00]
 [-1.69406589e-21  0.00000000e+00  0.00000000e+00  4.85524279e-22]]

*Scipy*
Scipy, P' - P:
 [[0.00000000e+00 8.47032947e-22 3.38813179e-21 3.38813179e-21]
 [3.38813179e-21 1.35525272e-20 1.52465931e-20 1.52465931e-20]
 [8.47032947e-22 3.38813179e-21 2.54109884e-21 3.39866995e-21]]
```

`OpenCV` results before
```
*OpenCV*
P (case 1): 
[[1. 0. 0. 0.]
 [0. 1. 0. 0.]
 [0. 0. 1. 0.]]
OpenCV, R.T @ R - I: 
[[0. 0. 0.]
 [0. 0. 0.]
 [0. 0. 0.]]
OpenCV scaled 1e-6, R.T @ R - I: 
[[0. 0. 0.]
 [0. 0. 0.]
 [0. 0. 0.]]
P (case 2): 
[[52. -7.  4. 12.]
 [-6. 49. 12.  8.]
 [ 4. 17.  1.  0.]]
OpenCV, R.T @ R - I: 
[[ 2.22044605e-16 -9.12253504e-17 -2.20527203e-19]
 [-9.12253504e-17  0.00000000e+00 -9.12405093e-18]
 [-2.20527203e-19 -9.12405093e-18  2.22044605e-16]]
OpenCV scaled 1e-6, R.T @ R - I: 
[[-1.28197013e-06  1.30450769e-07  7.67357467e-09]
 [ 1.30450769e-07 -1.52141637e-06 -9.92455574e-09]
 [ 7.67357467e-09 -9.92455574e-09 -1.35328272e-06]]

*OpenCV*
OpenCV, P' - P:
 [[-6.75449076e-11  1.73936564e-11 -4.94463312e-12 -1.61106020e-11]
 [ 1.41759913e-11 -7.54512016e-11 -1.67717374e-11 -9.74644390e-12]
 [-2.90254385e-12 -2.53521998e-11 -1.49130587e-12  4.00724440e-13]]
```

`OpenCV` results after
```
*OpenCV*
P (case 1): 
[[1. 0. 0. 0.]
 [0. 1. 0. 0.]
 [0. 0. 1. 0.]]
OpenCV, R.T @ R - I: 
[[0. 0. 0.]
 [0. 0. 0.]
 [0. 0. 0.]]
OpenCV scaled 1e-6, R.T @ R - I: 
[[0. 0. 0.]
 [0. 0. 0.]
 [0. 0. 0.]]
P (case 2): 
[[52. -7.  4. 12.]
 [-6. 49. 12.  8.]
 [ 4. 17.  1.  0.]]
OpenCV, R.T @ R - I: 
[[ 2.22044605e-16 -9.12253504e-17 -2.20527203e-19]
 [-9.12253504e-17  0.00000000e+00 -9.12405093e-18]
 [-2.20527203e-19 -9.12405093e-18  2.22044605e-16]]
OpenCV scaled 1e-6, R.T @ R - I: 
[[ 0.00000000e+00  4.36198333e-17 -2.66855078e-17]
 [ 4.36198333e-17  2.22044605e-16  3.17216400e-17]
 [-2.66855078e-17  3.17216400e-17 -2.22044605e-16]]

*OpenCV*
OpenCV, P' - P:
 [[ 6.77626358e-21  0.00000000e+00 -8.47032947e-22  1.69406589e-21]
 [-1.69406589e-21  6.77626358e-21  0.00000000e+00  3.38813179e-21]
 [ 8.47032947e-22  3.38813179e-21  2.11758237e-22  1.45657284e-21]]
```

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
- [ ] There is accuracy test, performance test and test data in opencv_extra repository, if applicable
      Patch to opencv_extra has the same branch name.
- [ ] The feature is well documented and sample code can be built with the project CMake
2024-03-27 12:43:56 +03:00

184 lines
5.5 KiB
C++

/*M///////////////////////////////////////////////////////////////////////////////////////
//
// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
// By downloading, copying, installing or using the software you agree to this license.
// If you do not agree to this license, do not download, install,
// copy or use the software.
//
//
// License Agreement
// For Open Source Computer Vision Library
//
// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
// Third party copyrights are property of their respective owners.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
// * Redistribution's of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
//
// * Redistribution's in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
//
// * The name of the copyright holders may not be used to endorse or promote products
// derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the Intel Corporation or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
// the use of this software, even if advised of the possibility of such damage.
//
//M*/
#include "test_precomp.hpp"
namespace opencv_test { namespace {
class CV_DecomposeProjectionMatrixTest : public cvtest::BaseTest
{
public:
CV_DecomposeProjectionMatrixTest();
protected:
void run(int);
};
CV_DecomposeProjectionMatrixTest::CV_DecomposeProjectionMatrixTest()
{
test_case_count = 30;
}
void CV_DecomposeProjectionMatrixTest::run(int start_from)
{
ts->set_failed_test_info(cvtest::TS::OK);
cv::RNG& rng = ts->get_rng();
int progress = 0;
for (int iter = start_from; iter < test_case_count; ++iter)
{
ts->update_context(this, iter, true);
progress = update_progress(progress, iter, test_case_count, 0);
// Create the original (and random) camera matrix, rotation, and translation
cv::Vec2d f, c;
rng.fill(f, cv::RNG::UNIFORM, 300, 1000);
rng.fill(c, cv::RNG::UNIFORM, 150, 600);
double alpha = 0.01*rng.gaussian(1);
cv::Matx33d origK(f(0), alpha*f(0), c(0),
0, f(1), c(1),
0, 0, 1);
cv::Vec3d rVec;
rng.fill(rVec, cv::RNG::UNIFORM, -CV_PI, CV_PI);
cv::Matx33d origR;
cv::Rodrigues(rVec, origR); // TODO cvtest
cv::Vec3d origT;
rng.fill(origT, cv::RNG::NORMAL, 0, 1);
// Compose the projection matrix
cv::Matx34d P(3,4);
hconcat(origK*origR, origK*origT, P);
// Decompose
cv::Matx33d K, R;
cv::Vec4d homogCameraCenter;
decomposeProjectionMatrix(P, K, R, homogCameraCenter);
// Recover translation from the camera center
cv::Vec3d cameraCenter(homogCameraCenter(0), homogCameraCenter(1), homogCameraCenter(2));
cameraCenter /= homogCameraCenter(3);
cv::Vec3d t = -R*cameraCenter;
const double thresh = 1e-6;
if (cv::norm(origK, K, cv::NORM_INF) > thresh)
{
ts->set_failed_test_info(cvtest::TS::FAIL_BAD_ACCURACY);
break;
}
if (cv::norm(origR, R, cv::NORM_INF) > thresh)
{
ts->set_failed_test_info(cvtest::TS::FAIL_BAD_ACCURACY);
break;
}
if (cv::norm(origT, t, cv::NORM_INF) > thresh)
{
ts->set_failed_test_info(cvtest::TS::FAIL_BAD_ACCURACY);
break;
}
}
}
TEST(Calib3d_DecomposeProjectionMatrix, accuracy)
{
CV_DecomposeProjectionMatrixTest test;
test.safe_run();
}
TEST(Calib3d_DecomposeProjectionMatrix, degenerate_cases)
{
for (int i = 0; i < 3; i++)
{
for (int j = 0; j < 2; j++)
{
cv::Matx34d P;
P(0, i) = 1;
P(1, (i + j + 1) % 3) = 1;
P(2, (i + 2 * j + 2) % 3) = 1;
cv::Matx33d K, R;
cv::Vec4d t;
decomposeProjectionMatrix(P, K, R, t);
EXPECT_LT(cv::norm(K * R, P.get_minor<3, 3>(0, 0), cv::NORM_INF), 1e-6);
}
}
}
TEST(Calib3d_DecomposeProjectionMatrix, bug_23733)
{
cv::Matx34d P(52, -7, 4, 12,
-6, 49, 12, 8,
4, 17, 1, 0);
P *= 1e-6;
cv::Matx33d K, R;
cv::Vec4d t;
decomposeProjectionMatrix(P, K, R, t);
EXPECT_LT(cv::norm(R.t() * R - cv::Matx33d::eye(), cv::NORM_INF), 1e-10);
cv::Matx34d M;
cv::hconcat(R, -R * cv::Vec3d(t[0] / t[3], t[1] / t[3], t[2] / t[3]), M);
cv::Matx34d P_recompose = K * M;
EXPECT_LT(cv::norm(P_recompose - P, cv::NORM_INF), 1e-16);
}
}} // namespace