The area moments of a polygon can be computed by generalizations of the shoelace formula for area.
In particular, the first and second order moments are given by https://en.wikipedia.org/wiki/Polygon#Centroid and https://en.wikipedia.org/wiki/Second_moment_of_area#Any_polygon.
I was wondering if there are simplified forms when the polygon is described by a Freeman chain (so that the coordinates vary by at most one unit from vertex to vertex).
I would solve the problem as follows. First, I need an analytic expression for an arbitrary polygon. Say that you have $n$ elements of length $S_k$ and angles $\Phi_k$ Then you can express the polygon in the complex plane with the fractal linkage (Gnomon: from Pharoahs to fractals, Michael J. Gazal$\acute{\text{e}}$, Princeton University Press, 1999) , that is
$$z=\text{cum}\sum_{k=1}^n S_ke^{i\Phi_k}$$
The perimeter, area, and centroid are given by
$$ s=\int |\dot z| du\\ A = {\textstyle{1 \over 2}}\int {\Im\left\{ {{z^{_*}}\dot z} \right\}} \,du\\\ z\left( G \right) = \frac{{\int {z\Im\left\{ {{z^{_*}}\dot z} \right\}} \,du}}{{3A}}\ $$
Zwikker, C. (1968) [The Advanced Geometry of Plane Curves and Their Applications, Dover Press] gives the derivations for the following moments:
$$ \text{Polar inertia}\\ J_0=\frac{1}{4}\int zz^{_*} {\Im\left\{ {{z^{_*}}\dot z} \right\}} du\\ \text{Binet's inertia moments}\\ J_x=\frac{1}{4}\int x^2 {\Im\left\{ {{z^{_*}}\dot z} \right\}} du\\ J_y=\frac{1}{4}\int y^2 {\Im\left\{ {{z^{_*}}\dot z} \right\}} du\\ \text{Centrifigal inertia moment}\\ J_{xy}=\frac{1}{4}\int xy {\Im\left\{ {{z^{_*}}\dot z} \right\}} du\\ $$
Notice that $J_0=J_x+J_y$.
At least, that's what I would do.