Une déformation de l'image est toujours present et cela est du a la lentille, et doit être traitée avant tout autre analyse, en particulier si l'on souhait estimer une position. La déformation provoque une détection incorrecte, et endommage la position des pixels, elle a un impact sur l'estimation de la position réelle. Un modèle de non-distorsion de l'objectif doit être appliqué et être le même pour les paramètres intrinsèques de la matrice de la caméra $\mathcal{A}$ (ne dépend pas de la vue de la scène). Ainsi, une fois estimé, il peut être réutilisé tant que la longueur focale est fixe (c'est-à-dire le zoom).

Modele de distortion

Je propose le modèle de distorsion de lentille suivant, qui est obtenu en ajoutant des coefficients de distorsion tangentielle au modèle de distorsion radiale, une distorsion optionnelle à prisme mince peut également être envisagée. La fonction ci-dessus $\Phi$ définit la cartographie polynomiale $2D\rightarrow2D$ de la correction de la distorsion.

$\left(\begin{matrix} x_u \\ y_u \end{matrix}\right) = \Phi \left(\begin{matrix} x_d \\ y_d \end{matrix}\right) = \left(\begin{matrix} x_d\,(1 + k_1 r^2 + k_2 r^4 + k_3 r^6) + 2 p_1 x_d y_d + p_2 (r^2 + 2 x_d^2) + s_1 r^2 + s_2 r^4 \\ y_d\,(1 + k_1 r^2 + k_2 r^4 + k_3 r^6) + 2 p_2 x_d y_d + p_1 (r^2 + 2 y_d^2) + s_3 r^2 + s_4 r^4 \end{matrix}\right)$

où $r^2 = x_d^2+y_d^2$, $(x_d,y_d)^t$ sont les coordonnées de pixels normalisées déformées, $(x_u,y_u)^t$ sont les coordonnées d'image normalisées non déformées, $r^2=x_d^2+y_d^2$. Les coordonnées normalisées de l'image sont calculées à partir des coordonnées des pixels en les traduisant au centre optique et en les divisant par la distance focale en pixels. Ainsi, $(x_u,y_u)^t$ et $(x_d,y_d)^t$ sont sans dimension. Nous pouvons réécrire l'équation pour séparer la partie radiale et tangentielle de l'équation :

$$\begin{array}{rlcccccc}x_u &= x_d &+& {x_d(k_1 r^2 + k_2 r^4 + k_3 r^6)} &+& {2 p_1 x_d y_d + p_2(r^2 + 2 x_d^2)} &+& {s_1 r^2 + s_2 r^4} \\  y_u &= y_d &+& {\underbrace{y_d(k_1 r^2 + k_2 r^4 + k_3 r^6)}} &+& {\underbrace{2 p_2 x_d y_d + p_1(r^2 + 2 y_d^2)}} &+& {\underbrace{s_3 r^2 + s_4 r^4}} \\ & & & {\Delta {radial}} &+& {\Delta {tangential}} &+& {\Delta {prism}} \end{array}$$

Chaque partie de l'équation est expliquée ci-dessus :

• Distorsion radiale : La distorsion radiale se produit lorsque les rayons lumineux se courbent plus près des bords d'une lentille qu'ils ne le font en son centre optique. Plus la lentille est petite, plus la distorsion est importante. Les coefficients de distorsion radiale modélisent ce type de distorsion comme une transformation 2D$\rightarrow$2D. Généralement, deux coefficients de distorsion radiale $k_1$,$k_2$ sont suffisants pour l'étalonnage. Pour les distorsions importantes, comme dans les objectifs grand angle, nous pouvons inclure k3. Les coefficients de distorsion d'ordre supérieur tels que $k_4$, $k_5$ et $k_6$ ne sont généralement pas pris en compte.
• Distorsion tangentielle : La distorsion tangentielle corrige les inclinaisons du plan de l'image après une distorsion radiale. Les coefficients de distorsion tangentielle modélisent ce type de distorsion sous la forme d'une transformation 2D$\rightarrow$2D.
• Distorsion du prisme : La distorsion à prisme mince est due à une légère inclinaison de la lentille ou du réseau de capteurs d'images et provoque également une distorsion radiale et tangentielle supplémentaire. Ce type de déformation n'est généralement pas pris en compte, mais nous l'avons également examiné.

Model d'estimation

class CameraModel(nn.Module):

    def __init__(self, e=0.1):
        super(CameraModel, self).__init__()
        self.K = torch.nn.Parameter(torch.rand(3, dtype=torch.float32)*e)
        self.P = torch.nn.Parameter(torch.rand(2, dtype=torch.float32)*e)
        self.S = torch.nn.Parameter(torch.rand(4, dtype=torch.float32)*e)
        self.fc = torch.nn.Linear(4*2, 2)

    def radial(self, r, xy):
        R = torch.stack([r, r**2, r**3], dim=1)
        K = self.K.repeat(len(r), 1)
        return xy * torch.sum(K * R, dim=1).unsqueeze(1)

    def tangential(self, r, xy):
        P = self.P.unsqueeze(0).repeat(len(xy), 1)
        px = torch.prod(xy, dim=1).unsqueeze(1)
        t1 = 2*P*px
        t2 = P.flip(1)*(r.unsqueeze(1)+xy**2)
        return t1 + t2

    def prismatic(self, r, xy):
        s1, s2, s3, s4 = self.S
        S1 = self.S[0] * r + self.S[1] * r**2
        S2 = self.S[2] * r + self.S[3] * r**2
        P = torch.stack([S1,S2], dim=1)
        return P

    def forward(self, xy):
        r = torch.sum(xy**2, dim=1)
        R = self.radial(r, xy).to(device=xy.device)
        T = self.tangential(r, xy).to(device=xy.device)
        P = self.prismatic(r, xy).to(device=xy.device)
        return xy + R + T + P

    def forward_numpy(self, mapx, mapy):
        ax = mapx.reshape(np.prod(mapx.shape))
        ay = mapy.reshape(np.prod(mapy.shape))

        xy = np.vstack([ax, ay]).T
        xy = torch.tensor(xy, dtype=torch.float32).cuda()
        xy = self.forward(xy)
        xy = xy.cpu().detach().numpy()

        ax = xy[:,0].reshape(mapx.shape)
        ay = xy[:,1].reshape(mapy.shape)
        return ax, ay

Estimation OpenCV

Pour annuler une image déformée, il faut évaluer le mappage inverse $\Phi^{-1}$, mais il n'existe pas d'expression de forme fermée pour ce mappage inverse. En pratique, l'évaluation est effectuée par un processus de table de recherche. Dans l'ensemble, ce modèle de caméra avec distorsion de l'objectif comprend jusqu'à onze paramètres intrinsèques en plus des paramètres extrinsèques, ou pose, avec six degrés de liberté, décrits par la matrice de rotation $R=[r_1 r_2 r_3]$ et le vecteur de translation $t$.

import numpy as np
import cv2

chessboard_shape = (9,6)
criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 30, 0.001)
objp = np.zeros((chessboard_shape[0] * chessboard_shape[1], 3), np.float32)
objp[:,:2] = np.mgrid[0:chessboard_shape[0],0:chessboard_shape[1]].T.reshape(-1,2)

