How can one intuitively think about quaternions? Quaternions came up while I was interning not too long ago and it seemed like no one really know how they worked. While eventually certain people were tracked down and were able to help with the issue, it piqued my interest in quaternions.
After reading many articles and a couple books on them, I began to know the formulas associated with them, but still have no clue how they work (why they allow rotations in 3D space to be specific). I back-tracked a little bit and looked at normal complex numbers with just one imaginary component and asked myself if I even understood how they allow rotations in 2D space. After a couple awesome moments of understanding, I understood it for imaginary numbers, but I'm still having trouble extending the thoughts to quaternions.
How can someone intuitively think about quaternions and how they allow for rotations in 3D space?
 A: I highly recommend that your read the presentation in Conway and Smith: On quaternions and octonions: Their geometry, arithmetic, and symmetry. Here's an excerpt from John Baez's very informative review:

It follows that the quaternions of norm 1 form a group under multiplication. This group is usually called SU(2), because people think of its elements as 2 × 2 unitary matrices with determinant 1. However, the quaternionic viewpoint is better adapted to seeing how this group describes rotations in 3 and 4 dimensions. The unit quaternions act via conjugation as rotations of the 3d space of "pure imaginary" quaternions, namely those with Re(q) = 0. This gives a homomorphism from SU(2) onto the 3d rotation group SO(3). The kernel of this homomorphism is {±1}, so we see SU(2) is a double cover of SO(3). The unit quaternions also act via left and right multiplication as rotations of the 4d space of all quaternions. This gives a homomorphism from SU(2) × SU(2) onto the 4d rotation group SO(4). The kernel of this homomorphism is {±(1, 1)}, so we see SU(2) × SU(2) is a double cover of SO(4).
These facts are incredibly important throughout mathematics and physics. With their help, Conway and Smith classify the finite subgroups of the 3d rotation group SO(3), its double cover SU(2), the 3d rotation/reflection group O(3), and the 4d rotation group SO(4). These classifications are all in principle "well known'. However, they seem hard to find in one place, so Conway and Smith's elegant treatment is very helpful

