Consider having ten equally sized bins for scalar values in the range [0, 10]. There are very intuitive and straightforward methods for assigning a scalar to its corresponding bin e.g. 1.7 belongs to bin 2 and 8.4 belongs to bin 9.
Are there any intuitive methods for the discretization of unit quaternions, such that I can retrieve a distinct bin for any quaternion?
It depends on just how efficient you want to be. I'll describe the method I would use for $S^2$; the approach for $S^3$ (or its quotient by inversion, $\mathbb{R}P^3$, which is the actual space of the quaternions) is similar.
Intuitively, think of a cube circumscribed around $S^2$; then the faces of this cube (or rather, their projections onto $S^2$ via lines through the origin) divide the sphere into six corresponding faces. You can then look at the 2d coordinates on each of these faces and bin into a grid based on those coordinates.
In practice, this means finding the coordinate of your unit vector (or quaternion) that's largest in absolute value; then binning based on (a) which coordinate is largest, (b) whether it's positive or negative (these two together represent the six possible 'faces' of the cube), and then (c) the bin it falls in based on the remaining two coordinates. Note that since those two coordinates are both smaller than the largest, their maximum possible value is $\frac{\sqrt{2}}{2}$; so multiplying both of them by $\sqrt{2}$ will give values in the range $[-1,1]$ for the binning.
The approach for $S^3$ is similar, though here since you're working with quaternions specifically you may be able to make the choice that your largest coordinate is always positive (by multiplying the rest of the quaternion by $-1$ as necessary); this reduces the number of possible projected faces from the 8 it would be down to 4. (Note that this has implications for things like animation, because while $q$ and $-q$ represent the same rotation, interpolating between $p$ and $q$ won't give the same result as interpolating between $p$ and $-q$.) Your bins will then be 3-dimensional based on the remaining three coordinates, and the exact same bound — that all of them are less than or equal to $\frac{\sqrt{2}}2$ in absolute value — holds.