public static void calcMoments(PointF[] contour, out double ixx, out double ixy, out double iyy)
{
// make poly circular / wrap-around
PointF[] contourWrapped = new PointF[contour.Length + 1];
Array.Copy(contour, contourWrapped, contour.Length);
contourWrapped[contour.Length] = contourWrapped.First();
contour = contourWrapped;
// calculate all moments
int n = contour.Length;
double areaSum = 0, sumMomentX = 0, sumMomentY = 0, sumMomentXY = 0, sumCentroidX = 0, sumCentroidY = 0;
for (int i = 1; i < n; i++)
{
double crossProd = contour[i - 1].X * contour[i].Y - contour[i].X * contour[i - 1].Y;
sumMomentX += (contour[i - 1].X * contour[i - 1].X + contour[i - 1].X * contour[i].X + contour[i].X * contour[i].X) * crossProd;
sumCentroidX += (contour[i - 1].X + contour[i].X) * crossProd;
sumMomentY += (contour[i - 1].Y * contour[i - 1].Y + contour[i - 1].Y * contour[i].Y + contour[i].Y * contour[i].Y) * crossProd;
sumCentroidY += (contour[i - 1].Y + contour[i].Y) * crossProd;
sumMomentXY += (2 * contour[i - 1].X * contour[i - 1].Y + contour[i - 1].X * contour[i].Y + contour[i].X * contour[i - 1].Y + 2 * contour[i].X * contour[i].Y) * crossProd;
areaSum += crossProd;
}
// centre and scale moments
double area = areaSum / 2, areaSq = area * area;
ixx = sumMomentX / 12 / area - sumCentroidX * sumCentroidX / 36 / areaSq;
ixy = sumMomentXY / 24 / area - sumCentroidX * sumCentroidY / 36 / areaSq;
iyy = sumMomentY / 12 / area - sumCentroidY * sumCentroidY / 36 / areaSq;
}