Acquire each camera's extrinsic martix

[yichuan1998]

Hi, Mac, it's me. Recently I was learning about multi-view geometry for human pose estimation. I managed to calibrate 2 cameras using cv2 but did't know how to calibrate 3 cameras. I could not find tutorial fromgoogle or github. I tried to learn it from the source code of pyxy3d and failed. Could you give me some advice? Thanks!

The code for calibrating two cameras:

def read_charucoboards(images):
    """
    Charuco base pose estimation.
    """
    allCorners = []
    allIds = []
    decimator = 0
    # SUB PIXEL CORNER DETECTION CRITERION
    criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 100, 0.00001)

    for im in images:
        frame = cv2.imread(im)
        gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)
        corners, ids, rejectedImgPoints = cv2.aruco.detectMarkers(gray, aruco_dict)

        if len(corners)>0:
            # SUB PIXEL DETECTION
            for corner in corners:
                cv2.cornerSubPix(gray, corner,
                                 winSize = (3,3),
                                 zeroZone = (-1,-1),
                                 criteria = criteria)
            res2 = cv2.aruco.interpolateCornersCharuco(corners,ids,gray,board)
            if res2[1] is not None and res2[2] is not None and len(res2[1])>3 and decimator%1==0:
                allCorners.append(res2[1])
                allIds.append(res2[2])

        decimator += 1

    imsize = gray.shape
    return allCorners,allIds,imsize


def calibrate_camera(allCorners,allIds,imsize):
    """
    Calibrates the camera using the dected corners.
    """

    cameraMatrixInit = np.array([[ 1000.,    0., imsize[0]/2.],
                                 [    0., 1000., imsize[1]/2.],
                                 [    0.,    0.,           1.]])

    distCoeffsInit = np.zeros((5,1))
    flags = (cv2.CALIB_USE_INTRINSIC_GUESS + cv2.CALIB_RATIONAL_MODEL + cv2.CALIB_FIX_ASPECT_RATIO)
    #flags = (cv2.CALIB_RATIONAL_MODEL)
    (ret, camera_matrix, distortion_coefficients0,
     rotation_vectors, translation_vectors,
     stdDeviationsIntrinsics, stdDeviationsExtrinsics,
     perViewErrors) = cv2.aruco.calibrateCameraCharucoExtended(
                      charucoCorners=allCorners,
                      charucoIds=allIds,
                      board=board,
                      imageSize=imsize,
                      cameraMatrix=cameraMatrixInit,
                      distCoeffs=distCoeffsInit,
                      flags=flags,
                      criteria=(cv2.TERM_CRITERIA_EPS & cv2.TERM_CRITERIA_COUNT, 10000, 1e-9))
    
    return camera_matrix, distortion_coefficients0


def calibrate_stereo(matrix0, dist0, matrix1, dist1):
    criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 100, 0.00001)
    imgs_cam0 = readimgs(readimgspath(config["cam0"]))
    imgs_cam1 = readimgs(readimgspath(config["cam1"]))
    

    objpoints = [] # 3d point in real world space
    imgpoints0 = [] # 2d points in image plane.
    imgpoints1 = [] # 2d points in image plane.

    for frame0, frame1 in zip(imgs_cam0, imgs_cam1):
        corners0, ids0, rejectedImgPoints0 = cv2.aruco.detectMarkers(frame0, aruco_dict)
        corners1, ids1, rejectedImgPoints1 = cv2.aruco.detectMarkers(frame1, aruco_dict)
        if not corners0 or not corners1: raise Exception("No markers detected")

        ret0, corners0, ids0 = cv2.aruco.interpolateCornersCharuco(corners0, ids0, frame0, board)
        ret1, corners1, ids1 = cv2.aruco.interpolateCornersCharuco(corners1, ids1, frame1, board)
        if not ret0 or not ret1: raise Exception("Can't interpolate corners")

        objpoints.append(board.getChessboardCorners())
        imgpoints0.append(corners0.reshape(-1, 2))
        imgpoints1.append(corners1.reshape(-1, 2))

    ret, CM0, dist0, CM1, dist1, R, T, E, F = cv2.stereoCalibrate(objpoints, imgpoints0, imgpoints1, matrix0, dist0,
                                                                 matrix1, dist1, frame_size, criteria = criteria, flags = cv2.CALIB_FIX_INTRINSIC)


    return R, T


