Calcul de la matrice d'homographie par la méthode des moindres carrés (méthode DLT)

introduction

En vision par ordinateur, une matrice d'homographie (matrice de conversion de projection) apparaît souvent. Il s'agit d'une matrice 3x3 qui peut projeter des coordonnées d'image de l'image originale vers l'image convertie lorsque la conversion de projection (agrandissement / réduction / rotation / mouvement parallèle, etc.) est effectuée sur une certaine image. A l'inverse, dans l'appariement d'images entre deux images, la matrice d'homographie peut être estimée lorsqu'il y a plusieurs points correspondants. En général, il y a une erreur dans les coordonnées des points correspondants, donc la matrice d'homographie qui minimise cette erreur doit être estimée par la méthode du carré minimum. Je peux.

Dans cet article, nous montrons la procédure d'estimation de la matrice d'homographie à l'aide de la méthode DLT à l'aide d'un exemple concret, et expliquons également en détail la théorie mathématique pour trouver la solution du carré minimum par décomposition de singularité (SVD), qui n'est pas bien expliquée. Je vais expliquer.

Matrice d'homographie

La matrice d'homographie et la transformation par projection des coordonnées de l'image sont exprimées comme suit. Dans la transformation d'homographie (transformation de projection), neuf composantes peuvent être déterminées librement, mais cette fois nous considérerons une transformation affine qui n'implique que la rotation et le mouvement parallèle.

H =  \lambda\left(
    \begin{array}{ccc}
      \cos \theta & -\sin\theta & t_1 \\
      \sin\theta & \cos\theta & t_2 \\
      0 & 0 & 1
    \end{array}
  \right)
\\
\vec{x}_1 = (x_1,y_1,1)^T,\ \vec{x}_2 = (x_2,y_2,1)^T \\
\vec{x}_2 = H\vec{x}_1

Ici, les coordonnées de l'image $ \ vec {x} _1, \ vec {x} _2 $ utilisent le système de coordonnées du même ordre. $ t_1 et t_2 $ représentent un mouvement parallèle, et le 2x2 supérieur gauche est la matrice de rotation. De plus, il existe une indéfinité multiple constante due à l'expansion et à la contraction, c'est-à-dire l'échelle, et $ \ lambda $ est un coefficient représentant cela. Il y a 9 composants dans la matrice d'homographie, mais le degré de liberté est de 8 en raison de l'indéfini d'un multiple constant. En d'autres termes, la matrice d'homographie peut être estimée s'il y a quatre points correspondants ou plus (deux équations x et y peuvent être établies pour chaque point correspondant).

exemple

Affinisez l'image d'origine pour générer deux images. La conversion Affin utilisait la fonction d'opencv. Il est possible d'estimer avec 4 points correspondants ou plus, mais cette fois nous préparerons 10 points. La transformation d'affine est effectuée par la matrice d'homographie pour obtenir les coordonnées. img1 src.png img2 affine_img.png

Coordonnées de point correspondantes

img1 (10, 3)
[[ 66 215   1]
 [135  32   1]
 [362 185   1]
 [479  40   1]
 [239 247   1]
 [407 147   1]
 [253  14   1]
 [ 11 103   1]
 [474 181   1]
 [ 57 222   1]]
img2 (10, 3)
[[ 57 229   1]
 [141  53   1]
 [354 225   1]
 [483  91   1]
 [226 276   1]
 [402 191   1]
 [260  45   1]
 [ 11 113   1]
 [466 231   1]
 [ 47 236   1]]

code


