Wie finde ich die Kovarianzmatrix eines Polygons?

9

Stellen Sie sich vor, Sie haben ein Polygon, das durch einen Satz von Koordinaten und dessen Schwerpunkt bei . Sie können das Polygon als gleichmäßige Verteilung mit einer polygonalen Grenze behandeln. (x1,y1)...(xn,yn)(0,0)Geben Sie hier die Bildbeschreibung ein

Ich bin nach einer Methode, die die Kovarianzmatrix eines Polygons findet .

Ich vermute, dass die Kovarianzmatrix eines Polygons eng mit dem zweiten Moment der Fläche zusammenhängt , aber ob sie gleichwertig sind, weiß ich nicht. Die Formeln in dem von mir verlinkten Wikipedia-Artikel scheinen (eine Vermutung hier ist mir aus dem Artikel nicht besonders klar) sich eher auf die Rotationsträgheit um die x-, y- und z-Achse als auf die Hauptachsen des Polygons zu beziehen.

(Wenn mich übrigens jemand darauf hinweisen kann, wie man die Hauptachsen eines Polygons berechnet, wäre das auch für mich nützlich.)

Es ist verlockend, nur PCA an den Koordinaten durchzuführen , aber dabei tritt das Problem auf, dass die Koordinaten nicht unbedingt gleichmäßig über das Polygon verteilt sind und daher nicht für die Dichte des Polygons repräsentativ sind. Ein extremes Beispiel ist der Umriss von North Dakota, dessen Polygon durch eine große Anzahl von Punkten definiert wird, die dem Roten Fluss folgen, sowie nur zwei weitere Punkte, die den westlichen Rand des Staates definieren.

Ingolifs
quelle
Mit "finden" gehe ich davon aus, dass Sie einfach aus dem Polygon abtasten und dann die Kovarianz der Stichproben berechnen. Ist das nicht das, was Sie vorhaben?
Stephan Kolassa
Können Sie Ihren Beitrag auch so bearbeiten, dass er Koordinaten für Ihr Polygon enthält, damit die Leute damit herumspielen können?
Stephan Kolassa
1
@StephanKolassa Ich meine, das Polygon als einheitliche bivariate Wahrscheinlichkeitsdichte mit polygonaler Grenze zu behandeln. Sicher, Sie können Punkte abtasten und das Limit wäre dasselbe, aber ich suche nach einer A-priori-Methode. Das Bild ist nur eine Illustration von Farbe, die ich verwendet habe. Die realen Daten, die ich verwenden möchte, sind die Umrisse von Staaten und Regionen.
Ingolifs
1
Sie haben Recht, dass der übliche Begriff für "Kovarianzmatrix" Trägheitsmoment oder zweites Moment ist. Die Hauptachsen sind in ihren Ausrichtungen ausgerichtet. Das Ausführen von PCA auf den Koordinaten ist falsch: Es kommt der Annahme gleich, dass sich die gesamte Masse an den Eckpunkten befindet. Die direktesten Methoden zur Berechnung des Schwerpunkts - der erste Moment - werden in meinem Beitrag unter gis.stackexchange.com/a/22744/664 erläutert . Die zweiten Momente werden auf die gleiche Weise mit geringfügigen Änderungen berechnet. Auf der Kugel sind besondere Überlegungen erforderlich.
whuber
2
Es funktioniert umgekehrt: Berechnen Sie den Trägheitstensor und ermitteln Sie daraus seine Hauptachsen. Die Technik in Ihrem Fall beinhaltet den Satz von Green, der zeigt, dass die erforderlichen Integrale kann als Konturintegrale um der Einform berechnet werden, wobeiSolche Formen sind leicht zu finden, da jede geeignete lineare Kombination von und funktioniert. Das Konturintegral ist eine Summe von Integralen über den Kanten.
μk,l(P)=Pxkyldxdy
Pωdω=xkyldxdy.xkyl+1dxxk+1yldy
whuber

Antworten:

10

Lassen Sie uns zuerst eine Analyse durchführen.

Angenommen, innerhalb des Polygons seine Wahrscheinlichkeitsdichte die proportionale Funktion Dann ist die Proportionalitätskonstante die Umkehrung des Integrals von über dem Polygon.Pp(x,y).p

