Evaluating a double integral of a complicated rational function 
Define the function $Q:\mathbb{C}^{2}\rightarrow\mathbb{C}$ to be the binary quadratic form,
$$Q{\left(z,w\right)}:=z^{2}+w^{2}.\tag{1a}$$
Also, define $P:\mathbb{C}^{4}\rightarrow\mathbb{C}$ to be the polynomial of degree $5$, in four variables,
$$\begin{align}
P{\left(a,b,x,y\right)}
&:=a\left[\left(a^{2}+1\right)\left(a^{2}+b^{2}+1\right)-4b^{2}\right]\\
&~~~~~+2\left[\left(a^{2}+1\right)\left(a^{2}+b^{2}+1\right)-2b^{2}\right]x\\
&~~~~~+a\left(a^{2}+b^{2}+1\right)x^{2}\\
&~~~~~+a\left(a^{2}+b^{2}+1\right)y^{2}.\tag{1b}\\
\end{align}$$
Note that $P{\left(a,b,x,y\right)}$ is obviously even in both $b$ and $y$.
Then, define the function $\mathcal{I}:\mathbb{R}^{2}\rightarrow\mathbb{R}$ via the double integral
$$\mathcal{I}{\left(a,b\right)}:=\int_{-\infty}^{\infty}\mathrm{d}x\int_{0}^{\infty}\mathrm{d}y\,\frac{2^{4}xy\,P{\left(a,b,x,y\right)}}{Q{\left(a+x,1-y\right)}\,Q{\left(a+x,1+y\right)}\,Q{\left(x,b-y\right)}\,Q{\left(x,b+y\right)}}.\tag{1c}$$
It's not hard to show then that $\mathcal{I}{\left(a,b\right)}$ is even in the second parameter $b$:
$$\mathcal{I}{\left(a,-b\right)}=\mathcal{I}{\left(a,b\right)};~~~\small{\left(a,b\right)\in\mathbb{R}^{2}}.$$


Problem: Given the pair of real parameters $\left(a,b\right)\in\mathbb{R}\times\mathbb{R}_{\ge0}$, find a closed form expression for the double integral $\mathcal{I}{\left(a,b\right)}$ in terms of elementary functions.


The obstacle in the way of solving this problem appears to be tedium more than anything else. Integrating the rational integrand in $(1c)$ over $x$ should in principle yield a piecewise rational function. Thus, subsequent integration over $y$ should lead to a function that is at the very least piecewise elementary, if not simpler.
However, attempting to solve the problem by brute force with partial fraction decompositions quickly leads to large numbers of cumbersome expressions, rending the integral quite unmanageable without a program such as Mathematica.
It is my hope that there is some cleverly efficient approach to this integral that I'm just not seeing at the moment. Any advice here would be welcome. Cheers!

 A: This is going to be a partial answer to this question. We will integrate the integrand over $y$ first. In this variable the integrand is clearly a rational function and as such it can be decomposed into partial fractions. Having done this we find the anti-derivative of all relevant terms using the following identity:
\begin{equation}
\int\frac{A+B y}{a+b y+y^2} dy = \frac{2(A-b B)}{\sqrt{4 a-b^2}} \arctan(\frac{b+2 y}{\sqrt{4 a-b^2}}) + \frac{B}{2} \log(a+b y+y^2)
\end{equation}
All we have to do now is to evaluate the resulting expression at infinity and at zero. Using Mathematica and substituting particular numbers for $a$ and $b$ we have checked that the value at infinity vanishes and as such the integral over $y$ boils down to the negative value of our expression at zero. We evaluated this value in Mathematica  and then tediously simplified it by hand. It appears that the resulting expression is very simple and neat. It reads:
\begin{eqnarray}
&&\int\limits_0^\infty \frac{2^{4}xy\,P{\left(a,b,x,y\right)}}{Q{\left(a+x,1-y\right)}\,Q{\left(a+x,1+y\right)}\,Q{\left(x,b-y\right)}\,Q{\left(x,b+y\right)}} dy=\\
&&\frac{2 \left(a^2-b+1\right) \left(\frac{a}{2}+x\right)+a b (b-1)}{b \left(\left(\frac{a}{2}+x\right)^2+\frac{1}{4} (1-b)^2\right)}\cdot\left(\arctan(\frac{b}{x})+\arctan(\frac{1}{x+a})\right)+\\
&&\frac{2 \left(a^2+b+1\right) \left(\frac{a}{2}+x\right)+a b (b+1)}{b \left(\left(\frac{a}{2}+x\right)^2+\frac{1}{4} (b+1)^2\right)}\cdot \left(\arctan(\frac{b}{x})-\arctan(\frac{1}{x+a})\right)+\\
&&\left( \frac{(1-b) \left(a^2-b+1\right)+2 a b \left(\frac{a}{2}+x\right)}{2 b \left(\left(\frac{a}{2}+x\right)^2+\frac{1}{4} (1-b)^2\right)}-\frac{(b+1) \left(a^2+b+1\right)-2 a b \left(\frac{a}{2}+x\right)}{2 b \left(\left(\frac{a}{2}+x\right)^2+\frac{1}{4} (b+1)^2\right)}\right)\cdot(\log(b^2+x^2) - \log(1+(a+x)^2))
\end{eqnarray}
valid for $x,a,b \in{\mathbb R}$. We have checked this equality by randomly sampling values of $a$,$b$ and $x$ and varifying that the right hand side equals the left hand side to high numerical precision.
It is clear that it is not hard to finish our calculation. All we need to do is to reduce the rational function into partial fractions and then find the respective anti-derivatives (which will involve elementary functions and di-logarithms only). We will try to finish this task as soon as possible.
