He jugado un poco con esto de Euler estilo (no preocuparse mucho acerca de establecer la convergencia, que creo que debe ser bastante simple aquí si uno necesita este) y encontró una manera de derivar a partir de cero.
En primer lugar hacer la sustitución de $t=\log(e^{x^2}-1)$ y, a continuación, hacer la segunda sustitución $w=-t$ $t<0$ parte como resultado de la integral para obtener
$$I = \frac{1}{2}\int_0^\infty \frac{dw}{w^2}\left(1 + e^{-w} - 2e^{-w/2}\right)\frac{1}{1+e^{-w}}$$
Ahora expandir la última parte en una serie de Taylor y la introducción (para ser capaces de manipular cada término de la integral de forma independiente) un factor de regularización $w^s$ en el integrando. Esto nos da
$$I(s) = \frac{1}{2}\sum_{n=0}^\infty(-1)^n\int_0^\infty dw\left(1 + e^{-w} - 2e^{-w/2}\right)w^{s-2}e^{-nw}$$
donde $I(0)$ es el valor que nos interesa. Ahora podemos fácilmente evaluar la integral en términos de la $\Gamma$-función,
$$I(s) = \frac{1}{2}\sum_{n=0}^\infty(-1)^n\left[\frac{1}{n^{s-1}} + \frac{1}{(n+1)^{s-1}} - \frac{2^s}{(2n+1)^{s-1}}\right]\Gamma(s-1)$$
Cada término de la suma está bien definido en el límite de $s\to 0$ y el uso de $\lim\limits_{s\to 0}\Gamma(s-1)s = -1$ junto con
$$\lim_{s\to 0} \frac{\frac{1}{n^{s-1}} + \frac{1}{(n+1)^{s-1}} - \frac{2^s}{(2n+1)^{s-1}}}{s} = n\left(2\log(n+1/2)-\log(n+1)-\log(n)\right)\\ + \log(n+1/2)-\log(n+1)$$
nos da
$$I = \frac{1}{2}\sum_{n=0}^\infty(-1)^{n+1}\left[n\left(2\log(n+1/2)-\log(n+1)-\log(n)\right)\\ + \log(n+1/2)-\log(n+1)\right]$$
que después de colapso de los registros, convirtiéndolo en un producto y la simplificación de los rendimientos de la identidad
$$e^{2I-\frac{1}{2}} = \lim_{n\to \infty} \frac{1}{(4n+1)^{2n}}\prod_{k=0}^{n-1}\frac{(4k+3)^{4k+3}}{(4k+1)^{4k+1}}$$
que es idéntica a la que se encuentra aquí (Eq. 53) con
$$I = \frac{K}{\pi}$$
Una prueba de la identidad del producto arriba se puede encontrar aquí.