The issue in the derivation from thermal QED is to evaluate the partition sum of the free electromagnetic field,
Z=\mathrm{Tr} \exp(-\beta \hat{H}).
The most simple way is to use the mode decomposition of the field operator in the totally radiation-gauge fixed formulation, which gives
\hat{H}=\sum_{\vec{p},\lambda} \omega(\vec{p}) \hat{N}(\vec{p},\lambda),
where \vec{p} runs through the values \vec{p} \in \frac{2 \pi}{L} \mathbb{Z}, where L is the length of a large cube within the thermal radiation field. For simplicity I've assumed periodic boundary conditions. The specific form of the boundary conditions is pretty unimportant in the thermodynamic limit, which we have to take at the end of the calculation. \omega(\vec{p})=|\vec{p}| is the energy of a photon of momentum \vec{p} and \lambda \in \{-1,1 \} are the two helicity values (polarizations) of a photon.
To get the Planck law for the occupation number distribution we evaluate the somewhat more general expression
Z[J]=\mathrm{Tr} \exp(-\sum_{\vec{p},\lambda} j(\vec{p},\lambda) \omega(\vec{p}) \hat{N}(\vec{p},\lambda).
To evaluate the trace we use the occupation-number eigenstates |\{N(\vec{p},\lambda) \}, where N(\vec{p},\lambda) \in \{0,1,2,\ldots\}. This gives
Z[j]=\prod_{\vec{p},\lambda} \sum_{N(\vec{p},\lambda)=0}^{\infty} \exp[-j(\vec{p},\lambda) \omega(\vec{p}) N(\vec{p},\lambda)] = \prod_{\vec{p},\lambda} \frac{1}{1-\exp[-j(\vec{p},\lambda) \omega(\vec{p})]}.
To proceed further we evaluate the logarithm, and convert the sum to an integral by taking the large-volume limit. Then in each volume in momentum space we have \frac{L^3}{(2 \pi)^3} \mathrm{d}^3 \vec{p} states:
\ln Z[j]=-V \sum_{\lambda} \int_{\mathbb{R}^3} \frac{\mathrm{d}^3<br />
\vec{p}}{(2 \pi)^3} \ln \left<br />
\{1-\exp[-j(\vec{p},\lambda)\omega(\vec{p})] \right \}.
The energy density of the radiation field is then given by
\mathrm{d} U=-\mathrm{d}^3 \vec{p} \frac{1}{V} \sum_{\lambda} \left ( \frac{\delta Z}{\delta j(\vec{p},\lambda)} \right )_{j(\vec{p},\lambda)=\beta}.
This gives
\mathrm{d} U=\frac{1}{4 \pi^3} \mathrm{d}^3 \vec{p} \frac{\omega(\vec{p})}{\exp[\beta \omega(\vec{p})]-1}.
Further we have
\vec{p}^2=\omega^2 \; \Rightarrow \; \mathrm{d}^3 \vec{p} = \mathrm{d}^2 \Omega \omega^2 \mathrm{d} \omega,
and thus finally the Planck radiation formula
\frac{\mathrm{d} U}{\mathrm{d} \omega}=\frac{\omega^3}{\pi^2} \frac{1}{\exp(\beta \omega)-1}.