Camera extrinsics calculation wrong

I am trying to find camera extrinsics from 2 images. I have the intrinsics from CameraCalibration and the scene has known dimensions (created with 3DSMAX).

The chessboard is 1000*1000, each square thus 125*125. The camera is at (0,0,3000) looking straight down at the chessboard which is centred at the origin. In the second image, the camera is translated (-1500, 0, -402) and rotated 30° on the Y axis to point again to the centre of the chessboard:

GoodFeaturesToTrack correctly identifies the 81 corners:

I create the 3d points of the chessboard corners, cvFindExtrinsicCameraParams2 to compute the intrinsics and cvRodrigues2 to get the rotation matrix. Here's the code

Imports Emgu.CV
Imports Emgu.CV.Structure
Imports Emgu.CV.CvInvoke
Imports Emgu.CV.CvEnum
Imports Emgu.CV.UI
Imports System.Drawing
Imports System.IO
Imports System.Diagnostics
Imports System.Math
Module main_

    Sub Main()

        Const MAXFEATURES As Integer = 100
        Dim featuresA(0)() As PointF
        Dim featuresB(0)() As PointF
        Dim features As Integer = 0
        Dim imgA As Emgu.CV.Image(Of Emgu.CV.Structure.Bgr, Byte)
        Dim imgB As Emgu.CV.Image(Of Emgu.CV.Structure.Bgr, Byte)
        Dim grayA As Emgu.CV.Image(Of Emgu.CV.Structure.Gray, Byte)
        Dim grayB As Emgu.CV.Image(Of Emgu.CV.Structure.Gray, Byte)
        Dim pyrBufferA As Emgu.CV.Image(Of Emgu.CV.Structure.Gray, Byte)
        Dim pyrBufferB As Emgu.CV.Image(Of Emgu.CV.Structure.Gray, Byte)
        Dim pointsA As Matrix(Of Single)
        Dim pointsB As Matrix(Of Single)
        Dim flags As Emgu.CV.CvEnum.LKFLOW_TYPE = Emgu.CV.CvEnum.LKFLOW_TYPE.DEFAULT
        Dim imagesize As Size
        Dim termcrit As New Emgu.CV.Structure.MCvTermCriteria(20, 0.03D)
        Dim status As Byte() = Nothing
        Dim errors As Single() = Nothing
        Dim red As Bgr = New Bgr(Color.Red)

        ' Load chessboards
        imgA = New Image(Of [Structure].Bgr, Byte)("chessboard centre.jpg")
        imgB = New Image(Of [Structure].Bgr, Byte)("chessboard left.jpg")
        grayA = imgA.Convert(Of Gray, Byte)()
        grayB = imgB.Convert(O开发者_JAVA技巧f Gray, Byte)()

        ' setup for feature detection
        imagesize = cvGetSize(grayA)
        pyrBufferA = New Emgu.CV.Image(Of Emgu.CV.Structure.Gray, Byte)(imagesize.Width + 8, imagesize.Height / 3)
        pyrBufferB = New Emgu.CV.Image(Of Emgu.CV.Structure.Gray, Byte)(imagesize.Width + 8, imagesize.Height / 3)
        features = MAXFEATURES

        ' Find features
        featuresA = grayA.GoodFeaturesToTrack(features, 0.01, 25, 3)
        grayA.FindCornerSubPix(featuresA, New System.Drawing.Size(10, 10), New System.Drawing.Size(-1, -1), termcrit)
        features = featuresA(0).Length

        ' Compute optical flow. Not necessary here but good to remember
        Emgu.CV.OpticalFlow.PyrLK(grayA, grayB, pyrBufferA, pyrBufferB, featuresA(0), New Size(25, 25), 3, termcrit, flags, featuresB(0), status, errors)
        Debug.Assert(featuresA(0).GetUpperBound(0) = featuresB(0).GetUpperBound(0))

        ' Copy features to an easier-to-use matrix and get min/max to create 3d points
        Dim minx As Double = Double.MaxValue
        Dim miny As Double = Double.MaxValue
        Dim maxx As Double = Double.MinValue
        Dim maxy As Double = Double.MinValue
        pointsA = New Matrix(Of Single)(features, 2)
        pointsB = New Matrix(Of Single)(features, 2)
        For i As Integer = 0 To features - 1
            pointsA(i, 0) = featuresA(0)(i).X
            pointsA(i, 1) = featuresA(0)(i).Y
            pointsB(i, 0) = featuresB(0)(i).X
            pointsB(i, 1) = featuresB(0)(i).Y
            If pointsA(i, 0) < minx Then
                minx = pointsA(i, 0)
            End If
            If pointsA(i, 1) < miny Then
                miny = pointsA(i, 1)
            End If
            If pointsA(i, 0) > maxx Then
                maxx = pointsA(i, 0)
            End If
            If pointsA(i, 1) > maxy Then
                maxy = pointsA(i, 1)
            End If
        Next

        ' Create 3D object points that correspond to chessboard corners
        ' (The chessboard is 1000*1000, squares are 125*125)
        Dim corners As Integer = Sqrt(features)
        Dim obj As New Matrix(Of Double)(features, 3)
        Dim squaresize2dx As Double = (maxx - minx) / 8 ' pixel width of a chessboard square
        Dim squaresize2dy As Double = (maxy - miny) / 8 ' pixel height of a chessboard square
        For i As Integer = 0 To features - 1
            obj(i, 0) = Math.Round((pointsA(i, 0) - minx) / squaresize2dx) * 125 ' X=0, 125, 250, 375 ... 1000
            obj(i, 1) = Math.Round((pointsA(i, 1) - miny) / squaresize2dy) * 125 ' idem in Y
            obj(i, 2) = 0
            ' Debug.WriteLine(pointsA(i, 0) & " " & pointsA(i, 1) & " " & obj(i, 0) & " " & obj(i, 1) & " " & obj(i, 2)) ' Just to verify
        Next

        ' These were calculated with CalibrateCamera using the same images
        Dim intrinsics As New Matrix(Of Double)({{889.1647, 0.0, 318.3721},
                                                 {0.0, 888.5134, 238.4254},
                                                 {0.0, 0.0, 1.0}})
        ' Find extrinsics
        Dim distortion As New Matrix(Of Double)({-0.036302, 2.008797, -29.674306, -29.674306})
        Dim translation As New Matrix(Of Double)(3, 1)
        Dim rotation As New Matrix(Of Double)(3, 1)
        cvFindExtrinsicCameraParams2(obj, pointsA, intrinsics, distortion, rotation, translation, False)

        ' Convert rotation vector to rotation matrix
        Dim rotmat As New Matrix(Of Double)(3, 3)
        Dim jacobian As New Matrix(Of Double)(9, 3)
        cvRodrigues2(rotation, rotmat, jacobian)

        ' From http://en.wikipedia.org/wiki/Rotation_representation paragraph "Conversion formulae between representations"
        Dim yr As Double = Asin(-rotmat(2, 0))
        Dim xr As Double = Asin(rotmat(2, 1) / Cos(yr))
        Dim zr As Double = Asin(rotmat(1, 0) / Cos(yr))

    End Sub