allCorners0, allIds0, imsize0 = read_charucoboards(readimgspath(r"D:\my_work\Only_test\charuco_stereo_cali\cam0"))
matrix0, dist0 = calibrate_camera(allCorners0, allIds0, imsize0)
calibrate_visualize(r"D:\my_work\Only_test\charuco_stereo_cali\cam0", matrix0, dist0, "cam0")

allCorners1, allIds1, imsize1 = read_charucoboards(readimgspath(r"D:\my_work\Only_test\charuco_stereo_cali\cam1"))
matrix1, dist1 = calibrate_camera(allCorners1, allIds1, imsize1)
calibrate_visualize(r"D:\my_work\Only_test\charuco_stereo_cali\cam1", matrix1, dist1, "cam1")

R, T = calibrate_stereo(matrix0, dist0, matrix1, dist1)

Sincerely

Murphy

[mprib]

Hi @yichuan1998 (Murphy)!

The phrase you want to google is "bundle adjustment." This is a process for estimating the point locations and the camera positions in such a way that it minimizes the reprojection error (or some cost function based on the reprojection error).

It seems intimidating, but really it's just a least_squares problem. This was the example I started with to get pyxy3d working:
https://scipy-cookbook.readthedocs.io/items/bundle_adjustment.html

There is also a lecturer on youtube whose videos really helped me get my head around it:
https://www.youtube.com/results?search_query=cyrill+stachniss+bundle+adjustment

In pyxy3d the bundle adjustment happens within the CaptureVolume class:
https://github.com/mprib/pyxy3d/blob/ec67a2f94134134174ad17aebea398b2bbe0b5e8/pyxy3d/calibration/capture_volume/capture_volume.py#L93

The capture volume is composed of the camera position estimates and the 3d point estimates of the charuco corners taken during the multicamera calibration. Basically, you start with your initial estimates of where the cameras are and where the board is, then figure out where the point would be projected on the camera view if your estimate were 100% accurate (including estimates of lens distortion). This will differ from where the point actually shows up in the image. The difference between the "projected point" and the actual point is the "reprojection error". With the reprojection error defined as a function of the 3d point positions and the camera positional estimates, the optimization basically runs through some kind of gradient descent to get the best fit. It nudges the points and the cameras around and rechecks if the error improves.

The thing that pyxy3d does that I'm proud of (and that results in a very fast, reliable convergence to an optimum) is that it is very careful with the initialization of the estimates. It basically follows these stages:

1. Stereocalibrate all possible pairs of cameras
2. rank camera calibration quality based on RMSE
3. Choose an anchor camera that has the best stereocalibrations
4. Construct a complete estimate of all camera positions relative to the anchor. If necessary, daisy-chain cameras to the anchor via a camera-in-common. This set of estimates is the "camera_array" part of the CaptureVolume
5. With all camera positions estitimated, I think it just does a simple stereotriangulation from all pairs and then averages to get the "point_estimates"
6. camera_array and point_estimates create the CaptureVolume, which can then be optimized via bundle adjustment.

Unfortunately (and surprisingly) there is not a stock cv2 approach to this. I think you kinda have to roll your own...

---

I realize that I'm throwing a lot at you at once...One thing I will say is that if you are going down this path to roll your own, prioritize a method for visualizing your current estimates because you want a way to quickly see when a small error has started to creep into your calculations (which often has a big impact and will yield weird results). Please feel free to keep questions coming!

[yichuan1998]

It's very kind of you for telling me these useful tips. According to my understanding, triangulation is needed before bundle adjustmen. I'm also confused about triangulation for 3 cameras(I only made it on two cameras). How can I obtain each camera's rotation matrix and translation vector？

Stereo calibration outputs the rotation matrix and translation vector from camera0 to camera1. In what way do you obtain rotation matrix and translation vector for each camera?

And why the size of rotation is 3x1, not 3x3?

Is it Rodrigues rotation？

I'm really a novice at camera calibration and 3d reconstruction. Thanks very much!

[mprib]

You are correct regarding the Rodrigues. The rotation within each camera is the 3x1 version and the rotation stored within the stereo_pair entries is the 3x3 equivalent. This is just an arbitrary artifact of the way things are stored. cv2.rodrigues will convert between the two different forms of the rotation and during the matrix math the 3x3 version is used.

