cv2.triangulatePoints method generates different decimal values of the reconstructed points.

[mxochicale]

Triangulating points using cv2.triangulatePoints() over _iter_triangulate_point_w_svd() is definitely quicker (0.91077 milliseconds vs 2.426845 milliseconds), however cv2.triangulatePoints seems to use a different computational methods to resolve such points. See below decimal values of reconstructed points.

• CODE SNIPPETS

FROM: def triangulate_points_hartley(...

# array shapes for input args _iter_triangulate_point_w_svd( [3, 4]; [3, 4]; [3, 1]; [3, 1] )
reconstructed_point = _iter_triangulate_point_w_svd(p1d, p2d, u1t, u2t)
#reconstructed_point = np.around(reconstructed_point, 3) #Truncate and discarding decimal digits (np.around)
print(f'\n {reconstructed_point}')

FROM: def triangulate_points_opencv(...
# array shapes for input args cv2.triangulatePoints( [3, 4]; [3, 4]; [2, 1]; [2, 1] )
reconstructed_point = cv2.triangulatePoints(p1mat, p2mat, u1t[:2], u2t[:2])
reconstructed_point /= reconstructed_point[3]  # Homogenize
#reconstructed_point = np.around(reconstructed_point, 3) #Truncate and discarding decimal digits (np.around)
print(f'\n {reconstructed_point}')

• LOGS

$ pytest -v -s tests/algorithms/test_triangulate.py #for individual tests

tests/algorithms/test_triangulate.py::test_triangulate_points_hartley
 [[  9.88256179]
 [-22.44358672]
 [127.9216597 ]]

 [[ 48.46473971]
 [-23.09607652]
 [119.95244623]]

 [[  6.94680311]
 [  1.97332699]
 [115.79217294]]

 [[ 44.77458535]
 [  1.76836361]
 [106.8097294 ]]

 Elapsed time for at.triangulate_points_hartley(): 2.426845 millisecs

 [[  9.88261104]
 [-22.44364881]
 [127.92238956]
 [  1.        ]]

 [[ 48.46887878]
 [-23.09824868]
 [119.96328148]
 [  1.        ]]

 [[  6.94682563]
 [  1.97329098]
 [115.79279036]
 [  1.        ]]

 [[ 44.77479814]
 [  1.76837639]
 [106.81026738]
 [  1.        ]]

 Elapsed time for at.triangulate_points_opencv(): 0.91077 millisecs

 rms_hartley: 1.2954729076604021 and test (rms_hartley < 1.5)

 rms_hartley_opencv:  1.3009661570992705 and test (rms_hartley < 1.5)

I guess, we might like to dig into the implementation of cv2.triangulatePoints to understand where the differences comes from.

[thompson318]

I'm not sure it's worth spending too much time on this. For our purposes both methods pass the regression test, so I suggest just use the faster method.
The other issue is that historically the implementations in OpenCV aren't that stable from version to version. I wouldn't be surprised if a different version of OpenCV results in a different rms_hartley_opencv. We could end up spending a lot of time on it for little gain.

[mxochicale]

triangulate_points_opencv is faster when computing rms_hartley values, where both rms_hartley and rms_hartley_opencv pass rms_hartley < 1.5). However, there are failures for other unit tests: assert recon_err < 6.4; assert stereo_recon_err < 4.5 and assert (np.abs(s_recon - 1.64274596) < 0.000001) (see (a) for logs). My guess is that it might be related to the use of different methods cv::SVD::solve(), cv::SVD::solveZ and convertTo.

The OpenCV triangulation method, based on HartleyZ00 12.2 pag.312 and appendix of Keir's thesis: A.2 N-View Triangulation, pp. 106, has not been changed much since the fist/second commit on Oct 2015/Jan2018 (see (b)).

(a)reconstruction errors with openCV produces different errors in unit tests

• test_stereo_davinci

   def test_stereo_davinci():
    
        left_images = []
        files = glob.glob('tests/data/ChAruco_LR_frames_Steve_Axis_Tests/ExtractedFrames_L/*.jpg')
        files.sort()
        for file in files:
            image = cv2.imread(file)
            print("Loaded:" + str(file))
            left_images.append(image)
        assert(len(left_images) == 59)
    
        right_images = []
        files = glob.glob('tests/data/ChAruco_LR_frames_Steve_Axis_Tests/ExtractedFrames_R/*.jpg')
        files.sort()
        for file in files:
            image = cv2.imread(file)
            print("Loaded:" + str(file))
            right_images.append(image)
        assert (len(right_images) == 59)
    
        ref_img = cv2.imread('tests/data/2020_01_20_storz/pattern_4x4_19x26_5_4_with_inset_9x14.png')
    
        minimum_number_of_points_per_image = 50
        detector = pd.CharucoPlusChessboardPointDetector(ref_img,
                                                         error_if_no_chessboard=False) # Try to accept as many as possible.
        calibrator = sc.StereoVideoCalibrationDriver(detector, detector, minimum_number_of_points_per_image)
        for i, _ in enumerate(left_images):
            try:
                number_left, number_right = calibrator.grab_data(left_images[i], right_images[i])
                if number_left < minimum_number_of_points_per_image:
                    print("Image pair:" + str(i) + ", left image, SKIPPED, due to not enough points")
                if number_right < minimum_number_of_points_per_image:
                    print("Image pair:" + str(i) + ", right image, SKIPPED, due to not enough points")
            except ValueError as ve:
                print("Image pair:" + str(i) + ", FAILED, due to:" + str(ve))
            except TypeError as te:
                print("Image pair:" + str(i) + ", FAILED, due to:" + str(te))
    
        reproj_err, recon_err, params = calibrator.calibrate()
        print("Reproj:" + str(reproj_err))
        print("Recon:" + str(recon_err))
        assert reproj_err < 1.1
>       assert recon_err < 6.4
E       assert 6.512972286520402 < 6.4

tests/video/test_charuco_plus_chessboard.py:56: AssertionError

• test_precalbration

    def test_precalbration():
        """ Use intrinsics from A to calibration B, currently failing. """
        left_intrinsics = np.loadtxt('tests/data/precalib/precalib_base_data/calib.left.intrinsics.txt')
        left_distortion = np.loadtxt('tests/data/precalib/precalib_base_data/calib.left.distortion.txt')
        right_intrinsics = np.loadtxt('tests/data/precalib/precalib_base_data/calib.right.intrinsics.txt')
        right_distortion = np.loadtxt('tests/data/precalib/precalib_base_data/calib.right.distortion.txt')
        l2r = np.loadtxt('tests/data/precalib/precalib_base_data/calib.l2r.txt')
        l2r_rmat = l2r[0:3, 0:3]
        l2r_tvec = l2r[0:3, 3]
    
        calib_dir = 'tests/data/precalib/data_moved_scope'
        calib_driver = get_calib_driver(calib_dir)
    
        stereo_reproj_err, stereo_recon_err, _ = \
            calib_driver.calibrate(
                override_left_intrinsics=left_intrinsics,
                override_left_distortion=left_distortion,
                override_right_intrinsics=right_intrinsics,
                override_right_distortion=right_distortion,
                override_l2r_rmat=l2r_rmat,
                override_l2r_tvec=l2r_tvec)
    
        tracked_reproj_err, tracked_recon_err, _ = \
            calib_driver.handeye_calibration()
    
        print(stereo_reproj_err, stereo_recon_err, tracked_reproj_err, tracked_recon_err)
        assert stereo_reproj_err < 4.5
>       assert stereo_recon_err < 4.5
E       assert 17.573385358544893 < 4.5

tests/video/test_precalib.py:125: AssertionError

• test_stereo_video_calibration

    def test_stereo_video_calibration():
    
        object_points, left_image_points, right_image_points, ids = \
            load_first_stereo_data()
    
        s_reproj, s_recon, \
            l_c, l_d, left_rvecs, left_tvecs, \
            r_c, r_d, right_rvecs, right_tvecs, \
            l2r_r, l2r_t, \
            essential, fundamental = \
            vc.stereo_video_calibration(ids,
                                        object_points,
                                        left_image_points,
                                        ids,
                                        object_points,
                                        right_image_points,
                                        (1920, 1080))
    
        assert (np.abs(s_reproj - 0.63022577) < 0.000001)
>       assert (np.abs(s_recon - 1.64274596) < 0.000001)
E       AssertionError: assert 0.0015900760471676545 < 1e-06
E        +  where 0.0015900760471676545 = <ufunc 'absolute'>((1.6411558839528324 - 1.64274596))
E        +    where <ufunc 'absolute'> = np.abs

tests/video/test_video_calibration.py:82: AssertionError

(b) triangulation.cpp implementation in OpenCV

triangulation.cpp for 4.x has not been changed much since the fist/second commit on Oct 2015/Jan2018:

• triangulation.cpp first commit on Oct 5, 2015 based on Libmv and here
  https://github.com/opencv/opencv_contrib/blob/2a14d220c70ad8c0d33c16766a810a890d91c4ef/modules/sfm/src/triangulation.cpp
• Fixed several warnings produced by clang-6.0.0 on triangulation.cpp and others on Jan 30, 2018: https://github.com/opencv/opencv_contrib/blob/c99d1c3b0446037045acd0bf54bcbe148e721b5e/modules/sfm/src/triangulation.cpp