def main():

    os.chdir(PATH_NAME)
    src = cv2.imread("./stk.png ", cv2.IMREAD_COLOR)
    src = cv2.resize(src,(512,256))
    cv2.imwrite("src.png ", src)

    print(src.shape)

    #Rotation 5deg, mouvement parallèle(10,10)
    theta = np.radians(10)
    t1 = 10
    t2 = 10
    #Matrice d'homographie
    H = np.array([[np.cos(theta),-np.sin(theta),t1]
                     ,[np.sin(theta),np.cos(theta),t2]
                     ,[0,0,1]])

    print("Homography matrix=\n" ,H)


    #Coordonnées de l'image dans l'image d'origine*4 points
    img_1 = []
    for x in range(4):
        v = random.randint(0,src.shape[0])
        u = random.randint(0,src.shape[1])
        img_1.append([u,v,1])
    img_1 = np.array(img_1)
    print("img1",img_1.shape)
    print(img_1)

    #Coordonnées de l'image N points dans l'image convertie (conversion de projection)
    img_2 = []
    for row in img_1:
        x2 = H @ row
        x2 = x2.astype('int')
        x2[0] = x2[0]
        x2[1] = x2[1]
        img_2.append(x2)
    img_2 = np.array(img_2)
    print("img2",img_2.shape)
    print(img_2)

    #M est une conversion affine d'image matricielle 2x3
    M = H[:2,]
    affine_img = cv2.warpAffine(src, M, (src.shape[1], src.shape[0]))

    cv2.imshow("color", affine_img)
    cv2.waitKey(0)
    cv2.destroyAllWindows()
    cv2.imwrite("affine_img.png ", affine_img)

    return

Estimation de la matrice d'homographie (méthode DLT)

Maintenant, entrons dans le sujet principal à partir d'ici. La méthode DLT est utilisée pour l'estimation. Il s'agit d'une équation linéaire qui est transformée en une forme qui peut utiliser la méthode des moindres carrés. La méthode des moindres carrés minimise l'erreur carrée de $ A \ boldsymbol {x} - \ boldsymbol {b} = 0 $ dans une équation linéaire simultanée qui peut être écrite par $ A \ boldsymbol {x} = \ boldsymbol {b} $. Estimez $ \ boldsymbol {x} $. La formule de conversion de projection est

 \boldsymbol{x}_2 = H\boldsymbol{x}_1

était. C'est sous la forme $ A \ boldsymbol {x} = \ boldsymbol {b} $, mais ce que vous voulez estimer n'est pas $ x $, mais $ H $ dans la matrice de coefficients. Ici, en utilisant le fait que les vecteurs projetés sont parallèles, $ \ boldsymbol {x} _2 \ times H \ boldsymbol {x} _1 = 0 $. Lorsque la formule est transformée, elle devient comme suit.

\left[
    \begin{array}{ccccccccc}
      0 & 0 & 0 & -x_1 & -y_1 & -1 & y_2 x_1 & y_2 y_1 & y_2 \\
      x_1 & y_1 & 1 & 0 & 0 & 0 & -x_2 x_1 & -x_2 y_1 & -x_2 \\
    \end{array}
  \right]
\left[
    \begin{array}{c}
      h_1 \\
      h_2 \\
      h_3 \\
      h_4 \\
      h_5 \\
      h_6 \\
      h_7 \\
      h_8 \\
      h_9 \\
    \end{array}
  \right]
=\boldsymbol{0}

Ici, le groupe de points des points correspondants de img_1 et img_2 est exprimé comme suit.

\boldsymbol{x}_A=\{\boldsymbol{x}_1^1,\boldsymbol{x}_1^2,\cdots,\boldsymbol{x}_1^k\}\\
\boldsymbol{x}_B=\{\boldsymbol{x}_2^1,\boldsymbol{x}_2^2,\cdots,\boldsymbol{x}_2^k\}

Si la formule ci-dessus est étendue par les points correspondants,

\left[
    \begin{array}{ccccccccc}
      0 & 0 & 0 & -x_A^1 & -y_A^1 & -1 & y_B^1 x_A^1 & y_B^1 y_A^1 & y_B^1 \\
      x_A^1 & y_A^1 & 1 & 0 & 0 & 0 & -x_B^1 x_A^1 & -x_B^1 y_A^1 & -x_B^1 \\
\vdots \\
      0 & 0 & 0 & -x_A^k & -y_A^k & -1 & y_B^k x_A^k & y_B^k y_A^k & y_B^k \\
      x_A^k & y_A^k & 1 & 0 & 0 & 0 & -x_B^k x_A^k & -x_B^k y_A^k & -x_B^k \\
    \end{array}
  \right]
