All the primitive Pythagorean Quadruples are known. This is Theorem 3 on page 176 and Theorem 4 on page 177 of Jones_Pall_1939.pdf, available at TERNARY as a pdf. The same information is on the first two pages of Pall_Automorphs_1940.pdf at the same site.
The short version is this: you have an odd number $ W $ Find all quadruples $a,b,c,d$ with $$ a^2 + b^2 + c^2 + d^2 = W, $$ where we are allowed to mix order, take the variables to be positive, negative, or zero. Then all the primitive quadruples, odd entry first, are given by
$$ \left( a^2 + b^2 - c^2 - d^2 \right)^2 + 4 \left( ad-bc \right)^2 + 4 \left(ac+bd \right)^2 = W^2. $$
Meanwhile, if $t$ is some (positive) divisor of $W,$ we can do the same process for $W/t$ and $W^2 / t^2,$ then multiply all $a,b,c,d$ by $t.$
July 2015: for another project, I decided to try to generate all the quadruples by the sort of three-parameter formulas one gets by stereographic projection to $\mathbb S^2.$ the results are really disappointing. I think I will stick with the four parameter thing. I found it in Jones and Pall, but it goes back at least to V. A. Lebesgue (and likely known to Euler), https://en.wikipedia.org/wiki/Pythagorean_quadruple#Parametrization_of_primitive_quadruples
and is a simple calculation using quaternions with ordinary integer coefficients. Alright, the (first correct) proof that the four parameter recipe gives all primitive solutions is from 1920, by Dickson. In 1941, Skolem gave a proof that can be adjusted to give a reasonably direct algorithm for taking a quadruple and reconstructing the four parameters. In 1956, one F. Steiger gave inequalities that make the mapping one-to-one. This is all reported in a 1962 article by Robert Spira, in the maa Monthly, May 1962, volume 69, number 5, pages 360-365, title The Diophantine Equation $x^2 + y^2 + z^2 = m^2.$ One thing I did not initially notice, Spira just discards the Pythagorean triples; in the quadruple setting, if one of the numbers is $0,$ we have a triple and can easily recover parameters.
Let's see, I was able to adapt the methods of Jones and Pall to $x^2 + y^2 + z^2 = 3 t^2,$ coefficients other than $3$ could also be handled this way. http://math.stackexchange.com/questions/1964607/when-will-a-parametric-solution-generate-all-possible-solutions/1965805#1965805
I will just give a sketch for now how to solve it, I can fill in more details if necessary.
First define vectors $v_i = (x_i-a,y_i-b),\, i = 1,2,3$. Then we can rewrite the equations as
\begin{align}
\langle v_1, v_2\rangle = a_{12},\ \langle v_2, v_3\rangle = a_{23},\ \langle v_3, v_1\rangle = a_{31},\\
\det[v_1^t\ v_2^t] = b_{12},\ \det[v_2^t\ v_3^t] = b_{23},\ \det[v_3^t\ v_1^t] = b_{31},\\
\end{align} where I changed dot's and det's to a's and b's.
Denote by $\vartheta_{ij}$ the (oriented) angle between vectors $v_i$ and $v_j$. Remember that dot product in $\mathbb R^2$ is given by formula $\|v\|\|w\|\cos\vartheta$ and that determinants above measure (signed) area of parallelogram formed by the vectors, which can also be expressed as $\|v\|\|w\|\sin\vartheta$, and thus the system becomes
\begin{align}
\|v_1\|\|v_2\|\cos\vartheta_{12} = a_{12},\ \|v_2\|\|v_3\|\cos\vartheta_{23} = a_{23},\ \|v_3\|\|v_1\|\cos\vartheta_{31} = a_{31},\\
\|v_1\|\|v_2\|\sin\vartheta_{12} = b_{12},\ \|v_2\|\|v_3\|\sin\vartheta_{23} = b_{23},\ \|v_3\|\|v_1\|\sin\vartheta_{31} = b_{31}.\\
\end{align}
Then, you can find the lengths $\|v_i\|$ by looking at the system $$\|v_1\|\|v_2\| = c_{12},\ \|v_2\|\|v_3\| = c_{23},\ \|v_3\|\|v_1\| = c_{31},$$ where $c_{ij} = \sqrt{a_{ij}^2+b_{ij}^2}$ obtained by squaring the above equations and using $\sin^2t+\cos^2t = 1$. Solving it gives you $$\|v_1\|=\sqrt{\frac{c_{12}c_{31}}{c_{23}}},\ \|v_2\|=\sqrt{\frac{c_{12}c_{23}}{c_{31}}},\ \|v_3\|=\sqrt{\frac{c_{23}c_{31}}{c_{12}}}.$$
The whole thing is, from geometric perspective, obviously rotationally invariant (rotating all of the vectors won't change the angles between them or the areas), so fix some angle $\vartheta$. Then, the solution can be represented as complex numbers as
$$v_1 = \|v_1\|e^{i\vartheta},\ v_2 = \|v_2\|e^{i(\vartheta +\vartheta_{12})},\ v_3 = \|v_3\|e^{i(\vartheta-\vartheta_{31})}.$$
More explicitly, $e^{i\vartheta_{12}} = \cos\vartheta_{12} + i\sin\vartheta_{12} = \frac{a_{12}}{c_{12}}+\frac{b_{12}}{c_{12}}i$ and $e^{i\vartheta_{31}} = \cos\vartheta_{31} + i\sin\vartheta_{31} = \frac{a_{31}}{c_{31}}+\frac{b_{31}}{c_{31}}i$, so we have
\begin{align}v_1 &= \sqrt{\frac{c_{12}c_{31}}{c_{23}}}(\cos\vartheta + i\sin\vartheta),\\
v_2 &= \sqrt{\frac{c_{12}c_{23}}{c_{31}}}(\cos\vartheta + i\sin\vartheta)(\frac{a_{12}}{c_{12}}+i\frac{b_{12}}{c_{12}}),\\
v_3 &= \sqrt{\frac{c_{23}c_{31}}{c_{12}}}(\cos\vartheta + i\sin\vartheta)(\frac{a_{31}}{c_{31}}-i\frac{b_{31}}{c_{31}}),\ \vartheta\in\mathbb R.\end{align}
All you have to do now is expand and $x$'s will be the real parts, while $y$'s the imaginary parts of the above complex number representation.
Note, however, that the system is overdetermined since knowing the angles between $v_1$ and $v_2$ and $v_2$ and $v_3$ will also give you the angle between $v_1$ and $v_3$.
Best Answer
There is a classical way of finding all rational solutions to a quadratic equation (such as the Pythagorean equation $x^2+y^2=1$). Fix one solution—in our case we'll fix $(x_1,x_2,x_3,y_1,y_2,y_3) = (-1,0,0,-1,0,0)$. Then for every other rational solution $(s_1,s_2,s_3,t_1,t_2,t_3)$, each of the "slopes" $m_1=\frac{s_1-(-1)}{t_3-0}$, $m_2=\frac{s_2-0}{t_3-0}$, $m_3=\frac{s_3-0}{t_3-0}$, $m_4=\frac{t_1-(-1)}{t_3-0}$, $m_5=\frac{t_2-0}{t_3-0}$ must be a rational number. On the other hand, given any five rational numbers $m_1,m_2,m_3,m_4,m_5$, the intersection of the quadratic surface $x_1^2+x_2^2+x_3^2=y_1^2+y_2^2+y_3^2$ with these five hyperplanes $m_1t_3=s_1+1$, $m_2t_3=s_2$, $m_3t_3=s_3$, $m_4t_3=t_1+1$, $m_5t_3=t_2$ yields a pair of points, one of which, namely $(-1,0,0,-1,0,0)$, we already know; solving the corresponding quadratic equation will give us the other point.
Therefore there is a bijective correspondence between $5$-tuples $(m_1,m_2,m_3,m_4,m_5)$ of rational numbers and rational solutions to $x_1^2+x_2^2+x_3^2=y_1^2+y_2^2+y_3^2$ other than $(-1,0,0,-1,0,0)$. Solving these equations yields \begin{align*} x_1 &= \frac{m_1^2-2 m_1m_4-m_2^2-m_3^2+m_4^2+m_5^2+1}{m_1^2+m_2^2+m_3^2-m_4^2-m_5^2-1} \\ x_2 &= \frac{2 m_2(m_1-m_4)}{m_1^2+m_2^2+m_3^2-m_4^2-m_5^2-1} \\ x_3 &= \frac{2 m_3(m_1 -m_4)}{m_1^2+m_2^2+m_3^2-m_4^2-m_5^2-1} \\ y_1 &= \frac{-m_1^2+2 m_1m_4-m_2^2-m_3^2-m_4^2+m_5^2+1}{m_1^2+m_2^2+m_3^2-m_4^2-m_5^2-1} \\ y_2 &= \frac{2 m_5(m_1-m_4)}{m_1^2+m_2^2+m_3^2-m_4^2-m_5^2-1} \\ y_3 &= \frac{2(m_1-m_4)}{m_1^2+m_2^2+m_3^2-m_4^2-m_5^2-1} \end{align*} If we write $m_j = \frac{n_j}d$ and clear all denominators, we get the parametrization in $6$ integer variables \begin{align*} x_1 &= n_1^2-2 n_1n_4-n_2^2-n_3^2+n_4^2+n_5^2 +d^2\\ x_2 &= 2 n_2(n_1-n_4) \\ x_3 &= 2 n_3 (n_1-n_4) \\ y_1 &= -n_1^2+2 n_1n_4-n_2^2-n_3^2-n_4^2+n_5^2+d^2 \\ y_2 &= 2 n_5(n_1-n_4) \\ y_3 &= 2 d (n_1-n_4). \end{align*} It would take extra work to produce a complete minimal parametrization (this one might hit solutions many times each). There are also some issues with gcds of the numbers. (Although since the original equation is homogeneous, finding integer and rational solutions is basically the same thing.)