3x4 camera matrix from blender camera
I wrote the function get_3x4_RT_matrix_from_blender to do this, listed below.
import bpy
import bpy_extras
from mathutils import Matrix
from mathutils import Vector

#---------------------------------------------------------------
# 3x4 P matrix from Blender camera
#---------------------------------------------------------------

# Build intrinsic camera parameters from Blender camera data
#
# See notes on this in 
# blender.stackexchange.com/questions/15102/what-is-blenders-camera-projection-matrix-model
def get_calibration_matrix_K_from_blender(camd):
    f_in_mm = camd.lens
    scene = bpy.context.scene
    resolution_x_in_px = scene.render.resolution_x
    resolution_y_in_px = scene.render.resolution_y
    scale = scene.render.resolution_percentage / 100
    sensor_width_in_mm = camd.sensor_width
    sensor_height_in_mm = camd.sensor_height
    pixel_aspect_ratio = scene.render.pixel_aspect_x / scene.render.pixel_aspect_y
    if (camd.sensor_fit == 'VERTICAL'):
        # the sensor height is fixed (sensor fit is horizontal), 
        # the sensor width is effectively changed with the pixel aspect ratio
        s_u = resolution_x_in_px * scale / sensor_width_in_mm / pixel_aspect_ratio 
        s_v = resolution_y_in_px * scale / sensor_height_in_mm
    else: # 'HORIZONTAL' and 'AUTO'
        # the sensor width is fixed (sensor fit is horizontal), 
        # the sensor height is effectively changed with the pixel aspect ratio
        pixel_aspect_ratio = scene.render.pixel_aspect_x / scene.render.pixel_aspect_y
        s_u = resolution_x_in_px * scale / sensor_width_in_mm
        s_v = resolution_y_in_px * scale * pixel_aspect_ratio / sensor_height_in_mm
    

    # Parameters of intrinsic calibration matrix K
    alpha_u = f_in_mm * s_u
    alpha_v = f_in_mm * s_v
    u_0 = resolution_x_in_px * scale / 2
    v_0 = resolution_y_in_px * scale / 2
    skew = 0 # only use rectangular pixels

    K = Matrix(
        ((alpha_u, skew,    u_0),
        (    0  , alpha_v, v_0),
        (    0  , 0,        1 )))
    return K

# Returns camera rotation and translation matrices from Blender.
# 
# There are 3 coordinate systems involved:
#    1. The World coordinates: &quot;world&quot;
#       - right-handed
#    2. The Blender camera coordinates: &quot;bcam&quot;
#       - x is horizontal
#       - y is up
#       - right-handed: negative z look-at direction
#    3. The desired computer vision camera coordinates: &quot;cv&quot;
#       - x is horizontal
#       - y is down (to align to the actual pixel coordinates 
#         used in digital images)
#       - right-handed: positive z look-at direction
def get_3x4_RT_matrix_from_blender(cam):
    # bcam stands for blender camera
    R_bcam2cv = Matrix(
        ((1, 0,  0),
         (0, -1, 0),
         (0, 0, -1)))

    # Transpose since the rotation is object rotation, 
    # and we want coordinate rotation
    # R_world2bcam = cam.rotation_euler.to_matrix().transposed()
    # T_world2bcam = -1*R_world2bcam * location
    #
    # Use matrix_world instead to account for all constraints
    location, rotation = cam.matrix_world.decompose()[0:2]
    R_world2bcam = rotation.to_matrix().transposed()

    # Convert camera location to translation vector used in coordinate changes
    # T_world2bcam = -1*R_world2bcam*cam.location
    # Use location from matrix_world to account for constraints:     
    T_world2bcam = -1*R_world2bcam @ location

    # Build the coordinate transform matrix from world to computer vision camera
    # NOTE: Use * instead of @ here for older versions of Blender
    # TODO: detect Blender version
    R_world2cv = R_bcam2cv@R_world2bcam
    T_world2cv = R_bcam2cv@T_world2bcam

    # put into 3x4 matrix
    RT = Matrix((
        R_world2cv[0][:] + (T_world2cv[0],),
        R_world2cv[1][:] + (T_world2cv[1],),
        R_world2cv[2][:] + (T_world2cv[2],)
         ))
    return RT

def get_3x4_P_matrix_from_blender(cam):
    K = get_calibration_matrix_K_from_blender(cam.data)
    RT = get_3x4_RT_matrix_from_blender(cam)
    return K@RT, K, RT

# ----------------------------------------------------------
# Alternate 3D coordinates to 2D pixel coordinate projection code
# adapted from https://blender.stackexchange.com/questions/882/how-to-find-image-coordinates-of-the-rendered-vertex?lq=1
# to have the y axes pointing up and origin at the top-left corner
def project_by_object_utils(cam, point):
    scene = bpy.context.scene
    co_2d = bpy_extras.object_utils.world_to_camera_view(scene, cam, point)
    render_scale = scene.render.resolution_percentage / 100
    render_size = (
            int(scene.render.resolution_x * render_scale),
            int(scene.render.resolution_y * render_scale),
            )
    return Vector((co_2d.x * render_size[0], render_size[1] - co_2d.y * render_size[1]))

