How does the calibration method by Toscani and Faugeras work?

Question

First, I tried a basic approach to get the parameters of the camera. In the matrix A, you have the two equations that a 3D point can give (multiplied by the number of points): Because we have uZ -uoZ -px =0 and vZ -v0Z - py=0 Then we have A*X=0 where X represents the unknows and A represents the observables.

To solve X, in an ideal scenario, we would perform the Singular Value Decomposition (SVD) of A. X is associated with the lowest value in D. However, there's an issue: X = 0 is a valid solution! So I think the SVD gives that solution, rather than the one we're looking for. I’m not sure if there is a problem in my code or if this is just "normal". objpoints is in the reference frame of the checkerboard, and imgpoints in the reference frame the camera. (they are computed below)

    A = []
    for i in range(len(objpoints)):
        X, Y, Z = objpoints[i]
        u, v = imgpoints[i][0]

        row1 = [u*X, u*Y, u*Z, -X, -Y, -Z, 0, 0, 0, -1, 0, u]
        row2 = [v*X, v*Y, v*Z, 0, 0, 0, -X, -Y, -Z, 0, -1, v]
        
        A.append(row1)
        A.append(row2)

    A = np.array(A)

    # Solve the homogeneous system using SVD
    _, _, V = np.linalg.svd(A)
    X = V[:, -1]
    print("Estimated Camera Parameters:", X)

    print("uo:",X[0]*X[3]+X[1]*X[4]+X[2]*X[6])

In the method by Toscani and Faugeras, it’s a bit different. This time we don't solve AX because 0 can be a solution. Instead, we try to minimize the following cost function C=||AX1 + B*X2||2 + (1-X1) ; In this approach, the problem has been decomposed.

A= [[uX, uY, uZ]**,[vX, vY, vZ]] and X1=[X1,X2,X3] (the 3 first parameters we obtained before) B=[[-X, -Y, -Z, 0, 0, 0, -1, 0, u],[ 0, 0, 0, -X, -Y, -Z, 0, -1, v]] and X2=[X4,..,X12] (the rest)

By differentiating C, we can obtain the **solution **to the problem. Derivative and then: Solution

As you can see, we need to compute the inverse of B.T @ B in order to obtain some parameters (not all of them!). However, I have a BIG problem: B.T @ B is supposed to be invertible. When I look at the method, there seems to be no reason why it shouldn’t be! There is no apparent dependency between rows and columns!

Yet, in the code I wrote (which you can find below), I get an** error stating that the matrix is singular**. I looked into it, and indeed, the 9x9 matrix is not of full rank (rank 8, to be precise). This makes it impossible to calculate X1 and X2. I’m not sure if I’ve made a mistake in my calculations, the matrix handling, or if this method is just meant to be used in a specific context. For now, I’m trying to make it work on a checkerboard image from the OpenCV calibration tutorial.

Here’s the code in question. It’s supposed to display:Cx= X[0]*X[3]+X[1]*X[4]+X[2]*X[6] which represents the center of the image along with Cy However, I ended up with a value around 10^(-3) , which is far from the expected 400 given by OpenCV.

import numpy as np
import cv2
import glob


CHECKERBOARD = (6, 9)  
square_size = 3 

def calibration(CHECKERBOARD,square_size):
    objp = np.zeros((CHECKERBOARD[0] * CHECKERBOARD[1], 3), np.float32)
    objp[:, :2] = np.mgrid[0:CHECKERBOARD[0], 0:CHECKERBOARD[1]].T.reshape(-1, 2) * square_size
    objpoints = []  
    imgpoints = [] 

    images = glob.glob('Folder/*.jpg') 
    #Only for one image. I have to remove the loop
    for fname in images: 
        img = cv2.imread(fname)
        gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
        ret, corners = cv2.findChessboardCorners(gray, CHECKERBOARD, cv2.CALIB_CB_ADAPTIVE_THRESH + cv2.CALIB_CB_FAST_CHECK + cv2.CALIB_CB_NORMALIZE_IMAGE)
        if ret:
            objpoints=objp
            imgpoints=corners
            cv2.drawChessboardCorners(img, CHECKERBOARD, corners, ret)
            cv2.imshow('Image', img)
            cv2.waitKey()  
        else:
            print(ret)

    cv2.destroyAllWindows()

    A = np.zeros((CHECKERBOARD[0] * CHECKERBOARD[1] * 2, 3), np.float32)  
    B = np.zeros((CHECKERBOARD[0] * CHECKERBOARD[1] * 2, 9), np.float32)
    #Algo 1: using Toscani Faugeras approach
    for i in range (CHECKERBOARD[0] * CHECKERBOARD[1]):
        X,Y,Z= objpoints[i]
        u,v=imgpoints[i][0]
        
        A[2*i][0]=u*X
        A[2*i][1]=u*Y
        A[2*i][2]=u*Z

        A[2*i+1][0]=v*X
        A[2*i+1][1]=v*Y
        A[2*i+1][2]=v*Z

        B[2*i][0]=-X
        B[2*i][1]=-Y
        B[2*i][2]=-Z
        #B[2*i][3]=0
        #B[2*i][4]=0
        #B[2*i][5]=0
        B[2*i][6]=-1
        B[2*i][7]=0
        B[2*i][8]=u
        
        #B[2*i+1][0]=0
        #B[2*i+1][1]=0
        #B[2*i+1][2]=0
        B[2*i+1][3]=-X
        B[2*i+1][4]=-Y
        B[2*i+1][5]=-Z
        B[2*i+1][6]=0
        B[2*i+1][7]=-1
        B[2*i+1][8]=v
    
    print(B) #Is there a problem? 
    E = A.T @ A - A.T @ B @ np.linalg.inv(B.T @ B) @ B.T @ A 
    eigenvalues, eigenvectors = np.linalg.eig(E)
    X1 = eigenvectors[:, np.argmin(eigenvalues)]
    X2 = -np.linalg.pinv(B.T @ B) @ B.T @ A @ X1
    print("X1",X1,"X2",X2)
    #print([i for i in range (9) if ((B.T @ B)[i,i]==0)] ) #2 erreurs
    #print(np.linalg.matrix_rank((B.T @ B)))

    print(X[0]*X[3]+X[1]*X[4]+X[2]*X[6])

calibration(CHECKERBOARD,square_size)

I couldn't find anything on the web on this. It seems that studying this is quite challenging.

If anyone sees an issue in the code or the approach, thanks! I know it's not the best method, but I wanted to try it out to compare the results with the OpenCV functions.