Is my perspective math correct?

Question

I've got a homework in which I have to calculate and plot some points using a pespective transformation, but I'm not sure my results are correct, since the 3d plot using Camera coordinates looks very different from the 2d plot using the image coordinates. Can you help me understand what's wrong?

This is what is given: The camera is at the point $_WT^C = [−1, 1, 5]^T$, specified in world coordinates (in meters). The camera coordinate system is rotated around the Y axis of the world reference by $\theta = 160^o$, so it's rotation matrix is $^wR_c = \begin{bmatrix}cos(\theta) & 0 & sin(\theta)\\ 0 & 1 & 0 \\ -sin(\theta) & 0 & cos(\theta)\end{bmatrix}$

Camera parameter are: $f = 16mm$, $s_x = s_y = 0.01 mm/px$, $o_x = 320 px$, $o_y = 240px$

Sample points (in world coordinates):

$^WP_1 = [1, 1, 0.5]^T$

$^WP_2 = [1, 1.5, 0.5]^T$

$^WP_3 = [1.5, 1.5, 0.5]^T$

$^WP_4 = [1.5, 1, 0.5]^T$

I have to calculate and plot the points in camera coordinates and in image coordinates, so I wrote the following code in Octave:

%camera intrinsic parameters
f = 16
Sx = 0.01
Sy = 0.01
Ox = 320
Oy = 240

%given points, in world coordinate
wP1 = transpose([1, 1, 0.5])
wP2 = transpose([1, 1.5, 0.5])
wP3 = transpose([1.5, 1.5, 0.5])
wP4 = transpose([1.5, 1, 0.5])

% camera translation matrix
wTc = transpose([-1, 1, 5])

% rotation angle converted to rad
theta = 160 / 180 * pi

%camera rotation matrix
wRc = transpose([cos(theta), 0, sin(theta); 0, 1, 0; -sin(theta), 0, cos(theta)])

%transform the points to homogeneous coordinates
wP1h = [wP1; 1]
wP2h = [wP2; 1]
wP3h = [wP3; 1]
wP4h = [wP4; 1]

%separate each line of the rotation matrix
R1 = transpose(wRc(1 , :))
R2 = transpose(wRc(2 , :))
R3 = transpose(wRc(3 , :))

%generate the extrinsic parameters matrix
Mext = [wRc, [-transpose(R1) * wTc; -transpose(R2) * wTc; -transpose(R3) * wTc]]

%intrinsic parameters matrix
Mint = [-f/Sx, 0, Ox; 0, -f/Sy, Oy; 0, 0, 1]

% calculate coordinates in camera coordinates
cP1 = wRc * (wP1 - wTc)
cP2 = wRc * (wP2 - wTc)
cP3 = wRc * (wP3 - wTc)
cP4 = wRc * (wP4 - wTc)

% put coordinates in a list for plotting

x = [cP1(1), cP2(1), cP3(1), cP4(1), cP1(1)]
y = [cP1(2), cP2(2), cP3(2), cP4(2), cP1(2)]
z = [cP1(3), cP2(3), cP3(3), cP4(3), cP1(3)]

%plot the points in 3D using camera coordinates
plot3(x, y, z, "o-r")

pause()

% calculate the points in image coordinates
iP1 = Mint * (Mext * wP1h)
iP2 = Mint * (Mext * wP2h)
iP3 = Mint * (Mext * wP3h)
iP4 = Mint * (Mext * wP4h)

%generate a list of points for plotting
x = [iP1(1) / iP1(3), iP2(1) / iP2(3), iP3(1) / iP3(3), iP4(1) / iP4(3), iP1(1) / iP1(3)]
y = [iP1(2) / iP1(3), iP2(2) / iP2(3), iP3(2) / iP3(3), iP4(2) / iP4(3), iP1(2) / iP1(3)]

plot(x, y, "o-r")

pause()

And these are the plots I've got from the script: I was expecting they were somewhat similar, but they don't look so.

Plot in camera coordinates

Plot in image coordinates

Answer 1

Identifying your axes in both figures and adding the camera position to your first figure would help you understand what's going on.

You could also have a single variables for all your points, generating a 2D matrix with the rows as each point and columns as the components $x$, $y$ and $z$. That way, you could handle the projection using a simple matrix multiplication instead of handling each row separately.

In the problem statement, would be interesting to know the direction of your rotation, and more importantly, the original camera direction and its up vector. I guess you camera is rotated by 160° counter-clockwise, the original camera direction is $[0, 0, 1]$ and its up vector is $[0, 1, 0]$. If any of those assumptions is wrong, the rest of the answer will be wrong.

There is no sign of mm to m conversion in your code. Either your focal should be $0.016$ and $S_x = S_y = 0.0001$ or your point coordinates should be multiplied by $0.00001$.

