[PYTHON] 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
img2
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é
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