Moment of Inertia of an arbitrary shape represented by a finite number of spheres

Quick Background

I’m writing software for my research (primarily focused around orbit determination and gravity field estimation around irregular bodies). For generating the truth reference models, I am using shape models of the asteroid in a polyhedron format, which I then pack full of spheres (or cubes, but to keep this problem description simple, I am going to limit my description to just talking about the spheres). This allows me to calculate a gravity field for irregular bodies very easily by simply considering it as a collection of packed spheres, and summing their gravitational potentials together.

Question

Using this collection of sphere though, I would also like to generate a moment of inertia about the center of mass, which will allow me to accurately model the rotation of the asteroid. Each sphere in the model has its own unique mass and radius (they are NOT necessarily all identical, but they WILL all be complete; so no partial spheres/cubes). For a simple representation of this, just so things are clear, please refer to the last section.

My first approach was to simply consider the moment of inertia of a sphere:

$$I = left(frac{2}{5}right)m r^2$$

and then apply the parallel axis theorem to all of the spheres like this:

$$I_{xx} = sum_{i}^{N} left(frac{2}{5}right) m_i r_i^2 + m_i(y_i^2 + z_i^2)$$
$$I_{yy} = sum_{i}^{N} left(frac{2}{5}right) m_i r_i^2 + m_i(x_i^2 + z_i^2)$$
$$I_{zz} = sum_{i}^{N} left(frac{2}{5}right) m_i r_i^2 + m_i(x_i^2 + y_i^2)$$

Where $m_i$, $r_i$, $x_i$, $y_i$,and $z_i$ are the mass, radius, and x,y,z coordinates of the $i^{th}$ sphere respectively.

However I quickly realized that this only yields me the diagonal terms of the inertia matrix. I’m unsure how to calculate the off-diagonal terms. This answer leads me to believe that my off diagonal terms would be calculated as:

$$I_{xy} = sum_{i}^{N} left(frac{2}{5}right) m_i r_i^2 + m_i x_i y_i$$
$$I_{xz} = sum_{i}^{N} left(frac{2}{5}right) m_i r_i^2 + m_i x_i z_i$$
$$I_{yz} = sum_{i}^{N} left(frac{2}{5}right) m_i r_i^2 + m_i y_i z_i$$

However it is framed as a continuous integral for a continuous medium. While I recognize that a continuous integral is analogous to a finite sum such as what I outlined above, I’m having a difficult time figuring out if this applies to a collection of spheres/cubes. Or if that second term in my sum ($m_i x_i y_i$) needs to be modified at all.

Any help would be greatly appreciated! I’d just like to make sure that I am proceeding correctly!

---

Sphere Packing Illustration

For reference of how the sphere packing and cube packing words, see the below image. It shows the initial geometry represented as a polygon, the geometry represented as a dense packing of spheres, and then the geometry represented as a packing of cubes. NOTE: I am using larger cubes/spheres than I would for the actual calculations above simply so that they are visible on this page at the resolution available to me.