I thought I would write out the quaternion proof, because it is really quite elegant if you have good algebraic terminology. One note: I won't be proving quite the right result. Let $r(n)$ be the number of quadruples $(a,b,c,d)$ with $n=a^2+b^2+c^2+d^2$ such that either $a$, $b$, $c$ and $d$ are all in $\mathbb{Z}$, or are all in $(1/2) + \mathbb{Z}$. The formula I'll be proving is
$$r(n) = 24 \sum_{\begin{matrix} d|n \\ d \ \mathrm{odd} \end{matrix}} d \quad \mathrm{for}\ n>0$$
The question of allocating the terms between $\mathbb{Z}$ and $(1/2) + \mathbb{Z}$ is a bit messier.
Let $H$ be the ring of quaternions of the form $a+bi+cj+dk$, with $(a,b,c,d)$ as above.
The following lemmas are proved in many sources on the quaternion proof of the Lagrange theorem:
$\bullet$ Every right ideal of $H$ is principal.
$\bullet$ The right ideal $(a+bi+cj+dk) H$ has index $(a^2+b^2+c^2+d^2)^2$ in $H$.
Let $q(n)$ be the number of right ideals of $H$ that have index $n^2$. Since the unit group of $H$ has size $24$, we have $r(n) = 24 q(n)$. So we concentrate on computing $q(n)$.
Let $I = (a+bi+cj+dk)H$ be a right ideal of $H$ with index $n^2$. Then $(a+bi+cj+dk)(a-bi-cj-dk) = n \in (a+bi+cj+dk)H$. So $I/nH$ is an index $n^2$ ideal in the ring $H/nH$. We see that $q(n)$ is the number of index $n^2$ ideals in $H/nH$.
Now, by the Chinese remainder theorem, if $m$ and $n$ relatively prime, then $H/mnH \cong H/mH \times H/nH$. So $q(n)$ is multiplicative, and we are reduced to proving the claim for $n$ a prime power.
We start in the case $n = p^r$ for $p$ odd. Let $\mathbb{Z}_p$ be the $p$-adic integers and let $H_p := H \otimes \mathbb{Z}_p$. We will start by showing that the number of index $p^{2r}$ ideals in $H_p$ is $p^r+p^{r-1} + \cdots + p+1$. We then check that all of those ideals contain $p^r H$, so this is also the number of ideals in $H/p^r H$.
The point is that $H_p \cong \mathrm{Mat}_{2 \times 2}(\mathbb{Z}_p)$ (again, for $p$ odd). By a standard pigeon hole argument, we can find $u$ and $v$ with $u^2+v^2 + 1 \equiv 0 \mod p$; by Hensel's lemma, we can find $u$ and $v$ in $\mathbb{Z}_p$ with $u^2+v^2+1=0$. Let
$$I = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \quad
J = \begin{pmatrix} u & v \\ v & -u \end{pmatrix} \quad
K = \begin{pmatrix} v & -u \\ -u & -v \end{pmatrix}.$$
Then $I$, $J$ and $K$ obey the quaternion relations, and we get a map $H_p \to \mathrm{Mat}_{2 \times 2}(\mathbb{Z}_p)$. We claim that this map is bijective. It is a map between two free $\mathbb{Z}_p$-modules of rank four, so we just need to check that the determinant of this map is a unit of $\mathbb{Z}_p$. This determinant is $-4(u^2+v^2)=4$. Since $p$ is odd, we are fine.
So, set $E_p := \mathrm{Mat}_{2 \times 2}(\mathbb{Z}_p)$. Let's count index $p^{2r}$ right ideals $E_p$. Let $M$ be the $E_p$-module $\mathbb{Z}_p^2$. As a module over itself, $E_p \cong M^{\oplus 2}$. Right ideals is another word for submodules, so we want to count index $p^{2r}$ submodules of $M^{\oplus 2}$. Since $E_p$ is Morita equivalent to $\mathbb{Z}_p$, we know that $E_p$ submodules of $M^{\oplus 2}$ are in bijection with $\mathbb{Z}_p$-submodules of $\mathbb{Z}_p^{\oplus 2}$. A quick check shows that being index $p^{2r}$ in $M^{\oplus 2r}$ corresponds to being index $p^r$ in $\mathbb{Z}_p^{\oplus 2}$. Moreover, every index $p^r$ submodule of $\mathbb{Z}_p^{\oplus 2}$ is contained in $p^r \mathbb{Z}_p^{\oplus 2}$, so every index $p^{2r}$ submodule of $M^{\oplus 2}$ is contained in $p^r M^{\oplus 2}$, fulfilling an earlier promise, and allowing us to shift our attention to the finite problem of counting index $p^r$ submodules in $(\mathbb{Z}/p^r)^{\oplus 2}$.
We are reduced to showing that the number of index $p^r$ submodules of $(\mathbb{Z}/p^r)^{\oplus 2}$ is $p^r + p^{r-1} + \cdots p+1$. Proof by induction on $r$. Let $N$ be such a submodule. There are two cases:
(1) $N$ is not contained in $p (\mathbb{Z}/p^r \mathbb{Z})^{\oplus 2}$. In this case, $N$ is generated by a single element. The number of elements of $(\mathbb{Z}/p^r)^{\oplus 2}$ which generate a submodule of index $p^r$ is $p^{2r} - p^{2r-2}$. Two such elements give the same submodule if their ratio is a unit of $\mathbb{Z}/p^r$; the number of such units is $p^r-p^{r-1}$. So the number of $N$ in this case is $(p^{2r}-p^{2r-2})/(p^r - p^{r-1}) = p^r+p^{r-1}$.
(2) $N$ is contained in $p (\mathbb{Z}/p^r \mathbb{Z})^{\oplus 2}$. In this case, $N$ is an index $p^{r-2}$ submodule of $p (\mathbb{Z}/p^r \mathbb{Z})^{\oplus 2} \cong (\mathbb{Z}/p^{r-1} \mathbb{Z})^{\oplus 2}$. By induction, the number of options for such an $r$ is $p^{r-2} + p^{r-3} + \cdots p+1$.
Adding the two cases together, we are done.
Okay, so what about $2^r$? What one wants to show is that there is only one index $2^{2r}$ ideal in $H$.
I haven't found a ring theory way to do this, but a direct proof isn't hard. Notice that, if $a^2+b^2+c^2+d^2$ is even, then $(a,b,c,d)$ must be integers, and an even number of $(a,b,c,d)$ must be odd. So we can write $a+b+cj+dk = \epsilon + 2 (a' + b'i+c'j+d'k)$ which $\epsilon$ is one of $0$, $1+i$, $1+j$, $1+k$, $i+j$, $i+k$, $j+k$ or $1+i+j+k$. Each option for $\epsilon$ is divisible by $(1+i)$, and we also have $2=(1+i)(1-i)$. So we have shown that, if $a^2+b^2+c^2+d^2$ is even, then $1+i$ divides $a+bi+cj+dk$. By induction, then, $a^2+b^2+c^2+d^2 = 2^r$, then $a+bi+cj+dk$ is a unit times $(1+i)^r$.
