Calculate moment of inertia of Koch snowflake
For concreteness sake, let's suppose we are computing the moment of inertia of the Koch Snowflake $K$ of uniform density and mass 1 shown on the left below. This snowflake is centered at the origin and has diameter 3 - that is, the maximum distance between two points is three. From here, it would be easy to compute the moment of another Koch flake of different size and uniform density, as long as it's centered at the origin.
The moment of inertia of this flake $K$, as I understand it, is simply $$\iint\limits_{K} (x^2+y^2) \, dA$$ divided by the area of the snowflake. There are any number of tools that can compute this numerically. Using Mathematica, I came up with 0.813719. (I could provide code, if desired.)
This value can also be computed exactly using the theory of self-similar integration as described in Bob Strichartz paper Evaluating Integrals Using Self-similarity. Applying this, I computed an exact value of $9/11$ or $0.\overline{81}$. The remainder of this answer will describe that process, though it does require a fairly sophisticated knowledge of self-similarity and measure theory.
We will use the fact that the Koch Snowflake is self-similar - composed of 7 copies of itself as shown in the figure above on the right. The scaling factor for the middle piece is $1/\sqrt{3}$ and it's $1/3$ for the others. Let us suppose that the functions in the IFS mapping $K$ onto the individual pieces are $T_0,T_1,\ldots,T_6$, with $T_0$ mapping onto the central piece. The exact functions in the IFS are
\begin{align} T_0(x,y) &= \frac{1}{\sqrt{3}}R\left(\frac{\pi}{2}\right)\left(\begin{array}{c}{x\\y}\end{array}\right) \\ T_{i}(x,y) &= \frac{1}{3}\left(\begin{array}{c}{x\\y}\end{array}\right) + \left(\begin{array}{c}{\cos\left(\frac{\pi}{6} + i\frac{\pi}{3}\right)\\ \sin\left(\frac{\pi}{6} + i\frac{\pi}{3}\right)}\end{array}\right), \end{align} where $\varphi$ ranges from $\pi/6$ to $11\pi/6$ in steps of $\pi/3$ and $R(\pi/2)$ represents a rotation through the angle $\pi/2$.
As explained by Strichartz, the integral of a function $f$ defined on a self-similar set can be computed with respect to any self-similar measure $\mu$. Specializing our notation somewhat to deal with this specific problem, a measure $\mu$ on $\mathbb R^2$ is called self-similar with respect to a list of weights $p_0,p_1,\ldots,p_6$, if it satisfies $$\mu(A) = \sum_{i=0}^6 p_i \mu(T_i^{-1}(A)$$ for every $A\subset\mathbb R^2$. As a result, any integral with respect to $\mu$ will satisfy $$\int f d\mu = \sum_{i=0}^6 p_i \int f \circ T_i \, d\mu.$$ As explained in Strichartz' paper there is a unique measure $\mu$ that satisfies the self-similarity condition for any given list of probabilities $p_0,p_1,\ldots,p_6$. We need to choose the probabilities so that the measure $\mu$ is uniform. Given an IFS with scaling factors $r_0, r_1, \ldots, r_6$, a uniform measure can be constructed by choosing $p_i = r_i^s$, where $s$ satisfies Moran's equation $$\sum_{i=0}^6 r_i^s = 1.$$ This yields the probability list $p_0=1/3$ and $p_i = 1/9$ for $i=1,\ldots,6$.
Now, since the total mass of $K$ with respect to $\mu$ is 1, we certainly have that the integral of any constant is just that constant, i.e. $$\int c \, d\mu = c.$$ We can use the self-similarity of the integral to compute the integrals of $f_1(x,y)=x$ and $f_2(x,y)=y$ as follows. Written down as a list of functions that return ordered pairs, the IFS is
\begin{align} T_0(x,y) &= \left(-\frac{y}{\sqrt{3}},\frac{x}{\sqrt{3}}\right) \\ T_1(x,y) &= \left(\frac{x}{3}+\frac{\sqrt{3}}{2},\frac{y}{3}+\frac{1}{2}\right) \\ T_2(x,y) &= \left(\frac{x}{3},\frac{y}{3}+1\right) \\ T_3(x,y) &= \left(\frac{x}{3}-\frac{\sqrt{3}}{2},\frac{y}{3}+\frac{1}{2}\right) \\ T_4(x,y) &= \left(\frac{x}{3}-\frac{\sqrt{3}}{2},\frac{y}{3}-\frac{1}{2}\right) \\ T_5(x,y) &= \left(\frac{x}{3},\frac{y}{3}-1\right) \\ T_6(x,y) &= \left(\frac{x}{3}+\frac{\sqrt{3}}{2},\frac{y}{3}-\frac{1}{2}\right) \end{align}
Note that $f_1$ returns just the first component and $f_2$ returns the second. Extracting those first components, multiplying them by the terms in the probability list and adding them up, to apply the self-similar integration identity, we get $$\int x \, d\mu = \frac{2}{9}\int x \, d\mu - \frac{1}{3\sqrt{3}} \int y \, d\mu.$$ Similarly, $$\int y \, d\mu = \frac{1}{3\sqrt{3}} \int x \, d\mu + \frac{2}{9}\int y \, d\mu.$$ This leads to a pair of equations in the unknowns $$\int x \, d\mu \: \text{ and } \int y \, d\mu,$$ with a unique solution that both are zero. From the symmetry of those functions and the domain, this is absolutely correct. A similar procedure may be carried out on the second order terms. We find that \begin{align} \int x^2 \, d\mu &= \int (1/3 + 2x^2/27 +y^2/9) \, d\mu \\ \int xy \, d\mu &= -\int xy/27 \, d\mu \\ \int y^2 \, d\mu &= \int (1/3 + x^2/9 + 2y^2/27) \, d\mu. \end{align} Solving these for the integrals of $x^2$ and $y^2$ and adding, we obtain the desired result.
An answer based off Mark McClure's but involving less theory and invoking the parallel axis theorem.
If we assume that the moment of inertia is always proportional to the mass, and proportional to the square of a characteristic length scale (in this case the length of the diameter as in Mark's answer), then $$ I = \gamma Ml^2, $$ where $\gamma$ is a constant. Using the principle of superposition, $$ I = I_{\text{centre}} + 6I_{\text{edge}}, $$ where the centre snowflake is in blue and the edge snowflake's are red/brown. Now $$\displaystyle I_{\text{centre}} = \gamma \frac{M}{3}\left(\frac{l}{\sqrt{3}}\right)^2 = \frac{\gamma}{9}Ml^2 = \frac{I}{9}.$$ Using the parallel axis theorem, $\displaystyle I_{\text{edge}} = I_{\text{COM}} + Md^2$, where $I_{\text{COM}}$ is the moment of inertia through the centre of mass, and $$\displaystyle I_{\text{COM}} = \gamma \frac{M}{9}\left(\frac{l}{3}\right)^2 = \frac{\gamma}{81}Ml^2 = \frac{I}{81},$$ and $d = l/3$, so $$\displaystyle I_{\text{edge}} = \frac{I}{81}+ \frac{Ml^2}{9},$$
So \begin{align*} I &= \frac{I}{9} + 6\left(\frac{I}{81}+ \frac{Ml^2}{9}\right), \\ I & = \frac{15I}{81} + \frac{2Ml^2}{3}, \\ \frac{22I}{27} & = \frac{2Ml^2}{3}, \\ I & = \frac{Ml^2}{11} \end{align*}
If we let $l=3$ we recover Mark's answer.
Here's the approach I would take.
(1) Moments of inertia are additive, so the idea would be to do a summation over all iterations of triangle-adding.
(2) The moment of inertia of the central triangle is straightforward to compute.
(3) At the $i_{th}$ iteration, there are $a_i$ triangles.
(4) The center of these triangles are a computable distance $d_i$ away from the axis. You can use the parallel axis theorem to give the moment of inertia about the center of the central triangle in terms of $d_i$ and the moment of inertia of the triangles at that iteration.
(5) Sum over all $i$.