Is is not that difficult to integrate, it is more difficult to find the perpendicular distance $r_\perp$ from a point $(x,y,z)$ of the cuboid to the space diagonal of the cuboid.
For elegance, denote the cuboid lengths as $a,b,c$ and place it into a Cartesian coordinate system with the origin at one vertex such that $x$-axis is along length $a$, $y$-axis is along length $b$ and $z$-axis is along length $c$.
Say we are interested in the space diagonal from the origin $(0,0,0)$ to the point $(a,b,c)$. Now we need to recall that the distance of point $A$ from the line from $B$ to $C$ is given by
$$d(A,BC) = \frac{\|\vec{BA} \times \vec{BC}\|}{\|\vec{BC}\|} = \frac{\|(A-B) \times (C-B)\|}{\|C-B\|}.$$
In our case, we are interested in the distance from $(x,y,z)$ to the space diagonal, which is
$$r_\perp = \frac{\|(x,y,z) \times (a,b,c)\|}{\|(a,b,c)\|} = \frac{(bx-ay)^2+(az-cx)^2+(cy-bz)^2}{\sqrt{a^2+b^2+c^2}}.$$
Now the moment of inertia is simply
\begin{align}
I_{\text{diag}} &= \int_{\text{cuboid}}r_\perp^2\,dm \\
&= \frac{m}{abc}\int_{\text{cuboid}} r_\perp^2\,dV \\
&= \frac{m}{abc(a^2+b^2+c^2)} \int_{x=0}^a \int_{y=0}^b\int_{z=0}^c (bx-ay)^2+(az-cx)^2+(cy-bz)^2\, dz\,dy\,dx\\
&= \frac16m\frac{a^2b^2+b^2c^2+c^2a^2}{a^2+b^2+c^2}
\end{align}
To verify that this result is indeed correct, we can proceed by the usual route as suggested in the comments. The moment inertia tensor of the cuboid around its center of mass is
$$\mathbf{I} = \frac{m}{12}\begin{bmatrix} b^2+c^2 & 0 & 0 \\ 0 & a^2+c^2 & 0 \\ 0 & 0 & a^2+b^2\end{bmatrix}$$
and the moment of inertia around an axis with unit vector $\hat{n}$ through the center of mass is given by $\mathbf{I}\hat{n}\cdot \hat{n}$. In our case, the space diagonal has vector $\hat{n} = \frac{(a,b,c)}{\|(a,b,c)\|} = \frac{(a,b,c)}{\sqrt{a^2+b^2+c^2}}$ so
$$I_{\text{diag}} = \mathbf{I}\hat{n}\cdot \hat{n} = \frac{m}{12}\frac1{a^2+b^2+c^2}((b^2+c^2)a^2+(a^2+c^2)b^2+(a^2+b^2)c^2) = \frac16m\frac{a^2b^2+b^2c^2+c^2a^2}{a^2+b^2+c^2}$$
which is of course the same result.
