# How to find the covariance matrix of a polygon?

Imagine you have a polygon defined by a set of coordinates $$(x_1,y_1)...(x_n,y_n)$$ and its centre of mass is at $$(0,0)$$. You can treat the polygon as a uniform distribution with a polygonal boundary.

I'm after a method that will find the covariance matrix of a polygon.

I suspect that the covariance matrix of a polygon is closely related to the second moment of area, but whether they are equivalent I'm not sure. The formulas found in the wikipedia article I linked seem (a guess here, it's not especially clear to me from the article) to refer to the rotational inertia around the x, y and z axes rather than the principal axes of the polygon.

(Incidentally, if anyone can point me to how to calculate the principal axes of a polygon, that would also be useful to me)

It is tempting to just perform PCA on the coordinates, but doing so runs into the issue that the coordinates are not necessarily evenly spread around the polygon, and are therefore not representative of the density of the polygon. An extreme example is the outline of North Dakota, whose polygon is defined by a large number of points following the Red river, plus only two more points defining the western edge of the state.

• By "find", I assume simply sampling from the polygon, then calculating the covariance of the samples, is not what you have in mind? – Stephan Kolassa Jul 4 '19 at 6:32
• Also, can you edit your post to include coordinates for your polygon, so people can play around with it? – Stephan Kolassa Jul 4 '19 at 6:33
• @StephanKolassa I mean treating the polygon as a uniform bivariate probability density with polygonal boundary. Sure, you can sample points and the limit would be the same thing, but i'm looking for an a-priori method. The picture is just an illustration from paint that I used. The real world data I intend to use are the outlines of states and regions. – Ingolifs Jul 4 '19 at 9:46
• You are correct that the usual term for "covariance matrix" is moment of inertia or second moment. The principal axes are oriented in its eigendirections. Running PCA on the coordinates is incorrect: it is tantamount to assuming all the mass is located at the vertices. The most direct methods of computation of the barycenter--the first moment--are discussed in my post at gis.stackexchange.com/a/22744/664. The second moments are computed in the same way with minor modifications. Special considerations are needed on the sphere. – whuber Jul 7 '19 at 13:21
• It works the other way: compute the inertial tensor and find its principal axes from that. The technique in your case involves Green's Theorem, which shows that the requisite integrals $$\mu_{k,l}(\mathcal{P})=\iint_{\mathcal{P}}x^ky^l\mathrm{d}x\mathrm{d}y$$ can be computed as contour integrals around $\partial\mathcal{P}$ of the one-form $\omega$ where $\mathrm{d}\omega=x^ky^l\mathrm{d}x\mathrm{d}y.$ Such forms are easy to find because any suitable linear combination of $x^ky^{l+1}\mathrm{d}x$ and $x^{k+1}y^l\mathrm{d}y$ will work. The contour integral is a sum of integrals over the edges. – whuber Jul 8 '19 at 20:22

Let's do some analysis first.

Suppose within the polygon $$\mathcal{P}$$ its probability density is proportional function $$p(x,y).$$ Then the constant of proportionality is the inverse of the integral of $$p$$ over the polygon,

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

The barycenter of the polygon is the point of average coordinates, computed as their first moments. The first one is

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

The inertial tensor can be represented as the symmetric array of second moments computed after translating the polygon to put its barycenter at the origin: that is, the matrix of central second moments

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

where $$(k,l)$$ range from $$(2,0)$$ to $$(1,1)$$ to $$(0,2).$$ The tensor itself--aka covariance matrix--is

$$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})}.$$

A PCA of $$I(\mathcal{P})$$ yields the principal axes of $$\mathcal{P}:$$ these are the unit eigenvectors scaled by their eigenvalues.

Next, let's work out how to do the calculations. Because the polygon is presented as a sequence of vertices describing its oriented boundary $$\partial\mathcal P,$$ it is natural to invoke

Green's Theorem: $$\iint_{\mathcal{P}} \mathrm{d}\omega = \oint_{\partial\mathcal{P}}\omega$$ where $$\omega = M(x,y)\mathrm{d}x + N(x,y)\mathrm{d}y$$ is a one-form defined in a neighborhood of $$\mathcal{P}$$ and $$\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.$$

For instance, with $$\mathrm{d}\omega = x^k y^l \mathrm{d}x\mathrm{d}y$$ and constant (i.e., uniform) density $$p,$$ we may (by inspection) select one of the many solutions, such as $$\omega(x,y) = \frac{-1}{l+1}x^k y^{l+1}\mathrm{d}x.$$

The point of this is that the contour integral follows the line segments determined by the sequence of vertices. Any line segment from vertex $$\mathbf{u}$$ to vertex $$\mathbf{v}$$ can be parameterized by a real variable $$t$$ in the form

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

