Как найти ковариационную матрицу многоугольника?

Представьте, что у вас есть многоугольник, определенный набором координат и его центр масс находится в . Вы можете рассматривать полигон как равномерное распределение с полигональной границей. $(x_1,y_1)...(x_n,y_n)$$(0,0)$

Мне нужен метод, который найдет ковариационную матрицу многоугольника .

Я подозреваю, что ковариационная матрица многоугольника тесно связана со вторым моментом площади , но я не уверен, эквивалентны ли они. Формулы, найденные в статье в Википедии, которую я связал, кажутся (наверное, из статьи мне это не особо понятно) относятся к инерции вращения вокруг осей x, y и z, а не к главным осям многоугольника.

(Кстати, если кто-нибудь может указать мне, как рассчитать главные оси многоугольника, это также будет полезно для меня)

Соблазнительно просто выполнить PCA по координатам , но при этом возникает проблема, заключающаяся в том, что координаты не обязательно равномерно распределены по многоугольнику и, следовательно, не отражают плотность многоугольника. Экстремальным примером является контур Северной Дакоты, полигон которого определяется большим количеством точек, следующих за Красной рекой, плюс только еще две точки, определяющие западный край штата.

Давайте сначала проведем некоторый анализ.

Предположим, что внутри многоугольника его плотность вероятности пропорциональна функции Тогда константа пропорциональности является обратной величиной от интеграла по многоугольнику,$\mathcal{P}$$p(x,y).$$p$

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

Барицентр многоугольника является точкой средних координат, вычисленных как их первых моментов. Первый из них

$$\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.$$

Тензор инерции может быть представлена в виде симметричной матрицы вторых моментов , вычисленных после перевода многоугольника поставить свой барицентр в начале координат: то есть матрица центральных моментов второго

$$\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$$

где составляет от до до Сам тензор - он же ковариационная матрица - это$(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})}.$$

РС из дает главную ось из эти единичные собственные векторы , масштабированные их собственных значений.$I(\mathcal{P})$П :$\mathcal{P}:$

---

Далее давайте разберемся, как делать расчеты. Поскольку многоугольник представлен в виде последовательности вершин, описывающих его ориентированную границу естественно вызывать$\partial\mathcal P,$

Теорема Грина: где - это единая форма, определенная в окрестности и

$$\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.$$

Например, с и постоянной ( то есть равномерной) плотностью мы можем (осмотром) выбрать один из множества решения, такие как$\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.$$

Дело в том, что контурный интеграл следует отрезкам, определяемым последовательностью вершин. Любой отрезок прямой от вершины до вершины может быть параметризован действительной переменной в форме$\mathbf{u}$$\mathbf{v}$$t$

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

где - это единичное нормальное направление от доСледовательно, значения варьируются от до При этой параметризации и являются линейными функциями от а и являются линейными функциями от Таким образом, подынтегральное выражение интеграла контура по каждому ребру становится полиномиальной функцией от которая легко вычисляется для малых и$\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.$

---

Реализация этого анализа так же проста, как и кодирование его компонентов. На нижнем уровне нам понадобится функция для интегрирования полиномиальной формы на отрезке прямой. Функции более высокого уровня агрегируют их для вычисления необработанных и центральных моментов, чтобы получить барицентр и инерционный тензор, и, наконец, мы можем оперировать этим тензором, чтобы найти главные оси (которые являются его масштабированными собственными векторами). RНиже код выполняет эту работу. Он не претендует на эффективность: он предназначен только для иллюстрации практического применения вышеизложенного анализа. Каждая функция проста, а соглашения об именах параллельны соглашениям анализа.

В код включена процедура для генерации допустимых замкнутых, односвязных, несамопересекающихся многоугольников (путем случайного деформирования точек вдоль окружности и включения начальной вершины в качестве конечной точки для создания замкнутого контура). Далее следуют несколько утверждений для построения многоугольника, отображения его вершин, примыкания к центру барицентра и построения главных осей в красном (наибольшем) и синем (наименьшем) виде, создавая ориентированную на многоугольник положительно ориентированную систему координат.

#
# 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 уже ответил. Я оставлю это как пример другого (возможно, менее изящного) подхода к проблеме.

Ковариационная матрица

Пусть случайная точка из равномерного распределения на многоугольник с площадью . Ковариационная матрица имеет вид:$(X,Y)$$P$$A$

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

где - дисперсия , - дисперсия , а - ковариация между и . Это предполагает нулевое среднее значение, так как центр масс многоугольника находится в начале координат. Равномерное распределение назначает постоянную плотность вероятности каждой точке в , поэтому:$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}$$

триангуляция

Вместо того, чтобы пытаться напрямую интегрироваться в сложную область, такую ​​как , мы можем упростить задачу, разбив на треугольных областей:$P$$P$$n$

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

В вашем примере одно из возможных разделов выглядит следующим образом:

Существуют различные способы создания триангуляции (см. Здесь ). Например, вы можете вычислить триангуляцию Делоне для вершин, а затем отбросить ребра, которые выходят за пределы (поскольку это может быть невыпуклым, как в примере).$P$

Интегралы по можно затем разбить на суммы интегралов по треугольникам:$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}$$

Треугольник имеет красивые простые границы, поэтому эти интегралы легче оценить.

Интегрирование по треугольникам

Существуют различные способы объединения по треугольникам. В этом случае я использовал трюк, который включает в себя отображение треугольника на единицу площади. Преобразование в барицентрические координаты может быть лучшим вариантом.

Ниже приведены решения для приведенных выше интегралов для произвольного треугольника определенного вершинами . Позволять:$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]$$

Затем:

$$\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}$$

Собираем все вместе

Пусть и содержат координаты x / y вершин для каждого треугольника , как указано выше. Вставьте в для каждого треугольника, отметив, что условия области отменяются. Это дает решение:$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}$$

Главные оси

Главные оси задаются собственными векторами ковариационной матрицы , как и в PCA. В отличие от PCA, у нас есть аналитическое выражение для , вместо того, чтобы оценивать его по выборочным точкам данных. Обратите внимание, что сами вершины не являются репрезентативной выборкой из равномерного распределения на , поэтому нельзя просто взять выборочную ковариационную матрицу вершин. Но * является * относительно простой функцией вершин, как видно из .$C$$C$$P$$C$$(4)$