mxochicale commented 1 year ago

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 `**](https://github.com/opencv/opencv_contrib/blob/4.x/modules/sfm/src/triangulation.cpp) to understand where the differences comes from.

thompson318 commented 1 year ago

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 1 year ago

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 ](https://github.com/opencv/opencv_contrib/blob/4.x/modules/sfm/src/triangulation.cpp)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](https://developer.blender.org/project/profile/59) and [here](https://developer.blender.org/diffusion/LMV/repository/master/)
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

SciKit-Surgery / scikit-surgerycalibration

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

array shapes for input args _iter_triangulate_point_w_svd( [3, 4]; [3, 4]; [3, 1]; [3, 1] )

reconstructed_point = np.around(reconstructed_point, 3) #Truncate and discarding decimal digits (np.around)

array shapes for input args cv2.triangulatePoints( [3, 4]; [3, 4]; [2, 1]; [2, 1] )

reconstructed_point = np.around(reconstructed_point, 3) #Truncate and discarding decimal digits (np.around)

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

SciKit-Surgery / scikit-surgerycalibration

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

array shapes for input args _iter_triangulate_point_w_svd( [3, 4]; [3, 4]; [3, 1]; [3, 1] )

reconstructed_point = np.around(reconstructed_point, 3) #Truncate and discarding decimal digits (np.around)

array shapes for input args cv2.triangulatePoints( [3, 4]; [3, 4]; [2, 1]; [2, 1] )

reconstructed_point = np.around(reconstructed_point, 3) #Truncate and discarding decimal digits (np.around)

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

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