The way to find the stars in some region is to loop through the stars and for each one evaluate a formula and do a test to see if it is in the region or not. There's not much point in using a cube first; on modern computers multiplications are cheaper than "if" statements!
So the important step is to find the formula for the shapes you are interested in.
The interior of a sphere centred on the origin with radius $r$ satisfies this inequality:
$$x^2+y^2+z^2 \leq r^2$$
This is just Pythagoras's theorem, applied to find the distance from the origin to the point $(x, y, z)$.
If you want the sphere centred elsewhere or with a different radius, you can just substitute for $x$, $y$, $z$:
$$x = x' - x_0$$
$$y = y' - y_0$$
$$z = z' - z_0$$
This gives:
$$(x'-x_0)^2 + (y'-y_0)^2 + (z'-z_0)^2 \leq r^2$$
In this formula, $(x_0, y_0, z_0)$ is the centre of the sphere.
The formula for the interior of a double cone with its point at the origin and its axis along the z-axis and a half-angle of 45 degrees is:
$$x^2 + y^2 \leq z^2$$
Again, this is just Pythagoras's theorem. The left-hand side is the squared distance from $(x, y, z)$ to the point on the z-axis at $(0, 0, z)$.
Note that this gives a double cone; if we want just one of the cones we need in addition something like:
$$z \geq 0$$
You can change the half-angle of the cone to $\phi$ by scaling $z$:
$$x^2 + y^2 \leq (z \tan(\phi))^2$$
If you want to change the axis of the cone, you can do it as follows. First, rewrite the equation as follows:
$$x^2 + y^2 + z^2 \leq z^2 (\tan^2(\phi) + 1)$$
Now the left-hand side looks like the formula for a sphere, and does not depend on the orientation of the cone. The right-hand side does. We need to replace $z$ on the right-hand side with a quantity $d$ that measures the distance along the desired axis of the cone.
Suppose you want the cone to point towards $(u, v, w)$. Then a measure of the distance in that direction is the dot-product:
$$d = \frac{u}{l}x + \frac{v}{l}y + \frac{w}{l}z$$
Where $l$ is the length of the vector $(u, v, w)$, which can be computed (again) using Pythagoras's theorem:
$$l = \sqrt{u^2 + v^2 + w^2}$$
So, putting this all together:
$$(x^2 + y^2 + z^2)(u^2 + v^2 + w^2) \leq (wx + vy + wz)^2 (\tan^2(\phi) + 1)$$
$$ux + vy + wz \geq 0$$
Finally, you can move the point of the cone around using a coordinate transform, just as for the sphere.
The "someone" you're describing was the Irish physicist and mathematician Willian Rowan Hamilton.
What he wanted was exactly a 3-dimensional number system -- not particularly for graphing as such, but more generally to get something that would make it possible to use algebraic techniques in space geometry as easy as the complex plane had by then already made it for plane geometry.
In modern terms one might say that what Hamilton was really looking for was what we know as three-dimensional vectors, but in those days there was a general feeling that you needed to have a well-behaved multiplication rule for your system before you were "allowed to" use algebraic notation for your calculations with it. The modern concept of a vector space (where we're quite comfortable with having things that can be added, but don't have a multiplication that satisfies the same rules as multiplication of real numbers) did not exist yet.
Hamilton later wrote that he had tried for a long time to find a multiplication rule that would work for three-dimensional quantities, but without luck (later it was proved to be impossible). Then in a sudden flash of inspiration he realized that a four-dimensional system would be possible. There's a plaque at the exact spot in Dublin where he said this insight occurred.
It is therefore not so much a matter of "deriving" the quaternions, as of thinking about the problem for a long time (and by that time Hamilton no doubt had a keen intuitive experience with the ways a multiplication rule can fail to work) and then suddenly noticing that this particular rule happens to work. Once it works, it doesn't need to have a neat story of how you found it -- though Frobenius later proved that the quaternions are the only finite-dimensional associative division algebra over $\mathbb R$ other than the real and complex numbers, so sooner or later someone would have found it.
Best Answer
In 3D space, initial coordinates are decided by the given Left-hand and Right-hand coordinates, sometimes called x-up, y-up, z-up conventions which can also be found (for example) in Wiki's
EulerAnglesitem.Initial Quaternions are defined corresponding to the intial x-y-z
For example, x-y-z-right hand coordinates, then 0,0,0,1 would be the initial quaternion.
Of course, the order also depends on the setting of program.
For example maybe DirectX and openGL are in different orders of the initial quaternion of the same coordinates.
And also your self programs could in one arbitrary initial order.
However, when one initial order was defined(corresponding to x-y-z), then the other 23 orders could also be well defined.
24 is the number of kinds of orders. When I decide which order was used in one 3D game engine which I cannot know the order by specification, I do experiments.
Maybe you could find a Matrix to Quaternion program, and then do rotations in Matrix, then convert the combined final matrix to Quaternions will lead you a good comprehension.