I will be using the method given by Phira in the question you linked. First let us assume that $x \ge y \ge z$
- Number of squares of size 1: $x*y*z$
- Number of squares of size 2: $(x-1) * (y-1) * (z-1)$
- ...
- Number of squares of size $z$: $(x-z+1)*(y-z+1) $
We can add these together to get the total number of cubes, this sum can be written as:
$$ f(x,y,z) = \sum_{i=0}^{z-1} (x-i)(y-i)(z-i) $$
However, I will re-index using $k=z-i$ and we get the sum:
$$ f(x,y,z) = \sum_{k=1}^z k(x-z+k)(y-z+k) $$
Distributing, we get:
$$ f(x,y,z) = \sum_{k=1}^z k^3 + k^2\big((y-z)+(x-z)\big) + k(x-z)(y-z) $$
Expressing this as 3 different sums each with a different degree of $k$, we get:
$$ f(x,y,z) = \sum_{k=1}^z k^3 + (x+y-2z) \sum_{k=1}^z k^2 + (x-z)(y-z) \sum_{k=1}^z k$$
These are all known sums now so we can simplify this greatly.
$$ f(x,y,z) = \frac{z^2(z+1)^2}{4} + (x+y-2z) \frac{z(z+1)(2z+1)}{6} + (x-z)(y-z)\frac{z(z+1)}{2} $$
As you can see this is a closed form solution. It is not pretty and the rest of my work will just be algebreaic manipulation to try to make it looks pretty. (You should check for mistakes youself, a lot of terms means a lot of places I could've made a mistake)
$$ f(x,y,z)= \frac{z(z+1)}{12}\left( 3z^2+3z + 2x+2y-4z + 4xz+4yx-8z^2 + 6xy - 6xz - 6yz +6z^2 \right) $$
$$ = \frac{z(z+1)}{12} \left( z^2 - z + 2x + 2y +6xy - 2xz - 2 yz \right) $$
$$ = \frac{z(z+1)xy}{2} + \frac{(z-1)z(z+1)}{12}(z-2x-2y) $$
