Obtain Projection matrix from Intrinsic and Homograpy matrix

[edwardnguyen1705]

Dear @JarnoRalli ,

Thanks for sharing your great work. I would like to ask for help, so thanks so much for your time.

Given camera Intrinsic (K matrix), I try to obtain Extrinsic (Rotation matrix and translation vector) to form a camera projection matrix. I have tried your lib, but I am not successful yet.

Here is the procedure

1. 4 corresponding points pairs

• uv
• points3D_obj

2. Find homo
3. Obtain R and t

If I use forwardprojectK to obtain uv_, this works, the obtained uv_ is similar to the predefined uv.

However, if I use forwardprojectP to obtain uv_, this does not work the obtained uv_ is very different from the predefined `uv'.

And I do not know why the later method does not work, the obtained points3D_cam of the two methods are very different.

points3D_cam = R @ points3D_obj + t

R = np.transpose(R)
world2cam = np.hstack((R, np.dot(-R, t).reshape(3,-1)))
points3D_h = np.vstack([points3D_h, np.ones([1, 4], dtype=np.float32)])
points3D_cam = world2cam@points3D_h

Here is the full src for reproducing.

import numpy as np
from python_camera_library import rectilinear_camera as cam
from python_camera_library import homography as homog
from python_camera_library import utils


K = np.array([[9.32545185e+02, 0.00000000e+00, 1.10520066e+03], 
     [0.00000000e+00, 9.32545185e+02, 4.84232068e+02], 
     [0.00000000e+00, 0.00000000e+00, 1.00000000e+00]], dtype=np.float32)

uv = np.array([[432, 1167, 550, 1524],
               [593, 495, 873, 571],
               [1, 1, 1, 1]], dtype=np.float32)

# mm
points3D_obj = np.array([[0, 15240,    0, 15240],
                        [0,      0, 8530,  8530],
                        [0,      0,    0,     0]], dtype=np.float32)

points3D_h = np.copy(points3D_obj)
points3D_h[-1, :] = np.ones([1, 4])

H = homog.dlt_homography(points3D_h, uv)
(R, t) = homog.extract_transformation(K, H)
scale = homog.transformation_scale(K, R, t, uv[:, 1], points3D_obj[:, 1])
t *= scale

points3D_cam = R @ points3D_obj + t
(pts3D, uv_, depth) = cam.forwardprojectK(points3D_cam, K, (6000, 6000))
print("points3D_cam = R @ points3D_obj + t")
print(points3D_cam)
print()
print("(pts3D, uv_, depth) = cam.forwardprojectK(points3D_cam, K, (6000, 6000))")
print(uv_)

R = np.transpose(R)
world2cam = np.hstack((R, np.dot(-R, t).reshape(3,-1)))
points3D_h = np.vstack([points3D_h, np.ones([1, 4], dtype=np.float32)])
points3D_cam = world2cam@points3D_h
print()
print("points3D_cam = world2cam@points3D_h")
print(points3D_cam)

P = np.dot(K, world2cam)
_, uv, _ = cam.forwardprojectP(points3D_obj, P, (6000, 6000))
print()
print("_, uv, _ = cam.forwardprojectP(points3D_obj, P, (6000, 6000))")
print(uv_)

[edwardnguyen1705]

Dear @JarnoRalli ,

I know how to fix the issue, just stacking R and t. Now, I try to project 2D points to 3D points, on the z plane.

world2cam = np.hstack((R, t))
world2cam = np.vstack([world2cam, np.array([0, 0, 0, 1], dtype=np.float32)])

points3D_h = np.vstack([points3D_h, np.ones([1, 4], dtype=np.float32)])
points3D_cam = world2cam@points3D_h
print()
print("points3D_cam = world2cam@points3D_h")
print(points3D_cam)

K_h = np.hstack([K, np.zeros((3, 1), dtype=np.float32)])
P = np.dot(K_h, world2cam)
print("Project matrix")
print(P)

_, uv, _ = cam.forwardprojectP(points3D_obj, P, (6000, 6000))
print()
print("_, uv, _ = cam.forwardprojectP(points3D_obj, P, (6000, 6000))")
print(uv_)

[JarnoRalli]

Hi @edwardnguyen1705 ,

I'm glad that you find the library useful. In the utils.py there is actually a function that is used for concatenating R (rotation) and t (translation) together:

def concatenateRt(R: np.array, t: np.ndarray) -> np.ndarray:
    """Concatenates a rotation matrix R and a translation vector t into a transformation matrix T.

    The outcome is as follows:
    [ R      | t]
    [[0 0 0] | 1]

    Parameters
    ----------
    R : np.array
        Rotation matrix, 3x3.
    t : np.ndarray
        Translation vector, 3x1.

    Returns
    -------
    np.ndarray
        Transformation matrix, 3x4.
    """

Unfortunately I'm quite busy the next week, but if you need more clarifications, please let me know.

[edwardnguyen1705]

Hi @JarnoRalli ,

Thanks for your reply.

I have obtained a projection matrix and I check it and it works. However, when I try a program with input is the projection matrix it does not work.
The doc of the program says that:
...the principal point (i.e., (Cx, Cy)) in the camera matrix is assumed to be at (0, 0) as image coordinates. But, the optical center (i.e., (Cx, Cy)) is located at the image center (i.e., (img_width/2, img_height/2)). Thus, to move the origin to the left-top of the camera image (i.e., the pixel coordinates), SV3DT internally adds (img_width/2, img_height/2) after the transformation using the camera matrix provided in projectionMatrix_3x4.
The 3x4 Camera Projection Matrix

I am confusing
"the principal point (i.e., (Cx, Cy)) in the camera matrix": is that in the intrinsic matrix the last column is [0, 0, 1]'? However, OpenCV cam calib always returns estimated (Cx, Cy)?

I see one of your examples:

python_camera_library/homography_example.py

Line 30 in b6154a1

K = np.array([[100, 0, 150], [0, 100, 150], [0, 0, 1]])

you also set the intrinsic matrix the last column be [0, 0, 1]'

If I set it for my intrinsic matrix, and I check that the projection matrix does not work.

If you know about my issue, please help me.
Thanks so much for your time.

[JarnoRalli]

Format of the camera calibration matrix K is as follows (assuming that the pixels are not slanted):

      [fx  0  cx]
K = [0  fy  cy]
      [0  0   1]

, where fx and fy are the focal lengths (typically fx=fy) and cx and cy is the principal point (typically close the centre of the image).
If you're trying to obtain the R and t from point correspondences, then you can use PnP from OpenCV: https://docs.opencv.org/4.x/d5/d1f/calib3d_solvePnP.html, those will give robust estimates.

[edwardnguyen1705]

Format of the camera calibration matrix K is as follows (assuming that the pixels are not slanted):

      [fx  0  cx]
K = [0  fy  cy]
      [0  0   1]

, where fx and fy are the focal lengths (typically fx=fy) and cx and cy is the principal point (typically close the centre of the image). If you're trying to obtain the R and t from point correspondences, then you can use PnP from OpenCV: https://docs.opencv.org/4.x/d5/d1f/calib3d_solvePnP.html, those will give robust estimates.

Hello @JarnoRalli ,

As you suggest, I use PnP and I obtain a more accurate result.
Thank you!