The geometric interpretation of quaternion multiplication is fundamentally 4-dimensional (unlike quaternion conjugation, which can be considered as an action on $\Bbb{R}^3$).
Let's start with an easy case. Say $q=a+bi$ with $b \neq 0$, $a^2+b^2=1$. That is, $q$ is a non-real unit quaternion in the subalgebra of $\Bbb{H}$ generated by $i$. What affect does multiplying by $q$ have on an arbitrary quaternion $r$?
First of all, if $r$ also lies in the subalgebra generated by $i$, then we can consider multiplication of $q$ and $r$ to be ordinary complex multiplication; that is, multiplication by $q$ rotates $r$ by $\theta$ in the $\{1, i\}$ plane.
Secondly, if $r$ lies in the orthogonal complement of that subalgebra, $r=cj+dk$, we can write $r=(c+di)j=j(c-di)$. The first of these representations can be used to left-multiply by $q$ via ordinary complex multiplication; the second one can be used to right-multiply. In either case, multiplying by $q$ rotates $r$ by $\theta$ in the $\{j, k\}$-plane; however, the sign difference means that the two multiplications rotate in opposite directions from each other.
We can then find the effect of multiplying $q$ by an arbitrary quaternion by projecting that quaternion into these two planes. That is, an arbitrary quaternion will have its $\{1, i\}$-projection and $\{j, k\}$-projection both rotated by $\theta$ when it is multiplied by $q$; the direction of the $\{j, k\}$-rotation, but not of the $\{1, i\}$-rotation, will be affected by whether we're left-multiplying or right-multiplying by $q$.
The general case works similarly. For any unreal unit quaternion $q$ that makes an angle $\theta$ with the real axis, multiplication by $q$ rotates by $\theta$ in the $\{1, q\}$-plane, and also rotates by $\theta$ in its orthogonal complement. The direction of the first rotation is fixed, but the direction of the second rotation depends on whether we're multiplying by $q$ on the left or on the right. You can see this just by noticing that any unreal quaternion generates a 2-dimensional subalgebra isomorphic to $\Bbb{C}$, making the previous few paragraphs work in general after some relabeling.
This also gives you a way to see why quaternion conjugation works the way it does on $\Bbb{R}^3$. If $q$ makes an angle $\theta$ with the real axis, then the map $r \mapsto qrq^{-1}$:
- is the identity map in the $\{1, q\}$-plane, since that plane forms a commutative subalgebra of $\Bbb{H}$
- involves a rotation by $2\theta$ in its orthogonal complement. Both $q$ and $q^{-1}$ rotate quaternions in $\{1, q\}^{\perp}$ by $\theta$. If they were multiplied on the same side, those rotations would have to cancel out; since left-multiplication behaves oppositely to right-multiplication, this means they must reinforce each other. Since the orthogonal complement of $\{1, q\}$ is orthogonal to $1$, it is pure imaginary, so we've reproduced the fact that quaternion conjugation corresponds to a double-angle rotation in $\Bbb{R}^3$ (when identified with $\Im(\Bbb{H})$).
Note also that multiplying by a general quaternion involves scaling by the norm of that quaternion, but conjugation conveniently causes the norms of $q$ and $q^{-1}$ to cancel.
For any finite-dimensional real vector space $V$, if $g$ is an inner product on $V$ (i.e $g:V\times V\to\Bbb{R}$, is bilinear, symmetric and positive definite), then we get an induced linear mapping $g^{\flat}:V\to V^*$ (read as 'gee-flat'), and the definition is
\begin{align}
g^{\flat}(v):=g(v,\cdot).
\end{align}
The meaning is that $g$ is a function of two variables, taking values in $\Bbb{R}$, so what we can do is fill up one of the entries to get $g(v,\cdot)$, and this can eat another vector $w$ to product a real number, so $g^{\flat}(v)=g(v,\cdot)\in V^*$.
Now, the mapping $g^{\flat}$ is injective, because if $v\in \ker(g^{\flat})$, then it means $g^{\flat}(v):=g(v,\cdot)=0$, which explicitly means for all $w\in V$, $g(v,w)=0$. So, in particular, $g(v,v)=0$, hence by positive-definiteness, $v=0$. Therefore, $\ker(g^{\flat})=\{0\}$, hence we have an injection. Recall that on a finite-dimensional space, $V$ and $V^*$ have the same dimension, so an injection is automatically an isomorphism. The inverse mapping $(g^{\flat})^{-1}:V^*\to V$ is typically denoted $g^{\sharp}:V^*\to V$.
We're going to apply this logic in the specific case $V=\Bbb{R}^3$ and $g(x,y)=\sum_{i=1}^3x_iy_i$ being the standard inner product. So, we have an isomorphism $g^{\flat}:\Bbb{R}^3\to(\Bbb{R}^3)^*$. Now, the definition of the cross product is that given $u,v\in\Bbb{R}^3$, we recall that the determinant is a multilinear function of the rows/columns, so the function $\det(u,v,\cdot):\Bbb{R}^3\to\Bbb{R}$ is linear. This means $\det(u,v,\cdot)\in (\Bbb{R}^3)^*$. Thus, the definition of the cross product is
\begin{align}
u\times v&=g^{\sharp}\bigg(\det(u,v,\cdot)\bigg)\in \Bbb{R}^3
\end{align}
So, it is the fact that $\Bbb{R}^3\cong (\Bbb{R}^3)^*$ via the inner product which makes it possible to define the cross product in this way.
Anyway, if all you want to know is how to calculate the cross product, then recall that in general, given any vector $\xi\in \Bbb{R}^3$, we can take any orthonormal basis, for example the standard basis $\{e_1,e_2, e_3\}$, and write $\xi=g(\xi,e_1)e_1+g(\xi,e_2)e_2+g(\xi,e_3)e_3$ (i.e dot the vector into the basis to get the projections onto the axis, and then multiply by the unit vector).
So, for the cross product, we can also do the same:
\begin{align}
u\times v&=\det(u,v,e_1)\cdot e_1+\det(u,v,e_2)\cdot e_2 + \det(u,v,e_3)\cdot e_3
\end{align}
You may have seen this summarized mnemonically as
\begin{align}
u\times v&=
\det
\begin{pmatrix}
u_1&u_2&u_3\\
v_1&v_2&v_3\\
\hat{x}&\hat{y}&\hat{z}
\end{pmatrix}
\end{align}
Best Answer
Before meeting the cross product, students will have already met the dot product. This is a form of multiplication that takes two vectors and gives you a scalar. It's natural to wonder whether there's also some kind of multiplication that takes two vectors and gives you another vector.
Obviously we could just write down any old function and call it "the cross product". But in order for it to actually be a nice form of multiplication there are a few properties that we would expect it to have:
$$\begin{align*} (\lambda\mathbf u) \times \mathbf v=\lambda(\mathbf u \times \mathbf v)=\mathbf u \times (\lambda\mathbf v)\\ \mathbf u \times (\mathbf v+\mathbf w)=\mathbf u \times \mathbf v+\mathbf u \times \mathbf w\\ (\mathbf u + \mathbf v)\times\mathbf w=\mathbf u \times \mathbf w+\mathbf v \times \mathbf w \end{align*}$$
(these properties are called "being bilinear").
There's one more property that we require in order for the cross product to make sense geometrically, which is that it shouldn't depend on which way you are looking at your problem. If we rotate our vectors and then take their cross product we should get the same answer as if we take their cross product and then rotate. Otherwise people looking at the same problem from different angles would get different answers! In symbols we represent the rotation by some orthogonal matrix $M$ with determinant $1$, and say that for each such matrix we want the following property:
$$M(\mathbf u\times\mathbf v)=(M\mathbf u)\times(M\mathbf v)$$
(This is called "invariance" or sometimes "covariance". Notice that the dot product also makes geometric sense in this way. If you rotate two vectors and then take their dot product you get the same answer as you would have gotten before the rotation. Or in other words $(M\mathbf u)\cdot(M\mathbf v)=u\cdot v$.)
Now here's the clever bit: I claim that the cross product is the only function with these properties. This explains why the cross product is interesting: it's the only form of multiplication that makes any sense at all! (Actually there are also the functions like $\lambda\mathbf u\times\mathbf v$ that are just scalings of the cross product by some factor, but each of these functions can be written in terms of any of the others and so we can just pick the usual cross product to be our favourite, and work in terms of that.)
Proof
I'll show that if we have a function, $\times$, that is bilinear and invariant (i.e. it obeys the four equations listed above) then it is in fact the cross product. We'll work in terms of the usual basis vectors $\mathbf i$, $\mathbf j$ and $\mathbf k$.
First we'll try to work out what $\mathbf i \times \mathbf j$ is. Let $M$ be the rotation of $180^\circ$ about the $k$-axis. Then $M\mathbf i=-\mathbf i$ and $M\mathbf j=-\mathbf j$. So we have
$$M(\mathbf i \times \mathbf j)=(M\mathbf i)\times(M\mathbf j)=(-\mathbf i)\times(-\mathbf j)=\mathbf i\times\mathbf j$$
(the last step used bilinearity to move the minus signs out so they could cancel). This means that $\mathbf i \times \mathbf j$ is fixed by the rotation $M$, and so must lie on the $\mathbf k$-axis. As I said before, we're going to allow ourselves to pick our favourite scaling, so since we know $\mathbf i \times \mathbf j$ is on the $\mathbf k$-axis we might as well assume that $\mathbf i \times \mathbf j=\mathbf k$.
There's a rotation (of $120^\circ$ about $\mathbf i+\mathbf j+\mathbf k$) that takes $\mathbf i$ to $\mathbf j$, $\mathbf j$ to $\mathbf k$, and $\mathbf k$ to $\mathbf i$. Applying invariance under this matrix to our equation $\mathbf i \times \mathbf j=\mathbf k$ gives us $\mathbf j \times \mathbf k=\mathbf i$. Applying it again gives $\mathbf k \times \mathbf i=\mathbf j$.
There's also a rotation (of $180^\circ$ about $\mathbf i+\mathbf j$) that takes $\mathbf i$ to $\mathbf j$, $\mathbf j$ to $\mathbf i$, and $\mathbf k$ to $-\mathbf k$. Applying invariance under this matrix to our equation $\mathbf i \times \mathbf j=\mathbf k$ gives us $\mathbf j \times \mathbf i=-\mathbf k$. Similarly we have $\mathbf k \times \mathbf j=-\mathbf i$ and $\mathbf i \times \mathbf k=-\mathbf j$.
Finally we want to know what $\mathbf i \times \mathbf i$ is. Let $M$ be the $180^\circ$ rotation about the $\mathbf k$-axis, as before. Then
$$M(\mathbf i \times \mathbf i)=(M\mathbf i)\times(M\mathbf i)=(-\mathbf i)\times(-\mathbf i)=\mathbf i\times\mathbf i$$
so $\mathbf i \times \mathbf i$ is fixed by $M$ and therefore lies on the $\mathbf k$-axis. But the same argument applied with a rotation about the $\mathbf j$-axis shows that $\mathbf i \times \mathbf i$ lies on the $\mathbf j$-axis too! These two axes only intersect at $\mathbf 0$. So $\mathbf i \times \mathbf i=\mathbf 0$ and by the same argument $\mathbf j \times \mathbf j=\mathbf 0$ and $\mathbf k \times \mathbf k=\mathbf 0$.
Now since we know how to cross product any two basis vectors we can calculate the cross product of any two vectors by multiplying out (using bilinearity):
$$\begin{align*}(u_i\mathbf i+u_j\mathbf j+u_k\mathbf k)\times(v_i\mathbf i+v_j\mathbf j+v_k\mathbf k)= &u_iv_i\mathbf i\times\mathbf i+u_iv_j\mathbf i\times\mathbf j+u_iv_k\mathbf i\times\mathbf k\\ +&u_jv_i\mathbf j\times\mathbf i+u_jv_j\mathbf j\times\mathbf j+u_jv_k\mathbf j\times\mathbf k\\ +&u_kv_i\mathbf k\times\mathbf i+u_kv_j\mathbf k\times\mathbf j+u_kv_k\mathbf k\times\mathbf k\\ =(u_jv_k-u_kv_j)\mathbf i+(u_kv_i-u_iv_k)&\mathbf j+(u_iv_j-u_jv_i)\mathbf k \end{align*}$$
This is the formula for the cross product.
The above was a rewrite of my original answer which said more or less the same thing as above but in more formal terms. I'll put my original answer here because I think some people reading this might like to see the technical details:
Given a $3$-dimensional oriented real inner-product space $V$ the group of symmetries preserving the inner-product and orientation is $\mathrm{SO}(V)$. The invariant tensors under $\mathrm{SO}(V)$ are $\delta_{ij}$, $\delta^{ij}$, and $\varepsilon_{ijk}$, along with the things they generate like $\delta^{ij}\varepsilon_{klm}$ and so on.
The tensors $\delta_{ij}$ and $\delta^{ij}$ are the inner-product and the inner-product induced on the dual space. These aren't very interesting because we defined $\mathrm{SO}(V)$ to preserve these, so we already knew that we were going to get them. But the tensor $\varepsilon_{ijk}$ is in some sense new. Therefore we are motivated to investigate $\varepsilon_{ijk}$ or equivalently the bilinear map $V\times V\rightarrow V$ given by $(v\times w)^i=\delta^{ij}\varepsilon_{jkl}v^kw^l$. This is the cross product.