A long long time ago I needed to come up with mass and angular momentums for some 2D polygons in a simple physics engine. Although I could have used arbitrary numbers I wanted the mass and momemtum to be as if the polygons were cut out of a thin sheet material. At the time, many online resources only had the solution for triangles or convex polygons. I wanted a more general solution.
If you want is the source code, you can just jump to PolyMath.js
Finding the mass and momemtum of these lamina was pretty cool because I was able to use some of the fancy calculus I learned in school. Some quick googling showed that this is a common problem given to undergrad students. This method is also extended to more general shapes involving Bezier curves and holes.DISCLAIMER: I am not a math or physics guy. None of this stuff is rigorous. There is plenty of handwaving. My goal was to come up with working code and a way to test to make sure my results were correct. This article assumes the reader is familiar with basic calculus and physics.
Finding the mass, center of mass and moment of inertia is very easy for a system of point masses. You can approximate a polygon with squares. The smaller the squares the better the approximation.
The mass of our polygon will be proportional to the area of the polygon. We can multiply the area by a constant if we want to make the overall polygon heavier or lighter. We'll call this constant our area density. From now on we will only talk about area since we can scale everything afterwards.
We will treat each of the little squares as a point mass. The position is the x and y of the pixel and the mass is the area of the square. \begin{align} A &= \sum_{i=1}^{N} a_i \\ \end{align}
The center of mass, [Cx, Cy], of N point masses is the weighted averages of the points: \begin{align} C_x &= \frac{1}{N} \sum_{i=1}^{N} a_i x_i \\ C_y &= \frac{1}{N} \sum_{i=1}^{N} a_i y_i \\ \end{align}
The polar moment of inertia, usually written as IO or just I, is a scalaer quantity and the rotational equivalent of mass. For a system of N point masses: \begin{align} I_O &= \sum_{i=1}^{N} a_i ( x_i^2 + y_i^2 ) \\ \end{align}
Here are the functions I used to implement our little squares method: PixelPolyMath.js The functions use the HTML5 Canvas ImageData to access the pixel data of an already filled polygon.
The more squares we use the more accurate our approximation. Anytime you see approximations getting better as a result of adding up smaller things, integrals are in order. However we have to implement the little squares method so that we have something to check our final analysis again. Our integrals are going to be over the area of the polygon. However our polygon data is just a list of points.
There is something called Green's Theorem that relates integrals over a bounded area to a line integrals around the perimeter of said area. Since our list of points describe the polygon's perimeter this is exactly what we need.
If we take P0 to be a point our polygon and P1 as the next point, we can describe the edge P0->P1 as a parametric line: \begin{align} x &= P0_x + (P1_x - P0_x) \, t \\ y &= P0_y + (P1_y - P0_y) \, t \\ \end{align} Where t is [0, 1]. The derivatives of our edge are also very easy: \begin{align} dx &= (P1_x - P0_x) \, dt \\ dy &= (P1_y - P0_y) \, dt \\ \end{align}
Using equations above we can describe our polygon as a list of lines equations. If our area's perimeter, C, is composed of multiple segments, C1, C2, ... CN. If C1 through CN are what is called piecewise-smooth, we can apply Green's Theorem to each segment, Ci, and add them up to the integral of the perimeter as a whole.
C1 through CN are piecewise-smooth if the end points are connected, and if each segment is differentiable along its entirety. Polygons meet both of these requirements so our list of edges are in fact piecewise-smooth. Note that that polygons at not differentiable at the vertexes since the derivative on either side is different. Here is Green's Theorem stated: \[ \oint_C P\, dx + Q\, dy = \iint_D \left( \frac{\partial Q}{\partial x} - \frac{\partial P}{\partial y} \right) dA \\ \] This says that we can express a double integral over an area(the dA) as a line integral provided we can find two functions P and Q where (dQ/dx - dP/dy) is equal to the integral we are originally trying to solve. Finding P and Q tends to be the difficult part. Once we have them we can just plug them into a TI-89 or Wolfram Alpha to solve.
We will work out the case of area step by step. Center of mass and polar moment of inertia are only sketched out. Area: \begin{align} A &= \iint_D 1 \, dA \end{align}
We have to choose the right P and Q equation to satisfy Green's Theorem: \begin{align} A &= \iint_D 1 \, dA = \iint_D \left( \frac{\partial Q}{\partial x} - \frac{\partial P}{\partial y} \right) dA = \oint_C P\, dx + Q\, dy \\ 1 &= \frac{\partial Q}{\partial x} - \frac{\partial P}{\partial y} \\ P &= - \frac{1}{2} y \\ Q &= \frac{1}{2} x \\ A &= \oint_C P\, dx + Q\, dy = \oint_C -\frac{1}{2} y \, dx + \frac{1}{2} x \, dy \\ A &= \frac{1}{2} \oint_C x \, dy - y \, dx \end{align}
If we plug in the parametric equations for our edges we get the area for line segment: \begin{align} A &= \frac{1}{2} \oint_C x \, dy - y \, dx \\ &= \frac{1}{2} \left( \oint_C x \, dy - \oint_C y \, dx \right) \\ &= \frac{1}{2} \left( \int_0^1 (P0_x + (P1_x - P0_x)t)(P1_y - P0_y) \, dt - \int_0^1 (P0_y + (P1_y - P0_y)t)(P1_x - P0_x) \, dt \right) \\ &= \frac{1}{2} \left( P0_x P1_y - P1_x P0_y \right) \end{align}
Area Centroid: \begin{align} C_x &= \frac{1}{A} \iint_D x \, dA = \frac{1}{2A} \oint_C x^2 \, dy \\ C_y &= \frac{1}{A} \iint_D y \, dA = \frac{1}{2A} \oint_C y^2 \, dx \\ \end{align}
Polar Moment of Inertia: \begin{align} I &= \iint_D (x^2 + y^2) \, dA = \frac{1}{3} \left( \oint_C x^3 \, dy - \oint_C y^3 \, dx \right) \\ \end{align}
A more thorough explanation of the integrals, as well as the choice of P and Q can be found at here[PDF]. Below is a comparison of the analytical solution to the "pixel" little squares approximation. You can add points to the polygon by clicking on the edges. Below that is the final source code used to compute our numbers. The results are negated because I made the map assuming that the Y axis goes "up" but on HTML canvases, as well as most screens, the Y axis goes down.
|
var PolyMath = {
calcArea: function(poly) {
var area = 0
for(var i = 0; i < poly.points.length; i++) {
var j = (i + 1)%poly.points.length
var p0 = poly.points[i]
var p1 = poly.points[j]
area += 0.5*(p0.x*p1.y - p0.y*p1.x)
}
return -area
},
calcCentroid: function(poly) {
var area = this.calcArea(poly)
var c = new Vect3(0, 0, 0)
for(var i = 0; i < poly.points.length; i++) {
var j = (i + 1)%poly.points.length
var p0 = poly.points[i]
var p1 = poly.points[j]
c.x += ( (p1.y - p0.y)*(p1.x*p1.x + p1.x*p0.x + p0.x*p0.x) )/(6*area)
c.y += -( (p1.x - p0.x)*(p1.y*p1.y + p1.y*p0.y + p0.y*p0.y) )/(6*area)
}
return c.neg()
},
calcMomentOfInertia: function(poly, com) {
var ixx = 0
var iyy = 0
for(var i = 0; i < poly.points.length; i++) {
var j = (i + 1)%poly.points.length
var p0 = poly.points[i]
var p1 = poly.points[j]
var x0 = p0.x - com.x
var y0 = p0.y - com.y
var x1 = p1.x - com.x
var y1 = p1.y - com.y
ixx += +(y1 - y0)*(x1*x1 + x0*x0)*(x1 + x0)/12
iyy += -(x1 - x0)*(y1*y1 + y0*y0)*(y1 + y0)/12
}
return -(ixx + iyy)
}
}
We can support fancier shapes if we use the equations for Bezier curves. Below are the equations to plug into the integrals. The equations also hold for "holes" if we represent the holes as curves going in the opposite direction. TrueType Fonts originally only used line segments and quadratic Bezier curves. parametric equations for quadratic Bézier curves: \begin{align} x &= (1 - t)^2 P0_x + 2(1 - t)t P1_x + t^2 P2_x \\ y &= (1 - t)^2 P0_y + 2(1 - t)t P1_y + t^2 P2_y \\ dx &= 2(1 - t)(P1_x - P0_x) + 2t(P2_x - P1_x) \, dt \\ dy &= 2(1 - t)(P1_y - P0_y) + 2t(P2_y - P1_y) \, dt \\ \end{align}
COMING SOON
|