μ0,0(P)=Pp(x,y)dxdy.

Der Schwerpunkt des Polygons ist der Punkt der Durchschnittskoordinaten, der als erste Momente berechnet wird. Der erste ist

μ1,0(P)=1μ0,0(P)Pxp(x,y)dxdy.

Der Trägheitstensor kann als die symmetrische Anordnung von zweiten Momenten dargestellt werden, die nach der Übersetzung des Polygons berechnet wird, um seinen Schwerpunkt auf den Ursprung zu setzen, dh die Matrix der zentralen zweiten Momente

μk,l(P)=1μ0,0(P)P(xμ1,0(P))k(yμ0,1(P))lp(x,y)dxdy

wobei von bis bis reichen Der Tensor selbst - auch bekannt als Kovarianzmatrix - ist(k,l)(2,0)(1,1)(0,2).

I(P)=(μ2,0(P)μ1,1(P)μ1,1(P)μ0,2(P)).

Eine PCA von ergibt die Hauptachsen von Dies sind die Einheitseigenvektoren, die durch ihre Eigenwerte skaliert werden.I(P)P:


Als nächstes wollen wir herausfinden, wie die Berechnungen durchgeführt werden. Da das Polygon als eine Folge von Eckpunkten dargestellt wird, die seine orientierte Grenze es natürlich, es aufzurufenP,

Green's Theorem: wobei ist eine Einform, die in einer Nachbarschaft von und

Pdω=Pω
ω=M(x,y)dx+N(x,y)dyP
dω=(xN(x,y)yM(x,y))dxdy.

Zum Beispiel können wir mit und konstanter ( dh einheitlicher) Dichte (durch Inspektion) eine der vielen auswählen Lösungen wiedω=xkyldxdyp,

ω(x,y)=1l+1xkyl+1dx.

Der Punkt dabei ist, dass das Konturintegral den Liniensegmenten folgt, die durch die Folge von Eckpunkten bestimmt werden. Jedes Liniensegment von Vertex zu Vertex kann durch eine reelle Variable im Formular parametrisiert werdenuvt

tu+tw

Dabei ist die Einheitsnormalrichtung von nachDie Werte von reichen daher von bis Unter dieser Parametrisierung sind und lineare Funktionen von und und sind lineare Funktionen von Somit wird der Integrand des Konturintegrals über jeder Kante eine Polynomfunktion von die leicht für kleine und ausgewertet werden kannwvuuv.t0|vu|.xytdxdydt.t,kl.


Die Implementierung dieser Analyse ist so einfach wie die Codierung ihrer Komponenten. Auf der untersten Ebene benötigen wir eine Funktion, um eine Polynom-Einform über ein Liniensegment zu integrieren. Übergeordnete Funktionen aggregieren diese, um die rohen und zentralen Momente zu berechnen, um den Schwerpunkt und den Trägheitstensor zu erhalten, und schließlich können wir diesen Tensor bearbeiten, um die Hauptachsen (die seine skalierten Eigenvektoren sind) zu finden. Der folgende RCode führt diese Arbeit aus. Es erhebt keinen Anspruch auf Effizienz: Es soll nur die praktische Anwendung der vorstehenden Analyse veranschaulichen. Jede Funktion ist unkompliziert und die Namenskonventionen entsprechen denen der Analyse.

Der Code enthält eine Prozedur zum Generieren gültiger geschlossener, einfach verbundener, sich nicht selbst schneidender Polygone (durch zufälliges Verformen von Punkten entlang eines Kreises und Einbeziehen des Startscheitelpunkts als Endpunkt, um eine geschlossene Schleife zu erstellen). Im Folgenden finden Sie einige Anweisungen zum Zeichnen des Polygons, Anzeigen seiner Eckpunkte, Angrenzen an das Schwerpunktzentrum und Zeichnen der Hauptachsen in Rot (am größten) und Blau (am kleinsten), wodurch ein polygonzentriertes positiv orientiertes Koordinatensystem erstellt wird.

Abbildung mit Polygon- und Hauptachse

