That can be seen as a consequence of a more general principle: if $\{a_n\}_{n\geq 0}$ is a sequence of complex numbers fulfilling weak constraints, $f(z)=\sum_{n\geq 0}a_n z^n$ is a holomorphic function in a neighbourhood of the origin. Viceversa, given a holomorphic function in a neighbourhood of the origin, its derivatives at the origin can be computed through suitable integrals, as a consequence of Cauchy's integral formula/ the residue theorem. Binomial coefficients are naturally related to shifted monomials, namely through
$$ (1+z)^n = \sum_{k\geq 0}\binom{n}{k}z^k, $$
where the LHS is an entire function. In particular
$$ \binom{n}{k}=\operatorname*{Res}_{z=0}\frac{(1+z)^n}{z^{k+1}}=\frac{1}{2\pi i}\oint_{|z|=1}\frac{(1+z)^n}{z^{k+1}}\,dz = \frac{1}{2\pi}\int_{0}^{2\pi}(1+e^{i\theta})^n e^{-ki\theta}\,d\theta $$
where the equality between the first term and the last one can also be seen as a direct consequence of the orthogonality relation $\int_{0}^{2\pi}e^{ai\theta}e^{-bi\theta}\,d\theta=2\pi\delta(a,b)$. With a bit of maquillage, this identity turns into the given one. In particular, by letting $\theta=2\varphi$ we get:
$$\binom{n}{k}=\frac{1}{\pi}\int_{0}^{\pi}(1+e^{2i\varphi})^n e^{-2ki\varphi}\,d\varphi=\frac{2^n}{\pi}\int_{0}^{\pi}\cos^n(\varphi)e^{(n-2k)i\varphi}\,d\varphi $$
and we are free to replace $e^{(n-2k)i\varphi}$ with its real part, since the LHS is clearly real:
$$\binom{n}{k}=\frac{2^n}{\pi}\int_{0}^{\pi}\cos^n(\varphi)\cos((n-2k)\varphi)\,d\varphi. $$
It is enough to exploit parity and we are done.