Q:

Issues with manual culling, using BMesh

I have a camera and a monkey (suzanne) in the scene.  I'm attempting to manually cull the back faces from the monkey. 
The problem is that certain faces get culled, despite its normal seen from the camera's view port. 

The monkey in object mode. 

The monkey after culling (using the function below). 

I transform the monkey's vertices into world space, then transform the vertices to face the camera (camera's inverse matrix).  This assumes the camera is at (0,0,0). 

I'm using the internal renderer.  Blender 2.78. 
def cull_manual(camera_obj, object_obj):
    matrix = camera_obj.matrix_local.inverted() * object_obj.matrix_local

    delete_faces = []

    mesh = object_obj.data
    new_bmesh = bmesh.from_edit_mesh(mesh)

    for face in new_bmesh.faces:
        verts = [(matrix * bmvert.co) for bmvert in face.verts]
        a = (verts[1] - verts[0])
        b = (verts[2] - verts[0]) if (len(verts) == 3) else (verts[3] - verts[0])
        normal = a.cross(b).normalized()

        v = verts[0].normalized()

        visibility = (-v).dot(normal)
        if visibility <= 0.0:
            delete_faces.append(face)

    bmesh.ops.delete(new_bmesh, geom = delete_faces, context = 5)
    bmesh.update_edit_mesh(mesh)

camera_obj = find_obj(scene, "Camera.05")
object_obj = find_obj(scene, "Suzanne")
cull_manual(camera_obj, object_obj)

I devised a way to do it the opposite way, where I move the object to the origin, and move the camera by applying the object's inverse transform:
def cull_from_camera(camera_obj, object_obj):
    matrix = object_obj.matrix_world.inverted()
    camera_matrix = matrix * camera_obj.matrix_world
    location, _, _ = camera_matrix.decompose()
    from_camera = location.normalized()
    delete_faces = []
    mesh = object_obj.data
    new_bmesh = bmesh.new()
    new_bmesh.from_mesh(mesh)
    for face in new_bmesh.faces:
        visibility = from_camera.dot(face.normal)
        E = 1.0e-05
        R = -0.5 # HERE
        if visibility <= (R + E):
            delete_faces.append(face)
    bmesh.ops.delete(new_bmesh, geom = delete_faces, context = 5)
    new_bmesh.to_mesh(mesh)
    mesh.update()

It works, but there are still some aggressive culling happening, unless I change the visibility check from 0.0 to -0.5. 
I can't think of a way to truly cull only the faces that are not visible to the camera. 

A:

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 < 0. 
        # For each 4D vertex, if w < 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 < 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 < 0.0:
        delete_faces.append(face)
bmesh.ops.delete(new_bmesh, geom = delete_faces, context = 5)
new_bmesh.to_mesh(mesh)
mesh.update()

