I give a brief answer. More details, in particular on the conjugated variables used for the electromagnetic field see my post
The origin of quantization
(Note that I use the convention $\hbar=1$ and the angular frequency $\omega = 2\pi f$ instead of simple frequency $f$ in that post).
The formula
$$ U =\frac{1}{2}\left(\epsilon_0 E_0^2 + \frac{B_0^2}{\mu_0^2}\right)\quad\text{respectively}\quad H(\equiv E) = \frac{1}{2}\int dV \left(\epsilon_0 E_0^2 + \frac{B_0^2}{\mu_0^2}\right)$$
is the energy density (respectively energy $E$ or Hamiltonian $H$ ) in classical electrodynamics. In order to get $$E=hf$$ the electromagnetic field has to be quantized.
Actually, not every electromagnetic field can be quantized, in particular static fields cannot be quantized. The field should have some oscillatory behaviour, or in other words be a (superposition of) electromagnetic wave(s).
In order to achieve that one first has to find the conjugated variables the Hamiltonian depends on. Once one has the Hamiltonian written in conjugated variables $\mathbf{Q_k}$ and $\mathbf{P_k}$ we can recognize that the form of the Hamiltonian looks like a summation of Hamiltonians of the harmonic oscillator which distinguish by different wave vectors $\mathbf{k}$. Up to now it is still classical physics.
Now Quantum physics:
We can study harmonic oscillators as a quantum system, if we do this --- well documented in every textbook on Quantum Mechanics -- we learn that the energy spectrum of a quantum harmonic oscillator is no longer continous, but can take on only discrete values ($N_\mathbf{k}$ is an integer number which can be 0 (no excitation) or 1 (first excitation), 2, .... 100 ("high" excitation),.... etc.)
$$ E_\mathbf{k} = hf_\mathbf{k} \left(N_\mathbf{k}+\frac{1}{2}\right)$$
$\mathbf{k}$ is the parameter which distinguishes the different harmonic oscillators since we have many of them in the electromagnetic field.
Therefore the total energy of quantum electromagnetic field can written as
$$ E = \sum_{\mathbf{k}} h f_\mathbf{k} \left(N_k+\frac{1}{2}\right) \quad\tag{1}$$
The last step is a conceptual one. For study just select a single harmonic oscillator characterized by the wave vector $\mathbf{k}$ from the ensemble of oscillators which represent the electromagnetic field. (you can also choose another oscillator with wave vector $\mathbf{h}$ -- it does not change anything) The energy excitations of a single quantum harmonic oscillator (and this course also valid for the one chosen in the selection) with a fixed frequency $f_\mathbf{k}$ are like a ladder which constant energy steps. One single step in the excitation of the energy of the electromagnetic field is called a photon with the energy (step):
$$ E_k = h f_k $$
Of course one can go all steps backwards from (1). For getting something which looks like the integral over the energy density one would need to define field strength operators of $\mathbf{E}$ and $\mathbf{B}$. They can actually be constructed from $\mathbf{Q_k}$ and $\mathbf{P_k}$, but right now I don't know the exact formula since in modern descriptions of the quantum electromagnetic field they are very rarely used. The discrete harmonic oscillator summation over different wave vectors $\mathbf{k}$ can always be written as a continous integral with corresponding math tools.
At some point one would to deprive $\mathbf{E}$ and $\mathbf{B}$ of the operator status in order to get the classical formula back. But it is a bit weird, because once you are quantum one would stay quantum. But if the quantum numbers $N_\mathbf{k}$ are of high value one can kind of neglect that the concerned operators fulfill (non-zero) commutator relations.
EDIT
This transition from the quantum world back to classical physics is based on the Correspondence Principle which says that in case of high quantum numbers there is no sensible difference between classical physics and quantum physics anymore.
In case of high quantum numbers, i.e. in case one deals with a high number of photons one would be able to describe such a situation with classical electrodynamics, i.e. also with a classical energy density to be computed from the field strengths.