Forward Projection From 3D To 2D And Triangulation
1- General Information
x=PX
• X∈R3: 3D point in Cartesian coordinates [X, Y, Z]
X∈R4: 3D point in homogeneous coordinates [X, Y, Z, 1]
• x∈R2: projected 2D image point [x, y]
x∈R3: projected 2D image point homogenous coordinate [x, y, 1]
• P3x4: Projection matrix
and:
P=K[R∣t]
K3x3: Intrinsic camera matrix
R3x3: Camera Orientation in the world
T3x1: Camera Position in the World
2- Calculating Intrinsic Camera Parameters
In my code, I used the checkerboard images to calculate the camera parameters using this
line of code:
ret, mtx, dist, rvecs, tvecs = [Link](objpoints, imgpoints, [Link][::-1],
None, None)
mtx: camera intrinsic parameters
dist: camera distortion coefficients
rvecs: camera rotation matrix
tvecs: camera position matrix
and we can use the following lines to construct the projection matrix:
rvec = rvecs[0]
tvec = tvecs[0]
R, _ = [Link](rvec)
Rt = [Link]((R, tvec))
P = K @ Rt # shape (3, 4)
However, in this code, since we did not have real-world object point coordinates
(checkerboard cell corners), the calculated rotation and translation matrices are not
accurate. We cannot use them for further 3D reconstruction, and only intrinsic camera
parameters will be used.
More importantly, since we moved the checkerboard in front of the camera, the camera
observed the checkerboard from a different viewpoint in each frame, and each image
provides information on its own pose of the checkerboard. Therefore, the camera’s position
and orientation (extrinsic parameters) relative to the checkerboard changed in each frame,
resulting in a distinct projection matrix for each image. As a result, the rotation (R) and
translation (T) matrices vary across images.
3- Estimating The Projection Matrix for Each Camera
As we were unable to directly calculate the projection matrix in the previous step, we need
a fitting model to estimate the mapping from 3D to 2D space. This method is known as the
Direct Linear Transformation (DLT) algorithm, which is a parameter estimation method
used to compute the forward projection matrix (P) for each camera based on known 2D
image points (checkerboard corner coordinates) observed by that camera.
For the purpose of this paper—and as requested by Dr. Dutta—we used the EasyWand
MATLAB tool, which uses a Sparse Bundle Adjustment (SBA) optimization method. But, in
my pose estimation project, I used its OpenCV equivalent: the solvePnP function, which
also performs optimization to estimate the camera pose from 2D-3D correspondences.
4- Estimating The Projection Matrix for Each Camera
After utilizing the DLT method to estimate the projection matrices P1 and P2 for the right
and left cameras, I developed the following triangulation function to reconstruct 3D points
from the 2D projections. However, the following command in OpenCV, cv2.
[Link], can also be used for this purpose.
def ReconstructButtonClicked(self):
file = open(DLTfile)
reader = [Link](file)
rows = []
for row in reader:
[Link](row)
LC1 = [float(rows[0][0]), float(rows[1][0]), float(rows[2][0]), float(rows[3][0]),
float(rows[4][0]),
float(rows[5][0]), float(rows[6][0]), float(rows[7][0]), float(rows[8][0]), float(rows[9][0]),
float(rows[10][0])]
LC2 = [float(rows[0][1]), float(rows[1][1]), float(rows[2][1]), float(rows[3][1]),
float(rows[4][1]),
float(rows[5][1]), float(rows[6][1]), float(rows[7][1]), float(rows[8][1]), float(rows[9][1]),
float(rows[10][1])]
RightCamPlanePoints = [[p1xR,p1yR],[p2xR,p2yR],[p3xR,p3yR]]
LeftCamPlanePoints = [[p1xL, p1yL], [p2xL, p2yL], [p3xL, p3yL]]
PlanePtsCameraView = []
for i in range(len(RightCamPlanePoints)):
u1 = RightCamPlanePoints[i][0]
v1 = RightCamPlanePoints[i][1]
u2 = LeftCamPlanePoints[i][0]
v2 = LeftCamPlanePoints[i][1]
u1C1 = u1
v1C1 = v1
u1C2 = u2
v1C2 = v2
Y = [u1C1 - LC1[3], v1C1 - LC1[7], u1C2 - LC2[3], v1C2 - LC2[7]]
A=[
[LC1[0] - u1C1 * LC1[8], LC1[1] - u1C1 * LC1[9], LC1[2] - u1C1 * LC1[10]],
[LC1[4] - v1C1 * LC1[8], LC1[5] - v1C1 * LC1[9], LC1[6] - v1C1 * LC1[10]],
[LC2[0] - u1C2 * LC2[8], LC2[1] - u1C2 * LC2[9], LC2[2] - u1C2 * LC2[10]],
[LC2[4] - v1C2 * LC2[8], LC2[5] - v1C2 * LC2[9], LC2[6] - v1C2 * LC2[10]]
]
A = [Link](A)
AT = [Link](A)
ATA = [Link](AT, A)
inv_ATA = [Link](ATA)
ATAinv_AT = [Link](inv_ATA, AT)
Y = [Link](Y)
X = [Link](ATAinv_AT, Y) # in Meters
[Link]()
[Link]([Link](X * 1000, 0)) # in mm
# print(PlanePtsCameraView)
if [Link](ProjectWorkingPath + '/files/[Link]'):
[Link](ProjectWorkingPath + '/files/[Link]', PlanePtsCameraView,
delimiter=',')
print("Pyramid heads plane is created!")
else:
[Link](ProjectWorkingPath + '/files/[Link]', PlanePtsCameraView,
delimiter=',')
print("Ground plane is created!")
[Link]("")
[Link]("")
[Link]("")
[Link]("")
[Link]("")
[Link]("")
return PlanePtsCameraView