We define the irradiance as the average density flux arriving at a surface with units $$\frac{W}{m^2}$$. So for a point light source, we have: $$E = \frac{\Phi}{4 \pi r^2}$$ since the area of a sphere is $$4 \pi r^2$$. Where $$\Phi$$ is the flux or power.

A (to me) similar concept is intensity which is the amount of power per angle. Again, for a sphere with a point light at the center, this is $$I = \frac{\Phi}{4 \pi r^2}$$ with the unit $$[\frac{W}{sr}]$$ (watt over steradian)

Now, the book defines radiance for a point $$p$$ as $$L = \frac{d\Phi}{d\omega dA^\perp}$$ in units $$[\frac{W}{sr\cdot m^2}]$$. Here, $$\omega$$ is the direction where the light comes from, $$A^\perp$$ is the projected are of $$A$$ as seen here: This means that practically, when I implement a point light source with a given power that shines at a point $$p$$, I need to do the following to arrive at radiance:

• Divide by $$4 \pi r^2$$ to convert power into $$[\frac{W}{sr}]$$, or in other words, intensity.
• Given intensity, I need to divide it by $$4 \pi r^2$$ and multiply by $$\cos \theta$$ to arrive at $$[\frac{W}{sr\cdot m^2}]$$, the final radiance. The multiplication by $$\cos \theta$$ is to project $$A$$ to $$A^\perp$$ and is the dot product of the surface normal $$n$$ with the direction $$w$$ (as both are normalized).

For both calculations, $$r$$ is the distance between the light source and my point $$p$$.

However, when I look at the source, this is not what happens. The point light returns intensity divided by $$r^2$$ as seen here:

return I / DistanceSquared(pLight, ref.p);


and the integrator then multiplies it with the dot product (and the brdf) in the whitted integrator

L += f * Li * AbsDot(wi, n) / pdf;


So what is wrong in my derivation? Why do we "only" divide once by $$4\pi r^2$$ (to get the Intensity I) and not twice? Aren't we missing either the power per area or the power per steradian?

Your definition for radiant intensity is wrong: it should be just $$\Phi / 4\pi$$. There are only $$4\pi$$ steradians in a sphere no matter how big it is, so $$r$$ doesn't come into it.
Also note that you can't calculate radiance for a point source—it would be infinite, due to the fact the point source emits a finite amount of flux compressed into zero size. It subtends zero solid angle, from the receiver's point of view, so the $$d\omega$$ factor in the denominator of radiance would be zero. The usable quantities are radiant intensity in a certain direction from the point source (which could vary with direction, for a non-omnidirectional light), and irradiance in a certain direction and distance. As you've seen, irradiance is obtained from radiant intensity by dividing by $$r^2$$—you could think of this as "area per steradian", as it's the conversion factor from $$4\pi$$ steradians to $$4\pi r^2$$ area of a sphere; then the units work out.
Where radiance would actually show up is when dealing with an area light rather than a point light. Then you would have the flux being distributed over a finite solid angle from the receiver's point of view, and you'd get the incident irradiance by integrating the light's radiance over that solid angle. The emitted radiance would be $$\Phi/(2\pi A_\text{light})$$, assuming it's emitted uniformly over the surface of the light and into all directions from each point. This is only $$2\pi$$ since it's only emitting into the outward-facing hemisphere. Also note that you don't do any division by $$r^2$$ for radiance—with area lights, the distance attenuation comes naturally as a result of the light subtending less solid angle from the receiver's point of view, when the receiver is farther away.