#
# Integrate a monomial one-form x^k*y^l*dx along the line segment given as an 
# origin, unit direction vector, and distance.
#
lintegrate <- function(k, l, origin, normal, distance) {
  # Binomial theorem expansion of (u + tw)^k
  expand <- function(k, u, w) {
    i <- seq_len(k+1)-1
    u^i * w^rev(i) * choose(k,i)
  }
  # Construction of the product of two polynomials times a constant.
  omega <- normal[1] * convolve(rev(expand(k, origin[1], normal[1])), 
                                expand(l, origin[2], normal[2]),
                                type="open")
  # Integrate the resulting polynomial from 0 to `distance`.
  sum(omega * distance^seq_along(omega) / seq_along(omega))
}
#
# Integrate monomials along a piecewise linear path given as a sequence of
# (x,y) vertices.
#
cintegrate <- function(xy, k, l) {
  n <- dim(xy)[1]-1 # Number of edges
  sum(sapply(1:n, function(i) {
    dv <- xy[i+1,] - xy[i,]               # The direction vector
    lambda <- sum(dv * dv)
    if (isTRUE(all.equal(lambda, 0.0))) {
      0.0
    } else {
      lambda <- sqrt(lambda)              # Length of the direction vector
      -lintegrate(k, l+1, xy[i,], dv/lambda, lambda) / (l+1)
    }
  }))
}
#
# Compute moments of inertia.
#
inertia <- function(xy) {
  mass <- cintegrate(xy, 0, 0)
  barycenter = c(cintegrate(xy, 1, 0), cintegrate(xy, 0, 1)) / mass
  uv <- t(t(xy) - barycenter)   # Recenter the polygon to obtain central moments
  i <- matrix(0.0, 2, 2)
  i[1,1] <- cintegrate(uv, 2, 0)
  i[1,2] <- i[2,1] <- cintegrate(uv, 1, 1)
  i[2,2] <- cintegrate(uv, 0, 2)
  list(Mass=mass,
       Barycenter=barycenter,
       Inertia=i / mass)
}
#
# Find principal axes of an inertial tensor.
#
principal.axes <- function(i.xy) {
  obj <- eigen(i.xy)
  t(t(obj$vectors) * obj$values)
}
#
# Construct a polygon.
#
circle <- t(sapply(seq(0, 2*pi, length.out=11), function(a) c(cos(a), sin(a))))
set.seed(17)
radii <- (1 + rgamma(dim(circle)[1]-1, 3, 3))
radii <- c(radii, radii[1])  # Closes the loop
xy <- circle * radii
#
# Compute principal axes.
#
i.xy <- inertia(xy)
axes <- principal.axes(i.xy$Inertia)
sign <- sign(det(axes))
#
# Plot barycenter and principal axes.
#
plot(xy, bty="n", xaxt="n", yaxt="n", asp=1, xlab="x", ylab="y",
     main="A random polygon\nand its principal axes", cex.main=0.75)
polygon(xy, col="#e0e0e080")
arrows(rep(i.xy$Barycenter[1], 2), 
       rep(i.xy$Barycenter[2], 2),
       -axes[1,] + i.xy$Barycenter[1],     # The -signs make the first axis .. 
       -axes[2,]*sign + i.xy$Barycenter[2],# .. point to the right or down.
       length=0.1, angle=15, col=c("#e02020", "#4040c0"), lwd=2)
points(matrix(i.xy$Barycenter, 1, 2), pch=21, bg="#404040")
whuber
quelle
+1 Wow, das ist eine großartige Antwort!
Amöbe
7

Edit: Habe nicht bemerkt, dass whuber schon geantwortet hat. Ich werde dies als Beispiel für eine andere (vielleicht weniger elegante) Herangehensweise an das Problem belassen.

Die Kovarianzmatrix

Let auf einem Polygon eine beliebige Stelle aus der Gleichverteilung mit Bereich . Die Kovarianzmatrix lautet:(X,Y)PA

C=[CXXCXYCXYCYY]

wobei die Varianz von , die Varianz von ist und die Kovarianz zwischen und . Dies setzt einen Mittelwert von Null voraus, da sich der Massenschwerpunkt des Polygons am Ursprung befindet. Die Gleichverteilung weist jedem Punkt in eine konstante Wahrscheinlichkeitsdichte zu , also:CXX=E[X2]XCYY=E[Y2]YCXY=E[XY]XY1AP

