This is an extended comment. (The answer is contained in the comment of Will Jagy:
indeed it is a very simple argument that every rational orthogonal matrix is obtained from a rational quaternion).

The paper of J. Cremona referred in the question gives an incorrect formula for the recovery of the quaternion. The correct formula in Cremona's notation is
$$\cos(\theta/2)+(b\mathbf{i}+c\mathbf{j}+d\mathbf{k})\sin(\theta/2),$$
where $(b,c,d)$ is the *unit* vector in the direction of the axis of rotation,
and $\theta$ is the angle of rotation. It is far from evident from this formula
that all rational rotations matrices arise from rational quaternions.
But the simple direct argument of Pall does not use this inversion formula.

The fact that all rotation matrices with *rational* entries can be parametrized by the formula in question
was discovered by Euler, Problema algebraicum ad affectiones prorsus singulares memorabile, Novi Commentarii Avcad. Sci. Imp. Petropolitanae, XV, 75-106 (Opera omnia, ser. 1, vol. 6, 287-315).

He also did the similar parametrization of rational $4\times 4$ rotations matrices.

All this was 72 years before the discovery of quaternions!

Let me cite Euler in English translation:

"This solution deserves full attention, as I arrived at it not using a definite method but rather by making guesses... If somebody can find a direct way to
this solution, then it will have to be admitted that he made an outstanding contribution to analysis."

According to Khrushchev (Orthogonal polynomials and continued fractions, Cambridge UP, 2008), the first proof is due to D. Grave (1937).

In the book of Vilenkin Special functions and the theory of group representations, one can find a generalization to arbitrary dimension.