Let's reason about where should the points end up in your image. For example, your camera center, without rotation, is pointing at the $[-1,1,x]$ line. Since the points are all on the $z=0.5$ plane, we can do the following analysis: if we focus on the $x$-axis, we can see that $tan(160°) \cdot (5 - 0.5) = 1.64...$, so the camera center will end up slightly to the left of the points (as the camera is at $x=-1$), so the center will end up at $\approx 0.64$, meaning that the points will appear on the right section of the image. Furthermore, the camera has the same $y$ coordinates as two point, and since the $y$ coordinates are not changed by the rotation, they should still end up on the same coordinates after the transform, meaning on the central row of the image.

A good way to check your answer is by using an existing 3D modeler like Blender: Be cautious with Blender' coordinates system, for example the default camera vector is [0, 0, -1]. Here's the render: Focal was set to another value to make the sphere more visible. So we see that the bottom two points are on the middle row of the image and the points are slightly to the right of the image.

I implemented your homework in Python:

import numpy as np

from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import axes3d, Axes3D


# Parameters
f_mm = 0.016
f_px = f_mm / 0.00001
t_cam = np.array([[-1., 1., 5.]]).T
t_cam_homogeneous = np.vstack((t_cam, np.array([[0]])))
theta = 160. * np.pi / 180.
ox = 320
oy = 240
# Rotation and points are in homogeneous coordinates
rot_cam = np.array([[np.cos(theta), 0, np.sin(theta)],
                    [0, 1, 0],
                    [-np.sin(theta), 0, np.cos(theta)]])
points = np.array([[1, 1, 0.5, 1],
                   [1, 1.5, 0.5, 1],
                   [1.5, 1.5, 0.5, 1],
                   [1.5, 1, 0.5, 1]]).T

# Compute projection matrix using intrinsics and extrinsics
intrinsics = np.array([[f_px, 0, ox],
                       [0, f_px, oy],
                       [0, 0, 1]])
extrinsics = np.hstack((rot_cam, rot_cam.dot(-t_cam)))

rot_cam2 = np.identity(4); rot_cam2[:3,:3] = rot_cam
camera_coordinates = rot_cam2.dot(points - t_cam_homogeneous)
camera_coordinates = camera_coordinates[:3,:] / camera_coordinates[3,:]

# Perform the projection
projected_points = intrinsics.dot(camera_coordinates)
projected_points = projected_points[:2,:] / projected_points[2,:]
projected_points[0,:] = -projected_points[0,:] # Inverted x-axis because camera is pointing toward [0, 0, 1]

fig = plt.figure()
ax = Axes3D(fig)
ax.scatter(points[0,:], points[1,:], points[2,:], label="Points")
ax.scatter(t_cam[0], t_cam[1], t_cam[2], c="red", label="Camera")
ax.set_xlabel("X axis"); ax.set_ylabel("Y axis"); ax.set_zlabel("Z axis")
plt.title("World coordinates")
plt.legend()
plt.savefig('world_coordinates.png', dpi=300, bbox_inches="tight")

fig = plt.figure()
ax = Axes3D(fig)
ax.scatter(camera_coordinates[0,:], camera_coordinates[1,:], camera_coordinates[2,:], label="Points")
ax.scatter(0, 0, 0, c="red", label="Camera")
ax.set_xlabel("X axis"); ax.set_ylabel("Y axis"); ax.set_zlabel("Z axis")
plt.title("Camera coordinates")
plt.legend()
plt.savefig('camera_coordinates.png', dpi=300, bbox_inches="tight")

plt.figure()
plt.scatter(projected_points[0,:], projected_points[1,:])
plt.xlabel("X axis"); plt.ylabel("Y axis")
plt.title("Image coordinates")
plt.savefig('image_coordinates.png', dpi=300, bbox_inches="tight")

plt.show()

This yield me those figures: Respectively: World coordinates, Camera coordinates, Camera coordinates rotated to fit slightly the camera orientation (note that here the camera vector goes toward the figure viewpoint, it does not "enter" the figure) and image coordinates.

So we see that the vertical coordinates for the bottom points are correctly on the middle row (240) and the points are on the right side of the image (horizontal value > 320).

I believe one bug you had is that you found negative X values, so you negated the focals (-f/Sxy) in the intrinsics matrix to compensate. The issue here is that we've assumed that the camera initially was pointing toward $[0, 0, 1]$ (otherwise the 160° rotation would not point toward the points). If you look at it that way, the $x$-axis increases when going toward the left, the inverse of this axis should be taken.

Both our results seems similar to me, just that you supposed a $[0, -1, 0]$ up vector for the camera (in fact both axes were mirrored because you negated both focals), and did the computation in mm instead of meters.