Stereo rectification done wrong :(

As others before me I am getting frustrated trying to calibrate a pair of stereo cameras.

I am using opencv 4.9.0 and a solid Charuco board (18x12) where each squares is about 6cm each. I took 30 images from different distances and position.

The calibration of each camera works fine and leads to a errors of 0.5. but the stereo calibration error is about 20.

left images:






























right images:




























my images look fine to me and i could check that the corners found between each frame correspond to the same ones and to the good Charuco IDs:


Obviously with such stereo calibration error the rectified maps are useless and end up zooming in a tiny portion of the images. I show here a stack of the original/ “rectified” images:
image

Last and not least, here is the code I am using:

import cv2 #version 4.9.0
import numpy as np
from glob import glob
import matplotlib.pyplot as plt



def read_chessboards(images,board):
    found_corners = []
    allIDs = []
    # Ssub pixel corner detection criteria
    criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 100, 0.00001)

    
    count=0
    img_ind=[]
    for gray in images:
        if len(gray.shape)==3:
            gray = cv2.cvtColor(gray, cv2.COLOR_BGR2GRAY)
        corners, ids, rejectedImgPoints = cv2.aruco.detectMarkers(gray, aruco_dictionary)

        if len(corners)>0:# if some corners are found we go through subpixel detection
            # Ssub pixel resolution of the corners
            for corner in corners:
                cv2.cornerSubPix(gray, corner,
                                 winSize = (3,3),
                                 zeroZone = (-1,-1),
                                 criteria = criteria)
            res2 = cv2.aruco.interpolateCornersCharuco(corners,ids,gray,board)
            if res2[1] is not None and res2[2] is not None and len(res2[1])>3 and count%1==0:
                found_corners.append(res2[1])
                allIDs.append(res2[2])
                img_ind.append(count)
        else:
            print('Not enough corners found in image ' +str(count) +'.')
        count+=1
        imsize = gray.shape
    return found_corners,allIDs,imsize

#allCorners,allIds,imsize=read_chessboards(images)

def calibrate_camera(allCorners,allIds,imsize,board):

    cameraMatrixInit = np.array([[ 1000.,    0., imsize[0]/2.],
                                 [    0., 1000., imsize[1]/2.],
                                 [    0.,    0.,           1.]])

    distCoeffsInit = np.zeros((5,1))
    flags = (cv2.CALIB_USE_INTRINSIC_GUESS + cv2.CALIB_RATIONAL_MODEL + cv2.CALIB_FIX_ASPECT_RATIO)
    #flags = (cv2.CALIB_RATIONAL_MODEL)
    (ret, camera_matrix, distortion_coefficients0,
     rotation_vectors, translation_vectors,
     stdDeviationsIntrinsics, stdDeviationsExtrinsics,
     perViewErrors) = cv2.aruco.calibrateCameraCharucoExtended(
                      charucoCorners=allCorners,
                      charucoIds=allIds,
                      board=board,
                      imageSize=imsize,
                      cameraMatrix=cameraMatrixInit,
                      distCoeffs=distCoeffsInit,
                      flags=flags,
                      criteria=(cv2.TERM_CRITERIA_EPS & cv2.TERM_CRITERIA_COUNT, 10000, 1e-9))
    print('Camera calibration error: ' ,ret )
    return ret, camera_matrix, distortion_coefficients0, rotation_vectors, translation_vectors