To link together more than two cameras in a common frame of reference, you can use a single camera as the anchor and stereocalibrate all other cameras relative to that. This is self-limiting because you can't always get a good view between an anchor and all other cameras.

So what to do?

You can then bridge stereopairs to put things in a common frame of reference. That happens in the code here:

https://github.com/mprib/pyxy3d/blob/ec67a2f94134134174ad17aebea398b2bbe0b5e8/pyxy3d/cameras/camera_array_initializer.py#L71

It involves using the transformation matrix which is a merge of the rotation and translation matrices, so that's what you'll want to start googling. In the simplest case you can invert the transformation matrix to swap the "anchor" camera in a stereopair (at least if memory serves).

I will say that you are walking down a topic where I found precious little relevant information online that I could interpret. Searching stack overflow for information about this seems to lead to many people with deep expertise in CV being....less than helpful. I kinda had to just mess around in conjunction with a good visualizer to make sure that I was understanding things correctly. So I'll underscore again: prioritize having a method to visualize your calculations in 3D before you start diving deep into doing the calculations. Pyqtgraph is an option. Python scripting in Blender may also be useful here as they have many camera /camera position features built in. I'm sure there are other high level 3d tools that can help with this, so probably worth poking around.

[mprib]

And on the topic of triangulating with 3 cameras, it remains fundamentally a least-squares problem where you are trying to minimize the reprojection error across all cameras. I employ the code solution developed by Lili Karashchuk (lambdaloop on github) within Anipose:

https://github.com/mprib/pyxy3d/blob/f0c69821c3a744d16b5bc44bbba0f9763a9a279c/pyxy3d/triangulate/triangulation.py#L63

This uses singular value decomposition (np.linalg.svd) for the least squares solution.

---

And just an edit to mention that for initializing the point estimates for the bundle adjustment I think I just average the stereopair triangulations once I have all the cameras in a common frame of reference.

[yichuan1998]

Thanks again for your detailed reply. I am gradually understanding these complex concepts gradually. for visualization, I use plt to draw all pose estimated 3d points in 3d figure.

This is, obviously, a wrong estimation of human pose, largely due to getting each camera's extrinsic matrix wrong. I will keep trying to figure out it.

[mprib]

@yichuan1998

I just wanted to check back in to let you know that a significant refactor has largely been completed (v0.3) and the workflow now fits in with pre-recorded video. It does require synchronized video for the multicamera calibration and triangulation, so that synchronization would have to be done outside of pyxy3d. While I'm still tracking down little issues here and there, I'm hoping to start getting feedback from others on the functionality of the package, so thought I'd reach out in case you cared to take a look.

My brain isn't yet ready to try to write up the documentation, but here are some video walk throughs of the process in case you are interested:

• description of basic file/folder structure: https://youtu.be/dvosFrFVDsU
• demo of core workflow: https://www.youtube.com/watch?v=XDfuQ_7m1oQ

Mac

[yichuan1998]

It is glad to hear that Pyxy3d has been making progress. I'd would love to test the new version of the package and provide you feedbacks. But recently I'm preparing the entrance examination of post graduate. I would be happy to test this after Christmas. Thanks for informing me this.

Murphy

[mprib]

I wish you great success in your exam!

[yichuan1998]

Hey, Mac, how was your Christmas day? I tested the new version and got some problems as below.

1. I created a new project in new folder but the log showed：
   
   Is it normal？
2. Nothing's there when I clicked 'camera' butten:
   
   Does it means that users have to record as many videos as the number of cameras in advance and renamed like 'port_1.mp4,port_2.mp4' ?

[mprib]

@yichuan1998 ,

Thank you for taking the time to take a look at this and share your experiences. It is incredibly helpful to me and much appreciated. I have fixed the bug that you described.

The workflow is now modified to decouple the recording from the post-processing so that a wider variety of hardware approaches could be used. There is a sample dataset here that illustrates the file and folder layout here in case you have any interest in poking around some more: https://utexas-my.sharepoint.com/:f:/g/personal/priblede_my_utexas_edu/EliBf-eSuxJMpltXfANlS4oB35S9muKtpCxg4VMKX9rt2Q?e=Ecwi62

Your experience helps me to realize where I need to improve the documentation and the feedback provided through the Workspace tab of the GUI. I will continue working on it so that hopefully it will be more self-explanatory to users in the future. Thank you for taking a look and sharing your feedback!

Mac