The azimuthal quantum number determines the shape of electron distribution around a nucleus ($s$ orbital has spherical distribution, $p$ orbital has dumbell-like distribution and so on).
But what determines the shape of those orbitals? Why doesn't the $s$ orbital look like a rectangle? Or a $p$ orbital look like a mobius strip?
The shapes of the orbitals are not random at all! They are the solutions to a particular differential equation called the hydrogen atom Schrödinger equation, \begin{gather*} - \frac{\hbar^2}{2m}\nabla^2 \psi + \frac{e^2}{4\pi \epsilon_0 r}\psi = E \psi \end{gather*} in a particular coordinate system called spherical coordinates. It turns out that the mathematics of this differential equation along with the condition of normalizability of the solutions $\psi$ restricts us to certain solutions for the energies $E$ and their corresponding wavefunctions $\psi$. These solutions get labeled by integers obeying the equations, \begin{align*} n &\in \mathbb{N} &&(n=1,2,3,\ldots) \\ l &\in [0,n-1] \subset \mathbb{Z} &&(l=0,1,2,\ldots, n-1)\\ m &\in [-l,l] \subset \mathbb{Z} &&(m=-l,-l+1,\ldots,0,\ldots,l-1,l) \end{align*} The wave functions themselves, which are related to the shapes we call "orbitals," have the extremely complicated form, \begin{gather*} \psi_{n,l,m} (r,\theta, \phi) \propto e^{-\frac{r}{na_0}} \left( \frac{2r}{na_0} \right)^l L_{n-l-1} ^{(2l+1)} \left( \frac{2r}{na_0} \right) P_l ^{m} (\cos(\theta)) e^{im\phi} \end{gather*} where the indexed functions $P_l^m$ and $L_{n-l-1} ^{(2l+1)}$ are examples of orthogonal polynomials called associated Legendre and associated Laguerre polynomials. As noted in the comments, the angular part of these functions get the special name of "spherical harmonic" functions. Notice that these functions are complex for every value of $m\neq 0$ as we have written them. If we take special combinations of these functions they will retain all the desirable properties but become completely real, as explained in a recent education article of mine. These are the orbital shapes that are drawn in chemistry books. The shapes of the orbitals are predominantly determined by the value of $l$, with $l=0$ orbitals having spherical symmetry, $l=1$ orbitals having the shapes of dumbbells, and $l=2$ orbitals mostly being shaped like clovers and one the shape of a dumbbell surrounded by a doughnut, and so forth. The value of $m$ mostly determines the "orientation" of the orbital. For instance, $l=1,m=\pm1$ corresponds to a dumbbell shaped orbital oriented along the $x$ or $y$ axis while $l=1,m=0$ is a dumbbell shaped orbital oriented along the $z$ axis. The quantum numbers are actually intimately related to the nodal structure of the wavefunctions, which means that the general shape can be predicted a priori, as my article explains.
The other two answers already here are correct and great. I'm going to build off the idea they proposed that the angular distribution is governed by spherical harmonics. But let me try to give a more abstract, but potentially more intuitive explanation.
There is a theorem from applying group theory to quantum mechanics that if the "Hamiltonian" of a system is symmetric under some group, then there must be an eigenbasis of the Hamiltonian which transforms under a representation of that symmetry.
Given the nature of the question, it's probably safe to assume that means nothing to the original asker, so I'll try to remove the jargon and explain it more simply. When an electron interacts with a proton to form Hydrogen, to a good approximation the electron feels the same force pointing toward the proton no matter which direction the electron is with respect to the proton. This means, for one thing, that the groundstate (the lowest energy state that atoms sit in most of the time at room temperature) of hydrogen is spherically symmetric, or that the electron has the same probability of being in any direction from the nucleus. Why should it be more likely to be above the proton if it feels no different force above the proton than below it?
But what about the excited states which are not spherically symmetric? Well it turns out that if they are not spherically symmetric, they must also be degenerate, which is a fancy way of saying that there is another quantum energy level with the same energy. And actually if you take all those degenerate states, and add up the probability distributions they represent, you once again have a spherically symmetric distribution. Also, this set of degenerate states, if you rotate your perspective, they transform into one another. For example, if I'm in a $L=1$, $m_L=+1$ state, and then I rotate my perspective by 180 degrees, now I must be in a $L=1$, $m_L=-1$ state. This is what I meant by "they must transform under a representation of the symmetry." In this case, the symmetry is rotations - the system is symmetric under rotations in the sense that if I rotate the system, it remains the same.
The spherical harmonics others have mentioned are all the different ways electron wavefunctions can transform under rotations.
There is a mathematical concept called the spherical harmonic that is useful in a variety of physics contexts. It turns out that it shows up when describing QM orbital shapes as well as, for example, three-dimensional electromagnetic radiation patterns from shortwave radio antennas. @Matt Hanson's answer has all the marvelous derivations and details of how this works!
I'm hardly a physicist, but this stuff has fascinated me all my life, and I'm now trying to help my daughter understand it. It's finally made sense. So here goes:
Watch FloatHeadPhysics video here (mentioned above) to understand how a waveform is shaped around nodes, at which the wave is "clamped" at zero. These nodes can be a point, as in a vibrating string, or for three dimensional orbitals, a flat plane, or a sphere, or ... a cone.
Now remember that an orbital uses spherical harmonics. Don't think XYZ. Think latitude and longitude. These angular nodes can be along a fixed latitude or longitude. (Plus of course the spherical nodes as we go to higher levels of n)
An s orbital has none of these angular nodes. The waveform is spherical.
A p orbital has one. It's convenient to think of it as following the X, Y or Z planes, but really it's following longitude and latitude. m=0 has a plane around the equator. m = -1 and 1 have a plane around the poles, and these are out of phase - consider one starting at longitude zero and the other at 90 degrees.
A d orbital has two. m=0 has both on a fixed latitude. Spread them evenly (I'm not sure of the exact angle) between the poles, and rotate around. You end up with two cones, one facing north and the other south. The gaps between thus create the pecular shape of the dumbell and doughnut. m= -1 and 1 have one node on latitude and one on longitude. The single longitude node just sits half way, on the equator. The latitude node cuts through the poles, and like the p orbitals can be in or out of phase. Hence you end up with the diagonal dumbells along XZ and YZ. For m=-2 and 2, both nodes are on the longitude. Thus you end up with the familiar XY orbital, plus another that's out of phase (by 45 degrees of the full sphere this time). Actually, I believe it's the XY that's out of phase.
f orbitals - same. m=0 has three latitude nodes and ends up with two doughnuts. m=-1 and 1 has two latitude and one longitude. m=-2 and 2 has one latitude and two longitude. m=-3 and 3 has three around the longitude. And again, the longitude nodes can be in or out of phase.
Now take a look at the orbital shapes again, but this time thinking about cutting them by latitude and longitude. Or both.
