One nice way to construct the 24-cell and 120-cell is to take a tetrahedron or icosahedron, look its rotational symmetries (of order 12 or 60, respectively) in $\textrm{SO}(3)$, and pull them back along the double cover $\Phi : \textrm{SU}(2) \to \textrm{SO}(3)$, identifying $\textrm{SU}(2)$ with the unit quaternions: in this way we get 24 or 120 points lying on $S^3$.
Is there a nice way to verify that these points form the vertices of a regular 4-polytope? It seems plausible to me, and I've seen it stated, but I've never seen a proof. I believe attempting this with the symmetries of the cube give you 48 points which do not form a regular 4-polytope, so ideally an explanation would explain why that fails without being too computational.
If necessary, I'd be ok assuming as given that the 24-cell and 600-cell exist, if that could help in verifying this construction produces them.
The double cover and identification with unit quaternions boils down to associating a rotation in 3-space, around an axis with unit vector $(x, y, z)$ by angle $\theta$, with the quaternion $q = c + xs\mathbf{i} + ys\mathbf{j} + zs\mathbf{k}$, where $c = \cos(\theta/2)$ and $s = \sin(\theta/2)$, and with its negation $-q$. The quaternion $q$ in turn is identified with the point $(c, xs, ys, zs)$ in $\mathbb{R}^4$.
Since we are dealing with the rotation axes of a regular polytope, they'll be symmetrically distributed (and in a "full-dimensional" way, i.e., not like a prism, whose symmetries are essentially lower-dimensional.)
Thus the corresponding points in $\mathbb{R}^4$ will also be "evenly" distributed with a high degree of symmetry. Since the points arise from a group, which acts transitively on its own elements, and since the group action is multiplication by unit quaternions, which are isometries of 4-space, the points will be acted on transitively by isometries of $\mathbb{R}^4$.
So I don't think there's really a fundamental reason that the convex hull of these points has to be a regular 4-polytope, but it does have to be a highly symmetric vertex-transitive one, which suggests that it's pretty likely to be regular.
Tetrahedral symmetry
Let's look at a tetrahedron, which has threefold axes of rotation through each of its vertices, and twofold axes of rotation through the midpoints of each opposite pair of edges. We'll need to pick some coordinates for a regular tetrahedron. Let me take alternate vertices of the cube $[-1,1]\times[-1,1]\times[-1,1]$, say the ones with evenly many minus signs, and scale them to unit length: that gives me vertices $$\Bigl(\frac{1}{\sqrt{3}},\frac{1}{\sqrt{3}},\frac{1}{\sqrt{3}}\Bigr), \Bigl(\frac{-1}{\sqrt{3}},\frac{-1}{\sqrt{3}},\frac{1}{\sqrt{3}}\Bigr), \Bigl(\frac{-1}{\sqrt{3}},\frac{1}{\sqrt{3}},\frac{-1}{\sqrt{3}}\Bigr), \Bigl(\frac{1}{\sqrt{3}},\frac{-1}{\sqrt{3}},\frac{-1}{\sqrt{3}}\Bigr).$$
Our rotations are:
The 16 vertices $\bigl(\frac{\pm 1}{2},\frac{\pm 1}{2},\frac{\pm 1}{2},\frac{\pm 1}{2}\bigr)$ taken by themselves form a tesseract. The remaining eight vertices create a pyramid raised above each facet of the tesseract. Each 2-face (square) $F$ is in two facets of the tesseract (cubes); the pyramids raised over those two facets each contain a square pyramid over $F$. It just so happens that these two square pyramids are actually "coplanar", that is, contained in the same 3-dimensional subspace, so that they join together in an octahedron cell, and $F$ is absorbed in its interior.
It might help to visualize an analogous process in 3 dimensions. Add a vertex above each face of a cube, forming pyramids above the faces. This will give you 24 triangular faces, but if the new vertices are at just the right height, the two triangles over each of the original edges will be coplanar, forming a rhombus, and the original edge will be absorbed in its interior. This gives you the rhombic dodecahedron:
That picture (from MathWorld) shows the skeleton of the rhombic dodecahedron balanced on the edges of the cube.
Since there's one of these octahedra for each 2-face of the tesseract, we get 24 octahedral cells. This is actually enough to know that it's a regular 4-polytope, given Theorem 4C4 of Peter McMullen's thesis:
Octahedral symmetry
If we consider a 3-cube with the unit-norm vertices $\bigl(\frac{\pm 1}{\sqrt{3}},\frac{\pm 1}{\sqrt{3}},\frac{\pm 1}{\sqrt{3}}\bigr)$, it has the same twelve rotations that we just discussed as symmetries, along with some more. There's a fourfold axis of rotation through the midpoints of each pair of opposite faces. These are the same axes $(0,0,1)$, $(0,1,0)$, $(1,0,0)$ that we saw for the tetrahedron, but now we can include rotations of $\pi/2$ and $3\pi/2$ as well as the rotation by $\pi$ that we had. This gives us six additional rotations, so twelve new points in $\mathbb{R}^4$:
$$ \Bigl(\frac{1}{\sqrt{2}},0,0,\frac{1}{\sqrt{2}}\Bigr),\quad \Bigl(\frac{-1}{\sqrt{2}},0,0,\frac{-1}{\sqrt{2}}\Bigr),\quad \Bigl(\frac{-1}{\sqrt{2}},0,0,\frac{1}{\sqrt{2}}\Bigr),\quad \Bigl(\frac{1}{\sqrt{2}},0,0,\frac{-1}{\sqrt{2}}\Bigr),\\ \Bigl(\frac{1}{\sqrt{2}},0,\frac{1}{\sqrt{2}},0\Bigr),\quad \Bigl(\frac{-1}{\sqrt{2}},0,\frac{-1}{\sqrt{2}},0\Bigr),\quad \Bigl(\frac{-1}{\sqrt{2}},0,\frac{1}{\sqrt{2}},0\Bigr),\quad \Bigl(\frac{1}{\sqrt{2}},0,\frac{-1}{\sqrt{2}},0\Bigr),\\ \Bigl(\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}},0,0\Bigr),\quad \Bigl(\frac{-1}{\sqrt{2}},\frac{-1}{\sqrt{2}},0,0\Bigr),\quad \Bigl(\frac{-1}{\sqrt{2}},\frac{1}{\sqrt{2}},0,0\Bigr),\quad \Bigl(\frac{1}{\sqrt{2}},\frac{-1}{\sqrt{2}},0,0\Bigr). $$
There are also twofold axes of rotation through the midpoints of each pair of opposite edges; there are six such pairs, giving us twelve more points in $\mathbb{R}^4$. The unit vectors are $\bigl(\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}},0\bigr)$, $\bigl(\frac{-1}{\sqrt{2}},\frac{1}{\sqrt{2}},0\bigr)$, $\bigl(\frac{1}{\sqrt{2}},0,\frac{1}{\sqrt{2}}\bigr)$, $\bigl(\frac{-1}{\sqrt{2}},0,\frac{1}{\sqrt{2}}\bigr)$, $\bigl(0, \frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}}\bigr)$, and $\bigl(0, \frac{-1}{\sqrt{2}},\frac{1}{\sqrt{2}}\bigr)$. The quaternions are $$ \Bigl(0,\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}},0\Bigr),\quad \Bigl(0,\frac{-1}{\sqrt{2}},\frac{-1}{\sqrt{2}},0\Bigr),\quad \Bigl(0,\frac{-1}{\sqrt{2}},\frac{1}{\sqrt{2}},0\Bigr),\quad \Bigl(0,\frac{1}{\sqrt{2}},\frac{-1}{\sqrt{2}},0\Bigr),\\ \Bigl(0,\frac{1}{\sqrt{2}},0,\frac{1}{\sqrt{2}}\Bigr),\quad \Bigl(0,\frac{-1}{\sqrt{2}},0,\frac{-1}{\sqrt{2}}\Bigr),\quad \Bigl(0,\frac{-1}{\sqrt{2}},0,\frac{1}{\sqrt{2}}\Bigr),\quad \Bigl(0,\frac{1}{\sqrt{2}},0,\frac{-1}{\sqrt{2}}\Bigr),\\ \Bigl(0,0,\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}}\Bigr),\quad \Bigl(0,0,\frac{-1}{\sqrt{2}},\frac{-1}{\sqrt{2}}\Bigr),\quad \Bigl(0,0,\frac{-1}{\sqrt{2}},\frac{1}{\sqrt{2}}\Bigr),\quad \Bigl(0,0,\frac{1}{\sqrt{2}},\frac{-1}{\sqrt{2}}\Bigr). $$
This gives us 48 points total.
If we consider only the 24 new points, we can identify the facet-defining hyperplanes: the sixteen variations of $\pm x_1 \pm x_2 \pm x_3 \pm x_4 = \sqrt{2}$, and the eight versions of $x_i = \frac{\pm1}{\sqrt{2}}$.
Each of these hyperplanes contains 6 of the new points. For instance, $x_1 + x_2 + x_3 + x_4 = \sqrt{2}$ matches $\bigl(\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}},0,0\bigr)$, $\bigl(\frac{1}{\sqrt{2}},0,\frac{1}{\sqrt{2}},0\bigr)$, $\bigl(\frac{1}{\sqrt{2}},0,0,\frac{1}{\sqrt{2}}\bigr)$, $\bigl(0,\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}},0\bigr)$, $\bigl(0,\frac{1}{\sqrt{2}},0,\frac{1}{\sqrt{2}}\bigr)$, and $\bigl(0,0,\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}}\bigr)$. These six points form an octahedron.
Thus the 24 new points form another 24-cell, rotated with respect to the first one; in fact each vertex of the first one is above the midpoint of a facet of the new one, e.g. $\bigl(\frac{1}{2},\frac{1}{2},\frac{1}{2},\frac{1}{2}\bigr)$ is above the octahedron we just considered.
It's not as easy to visualize the structure of the convex hull of the union of these two 24-cells. We can find a facet-defining hyperplane, though, by verifying that all 48 points satisfy $x_1 + x_2 + \sqrt{2}x_3 \leq 1 + 1/\sqrt{2}$. The four points meeting the equality condition are $\bigl(\frac{1}{2},\frac{1}{2},\frac{1}{2},\frac{\pm 1}{2}\bigr)$, $\bigl(\frac{1}{\sqrt{2}},0,\frac{1}{\sqrt{2}},0\bigr)$, and $\bigl(0,\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}},0\bigr)$.
These form a disphenoid tetrahedron; two of the edges have length 1 and the others have length $\sqrt{2-\sqrt{2}}$.
So we know this 4-polytope is not regular. (In fact, it is the Disphenoidal 288-cell.)
Icosahedral symmetry
Considering an icosahedron, its twelve vertices are in six opposite pairs, with an axis of five-fold symmetry through each pair: this gives us 24 rotations. Its twenty faces are in ten opposite pairs, with an axis of three-fold symmetry through each pair: this gives us 20 more rotations. Its thirty edges are in fifteen opposite pairs, with an axis of two-fold symmetry through each pair: this gives us 15 more rotations. Together with the identity, we have 24 + 20 + 15 + 1 = 60 rotations, and so 120 corresponding points in $\mathbb{R}^4$.
The identity rotation gives us the two points $(1,0,0,0)$ and $(-1,0,0,0)$ in $\mathbb{R}^4$.
We can take unit vectors for the six axes of five-fold symmetry to be $\frac{1}{\sqrt{10}}(0, \pm \sqrt{5-\sqrt{5}}, \sqrt{5+\sqrt{5}})$, and the rotations of these coordinates. Let $\gamma$ be $\frac{1}{4}(\sqrt{5} + 1)$ (half the golden ratio) and $\delta$ be $\frac{1}{4}(\sqrt{5}-1)$ (half the golden ratio's inverse). Using the facts that $\cos(\pi/5) = \gamma$, $\sin(\pi/5) = \frac{1}{4}\sqrt{10-2\sqrt{5}}$, $\cos(2\pi/5) = \delta$, and $\sin(2\pi/5) = \frac{1}{4}\sqrt{10+2\sqrt{5}}$, we get the points
$$ \Bigl(\gamma, 0, \delta, \frac{1}{2}\Bigr), \Bigl(\delta, 0, \frac{1}{2}, \gamma\Bigr), \Bigl(-\delta, 0, \frac{1}{2}, \gamma\Bigr), \Bigl(-\gamma, 0, \delta, \frac{1}{2}\Bigr),\\ \Bigl(\gamma, 0, -\delta, \frac{1}{2}\Bigr), \Bigl(\delta, 0, \frac{-1}{2}, \gamma\Bigr), \Bigl(-\delta, 0, \frac{-1}{2}, \gamma\Bigr), \Bigl(-\gamma, 0, -\delta, \frac{1}{2}\Bigr),\\ \Bigl(\gamma, \delta, \frac{1}{2}, 0\Bigr), \Bigl(\delta, \frac{1}{2}, \gamma, 0\Bigr), \Bigl(-\delta, \frac{1}{2}, \gamma, 0\Bigr), \Bigl(-\gamma, \delta, \frac{1}{2}, 0\Bigr),\\ \Bigl(\gamma, -\delta, \frac{1}{2}, 0\Bigr), \Bigl(\delta, \frac{-1}{2}, \gamma, 0\Bigr), \Bigl(-\delta, \frac{-1}{2}, \gamma, 0\Bigr), \Bigl(-\gamma, -\delta, \frac{1}{2}, 0\Bigr),\\ \Bigl(\gamma, \frac{1}{2}, 0, \delta\Bigr), \Bigl(\delta, \gamma, 0, \frac{1}{2}\Bigr), \Bigl(-\delta, \gamma, 0, \frac{1}{2}\Bigr), \Bigl(-\gamma, \frac{1}{2}, 0, \delta\Bigr),\\ \Bigl(\gamma, \frac{1}{2}, 0, -\delta\Bigr), \Bigl(\delta, \gamma, 0, \frac{-1}{2}\Bigr), \Bigl(-\delta, \gamma, 0, \frac{-1}{2}\Bigr), \Bigl(-\gamma, \frac{1}{2}, 0, -\delta\Bigr) $$ and all their negations.
We can take unit vectors for the ten axes of three-fold symmetry to be $\frac{1}{\sqrt{3}}(1,1,1)$, $\frac{1}{\sqrt{3}}(-1,1,1)$ and its permutations, and $(0, \pm \frac{\sqrt{5}+1}{2\sqrt{3}}, \frac{\sqrt{5}-1}{2\sqrt{3}})$ and the rotations of its coordinates. Rotating about these through $2\pi/3$ or $4\pi/3$ corresponds to the points
$$ \Bigl(\frac{1}{2}, \frac{1}{2}, \frac{1}{2}, \frac{1}{2} \Bigr), \Bigl(\frac{-1}{2}, \frac{1}{2}, \frac{1}{2}, \frac{1}{2} \Bigr),\\ \Bigl(\frac{1}{2}, \frac{-1}{2}, \frac{1}{2}, \frac{1}{2} \Bigr), \Bigl(\frac{-1}{2}, \frac{-1}{2}, \frac{1}{2}, \frac{1}{2} \Bigr),\\ \Bigl(\frac{1}{2}, \frac{1}{2}, \frac{-1}{2}, \frac{1}{2} \Bigr), \Bigl(\frac{-1}{2}, \frac{1}{2}, \frac{-1}{2}, \frac{1}{2} \Bigr),\\ \Bigl(\frac{1}{2}, \frac{1}{2}, \frac{1}{2}, \frac{-1}{2} \Bigr), \Bigl(\frac{-1}{2}, \frac{1}{2}, \frac{1}{2}, \frac{-1}{2} \Bigr),\\ \Bigl(\frac{1}{2}, 0, \gamma, \delta\Bigr), \Bigl(\frac{-1}{2}, 0, \gamma, \delta \Bigr),\\ \Bigl(\frac{1}{2}, 0, -\gamma, \delta\Bigr), \Bigl(\frac{-1}{2}, 0, -\gamma, \delta \Bigr),\\ \Bigl(\frac{1}{2}, \gamma, \delta, 0 \Bigr), \Bigl(\frac{-1}{2}, \gamma, \delta, 0 \Bigr),\\ \Bigl(\frac{1}{2}, -\gamma, \delta, 0 \Bigr), \Bigl(\frac{-1}{2}, -\gamma, \delta, 0 \Bigr),\\ \Bigl(\frac{1}{2}, \delta, 0, \gamma \Bigr), \Bigl(\frac{-1}{2}, \delta, 0, \gamma \Bigr),\\ \Bigl(\frac{1}{2}, \delta, 0, -\gamma \Bigr), \Bigl(\frac{-1}{2}, \delta, 0, -\gamma \Bigr) $$ and all their negations.
With the coordinates I chose for the icosahedron's vertices (seen above as the axes of five-fold rotation), the edge midpoints end up being $\frac{1}{\sqrt{10}}(0,0,\sqrt{5+\sqrt{5}})$, $\frac{1}{2\sqrt{10}}(\sqrt{5-\sqrt{5}},\sqrt{10+4\sqrt{5}},\sqrt{5+\sqrt{5}})$, and the rotations of these coordinates (with all assignments of +/- signs.) After scaling to unit vectors, and choosing one from each opposite pair, we get $(1,0,0)$, $(0,1,0)$, $(0,0,1)$ and all the rotations of $(\delta, \gamma, \frac{1}{2})$ with at most one minus sign. Rotating by an angle of $\pi$ about these vectors gives us the quaternions $$ (0,1,0,0), \quad (0,0,1,0), \quad (0, 0, 0, 1), \\ (0, \delta, \gamma, \frac{1}{2}), \quad (0, -\delta, \gamma, \frac{1}{2}), \quad (0, \delta, -\gamma, \frac{1}{2}), \quad (0, \delta, \gamma, -\frac{1}{2}), \\ (0, \gamma, \frac{1}{2}, \delta), \quad (0, -\gamma, \frac{1}{2}, \delta), \quad (0, \gamma, -\frac{1}{2}, \delta), \quad (0, \gamma, \frac{1}{2}, -\delta), \\ (0, \frac{1}{2}, \delta, \gamma), \quad (0, -\frac{1}{2}, \delta, \gamma), \quad (0, \frac{1}{2}, -\delta, \gamma), \quad (0, \frac{1}{2}, \delta, -\gamma) $$ and all their negations.
So after going through all the rotations of the icosahedron, we can summarize all the corresponding points in $\mathbb{R}^4$ the same way that Wikipedia puts it: the eight permutations of $(\pm 1, 0, 0, 0)$, the sixteen points $\bigl(\frac{\pm 1}{2},\frac{\pm 1}{2},\frac{\pm 1}{2},\frac{\pm 1}{2}\bigr)$, and the 96 even permutations of $\bigl(0,\frac{\pm 1}{2},\pm \delta, \pm \gamma\bigr)$.
So what is the convex hull of these points?
We can find a facet-defining hyperplane by verifying that all the points satisfy $x_1 + x_2 \leq \gamma + \frac{1}{2}$, and exactly four of the points are in that hyperplane: $$ \Bigl(\frac{1}{2}, \gamma, \delta, 0\Bigr), \Bigl(\frac{1}{2}, \gamma, -\delta, 0\Bigr), \Bigl(\gamma, \frac{1}{2}, 0, \delta\Bigr), \Bigl(\gamma, \frac{1}{2}, 0, -\delta\Bigr). $$
Since these points are all at a distance of $2\delta$ from eachother, they form a regular tetrahedron. If we trust that all the facets are regular tetrahedra (the vertex-transitivity will help a lot if you want to rigorously check this), then McMullen's Theorem 4C4 tells us this is a regular 4-polytope.