Here is an algebraic theoretical proof parallel to the "gaussian" proof of Fermat's theorem on sums of squares. You may find it too complicated, but actually the general trend in math is to try to embed particular questions into larger and larger general theories in order to solve in one stroke whole families of problems which are of the same nature (*)
First replace the Gauss ring $\mathbf Z[i]$ by the Eisenstein ring $E=\mathbf Z[j]$, which is the ring of integers of the quadratic field $\mathbf Q(j)$, where $j$ is a primitive 3rd root of unity. It is classically known that $E$ is a principal ideal ring, so the theory of division in $E$ for ideals and for elements is the same up to the units of $E$, which are plainly the powers of $\pm j$. The prime elements (ideals) of $E$ are deduced from those of $\mathbf Z$ in the following manner : a prime $p$ of $\mathbf Z$ does not necessarily remain prime in $E$, but it can decompose only in 3 possible ways : (1) $p$ is inert, i.e. remains prime, iff $p\equiv 2$ mod $3$ ; (2) $p$ splits, i.e. $p=\pi\bar\pi$ (where $\pi$ is a prime of $E$) iff $p\equiv 1$ mod $3$ ; (3) $p$ ramifies in $E$, i.e. $p=\pi^2$ iff $p=3$ (see any textbook in ANT).
Now the equation $a^2-ab+b^2=c^2$ in $\mathbf Z$ reads $(a+bj)(a+bj^2)=c^2$ in $E$, where the UFD property implies that the prime divisors $\pi$ of $c$ are exactly such that $\pi$ divides $a+bj$, say, and the conjugate $\bar \pi$ divides $a+bj^2$. By uniqueness again, $\pi\bar \pi$ will divide $c$, i.e. the rational primes $p$ which divide $c$ must split.
(*) An archetypical example lies in the attempts to solve FLT by the same approach as above. Actually Kummer showed FLT along these lines - but of course with more involved arguments - in the so called regular case, i.e. when the ring $\mathbf Z[\zeta_p]$ (where $\zeta_p$ is a primitive $p$-th root of unity) is principal.