Comment trouver la matrice de covariance d'un polygone?

Imaginez que vous ayez un polygone défini par un ensemble de coordonnées et son centre de masse est à . Vous pouvez traiter le polygone comme une distribution uniforme avec une limite polygonale. $(x_1,y_1)...(x_n,y_n)$$(0,0)$

Je recherche une méthode qui trouvera la matrice de covariance d'un polygone .

Je soupçonne que la matrice de covariance d'un polygone est étroitement liée au deuxième moment de l'aire , mais si elles sont équivalentes, je ne suis pas sûr. Les formules trouvées dans l'article wikipedia que j'ai lié semblent (une supposition ici, ce n'est pas particulièrement clair pour moi d'après l'article) pour faire référence à l'inertie rotationnelle autour des axes x, y et z plutôt que des axes principaux du polygone.

(Soit dit en passant, si quelqu'un peut me montrer comment calculer les axes principaux d'un polygone, cela me serait également utile)

Il est tentant de simplement effectuer une PCA sur les coordonnées , mais cela pose le problème que les coordonnées ne sont pas nécessairement réparties uniformément autour du polygone et ne sont donc pas représentatives de la densité du polygone. Un exemple extrême est le contour du Dakota du Nord, dont le polygone est défini par un grand nombre de points suivant la rivière Rouge, plus seulement deux autres points définissant le bord ouest de l'État.

Faisons d'abord une analyse.

Supposons que dans le polygone sa densité de probabilité soit la fonction proportionnelle Alors la constante de proportionnalité est l'inverse de l'intégrale de sur le polygone,$\mathcal{P}$$p(x,y).$$p$

$$\mu_{0,0}(\mathcal{P})=\iint_{\mathcal P} p(x,y) \mathrm{d}x\,\mathrm{d}y.$$

Le barycentre du polygone est le point des coordonnées moyennes, calculées comme leurs premiers moments. Le premier est

$$\mu_{1,0}(\mathcal{P})=\frac{1}{\mu_{0,0}(\mathcal{P})} \iint_{\mathcal P} x\,p(x,y)\mathrm{d}x\,\mathrm{d}y.$$

Le tenseur inertiel peut être représenté comme le réseau symétrique de seconds moments calculé après la translation du polygone pour mettre son barycentre à l'origine: c'est-à-dire la matrice des seconds moments centraux

$$\mu^\prime_{k,l}(\mathcal{P}) = \frac{1}{\mu_{0,0}(\mathcal{P})} \iint_{\mathcal P} \left(x - \mu_{1,0}(\mathcal{P})\right)^k\,\left(y - \mu_{0,1}(\mathcal{P})\right)^l\,p(x,y)\mathrm{d}x\,\mathrm{d}y$$

où varie de à à Le tenseur lui-même - alias matrice de covariance - est$(k,l)$$(2,0)$$(1,1)$$(0,2).$

$$I(\mathcal{P}) = \pmatrix{\mu^\prime_{2,0}(\mathcal{P}) & \mu^\prime_{1,1}(\mathcal{P}) \\ \mu^\prime_{1,1}(\mathcal{P}) & \mu^\prime_{0,2}(\mathcal{P})}.$$

Une PCA de donne les axes principaux de ce sont les vecteurs propres unitaires mis à l'échelle par leurs valeurs propres.$I(\mathcal{P})$$\mathcal{P}:$

---

Ensuite, découvrons comment faire les calculs. Parce que le polygone est présenté comme une séquence de sommets décrivant sa frontière orientée il est naturel d'invoquer$\partial\mathcal P,$

Théorème de Green: où est une forme unique définie dans un voisinage de et

$$\iint_{\mathcal{P}} \mathrm{d}\omega = \oint_{\partial\mathcal{P}}\omega$$ $\omega = M(x,y)\mathrm{d}x + N(x,y)\mathrm{d}y$$\mathcal{P}$ $$\mathrm{d}\omega = \left(\frac{\partial}{\partial x}N(x,y) - \frac{\partial}{\partial y}M(x,y)\right)\mathrm{d}x\,\mathrm{d}y.$$

Par exemple, avec et une densité constante ( c'est -à- dire uniforme) nous pouvons (par inspection) sélectionner l'un des nombreux telles que$\mathrm{d}\omega = x^k y^l \mathrm{d}x\mathrm{d}y$$p,$

$$\omega(x,y) = \frac{-1}{l+1}x^k y^{l+1}\mathrm{d}x.$$

Le fait est que l'intégrale de contour suit les segments de ligne déterminés par la séquence de sommets. Tout segment de ligne allant du sommet au sommet peut être paramétré par une variable réelle sous la forme$\mathbf{u}$$\mathbf{v}$$t$

$$t \to \mathbf{u} + t\mathbf{w}$$

où est la direction normale d'unité de àLes valeurs de vont donc de à Sous ce paramétrage, et sont des fonctions linéaires de et et sont des fonctions linéaires de Ainsi l'intégrande de l'intégrale de contour sur chaque bord devient une fonction polynomiale de qui est facilement évaluée pour les petits et$\mathbf{w} \propto \mathbf{v}-\mathbf{u}$$\mathbf{u}$$\mathbf{v}.$$t$$0$$|\mathbf{v}-\mathbf{u}|.$$x$$y$$t$$\mathrm{d}x$$\mathrm{d}y$$\mathrm{d}t.$$t,$$k$$l.$

---

La mise en œuvre de cette analyse est aussi simple que le codage de ses composants. Au niveau le plus bas, nous aurons besoin d'une fonction pour intégrer une forme polynomiale sur un segment de ligne. Les fonctions de niveau supérieur les agrégeront pour calculer les moments bruts et centraux pour obtenir le barycentre et le tenseur inertiel, et enfin nous pouvons opérer sur ce tenseur pour trouver les axes principaux (qui sont ses vecteurs propres à l'échelle). Le Rcode ci-dessous effectue ce travail. Il ne fait aucune prétention d'efficacité: il ne vise qu'à illustrer l'application pratique de l'analyse qui précède. Chaque fonction est simple et les conventions de dénomination sont parallèles à celles de l'analyse.

Le code comprend une procédure pour générer des polygones fermés valides, simplement connectés et non auto-entrecroisés (en déformant aléatoirement des points le long d'un cercle et en incluant le sommet de départ comme point final afin de créer une boucle fermée). Voici quelques instructions pour tracer le polygone, afficher ses sommets, jouxter le barycentre et tracer les axes principaux en rouge (le plus grand) et bleu (le plus petit), créant un système de coordonnées orienté positivement centré sur le polygone.

#
# 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")

Edit: Je n'ai pas remarqué que Whuber avait déjà répondu. Je vais laisser cela comme un exemple d'une autre approche (peut-être moins élégante) du problème.

La matrice de covariance

Soit un point aléatoire de la distribution uniforme sur un polygone en zone . La matrice de covariance est:$(X,Y)$$P$$A$

$$C = \begin{bmatrix} C_{XX} & C_{XY} \\ C_{XY} & C_{YY} \end{bmatrix}$$

où est la variance de , est la variance de et est la covariance entre et . Cela suppose une moyenne nulle, car le centre de masse du polygone est situé à l'origine. La distribution uniforme attribue une densité de probabilité constante à chaque point de , donc:$C_{XX} = E[X^2]$$X$$C_{YY} = E[Y^2]$$Y$$C_{XY} = E[XY]$$X$$Y$$\frac{1}{A}$$P$

$$C_{XX} = \frac{1}{A} \underset{P}{\iint} x^2 dV \quad C_{YY} = \frac{1}{A} \underset{P}{\iint} y^2 dV \quad C_{XY} = \frac{1}{A} \underset{P}{\iint} x y dV \tag{1}$$

Triangulation

Au lieu d'essayer de s'intégrer directement sur une région compliquée comme , nous pouvons simplifier le problème en partitionnant en sous-régions triangulaires:$P$$P$$n$

$$P = T_1 \cup \cdots \cup T_n$$

Dans votre exemple, un partitionnement possible ressemble à ceci:

Il existe différentes façons de produire une triangulation (voir ici ). Par exemple, vous pouvez calculer la triangulation de Delaunay des sommets, puis éliminer les arêtes qui se trouvent en dehors de (car elle peut être non convexe comme dans l'exemple).$P$

Les intégrales sur peuvent ensuite être divisées en sommes d'intégrales sur les triangles:$P$

$$C_{XX} = \frac{1}{A} \sum_{i=1}^n \underset{T_i}{\iint} x^2 dV \quad C_{YY} = \frac{1}{A} \sum_{i=1}^n \underset{T_i}{\iint} y^2 dV \quad C_{XY} = \frac{1}{A} \sum_{i=1}^n \underset{T_i}{\iint} x y dV \tag{2}$$

Un triangle a de belles frontières simples, donc ces intégrales sont plus faciles à évaluer.

Intégration sur des triangles

Il existe différentes façons d'intégrer des triangles. Dans ce cas, j'ai utilisé une astuce qui consiste à mapper un triangle sur le carré de l'unité. La transformation en coordonnées barycentriques pourrait être une meilleure option.

Voici des solutions pour les intégrales ci-dessus, pour un triangle arbitraire défini par des sommets . Laisser:$T$$(x_1,y_1), (x_2,y_2), (x_3,y_3)$

$$v_x = \left[ \begin{smallmatrix} x_1 \\ x_2 \\ x_3 \end{smallmatrix} \right] \quad v_y = \left[ \begin{smallmatrix} y_1 \\ y_2 \\ y_3 \end{smallmatrix} \right] \quad \vec{1} = \left[ \begin{smallmatrix} 1 \\ 1 \\ 1 \end{smallmatrix} \right] \quad L = \left[ \begin{smallmatrix} 1 & 0 & 0 \\ 1 & 1 & 0 \\ 1 & 1 & 1 \end{smallmatrix} \right]$$

Alors:

$$\underset{T}{\iint} x^2 dV = \frac{A}{6} \text{Tr}(v_x v_x^T L) \quad \underset{T}{\iint} y^2 dV = \frac{A}{6} \text{Tr}(v_y v_y^T L) \quad \underset{T}{\iint} x y dV = \frac{A}{12} (\vec{1}^T v_x v_y^T \vec{1} + v_x^T v_y) \tag{3}$$

Tout mettre ensemble

Soit et les coordonnées x / y des sommets pour chaque triangle , comme ci-dessus. Branchez dans pour chaque triangle, en notant que les termes de zone s'annulent. Cela donne la solution:$v_x^i$$v_y^i$$T_i$$(3)$$(2)$

$$C_{XX} = \frac{1}{6} \sum_{i=1}^n \text{Tr} \big( v_x^i (v_x^i)^T L \big) \quad C_{YY} = \frac{1}{6} \sum_{i=1}^n \text{Tr} \big( v_y^i (v_y^i)^T L \big) \quad C_{XY} = \frac{1}{12} \sum_{i=1}^n \big( \vec{1}^T v_x^i (v_y^i)^T \vec{1} + (v_x^i)^T v_y^i \big) \tag{4}$$

Axes principaux

Les axes principaux sont donnés par les vecteurs propres de la matrice de covariance , tout comme dans l'ACP. Contrairement à l'ACP, nous avons une expression analytique pour , plutôt que d'avoir à l'estimer à partir de points de données échantillonnés. Notez que les sommets eux-mêmes ne sont pas un échantillon représentatif de la distribution uniforme sur , donc on ne peut pas simplement prendre l'échantillon de matrice de covariance des sommets. Mais, * est * une fonction relativement simple des sommets, comme vu dans .$C$$C$$P$$C$$(4)$