$\newcommand{\d}{\,\mathrm{d}}\newcommand{\res}{\operatorname{Res}}$This answer, though long, addresses the keyhole technique(s) and the handling of the original integrand that OP found problematic.
1 - This reasoning is not correct, since the sum of integrals that you provide is: $$\begin{align}\int_{-\infty}^1f(x)\d x+\int_{-1}^\infty f(x)\d x&=\int_{-\infty}^1f(x)\d x+\int_{-1}^1f(x)\d x+\int_1^\infty f(x)\d x\\&=\int_{-\infty}^\infty f(x)\d x+\int_{-1}^1f(x)\d x\\&\overset{?}{=}\int_{-1}^1f(x)\d x\end{align}$$If and only if the integral along $(-\infty,\infty)$ is zero, which is in general going to be false (and I'm pretty sure is false in this particular case).
2 - You mentioned the keyhole contour. There is a lesser-known approach, which I call the "double keyhole contour", which works similarly and is very effective for situations with two poles and a branch cut. It is more appropriate then trying to combine two separate classical keyholes here, in my opinion. Here is a solution using this technique (the method is inspired by the technique of Example 6 from this demonstration). Although I aim to be clear in each step, in all limit computations I am skipping several steps of reasoning, so be sure to read carefully and make sure you agree with me. The core idea is to take two small circular integrals around the poles at $-1,1$, as you might for normal keyhole integration, but you include the branch points and then consider the “residue at infinity”, that is, $-1$ times the residues outside the contour, where the integrand is holomorphic (bar some simple poles) as the branch cut has been trapped inside the contour.
Define $f:\Bbb C\to\Bbb C$ by (for $z$ not a singularity): $$z\mapsto\frac{1}{z^4+1}\exp\left(-\frac{1}{2}\left[\log_1(1-z)+\log_2(1+z)\right]\right)$$Where $\log_1$ is a logarithm defined by $\arg_1(z)\in(-\pi,\pi]$ (the principal branch) and $\log_2$ is a logarithm defined by $\arg_2(z)\in[0,2\pi)$. Individually, $\log_1(1-z)$ and $\log_2(1+z)$ have branch cuts at $[1,\infty)$ and $[-1,\infty)$, and in combination they have only a branch cut of $[-1,1]$ since $f(z)$ is continuous along $(1,\infty)$ - for $x\gt1$ fixed: $$\begin{align}\lim_{y\to0^+}f(x+iy)&=\frac{1}{x^4+1}\cdot\frac{i}{\sqrt{x-1}}\cdot\frac{1}{\sqrt{x+1}}\\\lim_{y\to0^-}f(x+iy)&=\frac{1}{x^4+1}\cdot\frac{-i}{\sqrt{x-1}}\cdot\frac{-1}{\sqrt{x+1}}\end{align}$$And the limits agree. Furthermore, for $x\in(-1,1)$ we have that $f(x)=\frac{1}{(1+x^4)\sqrt{1-x^2}}$.
Take a suitably small $0\lt\varepsilon$. Define $\gamma_{1,\varepsilon}:[\pi/2,3\pi/2]\to\Bbb C,\,t\mapsto-1+\varepsilon\cdot e^{it}$ and define $\gamma_{2,\varepsilon}:[-\pi/2,\pi/2]\to\Bbb C,\,t\mapsto1+\varepsilon\cdot e^{it}$. Let $\gamma_\varepsilon$ be the closed anticlockwise contour taken by the union of $\gamma_{1,2}$ joined with straight line segments between their endpoints. If, for $1\le n\le4$, $\omega_n:=\exp\left(\frac{\pi n}{2}i-\frac{\pi}{4}i\right)$ the fourth roots of $-1$, then by the residue theorem: $$\int_{\gamma_{1,\varepsilon}+\gamma_{2,\varepsilon}}f(z)\d z+\int_{-1-i\varepsilon}^{1-i\varepsilon}f(z)\d z-\int_{-1+i\varepsilon}^{1+i\varepsilon}f(z)\d z=-2\pi i\cdot\left[\res_{z=\infty}+\sum_{n=1}^4\res_{z=\omega_n}\right]$$
Since the outside of the contour does not contain any branch cut - the integrand is meromorphic - and there are only four simple poles at the $\omega_n$.
A simple analysis of the semicircular integrals around $\gamma_{1,2}$ shows they are bounded by at least $O(\sqrt{\varepsilon})$ and vanish as $\varepsilon\to0^+$. We also see quite straightforwardly that: $$\int_{-1-i\varepsilon}^{1-i\varepsilon}f(z)\d z-\int_{-1+i\varepsilon}^{1+i\varepsilon}f(z)\d z\overset{\varepsilon\to0^+}{\longrightarrow}-2I$$From the choices of logarithm branch. I don't know much about the mysterious "residue at infinity", but since $\lim_{|z|\to\infty}z\cdot f(z)=0$ I am fairly certain the residue is here equal to $0$.
I compiled some tables to help with the mass computation of the residues at each $\omega_n$. This cleans things up a lot, but makes the text longer.
Since these are all simples poles and the exponential/logarithm is analytic near each $\omega_n$, we have: $$\begin{align}\res_{z=\omega_n}f(z)&=\exp\left(-\frac{1}{2}\left[\log_1(1-z)+\log_2(1+z)\right]\right)\cdot\res_{z=\omega_n}\frac{1}{1+z^4}\\&=\exp\left(-\frac{1}{4}\cdot\ln(|1-\omega_n^2|^2)\right)\exp\left(-\frac{i}{2}[\arg_1(1-\omega_n)+\arg_2(1+\omega_n)]\right)\cdot\frac{1}{4\omega_n^3}\end{align}$$We have: $$\begin{align}|1-\omega_n^2|^2&=\begin{cases}|1-i|^2=2&n=1,3\\|1+i|^2=2&n=2,4\end{cases}\\\\\arg_1(1-\omega_n)&=\begin{cases}-\arctan(\sqrt{2}+1)&n=1\\-\arctan(\sqrt{2}-1)&n=2\\\arctan(\sqrt{2}-1)&n=3\\\arctan(\sqrt{2}+1)&n=4\end{cases}\\\\\arg_2(1+\omega_n)&=\begin{cases}\arctan(\sqrt{2}-1)&n=1\\\arctan(\sqrt{2}+1)&n=2\\2\pi-\arctan(\sqrt{2}+1)&n=3\\2\pi-\arctan(\sqrt{2}-1)&n=4\end{cases}\\\\\frac{1}{4\omega_n^3}&=\begin{cases}\frac{-1-i}{4\sqrt{2}}&n=1\\\frac{1-i}{4\sqrt{2}}&n=2\\\frac{1+i}{4\sqrt{2}}&n=3\\\frac{-1+i}{4\sqrt{2}}&n=4\end{cases}\end{align}$$If we further note that $\arctan(\sqrt{2}+1)-\arctan(\sqrt{2}-1)=\arctan(2/2)=\pi/4$ by the arctangent difference formulae, we see that the residues are: $$\begin{align}\res_{z=\omega_1}f(z)&=\frac{-1-i}{4\sqrt{2}\cdot\sqrt[4]{2}}\cdot\exp(\pi i/8)\\\res_{z=\omega_2}f(z)&=\frac{1-i}{4\sqrt{2}\cdot\sqrt[4]{2}}\cdot\exp(-\pi i/8)\\\res_{z=\omega_3}f(z)&=\frac{1+i}{4\sqrt{2}\cdot\sqrt[4]{2}}\cdot\exp(\pi i/8)(-1)\\\res_{z=\omega_4}f(z)&=\frac{-1+i}{4\sqrt{2}\cdot\sqrt[4]{2}}\cdot\exp(-\pi i/8)(-1)\end{align}$$
The sum of residues is then: $$\begin{align}\frac{1}{4\sqrt{2}\cdot\sqrt[4]{2}}\cdot(-2(1+i)\exp(\pi i/8)+2(1-i)\exp(-\pi i/8))&=\frac{1}{2\sqrt{2}\cdot\sqrt[4]{2}}\cdot(-2i\sin(\pi/8)-2i\cos(\pi/8))\\&=\frac{-i}{2^{3/4}}(\sin(\pi/8)+\cos(\pi/8))\end{align}$$
Therefore we get to conclude that: $$I=\frac{\pi}{2^{3/4}}(\sin(\pi/8)+\cos(\pi/8))=\frac{\pi}{2}\sqrt{1+\sqrt{2}}$$Which is the correct answer.
