Finden Sie die Drehung zwischen Punktmengen

Ich habe zwei Sätze ( sourcund target) von Punkten (x,y), die ich ausrichten möchte. Was ich bisher gemacht habe ist:

• Finden Sie den Schwerpunkt jedes Punktesatzes
• Verwenden Sie den Unterschied zwischen den Zentroiden-Übersetzungen, den Punkt in xundy

Ich möchte die beste Drehung (in Grad) finden, um die Punkte auszurichten .

Irgendeine Idee?

Der M-Code befindet sich unten (mit Darstellungen zur Visualisierung der Änderungen):

# Raw data
## Source data
sourc = matrix( 
     c(712,960,968,1200,360,644,84,360), # the data elements 
     nrow=2, byrow = TRUE)

## Target data
target = matrix( 
  c(744,996,980,1220,364,644,68,336), # the data elements 
  nrow=2, byrow = TRUE)

# Get the centroids
sCentroid <- c(mean(sourc[1,]), mean(sourc[2,])) # Source centroid
tCentroid <- c(mean(target[1,]), mean(target[2,])) # Target centroid

# Visualize the points
par(mfrow=c(2,2))
plot(sourc[1,], sourc[2,], col="green", pch=20, main="Raw Data",
     lwd=5, xlim=range(sourceX, targetX),
     ylim=range(sourceY, targetY))
points(target[1,], target[2,], col="red", pch=20, lwd=5)
points(sCentroid[1], sCentroid[2], col="green", pch=4, lwd=2)
points(tCentroid[1], tCentroid[2], col="red", pch=4, lwd=2)

# Find the translation
translation <- tCentroid - sCentroid
target[1,] <- target[1,] - translation[1]
target[2,] <- target[2,] - translation[2]

# Get the translated centroids
tCentroid <- c(mean(target[1,]), mean(target[2,])) # Target centroid

# Visualize the translation
plot(sourc[1,], sourc[2,], col="green", pch=20, main="After Translation",
     lwd=5, xlim=range(sourceX, targetX),
     ylim=range(sourceY, targetY))
points(target[1,], target[2,], col="red", pch=20, lwd=5)
points(sCentroid[1], sCentroid[2], col="green", pch=4, lwd=2)
points(tCentroid[1], tCentroid[2], col="red", pch=4, lwd=2)

Dies kann mit dem Kabsch-Algorithmus erfolgen . Der Algorithmus findet die beste Schätzung der kleinsten Quadrate für die Rotation von$RX-Y$ wo $R$ ist Rotationsmatrix, $X$ und $Y$ sind Ihre Ziel- und Quellmatrizen mit 2 Zeilen und n Spalten.

In [ 1 ] wird gezeigt, dass dieses Problem durch Singularwertzerlegung gelöst werden kann. Der Algorithmus ist wie folgt:

1. Zentrieren Sie die Datensätze so, dass ihre Schwerpunkte auf dem Ursprung liegen.
2. Berechnen Sie die "Kovarianz" -Matrix $C$=$XY^T$.
3. Erhalten Sie die Singularwertzerlegung von $C=UDV^T$.
4. Richtungseinstellung $d=sign(det(C))$.
5. Dann die optimale Rotation $R=V\left( \begin{array}{ccc} 1 & 0 \\ 0 & d \\ \end{array} \right)U^T$

Ich kenne keine Implementierung in R, deshalb habe ich unten eine kleine Funktion geschrieben.

Ihre ersten Punkte:

src <- matrix(c(712,960,968,1200,360,644,84,360), nrow=2, byrow=TRUE)
trg <- matrix(c(744,996,980,1220,364,644,68,336), nrow=2, byrow=TRUE)

Kabsch-Algorithmus in einer R-Funktion:

kabsch2d <- function(Y, X) {
  X   <- X-rowMeans(X)
  Y   <- Y-rowMeans(Y)
  C   <- X %*% t(Y)
  SVD <- svd(C)
  D   <- diag(c(1, sign(det(C))))
  t(SVD$v) %*% D %*% t(SVD$u)
}

Zentrieren Sie die Punkte:

src <- src-rowMeans(src)
trg <- trg-rowMeans(trg)

Rotation erhalten:

rot <- kabsch2d(src, trg)

Ergebnis (schwarz - Originalquelle, rot - ursprüngliches Ziel, grün - gedrehtes Ziel)

plot(t(src), col="black", pch=19)
points(t(trg), col="red", pch=19)
points(t(rot %*% trg), col="green", pch=19)

[1] http://www.math.pku.edu.cn/teachers/yaoy/Fall2011/arun.pdf

Ich habe dies mit einer iterativen Optimalsuche gemacht und 2 Versionen getestet.
Ich habe die ursprünglichen Arrays genommen und sie zentriert, indem ich diese Arrays cSRCund nannte cTAR. Dann habe ich eine Schleife mit Winkeln gemacht$\varphi$ zwischen $0$ und $2 \pi$ und für jeden Winkel berechnete ich das Fehlerkriterium unter Verwendung der Differenz zwischen den gedrehten $\small D=rot(\text{cSRC} ,\varphi)- \text{cTAR}$.

1. In Version 1) habe ich als Kriterium die Quadratsumme aller Einträge in genommen $\small D$ wie

   $$err_1 = \small \sum_{k=1}^4 \small((D_{k,1})^2+(D_{k,2})^2)$$ und der Winkel $\small \varphi$bei dem der minimale Fehler aufgetreten ist, entspricht der kabsch2dProzedur in der Antwort von @Karolis.

2. In Version 2) habe ich als Kriterium die Summe der absoluten Entfernungen genommen, dh die Summe

   $$err_2=\small \sum_{k=1}^4 \small\sqrt{(D_{k,1})^2+(D_{k,2})^2}$$ und bekam einen etwas anderen Drehwinkel $\small \varphi$ für den kleinsten Fehler.

Ich weiß nicht, welches Kriterium besser zu Ihren Bedürfnissen passt.

Hier sind einige Ergebnisse aus dem Protokoll.

$$\small \begin{array} {r|cc} & \text{version } 1 & \text{version } 2\\ \hline \varphi & -0.04895304& -0.05093647 \\ \text{rotation} & \begin{bmatrix} 0.99880204& -0.04893349\\ 0.04893349& 0.99880204\\ \end{bmatrix} & \begin{bmatrix} 0.99870302 & -0.05091444\\ 0.05091444 & 0.99870302\\ \end{bmatrix} \\ \text{distances} & \begin{bmatrix} -6.80077266 & -0.86209739\\ 2.79924551 & -9.33782500\\ -0.61309522 & 6.94156520\\ 4.61462237 & 3.25835719\\ \end{bmatrix} & \begin{bmatrix} -6.78017751& -0.37062404 \\ 3.35787307 & -9.36574874 \\ -1.16459115 & 6.95324527 \\ 4.58689559 & 2.78312752 \\ \end{bmatrix} \end{array}$$