\left[
    \begin{array}{c}
      h_1 \\
      h_2 \\
      h_3 \\
      h_4 \\
      h_5 \\
      h_6 \\
      h_7 \\
      h_8 \\
      1 \\
    \end{array}
  \right]
=\boldsymbol{0}

Ici, puisqu'il y a indéfini d'un multiple constant, comme condition de contrainte pour réduire le degré de libertéh_9=1C'est dit. Matrice de coefficient 2k × 9 sur le côté gaucheAEt la composante de la matrice d'homographie exprimée en vecteur colonne\boldsymbol{h}PuisA\boldsymbol{h}=\boldsymbol{0}Ce sera. Vous êtes maintenant prêt à résoudre la méthode des moindres carrés.A\boldsymbol{x}=\boldsymbol{b}Correspond au cas où le côté droit de est 0 vecteur.A\boldsymbol{h}=\boldsymbol{0}Minimisez l'erreur de\boldsymbol{h}La matrice d'homographie est-elle requise par. La condition de contrainte est|\boldsymbol{h}|=0Dans certains cas.

Méthode du carré minimum

À propos, l'équation normale de la méthode des moindres carrés était la suivante.

A^{T}A\boldsymbol{x} = A^{T}\boldsymbol{b} 

Cette fois, le côté droit sera 0, donc

A^{T}A\boldsymbol{x} = 0

Ce sera. Cela correspond à cela.

Maintenant, Dans la méthode des moindres carrés

J = \frac{1}{2}  \sum_{i=1}^{m}(\sum_{j=1}^{n} a_{ij}x_{j}-b_i)^2 \rightarrow min

Avant de régler le différentiel partiel = 0 pour trouver la valeur minimale de

J = \frac{1}{2}( x^{T}A^{T}Ax-x^{T}A^{T}b-b^{T}Ax+b^{T}b )\\

C'était. Cette fois, c'est $ \ boldsymbol {b = 0} $, donc

J = \frac{1}{2} x^{T}A^{T}Ax

Ce sera. Ici, la matrice $ A ^ TA $ est connue pour être une matrice symétrique à valeur constante semi-régulière (matrice symétrique avec des valeurs propres de 0 ou positives), et à ce moment, elle est sous la forme quadratique $ \ boldsymbol {x} ^. Il existe un théorème selon lequel $ \ boldsymbol {x} $, qui minimise TA ^ TA \ boldsymbol {x} $, est le vecteur propre correspondant à la plus petite valeur propre de la matrice $ A ^ TA $. (→ Algèbre linéaire) Cette fois, $ \ boldsymbol {h} $ correspond à $ \ boldsymbol {x} $.

Décomposition de singularité

Le vecteur propre correspondant à la valeur propre minimale de la matrice $ A ^ TA $ peut être obtenu par décomposition de singularité (SVD). La décomposition en valeur singulière de $ A $

A = U\Sigma V^T \\
\Sigma = 
\left[
    \begin{array}{ccc}
  \sigma_1   &  & \\
      & \sigma_2 & \\
      & \cdots & \\
      & \ &  \sigma_r  \\
    \end{array}
  \right]

est. $ U $ sont les vecteurs propres de $ AA ^ T $ disposés par ordre décroissant de valeurs propres positives, et $ V $ sont les vecteurs propres de $ A ^ TA $ disposés par ordre décroissant de valeurs propres positives. $ \ Sigma_1, \ sigma_2, \ cdots, \ sigma_r $ sont appelés valeurs singulières et sont les racines carrées (> = 0) des valeurs propres de $ A ^ TA $. Ce que nous voulons ici, c'est le vecteur propre correspondant à la moindre valeur propre de $ A ^ TA $, qui est la solution des moindres carrés, donc la dernière ligne de $ V $ lui correspond. La valeur singulière peut être zéro. Cela signifie que $ A $ n'est pas un rang complet, mais même dans ce cas, la dernière ligne de $ V $, le vecteur propre correspondant à la singularité zéro, est la solution la moins carrée. En d'autres termes, si la valeur unique minimale est $ \ lambda_ {min} $,