See also the Wikipedia article Quaternions and spatial rotation.
A: I don't know if the following counts as "intuition", but perhaps it
can put things in some perspective.
It's a theorem that any rotation in $R^n$ can be written as the
composition of an even number of reflections through hyperplanes.
To construct a way of doing rotions by algebra (not matrix algebra,
but rather by multiplying vectors and similar geometric objects), it
is therefore enough to figure out two things:
(1) Given a vector $f$, how do we represent the reflection $R_f$
through the hyperplane orthogonal to $f$? We want $R_f(x)$ to be some
kind of product involving the vectors $x$ and $f$.
(2) The correspondence constructed in (1) should be such that for two
vectors $f$ and $g$, the algebraic object $fg$ (whatever that is)
should represent the composition $R_f \circ R_g$. Then any rotation
can be written as
$R_{f_1} \circ R_{f_2} \circ \cdots \circ R_{f_{2m}}$ and therefore represented by the object
$f_1 f_2 \cdots f_{2m}$.
For this to work, we have to give up the idea that the product of two
vectors must be a vector (as Hamilton famously struggled with for a
long time). Instead we let the vectors sit inside some bigger space of
"complex" (compound) objects that can be multiplied with each other
according to suitable rules (associative, but not commutative).
In the case of quaternions, these complex objects consist of a "scalar
part" plus a "vector part", and the construction (1) which makes
everything work is $R_f(x)=-fxf^{-1}$. To see that this is a
reflection, note that it's linear in $x$, it maps $f$ to $-f$, and it
maps vectors orthogonal to $f$ to themselves (since they satisfy
$fx=-xf$ because of the way quaternion multiplication is defined).
To do rotations in $R^n$ one constructs a so-called Clifford algebra,
consisting of complex objects with a "scalar part", a "vector part", a
"bivector part", a "trivector part", etc., up to an "n-vector part",
and the rules to multiply these objects are generated by the
requirement that for vectors $x$ and $y$, one wants $xy+yx=2 x\cdot y$
(or the negative of that, depending on convention). This makes
orthonal vectors anticommute, so that $R_f(x)=-fxf^{-1}$ gives a
hyperplane reflection, and everything works as desired.
I won't go into what a "k-vector" is, but the space of k-vectors
is $\binom{n}{k}$-dimensional, so that the whole Clifford algebra is
$2^n$-dimensional.
The Clifford product of two pure vectors is a scalar plus a bivector,
which is a little bit different from quaternions where vector times
vector equals a scalar plus a vector.
The full Clifford algebra for $R^3$ is 8-dimensional, and it's a bit
of a "coincidence" that one can actually manage just as well with the
4-dimensional quaternion algebra if one only wants to work with
rotations in $R^3$.
A: One way to see the analogy between the transition from real to complex numbers and from complex numbers to quaternions
is by means of the fact that the spaces of unit complex numbers and the unit quaternions, namely,
the unit circle $S^1$ and the unit 3-sphere $S^3$ respectively, are Lie groups.
The analogy can be seen in the group laws, which can be casted in a very similar form as follows:
A point on $S^1$ can be represented by the complex number $z = (x, y)$ , $x^2+y^2= 1$.
A point on $S^3$ can be represented by a quaternion $q = (z, w)$ , where $z$ and $w$ are complex numbers and
$|z|^2+|w|^2= 1$.
The $S^1$ multiplication law has the form:
$z_1 z_2 = (x_1x_2 - y_1 y_2, x_1y_2 + y_1 x_2)$
The $S^3$ multiplication law is:
$q_1 q_2 = (z_1z_2 - w_1 \bar{w_2}, z_1w_2 + w_1 \bar{z_2})$
It easy to verify that these multiplication laws satisfy the group laws. 
Of course the main difference is that the quaternion group law is no longer commutative: since
$q_2 q_1 = (z_1z_2 - \bar{w_1} w_2, z_1w_2 + \bar{w_1} z_2)$
A: Re the last answer, in Geometric Algebra (a Clifford Algebra interpreted geometrically), this is very intuitive, you can even visualize it.
See http://www.geometricalgebra.net/quaternions.html
A: Here's one way.  The group of unit quaternions is isomorphic to the special unitary group $\text{SU}(2)$, the group of $2 \times 2$ unitary complex matrices with determinant $1$.  This group acts on $\mathbb{C}^2$ in the obvious way, and so it also acts on lines in $\mathbb{C}^2$.  (These are complex lines, so they have real dimension $2$.)  The space of lines in $\mathbb{C}^2$ is the complex projective line $\mathbb{CP}^1$, and it turns out there is a natural way to think about this space as a sphere - namely, the Riemann sphere.  There is a beautiful projection which is pictured at the Wikipedia article which shows this; essentially one thinks of $\mathbb{CP}^1$ as $\mathbb{C}$ plus a "point at infinity" and then projects the latter onto the former in a way which misses one point.
So $\text{SU}(2)$ naturally acts on a sphere, and as it turns out it naturally acts by rotations.  This describes the famous 2-to-1 map $\text{SU}(2) \to \text{SO}(3)$ which allows quaternions to describe 3D rotations.
A: I find the conversion between quaternions and the axis-angle representation quite instructive.
In the axis-angle representation, you describe a rotation by specifying the axis of rotation as a unit vector $\vec\omega$ and an angle $\theta$ about which to rotate around this axis. An interesting fact is that any possible rotation can be described in this way.
The corresponding quaternion is given simply by $\left(\cos(\theta/2), \vec\omega\sin(\theta/2)\right)$. Here the notation $(a, \vec v)$, where $a$ is a scalar and $\vec v$ a real vector, denotes the quaternion $a + v_xi + v_yj + v_zk$, or $(a,v_x,v_y,v_z)$.
Here is the intuitive interpretation of this. Given a particular rotation axis $\omega$, if you restrict the 4D quaternion space to the 2D plane containing $(1,0,0,0)$ and $(0,\omega_x,\omega_y,\omega_z)$, the unit quaternions representing all possible rotations about the axis $\vec \omega$ form the unit circle in that plane. A rotation of $\theta$ about the axis $\vec \omega$ is the point at an angle $\theta/2$ from $(1,0,0,0)$ on that circle. For example, not rotating at all is $(1,0,0,0)$, rotating 180° is $(0,\omega_x,\omega_y,\omega_z)$, and rotating 360° is $(-1,0,0,0)$, which is the same as not rotating at all (see final paragraph).
Multiplying two quaternions is unintuitive, but I'm not bothered by that, because the composition of two rotations in real life is quite unintuitive in the first place.
(Note that the quaternions are a "double cover" of the space of rotations, in that any rotation actually has two quaternions, say $q$ and $-q$, that represent it: $\theta$ and $\theta + 2\pi$ are the same angle, but $\theta/2$ and $(\theta+2\pi)/2$ are not. This is the only "glitch" in the quaternion representation of rotations.)
A: Thinking about quaternions as 4D is misleading. Quaternions are the union of a scalar and a 3-vector. Think: time and space. Space is a 3-vector. You can point in directions in space. Time is a scalar. There is a past (negative time) and a future (positive time) and now (0), but no ability to point in the direction of time.
Think of a blinking light on a train, each event having a time and location. These events can be written as quaternions. If the train travels at a constant velocity, you are seeing the addition of quaternions.
One can make movies out of quaternions. Examples are available on my web site, http://visualphysics.org
The rotations in 3D space calculations ignore time.