images = glob.glob('../data/calibration/chessboard/*.jpg')
objpoints = [] # 3d point in real world space
imgpoints = [] # 2d points in image plane.

print('finding chesboard parameters : ')

for i in tqdm(images):
    image = cv2.imread(i)
    gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
    
    ret, corners = cv2.findChessboardCorners(cv2.cvtColor(image,cv2.COLOR_BGR2GRAY), chessboard_shape, None)
 
    if ret == True:
        objpoints.append(objp)
        corners2 = cv2.cornerSubPix(gray, corners, (11,11), (-1,-1), criteria)
        imgpoints.append(corners2)

print('computing calibration ... ')

h, w = image.shape[:2]
calibration = cv2.calibrateCamera(
    objpoints, imgpoints, gray.shape[::-1], None, None,
    flags = cv2.CALIB_FIX_INTRINSIC + cv2.CALIB_THIN_PRISM_MODEL
)

ret, mtx, dist, rvecs, tvecs = calibration
cameramtx, roi = cv2.getOptimalNewCameraMatrix(mtx, dist, (w,h), 1, (w,h))

np.save('../result/mtx.npy', mtx)
np.save('../result/dist.npy', dist)
np.save('../result/cameramtx.npy', cameramtx)

Resulstat

La calibration est estimée avec la calibration de la caméra opencv en utilisant 77 images d'échiquier à différents angles et distances du capteur. Nous pouvons trouver les chiffres de correction $\Delta$ appliqués ci-dessous, de radial, tangentiel, prisme et les deux échantillonnés à chaque 30 pixel.