blender camera from 3x4 matrix

Question

In computer vision, the transformation from 3D world coordinates to pixel coordinates is often represented by a 3x4 (3 rows by 4 cols) matrix P as detailed below. Given P I need code to return a camera in Blender.

This is the inverse of my previous question 3x4 camera matrix from blender camera

Context where this shows up

This is useful for e.g. Augmented reality where 3x4 cameras are computed from real imagery using computer vision / structure from motion algorithms which are then used in CG to render registered synthetic models. Right now I want to employ multiple real images to be used as multiview reference for detailed 3D modeling.

Some Details

An image pixel $(u,v)$ is generated from world $(x,y,z)$ coordinates through a 3x4 matrix using projective coordinates:

$$kx = PX$$

where $k$ is a constant, $x = (u,v,1)^t$, $X = (X,Y,Z,W)^t$, $P$ can be decomposed as

$$P=K[I|0] \begin{bmatrix} R & T\\ 0^\text{T} & 1 \end{bmatrix} = K[R|T]$$

and $K$ can be decomposed as

$$K=\begin{bmatrix} \alpha_u & s & u_o\\ 0 & \alpha_v & v_o\\ 0 & 0 & 1 \end{bmatrix}$$

given that:

$$f=\text{focal length}$$ $$\alpha_u=\frac{f\times u\text{ pixels}}{\text{unit length}}=\frac{f}{\text{width of a pixel in world units}}$$ $$\alpha_v=\frac{f\times y\text{ pixels}}{\text{unit length}}=\frac{f}{\text{height of a pixel in world units}}$$ $$u_o=u\text{ coordinate of principle point}$$ $$v_o=v\text{ coordinate of principle point}$$ $$s_\theta=\text{skew factor}$$

I need to generate a Blender camera from this.

Preliminary Research

I could certainly detail this up further in terms of all the coordinate systems and work out the solution myself if need be, but I'm looking to reuse well-tested shared code.

Libmv conversion code from internal to blender camera: libmv/src/ui/tvr/tvr_document.h

Directly related to: What is blender's camera projection matrix model? 3x4 camera matrix from blender camera

Other related questions: How can I get the camera's projection matrix?, How to find image coordinates of the rendered vertex?

Answer 1