A^TA\boldsymbol{h}=\lambda_{min} \boldsymbol{h}

La solution qui satisfait est le $ \ boldsymbol {h} $ souhaité. Divisez le tout pour que finalement $ h_ {9} $ devienne 1.

Expérience

L'estimation de $ H $ avec SVD, comme indiqué ci-dessous, montre une correspondance étroite. Si le score est faible, une solution sera obtenue s'il y en a 4 ou plus, mais ils ne correspondent pas très bien. Il semble y avoir diverses techniques telles que la prise du centre de gravité des coordonnées de l'image et leur normalisation afin de calculer précisément.

Homography matrix=
 [[ 0.9961947  -0.08715574 10.        ]
 [ 0.08715574  0.9961947  10.        ]
 [ 0.          0.          1.        ]]

Estimated Homography matrix=
 [[ 9.97456897e-01 -8.38633226e-02  8.92359681e+00]
 [ 8.68971782e-02  9.98203988e-01  9.29572691e+00]
 [ 7.52034139e-07  6.24553736e-06  1.00000000e+00]]

En utilisant cette estimation $ H $ pour transformer les coordonnées de l'image, c'est presque la même chose que l'utilisation du $ H $ original.

[382 105   1]
[382 105   1]

code:

    ##Estimation de la matrice H
    A = np.zeros((img_1.shape[0]*2,9))
    for i, (a, b) in enumerate(zip(img_1,img_2)):
        A[i * 2] = np.array([a[0], a[1], 1, 0, 0, 0, -b[0] * a[0], -b[0] * a[1], -b[0]])
        A[i * 2 + 1] = np.array([0, 0, 0, a[0], a[1], 1, -b[1] * a[0], -b[1] * a[1], -b[1]])
    print(A.shape)
    print("A=",A)

    print(H.flatten())
    print(np.cross(img_2[0],(H @ img_1[0])))
    print(A @ H.flatten())

    # singular value decomposition
    u, s, vh = svd(A)
    print('\nSVD result')
    print('shape of u, s, vh:', u.shape, s.shape, vh.shape)
    print('singular values:', s.round(2))
    min = 8
    print('minimum eigen vector:', vh[min])

    Hest = vh[min].reshape((3,3))
    Hest = Hest / Hest[2,2]
    print("Estimated Homography matrix=\n" ,Hest)

    print((H @ img_1[0] / (H @ img_1[0])[2]).astype('int'))
    print((Hest @ img_1[0] / (Hest @ img_1[0])[2]).astype('int'))

Recommended Posts

Calcul de la matrice d'homographie par la méthode des moindres carrés (méthode DLT)
Calcul de la matrice d'homographie par la méthode des moindres carrés (méthode DLT)
Analyse de régression simple par la méthode des moindres carrés
Comment calculer le jour Julius
[Python] Méthode de calcul avec numpy
Méthode du carré minimum et méthode d'estimation la plus probable (comparaison par ajustement du modèle)
Analyse de régression simple par la méthode des moindres carrés
Approximation par la méthode des moindres carrés d'un cercle à deux points fixes
Calcul numérique du fluide compressible par la méthode des volumes finis
Calcul de similitude par MinHash
[Calcul scientifique / technique par Python] Solution numérique du problème des valeurs propres de la matrice par multiplication de puissance, algèbre linéaire numérique
Visualisation de la matrice créée par numpy
Méthode du carré minimum et méthode d'estimation la plus probable (comparaison par ajustement du modèle)
[Calcul scientifique / technique par Python] Calcul du produit de la matrice par l'opérateur @, python3.5 ou supérieur, numpy
Résumé de la méthode de connexion par DB de SQL Alchemy
Trouvez le coefficient du polypole le moins carré
[Calcul scientifique / technique par Python] Calcul de matrice inverse, numpy
À propos de la précision de la méthode de calcul du rapport de circonférence d'Archimède
[Python] Calcul numérique d'une équation de diffusion thermique bidimensionnelle par méthode ADI (méthode implicite de direction alternative)