Without loss of generality, we can create spherical coordinates $(p, \theta, \phi)$ such that $\mathbf{r}$ is along the $\theta=0$ axis. This may seem a little backwards, since we usually have $\mathbf{r}$ as the position vector, which naturally varies. But here we leave $\mathbf{r}$ constant, and vary in $\mathbf{p}$-space, since that's what we're integrating over. In particular, variations in $\theta$ and $\phi$ correspond to variations in $\mathbf{p}$ rather than $\mathbf{r}$. Then, we can transform the integral into those coordinates, perform the trivial integral over $\phi$, then change coordinates again using $\sin\theta\, d\theta = -d\cos\theta$, and finally perform that integral too, leaving only the integral over $p$.
\begin{align}
D(x-y) &=
\int_0^\infty \int_0^\pi \int_0^{2\pi} \frac{1}{(2\pi)^3}\,
\frac{1}{2E_{\mathbf{p}}} e^{i\, p\, r\, \cos\theta}\, p^2\, \sin\theta\,
d\phi\, d\theta\, dp \\
&=
\frac{1}{(2\pi)^2}\,
\int_0^\infty \frac{p^2}{2E_{\mathbf{p}}} \left[\int_0^\pi
e^{i\, p\, r\, \cos\theta}\, \sin\theta\,
d\theta\right]\, dp \\
&=
\frac{1}{(2\pi)^2}\,
\int_0^\infty \frac{p^2}{2E_{\mathbf{p}}} \left[\int_1^{-1}
-e^{i\, p\, r\, x}\, dx\right]\, dp \\
&=
\frac{1}{(2\pi)^2}\,
\int_0^\infty \frac{p^2}{2E_{\mathbf{p}}}
\frac{-e^{-i\, p\, r} - -e^{i\, p\, r}}{i\, p\, r}\, dp.
\end{align}
Note that this assumes $E_{\mathbf{p}}$ is independent of $\theta$ and $\phi$, but that's fine because we usually have $E_{\mathbf{p}} = \sqrt{ \lvert \mathbf{p} \rvert^2 + m^2} = \sqrt{ p^2 + m^2}$.