where $$\mathbf{w} \propto \mathbf{v}-\mathbf{u}$$ is the unit normal direction from $$\mathbf{u}$$ to $$\mathbf{v}.$$ The values of $$t$$ therefore range from $$0$$ to $$|\mathbf{v}-\mathbf{u}|.$$ Under this parameterization $$x$$ and $$y$$ are linear functions of $$t$$ and $$\mathrm{d}x$$ and $$\mathrm{d}y$$ are linear functions of $$\mathrm{d}t.$$ Thus the integrand of the contour integral over each edge becomes a polynomial function of $$t,$$ which is easily evaluated for small $$k$$ and $$l.$$

Implementing this analysis is as straightforward as coding its components. At the lowest level we will need a function to integrate a polynomial one-form over a line segment. Higher level functions will aggregate these to compute the raw and central moments to obtain the barycenter and inertial tensor, and finally we can operate on that tensor to find the principal axes (which are its scaled eigenvectors). The R code below performs this work. It makes no pretensions of efficiency: it is intended only to illustrate the practical application of the foregoing analysis. Each function is straightforward and the naming conventions parallel those of the analysis.

Included in the code is a procedure to generate valid closed, simply connected, non-self-intersecting polygons (by randomly deforming points along a circle and including the starting vertex as its final point in order to create a closed loop). Following this are a few statements to plot the polygon, display its vertices, adjoin the barycenter, and plot the principal axes in red (largest) and blue (smallest), creating a polygon-centric positively-oriented coordinate system.

#
# 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))
#
# 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: Didn't notice that whuber had already answered. I'll leave this up as an example of another (perhaps less elegant) approach to the problem.

### The covariance matrix

Let $$(X,Y)$$ be a random point from the uniform distribution on a polygon $$P$$ with area $$A$$. The covariance matrix is:

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

where $$C_{XX} = E[X^2]$$ is the variance of $$X$$, $$C_{YY} = E[Y^2]$$ is the variance of $$Y$$, and $$C_{XY} = E[XY]$$ is the covariance between $$X$$ and $$Y$$. This assumes zero mean, since the polygon's center of mass is located at the origin. The uniform distribution assigns constant probability density $$\frac{1}{A}$$ to every point in $$P$$, so:

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

Instead of trying to directly integrate over a complicated region like $$P$$, we can simplify the problem by partitioning $$P$$ into $$n$$ triangular subregions:

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

In your example, one possible partitioning looks like this:

There are various ways to produce a triangulation (see here). For example, you could compute the Delaunay triangulation of the vertices, then discard edges that fall outside $$P$$ (since it may be nonconvex as in the example).

Integrals over $$P$$ can then be split into sums of integrals over the triangles:

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

A triangle has nice, simple boundaries so these integrals are easier to evaluate.

### Integrating over triangles

There are various ways to integrate over triangles. In this case, I used a trick that involves mapping a triangle to the unit square. Transforming to barycentric coordintes might be a better option.

Here are solutions for the integrals above, for an arbitrary triangle $$T$$ defined by vertices $$(x_1,y_1), (x_2,y_2), (x_3,y_3)$$. Let:

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

Then:

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

### Putting everything together

Let $$v_x^i$$ and $$v_y^i$$ contain the x/y coordinates of the vertices for each triangle $$T_i$$, as above. Plug $$(3)$$ into $$(2)$$ for each triangle, noting that the area terms cancel out. This gives the solution:

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

### Principal axes

The principal axes are given by the eigenvectors of the covariance matrix $$C$$, just as in PCA. Unlike PCA, we have an analytic expression for $$C$$, rather than having to estimate it from sampled data points. Note that the vertices themselves are not a representative sample from the uniform distribution on $$P$$, so one can't simply take the sample covariance matrix of the vertices. But, $$C$$ *is* a relatively simple function of the vertices, as seen in $$(4)$$.

• +1 This can be simplified by allowing oriented triangles, thereby eliminating the need for a proper triangulation. Instead, you can just establish an arbitrary center $O$ and sum the (signed) values over the triangles $OP_iP_{i+1}:$ this is how it's often done because it's much less fussy. It's easy to see that such a summation is essentially the same thing as applying Green's Theorem, because each term in the summation ultimately is a function of the edge $P_iP_{i+1}.$ This approach is illustrated in the "Area" section at quantdec.com/SYSEN597/GTKAV/section2/chapter_11.htm. – whuber Jul 9 '19 at 15:46
• @whuber Interesting, thanks for pointing this out – user20160 Jul 9 '19 at 16:00
• Both of these answers are good, albeit a bit over my education level. Once I'm sure I fully understand them I'll try to figure out who gets the bounty. – Ingolifs Jul 11 '19 at 2:21