[PYTHON] Berechnung der Homographiematrix nach der Methode der kleinsten Quadrate (DLT-Methode)
Einführung
In der Bildverarbeitung erscheint häufig eine Homographiematrix (Projektionsumwandlungsmatrix). Dies ist eine 3x3-Matrix, die von den Bildkoordinaten des Originalbilds auf das konvertierte Bild projizieren kann, wenn eine Projektionskonvertierung (Vergrößerung / Verkleinerung / Drehung / Parallelbewegung usw.) für ein bestimmtes Bild durchgeführt wird. Umgekehrt kann bei der Bildanpassung zwischen zwei Bildern die Homographiematrix geschätzt werden, wenn mehrere entsprechende Punkte vorhanden sind. Im Allgemeinen liegt ein Fehler in den Koordinaten der entsprechenden Punkte vor. Daher sollte die Homografiematrix, die diesen Fehler minimiert, nach der Minimum-Square-Methode geschätzt werden. Ich kann.
In diesem Artikel zeigen wir das Verfahren zur Schätzung der Homographiematrix unter Verwendung der DLT-Methode anhand eines konkreten Beispiels und erläutern ausführlich die mathematische Theorie zum Finden der minimalen quadratischen Lösung durch Singularitätszerlegung (SVD), die nicht gut erklärt ist. Ich werde erklären.
Homographiematrix
Die Homographiematrix und die Projektionstransformation der Bildkoordinaten werden wie folgt ausgedrückt. Bei der Homographietransformation (Projektionstransformation) können neun Komponenten frei bestimmt werden. Diesmal werden wir jedoch eine affine Transformation betrachten, die nur Rotation und parallele Bewegung umfasst.
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
Hier verwenden die Bildkoordinaten $ \ vec {x} _1, \ vec {x} _2 $ das Koordinatensystem gleicher Ordnung. $ t_1 und t_2 $ stehen für eine parallele Bewegung, und oben links 2x2 ist die Rotationsmatrix. Zusätzlich gibt es eine konstante mehrfache Unbestimmtheit aufgrund von Expansion und Kontraktion, dh Skalierung, und $ \ lambda $ ist ein Koeffizient, der dies darstellt. Die Homographiematrix enthält 9 Komponenten, der Freiheitsgrad beträgt jedoch 8, da ein konstantes Vielfaches unbestimmt ist. Mit anderen Worten kann die Homographiematrix geschätzt werden, wenn es vier oder mehr entsprechende Punkte gibt (zwei Gleichungen x und y können für jeden entsprechenden Punkt aufgestellt werden).
Beispiel
Affinieren Sie das Originalbild, um zwei Bilder zu generieren. Die Affin-Konvertierung verwendete die Funktion von opencv. Es ist möglich, mit 4 oder mehr entsprechenden Punkten zu schätzen, aber dieses Mal werden wir 10 Punkte vorbereiten. Die Affintransformation wird von der Homographiematrix durchgeführt, um die Koordinaten zu erhalten.
img1
img2
Entsprechende Punktkoordinaten
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)
#Drehung 5 Grad, parallele Bewegung(10,10)
theta = np.radians(10)
t1 = 10
t2 = 10
#Homographiematrix
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)
#Bildkoordinaten im Originalbild*4 Punkte
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)
#Bildkoordinaten N Punkte im konvertierten Bild (Projektionskonvertierung)
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 ist eine affine Umwandlung eines 2 × 3-Matrixbildes
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
Schätzung der Homographiematrix (DLT-Methode)
Kommen wir nun von hier aus zum Hauptthema. Die DLT-Methode wird zur Schätzung verwendet. Dies ist eine lineare Gleichung, die in eine Form umgewandelt wird, die die Methode der kleinsten Quadrate verwenden kann. Die Methode der kleinsten Quadrate minimiert den quadratischen Fehler von $ A \ boldsymbol {x} - \ boldsymbol {b} = 0 $ in einer linearen Simultangleichung, die durch $ A \ boldsymbol {x} = \ boldsymbol {b} $ geschrieben werden kann. Schätzen Sie $ \ boldsymbol {x} $. Die Projektionsumrechnungsformel lautet
\boldsymbol{x}_2 = H\boldsymbol{x}_1
war. Es hat die Form von $ A \ boldsymbol {x} = \ boldsymbol {b} $, aber was Sie schätzen möchten, ist nicht $ x $, sondern $ H $ in der Koeffizientenmatrix. Unter Verwendung der Tatsache, dass die projizierten Vektoren parallel sind, ist $ \ boldsymbol {x} _2 \ times H \ boldsymbol {x} _1 = 0 $. Wenn die Formel transformiert wird, wird sie wie folgt.
\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}
Hier wird die Punktgruppe der entsprechenden Punkte von img_1 und img_2 wie folgt ausgedrückt.
\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\}
Wenn die obige Formel um die entsprechenden Punkte erweitert wird,
\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}
Hier gibt es, da es eine Unbestimmtheit eines konstanten Vielfachen gibt, als Randbedingung, um den Freiheitsgrad zu verringern
Minimum-Quadrat-Methode
Übrigens war die Normalgleichung der Methode der kleinsten Quadrate wie folgt.
A^{T}A\boldsymbol{x} = A^{T}\boldsymbol{b}
Dieses Mal ist die rechte Seite 0, also
A^{T}A\boldsymbol{x} = 0
Es wird sein. Es entspricht dem.
Jetzt, In der Methode der kleinsten Quadrate
J = \frac{1}{2} \sum_{i=1}^{m}(\sum_{j=1}^{n} a_{ij}x_{j}-b_i)^2 \rightarrow min
Vor dem Setzen von Partial Differential = 0, um den Minimalwert von zu finden
J = \frac{1}{2}( x^{T}A^{T}Ax-x^{T}A^{T}b-b^{T}Ax+b^{T}b )\\
Es war. Diesmal ist es $ \ boldsymbol {b = 0} $, also
J = \frac{1}{2} x^{T}A^{T}Ax
Es wird sein. Hier ist bekannt, dass die Matrix $ A ^ TA $ eine halbregelmäßige symmetrische Matrix mit konstantem Wert ist (symmetrische Matrix mit Eigenwerten von 0 oder positiv), und zu diesem Zeitpunkt liegt sie in der quadratischen Form $ \ boldsymbol {x} ^ vor. Es gibt einen Satz, dass $ \ boldsymbol {x} $, der TA ^ TA \ boldsymbol {x} $ minimiert, der Eigenvektor ist, der dem kleinsten Eigenwert der Matrix $ A ^ TA $ entspricht. (→ Lineare Algebra) Dieses Mal entspricht $ \ boldsymbol {h} $ $ \ boldsymbol {x} $.
Singularitätszerlegung
Der Eigenvektor, der dem minimalen Eigenwert der Matrix $ A ^ TA $ entspricht, kann durch Singularitätszerlegung (SVD) erhalten werden. Die Singularwertzerlegung von $ A $
A = U\Sigma V^T \\
\Sigma =
\left[
\begin{array}{ccc}
\sigma_1 & & \\
& \sigma_2 & \\
& \cdots & \\
& \ & \sigma_r \\
\end{array}
\right]
ist. $ U $ sind die Eigenvektoren von $ AA ^ T $, die in absteigender Reihenfolge positiver Eigenwerte angeordnet sind, und $ V $ sind die Eigenvektoren von $ A ^ TA $, die in absteigender Reihenfolge positiver Eigenwerte angeordnet sind. $ \ Sigma_1, \ sigma_2, \ cdots, \ sigma_r $ heißen Singularwerte und sind die Quadratwurzeln (> = 0) der Eigenwerte von $ A ^ TA $. Was wir hier wollen, ist der Eigenvektor, der dem kleinsten Eigenwert von $ A ^ TA $ entspricht, der die Lösung mit dem kleinsten Quadrat ist, also entspricht die letzte Zeile von $ V $ diesem. Der Singularwert kann Null sein. Dies bedeutet, dass $ A $ nicht den vollen Rang hat, aber selbst dann ist die letzte Zeile von $ V $, der Eigenvektor, der der Singularität von Null entspricht, die Lösung mit dem kleinsten Quadrat. Mit anderen Worten, wenn der minimale eindeutige Wert $ \ lambda_ {min} $ ist,
A^TA\boldsymbol{h}=\lambda_{min} \boldsymbol{h}
Die Lösung, die erfüllt, ist das gewünschte $ \ boldsymbol {h} $. Teilen Sie das Ganze so, dass schließlich $ h_ {9} $ 1 wird.
Experiment
Die Schätzung von $ H $ mit SVD, wie unten gezeigt, zeigt eine enge Übereinstimmung. Wenn die Punktzahl niedrig ist, wird eine Lösung erhalten, wenn 4 oder mehr vorhanden sind, diese jedoch nicht sehr gut übereinstimmen. Es scheint verschiedene Techniken zu geben, z. B. den Schwerpunkt der Bildkoordinaten zu nehmen und sie zu normalisieren, um genau zu berechnen.
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]]
Die Verwendung dieses geschätzten $ H $ zum Transformieren der Bildkoordinaten entspricht fast der Verwendung des ursprünglichen $ H $.
[382 105 1]
[382 105 1]
Code:
##H-Matrixschätzung
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