How was the integral formula for the binomial coefficient discovered?

The trig identity $$\cos(x)^n = 2^{-n}\sum_{k=0}^n {n \choose k}\cos((n-2k)x)$$ combined with the result $$ \int_{-\pi}^\pi \cos(nx)\cos(mx) dx = \delta_{n,m}\pi,$$ if $n$ or $m\ne 0$ and $2\pi$ otherwise, easily gives the integral formula you want.

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.