Think of the energy density and the pressure as the diagonal elements of the stress-energy tensor of an isotropic gas, for an observer locally at rest. Each element $T^{\alpha\alpha}$ represents the flux of the $\alpha$-component of the energy-momentum vectors of all particles in the $\alpha$-direction (see also this question). In other words, we can write
$$
T^{\alpha\alpha}(\vec{x}) = \int v^\alpha p^\alpha\, n(\vec{x},\vec{p})\,\text{d}^3\vec{p},
$$
where $v^\alpha = (c,\vec{v})$ is the velocity vector, $p^\alpha = (E/c,\vec{p})$ is the energy-momentum vector, and
$$n(\vec{x},\vec{p}) = \frac{g}{(2\pi\hbar)^3}\frac{1}{e^{[E(p)-\mu]/kT}\pm 1} =
\frac{g}{(2\pi\hbar)^3}f(\vec{x},p)$$
is the phase-space number density. Therefore,
$$\begin{align}
T^{00} &= \rho c^2 = \frac{g}{(2\pi\hbar)^3}\int E(p)\, f(\vec{x},p)\,\text{d}^3\vec{p},\\
T^{11} &= \frac{g}{(2\pi\hbar)^3}\int v_xp_x\, f(\vec{x},p)\,\text{d}^3\vec{p},\\
T^{22} &= \frac{g}{(2\pi\hbar)^3}\int v_yp_y\, f(\vec{x},p)\,\text{d}^3\vec{p},\\
T^{33} &= \frac{g}{(2\pi\hbar)^3}\int v_zp_z\, f(\vec{x},p)\,\text{d}^3\vec{p}.
\end{align}$$
Since the distribution is isotropic in $p$, we have $T^{11} = T^{22} = T^{33} = P$, so that
$$
P = \frac{1}{3}\left(T^{11} + T^{22} + T^{33} \right) =
\frac{g}{(2\pi\hbar)^3}\int \frac{1}{3}\vec{v}\cdot\vec{p}\; f(\vec{x},p)\,\text{d}^3\vec{p}.
$$
Finally, from
$$
E = \frac{mc^2}{\sqrt{1-v^2/c^2}},\qquad \vec{p} = \frac{m\vec{v}}{\sqrt{1-v^2/c^2}},
$$
we find $\vec{v} = \vec{p}c^2/E$, and
$$
P = \frac{g}{(2\pi\hbar)^3}\int \frac{p^2c^2}{3E(p)}\, f(\vec{x},p)\,\text{d}^3\vec{p}.
$$
