Manually modeling and rendering a 3D object to 2D using numpy with extrinsic and intrinsic params

Question

My current goal is to learn the math behind classic computer vision by attempting to re-create from scratch the animation produced here via modeling the paper in world coordinates and then computing the image points using intrinsic and extrinsic parameters, and to do so with multiple types of curves.

I started by rendering a grid of points in the world coordinate system. I am using a pinhole model with a principle point at the c enter of the film/sensor. But I can't succeed to rotate the perspective with this code as once I use any angle other than 0 nothing renders anymore. More importantly, I am generally unconvinced that there isn't something fundamentally missing with my understanding or approach. I do not want to proceed to add curvature before I get rotation working. Thank you so much for your time to help.

import cv2, math
import numpy as np

def rotationMatrix(roll, yaw, pitch):
    # Convert degrees to radians
    roll = roll * (math.pi/180)
    yaw = yaw * (math.pi/180)
    pitch = pitch * (math.pi/180)

    # pitch is beta, gamma is roll, alpha must be yaw
    alpha = yaw
    beta = pitch
    gamma = roll
   
    cos, sin = math.cos, math.sin

    # Calculate rotation matrix using https://en.wikipedia.org/wiki/Rotation_matrix#General_rotations
    rotationMatrix = [[cos(alpha)*cos(beta), cos(alpha)*sin(beta)*sin(gamma)-sin(alpha)*cos(gamma), cos(alpha)*sin(beta)*cos(gamma)+sin(alpha)*sin(gamma)],
                      [sin(alpha)*cos(beta), sin(alpha)*sin(beta)*sin(gamma)+cos(alpha)*cos(gamma), sin(alpha)*sin(beta)*cos(gamma)-cos(alpha)*sin(gamma)],
                      [-sin(beta), cos(beta)*sin(gamma), cos(beta)*cos(gamma)]]

    # Print rotation matrix
    print("Rotation Matrix:")
    for i in range(3):
        for j in range(3):
            print(rotationMatrix[i][j], end=" ")
        print()
    return rotationMatrix

w, h = 100, 200
f = 200
img = np.zeros((h, w, 3))
world_points = []
for x in range(0, w, 2):
  for y in range(0, h, 2):
    p = (x/f - (w//2), y/f - h//2)
    world_points.append(p)

K = np.array([[f, 0, w//2, 0],
              [0, f, h//2, 0],
              [0, 0, 1, 0]], dtype=np.float32)
R = rotationMatrix(0,0,0)
tvec = np.array([0, 0, f], dtype=np.float32)
Mex = np.array([[R[0][0], R[0][1], R[0][2], tvec[0]],
                [R[1][0], R[1][1], R[1][2], tvec[1]],
                [R[2][0], R[2][1], R[2][2], tvec[2]],
                [0, 0, 0, 1]], dtype=np.float32)
print('K', K)
print('Mex', Mex)
for p in world_points:
    pnew = K @ Mex @ np.array([p[0], p[1], 0, 1])
    cv2.circle(img, (int(pnew[0]), int(pnew[1])), 0, (255, 255, 0), -1)

plt.imshow(img)
plt.show()