Given a measurable region $S\subset\mathbb R^n$, we define a symmetric linear transformation $M$ by a volume integral
$$M(v)=\frac{\int_{x\in S}\,x(x\cdot v)\,dV}{\int_{x\in S}\,dV}$$
which has coordinates, with respect to an orthonormal basis $\{e_i\}$,
$$M_{ij}=e_i\cdot M(e_j)=\frac{\int\,(x\cdot e_i)(x\cdot e_j)\,dV}{\int dV}=\frac{\int\,x_ix_j\,dV}{V}.$$
The rotational inertia tensor used in physics is related to $M$; for rotating $S$ in the $e_1e_2$ plane (or around the $e_3$ axis, assuming $n=3$), the inertia (normalized by the total volume and a constant density factor) is $M_{11}+M_{22}$.
Since any region can be approximated(?) as a polytope, which can be broken into simplices, let's assume that $S$ is a simplex, with vertices $a_0,a_1,a_2,\cdots,a_n$. Now how can $M(v)$ be expressed in terms of $a_i$?
The centroid of $S$ is $c=\tfrac1V\int x\,dV=\tfrac1{n+1}(a_0+a_1+a_2+\cdots+a_n)$. If we define
$$M_c(v)=\frac{\int\,(x-c)\big((x-c)\cdot v\big)\,dV}{V},$$
then we have a form of the parallel axis theorem:
$$M(v)=M_c(v)+c(c\cdot v)$$
so without loss of generality $c=0$, that is, $a_0=-a_1-a_2-\cdots-a_n$. (But we may instead want to set $a_0=0$.)
The volume of $S$ is
$$V=\tfrac1{n!}\lVert(a_1-a_0)\wedge(a_2-a_0)\wedge\cdots\wedge(a_n-a_0)\rVert$$
$$=\tfrac1{n!}\lVert(a_1\wedge a_2\wedge\cdots\wedge a_n)-(a_0\wedge a_2\wedge\cdots\wedge a_n)-(a_1\wedge a_0\wedge\cdots\wedge a_n)-\cdots-(a_1\wedge a_2\wedge\cdots\wedge a_0)\rVert$$
$$=\tfrac1{n!}\big\lVert\big(a_1\wedge a_2\wedge\cdots\wedge a_n\big)-\big((-a_1)\wedge a_2\wedge\cdots\wedge a_n\big)-\big(a_1\wedge(-a_2)\wedge\cdots\wedge a_n\big)-\cdots-\big(a_1\wedge a_2\wedge\cdots\wedge(-a_n)\big)\big\rVert$$
$$=\tfrac{n+1}{n!}\lVert a_1\wedge a_2\wedge\cdots\wedge a_n\rVert.$$
The case $n=1$ is trivial; $S$ is a line segment with endpoints $a_0,a_1$, so (treating vectors as scalars)
$$M=\frac{\int_{a_0}^{a_1}x^2\,dx}{\int_{a_0}^{a_1}dx}=\tfrac13(a_0^2+a_0a_1+a_1^2).$$
For $n=2$ ($S$ is a triangle), I chose coordinates where $a_0=0$ and $a_1$ is on the $e_1$ axis, and calculated
$$M_{11}=\tfrac16(a_{11}\!^2+a_{11}a_{21}+a_{21}\!^2)$$
$$M_{12}=\tfrac1{12}(a_{11}+2a_{21})a_{22}$$
$$M_{22}=\tfrac16a_{22}\!^2.$$
For higher $n$, these coordinate-based expressions would be more complicated. Is there a simple, general, coordinate-free expression?
Consider what happens to $M$ when a linear transformation $L$ is applied to $S$:
$$S'=L(S)=\{L(x)\mid x\in S\}$$
$$M'(v)=\frac{\int_{x\in S}\,L(x)\big(L(x)\cdot v\big)\,(\det L)dV}{\int_{x\in S}\,(\det L)dV}$$
$$=\frac{\int\,L(x)\big(L(x)\cdot v\big)\,dV}{\int dV}$$
$$=L\left(\frac{\int\,x\big(L(x)\cdot v\big)\,dV}{\int dV}\right)$$
$$=L\left(\frac{\int\,x\big(x\cdot L^*(v)\big)\,dV}{\int dV}\right)$$
$$=L\Big(M\big(L^*(v)\big)\Big)$$
where $L^*$ is the adjoint (or transpose) of $L$. Now we can take $S$ to be a regular simplex with circumradius $\lVert a_i\rVert=1$, since any two simplices $S,S'$ are related by some linear transformation.
By symmetry, all vertices $a_0,a_1,\cdots,a_n$ are eigenvectors of $M$ with the same eigenvalue $\lambda$, so $M(v)=\lambda v$.
We need a lemma:
$$a_0(a_0\cdot v)+a_1(a_1\cdot v)+a_2(a_2\cdot v)+\cdots+a_n(a_n\cdot v)=\frac{n+1}nv,$$
which can be verified by evaluating it on $a_i$ (since they span $\mathbb R^n$):
$$a_0(a_0\cdot a_0)+a_1(a_1\cdot a_0)+a_2(a_2\cdot a_0)+\cdots+a_n(a_n\cdot a_0)$$
$$=a_0(1)+a_1(-\tfrac1n)+a_2(-\tfrac1n)+\cdots+a_n(-\tfrac1n)$$
$$=a_0+\tfrac1n(-a_1-a_2-\cdots-a_n)$$
$$=a_0+\tfrac1n(a_0)=\frac{n+1}na_0.$$
And $a_i\cdot a_j=-\tfrac1n$ for $i\neq j$ because, for example,
$$1=a_0\cdot a_0=a_0\cdot(-a_1-a_2-\cdots-a_n)=-a_0\cdot a_1-a_0\cdot a_2-\cdots-a_0\cdot a_n$$
and these dot products are all the same by symmetry;
$$=-n\,a_0\cdot a_1.$$
Therefore,
$$M(v)=\lambda\frac n{n+1}\big(a_0(a_0\cdot v)+a_1(a_1\cdot v)+\cdots+a_n(a_n\cdot v)\big).$$
Now, for a general simplex $S'$ with vertices $a_i'=L(a_i)$, we have
$$M'(v)=L\Big(M\big(L^*(v)\big)\Big)$$
$$=L\left(\frac{\lambda n}{n+1}\Big(a_0\big(a_0\cdot L^*(v)\big)+a_1\big(a_1\cdot L^*(v)\big)+\cdots+a_n\big(a_n\cdot L^*(v)\big)\Big)\right)$$
$$=\frac{\lambda n}{n+1}\Big(L(a_0)\big(L(a_0)\cdot v\big)+L(a_1)\big(L(a_1)\cdot v\big)+\cdots+L(a_n)\big(L(a_n)\cdot v\big)\Big)$$
$$=\frac{\lambda n}{n+1}\big(a_0'(a_0'\cdot v)+a_1'(a_1'\cdot v)+\cdots+a_n'(a_n'\cdot v)\big).$$
This looks good, but we still don't know what $\lambda$ is.
In preparation for shifting the origin from the centroid to the vertex $a_0'$, define $b_i'=a_i'-a_0'$ for $i\neq0$, as the edges connecting to the vertex, and $b_0'=-b_1'-b_2'-\cdots-b_n'=(n+1)a_0'$.
$$M'(v)=\frac{\lambda n}{n+1}\Big(a_0'(a_0'\cdot v)+(b_1'+a_0')\big((b_1'+a_0')\cdot v\big)+\cdots+(b_n'+a_0')\big((b_n'+a_0')\cdot v\big)\Big)$$
$$=\frac{\lambda n}{n+1}\Big(a_0'(a_0'\cdot v)+b_1'(b_1'\cdot v)+\cdots+b_n'(b_n'\cdot v)+(b_1'+\cdots+b_n')(a_0'\cdot v)+a_0'\big((b_1'+\cdots+b_n')\cdot v\big)+n\,a_0'(a_0'\cdot v)\Big)$$
$$=\frac{\lambda n}{n+1}\Big(b_1'(b_1'\cdot v)+\cdots+b_n'(b_n'\cdot v)-(n+1)a_0'(a_0'\cdot v)\Big).$$
Let's shift and apply the parallel axis theorem:
$$M_{a_0'}'(v)=M'(v)+(-a_0')\big((-a_0')\cdot v\big)$$
$$=\frac{\lambda n}{n+1}\big(b_1'(b_1'\cdot v)+\cdots+b_n'(b_n'\cdot v)\big)+(1-\lambda n)a_0'(a_0'\cdot v)$$
$$=\frac{\lambda n}{n+1}\big(b_1'(b_1'\cdot v)+\cdots+b_n'(b_n'\cdot v)\big)+\frac{1-\lambda n}{(n+1)^2}b_0'(b_0'\cdot v).$$
To find $\lambda$, we'll cut the regular simplex $S$ into $n+1$ smaller simplices $S_0,S_1,\cdots,S_n$ (using the centroid as a new vertex), each of which has volume $\tfrac1{n+1}V$, and apply the above result to each one:
$$M(v)=\frac{\int_{x\in S}\,x(x\cdot v)\,dV}{V}$$
$$=\frac1V\int_{x\in S_0}x(x\cdot v)\,dV+\frac1V\int_{x\in S_1}x(x\cdot v)\,dV+\cdots+\frac1V\int_{x\in S_n}x(x\cdot v)\,dV$$
$$=\tfrac1{n+1}M^{(S_0)}(v)+\tfrac1{n+1}M^{(S_1)}(v)+\cdots+\tfrac1{n+1}M^{(S_n)}(v)$$
(I could have picked a better notation for $M$'s dependence on $S$ and the choice of origin)
$$=\frac1{n+1}\left(\frac{1-\lambda n}{(n+1)^2}a_0(a_0\cdot v)+\frac{\lambda n}{n+1}\big(a_1(a_1\cdot v)+a_2(a_2\cdot v)+\cdots+a_n(a_n\cdot v)\big)\right)$$
$$+\frac1{n+1}\left(\frac{1-\lambda n}{(n+1)^2}a_1(a_1\cdot v)+\frac{\lambda n}{n+1}\big(a_0(a_0\cdot v)+a_2(a_2\cdot v)+\cdots+a_n(a_n\cdot v)\big)\right)$$
$$+\frac1{n+1}\left(\frac{1-\lambda n}{(n+1)^2}a_2(a_2\cdot v)+\frac{\lambda n}{n+1}\big(a_0(a_0\cdot v)+a_1(a_1\cdot v)+\cdots+a_n(a_n\cdot v)\big)\right)$$
$$\vdots$$
$$+\frac1{n+1}\left(\frac{1-\lambda n}{(n+1)^2}a_n(a_n\cdot v)+\frac{\lambda n}{n+1}\big(a_0(a_0\cdot v)+a_1(a_1\cdot v)+\cdots+a_{n-1}(a_{n-1}\cdot v)\big)\right)$$
$\,$
$$=\frac1{n+1}\left(\frac{1-\lambda n}{(n+1)^2}\big(a_0(a_0\cdot v)+a_1(a_1\cdot v)+\cdots+a_n(a_n\cdot v)\big)+\frac{\lambda n}{n+1}\big(n\,a_0(a_0\cdot v)+n\,a_1(a_1\cdot v)+\cdots+n\,a_n(a_n\cdot v)\big)\right)$$
$\,$
$$=\frac1{n+1}\left(\frac{1-\lambda n+\lambda n^3+\lambda n^2}{(n+1)^2}\big(a_0(a_0\cdot v)+a_1(a_1\cdot v)+\cdots+a_n(a_n\cdot v)\big)\right)$$
$$=\frac{1+\lambda n(n^2+n-1)}{(n+1)^2n}\,v.$$
Thus
$$\lambda=\frac{1+\lambda n(n^2+n-1)}{(n+1)^2n}$$
which we can solve to find
$$\lambda=\frac1{n(n+2)}.$$
Finally, we get a formula for a general simplex with centroid at the origin:
$$M'(v)=\frac1{(n+1)(n+2)}\big(a_0'(a_0'\cdot v)+a_1'(a_1'\cdot v)+\cdots+a_n'(a_n'\cdot v)\big).$$
Letting $d_i'=a_i'-h$ be the vectors from an arbitrary point $h$ to the vertices, and going through the "shift" calculations again, we get a more general formula:
$$M_h'(v)=\frac1{(n+1)(n+2)}\big(d_0'(d_0'\cdot v)+d_1'(d_1'\cdot v)+\cdots+d_n'(d_n'\cdot v)\big)+\frac{n+1}{n+2}h(h\cdot v)$$
or, reverting to the notation from the OP, with $h=-c$,
Sloane and Conway wrote about this, but they only considered $I=\text{tr }M$. Their equation (15) agrees with mine.