# ----------------------------------------------------------
if __name__ == &quot;__main__&quot;:
    # Insert your camera name here
    cam = bpy.data.objects['Camera.001']
    P, K, RT = get_3x4_P_matrix_from_blender(cam)
    print(&quot;K&quot;)
    print(K)
    print(&quot;RT&quot;)
    print(RT)
    print(&quot;P&quot;)
    print(P)

    print(&quot;==== Tests ====&quot;)
    e1 = Vector((1, 0,    0, 1))
    e2 = Vector((0, 1,    0, 1))
    e3 = Vector((0, 0,    1, 1))
    O  = Vector((0, 0,    0, 1))

    p1 = P @ e1
    p1 /= p1[2]
    print(&quot;Projected e1&quot;)
    print(p1)
    print(&quot;proj by object_utils&quot;)
    print(project_by_object_utils(cam, Vector(e1[0:3])))

    p2 = P @ e2
    p2 /= p2[2]
    print(&quot;Projected e2&quot;)
    print(p2)
    print(&quot;proj by object_utils&quot;)
    print(project_by_object_utils(cam, Vector(e2[0:3])))

    p3 = P @ e3
    p3 /= p3[2]
    print(&quot;Projected e3&quot;)
    print(p3)
    print(&quot;proj by object_utils&quot;)
    print(project_by_object_utils(cam, Vector(e3[0:3])))

    pO = P @ O
    pO /= pO[2]
    print(&quot;Projected world origin&quot;)
    print(pO)
    print(&quot;proj by object_utils&quot;)
    print(project_by_object_utils(cam, Vector(O[0:3])))
    
    # Bonus code: save the 3x4 P matrix into a plain text file
    # Don't forget to import numpy for this
    nP = numpy.matrix(P)
    numpy.savetxt(&quot;/tmp/P3x4.txt&quot;, nP)  # to select precision, use e.g. fmt='%.2f'

Detailed Tests

I included some helper functions and tests to cross-check the 3D-2D projection process with another approach using world_to_camera_view from @ideasman42's answer to How to find image coordinates of the rendered vertex? The floating point coordinates match up perfectly up to many decimal places. To get the integer pixel coordinates after projection, just round them.

I tested this visually and analytically with a few scenes, to make sure both routines were exactly correct on actual rendered images.

I represented e1, e2, e3 and O with small cone objects whose tips were at the desired positions
and coded with different colors: Red for X, Green for Y, Blue for Z
I placed the camera at a convenient location at first,
R = identity, T = (0,0,5)
I used default Blender parameters for the camera intrinsics
I computed the P matrix by hand and computed it using my rountines and they match perfectly
Visually comparing both the Cycles render and the OpenGL render in GIMP show that the produced coordinates match well visually (up to perhaps +- 1 pixel error, which is the sort of thing this careful testing was carried out anyways - to be investigated).
I then arbitrarily rotated and translated the camera such that the world e1, e2, e3 remain visible and at general locations. The projections match visually in the rendered images, and also match to the ones given by the object utils routine given above.



Here is an example of a test output:
  K
&lt;Matrix 3x3 (2100.0000,    0.0000, 960.0000)
            (   0.0000, 2100.0000, 540.0000)
            (   0.0000,    0.0000,   1.0000)&gt;
RT
&lt;Matrix 3x4 (-0.0594, -0.9483, -0.3118, -0.6837)
            ( 0.6234,  0.2087, -0.7536, -0.1887)
            ( 0.7797, -0.2391,  0.5787,  4.2599)&gt;
P
&lt;Matrix 3x4 ( 623.8236, -2220.9482,   -99.1637, 2653.7441)
            (1730.0916,   309.2392, -1269.9426, 1904.1343)
            (   0.7797,    -0.2391,     0.5787,    4.2599)&gt;
==== Tests ====
Projected e1
&lt;Vector (650.3716, 721.1436, 1.0000)&gt;
proj by object_utils
&lt;Vector (650.3716, 721.1436)&gt;
Projected e2
&lt;Vector (107.6402, 550.4857, 1.0000)&gt;
proj by object_utils
&lt;Vector (107.6403, 550.4857)&gt;
Projected e3
&lt;Vector (527.9586, 131.0692, 1.0000)&gt;
proj by object_utils
&lt;Vector (527.9586, 131.0693)&gt;
Projected world origin
&lt;Vector (622.9655, 446.9949, 1.0000)&gt;
proj by object_utils
&lt;Vector (622.9656, 446.9949)&gt;

Here is the test set with the blender rotation and location for the camera.

Here is the OpenGL render. The coordinates output above match the tip of the cones visually (up to +- 1 px error, perhaps due to aliasing).

Limitations: this code currently does not support certain intrinsic camera parameter configurations, see my answer to What is blender&#39;s camera projection matrix model? for a list of the limitations to get_calibration_matrix_K_from_blender.
Object coordinates: to deal with object coordinates, first transform to world coordinates, see my answer to How to find image coordinates of the rendered vertex?
Matrix multiplication used to be * before Blender 2.8, newer versions use @ instead