End Module

The results don't seem correct, I was expecting the translation/rotation but i get this:

translation
208.394425348956 
-169.447506344527 
-654.273807995629 
rotation
-0.0224937226554797 
-2.13660350939653 
-1.10542281290682 
rotmat
-0.741100224945266 0.322885083546921 -0.588655824237707 
-0.293966101915684 0.632206237134128 0.716867633983572 
0.603617749499279 0.704315622822328 -0.373610915174894 
xr=1.08307908108382 yr=-0.648031006135158 zr=-0.377625254910525
xr=62.0558602250101° yr=-37.1294416451609° zr=-21.636333343925°

Would someone have an idea what I'm doing wrong? Thanks!

Found it. The distortion coefficients are

k1, k2, p1, p2, k3

and not

k1, k2, k3, k4, k5

as I had coded. When I set them to

-0.05716, 2.519006, -0.001674, -0.001021, -33.372798

the answer is (roughly) correct

Please look at the pin hole camera equation:

[u, v, 1]' ~ A * R|T * [x, y, z, 1]', where extrinsic matrix R|T has dimensions 3x4. Note that multiplying it with the 'ideal' or vanishing point at infinity such as [1,0,0,0] will give you a corresponding column of R|T. This simply means that to find all columns of R|T you have to have at least three vanishing points that can be easily found from your checkerboard patterns.

For now you have only 1, so you simply got lucky if your results look reasonable. Try again and select minimum 10-20 different slants, distances and tilts of the calibration rig.