#%%
def calibrate_stereo(board,left_Corners,left_IDs,left_imsize,left_camera_matrix, left_distortion_coefficients,right_Corners,right_IDs,right_imsize,right_camera_matrix, right_distortion_coefficients,left_img=None,right_img=None):
## calibrate stereo paright
 
    allCorners=board.getChessboardCorners() #array of realworld point 3D

    #sort and keep corners onl if they appear in the 2 set of images right and left
    left_corners_sampled = []
    right_corners_sampled = []
    world_corners_sampled = []
    left_IDs_sampled=[]
    right_IDs_sampled=[]
    # loop over all left ID/corners: to sort the corners found in each frames
    for i in range(len(left_Corners)):# loop over pair of images
        frame_left_corners=[]
        frame_right_corners=[]
        frame_world_corners=[]
        frame_left_IDs=[]
        frame_right_IDs=[]
        for ii in range(len(left_IDs[i])):# loop over left corners of frame i
            for jj in range (len(right_IDs[i])): # loop over right corners of frame j
                if left_IDs[i][ii]==right_IDs[i][jj]: # if IDs are the same for right and left we add the corners to our list of points for extrinsics calib.
                    frame_left_IDs.append(left_IDs[i][ii])
                    frame_right_IDs.append(right_IDs[i][jj])
                    frame_world_corners.append(allCorners[left_IDs[i][ii]]) #real world corners
                    frame_left_corners.append(left_Corners[i][ii]) # left corners
                    frame_right_corners.append(right_Corners[i][jj])  # right corners
                    continue
        left_IDs_sampled.append(np.array(frame_left_IDs))
        right_IDs_sampled.append(np.array(frame_right_IDs))
        left_corners_sampled.append(np.array(frame_left_corners, dtype=np.float32))
        right_corners_sampled.append(np.array(frame_right_corners, dtype=np.float32))
        world_corners_sampled.append(np.array( frame_world_corners,dtype=np.float32).squeeze())

        if left_img and right_img:# to plot images with corners and their IDs
            left=cv2.aruco.drawDetectedCornersCharuco(left_img[i].copy(),left_corners_sampled[i], left_IDs_sampled[i])
            right=cv2.aruco.drawDetectedCornersCharuco(right_img[i].copy(),right_corners_sampled[i], right_IDs_sampled[i])
            stack=np.hstack((left,right))
            cv2.imshow('left and right frame #' + str(i)+  'IDs',stack)

            cv2.waitKey(0)
            cv2.destroyAllWindows()
            cv2.waitKey(1)
    stereo_criteria = (cv2.TERM_CRITERIA_COUNT +cv2.TERM_CRITERIA_EPS, 1000, 1e-9)
    #perform the real calibration
    ret, M1, d1, M2, d2, R, T, E, F = cv2.stereoCalibrate(objectPoints=world_corners_sampled,
                                         imagePoints1=left_corners_sampled, 
                                         imagePoints2=right_corners_sampled,
                                         cameraMatrix1=left_camera_matrix, distCoeffs1=left_distortion_coefficients,
                                         cameraMatrix2=right_camera_matrix, distCoeffs2=right_distortion_coefficients,
                                         imageSize= left_imsize,
                                         criteria=stereo_criteria,
                                         flags=cv2.CALIB_FIX_INTRINSIC)
    print('stereo calibration error: ' ,ret)
    return ret,R,T,E,F

# define chAruco board
nx=18
ny=12
aruco_dictionary = cv2.aruco.getPredefinedDictionary(cv2.aruco.DICT_5X5_250)# this is the dictionary used to generate the boards
board = cv2.aruco.CharucoBoard((nx, ny),1, 0.75, aruco_dictionary) # this is the board used for calibration nx, by, squarelength, markerlength,dictionary)


# load left and right images
leftnames=sorted(glob('left*.png'))
rightnames=sorted(glob('right*.png'))
left_im=[]
right_im=[]
for im1, im2 in zip(leftnames, rightnames):
    left_im.append(cv2.imread(im1, 1))
    right_im.append(cv2.imread(im2, 1))

# perforrm lens distortion - left camera
left_Corners,left_IDs,left_imsize=read_chessboards(left_im,board)
left_ret, left_camera_matrix, left_distortion_coefficients, left_rotation_vectors, left_translation_vectors=calibrate_camera(left_Corners,left_IDs,left_imsize[::-1],board)
   
# perforrm lens distortion - rightcamera
right_Corners,right_IDs,right_imsize=read_chessboards(right_im,board)
right_ret, right_camera_matrix, right_distortion_coefficients, right_rotation_vectors, right_translation_vectors=calibrate_camera(right_Corners,right_IDs,right_imsize[::-1],board)
   
#perform stereo calibration
stereo_err,R,T,E,F=calibrate_stereo(board,left_Corners,left_IDs,left_imsize[::-1],left_camera_matrix, left_distortion_coefficients,right_Corners,right_IDs,right_imsize[::-1],right_camera_matrix, right_distortion_coefficients,left_im,right_im)


R_left,R_right,P_left,P_right,Q,left_ROI,right_ROI=cv2.stereoRectify(left_camera_matrix, left_distortion_coefficients,
                                right_camera_matrix, right_distortion_coefficients,
                                left_imsize[::-1], R, T)

leftmap = cv2.initUndistortRectifyMap(left_camera_matrix,left_distortion_coefficients, R_left,P_left,left_imsize[::-1],cv2.CV_32FC1)
rightmap = cv2.initUndistortRectifyMap(right_camera_matrix,right_distortion_coefficients, R_right,P_right,right_imsize[::-1],cv2.CV_32FC1)
count=0
#plot the raw images and the rectified ones:
for left, right in zip(left_im, right_im):

    left_rect= cv2.remap(left, leftmap[0], leftmap[1], cv2.INTER_LANCZOS4, cv2.BORDER_CONSTANT,0) 
 
    right_rect= cv2.remap(right, rightmap[0], rightmap[1], cv2.INTER_LANCZOS4, cv2.BORDER_CONSTANT,0) 
    plt.figure()
    plt.subplot(2,2,1)
    plt.imshow(left),plt.title(str(count) +' raw left')
    plt.subplot(2,2,2)
    plt.imshow(right),plt.title(str(count) +' raw right')
    plt.subplot(2,2,3)
    plt.imshow(left_rect),plt.title('rectified left')
    plt.subplot(2,2,4)
    plt.imshow(right_rect),plt.title('rectified right')

    plt.show()    
    count+=1     

   

   

Is there anything obvious that I am doing wrong ?
Thanks for any help !