(1)CXX=1APx2dVCYY=1APy2dVCXY=1APxydV

Triangulation

Anstatt zu versuchen, eine komplizierte Region wie direkt zu integrieren , können wir das Problem vereinfachen, indem wir in dreieckige Unterregionen unterteilen:PPn

P=T1Tn

In Ihrem Beispiel sieht eine mögliche Partitionierung folgendermaßen aus:

Geben Sie hier die Bildbeschreibung ein

Es gibt verschiedene Möglichkeiten, eine Triangulation zu erzeugen (siehe hier ). Sie können beispielsweise die Delaunay-Triangulation der Scheitelpunkte berechnen und dann Kanten verwerfen, die außerhalb von (da diese wie im Beispiel möglicherweise nicht konvex sind ).P

Integrale über können dann in Summen von Integralen über den Dreiecken aufgeteilt werden:P

(2)CXX=1Ai=1nTix2dVCYY=1Ai=1nTiy2dVCXY=1Ai=1nTixydV

Ein Dreieck hat schöne, einfache Grenzen, sodass diese Integrale leichter zu bewerten sind.

Über Dreiecke integrieren

Es gibt verschiedene Möglichkeiten, über Dreiecke zu integrieren. In diesem Fall habe ich einen Trick verwendet , bei dem ein Dreieck dem Einheitsquadrat zugeordnet wird. Die Umwandlung in baryzentrische Koordinaten könnte eine bessere Option sein.

Hier sind Lösungen für die obigen Integrale für ein beliebiges Dreieck das durch Eckpunkte . Lassen:T(x1,y1),(x2,y2),(x3,y3)

vx=[x1x2x3]vy=[y1y2y3]1=[111]L=[100110111]

Dann:

(3)Tx2dV=A6Tr(vxvxTL)Ty2dV=A6Tr(vyvyTL)TxydV=A12(1TvxvyT1+vxTvy)

Alles zusammenfügen

Lassen und die x / y - Koordinaten der Eckpunkte für jedes Dreieck enthalten , wie oben. Stecken Sie in für jedes Dreieck und beachten Sie, dass sich die Flächenbegriffe aufheben. Dies ergibt die Lösung:vxivyiTi(3)(2)

(4)CXX=16i=1nTr(vxi(vxi)TL)CYY=16i=1nTr(vyi(vyi)TL)CXY=112i=1n(1Tvxi(vyi)T1+(vxi)Tvyi)

Hauptachsen

Die Hauptachsen sind wie bei PCA durch die Eigenvektoren der Kovarianzmatrix . Im Gegensatz zu PCA haben wir einen analytischen Ausdruck für , anstatt ihn aus abgetasteten Datenpunkten schätzen zu müssen. Beachten Sie, dass die Eckpunkte selbst keine repräsentative Stichprobe aus der gleichmäßigen Verteilung auf , so dass man nicht einfach die Stichproben-Kovarianzmatrix der Eckpunkte nehmen kann. Aber, * * ist eine relativ einfache Funktion der Eckpunkt, wie in zu sehen .CCPC(4)

user20160
quelle
2
+1 Dies kann vereinfacht werden, indem orientierte Dreiecke zugelassen werden, wodurch die Notwendigkeit einer ordnungsgemäßen Triangulation entfällt. Stattdessen können Sie einfach ein beliebiges Zentrum und die (vorzeichenbehafteten) Werte über den Dreiecken wird es oft gemacht, weil es viel weniger pingelig ist. Es ist leicht zu erkennen, dass eine solche Summierung im Wesentlichen dasselbe ist wie die Anwendung des Greenschen Theorems, da jeder Term in der Summation letztendlich eine Funktion der KanteDieser Ansatz wird im Abschnitt "Bereich" unter quantdec.com/SYSEN597/GTKAV/section2/chapter_11.htm veranschaulicht . OOPiPi+1:PiPi+1.
whuber
@whuber Interessant, danke für den Hinweis
user20160
Beide Antworten sind gut, wenn auch etwas über meinem Bildungsniveau. Sobald ich sicher bin, dass ich sie vollständig verstehe, werde ich versuchen herauszufinden, wer das Kopfgeld bekommt.
Ingolifs