Calculating the value of a triple integral I would like to evaluate the following integral
$$
I=\int_{\mathbb{R}^3}\frac{x^2+y^2+1-z^2}{|(x,y,z-1)|^3|(x,y,z+1)|^3}\,dx\,dy\,dz.
$$
My attempt: with $R^2=x^2+y^2$ we have
$$
I=2\pi\int_{\mathbb{R}_+\times\mathbb{R}}\frac{R^2+1-z^2}{[(R^2+(z-1)^2)(R^2+(z+1)^2)]^{3/2}}R\,dR\,dz.
$$
The $[\cdots]$ expression in the denominator can be rearranged:
$$
\cdots=R^4+(1-z^2)^2+2R^2z^2+2R^2=(R^2+1-z^2)^2+4R^2z^2.
$$
In addition, one can do another change of variables $D=R^2$, and so
$$
I=\pi\int_{\mathbb{R}_+\times\mathbb{R}}
\frac{D+1-z^2}{[(D+1-z^2)^2+4Dz^2]^{3/2}}\,dD\,dz.
$$
Any thoughts on how to proceed? 

Probably the fastest method: Do not rewrite the denominator and integrate with respect to $D$ directly:
\begin{align}
I &= \pi \int \limits_{-\infty}^{\infty} \int \limits_0^\infty \frac{D + 1 - z^2}{[(D + (1-z)^2)(D + (1 + z)^2)]^{3/2}} \, \mathrm{d} D \, \mathrm{d} z \\
&= \pi \int \limits_{-\infty}^{\infty} \left[\frac{D - (1 - z^2)}{2\sqrt{(D + (1-z)^2)(D + (1 + z)^2)}}\right]_{D=0}^{D=\infty} \, \mathrm{d} z \\
&= \pi \int \limits_{-\infty}^{\infty} \frac{1 + \operatorname{sgn}(1-z^2)}{2} \, \mathrm{d} z = \pi \int \limits_{-1}^1  \, \mathrm{d} z = 2 \pi \, .
\end{align}
However, it seems difficult to come up with this antiderivative without a CAS or a lucky guess.

Alternatively, we can use Feynman parameters to write
\begin{align} I &= \pi \int \limits_{-\infty}^{\infty} \int \limits_0^\infty (D + 1 -z^2)\left[\frac{8}{\pi} \int \limits_0^1 \frac{\sqrt{u(1-u)}}{[(D+(1-z)^2)u + (D+(1+z)^2)(1-u)]^3} \, \mathrm{d} u\right] \, \mathrm{d} D \, \mathrm{d} z \\
&= 8 \int \limits_{-\infty}^{\infty} \int \limits_0^\infty \int \limits_0^1 \frac{(D + 1 -z^2)\sqrt{u(1-u)}}{[D + (1+z)^2 - 4 z u]^3} \, \mathrm{d} u \, \mathrm{d} D \, \mathrm{d} z \, .
\end{align}
Fubini's theorem allows us to integrate with respect to $D$ first, which is now straightforward (partial fractions). We find
\begin{align}
I &= 8 \int \limits_0^1 \sqrt{u(1-u)} \int \limits_{-\infty}^\infty \frac{1 + z - 2 z u}{[(1+z)^2 - 4 z u]^2} \, \mathrm{d} z \, \mathrm{d} u \\
&= 8 \int \limits_0^1 \sqrt{u(1-u)} \int \limits_{-\infty}^\infty \frac{1 + (1-2u) z}{[(z + 1-2u)^2 + 4 u (1-u)]^2} \, \mathrm{d} z \, \mathrm{d} u \\
&= 8 \int \limits_0^1 \sqrt{u(1-u)} \int \limits_{-\infty}^\infty \frac{4 u (1-u) \color{red}{+ (1-2u) \zeta}}{[\zeta^2 + 4 u (1-u)]^2} \, \mathrm{d} \zeta \, \mathrm{d} u 
\end{align}
after introducing $\zeta = z + 1 - 2u$ . The integral of the red term vanishes by symmetry, so the substitution $\zeta = 2 \sqrt{u(1-u)} \tau$ yields
$$ I = 4 \int \limits_0^1 \int \limits_{-\infty}^\infty \frac{\mathrm{d} \tau}{(1+\tau^2)^2} \, \mathrm{d} u = 4 \int \limits_{-\infty}^\infty \frac{\mathrm{d} \tau}{(1+\tau^2)^2} \, .$$
The final integral is a typical exercise in residue calculus, but integration by parts also does the trick:
$$ I = 4 \int \limits_{-\infty}^\infty \left[\frac{1}{1+\tau^2} - \tau \frac{\tau}{(1+\tau^2)^2}\right] \, \mathrm{d} \tau = 4 \left[\pi - \frac{1}{2} \int \limits_{-\infty}^\infty \frac{\mathrm{d}\tau}{1+\tau^2}\right] = 2 \pi \, .$$ 
