SciKit-Surgery/scikit-surgerycalibration

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

Opened this issue · 2 comments

mxochicale commented

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 commented

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 commented

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: