Issues with manual culling, using BMesh
I have tried two implementations, and I wasn't too satisfied. I should say that the second implementation is definitely the better of the two since I need to simply cull faces that block the camera (think camera outside facing a wall of a room).

The most accurate seems to be in NDC space. For every face, project each vertex, transform to NDC space, then construct a normal from the 2D triangle for that ace. If the z component of the normal is negative, cull face.

The important bit here is that we do not clip in either clip or NDC space. Clipping is necessary when vertices are behind the camera (in homogeneous coordinates, w ≤ 0). The result is the vertex being "mirrored" in front of the camera, resulting in normals flipping! So we take a shortcut, forget clipping altogether, and simply move the vertices as close to w as possible (1.0e-05).

projection_matrix = camera_obj.calc_matrix_camera(scene.render.resolution_x,
                                                  scene.render.resolution_y,
                                                  scene.render.pixel_aspect_x,
                                                  scene.render.pixel_aspect_y)
project_matrix = projection_matrix * (camera_obj.matrix_world.inverted() * (object_obj.matrix_world))
delete_faces = []
mesh = object_obj.data
new_bmesh = bmesh.new()
new_bmesh.from_mesh(mesh)
for face in new_bmesh.faces:
    face_ndc = []
    for vert in face.verts:
        # Project point
        pvert_4d = project_matrix * vert.co.to_4d()
        # Perspective divide; transform to viewport space (NDC).

        # We need to consider the case where where w &lt; 0.
        # For each 4D vertex, if w &lt; 0, the vertex will be
        # flipped when performing the perspective divide,
        # resulting in the normal being flipped.

        # We could do go through and implement clipping using
        # the Liang-Barsky line clipping algorithm, but we only
        # care about backface culling. Clamping w to 0.0 (close
        # to) does result in distorting the triangle, but the
        # normal will remain the same as the normal's magnitude
        # is proportional to the area of the triangle, but
        # independent of the direction of the normal.
        w = 1.0e-05 if (pvert_4d.w &lt; 0.0) else pvert_4d.w
        ndc_vert = pvert_4d / w
        face_ndc.append(ndc_vert)
    # Calculate the normal
    z = calc_normal(face_ndc)
    if z &lt; 0.0:
        delete_faces.append(face)
bmesh.ops.delete(new_bmesh, geom = delete_faces, context = 5)
new_bmesh.to_mesh(mesh)
mesh.update()