I wrote the function get_blender_camera_from_3x4_P to do this, listed below.

    # Input: P 3x4 numpy matrix
    # Output: K, R, T such that P = K*[R | T], det(R) positive and K has positive diagonal
    #
    # Reference implementations: 
    #   - Oxford's visual geometry group matlab toolbox 
    #   - Scilab Image Processing toolbox
    def KRT_from_P(P):
        N = 3
        H = P[:,0:N]  # if not numpy,  H = P.to_3x3()
    
        [K,R] = rf_rq(H)
          
        K /= K[-1,-1]
        
        # from http://ksimek.github.io/2012/08/14/decompose/
        # make the diagonal of K positive
        sg = numpy.diag(numpy.sign(numpy.diag(K)))
    
        K = K * sg
        R = sg * R
        # det(R) negative, just invert; the proj equation remains same:
        if (numpy.linalg.det(R) < 0):
           R = -R
        # C = -H\P[:,-1]
        C = numpy.linalg.lstsq(-H, P[:,-1])[0]
        T = -R*C
        return K, R, T
    
    # RQ decomposition of a numpy matrix, using only libs that already come with
    # blender by default
    #
    # Author: Ricardo Fabbri
    # Reference implementations: 
    #   Oxford's visual geometry group matlab toolbox 
    #   Scilab Image Processing toolbox
    #
    # Input: 3x4 numpy matrix P
    # Returns: numpy matrices r,q
    def rf_rq(P):
        P = P.T
        # numpy only provides qr. Scipy has rq but doesn't ship with blender
        q, r = numpy.linalg.qr(P[ ::-1, ::-1], 'complete')
        q = q.T
        q = q[ ::-1, ::-1]
        r = r.T
        r = r[ ::-1, ::-1]
    
        if (numpy.linalg.det(q) < 0):
            r[:,0] *= -1
            q[0,:] *= -1
        return r, q
    
    # Creates a blender camera consistent with a given 3x4 computer vision P matrix
    # Run this in Object Mode
    # scale: resolution scale percentage as in GUI, known a priori
    # P: numpy 3x4
    def get_blender_camera_from_3x4_P(P, scale):
        # get krt
        K, R_world2cv, T_world2cv = KRT_from_P(numpy.matrix(P))
    
        scene = bpy.context.scene
        sensor_width_in_mm = K[1,1]*K[0,2] / (K[0,0]*K[1,2])
        sensor_height_in_mm = 1  # doesn't matter
        resolution_x_in_px = K[0,2]*2  # principal point assumed at the center
        resolution_y_in_px = K[1,2]*2  # principal point assumed at the center
    
        s_u = resolution_x_in_px / sensor_width_in_mm
        s_v = resolution_y_in_px / sensor_height_in_mm
        # TODO include aspect ratio
        f_in_mm = K[0,0] / s_u
        # recover original resolution
        scene.render.resolution_x = resolution_x_in_px / scale
        scene.render.resolution_y = resolution_y_in_px / scale
        scene.render.resolution_percentage = scale * 100
    
        # Use this if the projection matrix follows the convention listed in my answer to
        # https://blender.stackexchange.com/questions/38009/3x4-camera-matrix-from-blender-camera
        R_bcam2cv = Matrix(
            ((1, 0,  0),
             (0, -1, 0),
             (0, 0, -1)))
    
        # Use this if the projection matrix follows the convention from e.g. the matlab calibration toolbox:
        # R_bcam2cv = Matrix(
        #     ((-1, 0,  0),
        #      (0, 1, 0),
        #      (0, 0, 1)))
    
        R_cv2world = R_world2cv.T
        rotation =  Matrix(R_cv2world.tolist()) @ R_bcam2cv
        location = -R_cv2world @ T_world2cv
    
        # create a new camera
        bpy.ops.object.add(
            type='CAMERA',
            location=location)
        ob = bpy.context.object
        ob.name = 'CamFrom3x4PObj'
        cam = ob.data
        cam.name = 'CamFrom3x4P'
     
        # Lens
        cam.type = 'PERSP'
        cam.lens = f_in_mm 
        cam.lens_unit = 'MILLIMETERS'
        cam.sensor_width  = sensor_width_in_mm
        ob.matrix_world = Matrix.Translation(location)@rotation.to_4x4()
    
        #     cam.shift_x = -0.05
        #     cam.shift_y = 0.1
        #     cam.clip_start = 10.0
        #     cam.clip_end = 250.0
        #     empty = bpy.data.objects.new('DofEmpty', None)
        #     empty.location = origin+Vector((0,10,0))
        #     cam.dof_object = empty
     
        # Display
        cam.show_name = True
        # Make this the current camera
        scene.camera = ob
        bpy.context.scene.update()
    
    def test2():
        P = Matrix([
        [2. ,  0. , - 10. ,   282.  ],
        [0. ,- 3. , - 14. ,   417.  ],
        [0. ,  0. , - 1.  , - 18.   ]
        ])
        # This test P was constructed as k*[r | t] where
        #     k = [2 0 10; 0 3 14; 0 0 1]
        #     r = [1 0 0; 0 -1 0; 0 0 -1]
        #     t = [231 223 -18]
        # k, r, t = KRT_from_P(numpy.matrix(P))
        get_blender_camera_from_3x4_P(P, 1)

#Tests

• I got an image with known 3x4 projection matrix P. This can be from a computer vision dataset (usually following the conventions of the Matlab camera calibration toolbox), tracked set, or by generating P from a blender camera as in 3x4 camera matrix from blender camera
• I ran the above algorithm for such P, obtaining a camera in Blender.
• The dataset I tested with has 3D reconstruction data as well
• Numpad 0 gives the camera view. I confirmed that it looks like the photo
• By adding a background image to the camera, I noticed that the 3D model projection is spot-on.

#Remarks

• This code is meant as a basis for a general solution. Modify this if you follow different conventions.
• Only unit aspect ratio (render option) is supported for now.
• The camera representation might sometimes look odd, but the only thing that matters is that it generates the right image.